Soil moisture–temperature coupling during extreme warm conditions in 2018 in Sweden: a case study with WRF-CTSM

Soil moisture–temperature coupling (SM–T) significantly influences the frequency and amplitude of heat extremes. It describes how variations in soil moisture affect surface air temperature conditions and vice versa. This study aims to determine the spatial extent and duration of SM–T in southern and central Sweden, an area increasingly recognized as a coupling hot spot, during the extreme warm conditions between May and August 2018 (MJJA 2018). The assessment of coupling is based on a multi-correlation overlay analysis of key coupling variables: surface soil moisture, evaporative fraction, and daily maximum 2 m temperature from four different simulations of the coupled regional climate model WRF-CTSM, along with a merged gridded GLEAM-E-OBS observational–reanalysis dataset. These datasets demonstrate robust precision in representing the magnitude and variability of the key coupling variables during the MJJA 2018 compared to in situ observations, though the precise timing and duration of the coupling are challenging to reproduce at the local scale. WRF-CTSM provides a more realistic depiction of the key coupling variables and their interactions when recent CTSM advancements are incorporated. On average, across the study region and all five datasets, SM–T persisted for 22 d throughout the MJJA period. The atmospheric leg alone (involving daily evaporative fraction and maximum 2 m temperature), averaged across datasets, contributed 92 % to the regional coupling duration.

1 Introduction

Given the high confidence in the increasing probability of compound hot and dry events accompanying the global mean temperature rise, as assessed by the IPCC AR6 Synthesis Report (Calvin et al., 2023), studies related to heat extremes are quickly expanding. From a physical science perspective, further research focused on understanding the underlying processes of heatwaves, including land-atmosphere (L–A) interactions, is still required (Barriopedro et al., 2023; Miralles et al., 2019), particularly on an individual event basis (Maraun et al., 2025). Here, an exceptionally warm period from May to August (MJJA) in 2018 over southern and central Sweden is chosen as a case study (Dirmeyer et al., 2021). This is due to its extensive anomalous warm and dry conditions spanning most of northern and central Europe. This prolonged heat period generated multifaceted impacts, making it a continuing topic of investigation within the research community (Bakke et al., 2020; Bastos et al., 2020; Dirmeyer et al., 2021; Hoy et al., 2020; Petch et al., 2020; Peters et al., 2020; Rousi et al., 2023; Salomón et al., 2022; Tak et al., 2024; Wilcke et al., 2020). The primary factor determining the onset of the warm period resulting in the heatwave conditions in 2018 was a sustained atmospheric blocking high-pressure system over the British Isles and the North Sea (Kueh and Lin, 2020). The pronounced effects were encountered over northern Europe, where drought persisted from late spring through autumn, spanning the entire growing and harvesting seasons and generating agricultural losses (C3S, 2024). Following the meteorological drought (precipitation deficit), the intensity of the event revealed rapid progression to agricultural drought (soil water deficit), and ultimately, to hydrological drought (reduced streamflow and groundwater storage) (Bakke et al., 2020). These series of circumstances are more common in mid-latitude regions and are all unusual in the generally moisture-abundant high latitudes (Dirmeyer et al., 2021; Koster et al., 2009). Nevertheless, the area of southern and central Sweden is increasingly seen as a L–A coupling hot spot region (Fu et al., 2024; Maraun et al., 2025; Yuan et al., 2023).

According to the Budyko framework (Fig. 1), high latitudes are considered part of the wet climate/soil moisture regime where the evaporation (or evapotranspiration, Miralles et al., 2020) is generally not driven by the wetness state of soil but by the available energy (Dirmeyer et al., 2009; Seneviratne et al., 2010). The mid-to-high latitudes can, however, experience a shift from the energy-limited to the soil moisture-limited evaporation regime at shorter periods when the enhanced warming decreases the soil moisture (SM) below the non-stationary critical SM threshold (Dirmeyer et al., 2021; Hsu and Dirmeyer, 2023b). This shift predominantly occurs during hot and dry summers, which are expected to become more common under a changing climate (Denissen et al., 2022; Hsu and Dirmeyer, 2023a, b). A key process in these L–A interactions is the soil moisture–temperature coupling (SM–T) illustrated in the uppermost part of Fig. 1. The respective illustrations refer to the effects of the 2 m temperature increase across different SM regimes. SM–T may occur when a region enters the transitional SM regime, while it is not associated with the dry and wet SM regimes (Maraun et al., 2025). In the dry regime, the soil water content is low and evaporation is minimal, so further drying does not induce an additional 2 m temperature rise. In the wet regime, the SM content is always sufficient to meet the atmospheric water demand. In the transitional regime, the initial temperature increase may result in a SM deficit that may generate a decrease in evaporation and thereby evaporative fraction, which may feed back to further 2 m temperature increase (Dirmeyer, 2011). This pathway of L–A interactions in the transitional SM regime can significantly magnify the intensity of the initial extreme temperature event (Maraun et al., 2025; Seneviratne et al., 2010). A shift from the wet to the transitional regime occurs when SM depletes and the typical positive relationship (Fig. 1, upper right) reverses, such that rising temperature leads to declining evaporation and subsequently, evaporative fraction.

Figure 1 Conceptual framework of the dependence of evaporative fraction on soil moisture in different evaporation and climate/soil moisture regimes [modified from Seneviratne et al., 2010]. The uppermost section of the figure illustrates the effects of the 2 m temperature increase across different regimes that may influence the occurrence of soil moisture–temperature coupling. Arrows indicate the direction of relationships between variables. For example, a rise in 2 m temperature leads to a decrease in soil moisture in the transitional climate/soil moisture regime.

Soil moisture–temperature interactions modify the surface energy balance, which in turn accelerates global warming in Earth System Models (Qiao et al., 2023). This underlines that the ability to represent the coupling will affect the future projections of heat extremes (Hsu and Dirmeyer, 2023b; Schwitalla et al., 2025; Zarakas et al., 2024). Our study employs the recently developed coupled regional climate model WRF-CTSM (UCAR, 2021; Mužić et al., 2025) for the extreme heat analysis. It combines the Weather Research and Forecasting model, WRF (Skamarock et al., 2021), and the Community Terrestrial Systems Model, CTSM, which at present corresponds to the Community Land Model, CLM (version CLM5) (Lawrence et al., 2019). WRF and CTSM are currently considered among the most sophisticated atmospheric and land surface models. This facilitates tackling some of the current knowledge gaps in heatwave research, indicated by Barriopedro et al. (2023) and Miralles et al. (2019). Accordingly, the two main objectives of this study are:

1. To evaluate the representation of the key variables for SM–T analysis during the warm MJJA 2018 period in a set of WRF-CTSM simulations and an additional independent gridded dataset derived from observations and reanalysis inputs (GLEAM-E-OBS).

2. To examine the spatial extent and duration of the SM–T and its components throughout the MJJA 2018 in a set of WRF-CTSM simulations and GLEAM-E-OBS.

These objectives follow the findings from the study by Dirmeyer et al. (2021), who, based on in situ and reanalysis datasets, found that SM–T amplified the strength of the drought and heatwave as part of the summer in 2018 in the same study region (i.e., southern and central Sweden). Here, four different WRF-CTSM simulations are employed to assess if the model captures the extremity of the event as well as the variability and magnitude of the key coupling variables throughout MJJA 2018; surface soil moisture (SM_SURF), evaporative fraction (EF), maximum 2 m temperature (T_MAX) as well as their drivers; net radiation (R_NET), sensible heat (SH), latent heat (LH), and precipitation (PR). To assess the local and regional spatial extent and duration of daily SM–T and its components, model simulations are accompanied by additional datasets described in the following.

2 Data and Methods

2.1 Datasets

Datasets at daily time resolution used to support the SM–T analysis during the MJJA 2018 include in situ observations, a gridded satellite-based dataset (ESA CCI), four WRF-CTSM simulations, and a gridded merged observational–reanalysis dataset (GLEAM-E-OBS). For clarity, observational products and model outputs are both referred to as “datasets” throughout this study. An overview of the respective datasets is listed in Table 1 and explained in the subsequent sections.

Table 1 An overview of the employed datasets.

2.1.1 Stations

ICOS ground observations from Drought 2018 Team and ICOS Ecosystem Thematic Centre (2020) (station Lanna), and from Warm Winter 2020 Team and ICOS Ecosystem Thematic Centre (2022) (stations Norunda, Hyltemossa, Rosinedal-3, and Degerö) are flux tower data products in FLUXNET2015 format (Pastorello et al., 2020). These stations provide daily data for the investigated variables: SM_SURF (corresponding to 0–0.1 m soil depth), SH and LH (for the calculation of EF), T_MAX, PR, and R_NET. Their locations within the study area are presented in Fig. 2.

Figure 2 An overview of the ground stations' locations (symbols) and the study area's elevation (shading). Lakes within the study area are denoted in white.

The study area and the selection of stations follow the methodology adopted from Dirmeyer et al. (2021). Here, three additional stations were considered (Norunda, Rosinedal-3, and Degerö), while Svartberget was omitted since the newest version of its dataset does not provide data for SM_SURF. Only the best-quality station data that correspond to the directly measured values are used in the analysis. In situ observations are considered as ground truth and used as an evaluation dataset in our study.

2.1.2 ESA CCI

As indicated above, the availability of stations providing SM data in the study region is generally sparse. Additionally, the WRF-CTSM SM_SURF model output is, among the variables used in this study, the only one not previously validated over Sweden (Mužić et al., 2025). Hence, the European Space Agency's Climate Change Initiative (ESA CCI) Soil Moisture (SM) product (Dorigo et al., 2017; Gruber et al., 2019; Preimesberger et al., 2021) is utilized as an additional evaluation dataset for the assessment of the regional-scale WRF-CTSM performance in simulating SM_SURF. It is a daily global satellite SM_SURF dataset available at 0.25° spatial resolution. The ESA CCI COMBINED product used in this study is derived from a range of active and passive satellite remote-sensing instruments measuring SM_SURF (the corresponding soil depth layer is 0–0.05 m) in the microwave range (Preimesberger et al., 2021). Here, only good-quality (no inconsistency flagged) data are selected, so all unreliable data points are masked out.

2.1.3 Model simulations

A set of four WRF-CTSM simulations for the year 2018 on 10.5 km horizontal resolution is performed over a domain covering Norway, Sweden, and Finland, shown in Fig. 2; the rectangle indicates the analysis region of southern and central Sweden. Employing four WRF-CTSM simulations generates more robust findings and provides insight into the role of model advancements over time for the SM–T assessment. All four WRF-CTSM simulations are based on the initial and boundary conditions from ERA5 (Hersbach et al., 2020), while the spectral nudging of temperature, horizontal winds, and humidity is performed above the planetary boundary layer (PBL). WRF parameterization schemes used are the Rapid Radiative Transfer Model for General Circulation Models (RRTMG) radiation scheme (Iacono et al., 2008), Kain-Fritsch (KF) convection scheme (Kain, 2004), WRF Single-Moment 6-class Microphysics Scheme (WSM6) (Hong and Lim, 2006), and the Yonsei University Scheme (YSU) PBL scheme (Hong et al., 2006). The leaf area index (LAI) is prescribed by the “satellite phenology” option in CTSM. Alongside the respective common settings, the main differences among the four simulations are summarized in Table 1 and described below.

Run 1 is considered a reference simulation since it follows the setup used in the WRF-CTSM evaluation study by Mužić et al. (2025). It includes a long simulation period between 2010 and 2022. The subsequent simulations (Runs 2–4) simulate only the year 2018 using a 2-year spin-up starting from 1 January 2016, of Run 1. Runs 2–4 use configurations with recent CTSM advancements (Lawrence et al., 2019), namely, a newer CTSM model version (including physics and parameterizations), the stomatal conductance method by Medlyn et al. (2011), 20 (non-bedrock hydrologically and biogeochemically active) soil layers spanning up to 8.6 m depth, and subgrid heterogeneity of plant functional types (PFTs). Run 3 additionally includes the plant hydraulic stress configuration (PHS) suitable for heatwave research, developed by Kennedy et al. (2019). Finally, Run 4 also incorporates the biomass heat storage configuration (BHS) introduced by Meier et al. (2019). Hence, Run 4 uses WRF coupled with CTSM 5.1 in its current default configuration. This way, the simulations' numbering order (Runs 1–4) follows the CTSM advancements tailored for improved representation of the land surface water and energy exchange.

2.1.4 GLEAM-E-OBS

Global Land Evaporation Amsterdam Model (GLEAM) (Hulsman et al., 2023; Koppa et al., 2022; Miralles et al., 2011, 2025; Zhong et al., 2022) and E-OBS (Cornes et al., 2018) datasets are in our study merged to supplement the WRF-CTSM MJJA 2018 analysis results with an independent gridded dataset originally provided at the analogous spatial resolution as the model output (0.1°). As it is a merged observational–reanalysis dataset, it is not considered as ground truth and is evaluated the same way as the WRF-CTSM simulations by benchmarking against in situ observations.

GLEAM 4.1a is a recently issued global dataset of surface water exchange variables on a 0.1° grid and with a daily temporal resolution (Miralles et al., 2025). It is based on satellite and reanalysis data and is included in this analysis as this version for the first time also provides SH, which is needed to calculate the EF. GLEAM variables employed in this study as part of the GLEAM-E-OBS dataset are the SM_SURF (soil depth layer corresponding to 0–0.1 m depth), actual evaporation (includes transpiration, bare soil evaporation, open water evaporation, interception loss, and snow sublimation), and SH. GLEAM 4.1a utilizes SM_SURF input data from the ESA CCI COMBINED 8.1 dataset.

E-OBS is a gridded observational dataset on 0.1° resolution over Europe produced based on the network of the European Climate Assessment and Data (ECA&D) stations (Cornes et al., 2018). E-OBS variables used in the blended GLEAM-E-OBS dataset include daily T_MAX and daily total PR.

2.1.5 Datasets integration

Gridded datasets from Table 1 are compared to the station data by selecting the nearest grid cells to each station and masking their data to reflect the temporal coverage of the respective stations. ESA CCI and GLEAM-E-OBS datasets are regridded to the spatial resolution and coordinate system of WRF-CTSM using bilinear interpolation. Spatio-temporal gaps in GLEAM-E-OBS are primarily associated with the discontinuities in the GLEAM SM_SURF variable as a result of different satellite passing cycles, complex topography, surface water, dense vegetation, urban structures, and frozen and snow-covered soil conditions (Dorigo et al., 2017). These gaps are masked from the rest of the GLEAM-E-OBS variables to ensure a consistent dataset. Furthermore, lake grid cells are masked from all datasets. The equations used for the unit conversions of SM_SURF and LH among datasets are given in Sect. S1. The evaluation metrics for comparing WRF-CTSM simulations and GLEAM-E-OBS with in situ observations and ESA CCI are the Pearson correlation coefficient (PCC), root mean square error (RMSE), and absolute bias (AB) described in Sect. S2.

2.2 Methods

The assessment of SM–T is in our study based on the two-legged coupling metric (Dirmeyer, 2011) with adjustments to facilitate the quantification of the coupling duration and the role of the coupling components throughout the exceptionally warm MJJA 2018 period. The approach consists of a correlation coefficient analysis (Dirmeyer et al., 2009) among daily SM_SURF, EF, and T_MAX. EF represents the ratio of LH to the available energy and is calculated as:

$$\begin{array}{}\text{(1)}& \mathrm{EF}={\displaystyle \frac{\mathrm{LH}}{\mathrm{SH}+\mathrm{LH}}}\end{array}$$

where EF is the evaporative fraction, and LH and SH are latent and sensible heat (W m⁻²), respectively (Shuttleworth et al., 1989).

The soil moisture–temperature coupling diagnostic is in our study defined as:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \mathrm{Coupling}\left({\mathrm{SM}}_{\mathrm{SURF}},\mathrm{EF},T_{\mathrm{MAX}}\right)\\ & =\left\{\begin{array}{ll}\mathrm{True},\mathrm{if}& r\left(T_{\mathrm{MAX}},{\mathrm{SM}}_{\mathrm{SURF}}\right)<\mathrm{0}\wedge \\ & r\left({\mathrm{SM}}_{\mathrm{SURF}},\mathrm{EF}\right)>\mathrm{0}\wedge \\ & r\left(\mathrm{EF},T_{\mathrm{MAX}}\right)<\mathrm{0}\wedge \\ & r\left(T_{\mathrm{MAX}},t\right)>\mathrm{0}\\ \mathrm{False},\mathrm{otherwise}& \end{array}\right.\end{array}\end{array}$$

where SM_SURF is surface soil moisture, EF is evaporative fraction, T_MAX is daily maximum temperature, t is time, and r is the Pearson correlation coefficient.

This definition (Eq. 2) employs variables used in the SM–T coupling assessment by Schwingshackl et al. (2017). Here, SM–T is broken down and assessed according to three coupling components: T_MAX:SM_SURF, SM_SURF:EF (terrestrial leg), and EF:T_MAX (atmospheric leg). Accordingly, the analysis includes three daily time series (one for each of the three coupling components) of the PCCs obtained from periods corresponding to a 92 d running window centered on the day in question. The respective running windows for the first and last 92 d of the MJJA period include days before May and after August. SM–T is considered present only in areas that simultaneously exhibit a positive T_MAX trend over the 92 d running window and a combination of correlations within the same window: negative T_MAX:SM_SURF, positive SM_SURF:EF, and negative EF:T_MAX. This combination is referred to as a complete coupling (Eq. 2). Selecting 92 d provides a robust statistic while not smoothing out the variability throughout the MJJA period (Findell et al., 2015). Building on the framework of Rousi et al. (2023), who examined SM:LH using 92 d running periods, this multi-correlation overlay approach serves as an alternative to the piecewise/segmented regression method used by, e.g., Dirmeyer et al. (2021), Schwingshackl et al. (2017), Ukkola et al. (2018), Tak et al. (2024), and Fu et al. (2024). Both methodologies indicate the direction and magnitude of the coupling, but the advantage of the multi-correlation overlay approach is that it does not rely on critical SM values. The critical SM is not stationary as the wilting point; it varies in magnitude throughout time at a certain location due to climate variability (Hsu and Dirmeyer, 2022)

SM_SURF is used as part of SM–T, as the surface soil layer shows the strongest coupling with the PBL (Dirmeyer et al., 2021) and is thus the most relevant in correlation analysis. A sensitivity analysis within our study (not shown) illustrates that the sign of the correlations, including SM_SURF (T_MAX:SM_SURF and SM_SURF:EF) remains qualitatively the same if replaced by the weighted mean SM across all soil layers, as also reported in Fu et al. (2024) and Schwingshackl et al. (2017). Although evaporation (or LH) can be adopted in the coupling analysis, EF is advantageous as its variability is not directly tied to net radiation and wind speed (Qiu et al., 2020), making it a more suitable measure of land surface dryness (Ukkola et al., 2018). Using evaporation alone can obscure the impact of SM on the L–A interactions (Qiu et al., 2020). Given the recent inclusion of SH in GLEAM 4.1a, which enables calculation of the EF, the GLEAM-E-OBS dataset is prepared to accompany the WRF-CTSM simulations in the coupling analysis.

3 Results & Discussion

3.1 WRF-CTSM and GLEAM-E-OBS evaluation

A station-based evaluation of WRF-CTSM simulations and GLEAM-EOBS for daily SM_SURF, EF, T_MAX, SH, LH, and total PR is presented in Fig. 3. The mean across the five stations and gridded datasets (Fig. 3) is computed after masking the daily values of each gridded dataset according to the missing observations of the corresponding station. Additionally, the mean of R_NET and its radiative fluxes is available for three out of five stations (Norunda, Hyltemossa, and Degerö) presented in Fig. S1 in the Supplement. These local-scale comparisons are addressed in the following two sections.

Figure 3 Grid cell level evaluation of WRF-CTSM simulations and GLEAM-E-OBS dataset against in situ observations averaged over five stations, illustrated as a 10 d centered running mean time series. Thin lines represent MJJA 2018, while thick lines (stations, ESA CCI, GLEAM-E-OBS, and Run 1) denote the mean of the preceding three years for which data are available at all considered stations (2015–2017). Key soil moisture–temperature coupling variables are evaluated in panels (a)–(c). Surface soil moisture depths are 0.05 m (ESA CCI) and 0.1 m (stations, GLEAM-E-OBS, and Runs 1–4). Evaporative fraction values are constrained between 0 and 1. Evaluation metrics in the bottom right are for daily values during MJJA 2018 and calculated relative to the station data. They include the Pearson correlation coefficient (PCC), root mean square error (RMSE), and absolute bias (AB) described in Sect. S2.

3.1.1 Surface soil moisture evaluation

This section includes a validation of SM_SURF from ESA CCI against station data to ensure its suitability in reproducing in situ observations and serving as a regional-level SM_SURF evaluation dataset despite initially being provided on a coarser resolution (0.25°) than the WRF-CTSM and GLEAM-E-OBS datasets (0.1°). Additionally, ESA-CCI represents a shallower soil depth layer (0.05 m) than the other datasets (0.1 m). In Fig. 3a, the evaluation metrics consistently indicate that the GLEAM-E-OBS SM_SURF is the most aligned with the in situ observations; the slightly higher absolute bias compared to ESA CCI is negligible. ESA CCI and Run 4 are also closely aligned with the station data. Run 1 SM_SURF is the driest among datasets at the selected locations (underestimation against the stations is approximately 0.1 mm³ mm⁻³) throughout the MJJA 2018 period. The Run 1 underestimation is also depicted over most of the study region in comparison to ESA CCI (Fig. 4b), whereas this is not a prevailing pattern in GLEAM-E-OBS and Run 4 (Fig. 4a–c). The corresponding mean regional MJJA 2018 absolute SM_SURF values in GLEAM-E-OBS, Run 1, and Run 4 are presented in Fig. S2. Figures 4 and S2 show that the bias pattern among datasets aligns with the spatial distribution of SM_SURF. Regions characterized by higher mean SM_SURF values (northwest, mountainous areas) exhibit positive biases, while regions with lower SM_SURF values (southeast) show negative biases relative to ESA CCI. These summer SM_SURF WRF-CTSM biases could be attributed to several sources, including biases in simulated precipitation, parameterizations in the land and atmosphere models, differences in soil depth, soil texture, freeze-thaw representation, vegetation type and properties, as well as the discrepancy between point scale observations and the model grid cell size (Li et al., 2020).

Figure 4 Mean regional MJJA 2018 surface soil moisture absolute bias for (a) GLEAM-E-OBS, (b) Run 1, and (c) Run 4 against ESA CCI. Surface soil moisture depths are 0.05 m for the ESA CCI product and 0.1 m for GLEAM-E-OBS, Run 1, and Run 4. Regional scale evaluation metrics in the upper left section of each panel include Pearson correlation coefficient (PCC), root mean square error (RMSE), and absolute bias (AB), described in Sect. S2.

3.1.2 Evaluation of near-surface meteorological and energy flux variables

WRF-CTSM simulations and GLEAM-E-OBS illustrate an overestimation of LH throughout most of MJJA 2018 compared to station data (Fig. 3e). On the contrary, SH is generally underestimated (Fig. 3d). The magnitude of the respective biases is the highest in July 2018, which, according to Dirmeyer et al. (2021), represents the month of the most pronounced L–A interactions over the study region. Run 1 and GLEAM-E-OBS commonly indicate distinct temporal patterns as opposed to Runs 2–4, which are more closely aligned, likely due to their similar setups. When all metrics from the LH and SH station-based MJJA 2018 evaluations (Figs. 3d–e and S4–S8, summarized in Table S2) are considered, Run 4 yields the best overall performance. Accounting for the MJJA period and all evaluation metrics, GLEAM-E-OBS has the lowest performance as it shows a more pronounced overestimation of LH and underestimation of SH throughout MJJA 2018 (Fig. 3d–e). This is combined with the highest MJJA T_MAX (Fig. 3c) among the datasets, underpinning that the key coupling variables in GLEAM-E-OBS are independent and derived from different sources. However, considering the SH and LH throughout July specifically (Fig. 3d–e), GLEAM-E-OBS shows improved performance as it reproduces the decreasing EF tendency observed in data from most stations (Figs. 3b and S4–S8). GLEAM–E-OBS does not consistently outperform the WRF-CTSM simulations in PCC for LH and SH at the station level during MJJA or July 2018. Nevertheless, it does show a higher PCC in MJJA 2018 when all key coupling variables are considered together (Tables S2 and S3). Correlation analysis is particularly relevant in our SM–T assessment, although in that context it primarily refers to the correlations among the key coupling variables within each dataset.

Accurately simulating EF is the most challenging among the key coupling variables. EF is sensitive to local surface SH and LH variations (including, e.g., intermittent cloudiness) and thus fluctuates throughout the warm period. Wang et al. (2019), in their station-based analysis during and outside heatwave periods, found that the EF reduction throughout heatwaves largely varies among study sites depending on local climate conditions, soil properties, and land cover types, but also within the same land cover type, indicating that responses are species-dependent. Hence, in contrast to the grid cells of the evaluated datasets, station observations are more sensitive to local fluctuations and show a higher magnitude in SH and a larger decrease in EF throughout July (Fig. 3b).

The high accuracy of the WRF-CTSM simulations and GLEAM-E-OBS in representing in situ daily T_MAX is further illustrated in Fig. 3c. The mean daily T_MAX MJJA 2018 PCC is 0.99 in all evaluated datasets. Considering also RMSE and AB, WRF-CTSM simulations are marginally more precise in representing the daily T_MAX than GLEAM-E-OBS (Fig. 3c). As expected, high surface warming is accompanied by high R_NET (Fig. S1). However, the R_NET is in Runs 2–4 represented through compensating fluxes, as the underlying processes differ from the observations. Runs 2–4 overestimate incoming solar radiation and underestimate longwave fluxes, implying reduced cloud cover (Fig. S1). Lower SH from a colder surface coincides with higher T_MAX, suggesting underestimated turbulent transport and associated cloud cover biases in the WRF-CTSM simulations. WRF has several options for PBL schemes, and a different choice could have potentially improved the results (e.g., Wang et al., 2021a), but it is beyond the scope of our study to assess a wide range of different schemes in WRF. Issues with cloud representation may have also contributed to a slight overestimation of PR (Fig. 3f) in Runs 2–4, likely driving the LH increase during July. Although not consistent across all sites, GLEAM-E-OBS PR generally aligns more closely with observations than Runs 2–4 (Table S4).

3.2 Anomalies in MJJA 2018

GLEAM-E-OBS and Run 1 generally agree on the areas of the most pronounced MJJA 2018 anomalies in SM_SURF, LH, SH, and T_MAX compared to the mean MJJA of other years throughout 2010–2022 (Fig. 5). The negative SM_SURF and positive T_MAX anomalies are spread over the entire study region in both datasets, as also shown by previous studies using ERA5 reanalysis on 0.25° resolution (Dirmeyer et al., 2021; Tak et al., 2024). Our study further emphasizes that a reduced EF (as depicted by anomalously low LH and high SH) persisted the longest in the southeastern part of Sweden (Fig. S9). This area thus reveals the most favourable conditions for the shift towards a transitional SM regime. This aligns with Sjökvist et al. (2019), who demonstrate that the confined areas of the southeastern part of Sweden had the longest number of consecutive heatwave days (up to 30 d), defined as days where the daily mean temperature exceeds 20 °C.

Figure 5 Mean MJJA 2018 anomalies in SM_SURF, LH, SH, and T_MAX in GLEAM-E-OBS (left column) and Run 1 (right column) based on the MJJA 2010–2022 reference period, excluding the year 2018.

3.3 Soil moisture–temperature coupling

3.3.1 Robust regional-scale patterns

A regional-scale SM–T during MJJA 2018 in GLEAM-E-OBS and the four WRF-CTSM simulations is presented in Fig. 6. The first column shows the spatio-temporal coverage of the complete coupling across all datasets, while the subsequent columns display the areas and durations of the individual coupling components. On average across the study area and all datasets, the complete SM–T persisted for 22 d. There is, however, a disparity in the total coverage and coupling duration among datasets (these differences are discussed in Sect. 3.3.2 and 3.3.3).

Figure 6 Fraction of days during MJJA 2018 that experienced soil moisture–temperature coupling. A time series of Pearson correlation coefficients (PCC) at each grid cell is based on 92 d running periods. Rows (a)–(e) denote regional datasets while columns demonstrate the parts of the coupling analysis, namely complete coupling, negative T_MAX:SM_SURF, positive SM_SURF:EF, and negative EF:T_MAX correlations, respectively, each co-occurring with increasing T_MAX. The mean duration over the study region is additionally indicated in each plot.

Overall, it is common for all datasets that the spatial pattern of the complete SM–T duration generally follows the pattern of the duration of the atmospheric leg (EF:T_MAX) (first and last columns in Fig. 6). This reflects a strong control of the atmosphere, particularly surface net radiation, temperature, and vapor pressure deficit (VPD) on EF (Wang et al., 2019) in high-latitude regions dominated by an energy-limited evaporation regime (Schwingshackl et al., 2017; Seneviratne et al., 2010). This is supported by the prevailing positive anomalies in LH throughout MJJA 2018 in station data, GLEAM-E-OBS, and Run 1 on in situ (Figs. 3e and S4e–S8e) and regional (Figs. 5c–d, S11, and S14c–d) scales. The respective regional-scale analyses further indicate that although the magnitude of the positive T_MAX anomaly is similar for May and July 2018 (locally exceeding 6 °C), the positive anomaly in evaporation is less pronounced and widespread in July, when most of the coupling occurred (Dirmeyer et al., 2021).

Given the strong control of changes in evaporation on coupling duration, the coupling analysis in our study disregards cases in which those changes are not part of the SM–T chain. By explicitly accounting for the condition of increasing T_MAX, we exclude cases when enhanced LH cooled the surface, leading to lower SH and T_MAX. Such cases typically occur during the early warm period when soils are still wet from snowmelt, or shortly after rainfall events during the warm season. This condition then rules out cases in which EF:T_MAX reveals a negative relation, but the T_MAX is dampening, or cases in which evaporation depletion is not relevant enough to trigger substantial changes in EF that would enhance the T_MAX. Furthermore, observational-based studies in Sweden and central Europe indicate instances of decreased evaporation during the warm 2018 period due to increased VPD and a reduction in stomatal conductance, even when SM is still not limited (Lindroth et al., 2020; Salomón et al., 2022). This plant stress response is a L–A interaction that further enhances warming, although it does not involve the effect of SM depletion (Fu et al., 2022; Grossiord et al., 2020; Schönbeck et al., 2022; Sulman et al., 2016). These interactions commonly occur on daily timescales and are excluded from our SM–T definition, as it requires a co-occurrence of all correlations, thereby ensuring the necessary transition to the soil moisture–limited evaporation regime. These constraints on the definition thus partly explain why each of the gridded datasets in Fig. 6 exhibits areas of non-existent SM–T throughout MJJA 2018. Additional factors that shape the spatio-temporal pattern of the coupling are discussed in the following sections.

3.3.2 Regional differences between GLEAM-E-OBS and WRF-CTSM runs

The complete SM–T in GLEAM-E-OBS occurred over most of the study region and lasted the longest, on average, 32 d throughout MJJA 2018 (Fig. 6a). Coupling in WRF-CTSM covered a smaller area and lasted for 17–23 d, depending on the simulation. GLEAM-E-OBS, Run 1 and Run 2 demonstrate that the coupling persisted the longest in the southeastern parts of Sweden (up to about half of the MJJA 2018 period, ∼ 60 d). The reduced sensitivity to coupling in WRF-CTSM runs in the respective area in comparison to GLEAM-E-OBS is related to the PR overestimation, as depicted in the regional-scale anomalies in Run 1 (Fig. S15). As LH inherently varies following changes in PR during the dry period, more frequent increases in LH than in the observations, specifically during July (Figs. S4–S8), resulted in a shorter coupling in WRF-CTSM.

An additional factor contributing to the differences in spatio-temporal coupling patterns stems from how datasets represent vegetation response to extended warm and dry periods. Both GLEAM-E-OBS and WRF-CTSM are equipped to represent evaporation over vegetated areas during water-limited conditions. In GLEAM-E-OBS, plant evaporation is estimated with a machine-learning approach that accounts for vegetation characteristics such as height, optical depth, and LAI, and incorporates plant access to groundwater (Hulsman et al., 2023; Miralles et al., 2025). WRF-CTSM additionally includes stomatal conductance and plant hydraulic parameterizations (Kennedy et al., 2019; Lawrence et al., 2019; Medlyn et al., 2011). As a result, WRF-CTSM includes more processes related to vegetation physiology when simulating evaporation. A lower complexity of factors leading to simultaneous correlations in GLEAM-E-OBS may yield a smoother signal and longer overlay of correlations than in the WRF-CTSM runs.

3.3.3 Regional differences among WRF-CTSM runs

The spatial extent and persistence of the complete SM–T decreases from Run 1 to Run 4 (Fig. 6b–e), along with the inclusion of more factors for physical representation of evaporation. Since Runs 2–4 show improved performance across LH and SH compared to Run 1 (Table S2) and include a single setup change per simulation (although with multiple parameters), their differences are examined in more detail in this section. Figure 7 supports these comparisons by illustrating the regional mean MJJA 2018 absolute difference in the complete coupling duration, SM_SURF, LH, SH, and T_MAX between Runs 3 and 2 (reflecting the effect of including the PHS configuration) and between Runs 4 and 3 (demonstrating the impact of including the BHS configuration). The analysis of the specific biogeophysical parameters (e.g., vegetation hydraulic conductivity parameters, soil physiographic parameters) and their relative importance in the coupled and uncoupled CTSM configurations is, however, an ongoing community effort that is beyond the scope of this study (Dagon et al., 2020; Gao et al., 2021; Jiménez-Rodríguez et al., 2024; Zarakas et al., 2024).

Figure 7 The impact of including plant hydraulic stress (PHS) and biomass heat storage (BHS) configurations illustrated as an MJJA 2018 absolute difference between Runs 3 and 2 (a–e) and between Runs 4 and 3 (f–j) in the number of days with complete coupling (a, f), daily surface soil moisture (b, g), latent heat (c, h), sensible heat (d, i), and maximum 2 m temperature (e, j).

While both the PHS configuration in Run 3 and the soil moisture stress (SMS) configuration used in Run 2 account for changes in stomatal conductance due to SM deficit and atmospheric VPD, Run 2 does not include the plant hydraulic failure, as it cannot resolve concurrent limitations to stomatal conductance associated with declining soil water or excessive plant xylem tension. In contrast, Run 3 with PHS simulates hydraulic limits on plant transpiration by resolving plant water potential at the leaf, stem, and root levels. A higher evaporation in Run 3 (Fig. 7c) primarily refers to transpiration (Fig. S16c) and is likely a result of a compensatory root water uptake from deeper soil layers during water stress conditions (Kennedy et al., 2019).

The impact of additionally including BHS (Run 4) in the PHS run (Run 3) leads to a reduction in T_MAX specifically over the areas with dense forest cover (Fig. S3), where the mean MJJA 2018 is reduced by up to 0.6 °C (Fig. 7j). This is expected since BHS reduces diurnal temperature variations in forests, owing to their high heat storage capacity of stems and leaves, to be more aligned with observations. In the BHS configuration, daytime T_MAX is lower while nighttime minimum 2 m temperature is higher, resulting in a slightly elevated mean daily temperature (Meier et al., 2019). Accordingly, no pronounced changes are found in LH and SH (Fig. 7h–i).

The evaporation rates in this study region are likely strongly influenced by the needleleaf evergreen forest, which is the dominant land cover type (Fig. S3). Figure 8 provides insight into the regional-scale development of evaporation components throughout the MJJA 2018. The prevalence of forests reveals that canopy transpiration accounts for the largest share among the evaporation components and determines the overall evaporation pattern. On average over the study region, soil evaporation and canopy transpiration demonstrate a decline in early July. Canopy evaporation (or interception) accounts for the smallest share among the evaporation components, and is, as opposed to the other two, not affected by SM availability (Lawrence et al., 2020) but fluctuates due to changes in atmospheric conditions throughout the warm period (Wang et al., 2021b). The fraction of evaporation components as part of the total evaporation is analogous in Run 4 (Fig. 8) and Runs 1–3 (Fig. S16); however, their magnitudes differ due to the parameterizations used in those simulations (e.g., the stomatal behavior parameterization affecting transpiration rates).

Figure 8 92 d centered running mean time series of evaporation components in Run 4 during MJJA 2018, averaged over the land grid cells in the study domain. Colored areas represent contributions of different evaporation components to the total evaporation, i.e., the sum of the colored areas at each time step represents total evaporation.

3.3.4 The role of coupling components

Since the complete SM–T closely follows the spatio-temporal coverage of the atmospheric leg (first and last columns in Fig. 6), Fig. 9 further demonstrates that the persistence of EF:T_MAX has a major role in determining whether periodic shifts to a soil moisture–limited evaporation regime occur in this study region (as also indicated in Fig. 1). Figure 9 shows how the stand-alone duration of each coupling component contributes to the complete coupling duration in MJJA 2018. When considering only the atmospheric leg, the probability of estimating the duration of the complete coupling in this study region ranges from 67 % to 98 %, depending on the dataset. This indicates that while atmospheric leg accounts for the smallest number of days within the complete coupling (Fig. 6), it predominantly occurs at times when the other two coupling components are also met (Fig. 9). The regional average across GLEAM-E-OBS and WRF-CTSM simulations is 92 %, suggesting a high, but not definitive, ability to determine the coupling duration based only on the number of days spent in the negative EF:T_MAX (with increasing T_MAX) during MJJA 2018. A slight discrepancy between the persistence of the atmospheric leg and the complete coupling period resembles a need to integrate all coupling components in the SM–T definition. Studies accounting for the atmospheric leg only commonly include additional metrics to capture SM–T (Gevaert et al., 2018; Knist et al., 2017; Maraun et al., 2025). A decrease in the number of days in T_MAX:SM_SURF and SM_SURF:EF as part of the complete coupling from Run 1 to Run 4 (Fig. 9b–e) likely indicates attenuation of the direct impact of SM_SURF on the complete coupling duration as the evaporation becomes affected by additional factors (discussed in Sect. 3.3.3).

Figure 9 The proportion of days spent in each of the three coupling components within the complete soil moisture–temperature coupling period. Rows (a)–(e) denote regional datasets while columns demonstrate the steps in the coupling analysis, namely negative T_MAX:SM_SURF, positive SM_SURF:EF, and negative EF:T_MAX correlations, each co-occurring with increasing T_MAX. The proportion is calculated at each grid cell as a ratio of the number of days exhibiting the investigated correlation during the complete coupling and the total number of days spent in the respective coupling component during MJJA 2018. The mean proportion over the study region is additionally indicated in each plot.

3.3.5 Coupling at the in situ scale

The SM–T duration throughout the MJJA 2018 period at each of the five stations (Fig. 2), as well as in the corresponding nearest grid cells of GLEAM-E-OBS and Runs 1–4, is presented in Figs. S17–S22. Due to the temporal discontinuities in the stations' data for selected variables, the multi-correlation overlay analysis is based on 92 d, or 10 d as a minimum, running periods. The complexity of factors at the local scale resembles a lack of consistency in the matching periods of complete coupling across datasets. As the analysis is based on daily timescales, it involves frequent fluctuations of the key coupling variables (Findell et al., 2024). There is also a measurement uncertainty; e.g., at the Hyltemossa site, the high-frequency variability of EF throughout a sustained period of anomalously negative EF in July (Fig. S4b) contributed to shortening running window correlations and reducing the complete SM–T duration (Fig. S17).

Besides model resolution and parameterizations, the subgrid hierarchical structure of CTSM is a caveat for SM–T analysis in WRF-CTSM at the local scale. Soil moisture is defined at the soil column level within the vegetated land unit in the grid cell. Multiple PFTs share the same soil column and compete for the available soil water. Surface turbulent fluxes and 2 m temperature are, however, defined on the PFT level (Lawrence et al., 2020). Given the influence of the activities of co-existing PFTs sharing the same SM resource, tracking changes in all key coupling variables at the single PFT level is not straightforward. In addition to vegetation, among the land surface characteristics that affect coupling at the in situ level are elevation and soil texture, defined at the grid cell level in CTSM. Although disentangling the role of specific factors for SM–T strength lies beyond the framework of our study, the land surface characteristics of the investigated stations and their corresponding WRF-CTSM grid cells are summarized in Table S1. Elevation, dominant soil texture, and vegetation types are all represented well in the model at the 0.1° grid spacing for all stations except the Degerö site. This is because Degerö is a boreal mire covering 6.5 km², making it a localized landscape feature surrounded by a pine forest (ICOS Sweden, 2025).

Selecting a single gridded dataset as the most suitable for representing SM–T is not straightforward, as insights about regional-scale patterns in our study primarily rely on a small number of available stations. In GLEAM-E-OBS (Fig. S18), the MJJA 2018 complete SM–T persists longer than at the stations (Fig. S17) or in WRF-CTSM runs (Figs. S19–S22). In contrast to the stations and Run 1, MJJA 2018 is not considered anomalous in GLEAM E-OBS, since the mean MJJA of other available years shows a similar persistence of complete coupling (Fig. S18). A high sensitivity to coupling in GLEAM-E-OBS is also depicted in the mean regional-scale complete coupling duration in MJJA 2010–2022 (excluding 2018), which is 56 d (not shown) as opposed to 32 d throughout the MJJA 2018 (Fig. 6a). One reason may be similar shapes of the key coupling variables over the years, following the seasonal climatology (Fig. S14). Accounting for all the stations, evaluation metrics, and key coupling variables, Run 4 stands out as the most aligned with the in situ observations (Tables S2 and S3). Run 4 incorporates additional parameterizations of the physical processes observed at the in situ scale during MJJA 2018. It includes compensatory root soil water uptake, which was observed even in shallow-rooted species, such as spruce, supported by hydraulic lift from below the rooting zone (Lindroth et al., 2020). Run 4 also reduces the direct impact of SM stress on vegetation (Fig. 9e), partly facilitated by leaf-level regulation of transpiration. It thus provides a more realistic representation of evaporation processes in this mid-to-high latitude environment. Further enhancements are, however, needed as plant water storage capacity and hydraulics are not an inherent part of the stomatal conductance parameterization (Jiménez-Rodríguez et al., 2024; Kennedy et al., 2019; Xu et al., 2023).

4 Conclusions

This study examines the extreme warm and dry conditions in southern and central Sweden during the MJJA 2018 period, where previous work demonstrated that SM–T amplified the heatwave intensity. Our analysis utilizes four WRF-CTSM simulations with varying setups based on recent CTSM developments (Runs 1–4). To complement the WRF-CTSM simulations, which are used here for the first time to analyze sustained warm and dry conditions, we also employ the GLEAM-E-OBS merged observational–reanalysis dataset, which provides all relevant variables at an analogous 0.1° grid spacing. The results of the two main study objectives, evaluating WRF-CTSM and GLEAM-E-OBS key coupling variables (SM_SURF, EF, T_MAX), and identifying the spatio-temporal extent of SM–T and its components during MJJA 2018, are addressed in the following.

Across all evaluation metrics, Runs 2–4 represent SM_SURF well (against station observations and the ESA CCI satellite-based dataset), while GLEAM-EOBS shows a slight overall advantage. EF (based on LH and SH) from the stations is more closely aligned with Runs 2–4, and all datasets virtually reproduce local T_MAX. The spatio-temporal extent of SM–T is assessed using a multi-correlation overlay approach based on 92 d running periods of combined negative T_MAX:SM_SURF, positive SM_SURF:EF (terrestrial leg), and negative EF:T_MAX (atmospheric leg) correlations under increasing T_MAX. On average across the study region and all five datasets, T_MAX:SM_SURF persisted for 72 d, SM_SURF:EF for 42 d, and EF:T_MAX for 25 d during MJJA 2018, while the complete SM–T lasted for 22 d. Therefore, there is a strong but not definitive basis for determining the complete coupling duration based solely on analysing the number of days spent in the atmospheric leg (92 % of the complete coupling duration). The differences in SM–T duration among datasets are most pronounced at the local scale, as land surface characteristics and climatic conditions during the warm period are largely determined by grid spacing and utilized parameterizations. In contrast to GLEAM-E-OBS, WRF-CTSM captures the anomalousness of SM–T in MJJA 2018. Run 4 appears to represent the complete coupling most closely, although the exact timing and duration are difficult to replicate at local scales.

By combining observations and modelling, our study demonstrates the potential of WRF-CTSM for studying extended periods of warm and dry conditions, as it includes a detailed representation of physical processes. This work also highlights the value of the multi-correlation overlay approach in disentangling and examining the interplay between soil moisture and 2 m temperature. However, this analysis is limited to one case study. To ensure a more robust assessment, the same model and approach should be applied to additional events. Finally, these findings underscore the need to address inconsistencies in how SM–T variables and their interactions are represented and parameterized at the in situ scale.