Soil Moisture-Temperature Coupling in a Set of Land Surface Models.

The land surface controls the partitioning of water and energy fluxes and therefore plays a crucial role in the climate system. The coupling between soil moisture and air temperature, in particular, has been shown to affect the severity and occurrence of temperature extremes and heat waves. Here we study soil moisture-temperature coupling in five land surface models, focusing on the terrestrial segment of the coupling in the warm season. All models are run off-line over a common period with identical atmospheric forcing data, in order to allow differences in the results to be attributed to the models' partitioning of energy and water fluxes. Coupling is calculated according to two semiempirical metrics, and results are compared to observational flux tower data. Results show that the locations of the global hot spots of soil moisture-temperature coupling are similar across all models and for both metrics. In agreement with previous studies, these areas are located in transitional climate regimes. The magnitude and local patterns of model coupling, however, can vary considerably. Model coupling fields are compared to tower data, bearing in mind the limitations in the geographical distribution of flux towers and the differences in representative area of models and in situ data. Nevertheless, model coupling correlates in space with the tower-based results (r = 0.5-0.7), with the multimodel mean performing similarly to the best-performing model. Intermodel differences are also found in the evaporative fractions and may relate to errors in model parameterizations and ancillary data of soil and vegetation characteristics.

Introduction

The interaction between the land surface and atmosphere is an important driver of climate, since it involves an exchange of energy, water, and chemical elements such as carbon. Soil moisture, in particular, links the water and energy cycles through its control on evaporation and has been shown to influence precipitation and temperature (e.g., Catalano et al., 2016; Koster et al., 2006; Seneviratne et al., 2010; Taylor et al., 2012). The link from soil moisture to temperature occurs when limited soil moisture availability restricts the amount of energy used for evaporation, leaving more energy available for heating. This effect has been related to the occurrence and persistence of extreme temperatures and heat waves (Hirschi et al., 2011; Lorenz et al., 2010; Miralles et al., 2014). Due to the severe socioeconomic impacts of heat waves (Easterling, 2000; García‐Herrera et al., 2010), and the expectation that their frequency and intensity may increase as we progress into the future (Meehl & Tebaldi, 2004; Quesada et al., 2012; Tebaldi et al., 2006), it is important to understand the strength of soil moisture‐temperature coupling in models and whether this coupling strength appears realistic in comparison to observations.

Recognition of the importance of the interaction between soil moisture and temperature has led to a number of studies at regional and global scales, using model simulations (Dirmeyer, 2011; Koster et al., 2006; Seneviratne et al., 2006), observational data sets (Hirschi et al., 2011, 2014; Miralles et al., 2014; Mueller & Seneviratne, 2012), and combinations of these (Li et al., 2017; Zscheischler et al., 2015). These studies have employed metrics of varying complexity, ranging from relatively simple correlation‐based indices (Knist et al., 2017; Li et al., 2017; Miralles et al., 2012; Seneviratne et al., 2006; Zscheischler et al., 2015) to more complex statistical inferences based on wavelets theory (Casagrande et al., 2015) or more specific model diagnostics (Fischer, Seneviratne, et al., 2007; Koster et al., 2006). These studies found that transitional regions between wet and dry climates are typically hot spots of soil moisture‐temperature coupling: in these regions, evaporation is both controlled by soil moisture availability and variable and large enough to have an impact on atmospheric dynamics. In wetter regions, soil moisture is not limiting and evaporation tends to be radiation driven; in drier regions, the variability of soil moisture and evaporation becomes too low to have an impact on the dynamics of air temperature. Despite this consistency in the location of global hot spots of soil moisture‐temperature coupling found by previous studies, there are considerable differences in the magnitudes and regional patterns of coupling that have been reported (Dirmeyer, 2011; Koster et al., 2006).

Understanding model differences in the representation of soil moisture‐temperature coupling is crucial, as inaccurate representation of such land‐atmosphere interactions may result in errors and biases in climate projections (Davin et al., 2016; Fischer et al., 2012). At the same time, a better representation of surface soil moisture in atmospheric models has shown to improve weather forecasts (e.g., Bisselink et al., 2011; Orth et al., 2016; Quesada et al., 2012; van den Hurk et al., 2012) and constrain predictions of future climate variability (Sippel et al., 2016; van den Hurk et al., 2016; Vogel et al., 2017). This has further implications as land surface models (LSMs) are increasingly used to assess drought conditions (Dai, 2013; Prudhomme et al., 2014; Ukkola et al., 2016) and develop early‐warning systems (McNally et al., 2017). However, intercomparing results from past soil moisture‐temperature coupling studies is not straightforward: they are not only commonly based on a single land surface model and coupling metric, and rarely contrasted against observational data, but are also affected by the particular choice of atmospheric forcing. These methodological inconsistencies imply that intermodel differences cannot be attributed to the skill of particular land surface models and highlight the need for a thorough intercomparison and validation of different models and coupling metrics based on common atmospheric forcing.

In this study, we focus on the terrestrial segment of the coupling. We apply two different metrics of soil moisture‐temperature coupling to a set of five land surface models and compare the results to observational data. These five models have been run off‐line with consistent forcing within the frame of the eartH2Observe project (http://www.earth2observe.eu), and thus intermodel differences can be attributed to the way in which particular models represent the link between atmospheric state and surface energy partitioning. In the next section, we describe the metrics and data in more detail. In section 3, we present and discuss the global coupling hot spots, an intermodel comparison of global coupling and a comparison of the model results to observational data. Finally, we conclude with a summary in section 4.

Methods

Metrics

The metrics applied in this study do not explicitly take soil moisture into account, but infer soil moisture's control on evaporation. The first metric is simply the correlation between latent heat flux and temperature (Seneviratne et al., 2006); to aid interpretability, this metric (here referred to as Χ) is presented as the negative of this correlation: where ρ is the Pearson correlation coefficient, E the surface latent heat flux, and T the near‐surface air temperature. The theoretical range of this metric is from −1 to 1. Negative values indicate areas where evaporation is energy limited. Positive values, on the other hand, indicate areas where high (low) evaporation corresponds to low (high) air temperature, thus where soil moisture is potentially limited and acting upon evaporation and temperature dynamics. The more positive the value, the stronger the soil moisture‐temperature coupling. This metric has previously been used to study soil moisture‐temperature coupling in regional and global circulation models and how coupling may change in the future as climate changes (Seneviratne et al., 2006).

The second metric is based on the surface radiation balance equation and was first described by Miralles et al. (2012). This metric focuses on the difference resulting from correlating air temperature to sensible heat, depending on whether the sensible heat is computed using actual evaporation or potential evaporation. In this way, the metric becomes less dependent on the potential confounding effects caused by radiation variability (due to, e.g., clouds) on both evaporation and temperature. The metric (П) is defined as where R n is surface net radiation, λE the surface latent heat flux, λE p the potential latent heat flux, and T the near‐surface air temperature. The ground heat flux is not included in this equation, as this term tends to be comparatively low at daily time scales and has been shown to have a limited effect on the metric (Miralles et al., 2012). In this way, the metric is positive when taking the actual soil water availability into account and helps explain a larger fraction of the air temperature variability, compared to the case in which no constraints due to water availability exist. Where water availability is not a limiting factor for evaporation, the first and second terms will become equal and the value of the metric will approach zero. The theoretical range for the metric is from −2 to 2, although typical values range from 0 to 1, with higher values indicating stronger coupling. The applicability of the metric, however, relies on the availability of λE p, which is not explicitly computed and reported by all land surface models (Schellekens et al., 2016). In addition, λE p is very sensitive to the method used to calculate it (Fisher et al., 2011; Sperna Weiland et al., 2015). Therefore, assuming that λE p variability is mostly dictated by the variability in R n (Priestley & Taylor, 1972), here we simplify equation (2) as follows:

Though the difference between the two terms is no longer only the availability of soil moisture, the basic rationale of the metric is still applicable. In the case of deserts, for example, where evaporation is negligible, the first term will be equal to the second term and the final result of the metric will be zero. On the other end of the spectrum, where soil moisture is high and evaporation is primarily limited by radiation, λE will be a close function of R , thus both correlations will have a similar magnitude and the metric will again approach zero. Where soil moisture limits evaporation, and this limitation leads to impacts on the variability of sensible heat flux that show an imprint on the dynamics of air temperature, we expect positive values of П in equation (3).

The Χ and Π metrics evaluate the strength of soil moisture‐temperature coupling throughout the period of interest, usually the summer season. In an additional analysis, we compare the global ensemble mean results of these metrics to a metric that focuses on extremes, namely the Vegetation‐Atmosphere‐Coupling (VAC) index. This metric was developed to evaluate the relationship between extremes in temperature and photosynthetic activity but has also been applied to evaporation (Zscheischler et al., 2015). Here we use the evaporation formulation of the VAC index, which identifies the concurrence of evaporation and temperature extremes based on the 70th percentiles of the absolute values of these variables. We sum the relative occurrences of extremes with opposite sign (i.e., high temperature with low evaporation, and vice versa) as a measure of soil moisture‐temperature coupling, as these are more likely to occur in water‐limited regimes where soil moisture controls evaporation and temperature. For more details on this metric, see Zscheischler et al. (2015).

Note that these semiempirical metrics are based on the assumption that coupling between surface fluxes and the atmosphere is dominated by soil moisture. While soil moisture is an important factor, in reality other factors can play a role. Depending on the region and time scales of study, factors such as land cover and land use change, vegetation variability, and variability in radiation (relevant for Χ, but largely accounted for in Π) will also impact surface fluxes.

Models and Data

Five models from the eartH2Observe project were included in this study. Four of them are LSMs, which are designed to simulate the exchange of water and energy between soil, vegetation, and atmosphere: the Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (HTESSEL) (Balsamo et al., 2009), the Joint UK Land Environment Simulator (JULES) (Best et al., 2011; Clark et al., 2011), the ORganizing Carbon and Hydrology in Dynamic EcosystEms (ORCHIDEE) (D'Orgeval et al., 2008), and the Interaction between Soil Biosphere Atmosphere model in the SURFace EXternalized (SURFEX) modeling platform (Decharme et al., 2013). The fifth model, the Global Land Evaporation Amsterdam Model (GLEAM) (Martens et al., 2016; Miralles et al., 2011), is not a traditional LSM in that its primary purpose is to derive global evaporation from satellite observations, but it includes many of the same processes as traditional LSMs. In the eartH2Observe project, all models (including several global hydrological models which are not designed to simulate land‐atmosphere interactions and which are therefore not included in this study) are run with a consistent atmospheric forcing database. Therefore, these five models are not coupled to an atmospheric model but they are run in off‐line mode. There was no prescribed set of static fields, such as land cover types and soil parameters, and each modeling group chose the most suitable data sets for their modeling system (see Table 1). For GLEAM, satellite forcing data that were not part of the atmospheric forcing data set—such as the surface soil moisture and vegetation optical depth microwave retrievals—were also retained. None of the studied models account for land cover or land use change throughout the study period.

Table 1

Overview of Important Model Characteristics

Model	Evaporation scheme	Land cover source	Vegetation classes	Soil layers	Soil depth (m)
HTESSEL	Penman‐Monteith	Global land cover characteristic (Loveland et al., 2000)	Tall and short vegetation	4	2.89
JULES	Penman‐Monteith	Global land cover characteristic v2	5 classes: broadleaf forest, needleleaf forest, C3 Grass, C4 grassland, shrubs	4	3
ORCHIDEE	Bulk method (Barella‐Ortiz et al., 2013)	IGBP map with Olson classification (de Rosnay & Polcher, 1998)	13 vegetation types for transpiration and interception loss, grouped into 3 ensembles (tall/short vegetation and bare soil) for throughfall and root uptake	11	2
SURFEX	Penman‐Monteith	ECOCLIMAP (Faroux et al., 2013)	12 plant functional types including broadleaf forest, needleleaf forest, C3 crops, C4 crops, grassland and bare soil	14	Variable, up to 12 m
GLEAM	Priestley‐Taylor	MODIS global vegetation continuous fields (Hansen et al., 2005)	Tall and short vegetation	3	2.5

The common forcing database is the WATCH Forcing Data methodology applied to ERA‐Interim reanalysis data (WFDEI) (Weedon et al., 2014). The forcing and model outputs have a common 0.5° resolution and cover the period 1979–2012. The T forcing and model‐derived R and λE are used to calculate the coupling metrics for each model using equations (1) and (3). Relevant to this study is that these models derive evaporation based on different methods: two surface energy balance approaches, namely, the Penman‐Monteith equation (HTESSEL, JULES, and SURFEX) and a modified Priestley and Taylor equation (GLEAM), and a turbulent diffusion approach (ORCHIDEE). More information on the ensemble of models in the eartH2Observe project can be found in Table 1 and Schellekens et al. (2016).

The fact that the models are run off‐line has consequences for the results of the coupling metrics. Soil moisture‐temperature coupling, and indeed other land‐atmosphere interactions, occur in two stages (Dirmeyer, 2011; Guo et al., 2006). In the first step, the land surface affects surface fluxes such as evaporation and sensible heat flux. In the second step, the surface fluxes affect atmospheric states and fluxes, such as precipitation and temperature. In this study, the prescribed temperature forcing means that errors in the simulation of evaporation do not propagate to the temperature; thus, the approach focuses on the terrestrial segment of soil moisture‐climate coupling (Dirmeyer, 2011). Nevertheless, previous work found that a large part of differences in model land‐atmosphere coupling occurs in this first step between soil moisture and evaporation (Guo et al., 2006). Thanks to this approach we are able to attribute any differences to the models' partitioning of energy and water fluxes—as opposed to differences in model forcing—while constraining the models with realistic, observation‐corrected, atmospheric dynamics.

Eddy covariance data from the FLUXNET2015 synthesis data set (http://fluxnet.fluxdata.org/) were used as a benchmark in this study. Only sites with at least 100 observations in the modeled period and an energy balance mismatch of less than 20% were included. The energy balance mismatch is calculated as the difference between R and the sum of latent, sensible, and ground heat fluxes measured at each tower. These criteria led to the selection of 59 sites spread over different continents and climates. The list of stations is presented in Table S1 (Amos et al., 2005; Anthoni et al., 2004; Ardö et al., 2008; Aubinet et al., 2001; Baldocchi et al., 2016; Beringer, Hacker, et al., 2011; Beringer, Hutley, et al., 2011; Bristow et al., 2016; Cleverly et al., 2013; Cook et al., 2004; Fischer, Billesbach, et al., 2007; Goldstein et al., 2000; Imer et al., 2013; Knohl et al., 2003; Leuning et al., 2005; López‐Ballesteros et al., 2017; Ma et al., 2007; Merbold et al., 2014; Meyer et al., 2015; Milyukova et al., 2002; Noormets et al., 2008; Scott et al., 2015; Serrano‐Ortiz et al., 2009; Suni et al., 2003; Verma et al., 2005; Zeeman et al., 2010). The metrics described in section 2.1 were calculated using both the measured λE as well as the λE calculated based on the energy balance residual at each tower (i.e., the difference between observed net radiation and the sum of the observed ground and sensible heat fluxes), to give an indication of the uncertainty in the eddy covariance data (i.e., Michel et al., 2016). For each of the tower sites, we calculated a representativeness measure based on the percentage of the 0.5° pixel covered by the same land cover classification as defined for the tower site. In this study, we classified sites by a tree/nontree cover classification and used the MODIS44b 250 m Vegetation Continuous Fields product (Hansen et al., 2005) to calculate the cover percentages at 0.5° resolution.

In this study, we focused on the 1979–2012 period covered by the eartH2Observe model simulations, and only on the warm season, because we assume soil moisture‐temperature coupling (and related hazards such as high‐temperature extremes) to be most relevant during summer. For each model pixel, the warm season is defined as the 3 month period centered on the month with the highest average temperature calculated over all years. This implies that the results shown in maps do not correspond to the same timestamp in every pixel, though the warm season is largely consistent with summertime outside the tropics (see Figure 1). The annual variability in the hottest month is generally low, as the standard deviation of the difference between the month used to define the warm season and the annual hottest months is smaller than 1 month in 92% of pixels. Higher‐standard deviations are found in tropical regions where soil moisture‐temperature coupling is not expected (Figure 1). Before calculating the metrics, all model and tower data were converted to anomaly time series by subtracting their climatological mean before applying equations (1) and (3). This climatological mean was calculated for each day as the average of the data falling within a 31 day window centered on that day of the year, over all available years of data. Finally, only days with no precipitation were considered to avoid the confounding effects of interception loss: evaporation of water intercepted by vegetation is independent from soil moisture and can introduce artifacts when assessing coupling using simple metrics (Guillod et al., 2014). Furthermore, eddy covariance data are unreliable under rainy conditions. In this study, precipitation days are defined as days with nonzero precipitation in either the model forcing data set or the observational tower data, or both. The percentage of rainy days excluded from the model analysis is shown in Figure S1.

Figure 1

(top) The month of maximum average air temperature (1979–2012) used to define the warm season and (bottom) the standard deviation of the difference between that month and annual months of maximum air temperature.

Results and Discussion

In this section, we first assess the global patterns of soil moisture‐temperature coupling according to the two semiempirical metrics described in section 2.1 (section 3.1). Then, the variability in model coupling strength is evaluated and the important regional differences are discussed (section 3.2). Finally, the model results are compared to eddy covariance data to evaluate whether the coupling strength in models is consistent with observational data (section 3.3).

Coupling Hot Spots

The global hot spots of soil moisture‐temperature coupling are shown in Figure 2. In this figure, areas with no coupling (i.e., values below zero in Χ and П) have been set to zero before calculating the multimodel mean to focus only on the areas where coupling is relevant. The location of these hot spots is consistent with previous studies (Dirmeyer, 2011; Koster et al., 2006; Miralles et al., 2012; Seneviratne et al., 2010) and highlights the strong coupling in transitional climate zones, such as northern Australia, India, southern Africa, the Great Plains in the United States, and, to a lesser extent, the Sahel.

Figure 2

The (top row) multimodel mean and (bottom row) standard deviation of soil moisture–temperature coupling according to (a) Χ, (b) П, and (c) the VAC index for the period 1979–2012.

The global patterns of soil moisture‐temperature coupling according to the Χ and П metrics are remarkably similar (spatial Pearson correlation of 0.94), though there are small regional differences. For example, coupling strength in subequatorial Africa is more concentrated in the south according to П, while it is more evenly distributed in Χ. It is important to note that, although both metrics are shown on the same scale, we cannot directly compare their magnitudes, as their theoretical basis and potential range are different. Specifically, the main difference between both metrics relates to their treatment of radiation, and the extent to which they manage to isolate the effect of soil moisture on temperature. Differences between the metrics are more apparent when we look at intermodel variability (spatial Pearson correlation of 0.76). The variability in Χ is highest in areas where coupling is relatively weak, such as the Sahel and Southeast Asia; this suggests that model variability is related to the spatial extent of the coupling hot spots rather than their magnitude. The variability in П, on the other hand, is also found in areas with strong coupling, such as in northern Australia and India.

To evaluate whether overall seasonal coupling strength is comparable to coupling strength focusing on extremes, we compare the results of Χ and П to the VAC index (Figure 2). Global patterns of coupling are similar to those of Χ and П, though the VAC index identifies strong coupling in Southeast Asia and along the northern coast of South America that is not observed in the other metrics. The agreement between the VAC index and Χ (spatial Pearson correlation of 0.89) is somewhat better than the agreement with П (spatial Pearson correlation of 0.78), which may be because the VAC index and Χ are both based on evaporation and temperature anomalies alone, while П also accounts for variability in net radiation. These results show that the global patterns of soil moisture‐temperature coupling strength are consistent whether the metrics focus on all available data or only on data extremes. In the rest of the manuscript, we focus on soil moisture‐temperature coupling strength using the analogous Χ and П metrics.

The hot spots illustrated here are similar to the regions highlighted by previous studies. Hot spots in the central United States, the Sahel, and India were also identified in studies investigating coupling in boreal summer (Koster et al., 2006; Miralles et al., 2012; Seneviratne et al., 2010), while hot spots in southern Africa and Australia are found in studies that explore the austral summer (Dirmeyer, 2011; Dirmeyer et al., 2013; Miralles et al., 2012). The hot spots of coupling are different from those presented in Zscheischler et al. (2015), mainly because the hot spots in the Northern Hemisphere are usually weaker or missing altogether in that study. The differences are likely because the VAC index was based on photosynthetic activity in the before‐mentioned study and the results were not limited to the dry days of the warm season. Overall, the coupling hot spots found in this study are most similar to Miralles et al. (2012) and Seneviratne et al. (2010), mainly due to the absence of strong coupling in eastern China (observed by Koster et al., 2006) and Argentina (observed by Dirmeyer, 2011). The greater similarity to Miralles et al. (2012) and Seneviratne et al. (2010) is not surprising considering that the metrics applied here were also used in those studies. Even so, there are differences in the extent of many of the hot spots and coupling in the Mediterranean region, for example, is notably weaker in this study. These differences can be attributed to differences in the study periods and the data sets analyzed in this study, as well as the modification of the П metric, and the distinction between online and fully coupled models.

Finally, one interesting feature in Figure 2 is the horizontal pattern visible in all metrics in the Sahel region. This pattern is caused by changes in the month of the highest temperature we used to define the warm season, which shifts from March in the southern Sahel to July in the Sahara (Figure 1). The fact that there is a north‐south gradient in coupling strength within each hottest month region suggests that the period of strongest coupling in the Sahel does not coincide with the warm season we defined by the month of maximum temperature, but is slightly lagged. This lag could be caused by the onset of the rainy season in the Sahel, which follows the month of maximum temperature by a month or two. The increase in soil moisture variability at the onset of the wet season increases soil moisture‐temperature coupling strength compared to hottest months where soil moisture and its variability are low. Thus, the coupling strength is lower when evaluating the dry season than when the study period includes the (onset of) the wet season. Previous studies have described a similar latitudinal shift in rainfall, temperature (Domínguez et al., 2010), and soil wetness (Taylor, 2008), as well as in the sensitivity of the latent heat flux to surface soil moisture (Dirmeyer, 2011) in the Sahel region. Note that there are likely more locations where the season with the highest soil moisture‐temperature coupling strength does not coincide with the warm season, but that we focus on this period because we assume it is the period that the related hazards such as temperature extremes are most relevant.

Intermodel Comparison

There are considerable intermodel differences in the magnitude and regional patterns of soil moisture‐temperature coupling. In Figure 3, we show global coupling based on Χ for each of the individual models, as well as the anomalies from the model mean. The latitudinal profiles show the mean coupling strength over the land pixels in each row, corresponding to 0.5° latitude. Figure 4 shows analogous results for П. Coupling strength in JULES according to both metrics is relatively high compared to the other models. Indeed, it is the only model in which part of the Amazon region is considered to be a coupling hot spot. The magnitudes of the coupling metrics for GLEAM, on the other hand, are relatively low. GLEAM points to the same global hot spots shown in Figure 2, but some areas with low to intermediate magnitudes of coupling are not shown, such as Argentina and southeastern United States. Overall, the patterns of positive and negative anomalies in Figures 3 and 4 indicate many intermetric similarities: an overestimation/underestimation of coupling (relative to the multimodel mean) according to one metric often corresponds to an overestimation/underestimation in the other metric. The most notable exception is HTESSEL, which showed relatively weak coupling globally based on Χ, but shows a patchwork of both positive and negative anomalies from the model mean according to П, which may relate to differences in sensitivity to radiation.

Figure 3

Patterns and latitudinal profiles of soil moisture‐temperature coupling based on (a) Χ for each model, and the same for (b) the deviation from the multimodel mean.

Figure 4

Same as Figure 3 but for П.

Explaining the observed differences in coupling strength between models is not straightforward. Since atmospheric forcing is the same for all models, the differences in coupling strength can be attributed to differences in the dynamics of the latent heat flux and their correspondence to the air temperature. In models, soil moisture is linked to transpiration and air temperature by the vegetation stress curve and/or through the parameterization of surface conductance. These are influenced by many factors in models, including ancillary data and choice of algorithms and parameterizations (Medlyn et al., 2015); in particular, the shape of the vegetation stress curve (Combe et al., 2016; Powell et al., 2013), simulated plant water use strategies (Kala et al., 2016), and soil and vegetation characteristics (e.g., soil texture, root depth, and land cover classification). The current experiment with consistent forcing is not sufficient to disentangle all these factors. A comprehensive set of experiments with different model parameterizations and/or structures would be required to understand which factors are at the root of the intermodel differences, such as done for ecosystem response to CO2 by Medlyn et al. (2015).

As an example, one factor that could be expected to play an important role is the type of evaporation algorithm used in each model. For instance, the low coupling strength in GLEAM in the humid temperate regions of southeastern United States and northern Argentina coincides with areas where potential evaporation rates are higher according to the Priestley‐Taylor than to Penman‐Monteith formulation (Fisher et al., 2011; Sperna Weiland et al., 2015). However, the differences between the two evaporation schemes are expected to be larger in dry regions where the vapor pressure deficit can be critical as a driver of evaporation, including the coupling hot spot regions of northern Australia and southwestern United States (Fisher et al., 2011; Sperna Weiland et al., 2015). In addition, the three models that are based on a Penman‐Monteith evaporation scheme (HTESSEL, JULES, and SURFEX) also show a large intermodel variability. JULES, for example, shows especially strong soil moisture‐temperature coupling for both metrics in tropical areas such as the Amazon, south‐east Asia and tropical Africa. This is in line with previous studies reporting an underestimation of evaporation rates in the Amazon in JULES (Blyth et al., 2011; Zulkafli et al., 2013), which was later attributed to an overestimation of soil water stress (Zulkafli et al., 2013). Alternatively, the fact that the overestimation of coupling is concentrated in the arc of deforestation in the Amazon may indicate that land cover parameterizations play an important role there. In contrast, SURFEX has relatively weak coupling in the Amazon, but exceptionally strong coupling in the colder climates of central Asia, where HTESSEL shows weak coupling. Finally, the ORCHIDEE model uses the bulk method to estimate evaporation, which uses a diffusive equation based on the air density, the humidity gradient, and the aerodynamic resistance (Barella‐Ortiz et al., 2013). Coupling strength in ORCHIDEE, falls mostly within the already described model variability, with weaker coupling in tropical climates and stronger coupling in parts of eastern Asia and the arctic regions. Thus, under the current experimental setup, it is not possible to identify a systematic difference in soil moisture‐temperature coupling strength depending on the model's choice of evaporation scheme.

In order to gain more insight into the differences in soil moisture‐temperature coupling, we compare the relation between coupling strength and evaporative fraction (EF), here defined as the ratio of latent heat over net radiation. Density plots of Χ and П against EF in Figure 5a show an asymmetrical, positively skewed relationship for all models, with the maximum coupling strength occurring at evaporative fractions ranging 0.1–0.2. The shape of the distribution near the maximum is rounder for Χ and sharper for Π. Though there appears to be a good relation between maximum coupling strength and EF for both metrics, there is considerable scatter toward lower Χ and П, as a result of both wet and arid regions presenting low coupling strength. Of the studied models, the distribution is most skewed for SURFEX and least skewed for GLEAM. For the latter, the ranges of the evaporative fraction and coupling strength are noticeably narrower than for the LSMs.

Figure 5

Density plots of coupling strength against evaporative fraction (a) using the actual values and (b) the difference from the multimodel mean. In Figure 5b, pixels are first selected by the criterion that the coupling strength should be larger than zero for at least three of five models. Then, the 5% driest pixels (red) and the 5% wettest pixels (blue) are illustrated separately; this classification is based on the mean climatology of precipitation at each pixel. Color shades are a measure of point density plotted on a natural log scale.

The shape of the relationship between coupling strength and EF in Figure 5a is consistent with what is expected based on the current state of knowledge, that is, that coupling is strongest in transitional climates. Furthermore, it means that wet and dry regions respond differently to errors in model EF. For example, an overestimation of EF in a dry region is an indication of lower evaporative stress, and thereby signifies a shift in model behavior toward a less dry and more transitional climate, which in turn corresponds to an increase in soil moisture‐temperature coupling in the absence of confounding factors. In a wet region, on the other hand, an overestimation of EF signifies a shift away from a transitional climate, and therefore a shift toward weaker soil moisture‐temperature coupling. Note that in these off‐line model simulations there is no actual shift in climate, but differences in EF indicate that model behavior in a certain region is more (or less) similar to a transitional climate. Figure 5b quantifies the strength of this relationship between EF and coupling strength in dry and wet climates for all models. First, based on the average precipitation in the warm season, we selected the 5% driest and wettest pixels where at least three out of five models identified positive values of soil moisture‐temperature coupling. Then, EF and the coupling strength of these pixels were plotted as anomalies relative to the multimodel mean EF and coupling strength, respectively. In drier regions, models with higher (lower) EF tend to present a higher (lower) soil moisture‐temperature coupling; the opposite occurs for wet regions. This relationship is stronger for the “wet” pixels than for the “dry” pixels, which could be a result of the relatively narrow range in evaporative fraction represented by the dry pixels.

The strength of the relationship between EF and Χ or П also varies from model to model. JULES and SURFEX, which generally show stronger coupling (Figures 3 and 4), present higher correlations between EF and the coupling metrics for dry (positive) and wet (negative) regions. On the other hand, GLEAM and HTESSEL, which generally show weaker coupling (Figures 3 and 4), present weaker correlations between EF and the coupling metrics. For ORCHIDEE, the negative relationship in wet climates is strong for both metrics, but the relationship between EF and coupling strength in dry climates is very weak for Χ, but strong for Π. It is important to note that only the tails of the distribution in rainfall are shown in Figure 5b; additional testing indicates that (weak) positive correlations persist only until the 10th rainfall percentile, then there is an abrupt transition to negative correlations above the 15th percentile (not shown). This is consistent with what would be expected based on the distributions in Figure 5a.

Tower Benchmarking

The model soil moisture‐temperature coupling is evaluated against coupling at 59 FLUXNET sites that have been sampled based on their availability of data (see section 2.2). At each pixel in which tower data are available, the coupling metrics for each model were recalculated for those days with in situ data. In other words, tower data were compared to model data from the 0.5° pixel containing the tower location. While thus far we focused on areas where soil moisture‐temperature coupling is relevant (i.e., positive values for the metrics), here we extend the analysis to include sites with negative values of the metrics; after all, these values can still indicate whether the dynamics of the model and tower data are similar. In addition, we study two characteristic sites in more detail, one with high coupling strength and one with low coupling strength (see below). Note that the level of agreement between model and tower data does not reflect the ability of the metrics to identify regions of soil moisture‐temperature coupling, but rather reflects whether the dynamics of the underlying variables are similar for models and measurements.

In Figure 6, model and tower soil moisture‐temperature coupling strength are compared. First, the coupling at the tower sites is compared to the multimodel mean coupling strength for the period 1979–2012 and illustrated in a global map, and second, each of the models is compared to the tower inferences separately. Performance metrics for this comparison are reported in Table 2. Model and tower Χ are consistent overall, with Pearson correlations ranging between 0.55 (ORCHIDEE) and 0.69 (HTESSEL), though there can be substantial differences at specific sites as shown in Figure 6a. The multimodel mean performs similarly well to the best‐performing model, HTESSEL, with a correlation of 0.69 and a root‐mean‐square error of 0.21 and a bias close to zero. Overall, the coupling strength is slightly overestimated in JULES and SURFEX when compared to the in situ data (bias of 0.11 and 0.12, respectively) and underestimated in GLEAM (bias of −0.11). The intermodel range at each respective site is generally larger than the estimate of the tower uncertainty represented by the shaded gray area in the figure.

Figure 6

Results of the comparison between model and tower (a) Χ and (b) Π. (top) Tower‐based coupling results as circles, with the multimodel mean for the 1979–2012 period in the background. (bottom) The magnitude of model‐ and tower‐based coupling at all sites, with sites ordered by decreasing coupling strength. The shaded gray area shows coupling results based on the measured tower latent heat flux and the latent heat flux calculated as the energy balance residual (see section 2.2). The darker circles (top) and the dashed line (bottom indicate the AU‐DaP (high coupling) and CH‐Fru (low coupling) sites that are further investigated in Figures 8 and 9.

Table 2

Pearson Correlation, Bias, and Root‐Mean‐Square Error Between Model and Tower Coupling

Model	r	Χ	RMSE	r	Π	RMSE
		Bias			Bias
HTESSEL	0.69	−0.01	0.21	0.67	−0.01	0.24
JULES	0.62	0.11	0.27	0.61	0.11	0.27
ORCHIDEE	0.55	0.04	0.28	0.50	0.02	0.32
SURFEX	0.62	0.12	0.26	0.59	0.10	0.30
GLEAM	0.58	−0.11	0.32	0.55	−0.08	0.30
Ensemble Mean	0.69	0.03	0.21	0.63	0.03	0.25

Figure 8

(a) Tower and (b) intermodel mean daily time series of net radiation, latent heat flux, and temperature at CH‐Fru (grassland, Switzerland). Only the warm season used to calculate the coupling metrics are shown. Top subpanels show actual values, bottom subpanels show anomaly time series.

Figure 9

Same as Figure 8 but for AU‐DaP (cleared savanna, Australia).

For П, the agreement between models and towers is similar to Χ (Figure 6b): Pearson correlations vary between 0.50 (ORCHIDEE) and 0.68 (HTESSEL) and the root‐mean‐square error varies between 0.24 (HTESSEL) and 0.32 (ORCHIDEE). With the highest Pearson correlation, lowest RMSE, and lowest absolute bias, HTESSEL is the best‐performing model for this metric. The performance of the model mean is slightly lower than the best‐performing model for this metric (r = 0.63, RMSE = 0.25, bias = 0.03), but still higher than the other four models. The model bias with respect to tower results is similar to Χ: overestimation by JULES (0.11) and SURFEX (0.10), and underestimation by GLEAM (−0.08). It is important to bear in mind that the order of sites in Figure 6a (bottom) is not exactly the same as the order in Figure 6b (bottom).

Even though it is encouraging that the models show reasonably good agreement with tower data, this comparison can only be considered to be a first and very general step toward validation of soil moisture‐temperature coupling in models. First, the model and tower footprints differ by several orders of magnitude; this means that the dominant climate controls at the tower site may not represent the most dominant controls at the scale of the model pixel, thus the two data sets would not show perfect agreement even under ideal circumstances. Second, the tower sites are not evenly distributed over land and climate types, with few sites located in coupling hot spots. Third, the off‐line simulations may not capture the full feedback loops occurring in nature and measured at the sites. Finally, as already mentioned above, the identified skill of each model in regard to its representation of soil moisture‐temperature coupling is conditioned on the suitability of the metrics to capture this interaction.

To address the issue of representativeness of the tower sites, together with the uncertainty in reanalysis forcing and in situ measurements, we calculated Pearson correlation coefficients between model and tower anomaly time series of air temperature, net radiation, and latent heat flux. The results are presented in Figure 7, where sites are sorted by descending soil moisture‐temperature coupling strength based on П, therefore, in the same order as Figure 6b. The propagation of errors through the models is evident from Figure 7. The level of agreement between model and tower temperature anomalies is generally quite high, with an average correlation of 0.89, but correlations decrease to approximately 0.47 and 0.32 for the anomalies in net radiation and the latent heat flux, respectively. The model spread also increases along this progression from forcing to simulated time series. Ideally, soil moisture would also be evaluated here. However, only a few sites provide root zone soil moisture data, and these are located in radiation‐limited climates where soil moisture‐temperature coupling is not relevant.

Figure 7

Correlations between model and tower anomaly time series of (a) air temperature, (b) net radiation, and (c) latent heat flux, as well as the absolute difference between model and tower coupling according to (d) Χ and (e) П at each of the tower sites. The sites have been sorted by increasing evaporative fraction.

One would expect the level of agreement between model and tower R n and λE to be linked to the agreement between the model and tower estimates of soil moisture‐temperature coupling (Χ, Π). This is supported by the results of HTESSEL, which showed the best agreement between model and tower coupling, and also shows high correlations between model and tower R n and λE; it also holds true for ORCHIDEE, which shows relatively low correlations between model and tower R n, λE, and coupling. Interestingly, however, GLEAM shows the best overall agreement with tower data in terms of λE in Figure 7, but had relatively poor performance in the soil moisture‐temperature coupling fields (Table 2). This shows that a realistic representation of the dynamics of latent heat flux and air temperature is not always accompanied by a realistic sensitivity of energy partitioning to air temperature and vice versa.

Figures 7d and 7e show the absolute difference between model and tower coupling. Per model, the errors tend to be higher at sites with lower agreement between tower and model R n and/or λE. For Χ, errors in λE and R n are approximately equally important, while errors in Π are more closely related to errors in R n than to errors in λE. Two unusual sites are ES‐Amo (site 27) and IT‐Noe (site 30), where correlations between model and tower temperature anomalies are substantially lower than average, and the error in model coupling is especially low; these towers are located in coastal areas and their results should be interpreted with caution. Also note the high similarity between the patterns of model coupling errors for both metrics (r = 0.6–0.8, depending on the model), which implies that the choice of one metric over the other does not alter the conclusions of the analysis.

We further examine errors in coupling by climate, land cover, and tower representativeness. The errors in coupling tend to be lower at towers located in arid/semiarid climates and higher in temperate and continental climates, respectively. Similarly, errors tend to be higher at forested sites than at sites covered by shorter vegetation types. However, these differences between climates and land cover types are generally not statistically significant. Finally, we evaluate the role tower representativeness here based on the portion of the 0.5° pixel that is represented by the tower land cover. Tower representativeness varies between 3 and 78% over the study sites, with a mean of 43%. Perhaps somewhat surprisingly, however, there is no clear relationship between tower representativeness and coupling error.

Finally, to further understand the behavior of models and coupling metrics over different climatic regimes, we analyzed daily time series of R n, λE, and T at two characteristic sites: (1) CH‐Fru, a grassland site in Switzerland (47.12°N, 8.54°E) where the values of the coupling metrics are low (Χ = −0.47, Π = −0.48, based on tower measurements) and (2) AU‐DaP, a cleared savanna site in Australia (14.06°S, 131.32°E) where coupling is strong (Χ = 0.45, Π = 0.64). The exact location of these towers is shown in Figure 6. At CH‐Fru, the latent heat flux accounts for a relatively large part of net radiation throughout the warm season (Figure 8). The anomaly time series show that the dynamics of both energy fluxes and temperature are very similar, as days with high/low net radiation correspond with high/low latent heat flux and temperatures. The energy available for heating (here R n – λE) can be inferred from Figure 8 and tends to increase with increasing net radiation. In terms of П, this means that R and R – λE are highly correlated, and therefore, the first and second terms of П will be of a similar magnitude (thereby the low values of this metric).

The time series at the site with strong soil moisture‐temperature coupling show a different pattern (Figure 9). At the beginning of the warm season, λE accounts for a small portion of R , but this increases toward the end of the season. The direction of the anomalies in λE and T tend to be opposite at the beginning of the warm season or in dry years, though this effect diminishes by the end of the season or in wetter years when λE accounts for a larger portion of R . This suggests that soil moisture availability plays an important role at this location. Since in the majority of the time series the direction of T anomalies is opposite to that of λE anomalies, coupling according to Χ is positive. Though R anomalies are positively correlated to T anomalies, the R n – λE anomalies are more correlated to T, leading to a positive value for Π as well.

Summary and Conclusions

The soil moisture‐temperature coupling strength in a set of land surface models was calculated based on two different correlation‐based metrics, focusing on the terrestrial segment of the coupling. The advantage of the metrics used in this study is that they are simple and can easily be applied to both observational and model data sets. The disadvantage is that they are simplifications of the actual physics and as such they cannot be used to infer actual causality. After a global evaluation of the hot spots of soil moisture‐temperature coupling, the model results were compared and discussed and evaluated against the eddy‐covariance tower data.

In agreement with previous studies, hot spots of soil moisture‐temperature coupling are found in transitional climates such as the Great Plains, southern Africa, India, and northern Australia. The location of these hot spots is similar across models and for both metrics, but at a regional scale the strength of the coupling can be highly variable. Of the models included in this study, soil moisture‐temperature coupling is comparatively stronger in JULES and SURFEX, and comparatively weaker in GLEAM. Despite the differences in representative area of the model and eddy‐covariance tower data, the errors in tower data and gridded forcing, and the limited number of sites available, there is relatively good agreement between models and observations (r = 0.5–0.7). The ensemble mean performs similar to the best‐performing model, suggesting that for studies of land‐atmospheric coupling an ensemble of models is recommended.

Since the models were run off‐line and based on common forcing, discrepancies between models, and between model and tower results, can be attributed to differences in the partitioning of incoming radiation and precipitation. Determining the cause of these differences is not straightforward, as they are a result of multiple factors, including model parameterizations and soil and vegetation characteristics. A set of additional experiments which isolate the effects of model structure and parameterization is necessary to quantify the importance of these factors and improve the representation of soil moisture‐temperature coupling in these models. Nevertheless, our results show that a realistic portrayal of the partitioning of incoming radiation does not always coincide with a realistic portrayal of the sensitivity of air temperature to this partitioning. This distinction is important, because accurate representation of soil moisture‐temperature coupling strength is crucial to predict future climate variability. In turn, accurate portrayal of climate variability is necessary to gain insight into how the occurrence and intensity of temperature extremes and droughts will change in the future, especially since the transitional climates where soil moisture‐temperature coupling occurs also tend to be regions that have been experiencing negative trends in water availability over the last decades (Mekonnen & Hoekstra, 2016; Wada et al., 2011).