Articles | Volume 25, issue 1
Research article
28 Jan 2021
Research article |  | 28 Jan 2021

Evaluating a land surface model at a water-limited site: implications for land surface contributions to droughts and heatwaves

Mengyuan Mu, Martin G. De Kauwe, Anna M. Ukkola, Andy J. Pitman, Teresa E. Gimeno, Belinda E. Medlyn, Dani Or, Jinyan Yang, and David S. Ellsworth

Land surface models underpin coupled climate model projections of droughts and heatwaves. However, the lack of simultaneous observations of individual components of evapotranspiration, concurrent with root-zone soil moisture, has limited previous model evaluations. Here, we use a comprehensive set of observations from a water-limited site in southeastern Australia including both evapotranspiration and soil moisture to a depth of 4.5 m to evaluate the Community Atmosphere-Biosphere Land Exchange (CABLE) land surface model. We demonstrate that alternative process representations within CABLE had the capacity to improve simulated evapotranspiration, but not necessarily soil moisture dynamics–highlighting problems of model evaluations against water fluxes alone. Our best simulation was achieved by resolving a soil evaporation bias, using a more realistic initialisation of the groundwater aquifer state and higher vertical soil resolution informed by observed soil properties, and further calibrating soil hydraulic conductivity. Despite these improvements, the role of the empirical soil moisture stress function in influencing the simulated water fluxes remained important: using a site-calibrated function reduced the soil water stress on plants by 36 % during drought and 23 % at other times. These changes in CABLE not only improve the seasonal cycle of evapotranspiration but also affect the latent and sensible heat fluxes during droughts and heatwaves. The range of parameterisations tested led to differences of ∼150W m−2 in the simulated latent heat flux during a heatwave, implying a strong impact of parameterisations on the capacity for evaporative cooling and feedbacks to the boundary layer (when coupled). Overall, our results highlight the opportunity to advance the capability of land surface models to capture water cycle processes, particularly during meteorological extremes, when sufficient observations of both evapotranspiration fluxes and soil moisture profiles are available.

1 Introduction

Droughts and heatwaves can have severe and long-lasting impacts on terrestrial ecosystems (Allen et al., 2015; Reichstein et al., 2013) and humans (Matthews et al., 2017; Pal and Eltahir, 2016). Global climate models are commonly used to project how anthropogenic climate change will affect the magnitude, frequency and intensity of droughts and heatwaves. Heatwaves are projected to increase in the future in response to climate change (Dosio et al., 2018; Zhao and Dai, 2017). The future of droughts is less clear: projections of an increase in future droughts are common in the literature (Ault, 2020), yet regional precipitation projections remain uncertain (Collins et al., 2013) and land surface processes relevant to drought are poorly represented in climate models (Ukkola et al., 2018a).

While there is no universal definition, drought can be classified into meteorological, agricultural, hydrological and socioeconomic drought. From a climate model perspective, drought is an anomalous lack of water at the land–atmosphere interface sustained over time. It begins with a reduction in precipitation (“meteorological” drought), and if this persists it can evolve into “agricultural” drought via low soil moisture or into “hydrological” drought through low streamflow or groundwater (Orth and Destouni, 2018). A critical feedback exists between low soil moisture availability and heatwaves (Seneviratne et al., 2010; Teuling et al., 2010; Vogel et al., 2017). As soil moisture becomes depleted, the surface energy partitioning becomes increasingly dominated by sensible heat fluxes (QH) relative to latent heat fluxes (QE). This can lead to a positive feedback whereby the high sensible heat fluxes warm the boundary layer, which, combined with the reduced evaporation, leads to increased atmospheric demand for moisture, exacerbating land desiccation (Miralles et al., 2019). A combination of drought and heatwaves leads to wide-ranging impacts on the functioning of terrestrial ecosystems (Reichstein et al., 2013; Schumacher et al., 2019). For example, during the European heatwave and drought in 2003, terrestrial carbon losses of up to 0.5 Pg C were reported, corresponding to roughly 4 years of European terrestrial net carbon uptake (Ciais et al., 2005).

Given projections of worsening heatwaves and potentially more droughts under future climate change, the importance of land surface models (LSMs) to capture land responses and feedbacks to the atmosphere during climate extremes is becoming increasingly recognised (Mazdiyasni and AghaKouchak, 2015; Schumacher et al., 2019; Yang et al., 2019). Despite many improvements to LSMs over the past decades, LSMs have remained poor at simulating water fluxes during water-stressed periods (Egea et al., 2011; De Kauwe et al., 2017; Powell et al., 2013; Trugman et al., 2018; Ukkola et al., 2016a), which likely contributes to biases in land–atmosphere feedbacks during heatwaves (Sippel et al., 2017). LSMs commonly underestimate interannual variations in terrestrial water storage (Humphrey et al., 2018), underestimate QE during droughts (Powell et al., 2013; Ukkola et al., 2016a) and lack “persistence” by responding too strongly to short-term precipitation variation (Tallaksen and Stahl, 2014). Poor representation of hydrological processes has been identified as a key reason for model biases. There is uncertainty around soil moisture dynamics, how soil texture information is translated to soil hydraulic properties through pedotransfer functions and how water fluxes are partitioned to different components of evapotranspiration and runoff (Clark et al., 2015; Lian et al., 2018; Van Looy et al., 2017). Various approaches have been adopted to improve LSM hydrology, such as the introduction of groundwater dynamics (Niu et al., 2007), alternative pedotransfer functions (Best et al., 2011) and subgrid-scale processes for runoff generation (Decker, 2015). By contrast, the functions used in LSMs to represent the effect of declining water availability on vegetation function are poorly constrained by data (Medlyn et al., 2016) and not consistently applied. Specifically, some models downregulate the maximum rate of Rubisco carboxylation, whilst others reduce stomatal parameters (De Kauwe et al., 2013). Models also do not account for differences in species-level sensitivity to drought (De Kauwe et al., 2015; Klein, 2014; Zhou et al., 2014). This model gap has driven a significant investment in new theoretical approaches (Dewar et al., 2018; Sperry et al., 2017; Wolf et al., 2016).

Despite model developments, it has remained difficult to disentangle the reasons behind poor model performance due to a lack of suitable observations. Root-zone soil moisture estimates are rare, and, whilst satellite estimates are available, they only cover the top few centimetres or are only available at coarse spatial resolution. Meanwhile, QE is routinely measured at the site scale, but gridded large-scale estimates remain highly uncertain (Pan et al., 2020). As such, many past model evaluations have focused on observed QH and QE from eddy covariance observations (Best et al., 2015) or near-surface soil moisture and evaporation from water balance sites (e.g. Schlosser et al., 2000). However, it is rare for evaluations of LSMs, designed for use in climate models, to utilise observations of soil moisture extending through the root zone with concurrent measurements of water fluxes at high temporal frequency. In this paper, we use a novel dataset from the water-limited Eucalyptus Free-Air CO2 Enrichment (EucFACE) experiment site in southeastern Australia to evaluate the Community Atmosphere-Biosphere Land Exchange (CABLE) LSM. At this site, frequent measurements of each component of the water balance were made coincident with soil moisture observations to a depth of 4.5 m. The highly variable rainfall at this site leads to extended dry-downs, and the heatwaves in summer commonly exceed 35 C. We use this high-quality dataset to assess multiple model assumptions commonly used across LSMs within a single model framework, evaluating both simulated fluxes and state variables at seasonal to annual scales and across weather (heatwaves) and climate (drought) phenomena.

Figure 1(a) Location of the experimental site in western Sydney, Australia (333659′′ S, 1504417′′ E) shown by the red star. (b) Distribution of six rings (© Google Maps, 2020, EucFACE experiment site, 1:50, available at:,150.7379194,356m/data=!3m1!1e3, last access: 2 July 2020). (c) Understorey vegetation and infrastructure inside a ring (photograph taken by Mengyuan Mu). (d) Canopy structure and central tower (photograph taken by Mengyuan Mu).

2 Methods and data

2.1 Site information

The EucFACE experiment is located on an ancient alluvial floodplain, 3.6 km from the Hawkesbury River in western Sydney, Australia (333659′′ S, 1504417′′ E) (Gimeno et al., 2018a; Fig. 1). The site has a temperate–subtropical transitional climate with a mean annual temperature of 17.8 C and the mean annual precipitation of 719.1 mm evenly distributed over the year. EucFACE is a water-limited site experiencing frequent droughts and low water availability. The site is in an open woodland with a canopy height of 18–23 m and a plant area index (including leaf and woody components) that varied between 1.3 and 2.2 m2 m−2 (mean =1.7m2 m−2) over the study period. The overstorey is dominated by a single species, Eucalyptus tereticornis Sm., with scattered individuals of Eucalyptus amplifolia Naudin. The upper soil layer is a loamy sand with a sand fraction >75 %; at 30–80 cm depth, there is a higher clay content layer (15–35 % clay), and below the clay layer sand–clay loam soil extends to a depth of 300 cm. Between 300–350 cm and 450 cm depth, the soil is >40 % clay (Gimeno et al., 2016). The observed water table is at ∼12m. The site is characterised as nutrient poor, especially lacking in available phosphorus (Crous et al., 2015; Ellsworth et al., 2017). In this paper we evaluate CABLE against the averaged data from rings 2, 3 and 6, which are exposed to the ambient atmospheric CO2.

2.2 Observation data

In our study, CABLE is driven by in situ meteorological data and observed leaf area index (LAI) from 2013 to 2019. The photosynthetically active radiation (PAR; LI-190, LI-COR, Inc., Lincoln, NE, USA), air temperature and relative humidity (HUMICAP® HMP 155, Vaisala, Vantaa, Finland) were measured every second, and 1 min averages were recorded on data loggers (CR3000, Campbell Scientific Australia, Townsville, Australia). Following Yang et al. (2020), the meteorological data were gap-filled (0.8 % of values, from 1 January 2013 to 31 December 2019) using linear interpolation, aggregated to 30 min averages and subsequently used to force CABLE at the 30 min resolution. LAI was calculated from the measurements of above- and below-canopy PAR at each ring following Duursma et al. (2016). Since the site LAI represents the plant area index (including both woody part and leaves), to reflect the actual leaves condition we follow Yang et al. (2020) and reduce the LAI by a constant branch and stem cover (0.8 m2 m−2) estimated by the lowest LAI when the canopy shed almost all leaves during November 2013. The CO2 concentration was measured every 5 min at each ring and then gap-filled and aggregated to 30 min averages.

To evaluate CABLE, we used measurements of transpiration (Etr), soil evaporation (Es) and volumetric water content (θ) at different soil depths (see below). Etr estimates are derived from tree sap flow velocities (three or four trees per experimental ring) using the heat pulse compensation technique (Gimeno et al., 2018a). Sap flow velocity is translated to Etr by multiplying the sapwood area estimated from basal area inside each ring and a correlation between sapwood and basal areas based on 35 trees adjacent to the experimental rings. Es is computed from the soil moisture change in the top 5 cm depth monitored at two locations in each of the three ambient rings. The Es data also include transpiration from the dynamic (flushes and wilts) understorey vegetation (Collins et al., 2018; Pathare et al., 2017). For Es, Gimeno et al. (2018a) excluded rainy days and days preceded by a day with >2mm d−1 of precipitation. To represent variability in Etr and Es across rings, we show the mean and the uncertainty within ring estimates in all figures.

We used two sets of observations for θ to evaluate CABLE's simulated soil hydrology. The first dataset is from neutron probe measurements monitored at two locations in each ring every 10 to 21 d (lower frequency in 2017), covering the period January 2013 to July 2019. These data are collected at 12 different depths: 25 cm intervals from 25 to 150 cm depth and 50 cm intervals from 150 to 450 cm depth. The neutron probe counts are converted to θ via the site-specific linear correlation between the raw reading of neutron probe and the lab-measured soil θ sampled at the same depth as probes (Gimeno et al. 2018a). The second dataset is daily derived measurements from frequency-domain reflectometers (CS650, Campbell Scientific Australia, Garbutt, Qld) at each ring, monitoring to a depth of 25 cm and covering the period January 2013 to December 2019.

2.3 Model description

CABLE is a LSM that can be used in stand-alone mode with prescribed meteorological forcing (Haverd et al., 2013; Ukkola et al., 2016b; Yang et al., 2020) or coupled to the Australian Community Climate and Earth System Simulator (ACCESS; Bi et al., 2013; Law et al., 2017) or the Weather and Research Forecasting (WRF) model (Decker et al., 2017; Hirsch et al., 2019b) to provide energy, water and momentum fluxes to the lower atmosphere. The standard version of CABLE has been widely evaluated (De Kauwe et al., 2015; Li et al., 2012; Lorenz et al., 2014; Ukkola et al., 2016b; Wang et al., 2011; Williams et al., 2009), and the model's overall performance in simulating energy, water and energy fluxes is in line with other LSMs (Best et al., 2015). A detailed description of model components can be found in Kowalczyk et al. (2006) and Wang et al. (2011). The version of CABLE used here includes multiple process updates (Decker, 2015; Decker et al., 2017; Kala et al., 2015).

2.3.1 Hydrology scheme

We use the hydrology scheme from Decker (2015) that includes an improved representation of subsurface hydrology similar to that implemented in the Community Land Model (CLM; Lawrence and Chase, 2007; Oleson et al., 2008). Saturation- and infiltration-excess runoff generation mechanisms are represented, and a dynamic groundwater component with aquifer water storage is included. CABLE uses six soil layers covering a depth to 4.6 m and allows for vertical heterogeneity in soil parameters. The scheme solves the vertical redistribution of soil water following the modified Richards equation (Decker and Zeng, 2009):

(1) θ t = - z K z Ψ - Ψ E - F soil ,

where θ is the volumetric water content of the soil (mm3 mm−3), K (mm s−1) is the hydraulic conductivity, Ψ (mm) is the soil matric potential, ΨE (mm) is the equilibrium soil matric potential, z (mm) is soil depth and Fsoil (mmmm-1s-1) is the sum of subsurface runoff and Etr (Decker, 2015). A 22.8 m deep unconfined aquifer is simulated below the six-layer soil column by incorporating a simple water balance model:

(2) d W aq d t = q re - q aq,sub ,

where Waq (mm) is the mass of water in the aquifer; qaq,sub (mm s−1) the subsurface runoff removed from the aquifer; and qre (mm s−1) the water flux between the aquifer and the bottom soil layer, computed by the modified Darcy's law as

(3) q re = K aq Ψ aq - Ψ n - Ψ E , aq Ψ E , n z wtd - z n ,

where Kaq (mm s−1) is the hydraulic conductivity within the aquifer, Ψaq and ΨE,aq (mm) are the soil matric potential and the equilibrium soil matric potential for the aquifer, and Ψn and ΨE,n (mm) are the soil matric potential and the equilibrium soil matric potential for the bottom soil layer. zwtd and zn (mm) are the depth of the water table and the lowest soil layer, respectively. The groundwater aquifer is assumed to sit above an impermeable layer of rock, giving a bottom boundary condition of

(4) q out = 0 .

Subsurface runoff (qsub, mm s−1) is calculated from

(5) q sub = sin d z d l q ^ sub e - z wtd f p ,

where dzdl is the mean subgrid-scale slope, q^sub (mm s−1) is the maximum rate of subsurface drainage assumed to be achieved when the whole soil column is saturated and fp is a tunable parameter. qsub is generated within the aquifer and for each saturated soil layer below the third soil layer.

2.3.2 Soil evaporation (Es)

The computation of Es (kgm-2s-1) considers the subgrid-scale soil moisture heterogeneity within a grid square (Decker, 2015) and is given as

(6) E s = F sat E s + ( 1 - F sat ) β s E s ,

where Fsat is the saturated fraction of a grid cell, Es (kgm-2s-1) is the potential evaporation without soil moisture stress and βs is an empirical soil moisture stress factor (see below) that limits evaporation as water becomes limiting in the top soil layer (Sakaguchi and Zeng, 2009). Es is given by

(7) E s = ρ a ( q sat ( T srf ) - q a ) r g ,

where ρa (kg m−3) is the air density, qsat(Tsrf) (kg kg−1) is the saturated specific humidity at the surface temperature, qa (kg kg−1) is the specific humidity of the air and rg (s m−1) is the aerodynamic resistance term.

βs is computed as

(8) β s = 0.25 1 - cos π θ unsat θ fc 2 ,

where θunsat (mm3 mm−3) is the volumetric water content in the unsaturated portion of the top soil layer (top 2 cm), and θfc (mm3 mm−3) is the field capacity in the top soil layer.

2.3.3 Transpiration (Etr)

CABLE's canopy is represented using a two-leaf model, which computes photosynthesis, stomatal conductance, Etr (kgm-2s-1) and leaf temperature separately for sunlit and shaded leaves. Etr (for each sunlit and shaded leaf) is calculated following the Penman–Monteith equation:

(9) E tr = Δ R n + C p M a D l g h + g r λ Δ + γ g h + g r g w ,

where λ (J kg−1) is the latent heat of vaporisation, Dl (Pa) is the vapour pressure deficit at the leaf surface, Cp (Jkg-1K-1) is the air heat capacity, Ma (kg mol−1) is the molar mass of air, Δ (Pa K−1) is the slope of the curve relating saturation vapour pressure to air temperature and γ (Pa K−1) is the psychrometric constant. gh, gr and gw (molm-2s-1) are the conductances for heat, radiation and water, respectively. Rn (W m−2) is the non-isothermal net radiation, calculated as

(10) R n = R n - C p M a ( T a - T l ) g r ,

where Rn (W m−2) is the net radiation under isothermal conditions, and Ta and Tl is the air and leaf temperature (K), respectively.

gw is calculated as

(11) g w - 1 = g a - 1 + g b - 1 + g s - 1 ,

where ga (molm-2s-1) is canopy aerodynamic conductance, and gb (molm-2s-1) is leaf boundary layer conductance for free and forced convection (Kowalczyk et al., 2006). gs (molm-2s-1) is the leaf stomatal conductance following Medlyn et al. (2011):

(12) g s = g 0 + 1.6 1 + g 1 β D l A C s ,

where A (µmolm-2s-1) is the photosynthetic rate; Cs (µmol mol−1) is the CO2 concentration at the leaf surface; β (unitless) is the soil moisture stress factor on plants; and g0 (molm-2s-1) and g1 (kPa0.5) are fitted parameters representing the residual stomatal conductance when A=0 and the sensitivity of conductance to the assimilation rate, respectively. g1 reflects the plant's water use strategy and was derived for each plant functional type in CABLE (De Kauwe et al., 2015) based on a global synthesis of stomatal behaviour (Lin et al., 2015). β is calculated as

(13) β = i = 1 n f root , i θ i - θ w , i θ f c , i - θ w , i ,

where θi, θfc,i and θw,i (mm3 mm−3) are the soil moisture content, the field capacity and the wilting point for soil layer i, respectively, and froot,i is the root mass fraction of soil layer i.

CABLE does not have the capacity to simulate interacting water fluxes between the understorey and overstorey vegetation. Instead, it uses a “tiling” approach (fractionally weights separate simulations). As a result, comparisons between CABLE's Es and data-derived Es during wetter periods would be expected to be an underestimate as we only consider the fluxes from the overstorey trees. To quantify the effect of the understorey transpiration on the water balance, we also ran an extra simulation for the grass understorey at this site with the same setting as Watr (see below) but using CABLE default grass physiology parameters and a fixed LAI (1 m2 m−2 – site average). The estimated multi-year mean transpiration of 0.94 mm d−1 can be regarded as an upper estimate since the simulation does not consider grass dynamics, overstorey rainfall interception, or water and energy competition between tree and grass. Not accounting for understorey transpiration will lead to an overestimate of moisture availability in the soil profile.

2.4 Experiment design

We conducted a series of model experiments based on weaknesses identified in previous LSM evaluation studies. In our experiments, we deliberately adopted a “layering” approach: sequentially resolving a key systematic model bias and then layering additional experiments to examine how much additional benefit each experiment added to model performance. We choose to first resolve a soil evaporation bias as it affects total evapotranspiration (ET) partitioning; however, its importance is limited to the top soil layers, particularly during the period following rain. We then modified the initial water table depth as this fundamentally affects the root-zone soil moisture state. Next, we explored assumptions related to soil column discretisation and parameters, and further optimised key hydraulic parameters to improve overall soil moisture biases. We chose to explore parameter assumptions at this point in the experimental set-up as the previous experiments aimed to resolve existing biases that affected the overall soil moisture availability. Finally, we explored alternative soil moisture stress functions as the last step because this factor integrates, and arises from, the soil moisture state. The experimental order allowed us to probe model biases in a systematic way, but it is important to note that there is no perfect experimental order and alternative permutations would lead to subtly different interpretation of results. In fact, this is what commonly happens in model evaluation that explore a single factor (e.g. the soil moisture stress factor). Each experiment is described in detail below, and a summary of all experiments is provided in Table 1.

Table 1The experiments conducted. Layers refers to the number of soil layers. Increase resistance refers to increasing surface resistance to soil evaporation. Soil heterogeneity indicates whether soil properties and hydraulic parameters change with depth. The adjustment of θw, θsat and Ksat and the method used to calculate β are the final two columns.

Download Print Version | Download XLSX

In all experiments, LAI and physiology parameters were prescribed based on site observations (Table S1 in the Supplement). We tested the difference of using the CABLE default evergreen broadleaf physiology parameters (Fig. S1 in the Supplement) compared to using the site physiology (Fig. 2) and found that using site parameters increases Etr (due to higher g1 and increased sensitivity of carbon fixation to temperature), in turn reducing Es and θ.

Figure 2Control simulation (Ctl). (a) Etr,Es and precipitation (P) between 2013 and 2015. The shaded areas represent uncertainty between three ambient rings. Both simulations and observations are smoothed with a 3 d window to aid visualisation. (b) θ in the top 0.25 m from 2013 to 2019. (c) The vertical distribution of θ measured at observed dates from 2013 to 2019. (d) The vertical distribution of θ in Ctl for observed dates from 2013 to 2019. (e) θ differences between CABLE and observations. Note the different time axis for panels (c–e) relative to panels (a, b) due to different sampling intervals for soil moisture and fluxes.


All experiments were spun up using an iterative process, recycling all years of the meteorological forcing until the change between two iterations was <0.001m3 m−3 for soil moisture, <0.01C for soil temperature and <0.0001m3 m−3 for aquifer moisture.

2.4.1 Control experiment

The control simulation (Ctl) uses the default version of CABLE with six soil layers (but with site-specific physiology and LAI). The soil hydraulic parameters are derived via the pedotransfer functions based on Cosby et al. (1984) using the global soil texture map from the Harmonized World Soil Database (Fischer et al., 2008). Soil parameters are the same throughout the six-layer soil column.

2.4.2 Increasing the resistance for soil evaporation (Sres)

Previous studies suggest LSMs vary widely in their simulation of Es. For example, De Kauwe et al. (2017) found that, in an ensemble of 10 models, 6 models simulated ∼2–3.5 times more Es than the other 4 models. LSMs also partition evapotranspiration between Etr and Es with a high degree of uncertainty (Lian et al., 2018). At many sites, high springtime evapotranspiration can be linked to excessive Es rather than Etr (Decker et al., 2017; Ukkola et al., 2016b) and can lead to biases in soil moisture availability later in the growing season.

We note that models have attempted to resolve this Es bias through different mechanisms, for example, via a litter layer (Haverd and Cuntz, 2010; Sakaguchi and Zeng, 2009) or by limiting Es via adding the resistances to vapour diffusion through the soil pores and the surface viscous sublayer (Decker et al., 2017; Haghighi and Or, 2015; Swenson and Lawrence, 2014). Here, we adopt a simple litter layer (Decker et al., 2017) which adds an additional surface resistance to vapour and heat fluxes but does not limit rainfall infiltration. After adding the additional resistance, Es is calculated as

(14) E s = ρ a ( q sat ( T srf ) - q a ) r g + r lit ,

where rlit is the resistance (s m−1) for diffusion via the litter layer of depth zl (m) (default value is 10 cm), given by

(15) r lit = z l d ,

where d is the diffusivity of water vapour in air (m2 s−1).

2.4.3 Water table initialisation experiment (Watr)

The parameters governing the groundwater aquifer saturation and water table depth are both highly uncertain and difficult to constrain from observations. We investigated the importance of a correct water table depth to the simulation soil moisture and water fluxes. To better match the observed water table depth at EucFACE, we changed the aquifer θsat from the model default value (0.235 m3 m−3) to θsat set based on the observed soil texture at 4.5 m depth (0.448 m3 m−3), which reduced the initial saturation of the aquifer from 100 to 52 %. This has the effect of lowering the water table to ∼12m, in line with observations (Gimeno et al. 2018a).

2.4.4 High-resolution soil experiment (Hi-Res)

Most LSMs assume that soil parameters are depth invariant through the soil profile. The number of layers typically ranges from a minimum of 2 through to 6 in CABLE and up to 20 in the Community Land Model (Lawrence et al., 2019). Here, we test the impact of increasing the number of discrete soil layers, informed by observations of the varying vertical soil texture at the EucFACE site. Recent soil maps (e.g. SoilGrids; Hengl et al., 2017) have begun to capture vertical variations in soil texture, so it is important to test the impact in LSMs.

We performed two sub-experiments in Hi-Res:

  1. the number of vertical soil layers was increased from 6 to 31 (to match the resolution of observed soil texture, which was sampled at 15 cm intervals) (Hi-Res-1);

  2. soil parameters were allowed to vary vertically based on observed soil texture (Hi-Res-2).

To implement vertically varying soil parameters, the observed fractions of sand, clay and silt; soil bulk density; and organic carbon fraction were taken from measurements at each ambient CO2 ring and interpolated into 31 layers using the ∼15 cm resolution of the observations. Soil hydraulic parameters are computed using the same pedotransfer functions as used in Ctl but allowed to vary with depth based on the vertical heterogeneity in soil properties. Since CABLE assumes that the Clapp and Hornberger parameter and the parameter for suction at saturation are identical for the aquifer and the bottom soil layer, adding depth-varying soil parameters in Hi-Res-2 also changes these two parameters for the aquifer.

2.4.5 Soil parameter optimisation experiment (Opt)

As it is impractical to measure soil hydraulic parameters at the global scale, pedotransfer functions are used to convert widely measured soil properties into global soil hydraulic parameter datasets (Dai et al., 2013; Kishné et al., 2017). However, most of the widely used pedotransfer functions are empirical equations derived from the limited experimental samples measured for the specific locations (Cosby et al., 1984; van Genuchten, 1980). The adaptability of these pedotransfer functions is always confined by their underrepresentation of some soil properties, such as soil aggregate stability or macroporosity (Puhlmann and von Wilpert, 2012), and can lead to a divergence in model parameters (Van Looy et al., 2017; Zhang and Schaap, 2019). As a result, parameter calibrations are commonly used to obtain more accurate representations.

First, we used the site observations to adjust the plant wilting point (θw) and volumetric water content at saturation (θsat). With each layer as θw is changed, the corresponding residual water content (θres) was also updated to ensure it was smaller than θw. θsat was set to the observed maximum from the daily data measured by frequency-domain reflectometers for the top 30 cm. θsat was not adjusted below 30 cm as the observed maximum θsat is unlikely to represent saturated conditions due to lower soil moisture variability at depth. θw and θres were adjusted for each 15 cm layer in the soil column using the observed minimum (OBSmin) in each layer. When OBSmin was below the default θres, θres was set to OBSmin and θw to OBSmin+0.0001m3 m−3. When θres <OBSmin <θw, θw was set to OBSmin. Otherwise θres and θw were not adjusted.

Second, we optimised Ksat to test whether allowing the soil column to drain faster or slower reduced model biases in the soil moisture profile. Ksat was optimised by concurrently minimising the root mean squared error (RMSE) between modelled and observed transpiration, soil evaporation, and soil moisture over the total column and over the top 0.25 m.

2.4.6 Soil water limitation on transpiration (β-hvrd and β-exp)

LSMs use different empirical functional forms to represent the effect of water stress on vegetation function (see Introduction). To explore the influence of different functional formulations, we compare CABLE's default function (Eq. 13) to two alternative parameterisations: (1) an alternative hypothesis that plants optimise their root water uptake to exploit resources, with the wettest soil layer determining soil water stress on plants (β-hvrd; Haverd et al., 2016) and (2) a site-calibrated function to observations at EucFACE over the top 1.5 m (β-exp; Yang et al., 2020). We note that a number of studies have tested different water stress formulations (e.g. Egea et al., 2011), but this process evaluation is often decoupled from analysis of other contributing errors (e.g. LAI and/or soil hydrology).

The β-hvrd method tends to predict less water stress than the default function (Eq. 13) in CABLE when the moisture is unevenly distributed within the soil column. This function takes the form

(16) β = max ( α i δ i , i = 1 , n ) ,


(17) α = θ - θ w θ s γ / ( θ - θ w ) , ( θ - θ w ) > 0 0 , ( θ - θ w ) 0 ,

where αi is proportional to the root “shut-down” function (Lai and Katul, 2000) in the ith soil layer, and δi=1 if there are roots at the ith soil layer; otherwise δi=0. n is the total number of soil layers.

In β-exp, β is an exponential function calibrated to the site observations. Yang et al. (2020) fitted a non-linear relationship between β and θ, based on a fitted exponent term q (0.425, Table S1) using measured soil moisture over the top 1.5 m from EucFACE:

(18) β = i = 1 n f root , i θ i - θ w , i θ fc , i - θ w , i q .

2.4.7 Evaluation metrics

We used five metrics to evaluate CABLE's performance compared to observations. RMSE and mean bias error (MBE) were used to evaluate overall performance, and Pearson's correlation coefficient (r) the temporal variability. The absolute differences in modelled and observed 5th (P5) and 95th (P95) percentile values were used to evaluate the lower and upper tails, respectively. As the observed data have gaps, the metrics were only calculated for days for which observations were available.

3 Results

3.1 Control experiment (Ctl)

We first evaluate the Ctl simulation by comparing to the observed Etr, Es and soil moisture (Fig. 2). Overall, CABLE simulates Etr similarly to the observed Etr (r=0.85, RMSE = 0.34 mm d−1; Table 2) but overestimates peak Etr, which is particularly evident in the austral summer of 2014, by 0.54 mm d−1 on average (P95 in Table 2). However, it is worth noting that during the summer of 2014 there was an outbreak of psyllids leading to canopy defoliation (Gherlenda et al., 2016), which may explain part of the model–data mismatch (CABLE only accounts for canopy defoliation via a decline in LAI but not other damage, e.g. to the phloem). Compared to Etr, CABLE simulates Es less well (r=0.65, RMSE = 0.70 mm d−1; Table 2, Fig. 2a). Whilst the observations exclude rainy days, when CABLE reaches its highest Es, CABLE systematically overestimates mean and peak Es during observed days by 0.12 and 1.22 mm d−1, respectively (MBE and P95 in Table 2). Figure 2b shows that CABLE has a significant wet bias in the top 0.25 m soil moisture and never falls to the observed values below 0.08 m3 m−3 during drier periods. Given the excessive Es (Fig. 2a), the failure of the top 25 cm to dry out is surprising and suggests a parameterisation error and/or the impact of not accounting for understorey transpiration (see Sect. 2). Figure 2e shows that the wet bias in soil moisture is systematic, extending throughout the soil column (particularly between 2.5 and 4.5 m).

Table 2Performance metrics for the different experiments. Bold numbers are the best value among these experiments.

Download Print Version | Download XLSX

Taken together, the evaluation of the Ctl simulation implies that a good simulation in one evaporative flux (Fig. 2a) can be achieved for the wrong physical reasons and is associated with major systematic biases in the simulation of near-surface and root-zone soil moisture (Fig. 2b–d).

3.2 Increasing the resistance to soil evaporation experiment (Sres)

Implementing a litter layer (a proxy for additional surface resistance to Es) in CABLE significantly reduces Es from 305 mm yr−1 in Ctl to 204 mm yr−1 in Sres (Fig. 3a, Table 3). The simulation of peak Es is significantly improved compared to Ctl, but CABLE still overestimated Es (MBE and P95 in Table 2); this is particularly evident during an observed dry period in late 2013. As a consequence of lower Es compared to Ctl, Etr is markedly increased (from 341 mm yr−1 in Ctl to 402 mm yr−1 in Sres, Table 3), which implies a reduction in soil moisture stress in the profile (lower β). This degrades the simulated Etr relative to the observations for all metrics, particularly from around October 2013 to March 2014 (Fig. 3b). With an overall reduction in evapotranspiration, CABLE displays a considerably worse soil moisture profile (cf. Figs. 3c and 2d) and a larger wet bias through most of the soil profile (cf. Figs. 3d and 2e). Thus, resolving the Es bias alone relocated the bias to other model components, where it was less easily identified using commonly available measurements.

Table 3Average values from each experiment. Precipitation (P), total evapotranspiration (ET), transpiration (Etr), soil evaporation (Es), canopy evaporation (Ec), total runoff (R) including surface and subsurface runoff, soil water drainage to aquifer (Dr), gross primary production (GPP), latent heat (QE), sensible heat (QH), and volumetric water content in the 4.6m soil column (θ).

Download Print Version | Download XLSX

Figure 3Increasing-soil-evaporation-resistance experiment (Sres). (a) Es between 2013 and 2015. (b) Etr between 2013 and 2015. In panel (a) and (b) the shaded areas represent uncertainty between three ambient rings, and both simulations and observations are smoothed with a 3 d window to aid visualisation. (c) The vertical distribution of θ in Sres at observed dates from 2013 to 2019. (d) θ difference between CABLE and observations. Note the different time axis for panels (c, d) relative to panels (a, b) due to different sampling intervals for soil moisture and fluxes.


3.3 Water table (Watr) and vertical soil structure (Hi-Res) experiments

Figure 4 shows that reconciling the parameterisation of the aquifer θsat with the bottom layer θsat based on observed soil properties (Watr) leads to a marked improvement in the simulated soil moisture profile. By increasing the point of saturation and initialising the aquifer to be drier relative to saturation, CABLE simulates a more negative water potential in the aquifer, which promotes vertical drainage and results in a realistic water table depth in line with observations (simulated and observed ∼12 m over 2013/14). The wet bias in the top 3 m is markedly reduced (cf. Figs. 4d and 2e); however, the model now has a clear dry bias between 3 and 4.6 m. The simulated moisture in the top 0.25 m (Fig. 4b) is now also in better agreement with the observations (0.06 m3 m−3 in Watr vs. 0.11 m3 m−3 in Sres, MBE in Table S2). Finally, the bias in both the simulated Es and Etr is reduced by >0.2mm d−1 (MBE in Table 2), particularly evident during the summer of 2014.

Figure 4Water table initialisation experiment (Watr). (a) Etr and Es between 2013 and 2015. The shaded areas represent uncertainty between three ambient rings. Both simulations and observations are smoothed with a 3 d window to aid visualisation. (b) θ in the top 0.25 m from 2013 to 2019. (c) The vertical distribution of θ in Watr at observed dates from 2013 to 2019. (d) θ difference between CABLE and observations. Note the different time axis for panels (c, d) relative to panels (a, b) due to different sampling intervals for soil moisture and fluxes.


Increasing the number of soil layers from 6 to 31 (Hi-Res-1; Fig. S2) leads to a small improvement in the simulated temporal correlation (0.78 in Watr vs. 0.83 in Hi-Res-1; Table 2) of soil moisture, without notably changing the fluxes. The higher vertical resolution in the soil enables the transition of the dry-down to be better captured, in contrast to the alternating wet and dry patterns associated with the coarse vertical resolution at depths between 0.5 and 3.0 m depth in Watr (cf. Figs. S2c and 4c).

Allowing the soil parameters to vary vertically based on observed soil texture (Hi-Res-2; Fig. 5) reduces the dry bias in the lower layers in Watr (Fig. 4) but leads to a greater wet bias throughout the upper soil profile (<3m). The error in soil moisture has reduced in the mean, low and high extremes compared to Ctl and Sres (MBE, P5 and P95 in Table 2). Overall, Fig. 5 highlights a simulation with CABLE where the fluxes of Etr, Es and soil moisture are all in reasonable agreement with the observations (Table 3), albeit with an overestimation of peak Etr.

Figure 5High-soil-resolution experiment (Hi-Res-2), which uses 31 soil layers with depth-varying hydraulic parameters informed by observed soil properties. (a) Etr and Es between 2013 and 2015. The shaded areas represent uncertainty between three ambient rings. Both simulations and observations are smoothed with a 3 d window to aid visualisation. (b) θ in the top 0.25 m from 2013 to 2019. (c) The vertical distribution of θ in Hi-Res-2 at observed dates from 2013 to 2019. (d) θ difference between CABLE and observations. Note the different time axis for panels (c–d) relative to panels (a–b) due to different sampling intervals for soil moisture and fluxes.


3.4 Soil parameter optimisation experiment (Opt)

To address the simulated wet bias in the soil moisture profile (Fig. 5), we used observations to prescribe the critical soil hydraulic parameters θw and θsat (Fig. S3) and to optimise Ksat (Figs. S4 and S5). Prescribing θw and θsat led to a much improved “operating range” of soil moisture in the top 0.25 cm (Fig. S3b) but did not reduce the wet bias in the soil profile or solve the slow drainage after rainfall events (cf. Figs. 5c and 2c). Overall, these changes had a limited effect on simulated Etr (344 vs. 327 mm yr−1 in Hi-Res-2 in Table 3) as might be expected because the profile was sufficiently wet as not to limit evapotranspiration, especially in the root zone of the top 1.5 m (Fig. S5d). A reduction of the simulated Es (138 vs. 165 mm yr−1 in Hi-Res-2; Table 3) was mainly associated with the drier shallow soil (Fig. S5b). The optimised Ksat increased drainage speed (cf. Figs. 5c and 3c) and reduced the overall wet biases (0.04 m3 m−3 in Opt vs. 0.07 m3 m−3 in Hi-Res-2, MBE in Table 2).

3.5 Soil water limitation on transpiration (β-hvrd and β-exp)

Replacing CABLE's default soil moisture stress function with an alternative hypothesis that plants maximise their root water uptake to exploit resources (β-hvrd) led to a substantial increase in Etr (Fig. 6a) relative to experiment Opt (from 344 to 403 mm yr−1, Table 3) because the function assumes that the soil water stress on plants is determined by the availability of water in the wettest soil layer. This overestimation of Etr led to a small reduction in the wet soil moisture bias (cf. Figs. S5d and 6d).

Figure 6Haverd water stress function experiment (β-hvrd). (a) Etr and Es between 2013 and 2015. The shaded areas represent uncertainty between three ambient rings. Both simulations and observations are smoothed with a 3 d window to aid visualisation. (b) θ in the top 0.25 m from 2013 to 2019. (c) The vertical distribution of θ in β-hvrd at observed dates from 2013 to 2019. (d) θ difference between CABLE and observations. Note the different time axis for panels (c, d) relative to panels (a, b) due to different sampling intervals for soil moisture and fluxes.


Figure 7 shows the impact of using a site-calibrated β function (β-exp) (Yang et al., 2020). Using this function also increased Etr relative to experiment Opt (from 344 to 373 mm yr−1, Table 3), degrading the simulation relative to the standard β (Opt). In both experiments, owing to the overall simulated wet bias in the soil profile, a decreased sensitivity to soil moisture availability (either using β-hvrd or β-exp) did not improve simulated evapotranspiration.

Figure 7Site-based water stress function experiment (β-exp). (a) Etr and Es between 2013 and 2015. The shaded areas represent uncertainty between three ambient rings. Both simulations and observations are smoothed with a 3 d window to aid visualisation. (b) θ in the top 0.25 m from 2013 to 2019. (c) The vertical distribution of θ in β-exp at observed dates from 2013 to 2019. (d) θ difference between CABLE and observations. Note the different time axis for panels (c, d) relative to panels (a, b) due to different sampling intervals for soil moisture and fluxes.


3.6 Implications for drought

Improving the simulation of Etr, Es and soil moisture in LSMs is important on the seasonal timescale, but the increasing use of models to simulate future drought highlights the value of examining how these improvements impact the expression of drought in LSMs. We focus on a period of extensive drought across southeastern Australia that begins in October 2017 and extends to the end of 2019. Due to rainfall data availability, we focus on the dry-down period between October 2017 and September 2018.

Figure 8Simulations for each experiment during the drought period (October 2017 to September 2018). (a) The root-zone soil moisture over the top 1.5 m (θ1.5 m) and rainfall (P; bars), with blue dots showing the observed soil moisture. (b) Etr, (c) water stress factor (β), (d) Es and (e) sensible heat (QH). All lines are smoothed with a 30 d window.


Figure 8 shows selected fluxes during the drought period, over which the soil slowly dries in the observations and in the models (Fig. 8a) and the shallow soil moisture is close to wilting point (e.g. Fig. 6b). The Sres experiment maintains the highest soil moisture throughout the drought period, and β-hvrd the lowest, with the range across all experiments exceeding 0.1 m3 m−3. These soil moisture variations lead to inconsistent behaviour in Etr (Fig. 8b) due to resulting differences in β (Fig. 8c). β-hvrd Etr is very high despite having the driest soil moisture (Fig. 8a). Because the calculation uses the wettest soil layer, this leads to notably muted temporal variation (Fig. 8c) relative to the other methods for β. The differences in soil moisture, and as a result β, lead to differences in Etr (Fig. 8b) of 20–50 mm per month until autumn/winter ( April–July), when lower evaporative demand leads to more similar simulations. Through summer ( November–February), Es varies markedly from around 10 mm month−1 (β-hvrd) to 35 mm month−1 (Ctl) (Fig. 8d). The differences in Etr and Es are mirrored by differences in QH (Fig. 8e), which varies by >30W m−2 between the experiments in October 2017 and March 2018.

Integrating the simulations over the drought period highlights the differences in simulating water stress (expressed as β) between experiments. Figure 9a shows that Sres and β-hvrd maintain a relatively high β during drought periods (median >0.7), while the remaining experiments are notably lower. The β-exp simulates median values of 0.63, which is notably higher than the Hi-Res-2 of 0.33 and Opt of 0.46. This difference originates from the calibrated functional form shown in Fig. 9b, where the exponent in the β-exp function leads to a delay in the onset (point of inflection) of moisture stress relative to the default linear function used in CABLE. Overall, in a single model, parameterisations led to a difference of 98 % relative to the averaged median of all the simulations for β simulated during drought.

Figure 9(a) Box plot of simulated β during a drought year (October 2017–September 2018) and all simulated years (2013–2019). The dashed line is the mean value of β in Ctl over the dry period. (b) β variance with root-zone soil moisture over the top 1.5 m (θ1.5 m) during all simulated years.


3.7 Implications for heatwaves

The link between soil moisture and heatwaves is well known (Teuling et al., 2010) and is usually examined in the context of a drying soil leading to higher QH relative to QE (as our simulations are uncoupled, we cannot examine the consequences of these changes on the boundary layer).

Figure 10 shows a heatwave that occurred on 19–22 January 2018, where the air temperatures exceeded 35 C for four consecutive days and exceeded 40 C on the last day (Fig. 10a). The evaporative fraction during the daytime (09:00–16:00 LT) is shown in Fig. 10b and highlights a remarkable range from ∼0.2 in Ctl to ∼0.7 in β-hvrd, suggesting much stronger evaporative cooling in β-hvrd. An obvious diurnal variation in evaporative fraction is characterised by a progressive decline from a peak at 09:00 LT. QE gradually declines through the four heatwave days (Fig. 10c) in all experiments. At the beginning of the heatwave (19 January) daytime QE ranges from >200W m−2 in β-hvrd and Sres to around 100 W m−2 in Ctl, Watr, Hi-Res-1, Hi-Res-2 and Opt. The differences in QE are mirrored by differences in QH (Fig. 10d), with daytime fluxes varying on the heatwave days by more than 150 W m−2.

Figure 10Simulations during an observed heatwave with relatively low soil moisture (19–22 January 2018). (a) Air temperature (Tair; in black) and soil moisture within root zone over the top 1.5 m (θ1.5 m). The black dashed line shows the 35 C threshold. (b) Evaporative fraction (EF; calculated for daytime conditions), (c) latent heat (QE) and (d) sensible heat (QH). The vertical dash lines in panels (b–d) are at a 6 h interval to assist plot reading. One day before the heatwave is also shown.


Figure 10c and d also highlight a key divergence in energy partitioning due to parameterisations and the emergent interactions with soil water availability. Models that show a pronounced midday depression in QE (e.g. Ctrl, Watr and Hi-Res-2) due to increasing diurnal vapour pressure deficit (D) and soil moisture stress show earlier diurnal peaks in QH (Fig. 10d). By contrast, parameterisations that are less limited by β (e.g. β-hvrd despite the lowest soil moisture, Fig. 10a) see an emergent shift in peak in QH to later in the afternoon. When coupled, these emergent differences due to the role of soil water availability – and importantly, how this is translated in canopy gas exchange via β – may have implications for surface interactions with the boundary layer.

Given the importance of the role of D during heat extremes, to further explore the role of high D on simulated Etr, we plotted modelled and measured transpiration as a function of binned D (Fig. 11). At high D (>2kPa), simulated Etr is overestimated. As the mismatch between simulated Etr and observed occurs at both low and high D (Fig. 11), the model improvements are unlikely to simply relate to an alternative parameterisation of the stomatal sensitivity to D but instead suggest a missing mechanism to limit canopy gas exchange with increasing D. The impact of this overestimation would likely have greater significance for summers with concurrent heatwaves and droughts (compound events that are common in Australia), as during heatwaves when D is higher the model would overestimate Etr, using up available soil moisture.

Figure 11Modelled hourly Etr compared with measured hourly Etr over 2013. The solid line represents the 1:1 line. The dashed line is the linear fit between modelled and measured Etr. Colours of dots indicate the range of vapour pressure deficit.


4 Discussion and conclusions

Land surface schemes used in climate models range in complexity, and different approaches translate into contrasting predictions of the exchange of carbon, energy and water (Fisher and Koven, 2020). Perhaps critically, how strongly the land is coupled to the atmosphere also varies widely and is typically attributed to soil moisture variability (Brantley et al., 2017; Dirmeyer, 2011; Guo et al., 2006). A key component of LSMs is how soil moisture availability impacts processes internal to the land model and, in turn, how these impact fluxes of carbon and water. In this paper we used a rich observational dataset from a water-limited site that experiences both high temperatures and pronounced periods of low rainfall to explore a range of alternative model-based assumptions at multiple timescales (daily to annual) within a single model framework. We focused on the capacity of the model to simulate both the state (soil moisture) and the fluxes (evapotranspiration and its components) at a water-limited site. Whilst our analysis is site specific, the issues indicated here have been reported to lead to systematic biases in LSMs across multiple sites (Ukkola et al., 2016a; Trugman et al., 2018). We demonstrated that the default simulation (Ctl, Fig. 2) was able to simulate good transpiration fluxes but for the wrong reasons: erroneously high soil evaporation with a marked wet soil moisture bias. Errors of this kind may not have been identified in previous LSM evaluations against eddy covariance data which mostly focus on QE (Best et al., 2015). Our results highlight a potential bias in model evaluations due to a limited capacity to assess soil moisture or the partitioning of evapotranspiration. We demonstrated that poor model behaviour could be overcome via four key steps: (i) reducing soil evaporation biases, (ii) correctly initialising the aquifer moisture content, (iii) adjusting soil parameters to match site conditions and (iv) replacing the function used to constrain transpiration as soil moisture becomes limiting. Since our study attempts to articulate the common issues in the simulated dry conditions in LSMs, we anticipate our findings would be applicable not only in many water-limited conditions but, equally, in more mesic systems too.

Given the critical role of drought-prone ecosystems in contributing to interannual variability in the land CO2 sink size (Ahlström et al., 2015), our approach has the potential to improve the representation of these systems in models. We note that despite these improvements we still simulated a persistent wet soil moisture bias (e.g. Fig. 5d). We think on balance this is unlikely to originate from not simulating a seasonal understorey transpiration as β-hvrd, which grossly overestimated overstorey transpiration and did not sufficiently dry out the profile (cf. Figs. S5d and 6d). Instead the soil moisture bias must relate to CABLE's representation of subsurface processes.

4.1 Soil evaporation

LSMs commonly overestimate soil evaporation especially under a sparse canopy or over bare land (De Kauwe et al., 2017; Swenson and Lawrence, 2014), suggesting this is a key model weakness. Errors in soil evaporation are rarely isolated in models and often contribute to errors in transpiration by limiting soil moisture availability later in the growing season (Ukkola et al., 2016b) as well as affecting the distribution of shallow versus deep soil moisture draw-down during drought. Here we used a simple approach that increased resistance to surface evaporation, approximating the role of surface litter (Decker et al., 2017). At this site, this increased resistance to surface evaporation improved agreement with observations (Sres; Fig. 3a) but did not resolve all biases. Notably, soil evaporation was not directly measured at the site, but instead derived from the change in observed soil moisture over the top 5 cm, while ignoring days following rain (when the soil evaporative flux would likely be largest) (Gimeno et al., 2018a). As such, it contains soil moisture changes due to transpiration from a seasonal grass understorey but ignores evaporation below the top 5 cm, complicating model evaluation. Nevertheless, the magnitude of CABLE biases points to systematic errors in soil evaporation that do not merely arise form observational uncertainty. Many soil evaporation schemes used in LSMs lack a physical basis, estimating soil evaporation from aerodynamic resistance and an empirical soil water stress and ignoring the role of soil pores. However, a number of studies using alternative process-based schemes have been shown to improve individual model simulations (Haverd and Cuntz, 2010; Lehmann et al., 2018; Or and Lehmann, 2019). For example, Swenson and Lawrence (2014) introduced a dry-surface-layer-based soil evaporation resistance into CLM to depict water diffusion from dry soil, reducing biases in evapotranspiration and total water storage relative to FLUXNET-MTE and GRACE datasets. Based on a pore-scale model (Haghighi and Or, 2015), Decker et al. (2017) added the resistances of capillary–viscous and boundary layer to the CABLE soil evaporation scheme, lowered the positive Es bias in springtime and improved seasonality of evapotranspiration. Hence, a focused intercomparison of competing approaches against data originating from different ecosystems would be a valuable area of future work.

4.2 Aquifer initialisation

Our results showed that the initialisation of the aquifer moisture store was critical to an improved simulation of the soil moisture profile. By default, CABLE equilibrates the aquifer state by assuming almost complete saturation at the start. If, as happened with the Ctl, the aquifer is initialised too wet, the simulated water table is too high and the water potential in the aquifer is unlikely to be below the lowest soil moisture layer, impeding vertical aquifer recharge. When we initialised from a drier starting position (Watr), the simulated soil moisture profile matched the observed better, with implications for other models using similar groundwater schemes (e.g. CLM4.5, Noah-MP, JULES and LEAFHYDRO). First, our results imply that LSMs that incorporate groundwater schemes need to be careful about aquifer initialisation because this strongly affects soil moisture dynamics. Second, there is no obvious solution to this initialisation and spin-up problem because drainage into the aquifer is a slow process, and it may take hundreds of years to reach a realistic equilibrium state. For global simulations, this suggests the need to a priori initialise the starting aquifer state and to assess against satellite-based products like GRACE (Döll et al., 2014; Niu et al., 2007) or implement offline spin-up using meteorological forcing consistent with the subsequent simulations. However, while spin-up with observations is attractive, when the resulting states are incorporated into a coupled global model, inconsistencies are inevitable. Third, CABLE currently assumes a spin-up approach for the aquifer identical to that of the soil moisture, iterating until state changes between sequences of years are smaller than some threshold. LSMs that employ similar iteration approaches (Gilbert et al., 2017) are likely to encounter similar problems to CABLE because the rate of drainage into the aquifer is very slow, leading to negligible changes between iterations and thus satisfying the criteria for equilibrium.

4.3 Soil layers and pedotransfer functions

LSMs typically define a fixed number of soil layers globally, anywhere up to 20 layers. Most LSMs assume constant parameters across the entire soil profile, either using an experimental look-up table based on soil classification or estimating parameters from empirical pedotransfer functions. We explored the implications of these assumptions by first increasing the number of soil layers to match the number of observed layers (Hi-Res-1; Fig. S2) and then implementing soil parameters that varied vertically based on site texture (Hi-Res-2; Fig. 5). Increasing the vertical resolution had a small impact on the soil moisture and fluxes but did improve the temporal variability in soil moisture compared to observations. The use of site soil texture better depicts the moisture distribution in the soil profile but led to a slightly degraded soil moisture simulation. These results again highlight uncertainties in the translation of soil texture information to hydraulic parameters for water retention and hydraulic conductivity via the empirical functions, which are derived from limited observations (Van Looy et al., 2017). The development in pedotransfer functions via machine learning or multi-model ensemble provides new avenues to reduce errors from parameters (Zhang and Schaap, 2017; Dai et al., 2019). High-resolution global soil datasets (e.g. SoilGrids, Hengl et al., 2017) covering multiple soil layers up to 2 m depth offer opportunities to improve LSM simulations of soil moisture by incorporating depth-varying soil parameters. It is noteworthy that these global datasets of soil hydraulic parameters (Montzka et al, 2017; Dai et al., 2019) have existed for several years but have not been widely used. Furthermore, at the EucFACE site, the observed soil texture information enabled the separation of parameter uncertainties from biases in process representations and model structural errors, a valuable step in better constraining LSM simulations.

4.4 Calibration of soil hydraulic parameters

A number of studies have used satellite-derived (passive and active microwave) estimates of soil moisture to optimise soil hydraulic parameters in the top few soil layers (Harrison et al., 2012). Clearly these approaches are a potential way to constrain LSMs globally given the plethora of satellite observations extending back to the 1970s. However, these approaches implicitly assume that improving near-surface soil moisture translates to improvements over the entire soil column, an assumption not supported by our results. Whilst the use of observation-constrained θw and θsat over the top 0.3 m improved the simulated dynamics of shallow soil, it did not result in a large reduction in the bias simulated in deeper soil moisture layers (Fig. S3). At this site, the inability to significantly improve soil moisture dynamics through calibration of soil hydraulic conductivity against observed soil moisture data likely relates to the complexity of the soil profile, which contains two clay layers at depth (30–80 and 300–450 cm). This vertical texture complexity meant that it was difficult to obtain unique parameter solutions that would sufficiently improve vertical drainage whilst simultaneously simulating moisture dynamics well (Fig. S5). On the other hand, the neutron probe measurements of soil moisture used for calibration also involve uncertainties (Gimeno et al. 2018a). The soil moisture estimates were derived by fitting two distinctive linear relationships between soil volumetric water content and raw neutron probe counts (see Fig. S6) for clay (below 3 m) and non-clay soil (above 3 m). As a result, the observation error would be greatest in layers where the soil type differs from the assumed soil type at that depth. However, the fitted relationships were robust, since clay soils largely dominated the deeper profile (below 3 m depth) and sand soils mostly dominated the shallow profile (above 3 m depth). Overall, our sensitivity experiments suggest that optimising soil properties alone is not sufficient, and calibration exercises should also account for vegetation information to reduce biases in subsurface processes.

4.5 Water stress functions

Studies commonly highlight the functions used to limit photosynthesis and stomatal conductance with water stress as a key weakness among models. The lack of theory in this space (Medlyn et al., 2016) has led to models employing a range of functions encompassing different shapes and sensitivities that are not constrained by data. Trugman et al. (2018) explored the role of soil moisture stress in simulated “potential” gross primary productivity (GPP) among CMIP5 models and argued that the functional form used to represent the effect of soil moisture stress was the major driver of carbon cycle uncertainty. Here we deliberately attempted to first resolve model biases through other avenues (e.g. soil evaporation, soil parameterisation), because it is likely that model biases originate from multiple sources (e.g. leaf area, soil moisture dynamics etc.). We were subsequently able to assess the capacity to then further improve model behaviour via the functional forms used to represent water stress.

We examined three alternative water stress functions: the linear θ-based function used in Ctl (common among models, e.g. SDGVM, Orchidee-CN and JULES), a non-linear θ-based β (β-exp) calibrated for this site from Yang et al. (2020) and a function based on Haverd et al. (2016) (β-hvrd). Haverd et al. (2016) hypothesised that plants optimise their root water uptake, only limiting function when water in the deepest accessible soil layer becomes limiting. They further argued that this behaviour did not vary among sites (and so species). De Kauwe et al. (2015) previously tested this hypothesis and demonstrated that it led to an underestimation of the effect of moisture stress, inconsistent with observations. Our results again show that this hypothesis is not supported by data and led to an overestimation of transpiration (Fig. 6) and little evidence of moisture stress (Fig. 9b). We found that, after reducing other model biases, the use of the calibrated β−exp function did reduce the simulated soil moisture stress (over the drought year, median β=0.63 vs. 0.33 in Hi-Res-2 and 0.46 in Opt; Fig. 9). Overall, the various experiments show markedly different median β (ranging from 0.67 to 0.99, considering all simulated years), consistent with previous evaluations that have highlighted differences in simulated β across models (De Kauwe et al., 2017; Medlyn et al., 2016; Powell et al., 2013; Trugman et al., 2018). However, our results highlight that differences originate as much from alternative model assumptions and biases (e.g. soil evaporation, soil parameters) as the functional forms themselves. Alternatives to the β functions have emerged to fill the theoretical gap, including plant hydraulic (Christoffersen et al., 2016; Xu et al., 2016) and stomatal optimality approaches (Sperry et al., 2017), but are yet to be widely adopted in LSMs (but see Eller et al., 2020; De Kauwe et al., 2020; Kennedy et al., 2019; Sabot et al., 2020). Replacing the empirical soil water stress factor by these plant physiology schemes reduces model arbitrariness associated with the representation of soil water stress and reduces the simulated biases in transpiration either over water deficit regions or in areas with obvious dry seasons (Bonan et al., 2014; De Kauwe et al., 2020; Sabot et al., 2020). We can envision that a wider application of these processes-based models will offer a chance to improve water stress representation in more LSMs.

4.6 Heatwaves

Differences between the versions of CABLE lead to a different initial soil moisture state at the beginning of a heatwave, ranging from ∼0.15m3 m−3 (β-hvrd) to 0.23 m3 m−3 (Sres) (Fig. 10). In addition to the impact of the initial state, differences between parameterisation also affect estimates of β, leading to large divergences in evaporative cooling during a heatwave. Consequently, some versions of CABLE respond to the heatwave with a reduction of QE and a higher QH during 12:00–16:00 LT, while other simulations maintain a high QE during this period and the high QH shifts to later afternoon (16:00–18:00 LT) (Fig. 10c and d). The magnitudes of QE and QH between simulations are also substantially different: Ctl would amplify a heatwave, warming and drying the boundary layer, while β-hvrd would tend to moisten and (relatively) cool the boundary layer. Many studies have shown that the land surface can play a key role in amplifying heatwaves (Hirsch et al., 2019a; Miralles et al., 2014; Teuling et al., 2010), and LSMs exhibit systematic biases in representing this feedback (Sippel et al., 2017; Ukkola et al., 2018b). For a mega-heatwave like the 2010 European heatwave, the contribution of the local surface to the sensible heat anomaly was 20 W m−2 (Schumacher et al., 2019). However, our results show the differences between parameterisations within a single LSM can result in a greater divergence than this value. Therefore, these feedbacks can be substantially changed through different parameterisations and, if coupled to an atmospheric model, may be large enough to change the frequency and magnitude of heatwaves within a model.

We also showed that at high D our model overestimated transpiration, which would have consequences for subsequent soil moisture availability. Renchon et al. (2018) recently highlighted this point at the Cumberland Plains eddy covariance site, which neighbours the EucFACE site. Yang et al. (2019) showed that the MAESPA canopy gas exchange model similarly overpredicted transpiration at high D, leading to an overprediction of annual transpiration by 19 %. By examining leaf gas exchange data, they demonstrated that the reduction of transpiration could be attributed to non-stomatal limitation of photosynthesis at high D. Although non-stomatal limitation is commonly observed under low soil moisture content (e.g. Zhou et al. 2013) and implemented in a number of LSMs (De Kauwe et al., 2015), non-stomatal limitation at high D has been much less commonly reported and is not, to our knowledge, implemented in any LSMs. To echo Yang et al. (2019), further data on non-stomatal limitation at high D should be a priority, to determine whether this mechanism is sufficiently widespread to warrant inclusion in LSMs.

4.7 Future directions

We have shown that improving a LSM for one water flux is achievable, but improving a model to capture individual components of evapotranspiration and the associated soil moisture state is more challenging. No single step is sufficient in isolation, and if observations only constrain one element of a model, biases can be transferred within a model. This can lead to a tendency to hide biases in seldom-observed states because soil moisture profiles are rarely measured along with aboveground fluxes. International observational networks (e.g. FLUXNET; Baldocchi et al., 2001) rarely report QE, QH and soil moisture through and below the root zone simultaneously, although soil moisture profiles do sometimes exist. Expanding observational networks to include soil moisture profiles could accelerate model development. The EucFACE dataset holds exceptional promise as a means of evaluating model simulations and refining new theory. It is freely available, contains observations of the complete water balance and captures responses to both droughts and heatwaves. More broadly, our results also speak for the importance of multi-variable model evaluation methods for LSMs (e.g. iLAMB (Hoffman et al., 2017) and HTESSEL (Orth et al., 2017)).

While the focus of our study has been primarily on the parameterisations of hydrology and subsurface processes, we did aim to minimise the uncertainties from the non-hydrological factors by using site characteristics, such as the aerodynamic conductance (determined by setting the canopy height) and photosynthesis parameters (e.g. maximum carboxylation rate and maximum rate of electron transport). However, due to the significance of these non-hydrological factors on evapotranspiration (e.g. Breil et al., 2020), further evaluation should be considered in future studies.

Finally, our results imply that caution is needed in the interpretation of simulated heatwaves and droughts in coupled climate models. The feedback via the land surface is a key component, and, as our model experiments show, a range of alternative approaches can produce very different coupling between the land and the atmosphere if embedded in a coupled model. Despite the difficulties in acquiring datasets of the complete water balance, as a community we need to find an avenue to better assess (coupled) model predictions. Critical zone observatory networks (Brantley et al., 2017) may be one means to better constrain models, but in all likelihood targeted field campaigns that collect observations of soil moisture, eddy covariance and the boundary layer are also needed.

Code and data availability

The CABLE code is available at (last access: 25 January 2021) (NCI, 2021) after registration. Here, we use CABLE revision r7278. Scripts for plotting and processing model outputs are available at (Mu, 2021). EucFACE observations are publicly available in Western Sydney University's archive at (Gimeno et al., 2018b) and at (Yang, 2019).


The supplement related to this article is available online at:

Author contributions

MGDK, MM, AJP and AMU put forward the general scientific questions, designed the model experiments, investigated the simulations and drafted the manuscript. TEG, BEM, JY and DSE endeavoured to collect, to process and to correct the EucFACE observations. All authors participated in the discussion and revision of the manuscript.

Competing interests

The authors declare that they have no conflict of interest.


Mengyuan Mu, Martin G. De Kauwe, Andy J. Pitman and Anna M. Ukkola acknowledge support from the Australian Research Council (ARC) Centre of Excellence for Climate Extremes (CE170100023). Mengyuan Mu acknowledges support from the UNSW University International Postgraduate Award (UIPA) scheme. Martin G. De Kauwe and Andy J. Pitman acknowledge support from the ARC Discovery Grant (DP190101823). Martin G. De Kauwe, Teresa E. Gimeno, Belinda E. Medlyn, Jinyan Yang and David S. Ellsworth acknowledge support from the NSW Research Attraction and Acceleration Program (RAAP). EucFACE is supported by the Australian Commonwealth Government in collaboration with Western Sydney University. The facility was constructed as part of the nation-building initiative of the Australian Government. We thank the National Computational Infrastructure at the Australian National University, an initiative of the Australian Government, for access to supercomputer resources. We sincerely appreciate Burhan Amiji and Vinod Kumar for the collection of the neutron probe measurements. We also thank Craig Barton and Craig McNamara for their excellent technical support and Mingkai Jiang for the suggestion on the EucFACE understorey.

Financial support

This research has been supported by the Australian Research Council Centre of Excellence for Climate Extremes (grant no. CE170100023), the Australian Research Council Discovery Grant (grant no. DP190101823) and the NSW Research Attraction and Acceleration Program (RAAP).

Review statement

This paper was edited by Ryan Teuling and reviewed by Rene Orth and one anonymous referee.


Ahlström, A., Raupach, M. R., Schurgers, G., Smith, B., Arneth, A., Jung, M., Reichstein, M., Canadell, J. G., Friedlingstein, P., Jain, A. K., Kato, E., Poulter, B., Sitch, S., Stocker, B. D., Viovy, N., Wang, Y. P., Wiltshire, A., Zaehle, S., and Zeng, N.: The dominant role of semi-arid ecosystems in the trend and variability of the land CO2 sink, J. Geophys. Res.-Space, 120, 4503–4518,, 2015. 

Allen, C. D., Breshears, D. D., and McDowell, N. G.: On underestimation of global vulnerability to tree mortality and forest die-off from hotter drought in the Anthropocene, Ecosphere, 6, 1–55,, 2015. 

Ault, T. R.: On the essentials of drought in a changing climate, Science, 368, 256–260,, 2020. 

Baldocchi, D., Falge, E., Gu, L., Olson, R., Hollinger, D., Running, S., Anthoni, P., Bernhofer, C., Davis, K., Evans, R., Fuentes, J., Goldstein, A., Katul, G., Law, B., Lee, X., Malhi, Y., Meyers, T., Munger, W., Oechel, W., Paw, K. T., Pilegaard, K., Schmid, H. P., Valentini, R., Verma, S., Vesala, T., Wilson, K., and Wofsy, S.: FLUXNET: A new tool to study the temporal and spatial variability of ecosystem–scale carbon dioxide, water vapor, and energy flux densities, B. Am. Meteorol. Soc., 82, 2415–2434,<2415:FANTTS>2.3.CO;2, 2001. 

Best, M. J., Pryor, M., Clark, D. B., Rooney, G. G., Essery, R. L. H., Ménard, C. B., Edwards, J. M., Hendry, M. A., Porson, A., Gedney, N., Mercado, L. M., Sitch, S., Blyth, E., Boucher, O., Cox, P. M., Grimmond, C. S. B., and Harding, R. J.: The Joint UK Land Environment Simulator (JULES), model description – Part 1: Energy and water fluxes, Geosci. Model Dev., 4, 677–699,, 2011. 

Best, M. J., Abramowitz, G., Johnson, H. R., Pitman, A. J., Balsamo, G., Boone, A., Cuntz, M., Decharme, B., Dirmeyer, P. A., Dong, J., Ek, M., Guo, Z., Haverd, V., van den Hurk, B. J. J., Nearing, G. S., Pak, B., Peters-Lidard, C., Santanello, J. A., Stevens, L., and Vuichard, N.: The plumbing of land surface models: benchmarking model performance, J. Hydrometeorol., 16, 1425–1442,, 2015. 

Bi, D., Dix, M., Marsland, S., O'Farrell, S., Rashid, H., Uotila, P., Hirst, A., Kowalczyk, E., Golebiewski, M., Sullivan, A., Yan, H., Hannah, N., Franklin, C., Sun, Z., Vohralik, P., Watterson, I., Zhou, X., Fiedler, R., Collier, M., Ma, Y., Noonan, J., Stevens, L., Uhe, P., Zhu, H., Griffies, S., Hill, R., Harris, C., and Puri, K.: The ACCESS coupled model: description, control climate and evaluation, Aust. Meteorol. Ocean., 63, 41–64,, 2013. 

Bonan, G. B., Williams, M., Fisher, R. A., and Oleson, K. W.: Modeling stomatal conductance in the earth system: linking leaf water-use efficiency and water transport along the soil-plant-atmosphere continuum, Geosci. Model Dev., 7, 2193–2222,, 2014. 

Brantley, S. L., McDowell, W. H., Dietrich, W. E., White, T. S., Kumar, P., Anderson, S. P., Chorover, J., Lohse, K. A., Bales, R. C., Richter, D. D., Grant, G., and Gaillardet, J.: Designing a network of critical zone observatories to explore the living skin of the terrestrial Earth, Earth Surf. Dynam., 5, 841–860,, 2017. 

Breil, M., Davin, E. L., and Rechid, D.: What determines the sign of the evapotranspiration response to afforestation in the European summer?, Biogeosciences Discuss. [preprint],, in review, 2020. 

Christoffersen, B. O., Gloor, M., Fauset, S., Fyllas, N. M., Galbraith, D. R., Baker, T. R., Kruijt, B., Rowland, L., Fisher, R. A., Binks, O. J., Sevanto, S., Xu, C., Jansen, S., Choat, B., Mencuccini, M., McDowell, N. G., and Meir, P.: Linking hydraulic traits to tropical forest function in a size-structured and trait-driven model (TFS v.1-Hydro), Geosci. Model Dev., 9, 4227–4255,, 2016. 

Ciais, P., Reichstein, M., Viovy, N., Granier, A., Ogée, J., Allard, V., Aubinet, M., Buchmann, N., Bernhofer, C., Carrara, A., Chevallier, F., De Noblet, N., Friend, A. D., Friedlingstein, P., Grünwald, T., Heinesch, B., Keronen, P., Knohl, A., Krinner, G., Loustau, D., Manca, G., Matteucci, G., Miglietta, F., Ourcival, J. M., Papale, D., Pilegaard, K., Rambal, S., Seufert, G., Soussana, J. F., Sanz, M. J., Schulze, E. D., Vesala, T., and Valentini, R.: Europe-wide reduction in primary productivity caused by the heat and drought in 2003, Nature, 437, 529–533,, 2005. 

Clark, M. P., Fan, Y., Lawrence, D. M., Adam, J. C., Bolster, D., Gochis, D. J., Hooper, R. P., Kumar, M., Leung, L. R., Mackay, D. S., Maxwell, R. M., Shen, C., Swenson, S. C., and Zeng, X.: Improving the representation of hydrologic processes in Earth System Models, Water Resour. Res., 51, 5929–5956,, 2015. 

Collins, L., Bradstock, R. A., Resco de Dios, V., Duursma, R. A., Velasco, S., and Boer, M. M.: Understorey productivity in temperate grassy woodland responds to soil water availability but not to elevated [CO2], Global Change Biol., 24, 2366–2376,, 2018. 

Collins, M., Knutti, R., Arblaster, J., Dufresne, J.-L., Fichefet, T., Friedlingstein, P., Gao, X., Gutowski, W. J., Johns, T., Krinner, G., Shongwe, M., Tebaldi, C., Weaver, A. J., and Wehner, M.: Long-term climate change: projections, commitments and irreversibility, in: Climate Change 2013: The Physical Science Basis, Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., Midgley, P. M., Cambridge University Press, Cambridge, UK and New York, NY, USA, 1029–1136, 2013. 

Cosby, B. J., Hornberger, G. M., Clapp, R. B., and Ginn, T. R.: A statistical exploration of the relationships of soil moisture characteristics to the physical properties of soils, Water Resour. Res., 20, 682–690,, 1984. 

Crous, K. Y., Ósvaldsson, A., and Ellsworth, D. S.: Is phosphorus limiting in a mature Eucalyptus woodland? Phosphorus fertilisation stimulates stem growth, Plant Soil, 391, 293–305,, 2015. 

Dai, Y., Shangguan, W., Duan, Q., Liu, B., Fu, S., and Niu, G.: Development of a China dataset of soil hydraulic parameters using pedotransfer functions for land surface modeling, J. Hydrometeorol., 14, 869–887,, 2013. 

Dai, Y., Xin, Q., Wei, N., Zhang, Y., Shangguan, W., Yuan, H., Zhang, S., Liu, S., and Lu, X.: A Global High-Resolution Data Set of Soil Hydraulic and Thermal Properties for Land Surface Modeling, J. Adv. Model. Earth Syst., 11, 2996–3023,, 2019. 

Decker, M.: Development and evaluation of a new soil moisture and runoff parameterization for the CABLE LSM including subgrid-scale processes, J. Adv. Model. Earth Syst., 7, 1788–1809,, 2015. 

Decker, M. and Zeng, X.: Impact of modified Richards equation on global soil moisture simulation in the Community Land Model (CLM3.5), J. Adv. Model. Earth Syst., 1, 1–22,, 2009. 

Decker, M., Or, D., Pitman, A., and Ukkola, A.: New turbulent resistance parameterization for soil evaporation based on a pore-scale model: Impact on surface fluxes in CABLE, J. Adv. Model. Earth Syst., 9, 220–238,, 2017. 

De Kauwe, M. G., Medlyn, B. E., Zaehle, S., Walker, A. P., Dietze, M. C., Hickler, T., Jain, A. K., Luo, Y., Parton, W. J., Prentice, I. C., Smith, B., Thornton, P. E., Wang, S., Wang, Y. P., Wårlind, D., Weng, E., Crous, K. Y., Ellsworth, D. S., Hanson, P. J., Seok Kim, H., Warren, J. M., Oren, R., and Norby, R. J.: Forest water use and water use efficiency at elevated CO2: A model-data intercomparison at two contrasting temperate forest FACE sites, Global Change Biol., 19, 1759–1779,, 2013. 

De Kauwe, M. G., Zhou, S.-X., Medlyn, B. E., Pitman, A. J., Wang, Y.-P., Duursma, R. A., and Prentice, I. C.: Do land surface models need to include differential plant species responses to drought? Examining model predictions across a mesic-xeric gradient in Europe, Biogeosciences, 12, 7503–7518,, 2015. 

De Kauwe, M. G., Medlyn, B. E., Walker, A. P., Zaehle, S., Asao, S., Guenet, B., Harper, A. B., Hickler, T., Jain, A. K., Luo, Y., Lu, X., Luus, K., Parton, W. J., Shu, S., Wang, Y. P., Werner, C., Xia, J., Pendall, E., Morgan, J. A., Ryan, E. M., Carrillo, Y., Dijkstra, F. A., Zelikova, T. J., and Norby, R. J.: Challenging terrestrial biosphere models with data from the long-term multifactor Prairie Heating and CO2 Enrichment experiment, Global Change Biol., 23, 3623–3645,, 2017. 

De Kauwe, M. G., Medlyn, B. E., Ukkola, A. M., Mu, M., Sabot, M. E. B., Pitman, A. J., Meir, P., Cernusak, L., Rifai, S. W., Choat, B., Tissue, D. T., Blackman, C. J., Li, X., Roderick, M., and Briggs, P. R.: Identifying areas at risk of drought-induced tree mortality across South-Eastern Australia, Global Change Biol., 26, 5716–5733,, 2020. 

Dewar, R., Mauranen, A., Mäkelä, A., Hölttä, T., Medlyn, B., and Vesala, T.: New insights into the covariation of stomatal, mesophyll and hydraulic conductances from optimization models incorporating nonstomatal limitations to photosynthesis, New Phytol., 217, 571–585,, 2018. 

Dirmeyer, P. A.: A history and review of the Global Soil Wetness Project (GSWP), J. Hydrometeorol., 12, 729–749,, 2011. 

Döll, P., Muller, H. S., Schuh, C., Portmann, F. T., Eicker, A., Schmied, H. M., Schuh, C., Portmann, F. T., and Eicker, A.: Global-scale assessment of groundwater depletion and related groundwater abstractions: Combining hydrological modeling with information from well observations and GRACE satellites, J. Am. Water Resour. Assoc., 5, 2–2,, 2014. 

Dosio, A., Mentaschi, L., Fischer, E. M., and Wyser, K.: Extreme heat waves under 1.5 C and 2 C global warming, Environ. Res. Lett., 13, 054006,, 2018. 

Duursma, R. A., Gimeno, T. E., Boer, M. M., Crous, K. Y., Tjoelker, M. G., and Ellsworth, D. S.: Canopy leaf area of a mature evergreen Eucalyptus woodland does not respond to elevated atmospheric [CO2] but tracks water availability, Global Change Biol., 22, 1666–1676,, 2016. 

Egea, G., Verhoef, A., and Vidale, P. L.: Towards an improved and more flexible representation of water stress in coupled photosynthesis–stomatal conductance models, Agr. Forest Meteorol., 151, 1370–1384,, 2011. 

Eller, C. B., Rowland, L., Mencuccini, M., Rosas, T., Williams, K., Harper, A., Medlyn, B. E., Wagner, Y., Klein, T., Teodoro, G. S., Oliveira, R. S., Matos, I. S., Rosado, B. H. P., Fuchs, K., Wohlfahrt, G., Montagnani, L., Meir, P., Sitch, S., and Cox, P. M.: Stomatal optimization based on xylem hydraulics (SOx) improves land surface model simulation of vegetation responses to climate, New Phytol., 226, 1622–1637,, 2020. 

Ellsworth, D. S., Anderson, I. C., Crous, K. Y., Cooke, J., Drake, J. E., Gherlenda, A. N., Gimeno, T. E., Macdonald, C. A., Medlyn, B. E., Powell, J. R., Tjoelker, M. G., and Reich, P. B.: Elevated CO2 does not increase eucalypt forest productivity on a low-phosphorus soil, Nat. Clim. Change, 7, 279–282,, 2017. 

Fischer, G., Nachtergaele, F., Prieler, S., van Velthuizen, H. T., Verelst, L., and Wiberg, D.: Global Agro-ecological Zones Assessment for Agriculture (GAEZ 2008), IIASA, Laxenburg, Austria and FAO, Rome, Italy, 2008. 

Fisher, R. A. and Koven, C. D.: Perspectives on the future of land surface models and the challenges of representing complex terrestrial systems, J. Adv. Model. Earth Syst., 12, e2018MS001453,, 2020. 

Gherlenda, A. N., Esveld, J. L., Hall, A. A. G., Duursma, R. A., and Riegler, M.: Boom and bust: rapid feedback responses between insect outbreak dynamics and canopy leaf area impacted by rainfall and CO2, Global Change Biol., 22, 3632–3641,, 2016. 

Gilbert, J. M., Maxwell, R. M., and Gochis, D. J.: Effects of water-table configuration on the planetary boundary layer over the San Joaquin River Watershed, California, J. Hydrometeorol., 18, 1471–1488,, 2017. 

Gimeno, T. E., Crous, K. Y., Cooke, J., O'Grady, A. P., Ósvaldsson, A., Medlyn, B. E., and Ellsworth, D. S.: Conserved stomatal behaviour under elevated CO2 and varying water availability in a mature woodland, edited by D. Whitehead, Funct. Ecol., 30, 700–709,, 2016. 

Gimeno, T. E., McVicar, T. R., O'Grady, A. P., Tissue, D. T., and Ellsworth, D. S.: Elevated CO2 did not affect the hydrological balance of a mature native Eucalyptus woodland, Global Change Biol., 24, 3010–3024,, 2018a. 

Gimeno, T. E., McVicar, T. R., O'Grady, A. P., Tissue, D. T., and Ellsworth, D. S.: EucFACE Hydrological and meteorological measurements from 2012-04-30 to 2014-11-15, Western Sydney University, Sydney, Australia,, 2018b. 

Guo, Z., Dirmeyer, P. A., Koster, R. D., Sud, Y. C., Bonan, G., Oleson, K. W., Chan, E., Verseghy, D., Cox, P., Gordon, C. T., McGregor, J. L., Kanae, S., Kowalczyk, E., Lawrence, D., Liu, P., Mocko, D., Lu, C.-H., Mitchell, K., Malyshev, S., McAvaney, B., Oki, T., Yamada, T., Pitman, A., Taylor, C. M., Vasic, R., and Xue, Y.: GLACE: the Global Land–Atmosphere Coupling Experiment. part II: analysis, J. Hydrometeorol., 7, 611–625,, 2006. 

Haghighi, E. and Or, D.: Linking evaporative fluxes from bare soil across surface viscous sublayer with the Monin–Obukhov atmospheric flux-profile estimates, J. Hydrol., 525, 684–693,, 2015. 

Harrison, K. W., Kumar, S. V., Peters-Lidard, C. D., and Santanello, J. A.: Quantifying the change in soil moisture modeling uncertainty from remote sensing observations using Bayesian inference techniques, Water Resour. Res., 48, W11514,, 2012. 

Haverd, V. and Cuntz, M.: Soil–Litter–Iso: A one-dimensional model for coupled transport of heat, water and stable isotopes in soil with a litter layer and root extraction, J. Hydrol., 388, 438–455,, 2010. 

Haverd, V., Raupach, M. R., Briggs, P. R., Canadell, J. G., Isaac, P., Pickett-Heaps, C., Roxburgh, S. H., van Gorsel, E., Viscarra Rossel, R. A., and Wang, Z.: Multiple observation types reduce uncertainty in Australia's terrestrial carbon and water cycles, Biogeosciences, 10, 2011–2040,, 2013. 

Haverd, V., Cuntz, M., Nieradzik, L. P., and Harman, I. N.: Improved representations of coupled soil-canopy processes in the CABLE land surface model (Subversion revision 3432), Geosci. Model Dev., 9, 3111–3122,, 2016. 

Hengl, T., Mendes de Jesus, J., Heuvelink, G. B. M., Ruiperez Gonzalez, M., Kilibarda, M., Blagotić, A., Shangguan, W., Wright, M. N., Geng, X., Bauer-Marschallinger, B., Guevara, M. A., Vargas, R., MacMillan, R. A., Batjes, N. H., Leenaars, J. G. B., Ribeiro, E., Wheeler, I., Mantel, S., and Kempen, B.: SoilGrids250m: Global gridded soil information based on machine learning, edited by B. Bond-Lamberty, PLoS One, 12, e0169748,, 2017. 

Hirsch, A. L., Evans, J. P., Di Virgilio, G., Perkins-Kirkpatrick, S. E., Argüeso, D., Pitman, A. J., Carouge, C. C., Kala, J., Andrys, J., Petrelli, P., and Rockel, B.: Amplification of Australian heatwaves via local land-atmosphere coupling, J. Geophys. Res.-Atmos., 124, 13625–13647,, 2019a. 

Hirsch, A. L., Kala, J., Carouge, C. C., De Kauwe, M. G., Di Virgilio, G., Ukkola, A. M., Evans, J. P., and Abramowitz, G.: Evaluation of the CABLEv2.3.4 Land Surface Model Coupled to NU-WRFv3.9.1.1 in simulating temperature and precipitation means and extremes Over CORDEX AustralAsia within a WRF physics ensemble, J. Adv. Model. Earth Syst., 11, 4466–4488,, 2019b. 

Hoffman, F. M., Koven, C. D., Keppel-Aleks, G., Lawrence, D. M., Riley, W. J., Randerson, J. T., Ahlström, A., Abramowitz, G., Baldocchi, D. D., Best, M. J., Bond-Lamberty, B., De Kauwe, M. G., Denning, A. S., Desai, A. R., Eyring, V., Fisher, J. B., Fisher, R. A., Gleckler, P. J., Huang, M., Hugelius, G., Jain, A. K., Kiang, N. Y., Kim, H., Koster, R. D., Kumar, S. V., Li, H., Luo, Y., Mao, J., McDowell, N. G., Mishra, U., Moorcroft, P. R., Pau, G. S. H., Ricciuto, D. M., Schaefer, K., Schwalm, C. R., Serbin, S. P., Shevliakova, E., Slater, A. G., Tang, J., Williams, M., Xia, J., Xu, C., Joseph, R., and Koch, D.: 2016 International Land Model Benchmarking (ILAMB) Workshop Report, USA,, 2017. 

Humphrey, V., Zscheischler, J., Ciais, P., Gudmundsson, L., Sitch, S., and Seneviratne, S. I.: Sensitivity of atmospheric CO2 growth rate to observed changes in terrestrial water storage, Nature, 560, 628–631,, 2018. 

Kala, J., De Kauwe, M. G., Pitman, A. J., Lorenz, R., Medlyn, B. E., Wang, Y.-P., Lin, Y.-S., and Abramowitz, G.: Implementation of an optimal stomatal conductance scheme in the Australian Community Climate Earth Systems Simulator (ACCESS1.3b), Geosci. Model Dev., 8, 3877–3889,, 2015. 

Kennedy, D., Swenson, S., Oleson, K. W., Lawrence, D. M., Fisher, R., Lola da Costa, A. C., and Gentine, P.: Implementing Plant Hydraulics in the Community Land Model, Version 5, J. Adv. Model. Earth Syst., 11, 485–513,, 2019. 

Kishné, A. S., Yimam, Y. T., Morgan, C. L. S., and Dornblaser, B. C.: Evaluation and improvement of the default soil hydraulic parameters for the Noah Land Surface Model, Geoderma, 285, 247–259,, 2017. 

Klein, T.: The variability of stomatal sensitivity to leaf water potential across tree species indicates a continuum between isohydric and anisohydric behaviours, Funct. Ecol., 28, 1313–1320,, 2014. 

Kowalczyk, E. A., Wang, Y. P., and Law, R. M.: The CSIRO Atmosphere Biosphere Land Exchange (CABLE) model for use in climate models and as an offline model, CSIRO Mar. Atmos. Res. Pap. 13, CSIRO, Australia, 42 pp.,, 2006. 

Lai, C. T. and Katul, G.: The dynamic role of root-water uptake in coupling potential to actual transpiration, Adv. Water Resour., 23, 427–439,, 2000. 

Law, R. M., Ziehn, T., Matear, R. J., Lenton, A., Chamberlain, M. A., Stevens, L. E., Wang, Y.-P., Srbinovsky, J., Bi, D., Yan, H., and Vohralik, P. F.: The carbon cycle in the Australian Community Climate and Earth System Simulator (ACCESS-ESM1) – Part 1: Model description and pre-industrial simulation, Geosci. Model Dev., 10, 2567–2590,, 2017. 

Lawrence, D. M., Fisher, R. A., Koven, C. D., Oleson, K. W., Swenson, S. C., Bonan, G., Collier, N., Ghimire, B., van Kampenhout, L., Kennedy, D., Kluzek, E., Lawrence, P. J., Li, F., Li, H., Lombardozzi, D., Riley, W. J., Sacks, W. J., Shi, M., Vertenstein, M., Wieder, W. R., Xu, C., Ali, A. A., Badger, A. M., Bisht, G., van den Broeke, M., Brunke, M. A., Burns, S. P., Buzan, J., Clark, M., Craig, A., Dahlin, K., Drewniak, B., Fisher, J. B., Flanner, M., Fox, A. M., Gentine, P., Hoffman, F., Keppel-Aleks, G., Knox, R., Kumar, S., Lenaerts, J., Leung, L. R., Lipscomb, W. H., Lu, Y., Pandey, A., Pelletier, J. D., Perket, J., Randerson, J. T., Ricciuto, D. M., Sanderson, B. M., Slater, A., Subin, Z. M., Tang, J., Thomas, R. Q., Val Martin, M., and Zeng, X.: The Community Land Model Version 5: description of new features, benchmarking, and impact of forcing uncertainty, J. Adv. Model. Earth Syst., 11, 4245–4287,, 2019. 

Lawrence, P. J. and Chase, T. N.: Representing a new MODIS consistent land surface in the Community Land Model (CLM 3.0), J. Geophys. Res., 112, G01023,, 2007. 

Lehmann, P., Merlin, O., Gentine, P., and Or, D.: Soil texture effects on surface resistance to bare-soil evaporation, Geophys. Res. Lett., 45, 10398–10405,, 2018. 

Li, L., Wang, Y. P., Yu, Q., Pak, B., Eamus, D., Yan, J., Van Gorsel, E., and Baker, I. T.: Improving the responses of the Australian community land surface model (CABLE) to seasonal drought, J. Geophys. Res-Biogeo, 117, 1–16,, 2012. 

Lian, X., Piao, S., Huntingford, C., Li, Y., Zeng, Z., Wang, X., Ciais, P., McVicar, T. R., Peng, S., Ottlé, C., Yang, H., Yang, Y., Zhang, Y., and Wang, T.: Partitioning global land evapotranspiration using CMIP5 models constrained by observations, Nat. Clim. Change, 8, 640–646,, 2018. 

Lin, Y.-S., Medlyn, B. E., Duursma, R. a., Prentice, I. C., Wang, H., Baig, S., Eamus, D., de Dios, V. R., Mitchell, P., Ellsworth, D. S., de Beeck, M. O., Wallin, G., Uddling, J., Tarvainen, L., Linderson, M.-L., Cernusak, L. A., Nippert, J. B., Ocheltree, T. W., Tissue, D. T., Martin-StPaul, N. K., Rogers, A., Warren, J. M., De Angelis, P., Hikosaka, K., Han, Q., Onoda, Y., Gimeno, T. E., Barton, C. V. M., Bennie, J., Bonal, D., Bosc, A., Löw, M., Macinins-Ng, C., Rey, A., Rowland, L., Setterfield, S. a., Tausz-Posch, S., Zaragoza-Castells, J., Broadmeadow, M. S. J., Drake, J. E., Freeman, M., Ghannoum, O., Hutley, L. B., Kelly, J. W., Kikuzawa, K., Kolari, P., Koyama, K., Limousin, J.-M., Meir, P., Lola da Costa, A. C., Mikkelsen, T. N., Salinas, N., Sun, W., and Wingate, L.: Optimal stomatal behaviour around the world, Nat. Clim. Change, 5, 459–464,, 2015. 

Lorenz, R., Pitman, A. J., Donat, M. G., Hirsch, A. L., Kala, J., Kowalczyk, E. A., Law, R. M., and Srbinovsky, J.: Representation of climate extreme indices in the ACCESS1.3b coupled atmosphere-land surface model, Geosci. Model Dev., 7, 545–567,, 2014. 

Matthews, T. K. R. R., Wilby, R. L., and Murphy, C.: Communicating the deadly consequences of global warming for human heat stress, P. Natl. Acad. Sci. USA, 114, 3861–3866,, 2017. 

Mazdiyasni, O. and AghaKouchak, A.: Substantial increase in concurrent droughts and heatwaves in the United States, P. Natl. Acad. Sci. USA, 112, 11484–11489,, 2015. 

Medlyn, B. E., Duursma, R. A., Eamus, D., Ellsworth, D. S., Prentice, I. C., Barton, C. V. M., Crous, K. Y., De Angelis, P., Freeman, M., and Wingate, L.: Reconciling the optimal and empirical approaches to modelling stomatal conductance, Global Change Biol., 17, 2134–2144,, 2011. 

Medlyn, B. E., De Kauwe, M. G., Zaehle, S., Walker, A. P., Duursma, R. A., Luus, K., Mishurov, M., Pak, B., Smith, B., Wang, Y. P., Yang, X., Crous, K. Y., Drake, J. E., Gimeno, T. E., Macdonald, C. A., Norby, R. J., Power, S. A., Tjoelker, M. G., and Ellsworth, D. S.: Using models to guide field experiments: a priori predictions for the CO2 response of a nutrient- and water-limited native Eucalypt woodland, Global Change Biol., 22, 2834–2851,, 2016. 

Miralles, D. G., Teuling, A. J., Van Heerwaarden, C. C., and De Arellano, J. V. G.: Mega-heatwave temperatures due to combined soil desiccation and atmospheric heat accumulation, Nat. Geosci., 7, 345–349,, 2014. 

Miralles, D. G., Gentine, P., Seneviratne, S. I., and Teuling, A. J.: Land-atmospheric feedbacks during droughts and heatwaves: state of the science and current challenges, Ann. NY Acad. Sci., 1436, 19–35,, 2019. 

Montzka, C., Herbst, M., Weihermüller, L., Verhoef, A., and Vereecken, H.: A global data set of soil hydraulic properties and sub-grid variability of soil water retention and hydraulic conductivity curves, Earth Syst. Sci. Data, 9, 529–543,, 2017. 

Mu, M.: bibivking/Evaluate_CABLE_EucFACE, Zenodo,, 2021. 

NCI: CABLE: The Community Atmosphere Biosphere Land Exchange Model, available at:, last access: 25 January 2021. 

Niu, G.-Y., Yang, Z.-L., Dickinson, R. E., Gulden, L. E., and Su, H.: Development of a simple groundwater model for use in climate models and evaluation with Gravity Recovery and Climate Experiment data, J. Geophys. Res., 112, D07103,, 2007. 

Oleson, K. W., Niu, G.-Y., Yang, Z.-L., Lawrence, D. M., Thornton, P. E., Lawrence, P. J., Stöckli, R., Dickinson, R. E., Bonan, G. B., Levis, S., Dai, A., and Qian, T.: Improvements to the Community Land Model and their impact on the hydrological cycle, J. Geophys. Res.-Biogeo., 113, 1–26,, 2008. 

Or, D. and Lehmann, P.: Surface evaporative capacitance: how soil type and rainfall characteristics affect global-scale surface evaporation, Water Resour. Res., 55, 519–539,, 2019. 

Orth, R. and Destouni, G.: Drought reduces blue-water fluxes more strongly than green-water fluxes in Europe, Nat. Commun., 9, 1–8,, 2018. 

Orth, R., Dutra, E., Trigo, I. F., and Balsamo, G.: Advancing land surface model development with satellite-based Earth observations, Hydrol. Earth Syst. Sci., 21, 2483–2495,, 2017. 

Pal, J. S. and Eltahir, E. A. B.: Future temperature in southwest Asia projected to exceed a threshold for human adaptability, Nat. Clim. Change, 6, 197–200,, 2016. 

Pan, S., Pan, N., Tian, H., Friedlingstein, P., Sitch, S., Shi, H., Arora, V. K., Haverd, V., Jain, A. K., Kato, E., Lienert, S., Lombardozzi, D., Nabel, J. E. M. S., Ottlé, C., Poulter, B., Zaehle, S., and Running, S. W.: Evaluation of global terrestrial evapotranspiration using state-of-the-art approaches in remote sensing, machine learning and land surface modeling, Hydrol. Earth Syst. Sci., 24, 1485–1509,, 2020. 

Pathare, V. S., Crous, K. Y., Cooke, J., Creek, D., Ghannoum, O., and Ellsworth, D. S.: Water availability affects seasonal CO2-induced photosynthetic enhancement in herbaceous species in a periodically dry woodland, Global Change Biol., 23, 5164–5178,, 2017. 

Powell, T. L., Galbraith, D. R., Christoffersen, B. O., Harper, A., Imbuzeiro, H. M. A., Rowland, L., Almeida, S., Brando, P. M., da Costa, A. C. L., Costa, M. H., Levine, N. M., Malhi, Y., Saleska, S. R., Sotta, E., Williams, M., Meir, P., and Moorcroft, P. R.: Confronting model predictions of carbon fluxes with measurements of Amazon forests subjected to experimental drought, New Phytol., 200, 350–365,, 2013. 

Puhlmann, H., and von Wilpert, K.: Pedotransfer functions for water retention and unsaturated hydraulic conductivity of forest soils, J. Plant Nutr. Soil Sci., 175, 221–235,, 2012. 

Reichstein, M., Bahn, M., Ciais, P., Frank, D., Mahecha, M. D., Seneviratne, S. I., Zscheischler, J., Beer, C., Buchmann, N., Frank, D. C., Papale, D., Rammig, A., Smith, P., Thonicke, K., van der Velde, M., Vicca, S., Walz, A., and Wattenbach, M.: Climate extremes and the carbon cycle, Nature, 500, 287–295,, 2013. 

Renchon, A. A., Griebel, A., Metzen, D., Williams, C. A., Medlyn, B., Duursma, R. A., Barton, C. V. M., Maier, C., Boer, M. M., Isaac, P., Tissue, D., Resco de Dios, V., and Pendall, E.: Upside-down fluxes Down Under: CO2 net sink in winter and net source in summer in a temperate evergreen broadleaf forest, Biogeosciences, 15, 3703–3716,, 2018. 

Sabot, M. E. B. B., De Kauwe, M. G., Pitman, A. J., Medlyn, B. E., Verhoef, A., Ukkola, A. M., and Abramowitz, G.: Plant profit maximization improves predictions of European forest responses to drought, New Phytol., 226, 1638–1655,, 2020. 

Sakaguchi, K. and Zeng, X.: Effects of soil wetness, plant litter, and under-canopy atmospheric stability on ground evaporation in the Community Land Model (CLM3.5), J. Geophys. Res.-Atmos., 114, 1–14,, 2009. 

Schlosser, C. A., Slater, A. G., Robock, A., Pitman, A. J., Vinnikov, K. Y., Henderson-Sellers, A., Speranskaya, N. A., Mitchell, K.: Simulations of a boreal grassland hydrology at Valdai, Russia: PILPS Phase 2(d), Mon. Weather Rev., 128, 301–321,<0301:SOABGH>2.0.CO;2, 2000. 

Schumacher, D. L., Keune, J., van Heerwaarden, C. C., Vilà-Guerau de Arellano, J., Teuling, A. J., and Miralles, D. G.: Amplification of mega-heatwaves through heat torrents fuelled by upwind drought, Nat. Geosci., 12, 712–717,, 2019. 

Seneviratne, S. I., Corti, T., Davin, E. L., Hirschi, M., Jaeger, E. B., Lehner, I., Orlowsky, B., and Teuling, A. J.: Investigating soil moisture-climate interactions in a changing climate: A review, Earth-Sci. Rev., 99, 125–161,, 2010. 

Sippel, S., Zscheischler, J., Mahecha, M. D., Orth, R., Reichstein, M., Vogel, M., and Seneviratne, S. I.: Refining multi-model projections of temperature extremes by evaluation against land-atmosphere coupling diagnostics, Earth Syst. Dynam., 8, 387–403,, 2017. 

Sperry, J. S., Venturas, M. D., Anderegg, W. R. L., Mencuccini, M., Mackay, D. S., Wang, Y., and Love, D. M.: Predicting stomatal responses to the environment from the optimization of photosynthetic gain and hydraulic cost, Plant. Cell Environ., 40, 816–830,, 2017. 

Swenson, S. C. and Lawrence, D. M.: Assessing a dry surface layer-based soil resistance parameterization for the Community Land Model using GRACE and FLUXNET-MTE data, J. Geophys. Res.-Atmos., 119, 10299–10312,, 2014. 

Tallaksen, L. M. and Stahl, K.: Spatial and temporal patterns of large-scale droughts in Europe: Model dispersion and performance, Geophys. Res. Lett., 41, 429–434,, 2014. 

Teuling, A. J., Seneviratne, S. I., Stöckli, R., Reichstein, M., Moors, E., Ciais, P., Luyssaert, S., van den Hurk, B., Ammann, C., Bernhofer, C., Dellwik, E., Gianelle, D., Gielen, B., Grünwald, T., Klumpp, K., Montagnani, L., Moureaux, C., Sottocornola, M., and Wohlfahrt, G.: Contrasting response of European forest and grassland energy exchange to heatwaves, Nat. Geosci., 3, 722–727,, 2010. 

Trugman, A. T., Medvigy, D., Mankin, J. S., and Anderegg, W. R. L.: Soil moisture stress as a major driver of carbon cycle uncertainty, Geophys. Res. Lett., 45, 6495–6503,, 2018. 

Ukkola, A. M., De Kauwe, M. G., Pitman, A. J., Best, M. J., Abramowitz, G., Haverd, V., Decker, M., and Haughton, N.: Land surface models systematically overestimate the intensity, duration and magnitude of seasonal-scale evaporative droughts, Environ. Res. Lett., 11, 104012,, 2016a. 

Ukkola, A. M., Pitman, A. J., Decker, M., De Kauwe, M. G., Abramowitz, G., Kala, J., and Wang, Y.-P.: Modelling evapotranspiration during precipitation deficits: identifying critical processes in a land surface model, Hydrol. Earth Syst. Sci., 20, 2403–2419,, 2016b. 

Ukkola, A. M., Pitman, A. J., De Kauwe, M. G., Abramowitz, G., Herger, N., Evans, J. P., and Decker, M.: Evaluating CMIP5 model agreement for multiple drought metrics, J. Hydrometeorol., 19, 969–988,, 2018a. 

Ukkola, A. M., Pitman, A. J., Donat, M. G., De Kauwe, M. G., and Angélil, O.: Evaluating the contribution of land-atmosphere coupling to heat extremes in CMIP5 models, Geophys. Res. Lett., 45, 9003–9012,, 2018b. 

van Genuchten, M. T.: A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44, 892–898,, 1980. 

Van Looy, K., Bouma, J., Herbst, M., Koestel, J., Minasny, B., Mishra, U., Montzka, C., Nemes, A., Pachepsky, Y. A., Padarian, J., Schaap, M. G., Tóth, B., Verhoef, A., Vanderborght, J., Ploeg, M. J., Weihermüller, L., Zacharias, S., Zhang, Y., and Vereecken, H.: Pedotransfer functions in earth system science: challenges and perspectives, Rev. Geophys., 55, 1199–1256,, 2017. 

Vogel, M. M., Orth, R., Cheruy, F., Hagemann, S., Lorenz, R., van den Hurk, B. J. J. M., and Seneviratne, S. I.: Regional amplification of projected changes in extreme temperatures strongly controlled by soil moisture-temperature feedbacks, Geophys. Res. Lett., 44, 1511–1519,, 2017. 

Wang, Y. P., Kowalczyk, E., Leuning, R., Abramowitz, G., Raupach, M. R., Pak, B., van Gorsel, E., and Luhar, A.: Diagnosing errors in a land surface model (CABLE) in the time and frequency domains, J. Geophys. Res., 116, G01034,, 2011. 

Williams, M., Richardson, A. D., Reichstein, M., Stoy, P. C., Peylin, P., Verbeeck, H., Carvalhais, N., Jung, M., Hollinger, D. Y., Kattge, J., Leuning, R., Luo, Y., Tomelleri, E., Trudinger, C. M., and Wang, Y.-P.: Improving land surface models with FLUXNET data, Biogeosciences, 6, 1341–1359,, 2009. 

Wolf, A., Anderegg, W. R. L., and Pacala, S. W.: Optimal stomatal behavior with competition for water and risk of hydraulic impairment, P. Natl. Acad. Sci. USA, 113, E7222–E7230,, 2016. 

Xu, X., Medvigy, D., Powers, J. S., Becknell, J. M., and Guan, K.: Diversity in plant hydraulic traits explains seasonal and inter-annual variations of vegetation dynamics in seasonally dry tropical forests, New Phytol., 212, 80–95,, 2016. 

Yang, J.: MAESPA_EUCFACE_PARAM: Low sensitivity of gross primary production to elevated CO2 in a mature eucalypt woodland, Zenodo,, 2019. 

Yang, J., Duursma, R. A., De Kauwe, M. G., Kumarathunge, D., Jiang, M., Mahmud, K., Gimeno, T. E., Crous, K. Y., Ellsworth, D. S., Peters, J., Choat, B., Eamus, D., and Medlyn, B. E.: Incorporating non-stomatal limitation improves the performance of leaf and canopy models at high vapour pressure deficit, Tree Physiol., 39, 1961–1974,, 2019. 

Yang, J., Medlyn, B. E., De Kauwe, M. G., Duursma, R. A., Jiang, M., Kumarathunge, D., Crous, K. Y., Gimeno, T. E., Wujeska-Klause, A., and Ellsworth, D. S.: Low sensitivity of gross primary production to elevated CO2 in a mature eucalypt woodland, Biogeosciences, 17, 265–279,, 2020. 

Zhang, Y. and Schaap, M. G.: Weighted recalibration of the Rosetta pedotransfer model with improved estimates of hydraulic parameter distributions and summary statistics (Rosetta3), J. Hydrol., 547, 39–53,, 2017. 

Zhang, Y. and Schaap, M. G.: Estimation of saturated hydraulic conductivity with pedotransfer functions: A review, J. Hydrol., 575, 1011–1030,, 2019. 

Zhao, T. and Dai, A.: Uncertainties in historical changes and future projections of drought. Part II: model-simulated historical and future drought changes, Climatic Change, 144, 535–548,, 2017.  

Zhou, S., Duursma, R. A., Medlyn, B. E., Kelly, J. W. G., and Prentice, I. C.: How should we model plant responses to drought? An analysis of stomatal and non-stomatal responses to water stress, Agr. Forest Meteorol., 182–183, 204–214,, 2013.  

Zhou, S., Medlyn, B., Sabaté, S., Sperlich, D., Prentice, I. C., and Whitehead, D.: Short-term water stress impacts on stomatal, mesophyll and biochemical limitations to photosynthesis differ consistently among tree species from contrasting climates, Tree Physiol., 34, 1035–1046,, 2014. 

Short summary
Land surface model (LSM) is a critical tool to study land responses to droughts and heatwaves, but lacking comprehensive observations limited past model evaluations. Here we use a novel dataset at a water-limited site, evaluate a typical LSM with a range of competing model hypotheses widely used in LSMs and identify marked uncertainty due to the differing process assumptions. We show the extensive observations constrain model processes and allow better simulated land responses to these extremes.