## Abstract

Evaporation plays an extremely important role in determining summertime surface temperature variability over land. Observations show the relationship between evaporation and soil moisture generally conforms to the Budyko “two regime” framework; namely, that evaporation is limited by available soil moisture in dry climates and by radiation in wet climates. This framework has led climate models to different parameterizations of the relationship between evaporation and soil moisture in wet and dry regions. We have developed the Simple Land–Atmosphere Model (SLAM) as a tool for studying land–atmosphere interaction in general, and summertime temperature variability in particular. We use the SLAM to show that a negative feedback between evaporation and surface temperature gives rise to the two apparent evaporation “regimes” and provide analytic solutions for evaporative cooling anomalies that demonstrate the nonlinear impact of soil moisture perturbations. Stemming from the temperature dependence of vapor pressure deficit, the feedback we identify has important implications for how transitions between wet and dry land surfaces may impact temperature variability as the climate warms. We also elucidate the impacts of surface moisture and insolation perturbations on latent and sensible heat fluxes and on surface temperature variability.

## 1. Introduction

As the climate warms, the impacts of increasing summertime temperatures are becoming more evident. Global warming has increased the severity of maximum monthly summertime temperatures (Kirtman et al. 2013; Diffenbaugh et al. 2017). The temperature of the hottest day of the year has increased nearly 1°C per decade over the past 30 years in both Houston and Moscow, while the average trend over Eurasia is 0.3°C per decade (Papalexiou et al. 2018). Small alterations to climate variability result in large changes to the probability of extreme events like heat waves and droughts that depend on temperature thresholds (Katz and Brown 1992), and constitute a major challenge to climate adaptation. These challenges and their consequences were on display in 2003, when an unprecedented seasonal heat wave killed 70 000 people in Europe (Robine et al. 2008). Schär et al. (2004) found that this seasonal anomaly could not be explained by mean global warming and invoked increased temperature variability to explain the extremely unlikely heat wave. Just seven years later, an even more extreme seasonal heat wave killed 55 000 people in Russia and likely broke 500-year temperature records over much of central and eastern Europe (Barriopedro et al. 2011; Tingley and Huybers 2013).

In addition to the mortality associated with heat waves, temperature extremes impact global food security by contributing to year-to-year uncertainty in crop yields. The 2010 Russian heat wave noted above caused a 25% reduction in wheat yield, and Tigchelaar et al. (2018) found that mean global warming dramatically increases crop yield variability. Specifically, a 4°C mean global warming makes the likelihood that each of the four largest maize-producing countries will experience a 10% or greater drop in yield during a particular year 87%, while in the current climate this kind of synchronized shock is extremely unlikely.

Given the impacts of extreme summertime temperatures and the dependence of statistical extremes on natural variability, understanding the physical processes that give rise to summertime temperature variability is extremely important. For over a decade, regional models forced with anthropogenic CO_{2} emissions have projected a robust increase in temperature variability over Europe (Vidale et al. 2007). In particular, models that simulate the modern European climate most accurately project a 20% increase in daily temperature variability over large swaths of southern and central Europe by the end of the twenty-first century under business as usual greenhouse gas emissions (Fischer et al. 2012). Global climate models forced by increasing greenhouse gas emissions project significant (up to 40%) increases in the standard deviation of monthly summertime temperatures, particularly over North and South America, Europe, and Southeast Asia (see Fig. 1 and Table 1 for a list of models in the ensemble).

Our confidence in these projections is tied to the ability of the models to represent the present-day summertime climate. Unfortunately, regional and global climate models have large biases in continental land surface temperatures, with summer mean temperatures that are too high and more variable than the historical record (Lenderink et al. 2007; Morcrette et al. 2018). Further, Donat et al. (2017) found that changes in extremely warm temperatures in climate model simulations of the past 50 years are much larger than corresponding changes in the observational record. The source of these errors differs across models, and due to the complexity of land–atmosphere interaction no singular cause has been identified. However, model representations of evaporation have been implicated as a primary cause of temperature biases in climate models (Mueller and Seneviratne 2014; Merrifield and Xie 2016; Ma et al. 2018). Variations in evaporation are linked to variations in soil moisture, which has become widely regarded as a fundamental control on summertime temperatures (Seneviratne et al. 2010).

While current observations cannot reveal the mechanistic relationship between soil moisture and evaporation (Koster et al. 2015), observational studies have shown two distinct patterns of evaporation behavior (Ryu et al. 2008; Teuling et al. 2009; Gentine et al. 2012). The conceptual framework underpinning these two patterns comes from Budyko (1961). Budyko proposed that *evaporation efficiency* (defined as the ratio of actual to potential evaporation) is linearly dependent on soil moisture up to a critical value, above which it is constant (see Fig. 2). Below this critical value evaporation is considered “moisture controlled” while above it is considered “climate (or energy) controlled” (Eagleson 1978). Even as models of the land surface have become more complex, the nonlinear connection between soil moisture and evaporation proposed by Budyko is commonly invoked to explain model output and observations of evaporation [see discussion in Seneviratne et al. (2010)]. However, no critical value of soil moisture has been observed (Koster et al. 2006), begging the question of what really distinguishes the two apparent soil moisture “regimes.”

Despite the conceptual power of Budyko’s framework, Dirmeyer et al. (2006) showed that an ensemble of climate models had no consistent representation of the connection between soil moisture and evaporation. The common invocation of the two soil moisture regimes proposed by Budyko contrasts starkly with the murky representation of evaporation behavior in climate models. Developing a process-based account of these apparent soil moisture regimes is therefore critical to understanding the thermodynamics that govern land–atmosphere interaction and summertime temperature variability.

In this paper, we aim to illuminate the essential physics underpinning the relationship between soil moisture and evaporation. To achieve this, we construct the one-dimensional Simple Land–Atmosphere Model (SLAM) that represents the land–atmosphere fluxes of energy and moisture with as few parameters as possible. In section 2, we will briefly describe the SLAM and show results from an evaluation exercise. An in-depth description of the model’s derivation and equations is given in the appendix. The SLAM relies on time series of forcing and boundary conditions that are not widely available on subseasonal time scales; in section 3 we develop a suite of synthetic forcings that we use to create an ensemble of model runs with varying surface moisture initializations. In section 4, we show that the two apparent soil moisture regimes are a fundamental feature of land–atmosphere interaction caused by a feedback between evaporative cooling and vapor pressure deficit modulated by available soil moisture. In section 5 we present a discussion of results and conclusions.

## 2. SLAM description and evaluation

### a. Model description

Figure 3a shows a schematic of the SLAM where all model layers and thermodynamic variables are labeled. All model layers have a temperature (K) and moisture variable that are assumed to be homogenous within the layer. The moisture variables are specific humidity *q* [kg H_{2}O (kg air)^{−1}] for the atmospheric layers and volumetric soil water content *m* [m^{3} H_{2}O (m dry soil)^{−3}] for the soil layers. The volumetric soil water content is defined as the volume occupied by liquid water in a unit volume of soil. Smith and Mullins (2001) use the mass ratio of liquid water to dry soil *w* [kg H_{2}O (kg dry soil)^{−1}] to define volumetric soil water content:

In Eq. (1), *ρ*_{d} and *ρ*_{l} are the bulk densities (kg m^{−3}) of dry soil and liquid water. In practice, *m* has an upper limit *m*_{sat} set by the volume of air-filled pore space in a particular dry soil sample.

The soil surface layer has thickness *d*_{2} = 10 cm; this encompasses the region where radiative heating is largest during the diurnal cycle. The depth of the total soil column *D*_{r} = 1 m contains the full rooting zone and is deep enough that we can assume a constant temperature below this level; the constant “cold abyss” temperature *T*_{D} is one of the model’s boundary conditions. The lower atmospheric boundary layer thickness *d*_{3} = 100 m incorporates nearly all possible vegetation types, while the total boundary layer height *D*_{a} = 1000 m is a representative value over land surfaces. Temperature *T*_{Top}, specific humidity *q*_{Top}, and cloud fraction *c*_{f} boundary conditions are prescribed above the upper boundary layer.

All model fluxes are shown in Fig. 3b; model forcings and fluxes controlled by boundary conditions are shown in purple. Within the four model layers, fluxes of moisture (kg H_{2}O m^{−2} s^{−1}) and energy (W m^{−2}) are shown in blue and orange respectively. The term $P$ is the precipitation into the surface soil layer, while transpiration from the root and surface layers are *η*_{1} and *η*_{2} respectively. The term *Q*_{1,2} is the infiltration of liquid water between the surface and root layers, while *Q*_{↓} is drainage of excess soil moisture out of the model. Terms *Q*_{3,4} and *Q*_{↑} are turbulent fluxes of water vapor between the lower and upper atmospheric boundary layers, and between the upper boundary layer and free troposphere, respectively. The term $R$ is the net shortwave radiation absorbed at the surface, while *F*_{2}, *F*_{3}, and *F*_{4} represent net longwave radiation absorbed in the layer denoted by the subscript. Also, *H*_{2,3}, *H*_{3,4}, and *H*_{↑} are the turbulent sensible heat fluxes from the soil surface to the lower boundary layer, between the lower and upper boundary layers, and between the upper boundary layer and free troposphere, respectively. Terms *H*_{1,2} and *H*_{↓} represent conductive heat fluxes that transfer energy from the surface to the root layer and from the root layer out of the model, respectively. Finally, *E* represents evaporation; because evaporation both cools and dries the surface, it is shown in both blue and orange in Fig. 3b.

A detailed description of the model’s derivation and all equations governing the model fluxes is given in the appendix, but the evaporation *E* (kg H_{2}O m^{−2} s^{−1}) equation is of fundamental importance, so we will present it here:

Here *r*_{s} (s m^{−1}) is the surface resistance that governs the rate at which energy is transferred from the land surface into the lower boundary layer by turbulent eddies, while *f*_{υ} is the fraction of the land surface covered by vegetation. The density of air is *ρ*_{a}, *q*_{s} [kg H_{2}O (kg air)^{−1}] is the saturation specific humidity evaluated at the surface temperature *T*_{2}, and *q*_{3} is the specific humidity in the lower boundary layer.

In Eq. (2), *X* is a coupling term that connects the soil water content in the surface layer *m* to evaporation. In the Budyko framework, this coupling term is nonlinear (see Fig. 2):

In Eq. (3), the wilting point *m*_{w} is the point at which evaporation becomes impossible, and *m*_{crit} is a parameter that controls the distinction between the apparent moisture-controlled and climate-controlled regimes.

To investigate whether the Budyko parameterization is crucial to the development of these apparent regimes, we assume a simpler form of the coupling term:

Under this assumption, *X* varies linearly between zero and one and can be thought of as the fractional soil saturation (it is expressed as a percentage in our figures). We constrain soil moisture in each layer so that it does not exceed saturation *m*_{sat} or go beneath the wilting point *m*_{w}. While certainly a simplification of the complex processes that relate soil moisture to evaporation, we will show that this simple coupling deployed in the SLAM generates realistic evaporation behavior across the soil moisture spectrum (see section 4).

### b. Model evaluation

The SLAM needs time series of net absorbed shortwave radiation $R$, precipitation $P$, and cloud fraction *c*_{f} as well as temperature and humidity boundary conditions *T*_{Top} and *q*_{Top} to generate output. The Atmospheric Radiation Measurement facility in the Southern Great Plains (SGP) provides data on solar radiation $R$, precipitation $P$, and cloud fraction *c*_{f} every minute; these data are shown in Figs. 4a–c (Riihimaki and Shi 1994; Long et al. 2014). Boundary conditions *T*_{Top} and *q*_{Top} come from interpolated radiosonde data for temperature and specific humidity at 1000 m that are provided at every minute and are shown in Figs. 4d and 4e (Troyan and Jensen 1998). Data for these five fields are available from 0000 CDT 1 June 2014 through 0000 CDT 31 August 2014.

In addition to the time series shown in Fig. 4 we used the parameters listed in Table 2 for our simulation. We refer the interested reader to the appendix for information on how these parameters are incorporated into model equations. All parameters on the left-hand side of the table are held constant and govern the turbulent energy and moisture fluxes. Vegetation fraction *f*_{υ}, and stomatal resistance *r*_{st} are taken as representative values for an average grassland (Wei et al. 2017; Garratt 1992). Surface resistance *r*_{s} and maximum resistance to turbulent heat transport $ra,T\xaf$ are estimates for a reasonably smooth land surface (Garratt 1992). The maximum resistance to turbulent moisture transport $ra,q\xaf$ is much smaller than the equivalent resistance to turbulent heat transport to keep water vapor well mixed in the boundary layer. We estimate the average insolation maximum *σ*_{max} from the radiation observations (see Fig. 4a).

The right-hand side of Table 2 contains soil parameters. Deep root fraction *f*_{r} for a grassland is estimated from Jackson et al. (1996). Bulk land surface density and heat capacity vary linearly between the values indicated in the table (see the appendix). Tong et al. (2016) present functional forms for how soil conductivity *λ* increases with soil moisture.

The SLAM output is insensitive to changes in initial conditions for temperature and specific humidity in the various layers. For simplicity, we take initial conditions directly from the observations. The constant temperature of the cold abyss *T*_{D} was set to 280 K, and the SLAM output is also insensitive to changes in this value.

Soil moisture initial conditions do not have a large impact on model output, but this is due to the 100-mm precipitation event that occurs during the first week of summer (hence, in the first week of the simulation; Fig. 4b). This event eliminates any influence of the soil moisture initial condition on the rest of the simulation, as it effectively saturates both soil layers in the model. Still, the soil moisture observations from the SGP site start near their minimum value, so we prescribe this “dry start” for our simulation in both the surface and root layers.

In Fig. 5, we show observations from the SGP site in blue and the SLAM output in orange. Tower observations of temperature and specific humidity at 60 m are sampled once per minute (Xie and Chen 2012) and compared to *T*_{3} and *q*_{3} output from the SLAM. Sensible and latent heat flux observations derived from eddy covariance estimates with 30-min time steps (McCoy et al. 2017) are compared to *H*_{2,3} and *L*(*E* + *η*_{1} + *η*_{2}) output from the SLAM. Volumetric soil water observations were gathered from 5-cm depth once per hour, normalized using Eq. (4) and the measured maxima and minima of the observations for the corresponding *m*_{sat} and *m*_{w} parameters, and compared to *X* output from the SLAM (Ermold et al. 1996). To quantify model performance, Fig. 5 also shows the variance in daily averaged values of each observed time series explained by daily averaged SLAM output (*r*^{2}), and the ratio of standard deviation in the daily averaged SLAM output to that observed (*σ*_{M}/*σ*_{O}).

The SLAM output explains more than 60% of the variance in observed daily averaged 60-m temperature, and the standard deviation in the SLAM’s daily averaged temperatures closely matches that observed. Similarly, the variance in daily averaged SGP 60-m specific humidity is largely explained by the SLAM output, while the standard deviation in SLAM’s daily averaged values is somewhat less than observed. Mean biases in the SLAM output for temperature and specific humidity are −0.5°C and −0.5 g kg^{−1}, respectively.

The SLAM’s representation of turbulent heat fluxes has some potentially important departures from the eddy covariance observations. The mean biases in the latent and sensible heat fluxes output from the SLAM are −27 and +26 W m^{−2}, respectively. In addition, the SLAM explains roughly half of the variance in the daily averaged values of both fluxes. These errors could be due to the model’s relative simplicity: the SLAM does not account for wind variability and advection that contribute to turbulent energy flux variability, or for vegetation dynamics that impact the partitioning between latent and sensible heat fluxes. Despite neglecting these complexities, the standard deviation in daily averaged latent heat flux values output by the SLAM is 76% of that observed, indicating that this simple model is able to generate somewhat realistic variability in latent heat flux.

The discrepancies between the observed turbulent fluxes and those output by the SLAM could also result from errors in the eddy covariance observations that have known problems with energy conservation (Franssen et al. 2010). The SLAM’s departures from the observed turbulent fluxes are not always reflected in the temperature errors. For example, in early August, the SLAM output has larger sensible heat fluxes and smaller latent heat fluxes than the observations. Both of these flux biases should generate a warm bias in atmospheric temperature, but the SLAM temperature is lower than observed during this period. Further, while the observations suggest that the SLAM has a low bias in latent heat flux, the soil moisture observations indicate that evaporation from the surface layer is well represented by the SLAM.

The surface saturation *X* simulated by SLAM agrees with the observations, even though the actual values of surface volumetric soil water *m* differ between the SLAM and the observations due to our choice *m*_{sat} and *m*_{w} parameters. However, changing these parameters produces almost identical time series for all five quantities shown in Fig. 5. Because the surface saturation expressed in Eq. (4), rather than the value of *m*, is used to regulate evaporation [see Eq. (2)], the importance of SLAM’s soil moisture values lies in their variability more than their mean; this is a general feature of land surface models (Koster et al. 2009).

## 3. Synthetic forcing

The skill of the SLAM to reproduce the SGP observations for one summer motivates us to use it to understand the processes that control summertime temperatures and temperature variability more generally. To do this, we need a large ensemble of experiments and thus a large ensemble of model forcing and boundary condition time series with high sampling frequency. Such time series are not available for a hydrologically diverse set of regions, or on time scales long enough to study temperature variability. Further, using observations to drive the SLAM makes separating atmospheric forcings from land–atmosphere feedbacks nearly impossible. We have therefore developed synthetic forcing and boundary condition time series that can substitute for observational data and allow us to investigate the relationship between soil moisture, evaporation, and summertime temperatures. These time series and a description of the forcing ensemble used to drive the SLAM are presented in this section, but readers interested primarily in the model results may skip to section 4.

### a. Forcing description

We begin with time series for $R$, *T*_{Top}, and *q*_{Top}. Power spectra of ERA-Interim reanalysis output demonstrate that net surface insolation, 850-hPa temperature, and 850-hPa specific humidity are a combination of red noise and a diurnal cycle during the summer (Dee et al. 2011). With these power spectra in mind, we write time series for *T*_{Top}, *q*_{Top}, and $R$ as

In Eqs. (5)–(7), the Ψ_{x}(*r*_{x}, *t*) terms are red noise time series controlled by a 6-h lag autocorrelation coefficient *r*_{x}. Each red noise component has a multiplicative constant *β*_{x} that controls the amount of red noise variability in each time series. In Eqs. (5) and (6), $TTop\xaf$ and $qTop\xaf$ are the mean temperature and specific humidity at the upper boundary. Although there are diurnal cycles in both temperature and specific humidity at 1000 m, we assume that those cycles are a response to surface processes and do not include them in our external forcing. In Eq. (7), *D*(*t*) is an imposed diurnal insolation cycle with maximum *σ*, while *D*_{max}(*t*) is a cloud-free diurnal cycle with a higher maximum value of peak insolation *σ*_{max}. The term $H$ represents the Heaviside function, which ensures that $R(t)$ never dips below zero.

Once we have the radiation time series [Eq. (7)], we can write the cloud fraction as a function of the insolation red noise time series:

The Heaviside function ensures that *c*_{f} never dips below zero, and we require that 0 < $cf<1$. When the red noise component of the insolation forcing is positive, the cloud fraction must be zero, while a negative value of Ψ_{R} indicates a positive cloud fraction. The unitless *α* term governs how sensitive the cloud fraction is to variations in the insolation forcing; higher values lead to higher cloud fractions for the same insolation forcing variations.

The precipitation time series is generated through stochastic processes that are initiated whenever the cloud fraction is greater than zero, linking precipitation to both cloud fraction and net insolation. At any time step when the sky is cloudy, a random number between 0 and 1 is generated and compared to a threshold value set to 0.9; if the random number is greater than the threshold, precipitation is triggered. While somewhat arbitrary, experiments with this synthetic forcing algorithm show that this threshold value is high enough that cloudy periods generally have less than a trace of precipitation but low enough that significant rainfall occurs on monthly time scales. The rain rate at each time step when precipitation occurs is given by a random value according to a lognormal distribution with specified mean $P$ and standard deviation $\sigma P$. This procedure generates rain rate probability density functions similar to those found in Sauvageot (1994).

### b. Synthetic summers

Using the equations from section 4a, we create an ensemble of forcing time series for SLAM experiments to investigate the impact of soil moisture on summertime temperature variability. The ensemble has 50 sets of the five time series $R$, $P$, *T*_{Top}, *q*_{Top}, and *c*_{f} required to drive the SLAM. Simulations are made to start on 1 June and extend 92 days (three months of summer in the midlatitudes).

The parameters chosen in the forcing algorithm were tuned so that climatological mean values and monthly standard deviations for the $R$, $P$, *T*_{Top}, and *q*_{Top} are similar to those from the central United States in the summer months of June–August (see circled region in Fig. 6). The central United States has been identified as a hot spot of land–atmosphere interaction because it is a transition zone between the wet climate of the American East Coast and the dry climate of the American West (Koster et al. 2004). The monthly means and standard deviations of the synthetic forcing were compared to satellite observations of shortwave radiation from the CERES satellite (CERES Science Team 2000), interpolated weather station precipitation data from the Earth Systems Research Laboratory (Matsuura 2001), and ERA-Interim temperature and specific humidity output at the 850-hPa level. Summertime monthly standard deviations for these four quantities from the years 2000–14 are shown in Fig. 6.

Parameters used to create the forcing ensemble, as well as ensemble monthly means and standard deviations, are shown in Table 3. The $R$ and $P$ monthly means and standard deviations match those found in observations over the central United States. The mean value of the upper-level boundary conditions were taken to be the climatology from ERA-Interim, but the standard deviations in both *T*_{Top} and *q*_{Top} used in the model are reduced compared to those from observations (see Figs. 6c and 6d): we interpret much of the variability in *T* and *q* at 850 hPa as a response to land surface processes rather than to external forcing.

To investigate soil moisture’s impact on summertime temperature variability, we created four ensembles of 50 simulations; the simulations in each ensemble were given the same 50 sets of forcings and boundary conditions. Ensembles differ in that they have progressively more saturated soil moisture initial conditions. These different initial conditions were applied to the *surface and root layers*. Each ensemble has 50 members, for a total of 200 three-month summer simulations. Model parameters are identical to those used to evaluate the SLAM’s performance in the evaluation exercise (see Table 2 and section 2b).

## 4. Evaporation and soil moisture

### a. The source of regime behavior

A logical starting point in our search to understand the relationship between evaporation and soil moisture is vapor pressure deficit *V*:

Like relative humidity, *V* is a measure of atmospheric water vapor demand and is a function of temperature and specific humidity.

Figure 7a shows daily composites of specific humidity observations at 60 and 1000 m from the SGP site during the summer of 2014, along with the daily composite of *q*_{3} output from the SLAM simulation driven by SGP forcings and boundary conditions (see section 2b). The phase relationship between observed *q*(60 m) and *q*(1000 m) indicates that turbulent mixing of water vapor through the boundary layer and dry air entrainment during the daytime are larger influences on near-surface specific humidity than evapotranspiration. Similar phasing of above and within boundary layer surface specific humidity has been found in other observational studies (van Heerwaarden et al. 2010). In contrast, the diurnal cycle in SLAM near-surface specific humidity suggests that evapotranspiration influences *q*_{3} more than turbulence. While the mean value of SLAM’s *q*_{3} and SGP’s *q*(60 m) are similar, the inconsistent phasing of the two signals suggests that a more sophisticated representation of turbulence and/or a variable depth boundary layer in the SLAM could make near-surface specific humidity output more realistic. In addition, the observation-based boundary condition *q*_{Top} = *q*(1000 m) has a diurnal cycle in phase with the SLAM output; the prescribed boundary condition may be unduly influencing the simulation of *q*_{3}. This issue could be remedied by prescribing a constant value of *q*_{Top} rather than the observed value.

Importantly, the inconsistency between the SLAM’s representation of near-surface specific humidity and the observed composite does not preclude accurate modeling of evaporation because *V* is largely controlled by temperature and not by the near-surface specific humidity. Figure 7b shows daily composites of *V* derived from SGP observations of 60-m temperature and specific humidity and the SLAM’s *V* derived from *T*_{3} and *q*_{3}. The SLAM simulates a strong diurnal cycle that is very similar to that observed: the diurnal cycle in *V* is largely determined by the diurnal cycle in temperature. That temperature is the main driver of *V* is not surprising because the Clausius–Clapeyron relationship that governs saturation specific humidity is a function of temperature [see Eq. (9)] and there is a large range in the diurnal cycle of *q*_{s}(*T*) relative to *q*. We expect some inconsistency between the SLAM and observed *V* due to the SLAM’s *q*_{3} errors, specifically the underprediction of *V* during the day when the SLAM overpredicts *q*_{3}. However, despite the inconsistencies in specific humidity between the SLAM and the observations, the two *V* signals agree very well; it is encouraging that the primary quantity driving evapotranspiration is well simulated by the SLAM.

The SLAM’s equation for evaporation [see Eq. (2)], where *E* ∝ *XV*, links soil moisture, vapor pressure deficit, and temperature, but we have not demonstrated the influence of each quantity on the others. Figure 7c shows daily composites of *V* taken from all synthetic model runs colored by the daily averaged surface soil moisture *X*; a uniform moisture increment (of *dX* = 0.1) separates each pair of lines. Brown composites indicate days when the surface is nearly dry, while dark green lines indicate a nearly saturated surface. Increasing soil moisture damps the diurnal cycle of *V* because more evaporative cooling on days with more available soil moisture drives down temperature, and therefore *V*.

As soil moisture increases, the tight grouping of the green lines compared to the brown lines in Fig. 7c suggests that *V* becomes less sensitive to increasing soil moisture. What is the source of the vapor pressure deficit’s decreased sensitivity to soil moisture on the wettest days? To address this question, we define evaporative cooling *T*_{E} (°C) as the cumulative cooling of the surface due to evaporation over the course of one day. We assume that an increment *dE* in evaporation rate produces a proportional change in evaporative cooling *dT*_{E}. Therefore, a Taylor series expansion of Eq. (2) yields

where *C* [K kg air (kg H_{2}O)^{−1}] is a constant given by

where *L* [J (kg H_{2}O)^{−1}] is the enthalpy of vaporization, *c*_{s} (J K^{−1} m^{−2}) is an effective soil heat capacity, and *τ* (s) is the time scale over which we integrate the evaporative cooling. The terms on the right-hand side of Eq. (10) represent two separate contributions to evaporative cooling: the first is a land surface control driven by a perturbation in soil moisture *dX* while the second is an atmospheric control driven by a perturbation in vapor pressure deficit *dV*. Our results (Fig. 7a) and those of others (e.g., van Heerwaarden et al. 2010; Byrne and O’Gorman 2016) have indicated that boundary layer specific humidity is relatively insensitive to evaporative cooling, so we will assume that all *V* perturbations are driven by the influence of evaporative cooling and radiative forcing on near-surface temperature and that specific humidity perturbations are negligible. We can then write a perturbation of vapor pressure deficit as

where *γ* = *dq*_{s}/*dT* [kg H_{2}O (kg air)^{−1} K^{−1}] evaluated at some mean temperature $T\xaf$, and *dT*_{F} is a radiatively induced temperature anomaly that is independent of evaporative cooling.

In the simplest case, we assume that *dT*_{F} = 0 and *dV* = −*γ dT*_{E}. Substituting this assumption into Eq. (10), we obtain

or

We can integrate Eq. (14) from *X* = 0 to *X* to obtain the following closed form solution for *V*(*X*):

where *V*(*X* = 0) is the vapor pressure deficit assuming that no evaporation occurs and some combination of net radiation and sensible heat flux balance the surface energy budget. Daily composites of *V*(*X*) computed by Eq. (15) using *V*(*X* = 0) from the composite shown in Fig. 7c are shown in Fig. 7d and agree remarkably well with the daily *V* composites output from the SLAM experiments (cf. Fig. 7c with Fig. 7d).

The solid black line in Fig. 8 shows the evaporative cooling anomaly *dT*_{E} generated over the course of one day by a soil moisture perturbation *dX* = 0.1 as a function of *X* using Eq. (13) and assuming summertime mean values of *γ* and *V* from the SGP site during the summer of 2014. The damping structure is clearly visible: the same moisture perturbation generates a nearly 6°C anomaly over a completely dry soil, compared to a roughly 1°C anomaly over a completely saturated land surface. In the real world, *γ*, *V*, and *X* change in time: the black dots in Fig. 8 show daily values of *dT*_{E} predicted from daily values of *γ*, *V*, and *X* from the SGP observations. Incorporating daily *γ*, *V*, and *X* variability into Eq. (13) adds very little additional information to this calculation, suggesting that the dominant physical relationship associated with this nonlinear damping structure is the feedback between evaporative cooling and vapor pressure deficit.

This feedback that couples evaporative cooling to soil moisture *is overwhelmingly due to the Clausius–Clapeyron relationship*. The Clausius–Clapeyron relationship gives vapor pressure deficit a strong temperature dependence that connects evaporative cooling, soil moisture, and atmospheric temperature through the *γ* factor. The right-hand side of Eq. (13) illustrates a tug-of-war between soil moisture anomalies that *increase* the land surface’s capacity for evaporative cooling, and *decrease* the atmosphere’s demand for water vapor. When the soil is dry, low mean evaporative cooling implies large a vapor pressure deficit and the only limitation on an evaporative cooling anomaly is the soil moisture perturbation. In the limit ⟨*X*⟩ → 0, this translates to “free” evaporation where *CV dX* = *dT*_{E}. As soil moisture increases, evaporative cooling anomalies decrease as (1 + *CγX*)^{−1}. Although *γ* itself has a temperature dependence, it is large enough at average summertime temperatures that this inverse relationship between *dT*_{E} and *X* asymptotes as *X* → 1. Hence, at high mean soil moisture, evaporative cooling anomalies are insensitive to soil moisture perturbations. This insensitivity is indicative of a strong negative feedback between evaporative cooling and vapor pressure deficit that is most active at high soil moisture.

Hence, instead of the two soil moisture regimes assumed by Budyko, Fig. 8 and Eq. (13) indicate that there is a continuous transition between two limits: (i) a high evaporative cooling sensitivity to soil moisture brought on by high vapor pressure deficit when the soil is dry, and (ii) a low evaporative cooling sensitivity to soil moisture brought on by low vapor pressure deficit when the soil is wet.

### b. Impact of insolation on evaporative cooling

In the real world, evaporation is not the only source of temperature anomalies. If we combine Eqs. (10) and (12) and include a nonzero radiatively forced temperature anomaly that impacts *dV*, we obtain

The first term on the right-hand side of Eq. (16) is identical to the right-hand side of Eq. (13) and describes the impact of a soil moisture perturbation on evaporative cooling. The second term incorporates radiatively forced perturbations to vapor pressure deficit. The shaded region of Fig. 8 shows how a range of temperature perturbations −1°C ≤ *dT*_{F} ≤ 1°C modulates evaporative cooling anomalies across the soil moisture spectrum. When *X* is low, forced temperature perturbations cannot impact evaporative cooling and soil moisture acts as the main constraint on *dT*_{E}. At high values of *X*, the shaded region becomes more substantial, indicating that available soil moisture gives land surfaces the capacity to translate radiatively forced temperature perturbations into evaporative cooling anomalies. We have already demonstrated that *dT*_{E} is insensitive to *dX* for highly saturated soils, thus, we expect temperature perturbations forced by net insolation to be the primary control on evaporative cooling anomalies when *X* is large.

So far, we have argued that a negative feedback exists between soil moisture, evaporative cooling, and vapor pressure deficit. Because of this feedback’s damping structure (Fig. 8), soil moisture anomalies are the primary control on evaporative cooling at low *X*, while at high *X* evaporative cooling is most sensitive to radiatively forced temperature perturbations. To quantify these two patterns of behavior using the SLAM, we define the cumulative daily evaporative cooling (*T*_{E,Σ}) (in °C) as follows:

The term *T*_{E,Σ} gauges the amount of evaporative cooling in the SLAM over one day. If we assume that soil moisture stays constant over the course of one day, we can write Eq. (17) as

Each of the composites in Fig. 7c shows a diurnal cycle of vapor pressure deficit associated with a particular soil moisture value *V*(*X*, *t*). We can integrate each of these composites according to Eq. (18): each integration gives a point in *X*, *T*_{E,Σ} space, and by performing the integration with each *V* composite we can generate a curve of *T*_{E,Σ} as a function of *X*. This exercise can also yield an analytic solution for *T*_{E,Σ} by using Eq. (15):

The red dashed lines in Fig. 9a show three (*X*, *T*_{E,Σ}) curves computed according to Eq. (18) using the composites in Fig. 7c for three different values of *r*_{s} [and therefore *C*; see Eq. (11)], while the black dashed curves show three of the same curves computed according to Eq. (19) with the same three values for *r*_{s}. The distinctive nonlinearity in both sets of curves arises from application of the governing equations *without definition of a critical soil moisture value separating two patterns of behavior*, and is relatively insensitive to the value of *r*_{s}. This is consistent with the discussion of evaporative cooling anomalies above: Fig. 8 shows that cooling anomalies are damped at high values of mean soil moisture, indicating that there is some upper limit on evaporative cooling that can only be modulated by radiative forcing on days when mean soil moisture is high. To explain the reduction in evaporative cooling that occurs at the largest value of soil moisture shown by the red dashed curves in Fig. 9a, we need to examine the connection between radiatively forced temperature perturbations and evaporative cooling.

We use daily evaporation time series output from each day in the SLAM ensembles and Eq. (17) to compute the *T*_{E,Σ} values shown as blue scatter points in Fig. 9a. To illustrate model behavior, each point in Fig. 9a is color coded by daily average insolation; dark blue points correspond to high insolation, while light blue points indicate cloudy days with low insolation. A clear pattern appears across the soil moisture spectrum: higher insolation allows for greater evaporative cooling, while lower insolation restricts the energy available for evaporation. Figure 9a also shows that at the extremely wet end of the soil moisture spectrum, days with reduced evaporative cooling are associated with low insolation (due to rainfall) that drives down *V*.

Figure 9b shows probability distribution functions of *T*_{E,Σ} taken from days in the ensemble with five different values of *X* indicated in the legend. At low values of mean soil moisture, evaporative cooling is tightly constrained by available soil moisture and radiatively forced temperature perturbations cannot generate much spread around the mean value, leading to a small *σ*(*T*_{E,Σ}). In contrast, high soil moisture amplifies the radiatively forced temperature perturbations and generates a large *σ*(*T*_{E,Σ}). This is consistent with our discussion of forced temperature perturbations that preferentially amplify cooling anomalies on days with high soil moisture.

We have argued that one physical law (the Clausius–Clapeyron relationship) governs the nonlinear relationship between evaporation and soil moisture, first noted by Budyko. We have shown in Eq. (13) that even a linearized version of the Clausius–Clapeyron relationship’s strong temperature dependence gives rise to constraints on evaporation that change across the soil moisture spectrum. On the dry end of the spectrum, soil moisture perturbations strongly amplify evaporative cooling, while on the wet end, evaporative cooling becomes insensitive to soil moisture perturbations and is driven primarily by radiative forcing. To investigate the impacts of these different constraints on seasonal time scales, we turn to monthly averaged model output from the SLAM. Figure 10 shows scatterplots of monthly averaged latent heat flux [LHF = *L*(*E* + *η*_{1} + *η*_{2})] as a function of monthly averaged surface saturation *X* from each of the four synthetic ensemble experiments. Experiment 1, where the land surface was initialized with almost no soil moisture, is shown in dark brown (Fig. 10a), while experiment 4, where the SLAM was initialized with an almost completely saturated soil column, is shown in dark green (Fig. 10d). Figures 10e–h show the same monthly averaged values of LHF as a function of net insolation $R$, also partitioned by experiment. Correlations *r* for each scatterplot, along with those for surface sensible heat flux (SSHF), are shown in Table 4. In contrast to LHF, correlations of soil moisture and net insolation with SSHF are nearly constant across the soil moisture experiments. The consistency in SSHF behavior across the soil moisture spectrum is a feature of global climate models and the ERA-Interim reanalysis (Tétreault-Pinard 2013).

In Figs. 10a–d, the correlation between LHF and *X* switches from positive to negative as soil moisture is increased. From our discussion above, we anticipate positive correlation at low soil moisture values because the amplification of evaporative cooling anomalies is highly sensitive to moisture perturbations when *X* is low. As we move to progressively more saturated initialization experiments, we expect the soil moisture control on LHF to diminish [Eq. (13)]. Figures 10c and 10d show an even more marked shift in behavior across the soil moisture spectrum; namely, the negative correlation between LHF and *X* as the land surface becomes increasingly saturated. To explain this behavior, we next examine the relationship between soil moisture, precipitation, and radiative forcing.

We have shown that at high *X*, variability in radiative forcing becomes the dominant source of evaporative cooling variability. In our forcing ensemble, the correlation between monthly insolation and precipitation is −0.58, implying that soil moisture and insolation are also anticorrelated on monthly time scales. At large *X*, the soil moisture control on evaporation diminishes and we expect *X* and LHF to become negatively correlated because positive soil moisture perturbations are associated with months with negative insolation anomalies: since the soil moisture perturbations cannot influence evaporation, insolation perturbations become the only drivers of the correlation. The increasing radiative control on latent heat flux is evident in the increased correlation between $R$ and *X* across the soil moisture initialization experiments. In the driest experiment (Fig. 10e), the correlation between $R$ and LHF is weak and negative. This weak correlation represents a tug-of-war between radiative forcing and soil moisture perturbations that are generated by precipitation. The negative correlation between monthly insolation and precipitation generates a negative correlation between $R$ and *X* that weakens the correlation between $R$ and LHF.

This shift in correlation on monthly time scales comes about because there is a nonlinear relationship between evaporative cooling and soil moisture on daily time scales (Fig. 9a): it is not a product of two distinct soil moisture regimes, but rather it is a consequence of the feedback between evaporative cooling and vapor pressure deficit that preferentially damps evaporation when the soil is wet. We stress again that we have not prescribed any nonlinear behavior in the model that would force this shift in correlation across the spectrum.

### c. Evaporative cooling and summertime temperature variability

We now investigate summertime temperature variability generated by SLAM in the synthetic forcing experiments. Figure 11 shows the distributions of daily averaged near-surface temperature *T*_{3} and surface soil moisture *X* for each of the four experiments. The *x* axes of both plots show the column soil saturation prescribed at the beginning of each summer simulation; the colors of the box plots are the same as those from Fig. 10. The most obvious changes across the four experiments are the mean cooling and surface saturation increase as the initial column moisture grows. However, impacts on variability are also evident in Fig. 11 and summarized in Table 5.

There is a monotonic *but nonlinear* decrease in the standard deviation in temperature *σ*(*T*_{3}) with increasing initial soil moisture. From a 16% reduction in *σ*(*T*_{3}) between experiments 1 and 2 to a 7% reduction between experiments 3 and 4, increasing the initial soil moisture has diminishing returns on decreases in *σ*(*T*_{3}) in the SLAM. We might expect *σ*(*T*_{3}) to be proportional to *σ*(*X*) because of the connection between evaporative cooling and soil moisture, but the muting of *σ*(*X*) with larger initial soil moisture is much more pronounced than the muting of *σ*(*T*_{3}). Changes in soil moisture variability alone cannot explain the way that temperature variability changes across these experiments.

Note that in Fig. 11a, the minimum daily averaged *T*_{3} remains nearly constant between the four experiments: increasing initial soil moisture does not impact minimum daily averaged temperatures across the experiments. We can explain this behavior in terms of the radiative control on evaporation on extremely wet days seen in Figs. 8–10: days with extremely low values of net insolation will drive down vapor pressure deficit, inhibiting any evaporative cooling regardless of available soil moisture. No amount of excess soil moisture can influence the radiative forcing that drives the minimum temperatures in our experiments.

In contrast, the warmest daily averaged *T*_{3} values realized in each experiment decrease significantly as the model is initialized with more soil moisture. The drop in warmest daily averaged temperatures is largest between experiments 1 and 2; we have demonstrated in Fig. 8 that evaporative cooling anomalies are most sensitive to soil moisture when *X* is low. At high *X* where we expect evaporative cooling anomalies to be less sensitive to soil moisture perturbations, we see that the warmest temperatures become less sensitive to increasing the soil moisture initialization. The decreased sensitivity of evaporative cooling anomalies to soil moisture anomalies manifests in a reduced sensitivity of temperature variability to soil moisture initialization across the four experiments.

## 5. Summary and conclusions

We have developed the Simple Land–Atmosphere Model (SLAM) and evaluated the model’s performance by comparing its output to observations of summertime surface climate variability at the Atmospheric Radiation Measurement site in the SGP. The SLAM was designed to include the processes most relevant to summertime temperature variability while limiting the number of parameters. Although we do not expect that the model’s simplicity compromises the essential results and conclusions that we draw, we suggest the following pathways to future work that could improve the model’s capacity to simulate summertime temperature variability:

Currently the SLAM has a fixed constant boundary layer depth that does not account for changes in effective heat capacity of the atmosphere’s lowest layer. While this could influence our model’s representation of boundary layer specific humidity, we consider it unlikely that this simplification impacts our conclusions because evaporation is overwhelmingly controlled by surface temperature rather than atmospheric moisture.

The SLAM’s treatment of vegetation dynamics is very simple. A more advanced vegetation scheme could improve the representation of transpiration and allow modeling of densely vegetated regions, particularly forests. While we restrict our analysis to evaporation in this study, the scatterplots in Fig. 10 include the effects of transpiration and suggest that our conclusions about the relationship between soil moisture and evaporation apply to transpiration as well.

The runoff and drainage fluxes are currently applicable only to flat or rolling terrain. A more sophisticated runoff scheme would be necessary to accurately model soil moisture in regions where runoff dynamics are critical to soil hydrology.

The column framework adopted by the SLAM does not include the impact of atmospheric dynamics on turbulent heat fluxes. In particular, we do not account for low-level wind variability that could modulate these fluxes.

Despite its simplicity, the SLAM displays skill in reproducing key features of the summertime climate at the SGP site. To understand the relationship between evaporation and soil moisture, we generated a synthetic forcing ensemble and used it to drive the SLAM. We created four experiments of model runs with varying soil moisture initial conditions that share the same forcings and boundary conditions.

Without prescribing a nonlinear parameterization that distinguishes between two apparent soil moisture regimes, the SLAM output features a nonlinear relationship between soil moisture and evaporation that very nearly corresponds to the one proposed by Budyko in his 1961 paper. We have shown that the nonlinearity arises from a feedback between evaporative cooling and atmospheric vapor pressure deficit, the strength of which is governed by the temperature dependence of the Clausius–Clapeyron relationship. A set of simple analytic equations demonstrates that this feedback preferentially damps the influence of soil moisture perturbations on evaporative cooling when mean soil moisture is high. For wet soils, the feedback makes radiative forcing the primary driver of evaporative cooling, while for dry soils evaporative cooling anomalies are highly sensitive to soil moisture perturbations.

The relationship between soil moisture and evaporative cooling is of paramount importance to the distribution of summertime temperatures. In our experiments, summertime temperature variability becomes progressively less sensitive to increasing initial soil moisture, a finding that is consistent with previous studies (e.g., Koster et al. 2006). The explanation relies only on the negative feedback between temperature and evaporation, and not on the existence of a critical value of soil moisture that distinguishes the two apparent regimes. The Clausius–Clapeyron relationship that connects temperature and soil moisture through evaporative cooling is thus a sufficient reason to expect the emergence of apparent soil moisture regimes over land surfaces. While other sources of nonlinearity between evaporation and soil moisture surely exist, the impacts of soil moisture perturbations on temperature variability across climatologically distinct wet and dry regimes that have been identified in observations require only this simple physical explanation.

Our results suggest that large-scale land surface drying would not only increase mean temperatures due to less evaporative cooling; it would also increase temperature variability on all time scales by extending the warm tail of the temperature distribution. Large-scale surface drying is projected in the CMIP5 ensemble (Berg et al. 2017), while relative humidity is projected to decrease over land surfaces (Byrne and O’Gorman 2016). Using simple models to understand the source of these changes, and how they may impact summertime temperature variability, is a vital strategy for both validating climate model projections and gaining insight into land–atmosphere interaction.

## Acknowledgments

The authors are grateful for groups participating in the Atmospheric Radiation Measurement (ARM) program, and those who contributed to making the data publicly available online. Similarly, we thank groups that participated in the compilation of the CMIP5 archive on the LLNL server. Abigail Swann, Marysa Laguë, and the Ecoclimate group at the University of Washington provided invaluable feedback at every stage of this research, which contributed to L.R.V.Z’s master’s thesis. We also thank three anonymous reviewers for their helpful feedback. L.R.V.Z. is funded by the NSF GRFP. D.S.B. is funded by the Tamakai Foundation.

### APPENDIX

#### Model Derivation and Equations

##### a. Enthalpy equations

We begin with a thermodynamic formulation of our simplified land–atmosphere system, then describe the SLAM’s fluxes of energy and moisture in terms of the model’s state variables. To obtain model equations that allow us to integrate our state variables forward in time, we need enthalpy equations for the atmospheric and soil layers. In the atmospheric layers, moist static energy *h*_{a} is the sum of enthalpy contributions from dry air and water vapor:

In Eq. (A1), *c*_{a} (J kg^{−1} K^{−1}) is the heat capacity of moist air (assumed constant), *L* (J kg^{−1}) is the enthalpy of vaporization, and *g* (m s^{−2}) is the gravitational potential. We assume that variations in height are negligible within each layer, giving us a specific enthalpy tendency equation for the atmospheric layers:

No single equation for soil enthalpy has been firmly established, largely due to the complex thermodynamics of porous media and the multiphase nature of any soil system [see Nitao and Bear (1996) for a more rigorous thermodynamic treatment.] To derive a soil enthalpy equation using only the thermodynamic variables we have defined in the SLAM, we consider a system that involves dry soil, liquid water, water vapor, and dry air. The total enthalpy of this system *H* is given by

where *M*_{x} terms are masses and *h*_{x} terms represent the specific enthalpy of each substance given by *h*_{x} = *c*_{p,x}*T* + *gz*. The subscripts *d*, *l*, *υ*, and *o* denote dry soil, liquid water, water vapor, and dry air respectively. We can differentiate Eq. (A3) under the assumptions that the mass of dry soil is constant, that changes in the mass of dry air have a negligible contribution to the enthalpy tendency, and that the height *z* of the system is constant:

The water vapor tendency in Eq. (A4) is driven entirely by evaporation of liquid water. However, not all changes in liquid water mass are due to evaporation of soil water; some are externally forced. We denote this by modifying Eq. (A4):

The *E* and *F* subscripts denote changes in liquid water mass associated with evaporation and external forcing respectively. Implicit in Eq. (A5) is the assumption that $dM\upsilon /dt=\u2212dMl/dt|E$. By factoring the mass of dry soil *M*_{d} out of Eq. (A5), we obtain an expression in terms of liquid water mass fraction *w* [kg H_{2}O (kg dry soil)^{−1}]:

where we have substituted *L* in for the specific enthalpy difference between water vapor and liquid water. In going from Eq. (A5) to (A6), we have also ignored *M*_{υ} and *M*_{o} because both of these terms are small compared to *M*_{d}. The temperature *T* of liquid water coming in and out of the system is assumed to be the same as that of the system itself, an assumption we will return to below. Dividing by the total mass of the system (*M*_{d} + *M*_{l}), we obtain

where *h* is the total soil enthalpy per unit mass. We can expand Eq. (A7) under the assumption that *M*_{l}/*M*_{d} = *w* ≪ 1:

Using the definition of volumetric soil moisture [Eq. (1)], we can write Eq. (A8) in terms of *m* instead of *w* after eliminating terms of *O*(*w*^{2}) and higher:

From Eq. (A8), the heat capacity of the soil system is seen to be linearly dependent on volumetric soil moisture: *c*_{s}(*m*) = *c*_{d} + (*c*_{l} − *c*_{d})(*ρ*_{l}/*ρ*_{d})*m*. Similarly, we write *ρ*_{s}(*m*) = *ρ*_{d} + *ρ*_{l}*m* to define a bulk land surface density that increases linearly with soil moisture.

The three terms on the right-hand side of Eq. (A9) demonstrate the triple role of liquid water in changes to soil enthalpy. The first term illustrates how liquid water increases the heat capacity of the soil. The second term illustrates the effect of evaporation; in a closed soil system that does not interact with its environment, *dh*/*dt* = 0 and Eq. (A9) mandates that evaporation of liquid water within the soil layer must be accompanied by a corresponding reduction in soil temperature. The third term accounts for externally forced changes to soil moisture (i.e., precipitation) that change the system’s enthalpy. The next step in model development is to couple the two enthalpy equations [Eqs. (A2) and (A9)] to moisture tendency equations.

##### b. Moisture budget and enthalpy closure

Since there are two unknowns (*T* and *m* or *T* and *q*) in each of the enthalpy tendency equations [Eqs. (A2) and (A9)], we need two equations in each model layer to fully describe our system. Water must be conserved, so we can write water budgets for each layer in terms of model fluxes. We can use the moisture fluxes shown in Fig. 3b to write a moisture budget for each layer from bottom to top with all terms defined in section 2b (given in kg H_{2}O m^{−2} s^{−1}):

In Eqs. (A10)–(A13), the *A*_{i} constants govern the change in *m* or *q* for a specific mass flux of liquid water or water vapor. For the soil layers *A*_{i} = *d*_{i}*ρ*_{l} (kg H_{2}O m^{−2}) where *d*_{i} is the layer thickness, while for the atmospheric layers *A*_{i} = *d*_{i}*ρ*_{a} (kg air m^{−2}).

By combining the moisture budgets and enthalpy tendency equations, we can close our system and derive the temperature tendency equations for each layer. As noted above, changes in enthalpy in the model layers must be due to some combination of external forcings (in W m^{−2}); the sum of these forcings on each model layer is shown below (all terms described in section 2b):

In Eqs. (A14)–(A17), the *B*_{i} (kg m^{−2}) terms govern the change in *h* for a specific enthalpy flux. For each layer, *B*_{i} = *d*_{i}*ρ*_{i} (kg m^{−2}) where *d*_{i} and *ρ*_{i} represent layer thickness and density.

By summing Eqs. (A14) through (A17), we define the model energy balance equation, which is illustrated in Fig. A1:

In Eq. (A18), *F*_{N} is the sum of all net longwave terms. Enthalpy is introduced into the model through the model forcings $R$ and $P$, while boundary fluxes *H*_{↑}, *Q*_{↑}, *H*_{↓}, and *Q*_{↓} act primarily as enthalpy sinks. We have neglected terms that involve the difference between the two soil layer temperatures, as this difference is small relative to mean *T*_{1} or *T*_{2}.

Transpiration moistens the atmosphere, removes liquid water from the soil, and acts as a net enthalpy source in the model because *L*(*η*_{1} + *η*_{2}) > *c*_{l}(*η*_{1}*T*_{1} + *η*_{2}*T*_{2}). This enthalpy input to the model, along with the energy required to transport liquid water from the surface to the vegetation tops, is supplied by plants. Plants also facilitate the phase change from liquid to vapor *without cooling the soil or atmosphere*.

Using conservation equations for moisture [Eqs. (A10)–(A13)] and enthalpy [Eqs. (A14)–(A17)], we can write down the temperature tendency equations for each layer with all terms (in W m^{−2}):

In Eqs. (A19)–(A22), *C*_{i} (J m^{−2} K^{−1}) acts as the heat capacity of the layer and is given by *C*_{i} = *c*_{p,i}*ρ*_{i}*d*_{i} where *c*_{p,i} is the specific heat of the layer. The moisture [Eqs. (A10)–(A13)] and temperature [Eqs. (A19)–(A22)] tendency equations form the backbone of the SLAM. Next, we write each of the fluxes in these tendency equations in terms of the SLAM’s state variables, *T*, *q*, and *m*.

##### c. Soil moisture flux

Moisture movement through porous media is a complex physical process that involves parameters like hydraulic conductivity and soil moisture diffusivity that vary nonlinearly with *m* (DeVries 1958; Libardi et al. 1982; Rawls et al. 1982). These dynamic parameters certainly regulate land surface moisture but we have chosen to neglect them in our model in favor of a simpler formulation of soil moisture movement. We assume that the two layers are in equilibrium unless the soil moisture passes the field capacity *m*_{sat}, at which point the moisture is transferred to the layer below by infiltration *Q*_{1,2} or drainage *Q*_{↓}:

Here, *dt* is the time step used by the model and $H$ is the Heaviside function. The infiltration and drainage fluxes operate only when this saturation value is exceeded. Since the SLAM is a one-dimensional model, we assume that all excess moisture flows downward out of the root layer and into the cold abyss (see Fig. 3b).

##### d. Soil heat flux

Temperature gradients in the soil column can be large near the surface; transfer of heat down these temperature gradients is usually modeled as a diffusive process. The SLAM follows this framework for the conductive soil heat fluxes:

The *λ* value represents soil thermal conductivity (W m^{−1} K^{−1}). We use layer thicknesses as the relevant depths for calculating temperature gradients in the soil.

##### e. Surface sensible heat flux and evapotranspiration

Energy fluxes from the land surface to the atmosphere are often modeled after current flowing through a circuit. In this analogy, the “voltage drop” is determined by the vertical temperature gradient or atmospheric vapor pressure deficit, while the “resistance” is a property of the land surface and the overlying atmosphere. This approach has been deployed to model land surface energy fluxes since the 1980s (Sellers et al. 1986), continues in the present decade (Best et al. 2011), and has proven useful in understanding processes and sources of error in global climate models (Hirsch et al. 2016).

We have discussed the SLAM’s evaporation formula in section 2b. Sensible heat flux *H*_{2,3} and transpiration from the surface *η*_{2} and root layers *η*_{1} are all determined as a function of vertical temperature gradient or vapor pressure deficit and a resistance parameter:

In complex models, *r*_{s} (s m^{−1}) is governed by stability, friction velocity, and roughness length (Garratt 1992). Similarly, complex models parameterize the stomatal resistance *r*_{st} (s m^{−1}) according to the plant species, ambient CO_{2} concentrations, and photosynthetic rate (Collatz et al. 1991). For simplicity, and to avoid the introduction of too many parameters, we hold values of these two quantities constant in the SLAM. The parameters *f*_{υ} and *f*_{r} influence the partitioning of moisture fluxes in the SLAM between evaporation and transpiration and between the two soil layers. The *f*_{υ} parameter is the fractional vegetation cover of the land surface, while *f*_{r} is the fraction of plant roots that penetrate the 10-cm layer. We will draw values of *f*_{r} from Jackson et al. (1996). For a discussion of the moisture-flux coupling term *X*, see section 2b.

There is an important distinction between evaporation and transpiration based on the level where each process is taking place. The vapor pressure deficit at the surface (driven by *T*_{2}) drives evaporation while the same deficit in the lower boundary layer (driven by *T*_{3}) drives transpiration from both root layers because we assume that vegetation responds to atmospheric temperature while evaporation is driven by the land surface temperature.

##### f. Turbulent heat and moisture fluxes within the atmosphere

Turbulent fluxes between the atmospheric layers are also formulated in the resistance framework. In the SLAM, boundary layer turbulence transports energy down temperature and humidity gradients. Turbulent fluxes are inhibited by a resistance that depends on the stability of the layer in question:

In Eqs. (A30)–(A33), we introduce two new resistance parameters, one for atmospheric heat transfer *r*_{a,T} (s m^{−1}) and another for atmospheric vapor transfer *r*_{a,q} (s m^{−1}). In general, water vapor is well mixed in the boundary layer while temperature has a distinct vertical structure; resistances governing vapor transfer in our model are smaller than those governing heat transfer. Since turbulence is largely dependent on atmospheric stability, we vary atmospheric turbulent resistance according to the sum of buoyancy fluxes:

In Eq. (A34), $ra,T\xaf$ and $ra,q\xaf$ are the maximum nighttime values of turbulent resistance for heat and water vapor, while *σ*_{max} is a constant equal to the daily maximum cloud-free insolation (in W m^{−2}). This formulation ensures that nighttime turbulence is weak, while daytime surface fluxes contribute to vigorous boundary layer mixing. These resistances are calculated at each time step.

##### g. Longwave radiation

Net absorbed longwave radiation in each model layer is a function of the vertical temperature and specific humidity profile. The emissivity of the atmospheric layers *ϵ*_{i} is logarithmic with specific humidity *q*_{i}:

The parameters in Eq. (A35) were chosen so that the sensitivity of downward longwave radiation to temperature and humidity in the SLAM broadly matches the sensitivity found in radiative kernels (see Previdi 2010). Equations for the net longwave absorption in each layer are given below:

An important component of Eqs. (A36)–(A38) is the longwave emissivity of the free tropospheric layer *ϵ*_{↑}, which is governed both by specific humidity and cloud fraction *c*_{f}:

The value of 0.4 was chosen to yield values of *ϵ*_{Top} near unity for typical values of boundary layer specific humidity when the sky is completely cloudy. The emissivity *ϵ*_{Top} is capped at one, which is the value of a perfect blackbody associated with a completely cloudy sky.

##### h. A note on numerical methods

So far, we have discussed computation of fluxes in the SLAM from state variables (*T*, *q*, and *m*). To integrate these equations forward, we use an explicit numerical method where the fluxes are computed from the state variables at time step *i*, then used to find the state variables at time step *i* + 1. A schematic set of equations for this explicit method is shown below where a flux *F* is computed using state variables *T*, *q*, and *m* from the *i*th time step and then used to integrate the state variables forward:

An important value in these equations is the time step *dt*, which must be small for this explicit method to successfully model the short time scale variability that the SLAM is designed to study. This method represents a departure from land surface schemes used in global climate models, which often use an implicit numerical method that can result in loss of energy conservation on subdaily time scales (Shultz et al. 2000). However, as long as a small enough time step is used, the explicit method is accurate.

## REFERENCES

*The Heat Balance of the Earth’s Surface.*U.S. Department of Commerce, 259 pp.

*The Atmospheric Boundary Layer*. Cambridge University Press, 334 pp.

*Climate Change 2013: The Physical Science Basis*, T. F. Stocker et al., Eds., Cambridge University Press, 953–1028.

*Soil and Environmental Analysis: Physical Methods.*Taylor and Francis, 656 pp.

## Footnotes

^{a}

Current affiliation: Department of Earth and Space Science, University of Washington, Seattle, Washington.

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).