Accurate assessment of land–atmosphere coupling in climate models requires high-frequency data output

Land–atmosphere (L–A) interactions are important for understanding convective processes, climate feedbacks, the development and perpetuation of droughts, heatwaves, pluvials, and other land-centered climate anomalies. Local L–A coupling (LoCo) metrics capture relevant L–A processes, highlighting the impact of soil and vegetation states on surface flux partitioning and the impact of surface fluxes on boundary layer (BL) growth and development and the entrainment of air above the BL. A primary goal of the Climate Process Team in the Coupling Land and Atmospheric Subgrid Parameterizations (CLASP) project is parameterizing and characterizing the impact of subgrid heterogeneity in global and regional Earth system models (ESMs) to improve the connection between land and atmospheric states and processes. A critical step in achieving that aim is the incorporation of L–A metrics, especially LoCo metrics, into climate model diagnostic process streams. However, because land–atmosphere interactions span timescales of minutes (e.g., turbulent fluxes), hours (e.g., BL growth and decay), days (e.g., soil moisture memory), and seasons (e.g., variability in behavioral regimes between soil moisture and latent heat flux), with multiple processes of interest happening in different geographic regions at different times of year, there is not a single metric that captures all the modes, means, and methods of interaction between the land and the atmosphere. And while monthly means of most of the LoCo-relevant variables are routinely saved from ESM simulations, data storage constraints typically preclude routine archival of the hourly data that would enable the calculation of all LoCo metrics.

Here, we outline a reasonable data request that would allow for adequate characterization of sub-daily coupling processes between the land and the atmosphere, preserving enough sub-daily output to describe, analyze, and better understand L–A coupling in modern climate models. A secondary request involves embedding calculations within the models to determine mean properties in and above the BL to further improve characterization of model behavior. Higher-frequency model output will (i) allow for more direct comparison with observational field campaigns on process-relevant timescales, (ii) enable demonstration of inter-model spread in L–A coupling processes, and (iii) aid in targeted identification of sources of deficiencies and opportunities for improvement of the models.

1 Introduction

Much progress has been made in understanding and characterizing land–atmosphere (L–A) interactions in recent years (for an overview of some advances, see Santanello et al., 2018). The importance of L–A interactions has been demonstrated in the initiation, perpetuation, propagation, and termination of droughts (e.g., Otkin et al., 2018; Roundy et al., 2013; Herrara-Estrada et al., 2019; Wu and Dirmeyer, 2020); in the exacerbation of heatwaves (Findell et al., 2017; Alizadeh et al., 2020; Petch et al., 2020; Selten et al., 2020; Seo et al., 2020; Dirmeyer et al., 2021; Benson and Dirmeyer, 2021); and in the timing of monsoon or rainy-season onset (e.g., west Africa: Berg et al., 2017; India: Tuinenberg et al., 2014; the Amazon: Wright et al., 2017). These and other studies collectively suggest the importance of accurately modeling processes at the heart of these feedbacks and interactions. However, output from climate model simulations is rarely saved at high-enough frequencies to capture the rapidly changing features and fluxes that are crucial to the proper characterization of the many links in the chain of L–A interactions (Santanello et al., 2011). These individual linkages include the following:

• the impact of surface temperature, soil moisture, and vegetation on turbulent fluxes at the L–A interface;

• the impact of those fluxes on boundary layer (BL) mixing and moist static energy (MSE);

• the impact of BL processes (e.g., growth rate and buoyancy) on entrainment of air above the BL; and

• their cumulative impact on

  • the BL height, temperature, and humidity and

  • the development of clouds and/or precipitation.

Figure 1 LoCo metrics assess interactions between different parts of the Earth system (y axis) at different temporal (x axis) scales. Yin et al. (2023) and Seo and Dirmeyer (2022) highlight the need to recognize that the two-legged metrics of Dirmeyer (2011) yield results that are dependent on the temporal frequency of the input data, thus requiring a separation between hourly-to-daily (H2D) and monthly-to-seasonal (M2S) versions of the two-legged metrics. Esit et al. (2021) show promising predictability benefits from soil moisture initialization, extending the scope of soil moisture memory into the seasonal-to-decadal time frame. Modified from Santanello et al. (2018).

Figure 1 schematically demonstrates that individual metrics of L–A coupling capture different aspects of these complex linkages. While some metrics focus on the physical processes that operate within the diurnal cycle (e.g., mixing diagrams, Santanello et al., 2009, 2011), others focus on the signal of L–A interactions emerging from long-term multi-variate statistics (e.g., the triggering feedback strength or TFS, Findell et al., 2011). Because of this complexity, we cannot select just one variable, metric, or timescale to assess the strength of a model's coupling between the land and the atmosphere.

The objects in Fig. 1 highlight the distinction between metrics that elucidate physical processes directly (within the diurnal cycle) and those that look at the statistical behavior in data aggregated into long time series using sub-daily, daily, or longer-term mean values in the statistical analyses. Both classes of metrics provide useful information about L–A coupling; when used to inform model development and improvement, the statistical metrics can reveal symptoms of model behavior, while the process-oriented metrics can potentially diagnose causes (see Neelin et al., 2023, for a detailed appreciation and application of process-oriented diagnostics to assess and improve model behavior). For the purposes of demonstrating some of the critical information that can be learned from analyzing observations and models at sub-daily timescales, we will focus on the use of mixing diagrams (Santanello et al., 2009, 2011), two-legged metrics at multiple scales (Dirmeyer, 2011; Yin et al., 2023; Seo and Dirmeyer, 2022), and the triggering feedback strength (TFS, Findell et al., 2011, 2015).

Recent observational field campaigns have included high-frequency observations that can be compared to output from models covering a wide range of purposes and scales (e.g., Earth system models or ESMs, regional climate models, large-eddy simulations, and single-column land–atmosphere models) to test assumptions about L–A behavior. These include the Land–Atmosphere Feedback Experiment (LAFE) at the Southern Great Plains (SGP) site near Lamont, Oklahoma, USA (Wulfmeyer et al., 2018); the Chequamegon Heterogeneous Ecosystem Energy-balance Study Enabled by a High-density Extensive Array of Detectors (CHEESEHEAD) in Wisconsin, USA (Butterworth et al., 2021); and the Land surface Interactions with the Atmosphere over the Iberian Semi-arid Environment (LIAISE) experiment in northeastern Spain (Boone et al., 2021). For example, using high-frequency data from three observational towers from LAFE, Wulfmeyer et al. (2022) demonstrate some of the shortcomings of Monin–Obukov similarity theory (MOST, Monin and Obukhov, 1954) in the estimation of surface fluxes of sensible heat, latent heat, and momentum in unstable conditions. The widespread use of MOST in many model parameterizations speaks to the progress enabled by its implementation. However, the recent acquisition of high-frequency observations like those from LAFE and longer-life-span land–atmosphere feedback observatories (LAFOs) with the same instrumentation (Späth et al., 2023) exposes model shortcomings which can only be evaluated with high-frequency model output. While high frequency in the context of GCMs means something different than in the context of boundary layer turbulence (typically on the order of seconds), the data request presented here will enable evaluation of processes occurring on hourly to 3-hourly timescales, enabling a leap forward in understanding both the processes themselves and ESM representations of those L–A coupling processes.

The spatial scales of individual grid cells in ESM simulations included in the most recent Coupled Model Intercomparison Project (CMIP6) typically range from 50 to 250 km, with models run at resolutions finer than 50 km being eligible for participation in the High-Resolution Model Intercomparison Project (HighResMIP; Haarsma et al., 2016). These resolutions suggest that the footprint sampled from in situ observations (ranging from centimeter-scale soil moisture probes to wind- and height-dependent flux tower sampling fetches on the order of hundreds of meters) is substantially smaller than individual ESM grid cells. This suggests that, when possible, observational comparisons should be made against sub-grid tiles representing fractional areas of differing land use types. However, saving tile-specific high-frequency data is likely not to be feasible for most modeling centers. Given that reality, the data request outlined here will enable the previously impossible assessment of grid cell mean behavior throughout the diurnal cycle. Future work motivated by the CLASP project can extend these lines of inquiry to issues centered on sub-grid spatial heterogeneity or to comparisons with global storm-resolving efforts like those of Stevens et al. (2019).

While short-term simulations saving high-frequency output would allow for a comparison of models with field campaigns, to accurately capture the long-term signal of L–A coupling characterized by the statistically based L–A metrics shown in Fig. 1, sub-daily output of fields at the L–A interface must be saved as part of the routine diagnostic output from long simulations. Furthermore, previous studies have demonstrated that metrics assessing interactions between directly observed variables (e.g., TFS is not directly observed but assesses the relationship between observed fluxes and precipitation) require longer datasets than the directly observed variables themselves (e.g., precipitation) to adequately sample the joint parameter space and compute a statistically robust climatology (Findell et al., 2015).

To assess the coupling strength and details of the interactions in different parts of the L–A system of a GCM, a comprehensive data request would include the following:

• hourly 3D atmospheric profiles of potential temperature (θ), humidity (q), and three-dimensional winds (u, v, w);

• hourly 3D soil profiles of moisture content (SM) and temperature (T_soil); and

• hourly 2D fields of surface pressure; BL height (h_PBL); precipitation (P); sensible heat flux (H); evapotranspiration (ET) and its component parts; near-surface (2 m) temperature, humidity, and winds; net radiation (R_net) fields (incoming and outgoing short- and longwave radiation – SW_down, SW_up, LW_down, and LW_up); and land surface temperature (LST).

The atmospheric profiles should cover the region from the surface to the mid-troposphere in order to capture characteristics of air entrained at the top of the BL. The soil profiles should span from the top of the soil column down through the root zone at minimum. These data would allow for calculation of a host of LoCo metrics, including all but one of the metrics displayed in Fig. 1, at the timescales that are most relevant to the daytime processes the metrics are meant to describe (the GLACE-Ω metric can only be determined with specific model simulations; see Koster et al., 2004). However, we recognize that this would require copious amounts of archive capacity. Here, we aim to reduce that request substantially and include only two-dimensional fields. Our goal is to define a data request which is reasonable in its storage requirements but still provides enough information to characterize the core aspects of the sub-diurnal processes central to L–A interactions. More specifically, the goal is to define a small but sufficient number of data samples per day from two-dimensional fields capturing the sub-diurnal evolution and variability of the following:

• boundary layer properties (BL height and vertically averaged or representative mixed-layer heat content, humidity, and advection)

• fluxes and radiation fields (precipitation, sensible and latent heat fluxes, and net radiation or individual components)

• a bulk measure of stability and humidity deficits above the BL

• root zone and/or surface soil moisture and temperature conditions.

In Sect. 2, we highlight the complexity of the L–A system, showing the many interaction pathways between individual component parts. In Sect. 3, we demonstrate why sub-daily data are required, use these results to provide substantive rationale for the minimum data frequency required to adequately characterize the sub-daily processes of interest, and share an example of the type of behavior that could be routinely assessed if the requested data were regularly made available for model development and/or evaluation. In Sect. 4, we put forth our data request proposal, followed by conclusions in Sect. 5.

2 Highlighting the complexity of the land–atmosphere system

The novel pipe diagrams in Fig. 2 compile linkages as coupling strength indices assessed from daily summertime (June–July–August, JJA) data at the AmeriFlux tower at the SGP field site, along with corresponding diagrams from two versions of the National Oceanic and Atmospheric Administration (NOAA) United Forecast System (UFS) model for the grid cell closest to the SGP site. These coupling strength indices are modeled after the two-legged metrics, named in recognition of the two phases of interaction: the terrestrial leg, which assesses the connection between soil moisture and surface fluxes, and the atmospheric leg, which focuses on the connection between surface fluxes and the BL (Fig. 1; Dirmeyer et al., 2011, 2014). Pipe diagrams from approximately 170 flux tower locations were used during recent model development, aiding the evaluation of UFS prototype 6 (P6) to prototype 7 (P7) (Seo et al., 2023). An advantage of these diagrams is the ability to visualize a host of different L–A linkages at once and thus to identify systematic model biases or behaviors.

Figure 2 Land–atmosphere coupling pipe diagrams calculated from data at the US-ARM Southern Great Plains site (36.61° N, 97.49° W) and UFS model grid cells containing that location, demonstrating the complexity of land–atmosphere interactions, the need for more than one measure to assess coupling, and some of the potential inadequacies of modeled coupling at this example location. Widths of pipes are proportional to coupling index magnitude: $\left|\mathit{\sigma }\right(T\left)r\right(S,T\left)\right|$. Where the sign of the correlation between two terms is the opposite of the expected coupling behavior, dashed blue links indicate severed feedbacks (see text for more information). The three narrow, dashed, faint lines in the bottom panel indicate weak, uncoupled correlations (faint-blue color) and weak coupling index magnitude (very thin lines) for those variable pairings. Daily standard deviations in boxes and coupling index pipes list magnitudes and units; coupling correlations are shown as r=. Adapted from http://cola.gmu.edu/dirmeyer/ufs/P6vP7_loco_chains_AMX.html (last access: 15 December 2023).

The individual coupling strength indices in Fig. 2 are all indicative of both the sensitivity of a target variable T (e.g., latent heat flux) to a source variable S (e.g., soil moisture) and the amount of observed variability:

$$\begin{array}{}\text{(1)}& \mathit{\sigma }\left(T\right)r\left(S,T\right),\end{array}$$

where σ(T) is the daily standard deviation of the target variable, and r(S,T) is the correlation between the two variables. In each pipe diagram, the absolute value of this index is proportional to the width of the link, the strength of the σ(T) term is listed under each variable name and is visually revealed through the intensity of the green color in the box around the variable name, and the strength of the correlation term is enumerated with r= in each link and is visually represented by the color intensity of the link. The physically expected sign of the correlation between each source and target variable is given by a red triangle in the link when a positive correlation is expected (e.g., high soil moisture is associated with high latent heat flux) and by a blue triangle when a negative correlation is expected (e.g., high sensible heat flux is associated with a low evaporative fraction: $\mathrm{EF}=\mathit{\lambda }E/(H+\mathit{\lambda }E)$, where H is sensible heat flux, and λE is latent heat flux, with E being the evaporation rate and λ being the latent heat of vaporization). When the calculated correlation is of the opposite sign compared to this expectation then the variability of the source term is not driving the variability of the target term; thus, the feedback is severed, and the link is represented with a dashed line.

Comparing Fig. 2a and b quickly reveals that, at the grid cell closest to the SGP site, the UFS P6 model exhibits stronger variability in surface fields and stronger coupling between the soils (both moisture and temperature) and the fluxes than is measured at the observational flux tower. Additionally, the modeled fluxes exhibit stronger coupling to 2 m humidity and temperature than the observations show. Though observations have inherent uncertainties from measurement error and issues associated with the representativeness of a single point of the broader region characterized by the model grid cell, this information was used during the model development process, with changes being made to both the land model (Noah in UFS P6 to Noah-MP in UFS P7) and the boundary layer parameterization to improve the full spectrum of coupling strengths manifesting in UFS P7 (details of changes between UFS prototype versions are provided in Stefanova et al., 2022). As a result, the UFS P7 pipe diagram in Fig. 2c is a better match to the observations of Fig. 2a than that of UFS P6 in Fig. 2b.

Pipe diagrams like Fig. 2 can be extended vertically to include additional physical fields and states, accounting for additional links in the LoCo process chain (Santanello et al., 2018). For example, BL properties could include average BL potential temperature, humidity, or moist enthalpy; BL height; and the height or pressure of the lifted condensation level (LCL). A final layer at the top of these pipe diagrams could include information about clouds and precipitation. The myriad of possible links in the process chains connecting individual elements within these pipe diagrams – and, indeed, within the physical land–atmosphere system – demonstrates the complexity of interactions between the land and the near-surface atmosphere. Figure 2 demonstrates that model parameterizations influence the modeled strength and connectivity of different parts of the L–A system and that confronting models with process-level observations from different climate regimes can help expose model deficiencies and limitations. For that to be possible, however, model output must be temporally equivalent to the observations on hand, and it must adequately sample the behavior of the physical processes of interest. While daily data were successfully leveraged to improve land–atmosphere coupling in the UFS model, the next section demonstrates some of the processes requiring sub-daily data.

3 Establishing the need for high-temporal-resolution data

The triggering feedback strength (TFS, Findell et al., 2011) is a measure of the sensitivity of afternoon rainfall occurrence to morning-time evaporative fraction ($\mathrm{EF}=\mathit{\lambda }E/(H+\mathit{\lambda }E)$). Using 3-hourly data from the North American Regional Reanalysis (NARR; Mesinger et al., 2006), Findell et al. (2011) showed that high morning EF enhances the probability of afternoon rainfall east of the Mississippi and in Mexico, with higher EF leading to increases in afternoon rainfall probability of between 10 % and 25 % in these regions. By contrast, the intensity of rainfall was shown to be largely insensitive to surface flux partitioning, as assessed by the amplification feedback strength (AFS; Findell et al., 2011).

A follow-up study by Berg et al. (2013) showed that the Geophysical Fluid Dynamics Laboratory (GFDL) model AM2.1 exhibited similar sensitivity of afternoon rainfall likelihood to morning surface flux partitioning in the eastern US and Mexico and a similar insensitivity of rainfall intensity to surface flux partitioning. However, the similar TFS results from AM2.1 and NARR occurred for different reasons. Like the two-legged metrics discussed above (Eq. 1), the TFS is computed with a sensitivity term (the sensitivity of the probability of afternoon rain to variations in morning-time EF) multiplied by a standard deviation term (σ_EF). In contrast to the two-legged metrics, however, the calculation of the TFS is a summation of purposefully binned or segmented data to account for the possibility of non-uniform sensitivities in different EF regimes; indeed, sensitivity strength is substantially larger at EF >0.6 than at smaller EF values (Findell et al., 2011). Berg et al. (2013) showed that the regions with high TFS values in AM2.1 were driven by larger EF variability (peak σ_EF values of 0.2 in NARR compared to 0.4 in AM2.1), while regions with high TFS values in NARR were driven by larger mean rainfall sensitivities (peak mean sensitivities above 2 in NARR compared to less than 1 in AM2.1). The large values of σ_EF in the AM2.1 results also explained an additional region of high TFS values in AM2.1 in the northern central Great Plains of the US, extending into adjacent areas in southern Canada.

Figure 3 (a) The triggering feedback strength (TFS; units of probability of afternoon (noon to 18:00 local time (LST)) rain) in summer (JJA) based on ERA5 hourly data from 1991 to 2020. The TFS algorithm follows Findell et al. (2015) but with a 10-bin segmentation of daily evaporative fraction (EF). Positive values indicate that the morning EF positively affects the probability of the occurrence of afternoon precipitation. (b) The mean value of the sensitivity term contributing to the TFS $\left(\stackrel{\mathrm{‾}}{\partial \mathrm{\Gamma }\left(r\right)/\partial \mathrm{EF}}\right)$ – thus the mean sensitivity of afternoon rainfall to morning-time surface flux partitioning. (c) The variability term contributing to the TFS: the standard deviation of EF.

Figure 3 shows the June–July–August TFS (panel a) and its two component parts (panels b and c) calculated from hourly European Centre for Medium-Range Weather Forecasts 5th reanalysis data (ERA5; Hersbach et al., 2018, 2020). Comparison with Findell et al. (2011) and Berg et al. (2013) shows that NARR, ERA5, and AM2.1 exhibit the same range of sensitivity of afternoon rainfall triggering to morning-time flux partitioning, but in the ERA5 data, the peak TFS values of 15 %–25 % only manifest in Mexico, with some extension into the southern part of the mountainous US southwest. While the eastern US region shows up with relatively elevated component contributions in ERA5 (Fig. 3b–c), the resultant TFS values are only 5 %–10 % in most of the eastern US and approach 15 % in much of Florida (Fig. 3a). The individual terms contributing to ERA5's TFS results have peak values matching the smaller EF variability of the NARR data rather than the high variability of AM2.1 (Fig. 3c here and Fig. 6 of Berg et al., 2013) and sensitivities matching the smaller AM2.1 values rather than those of the NARR data (Fig. 3b here and Fig. 7 of Berg et al., 2013). These differences across reanalysis datasets are likely impacted by differences in the data assimilation protocols and observational datasets ingested by ERA5 and NARR. In addition, the TFS may also be highly sensitive to each system's parameterizations of the surface layer, boundary layer, and convection since the surface fluxes at the heart of the TFS are not assimilated variables but are wholly model dependent (Kalnay et al., 1996). Additional investigation is necessary to better understand these differences between the reanalyses and the model, but this behavior can only be exposed with analysis of sufficiently high-frequency data. Here, data frequencies of at least 3 h were essential to enable the separation of morning-time fluxes and afternoon precipitation events.

Figure 4 Two-legged metric analysis demonstrating the impact of different averaging periods on the assessment of coupling strength in summer (JJA and DJF for the Northern Hemisphere and the Southern Hemisphere, respectively). The diagnoses are based on the ERA5 (ECMWF reanalysis 5) reanalysis data from 1991 to 2020. The coupling strength between the sensible heat flux and p_LCL is estimated by the TLM algorithm (Dirmeyer et al., 2006). Strongly coupled regions (top 15 % of land grid cells) are diagnosed by using different time series (i.e., D: daytime-only mean, E: 24 h entire-day mean, and M: monthly mean). The Euler diagram is employed to illustrate the spatial differences between the three diagnoses. The areas of colored components in the Euler diagram are proportional to the sizes of specific sets. Modified from Yin et al., 2023.

While a paucity of high-frequency data has forced many previous analyses of two-legged metrics (Eq. 1) to rely on monthly mean data (e.g., Dirmeyer et al., 2014; Hu et al., 2021; Lorenz et al., 2015), Yin et al. (2023) highlight the need to recognize that the two-legged metrics yield results that are dependent on the temporal frequency of the input data (Fig. 4 and the H2D-M2S distinctions in Fig. 1), in part because the magnitude of variability is dependent on the averaging period of the data being analyzed and in part because the inclusion of nighttime hours can mask the daytime feedbacks that are at the heart of the sensitivity between the variables of interest. Figure 4 shows that the assessment of the strength of the atmospheric leg measuring the impact of sensible heat flux H on BL growth (as assessed by the pressure of the lifting condensation level p_LCL) can be very different when using monthly (M), 24 h entire-day (E), or daytime-only (D, 07:00 to 15:00 local time (LST)) time series. Different averaging periods of the input data effectively allow one to ask different questions about coupling: monthly averaged data tell us about the seasonal variability of the terms being assessed and their coupling, while daytime-only data are needed to tell us about the direct impact of surface fluxes on BL properties, for example. In regions where the month-to-month variability is small (e.g., where mean H and p_LCL values are similar for all summer months), substantial day-to-day variability in these terms will not be captured by monthly mean values (e.g., orange regions in Fig. 4). However, in regions where the progression into deeper days of summer tends to bring drier and drier conditions, differences across summer months (e.g., June compared to August) can be substantial, so monthly mean time series will still show high variability and potentially result in a diagnosis of a large coupling strength (e.g., blue regions in Fig. 4). Comparing daily to sub-daily scales, Fig. 4 shows about 30 % disagreement in the highlighted regions with strong H–p_LCL coupling determined from E versus D time series. The nighttime component of E was shown to obscure the diurnal coupling signal in some areas, with complications caused by regionally specific mechanisms (particularly in the very arid regions adjacent to the Mediterranean Sea) or UTC-based time smoothing (Yin et al., 2023). These differences highlight the need for sub-daily data to accurately capture the process-level connections between surface fluxes and the BL response.

Figure 5 (a) The hourly mixing diagrams at water-limited (red) and energy-limited (blue) flux tower sites exhibit the coevolution of moisture (x axis) and thermal (y axis) energy content per unit mass within the BL (modified from Fig. 5d in Seo and Dirmeyer, 2022). The marks are shaded by the color determined by two-legged couplings corresponding to the local hour (referring to Fig. 5a in Seo and Dirmeyer, 2022). The black-edged circle and square are the mean of the 24 h values in water- and energy-limited regimes, respectively. The colored numbers are the area within the curves (multiply displayed value by 10⁶; units: J kg⁻²); these values quantify the diurnal energetic asymmetry captured by each mixing diagram. Dashed black lines are lines of constant relative humidity. Note that x- and y-axis ranges differ. (b) The 3 h mixing diagrams in both climate regimes, computed with three different starting times: hour 0:00 (S0: solid), hour 01:00 (S1: dashed dot), and hour 02:00 (S2: dashed) LST. (c) The 6 h mixing diagrams as in (b) but with starting times at hours 00:00, 02:00, and 04:00 LST.

Seo and Dirmeyer's (2022) thorough evaluation of the hourly evolution of BL temperature and humidity at flux tower observational sites can be leveraged to determine the minimum number of data points needed per day to adequately capture both the thermal and the moisture evolution of the BL. Figure 5a shows hourly mixing diagrams spanning all hours of the day based on Seo and Dirmeyer (2022), plotting moisture (x axis) and heat (y axis) energy content per unit mass within the mixed layer, averaged across the 10 % of the 230 stations that were the most moisture-limited (red circles) and the most energy-limited (blue squares; see Fig. S1 in the Supplement for a global map with station locations). Through their detailed analysis, Seo and Dirmeyer (2022) highlight differences in the timing of the BL response to moisture fluxes compared to heat fluxes, with the thermal process chain often leading ahead of the moist process chain by 2–3 h during the day and rapid thermal decoupling in the late afternoon contrasted with a gradual decline in moist coupling throughout the evening hours. They also highlight dependence of the timing of humidity minimums on moisture availability: Fig. 5a shows that the driest time for the BL is during the early afternoon in moisture-limited regimes but before sunrise in energy-limited regimes. Both moisture- and energy-limited regions show a morning-time peak in BL humidity (07:00–09:00 LST).

Findell et al. (2017) showed that some of these behaviors can be captured in a statistical sense using monthly mean diurnal cycles of temperature and moisture, but a full step-by-step understanding of these detailed processes and interactions requires many data points per day. Figure 5b shows that 3 h output generally captures the critical phases and the maximum extent of the diurnal excursions in T−q phase space, as well as the bulk of the diurnal asymmetry of the T−q evolution. The numbers to the right of each mixing diagram quantify the area within the curve (e.g., 8.26×10⁶ J² kg⁻² in the water-limited diagram of Fig. 5a compared to 1.23×10⁶ J² kg⁻² in the energy-limited composite) and make it clear that, while the 3 h mixing diagrams underestimate the diurnal asymmetry, the process-relevant distinction of small asymmetry in energy-limited regimes compared to large asymmetry in water-limited regimes remains clear. While 6-hourly data (Fig. 5c) can capture the approximate timing of the humidity minimums (late afternoon versus early morning), such infrequent sampling can miss the most rapidly changing portions of the daytime T−q evolution (e.g., samples beginning at 00:00 LST), leading to inaccurate assessments of the extent of the diurnal asymmetry in T−q energetic phase space.

4 Justifying our choices on how to reduce the data request

4.1 Strategy for the reduction in time frequency

To determine the optimal strategy for reducing the time frequency of the data request while still achieving the coupling-assessment goals discussed above, we consider two possible strategies: (i) regular, gridded time intervals or (ii) time intervals based on the local solar day (e.g., values for nighttime, morning, and afternoon). Positive arguments for the first approach include the lack of subjectivity and the ease of implementation. Counterarguments center around geographic differences imposed by the gridded approach. For instance, in the summertime, sunrise times along one longitudinal band differ by about 2 h between high-latitude regions and the Equator. Thus, a 05:00 LST data point on the summer equinox would be after sunrise in Sweden but before sunrise in the Congo, even though both have a longitude of 18° E. This poses difficulties for investigations of, for example, the triggering feedback strength (TFS), meant to capture the impact of early-morning evaporative fraction on afternoon precipitation (Findell et al., 2011). Capturing specific times of day becomes more complicated with a reduced frequency of data collection or archival. Since hourly data represent 15° longitudinal bands around the globe, coarser-frequency data inherently require grouping broad longitudinal slices into common time points. For 6 h data, an attempt to capture early-morning conditions within one 90° longitudinal slice would produce local times that might span six time zones, potentially ranging from, for example, 03:00 LST at the western edge to 08:00 LST on the eastern edge. Clearly, processes at the land surface and within the boundary layer differ substantially between these times of day. Higher-frequency data would reduce the severity of these issues, albeit with more archive space required.

Positive arguments for a data-archiving scheme linked to the solar day include fewer data points (and thus less archive capacity) being needed to capture the three main phases of BL behavioral regimes (nighttime, morning, afternoon) and a more uniform understanding of the solar conditions associated with each data point. However, any sub-daily selection based on the solar day requires a priori decisions that might be appropriate for one purpose but which would restrict the appropriateness of further study. For example, mixing diagrams are useful tools to understand BL evolution within each of the three solar-day phases mentioned above. Saving average values within these three phases would eliminate the possibility of any sort of mixing diagram analysis of model behavior. Additionally, interpretation of solar-day-based data would be complicated by each archived data point representing different numbers of hours, both from day to day at one location and from location to location on each day. Furthermore, this strategy would require additional code being written and implemented at each climate modeling center, and thus, the possibility of differences in implementation quickly emerges.

Here, we opt to make a request for regular, temporally gridded data to avoid the complications of solar-day-based archiving and to maintain flexibility for future data usage. The negative features of the regular, gridded temporal data requests can be reduced with increased frequency of data storage. We propose 3 h data as a minimum request, with 1 or 2 h as improvements to that minimum. If this data request is still too cumbersome, a mask of oceanic regions can potentially be used to reduce the data volume by up to two-thirds, though these data may be useful for the study of ocean–atmospheric boundary layer coupling processes.

4.2 Other issues to confront

In addition to decisions related to the reduction in the time frequency of data archiving, our data request must tackle difficult decisions related to (i) capturing mean BL properties while the height of the BL is changing, (ii) capturing adequate measures of the temperature and humidity gradients above the BL, and (iii) capturing soil conditions (moisture and temperature) most relevant to the partitioning of energy into surface fluxes of latent and sensible heat.

Determining average properties within the BL at any given time requires knowledge of the height of the BL (h_PBL). Model-computed h_PBL is determined using different methods in different models, producing values which are self-consistent within each model's framework, and therefore should adequately capture the time evolution of the BL height at a given location and the relative BL heights at different locations. However, nighttime values of BL average properties will necessarily represent something different than daytime values, and values during the transition times of day will be tricky to compute and difficult to rely on. In addition, these times of day will change throughout the year. All of these issues suggest that care is needed in implementing these calculations and interpreting the results.

Finally, for characterization of soil conditions most relevant to surface energy partitioning, a root zone soil moisture would be most appropriate. However, since the root zone is both dynamic and dependent on vegetation type, no single depth can adequately capture the true root zone. Here, we opt for a near-surface measure of the top 10 cm plus a slightly deeper measure averaged over the 10–100 cm interval. In both cases, we recognize that these are characterizations of the model's soil wetness but that this variable is a model-specific quantity that is different from in situ or remotely sensed measures of soil wetness and should be interpreted with recognition of the model value's mean and variability (e.g., Koster et al., 2009; Benson and Dirmeyer, 2023).

Table 1 Specifics of request A. Grid cell average values are either 3 h time means (3 h) or are at an instantaneous point in time at the end of the time interval (3 hPt).

The ^* symbol indicates variables without standard Climate Model Output Rewriter (CMOR) names.

Table 2 Specifics of request B. The six requested pressure levels are every 50 hPa between 950 and 700 hPa. Grid cell average values are instantaneous in time at the end of the time interval (3 hPt).

5 The data request

Here, we present a concrete data request, dividing the request into three categories based on the analyses that would be enabled and by the work required by model developers. Request A is the highest-priority request and focuses on standard model output of surface fields saved at higher-frequency intervals than is currently routinely practiced, thus requiring no additional work by model developers, just additional archive space. This request includes both tier-1 and tier-2 variables. Request B focuses on the archival of variables in the lowest 300 mb of the troposphere and includes only tier-2 variables. Like request A, request B requires no additional work by model developers, just additional archive space, while request C requires in-code modifications to calculate average properties within and above the BL (tier-3 variables). After each request, we briefly mention which metrics (mostly from Fig. 1) and analyses would become possible with these additional data.

The data length requirements of Findell et al. (2015) suggest that a minimum of 10 years of data should provide for robust statistical analyses. Thus, for any simulation and/or time period of climatological interest, we request that these data are saved for at least a 10-year block of time. For historical and future scenario runs, it would be advantageous to have 10-year blocks saved at the beginning and end of the simulations.

5.1 Request A: high-frequency archival of surface variables already included in standard model output

Table 1 details the variables included in request A. The 10 tier-1 variables would allow for the computation of several two-legged metrics at sub-daily timescales (including all of those included in Fig. 2), soil moisture memory, TFS and AFS, basic mixing diagrams, and the percentile soil moisture–aridity index framework of Duan et al. (2023). Assuming a 1° grid (for reference) without data compression, archiving of tier-1 variables would require approximately 13 GB yr⁻¹.

Request A also includes several tier-2 priority variables: deeper soil moisture information and component terms of net radiation and evapotranspiration. These additional terms would allow for a more in-depth understanding of model depictions of radiative processes and of the role of vegetation in driving evaporative fluxes and feedbacks. However, they would nearly double the required archival requirements and, thus, have been deemed to be tier-2 priority variables.

Of the 10 tier-1 variables listed in Table 1, the first eight were included at a 3 h frequency in the HighResMIP data protocol (Haarsma et al., 2016), with soil temperature (tsl) saved at a 6 h frequency and boundary layer depth (bldep) saved monthly. HighResMIP also included 16 other variables in their 3 h data request (for a total of 24 3 h variables), indicating that saving all of the request-A variables is not an insurmountable challenge.

Table 3 Specifics of request C. The BL mean properties should be vertically integrated from 0.1×h_PBL to 0.8×h_PBL, while the gradients across the BL top should be calculated over the interval 0.8×h_PBL to 1.2×h_PBL. CMOR names are not currently available for these quantities.

5.2 Request B: high-frequency archival of data at several specified lower-tropospheric pressure levels

Table 2 details the five variables included in request B for archival of select lower-tropospheric pressure levels, specifically temperature, humidity, and three-dimensional winds. The priority here is to enable systematic exploration of BL processes throughout various stages of growth, development, and decay. Saving high-frequency data of full atmospheric profiles is not realistic, but saving a few select pressure levels would allow for the computation of atmospheric stability and humidity deficits in the early-morning hours (i.e., metrics like CTP and HI_low); mean properties within the BL; d$\mathit{\theta }/\mathrm{d}z$ and $\mathrm{d}q/\mathrm{d}z$ above the BL; the heated condensation framework; and more complex mixing diagrams than request A would enable, including identification of advection and entrainment terms during multiple phases of BL growth and development. The six specific pressure levels requested are every 50 hPa between 950 and 700 hPa.

5.3 Request C: variables requiring code modifications for internal computation

With request C, we aim to enable more accurate mixing diagram work than is possible with request B while simultaneously reducing the archive requirements needed to assess mean properties within and above the BL. Request C entails code modifications to determine, at each time step, the BL mean thermal and moist energy content per unit mass (c_pθ and λq, respectively), changes in these terms due to advection, and the mean potential temperature and humidity gradients across the top of the BL given by h_PBL (or the CMOR variable name bldep in Table 1). For a standard definition of h_PBL, we suggest the bulk Richardson number definition of Seidel et al. (2012), consistent with the data available in reanalyses such as ERA5 and MERRA2. Specifically, we request the mean BL properties vertically integrated from 0.1×h_PBL to 0.8×h_PBL and the mean θ and q gradients over the interval closest to 0.8×h_PBL to 1.2×h_PBL, given model-level constraints (see Turner et al., 2014, for the selection of these vertical bounds). These properties should be saved every 3 h.

While request C would reduce the archive requirements for mixing diagram work and provide a fuller picture of mean mixed-layer behavior, it would not allow for some of the other metric calculations that request B does cover. Thus, these are complementary requests rather than substitutes for each other.

6 Conclusions

Increasing the time resolution of model output describing components of land–atmosphere coupling and processes within the land–atmosphere interface is essential to fully and accurately model, understand, and predict these processes and to compare modeled processes with observational datasets. The data request described here will allow us to compare coupled Earth system and climate models with observations from field campaigns and to compare both diurnal and long-term properties of L–A interactions in different models and during model development. These sorts of comparisons are essential to fully assess the land–atmosphere coupling behaviors of different GCMs. Furthermore, these improvements to our understanding of processes at the land surface are essential to understanding the vulnerability of humans and ecosystems to changing climatic conditions and to improving our resilience in the face of a likely increase in extremes.