Exploring how groundwater buffers the influence of heatwaves on vegetation function during multi-year droughts

The co-occurrence of droughts and heatwaves can have significant impacts on many socioeconomic and environmental systems. Groundwater has the potential to moderate the impact of droughts and heatwaves by moistening the soil and enabling vegetation to maintain higher evaporation, thereby cooling the canopy. We use the Community Atmosphere Biosphere Land Exchange (CABLE) land surface model, coupled to a groundwater scheme, to examine how groundwater influences ecosystems under conditions of co-occurring droughts and heatwaves. We focus specifically on south-east Australia for the period 2000–2019, when two significant droughts and multiple extreme heatwave events occurred. We found groundwater plays an important role in helping vegetation maintain transpiration, particularly in the first 1–2 years of a multi-year drought. Groundwater impedes gravity-driven drainage and moistens the root zone via capillary rise. These mechanisms reduced forest canopy temperatures by up to 5 ^∘C during individual heatwaves, particularly where the water table depth is shallow. The role of groundwater diminishes as the drought lengthens beyond 2 years and soil water reserves are depleted. Further, the lack of deep roots or stomatal closure caused by high vapour pressure deficit or high temperatures can reduce the additional transpiration induced by groundwater. The capacity of groundwater to moderate both water and heat stress on ecosystems during simultaneous droughts and heatwaves is not represented in most global climate models, suggesting that model projections may overestimate the risk of these events in the future.

1 Introduction

Droughts and heatwaves are important socioeconomic and environmental phenomena, impacting regional food production (Kim et al., 2019; Lesk et al., 2016), water resources (Leblanc et al., 2009; Orth and Destouni, 2018), and the resilience of ecosystems (Ibáñez et al., 2019; Ruehr et al., 2019; Sandi et al., 2020). When droughts and heatwaves co-occur (a “compound event”) the consequences can be particularly severe, reducing the terrestrial carbon sink (Ciais et al., 2005), potentially accelerating tree die-off (Allen et al., 2010, 2015; Birami et al., 2018), and setting conditions conducive to wildfires (Jyoteeshkumar reddy et al., 2021). One region experiencing severe coincident heatwaves and drought is Australia (Mitchell et al., 2014). Drought in Australia is associated with large-scale modes of variability, including the El Niño–Southern Oscillation and the Indian Ocean Dipole (van Dijk et al., 2013), and periods of below-average rainfall can extend for multiple years (Verdon-Kidd and Kiem, 2009). Heatwaves are commonly synoptically driven, associated with blocking events that can be sustained over many days (Perkins-Kirkpatrick et al., 2016; Perkins, 2015). Modes of variability and synoptic situations are important in setting up conditions conducive to drought and heatwave. However, once a heatwave or drought has become established, land–atmosphere interactions can intensify and prolong both heatwaves and droughts (Miralles et al., 2019), affect their intensity, and influence the risk of their co-occurrence (Mukherjee et al., 2020). The role of the land surface in amplifying or dampening heatwaves and droughts is associated with the partitioning of available energy between sensible and latent heat (Fischer et al., 2007; Hirsch et al., 2019) and is regulated by subsurface water availability (Teuling et al., 2013; Zhou et al., 2019). As soil moisture becomes more limiting, more of the available energy is converted into sensible heat, reducing evaporative cooling via latent heat. Changes in the surface turbulent energy fluxes influence the humidity in the boundary layer, the formation of clouds, incoming solar radiation, and the generation of rainfall (D'Odorico and Porporato, 2004; Seneviratne et al., 2010; Zhou et al., 2019). The sensible heat fluxes warm the boundary layer, leading to heat that can accumulate over several days and exacerbate heat extremes (Miralles et al., 2014), which can in turn increase the atmospheric demand for water and intensify drought (Miralles et al., 2019; Schumacher et al., 2019).

Vegetation access to groundwater has the potential to alter these land–atmosphere feedbacks by maintaining vegetation function during extended dry periods, supporting transpiration and moderating the impact of droughts and heatwaves (Marchionni et al., 2020; Miller et al., 2010). Where the water table is relatively shallow, capillarity may bring water from the groundwater towards the surface root zone, increasing plant water availability. Where the water table is deeper, phreatophytic vegetation with tap roots can directly access groundwater (Zencich et al., 2002). The presence of groundwater and the access to groundwater by vegetation are therefore likely to buffer vegetation drought and heatwave stress. For example, groundwater may help vegetation sustain transpiration and consequently cool plant canopies via evaporation. This is particularly critical during compound events where cessation of transpiration would increase the risk of impaired physiological function and the likelihood that plants would exceed thermal limits and risk mortality (Geange et al., 2021; O'sullivan et al., 2017; Sandi et al., 2020).

Quantifying the influence of groundwater on vegetation function has remained challenging as concurrent observations of groundwater dynamics, soil moisture, and energy and water fluxes are generally lacking over most of Australia and indeed many parts of the world. Land surface models (LSMs) provide an alternative tool for studying the interactions between groundwater, vegetation, and surface fluxes in the context of heatwaves and droughts (Gilbert et al., 2017; Martinez et al., 2016a; Maxwell et al., 2011; Shrestha et al., 2014). However, there has been very little work focused on the influence of groundwater on droughts and heatwaves occurring at the same time (Keune et al., 2016; Zipper et al., 2019). Our key goal in this paper is therefore to examine the timescales and extent to which vegetation utilises groundwater during drought and heatwaves and determine the degree to which groundwater can mitigate the impacts of compound extremes. We focus on droughts and heatwaves occurring over south-eastern (SE) Australia during 2000–2019 using the Community Atmosphere Biosphere Land Exchange (CABLE) LSM. SE Australia is an ideal case study since its forest and woodland ecosystems are known to be dependent on groundwater (Eamus and Froend, 2006; Kuginis et al., 2016; Zencich et al., 2002), and it has experienced two multi-year droughts and record-breaking heatwaves over the last 2 decades. By examining the role of groundwater in influencing droughts and heatwaves and by understanding how well CABLE can capture the relevant processes, we aim to build confidence in the simulations of land–atmosphere interactions for future droughts and heatwaves.

2 Methods

2.1 Study area

The climate over SE Australia varies from humid temperate near the coast to semi-arid in the interior. In the last 20 years, SE Australia experienced the 9-year Millennium Drought during 2001–2009 (van Dijk et al., 2013), where rainfall dropped from a climatological average (1970–1999) of 542 to 449 mm yr⁻¹, and a 3-year intense recent drought during 2017–2019 where rainfall dropped to 354 mm yr⁻¹ (Fig. S1). It has also suffered record-breaking summer heatwaves in 2009, 2013, 2017, and 2019 (Bureau of Meteorology, 2013, 2017, 2019; National Climate Centre, 2009). Here we investigate groundwater interactions during the period 2000–2019, focusing on the Millennium Drought (MD, 2001–2009) and the recent drought (RD, 2017–2019).

2.2 Overview of CABLE

CABLE is a process-based LSM that simulates the interactions between climate, plant physiology, and hydrology (Wang et al., 2011). Above ground, CABLE simulates the exchange of carbon, energy, and water fluxes, using a single-layer, two-leaf (sunlit and shaded) canopy model (Wang and Leuning, 1998), with a treatment of within-canopy turbulence (Raupach, 1994; Raupach et al., 1997). CABLE includes a six-layer soil model (down to 4.6 m) with soil hydraulic and thermal characteristics dependent on the soil type and soil moisture content. CABLE has been extensively evaluated (e.g. Abramowitz et al., 2008; Wang et al., 2011; Zhang et al., 2013) and benchmarked (Abramowitz, 2012; Best et al., 2015) at global and regional scales.

Here we adopt a version of CABLE (Decker, 2015; Decker et al., 2017) which includes a dynamic groundwater component with aquifer water storage. This version, CABLE-GW, has been previously evaluated by Decker (2015), Ukkola et al. (2016b), and Mu et al. (2021) and shown to perform well for simulating water fluxes. CABLE code is freely available upon registration (https://trac.nci.org.au/trac/cable/wiki, last access: 29 August 2021); here we use CABLE SVN revision 7765.

2.3 Hydrology in CABLE-GW

The hydrology scheme in CABLE-GW solves the vertical redistribution of soil water via a modified Richards equation (Zeng and Decker, 2009):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}=-{\displaystyle \frac{\partial }{\partial z}}K{\displaystyle \frac{\partial }{\partial z}}\left(\mathrm{\Psi }-{\mathrm{\Psi }}_E\right)-F_{\mathrm{soil}},\end{array}$$

where θ is the volumetric water content of the soil (mm³ mm⁻³), K is the hydraulic conductivity (mm s⁻¹), z is the soil depth (mm), Ψ and Ψ_E are the soil matric potential (mm) and the equilibrium soil matric potential (mm), and F_soil is a sink term related to subsurface runoff and transpiration (s⁻¹) (Zeng and Decker, 2009; Decker, 2015). To simulate groundwater dynamics, an unconfined aquifer is added to the bottom of the soil column with a simple water balance model:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}{W}_{\mathrm{aq}}}{\mathrm{d}t}}=q_{\mathrm{re}}-q_{\mathrm{aq},\mathrm{sub}},\end{array}$$

where W_aq is the mass of water in the aquifer (mm), q_aq,sub is the subsurface runoff in the aquifer (mm s⁻¹), and q_re is the water flux between the aquifer and the bottom soil layer (mm s⁻¹) computed by the modified Darcy's law:

$$\begin{array}{}\text{(3)}& q_{\mathrm{re}}={\displaystyle \frac{(K_{\mathrm{aq}}+K_{\mathrm{bot}})}{\mathrm{2}}}{\displaystyle \frac{\left({\mathrm{\Psi }}_{\mathrm{aq}}-{\mathrm{\Psi }}_n\right)-\left({\mathrm{\Psi }}_{E,\mathrm{aq}}-{\mathrm{\Psi }}_{E,n}\right)}{z_{\mathrm{wtd}}-z_n}},\end{array}$$

where K_aq and K_bot are the hydraulic conductivity in the aquifer and in the bottom soil layer (mm s⁻¹), Ψ_aq and Ψ_E,aq are the soil matric potentials for the aquifer (mm), and Ψ_n and Ψ_E,n are the soil matric potentials for the bottom soil layer (mm). z_wtd and z_n are the depth of the water table (mm) and the lowest soil layer (mm), respectively. The positive q_re refers to the downward water flow from soil column to aquifer (i.e. vertical drainage, D_r), and the negative q_re is the upward water movement from the aquifer to the soil column (i.e. recharge, Q_rec). CABLE-GW assumes that the groundwater aquifer sits above impermeable bedrock, giving a bottom boundary condition of

$$\begin{array}{}\text{(4)}& q_{\mathrm{out}}=\mathrm{0}.\end{array}$$

CABLE-GW computes the subsurface runoff (q_sub, mm s⁻¹) using

$$\begin{array}{}\text{(5)}& q_{\mathrm{sub}}=\sin \stackrel{\mathrm{‾}}{{\displaystyle \frac{d_z}{d_l}}}{\widehat{q}}_{\mathrm{sub}}{e}^{-\frac{z_{\mathrm{wtd}}}{f_{\mathrm{p}}}},\end{array}$$

where $\stackrel{\mathrm{‾}}{\frac{d_z}{d_l}}$ is the mean subgrid-scale slope, ${\widehat{q}}_{\mathrm{sub}}$ is the maximum rate of subsurface drainage (mm s⁻¹), and f_p is a tunable parameter. q_sub is generated from the aquifer and the saturated deep soil layers (below the third soil layer).

2.4 Experiment design

To explore how groundwater influences droughts and heatwaves, we designed two experiments, with and without groundwater dynamics, driven by the same 3 h time-evolving meteorology forcing and the same land surface properties (see Sect. 2.5 for datasets) for the period 1970–2019 at a 0.05^∘ spatial resolution with a 3 h time step. To correct a tendency for high soil evaporation, we implemented a parameterisation of soil evaporation resistance that has previously been shown to improve the model (Decker et al., 2017; Mu et al., 2021).

2.4.1 Groundwater experiment (GW)

This simulation uses the default CABLE-GW model, which includes the unconfined aquifer to hold the groundwater storage and simulates the water flux between the bottom soil layer and the aquifer. We first ran the default CABLE-GW with fixed CO₂ concentrations at 1969 levels for 90 years by looping the time-evolving meteorology forcing over 1970–1999. At the end of the 90-year spin-up, moisture in both the soil column and the groundwater aquifer reached an effective equilibrium when averaged over the study area. We then ran the model from 1970 to 2019 with time varying CO₂. We omit the first 30 years of this period and analyse the period 2000–2019 to allow for further equilibrium with the time-evolving CO₂.

2.4.2 Free drainage experiment (FD)

Many LSMs, including those used in the Coupled Model Intercomparison Project 5 (CMIP5), still use a free drainage assumption and neglect the parameterisation of the unconfined aquifer. To test the impact of this assumption we decoupled the aquifer from the bottom soil layer and thus removed the influence of groundwater dynamics (experiment FD). In FD, at the interface between the bottom soil layer and the aquifer, soil water can only move downwards as vertical drainage at the rate defined by the bottom soil layer's hydraulic conductivity:

$$\begin{array}{}\text{(6)}& q_{\mathrm{re}}=K_{\mathrm{bot}}.\end{array}$$

This vertical drainage is added to the subsurface runoff flux:

$$\begin{array}{}\text{(7)}& q_{\mathrm{sub}}=q_{\mathrm{sub}}+q_{\mathrm{re}}.\end{array}$$

The simulated water table depth (WTD, i.e. z_wtd) in CABLE-GW affects the water potential gradient between the soil layers via Ψ_E (Zeng and Decker, 2009) and impacts q_sub (Eq. 5). However, in FD, decoupling the soil column from the aquifer and adding vertical drainage directly to subsurface runoff causes an artificial and unrealistic decline in WTD. To solve this problem, we assume a fixed WTD in the FD simulations at 10 m in order to remove this artefact from the simulation of Ψ_E and q_sub.

The FD simulations are initialised from the near-equilibrated state at the end of the 90-year spin-up used in GW. The period 1970–2019 is then simulated using varying CO₂ and the last 20 years are used for analysis.

2.4.3 Deep root experiment (DR)

The parameterisation of roots, including the prescription of root parameters in LSMs, is very uncertain (Arora and Boer, 2003; Drewniak, 2019) and LSMs commonly employ root distributions that are too shallow (Wang and Dickinson, 2012). The vertical distribution of roots influences the degree to which plants can utilise groundwater, and potentially the role groundwater plays in influencing droughts and heatwaves. To explore the uncertainty associated with root distribution, we added a “deep root” (DR) experiment by increasing the effective rooting depth in CABLE for tree areas. In common with many LSMs, CABLE-GW defines the root distribution following Gale and Grigal (1987):

$$\begin{array}{}\text{(8)}& f_{\mathrm{root}}=\mathrm{1}-{\mathit{\beta }}_{\mathrm{root}}^z,\end{array}$$

where f_root is the cumulative root fraction (between 0 and 1) from the soil surface to depth z (cm), and β_root is a fitted parameter specified for each plant functional type (PFT) (Jackson et al., 1996). In CABLE, the tree areas in our study region are simulated as evergreen broadleaf PFT (Fig. S2a) with a β_root= 0.962, implying that only 8 % of the simulated roots are located below the soil depth of 0.64 m (Fig. S2b). However, field observations (Canadell et al., 1996; Eberbach and Burrows, 2006; Fan et al., 2017; Griffith et al., 2008) suggest that the local trees tend to have a far deeper root system, possibly to help cope with the high climate variability. We therefore increased β_root for the evergreen broadleaf PFT to 0.99, which assumes 56 % of roots are located at depths below 0.64 m and 21 % of roots below 1.7 m (Fig. S2b). This enables the roots to extract larger quantities of deep soil water moisture, which is more strongly influenced by groundwater.

Given that we lack the detailed observations to set root distributions across SE Australia, we undertake the DR experiment as a simple sensitivity study. We only run this experiment during January 2019, when the record-breaking heatwaves compound with the severe recent drought. The DR experiment uses identical meteorology forcing and land surface properties as GW and FD and is initialised by the state of the land surface on the 31 December 2018 from the GW experiment.

2.5 Datasets

Our simulations are driven by the atmospheric forcing from the Australian Water Availability Project (AWAP), which provides daily gridded data covering Australia at 0.05^∘ spatial resolution (Jones et al., 2009). This dataset has been widely used to force LSMs for analysing the water and carbon balances in Australia (Haverd et al., 2013; De Kauwe et al., 2020; Raupach et al., 2013; Trudinger et al., 2016). The AWAP forcing data include observed fields of precipitation, solar radiation, minimum and maximum daily temperatures, and vapour pressure at 09:00 and 15:00. Since AWAP forcing does not include wind and air pressure we adopted the near-surface wind speed data from McVicar et al. (2008) and assume a fixed air pressure of 1000 hPa. Due to missing observations before 1990, the solar radiation input for 1970–1989 was built from the 1990–1999 daily climatology. Similarly, wind speeds for 1970–1974 are built from the 30-year climatology from 1975 to 2004. We translated the daily data into 3-hourly resolution using a weather generator (Haverd et al., 2013).

The land surface properties for our simulations are prescribed based on observational datasets. Land cover type is derived from the National Dynamic Land Cover Data of Australia (DLCD) (https://www.ga.gov.au/scientific-topics/earth-obs/accessing-satellite-imagery/landcover, last access: 29 August 2021). We classify DLCD's land cover types to five CABLE PFTs: crop (irrigated/rainfed crop, pasture, and sugar DLCD classes), broadleaf evergreen forest (closed/open/scattered/sparse tree), shrub (closed/open/scattered/sparse shrubs and open/scattered/sparse chenopod shrubs), grassland (open/sparse herbaceous), and barren land (bare areas). We then resample the DLCD dataset from the 250 m resolution to the 0.05^∘ resolution. The leaf area index (LAI) in CABLE is prescribed using a monthly climatology derived from the Copernicus Global Land Service product (https://land.copernicus.eu/global/products/lai, last access: 29 August 2021). The climatology was constructed by first creating a monthly time series by taking the maximum of the 10 daily time steps each month and then calculating a climatology from the monthly data over the period 1999–2017. The LAI data were resampled from the original 1 km resolution to the 0.05^∘ resolution following De Kauwe et al. (2020). Soil parameters are derived from the soil texture information (sand, clay, and silt fraction) from SoilGrids (Hengl et al., 2017) via the pedotransfer functions in Cosby et al. (1984) and resampled from 250 m to 0.05^∘ resolution.

To evaluate the model simulations, we use monthly total water storage anomaly (TWSA) at 0.5^∘ spatial resolution from the Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow On products (Landerer et al., 2020; Watkins et al., 2015; Wiese et al., 2016, 2018). The RLM06M release is used for February 2002–June 2017 and for June 2018–December 2019. We also use total land evapotranspiration from the 2000–2018 monthly Derived Optimal Linear Combination Evapotranspiration (DOLCE version 2, Hobeichi et al., 2021a) at 0.25^∘ resolution, as well as the 2000–2019 daily Global Land Evaporation Amsterdam Model (GLEAM version 3.5, https://www.gleam.eu/; Martens et al., 2017; Miralles et al., 2011) at 0.5^∘ spatial resolution. For daytime land surface temperature (LST) we use the Moderate Resolution Imaging Spectroradiometer (MODIS) datasets from Terra and Aqua satellites (products MOD11A1 and MYD11A1, Wan and Li, 1997; Wan et al., 2015a, b) at 1 km spatial resolution. We only consider pixels and time steps identified as good quality (QC flags 0). Only the daytime LST values are used due to the lack of good-quality nighttime LST data. The Terra overpass occurs at 10:00 and Aqua at 14:00 local time. To analyse the compound events in January 2019, we linearly interpolate the 3-hourly model outputs to 14:00 to match the overpass time of the Aqua LST. The GRACE, GLEAM, and MODIS datasets were resampled to the AWAP resolution using bilinear interpolation.

To evaluate model performance during heatwaves, we identify heatwave events using the excess heat factor index (EHF, Nairn and Fawcett, 2014). EHF is calculated using the daily AWAP maximum temperature, as the product of the difference of the previous 3 d mean to the 90th percentile of the 1970–1999 climatology and the difference of the previous 3 d mean to the preceding 30 d mean. A heatwave occurs when the EHF index is greater than 0 for at least 3 consecutive days. We only focus on summer heatwaves occurring between December and February of the following year.

3 Results

3.1 Simulations for the Millennium Drought and the recent drought

Previous studies have shown that simulations by LSMs diverge as the soil dries (Ukkola et al., 2016a), associated with systematic biases in evaporative fluxes and soil moisture states in the models (Mu et al., 2021; Swenson and Lawrence, 2014; Trugman et al., 2018). We therefore first evaluate how well CABLE-GW captures the evolution of terrestrial water variability during two recent major droughts.

Figure 1a shows the total water storage anomaly during 2000–2019 observed by GRACE and simulated in GW and FD. Both GW and FD accurately capture the interannual variability in total water storage for SE Australia (r= 0.96 in GW and 0.90 in FD). Both model configurations simulate a decline in TWSA through the first drought period (up to 2009, see Fig. S1), the rapid increase in TWSA from 2010 associated with higher rainfall, a decline from around 2012 due to the re-emergence of drought conditions, and the rapid decline during the recent drought after conditions had eased in 2016 (Fig. S1). FD underestimates the magnitude of monthly TWSA variance (standard deviation, SD = 37.18 mm) compared to GRACE (47.74 mm) or GW (47.67 mm), particularly during the wetter periods (2000, 2011–2016) and the first ∼2 years of the droughts (2001–2002, 2017–2018) (Fig. 1a). This underestimation in FD compared to GW is linked to the lack of aquifer water storage in the FD simulations, which provides a reservoir of water that changes slowly and has a memory of previous wet/dry climate conditions (Fig. 1a).

Figure 1 (a) Total water storage anomaly (TWSA) during 2000–2019 and (b) accumulated P–E for the two droughts over SE Australia. In panel (a), observations from GRACE are shown in black, the GW simulation in green, FD in orange, and the aquifer water storage anomaly in GW in blue. The shading in panel (a) highlights the two drought periods. The left top corner of panel (a) displays the correlation (r) between GRACE and GW/FD, as well as the standard deviation (SD, mm) of GRACE, GW, and FD over the periods when GRACE and the simulations coincide. Panel (b) shows the accumulated P–E for two periods; the dark lines show the 2001–2009 Millennium Drought (MD) and the light lines show the 2017–2019 recent drought (RD). The correlation (r) between the P and E is shown in the legend of panel (b).

Figure 1b shows the accumulated precipitation (P) minus evapotranspiration (E) over the two drought periods. GW increases the evapotranspiration relative to FD such that the accumulated P–E decreases from about 786 to 455 mm during the Millennium Drought, which is within the range of DOLCE (460 mm) and GLEAM (97 mm) estimates. A similar result, although over a much shorter period, is also apparent for the recent drought (Fig. 1b). The lower P–E in GW suggests that the presence of groundwater storage can alleviate the vegetation water stress during droughts and reduces the reliance of E on P, indicated by a small reduction in the correlation (r) between E and P from 0.28 in FD to 0.24 in GW for MD and a reduction from 0.42 to 0.37 for RD (Fig. 1b). Although the evapotranspiration products display some differences, the GW simulations are closer overall to both the DOLCE and the GLEAM observationally constrained estimates. The better match of GW than FD to the two evapotranspiration products implies that adding groundwater improves the simulations during droughts, whilst the remaining mismatch would tend to suggest further biases in simulated evapotranspiration arising from multiple sources (e.g. a mismatch in leaf area index or contributions from the understorey). The difference in E is also demonstrated spatially in Fig. S3. During the Millennium Drought, the GW simulations show a clear improvement over FD in two aspects. GW shows smaller biases in E along the coast where FD underestimates E strongly (Fig. S3b–c). The areas where E is underestimated are also smaller in extent in GW, suggesting that GW overall reduces the dry bias. The magnitude of the bias in GW reaches around 300 mm over small areas of SE Australia, while in the FD simulations biases are larger, reaching 400 mm over a larger area. Plant photosynthesis assimilation rates are associated with transpiration via stomata conductance. Figure S4 presents the spatial maps of gross primary productivity (GPP) during the two droughts. GW simulations increase carbon uptake by 50–300 g C yr⁻¹ along the coasts (Fig. S4c, f). However, since CABLE uses a prescribed LAI and does not simulate any feedback between water availability and plant growth (e.g. defoliation) and its impact on GPP, we only focus on how GW influences evapotranspiration and the surface energy balance in the subsequent sections.

Overall, Figs. 1 and S3 indicates that representing groundwater in the model improves the simulation of the interannual variability in the terrestrial water cycle and storage, particularly during droughts.

3.2 The role of groundwater in sustaining evapotranspiration during droughts

We next explore the mechanisms by which including groundwater modifies the simulation of evapotranspiration. Figure 2 displays the overall influence of groundwater on water fluxes during the recent drought. GW simulates 50–200 mm yr⁻¹ more E over coastal regions where there is high tree cover (Fig. 2a; see Fig. S2a for land cover). Adding groundwater also increases E in most other regions, although the impact is negligible in many inland and non-forested regions (i.e. west of 145^∘ E). We identified a clear connection between E (Fig. 2a) and the simulated WTD in the GW simulations (Fig. S5). GW simulates 110 mm yr⁻¹ more E when the WTD is shallower than 5 m deep, 22 mm yr⁻¹ when the WTD is 5–10 m deep, but only 3 mm yr⁻¹ more when the WTD is below 10 m. Higher transpiration (E_t) in GW explains 78 % of the evapotranspiration difference between GW and FD where WTD is shallower than 5 m (Fig. 2b). This is confirmed by the change in the soil evaporation (ΔE_s) (Fig. 2c) where adding groundwater increases E_s by negligible amounts over most of SE Australia but by up to 25 mm yr⁻¹ in regions underlain by shallow groundwater (Fig. S5), which is consistent with field observations that indicate that E_s can be substantial under conditions of a very shallow water table (Thorburn et al., 1992). In the very shallow WTD areas, the excess E_s in GW results from the capillary rise of moisture from the shallow groundwater to the surface.

Figure 2 The overall influence of groundwater during the recent drought. Panels (a)–(e) are the difference (GW–FD) in evapotranspiration (ΔE), transpiration (ΔE_t), soil evaporation (ΔE_s), vertical drainage (ΔD_r) and recharge from the aquifer to soil column (ΔQ_rec), respectively. In the bottom right of panels (a)–(e), the average of each variable over selected water table depths (WTDs) is provided. Panel (f) is the maximum nighttime water stress factor difference (Δβ) between 03:00 (i.e. predawn when the soil is relatively moist following capillary lift overnight) and 21:00 the previous day. We only include rainless nights in January 2019 to calculate Δβ to remove any influence of overnight rainfall.

A significant factor in explaining how groundwater influences E is through changes in vertical drainage and recharge from the aquifer to the soil column. Figure 2d shows that the vertical drainage (D_r) both increases and decreases depending on the location. The addition of groundwater reduces vertical drainage by 74 mm yr⁻¹ where WTD is shallower than 5 m. In some regions, the drainage increases with the inclusion of groundwater by up to 100 mm yr⁻¹, especially in the areas where WTD is ∼5 m. This is associated with the WTD being slightly below the bottom of the soil column (4.6 m). When the groundwater aquifer is nearly full in GW, the wetter soil in the bottom layer leads to a much higher hydraulic conductivity in GW than in FD, leading to higher vertical drainage in GW and a positive ΔD_r. Inland, where the WTD tends to be much deeper there is no significant difference in D_r between GW and FD.

Figure 2e shows the difference in recharge into the upper soil column (ΔQ_rec) between GW and FD. The recharge from the aquifer into the bottom soil layer provides 17 mm yr⁻¹ of extra moisture in GW, where the WTD is between 5–10 m, and 10 mm yr⁻¹ where the WTD is < 10 m, partially explaining the changes in E and E_t in areas with a deep WTD. However, there is no significant ΔQ_rec in regions with a shallow WTD (∼5 mm yr⁻¹), suggesting the influence of groundwater is mainly via reduced drainage in these locations. Recharge from the aquifer to the soil column can only occur when WTD is below the soil column (bottom boundary at 4.6 m depth). If WTD is shallow and within the soil column, the interface is saturated and no recharge from the aquifer to the soil column can occur and water only moves downwards by gravity.

The combined impact of reduced drainage in GW (Fig. 2d) and recharge from the aquifer into the root zone (Fig. 2e) is an increased water potential gradient between the drier top soil layers and the wetter deep soil layers, encouraging overnight capillary rise. Taking the hot and dry January 2019 as an example, when the compound events occurred, Fig. 2f shows the maximum water stress factor difference (Δβ) overnight (between 21:00 and 03:00, i.e. predawn when soil is relatively moist following capillary lift overnight). We only consider rainless nights to exclude the impact of drainage induced by precipitation. The water stress factor (β) is based on the root distribution and moisture availability in each soil layer and represents the soil water stress on transpiration as water becomes limiting. Figure 2f implies that while the redistribution of moisture is small overall, in some locations it can reduce moisture stress by up to 4 %–6 %.

3.3 The impact of groundwater during heatwaves

We next explore whether the higher available moisture due to the inclusion of groundwater enables the canopy to cool itself via evapotranspiration during heatwaves by examining the temperature difference between the simulated canopy temperature (T_canopy, ^∘C) and the forced air temperature (T_air, ^∘C). We focus on the forested regions (Fig. S2a) as the role of groundwater in enhancing plant water availability was shown to be largest in these regions (Fig. 2).

Figure 3a shows the average T_canopy−T_air (ΔT, ^∘C) over the forested regions for summer heatwaves from the GW and FD simulations, with the grey line indicating the median ΔT difference. During heatwaves, the inclusion of groundwater moistens the soil and supports higher transpiration, cooling the canopy and reducing ΔT relative to FD by up to 0.76 ^∘C (e.g. summer heatwaves in 2013). As the drought lengthens in time, the depletion of moisture gradually reduces this effect, from an average reduction of 0.52 ^∘C of the first 3 years to 0.16 ^∘C of the last 3 years in the Millennium Drought (Fig. 3a). The impact of groundwater is clear in the evaporative fraction (Fig. 3b) where in periods of higher rainfall (e.g. 2010–2011; Fig. S1) and at the beginning of a drought (2001, 2017) the evaporative fraction (EF) is higher (0.03 to 0.18). This implies that more of the available energy is exchanged with the atmosphere in the form of latent rather than sensible heat. However, the strength of the cooling effect decreases as the droughts extend and the transpiration difference (ΔE_t, mm d⁻¹) diminishes quickly (Fig. 3c) because the vegetation becomes increasingly water-stressed (Fig. 3d) which consequently limits transpiration. For all variables (ΔT, EF, E_t, and β), the difference between GW and FD is greatest during the wetter periods (e.g. 2013) and in the first 1–2 years of the multi-year drought (2001–2002 for the Millennium Drought or 2017–2018 for the recent drought). After the drought becomes well established, the FD and GW simulations converge as depleting soil moisture reservoirs reduce the impact of groundwater on canopy cooling and evaporative fluxes.

Figure 3 Groundwater-induced differences in (a) T_canopy−T_air (ΔT), (b) evaporative fraction (EF), (c) transpiration (E_t), and (d) water stress factor (β) during 2000–2019 summer heatwaves over forested areas (the green region in Fig. S2a). The left y axis is the scale for boxes. The blue boxes refer to the GW experiment and the red boxes to FD. For each box, the middle line is the median, the upper border is the 75th percentile, and the lower border is 25th percentile. The right y axis is the scale for the grey lines which display the difference in the medians (GW–FD). The shadings highlight the two drought periods.

Figure 4 Land response to heatwaves during the recent drought. Panels (a)–(c) are the mean T_canopy−T_air (ΔT), evaporative fraction (EF), and soil water stress factor (β) in GW, respectively, during 2017–2019 summer heatwaves. Panels (d)–(f) are the difference (GW–FD) of T_canopy, EF, and β. In the bottom right of each plot, the average of each variable over selected water table depths (WTDs) is provided. Note that the colour bar is switched between (d) and (e)–(f).

Figure 4a shows the spatial map of ΔT simulated in GW during heatwaves in the 2017–2019 drought. It indicates that both land cover type (Fig. S2a) and WTD (Fig. S5) contribute to the ΔT pattern. The evaporative cooling via transpiration is stronger over the forested areas compared to cropland or grassland and stronger in the regions with a wetter soil associated with a shallower WTD. However, EF is mainly determined by WTD (compare Figs. 4b and S5). Inland, where the WTD is deeper and the soil is drier, most of the net radiation absorbed by the land surface is partitioned into sensible rather than latent heat (Fig. 4b). However, in the coastal regions with a shallow WTD, the wetter soil reduces the water stress (Fig. 4c), enables a higher EF (Fig. 4b), and alleviates heat stress on the leaves (Fig. 4a). Along the coast where WTD is shallow, GW simulates a cooler canopy temperature due to the higher evaporative cooling (Fig. 4e), which is the consequence of a lower soil water stress (Fig. 4f) linked to the influence of groundwater (Fig. S5).

Figure 5 shows the density scatter plot of ΔT versus WTD in SE Australia forested areas during heatwaves in 2000–2019. A shallow WTD moderates the temperature difference between the canopy and the ambient air during heatwaves leading to a smaller temperature difference. Meanwhile, as the WTD increases, due to the limited rooting depth in the model, the ability of the groundwater to support transpiration and offset the impact of high air temperatures is reduced. Figure 5 shows a large amount of variations but nonetheless implies a threshold of ∼6 m, whereafter there is a decoupling and little influence from groundwater during heatwaves. However, the absolute value of the threshold is likely CABLE-specific and associated with the assumption of a 4.6 m soil depth, which also sets the maximum rooting depth (roots can only extend to the bottom of the soil and cannot directly access the groundwater aquifer in CABLE). The CABLE soil depth comes from observational evidence of most roots being situated within the top 4.6 m (Canadell et al. 1996). Since the model assumes no roots exist in the groundwater aquifer, when the water table is below this depth, the water fluxes become largely uncoupled between the soil column and the groundwater aquifer, leading to a negligible impact of GW below ∼6 m depth.

Figure 5 A density scatter plot of T_canopy−T_air (ΔT) versus water table depth (WTD) in GW simulations over forested areas in all heatwaves during 2000–2019. Every tree pixel on each heatwave day accounts for one record, and the darker colours show higher recorded densities.

3.4 The impact of groundwater during the drought and heatwave compound events

To examine the influence of groundwater on heatwaves occurring simultaneously with drought, we focus on a case study of the record-breaking heatwaves in January 2019, which is the hottest month on record for the study region (Bureau of Meteorology, 2019). The unprecedented prolonged heatwave period started in early December 2018 and continued through January 2019 with three peaks. We select 2 d (15 and 25 January 2019), when heatwaves spread across the study region, from the second and third heatwave phases (Fig. S6).

We evaluate CABLE T_canopy against MODIS LST observations, concentrating on forested areas where MODIS LST should more closely reflect vegetation canopy temperatures, but note that this comparison is not direct as the satellite estimate will contain contributions from the understorey and soil. Figure 6a–b show the good-quality MODIS LST minus T_air at 14:00 (ΔT_{MOD_14:00}) over forested regions on the 15 and 25 January 2019, and Fig. 6c–d display the matching GW-simulated ΔT at 14:00 (ΔT_{GW_14:00}). Overall, ΔT_{GW_14:00} increases from the coast to the interior in both heatwaves, consistent with the ΔT_{MOD_14:00} pattern in both heatwaves, although ΔT_{GW_14:00} appears to be biased high relative to ΔT_{MOD_14:00} along the coastal forests (Fig. S7a–b).

Figure 6 The simulation of two heatwaves on 15 (left column) and 25 January 2019 (right column). The first row shows the difference between MODIS land surface temperature (LST) and T_air at 14:00 (ΔT_{MOD_14:00}) (only forested areas with good LST quality data are displayed). The second row is the GW simulation of ΔT at 14:00 (ΔT_{GW_14:00}). The third row is the difference of ΔT at 14:00 between GW and FD simulations ($\mathrm{\Delta }{T}_{\mathrm{GW}\mathrm{\_}\mathrm{14}:\mathrm{00}}-\mathrm{\Delta }{T}_{\mathrm{FD}\mathrm{\_}\mathrm{14}:\mathrm{00}}$). The last row is the same as the third row but for the difference between the DR and GW simulations ($\mathrm{\Delta }{T}_{\mathrm{DR}\mathrm{\_}\mathrm{14}:\mathrm{00}}-\mathrm{\Delta }{T}_{\mathrm{GW}\mathrm{\_}\mathrm{14}:\mathrm{00}}$). Note that the comparison between GW/FD/DR and MODIS LST is shown in Fig. S7.

Figure 6e–f show the ΔT_14:00 difference between GW and FD. Access to groundwater can reduce canopy temperature by up to 5 ^∘C, in particular where the WTD is shallow. While reductions of 5 ^∘C are clearly limited in spatial extent, the overall pattern of cooling is quite widespread and coincident with the groundwater-induced E_t increase (Fig. S8a–b), implying a reduction in heat stress along coastal regions with a shallow WTD during heatwaves. Generally, GW matches MODIS LST better than FD despite the bias in both simulations (compare Figs. S7a–b and S7c–d). Nevertheless, the temperature reduction between GW and FD is still modest (< 1 ^∘C) for most of the forested regions. This may be related to the shallow root distribution assumed in many LSMs, which prevents roots from directly accessing the moisture stored in the deeper soil (note, CABLE assumes 92 % of a forest's roots are in the top 0.64 m, Fig. S2b). To examine this possibility, we performed the deep root (DR) sensitivity experiment which prescribed more roots in the deeper soil for forests (56 % below 0.64 m depth). Figure 6g–h illustrate the difference between ΔT_14:00 in DR and ΔT_{GW_14:00}. By enabling access to moisture in the deeper soil, the LSM simulates further cooling by 0.5–5 ^∘C across the forests associated with an E_t increase of 25–250 W m⁻² (Fig. S8c–d). The prescribed deeper roots also lead to an overall better simulation of ΔT at 20:00 relative to the MODIS LST (Figs. S7e–f vs. Fig. S7a–b).

Figure 7 shows the diurnal cycles of ΔT for the two selected regions (red boxes in Fig. 6) compared with the MODIS LST estimates. The region highlighted for the 15 January (Fig. 7a) has a WTD of 4–7 m, while the region highlighted for the 25 January (Fig. 7b) has a WTD < 4 m (Fig. S5). In both regions, the simulated ΔT is highest in FD, lower in GW and lowest in DR. Where the WTD is 4–7 m (Fig. 7a), the three simulated ΔT are slightly lower than ΔT calculated by MODIS LST (red squares). However, in the shallower WTD region (Fig. 7b), the simulated ΔT between experiments is more dispersed across experiments and exceeds the MODIS ΔT at both time points, implying that neglecting groundwater dynamics and deep roots is more likely to cause an overestimation of heat stress in the shallower WTD region. The shallower WTD region (Fig. 7b) tends to have a high LAI coverage, implying that the MODIS LST represents a good approximation of the canopy temperature over this region. Consequently, the lower MODIS ΔT implies that CABLE is likely underestimating transpiration, leading to an overestimation of ΔT in all three simulations.

Figure 7 Diurnal cycle of ΔT on 15 (left column) and 25 January 2019 (right column) over the selected regions shown in Fig. 6. The shadings show the uncertainty in every simulation defined as one standard deviation (SD) among the selected pixels. The red dots are MODIS LST minus T_air with the uncertainty shown by the red error bars. For both regions, only pixels available in MODIS are shown.

3.5 Constraints on groundwater mediation during the compound events

We finally probe the reasons for the apparent contradiction between the large impact of groundwater on E during drought (Fig. 2a) but a smaller impact on ΔT during the compound events (Fig. 7). Figure 8 shows three factors (β, vapour pressure deficit (D), and T_air) that constrain the impact of groundwater on ΔT in CABLE during heatwaves in January 2019. Figure 8a shows the difference in ΔT between GW and FD as a function of Δβ, suggesting that the inclusion of groundwater has a large impact on ΔT when there is a coincidental and large difference in β between the GW and FD simulations. Figure 8b indicates a clear threshold at D= 3 kPa where GW and FD converge, while Fig. 8c shows a convergence threshold when T_air exceeds 32 ^∘C. Above these two thresholds, access to groundwater seemingly becomes less important in mitigating plant heat stress. There are two mechanisms in CABLE that explain this behaviour. First, as D increases, CABLE predicts that stomata begin to close following a square root dependence (De Kauwe et al., 2015; Medlyn et al., 2011). Second, as T_air increases, photosynthesis becomes inhibited as the temperature exceeds the optimum for photosynthesis. In both instances, evaporative cooling is reduced, regardless of the root zone moisture state dictated by groundwater access. That is to say, access to groundwater has limited capacity to directly mediate the heat stress on plants during a compound event when the air is very dry or very hot.

Figure 8 Density scatter plots showing the three factors that influence the difference in T_canopy−T_air between GW and FD (ΔT, expressed as GW–FD difference). Panel (a) is ΔT difference against the β difference (GW–FD) (Δβ), (b) is ΔT difference against vapour pressure deficit (D), and (c) is ΔT difference against T_air. Each point corresponds to a tree pixel on a heatwave day in January 2019. The darker colours illustrate where the records are denser. The correlation (r) between the x and y axes is shown in the bottom left of each panel.

4 Discussion

In the absence of direct measurements, we used the CABLE-GW LSM, constrained by satellite observations to investigate how groundwater influences ecosystems under conditions of co-occurring droughts and heatwaves. We found that the influence of groundwater was most important during the wetter periods and the first ∼ 2 years of a multi-year drought (∼2001–2002 and 2017–2018; Figs. 1 and 3). This primarily occurred via impedance of gravity-driven drainage (Fig. 2d) but also via capillary rise from the groundwater aquifer (Fig. 2e). This moistening enabled the vegetation to sustain higher E for at least a year (Fig. S9). As the droughts progressed into multi-year events, the impact of groundwater diminished due to a depletion of soil moisture stores regardless of whether groundwater dynamics were simulated.

When a heatwave occurs during a drought, and in particular early in a drought, the extra transpiration enabled by representing groundwater dynamics helps reduce the heat stress on vegetation (e.g. the reduction of 0.64 ^∘C of ΔT over the forests in 2002, Fig. 3a). This effect is particularly pronounced in regions with a shallower WTD (e.g. where the groundwater was within the first 5 m, there was a 0.5 ^∘C mean reduction in ΔT in the recent drought; Fig. 4d). Importantly, the role played by groundwater diminishes as the drought lengthens beyond 2 years (Fig. 3). Additionally, either the lack of deep roots or stomatal closure caused by high $D/T_{\mathrm{air}}$ can reduce the additional transpiration induced by groundwater. The latter plant physiology feedback dominates during heatwaves co-occurring with drought, even if the groundwater's influence has increased root zone water availability.

Our results highlight the impact of groundwater on both land surface states (e.g. soil moisture) and on surface fluxes and how this impact varies with the length and intensity of droughts and heatwaves. The results imply that the dominant mechanism by which groundwater buffered transpiration was through impeding gravity-driven drainage. We found a limited role for upward water movement from the aquifer due to simulated shallow WTD (which was broadly consistent with the observations in Fan et al., 2013). Further work will be necessary to understand how groundwater interacts with droughts and heatwaves and what these interactions mean for terrestrial ecosystems and the occurrence of the compound extreme events, particularly under the projection of intensifying droughts (Ukkola et al., 2020) and heatwaves (Cowan et al., 2014).

4.1 Changes in the role of groundwater in multi-year droughts

Groundwater is the slowest part of the terrestrial water cycle to change (Condon et al., 2020) and can have a memory of multi-year variations in rainfall (Martínez-de la Torre and Miguez-Macho, 2019; Martinez et al., 2016a). Our results show that the effect of groundwater on the partitioning of available energy between latent and sensible heat fluxes is influenced by the length of drought. As the drought extends in time, the extra E sustained by groundwater decreases (e.g. during the Millennium Drought, Fig. S9). The role of a drying landscape in modifying the partitioning of available energy between latent and sensible heat fluxes is well known and has been extensively studied (Fan, 2015; Miralles et al., 2019; Seneviratne et al., 2010). Our results add to the knowledge by quantifying the extent of the groundwater control, and eliciting the timescales of influence and the mechanisms at play. The importance of vegetation–groundwater interactions on multi-year timescales has been identified previously. Humphrey et al. (2018) hypothesised that climate models may underestimate the amplitude of global net ecosystem exchange because of a lack of deep-water access. Our regionally based results support this hypothesis and in particular highlight the importance of groundwater for explaining the amplitude of fluxes in wet periods as well as sustaining evapotranspiration during drought (Fig. 1).

4.2 Implications for land–atmosphere feedbacks during compound events

Our results show that during drought–heatwave compound events, the existence of groundwater eases the heat stress on the forest canopy and reduces the sensible heat flux to the atmosphere. This has the potential to reduce heat accumulating in the boundary layer and help ameliorate the intensity of a heatwave (Keune et al., 2016; Zipper et al., 2019). The presence of groundwater helps dampen a positive feedback loop whereby during drought–heatwave compound events, the high exchange of sensible and low exchange of latent heat can heat the atmosphere and increase the atmospheric demand for water (De Boeck et al., 2010; Massmann et al., 2019), intensifying drying (Miralles et al., 2014). The lack of groundwater in many LSMs suggests a lack of this moderating process and consequently a risk of overestimating the positive feedback on the boundary layer in coupled climate simulations. Our results show that neglecting groundwater leads to an average overestimate in canopy temperature by 0.2–1 ^∘C where the WTD is shallow (Fig. 4d) but as much as 5 ^∘C in single heatwave events (Fig. 6e–f), leading to an increase in the sensible heat flux (Fig. 4e).

The capacity of groundwater to moderate this positive land–atmosphere feedback is via modifying soil water availability. Firstly, soil water availability influenced by WTD affects how much water is available for E. In the shallow WTD regions, the higher soil water is likely to suppress the mutual enhancement of droughts and heatwaves (Keune et al., 2016; Zipper et al., 2019), particularly early in a drought. However, this suppression becomes weaker as the WTD deepens, in particular at depths beneath the root zone (e.g. 4.6 m in CABLE-GW) or as a drought lengthens. Our results imply that the land amplification of heatwaves is likely stronger in the inland regions (Hirsch et al., 2019), where the WTD is lower than 5 m and the influence of groundwater diminishes (Fig. S5), and once a drought has intensified significantly.

On a dry and hot heatwave afternoon, plant physiology feedbacks to high D and high T_air dominate transpiration and reduce the influence of groundwater in moderating heatwaves. In CABLE, stomatal closure occurs either directly due to high D (> 3 kPa) (De Kauwe et al. 2015) or indirectly due to biochemical feedbacks on photosynthesis at high T_air (> 32 ^∘C) (Kowalczyk et al., 2006); both processes reduce transpiration to near zero, eliminating the buffering effect of groundwater on canopy temperatures. While the timing of the onset of these physiology feedbacks varies across LSMs due to different parameterised sensitivities of stomatal conductance to atmospheric demand (Ball et al., 1987; Leuning et al., 1995) and different temperature dependence parameterisations (Badger and Collatz, 1977; Bernacchi et al., 2001; Crous et al., 2013), importantly, stomatal closure during heat extremes would be model invariant.

4.3 Uncertainties and future directions

Our study uses a single LSM, and consequently the parameterisations included in CABLE-GW influence the quantification of the role of groundwater on droughts and heatwaves. We note CABLE-GW has been extensively evaluated for water cycle processes (Decker, 2015; Decker et al., 2017; Mu et al., 2021; Ukkola et al., 2016b), but evaluation for groundwater interactions remains limited due to the lack of suitable observations (e.g. regional WTD monitoring or detailed knowledge of the distribution of root depths). Figure 1 gives us confidence that CABLE-GW is performing well, based on the evaluation against the GRACE, DOLCE, and GLEAM products, as well as previous work that showed the capacity of CABLE-GW to simulate E well (Decker, 2015; Decker et al., 2017). However, we also note that key model parameterisations that may influence the role of groundwater are particularly uncertain.

We need to be cautious about the “small” groundwater impact on the canopy temperature and associated turbulent energy fluxes during high D or high T_air (Figs. 3, 4, 6 and 7). The thresholds of D and T_air currently assumed by LSMs are in fact likely to be species specific. Australian trees in particular have evolved a series of physiological adaptations to reduce the negative impact of heat extremes. It is important to note that most LSMs parameterise their stomatal response to D for moderate ranges (< 2 kPa), which leads to significant biases at high D (Yang et al., 2019), a feature common in Australia and during heatwaves in general. New theory is needed to ensure that models adequately capture the full range of stomatal response to variability in D (low and high ranges). Similarly, while there is strong evidence to suggest that the optimum temperature for photosynthesis does not vary predictably with the climate of species origin (Kumarathunge et al., 2019) (implying model parameterisations do not need to vary with species), findings from studies do vary (Cunningham and Reed, 2002; Reich et al., 2015). Moreover, evidence that plants acclimate their photosynthetic temperature response is strong (Kattge and Knorr, 2007; Kumarathunge et al., 2019; Mercado et al., 2018; Smith et al., 2016; Smith and Dukes, 2013). As a result, it is likely that LSMs currently underestimate groundwater influence during heatwaves due to the interaction with plant physiology feedbacks. This is a key area requiring further investigation. For example, Drake et al. (2018) demonstrated that during a 4 d heatwave > 43 ^∘C, Australian Eucalyptus parramattensis trees did not reduce transpiration to zero as models would commonly predict, allowing the trees to persist unharmed in a whole-tree chamber experiment. Although De Kauwe et al. (2019) did not find strong support for this phenomenon across eddy covariance sites, if this physiological response is common across Australian woodlands, it would change our view on the importance of soil water availability (therefore groundwater) for the evolution of heatwave or even compound events. Coupled model sensitivity experiments may be important to determine the magnitude that such a physiological feedback would present and could guide the direction of future field/manipulation experiments.

Root distribution and root function and thereby how roots utilise groundwater are uncertain in models (Arora and Boer, 2003; Drewniak, 2019; Wang et al., 2018; Warren et al., 2015) and indeed in observations (Fan et al., 2017; Jackson et al., 1996; Schenk and Jackson, 2002). Models often ignore how roots forage for water and respond to moisture heterogeneity, limiting the model's ability to accurately reflect the plant usage of groundwater (Warren et al., 2015). In LSMs, roots are typically parameterised using a fixed distribution and normally ignore water uptake from deep roots. This assumption neglects any climatological impact of root distribution and the differentiation in root morphology and function (fine roots vs. tap roots), leading to a potential underestimation of groundwater utilisation in LSMs (see our deep root experiment, Fig. 6g–h). This assumption may be particularly problematic in Australia where vegetation has developed significant adaptation strategies to cope with both extreme heat and drought, including deeply rooted vegetation that can access groundwater (Bartle et al., 1980; Dawson and Pate, 1996; Eamus et al., 2015; Eberbach and Burrows, 2006; Fan et al., 2017). We also note that CABLE does not directly consider hydraulic redistribution, defined as the passive water movement via plant roots from moister to drier soil layers (Burgess et al., 1998; Richards and Caldwell, 1987). Neglecting hydraulic redistribution has the potential to underestimate the groundwater transported upwards and understate the importance of groundwater on ecosystems.

On the atmosphere side, the existence of groundwater increases the water flux from the land to the atmosphere, particularly in regions of shallow WTD, during the first 1–2 years of a drought. This has the potential to moisten the lower atmosphere and may encourage precipitation (Anyah et al., 2008; Jiang et al., 2009; Martinez et al., 2016b; Maxwell et al., 2011). However, our experiments are uncoupled from the atmosphere so while there is the potential for the higher E to affect the boundary layer moisture (Bonetti et al., 2015; Gilbert et al., 2017; Maxwell et al., 2007), clouds, and precipitation, we cannot conclude that it would until we undertake future coupled simulations.

Finally, we note we have focused on the role of groundwater in a natural environment. Humans extract large quantities of groundwater in many regions (Döll et al., 2014; Wada, 2016). Adding human management of groundwater into LSMs enables an examination of how this affects the vulnerability of ecosystems to heatwaves and drought and may ultimately identify those vulnerable ecosystems close to tipping points that are priorities for protection.

5 Conclusions

In conclusion, we used the CABLE LSM, constrained by satellite observations, to explore the timescales and extent to which groundwater influences vegetation function and turbulent energy fluxes during multi-year droughts. We showed that groundwater moistened the soil during the first ∼ 2 years of a multi-year drought, which enabled the vegetation to sustain higher evaporation (50–200 mm yr⁻¹ over the coastal forest regions) during drought onset. This cooled the forest canopy on average by 0.03–0.76 ^∘C in heatwaves during 2001–2019 and by as much as 5 ^∘C in regions of shallow water table depths in the heatwave in January 2019, helping to moderate the heat stress on vegetation during heatwaves. However, the ability of groundwater to buffer vegetation function varied with the length and intensity of droughts and heatwaves, with its influence decreasing with prolonged drought conditions. Importantly, we also demonstrated that the capacity of the groundwater to buffer evaporative fluxes during heatwaves is constrained by plant physiology feedbacks which regulate stomatal control, irrespective of soil water status. Given the increased risk of regional heatwaves and droughts in the future, the role of groundwater on land–atmosphere feedbacks and on terrestrial ecosystems needs to be better understood in order to constrain future projections.