Relative humidity gradients as a key constraint on terrestrial water and energy fluxes

Earth's climate and water cycle are highly dependent on terrestrial evapotranspiration and the associated flux of latent heat. Although it has been hypothesized for over 50 years that land dryness becomes embedded in atmospheric conditions through evaporation, underlying physical mechanisms for this land–atmosphere coupling remain elusive. Here, we use a novel physically based evaporation model to demonstrate that near-surface atmospheric relative humidity (RH) fundamentally coevolves with RH at the land surface. The new model expresses the latent heat flux as a combination of thermodynamic processes in the atmospheric surface layer. Our approach is similar to the Penman–Monteith equation but uses only routinely measured abiotic variables, avoiding the need to parameterize surface resistance. We applied our new model to 212 in situ eddy covariance sites around the globe and to the FLUXCOM global-scale evaporation product to partition observed evaporation into diabatic vs. adiabatic thermodynamic processes. Vertical RH gradients were widely observed to be near zero on daily to yearly timescales for local as well as global scales, implying an emergent land–atmosphere equilibrium. This equilibrium allows for accurate evaporation estimates using only the atmospheric state and radiative energy, regardless of land surface conditions and vegetation controls. Our results also demonstrate that the latent heat portion of available energy (i.e., evaporative fraction) at local scales is mainly controlled by the vertical RH gradient. By demonstrating how land surface conditions become encoded in the atmospheric state, this study will improve our fundamental understanding of Earth's climate and the terrestrial water cycle.

1 Introduction

Latent heat flux (LE) associated with plant transpiration and evaporation from soil and intercepted water (i.e., evapotranspiration, ET) links the water cycle with the terrestrial energy budget. More than half of the incoming radiation energy at the land surface is consumed as LE, making ET the second largest flux in the terrestrial water balance after precipitation (Oki and Kanae, 2006). Also, LE is a controlling factor for near-surface climatic conditions such as temperature and humidity (Ma et al., 2018; Byrne and O'Gorman, 2016). While most research has been devoted to developing and improving rate-limiting parameters constraining LE (e.g., García et al., 2013; Martens et al., 2017), exploring the governing physics of LE has received less attention following earlier pioneering work (Schmidt, 1915; Penman, 1948; Bouchet, 1963; Monteith, 1965; Priestley and Taylor, 1972). Nevertheless, improvement of the theoretical understanding of LE still remains an essential cornerstone to correctly simulate and predict climate and hydrological cycles (Emanuel, 2020).

Climatic conditions over the land surface have been getting not only warmer but also drier in recent decades (i.e., decrease in relative humidity) (Sherwood and Fu, 2014; Willett et al., 2014; Byrne and O'Gorman, 2018), but land–atmosphere feedback processes shaping the near-surface atmospheric state are still not well understood. In the early 1960s, Bouchet (1963) hypothesized that land surface dryness is coupled to the atmospheric state through LE, with the Bouchet hypothesis now widely accepted (Ramírez et al., 2005; Fisher et al., 2008; Mallick et al., 2014). However, the underlying physical mechanisms for this land–atmosphere coupling still remain elusive (McNaughton and Spriggs, 1989). Recently, McColl et al. (2019) introduced a novel theoretical perspective on land–atmosphere coupling which is referred to as “surface flux equilibrium (SFE)”. They hypothesized that relative humidity (RH) reaches a steady-state value in an idealized atmospheric boundary layer at daily to monthly timescales. Under steady RH conditions (i.e., the SFE state), LE can be determined using only the atmospheric state and radiative energy. Although this method performed well compared to actual LE observations for inland continental sites (McColl and Rigden, 2020; Chen et al., 2021), a further investigation is needed to understand how dynamics of turbulent heat fluxes in the atmospheric surface layer at sub-daily timescales evolve to the SFE state.

A traditional way to express the atmospheric surface layer processes is to partition LE into diabatic and adiabatic processes using the Penman–Monteith (PM) equation (Monteith, 1965), as proposed by Monteith (1981). The PM equation combines the energy balance equation with mass-transfer theory for water vapour and sensible heat, resulting in diabatic (radiative energy-related) and adiabatic (vapour pressure deficit-related) processes for a parcel of air in contact with a saturated surface (Monteith, 1981). ​​​​​​​

$$\begin{array}{}\text{(1)}& \mathit{LE}=\underset{\text{Diabatic process}}{\underbrace{{\displaystyle \frac{S}{S+\mathit{\gamma }\left(\frac{r_{\mathrm{a}}+r_{\mathrm{s}}}{r_{\mathrm{a}}}\right)}}\cdot Q}}+\underset{\text{Adiabatic process}}{\underbrace{{\displaystyle \frac{\mathit{\rho }{c}_p}{S+\mathit{\gamma }\left(\frac{r_{\mathrm{a}}+r_{\mathrm{s}}}{r_{\mathrm{a}}}\right)}}\cdot {\displaystyle \frac{e^{\ast }\left(T_{\mathrm{a}}\right)-e_{\mathrm{a}}}{r_{\mathrm{a}}}}}},\end{array}$$

where S is the linearized slope of saturation vapour pressure versus temperature (hPa K⁻¹), γ is the psychrometric constant (hPa K⁻¹), ρ is the air density (kg m⁻³), c_p is the specific heat capacity of air at constant pressure (MJ kg⁻¹ K⁻¹), and Q is available radiative energy (i.e., the difference between net radiation (R_n) and soil heat flux (G) expressed in units of W m⁻²). e^∗(T_a) is the saturation vapour pressure (hPa) corresponding to the air temperature (T_a) measured at a reference height (typically 2 m or eddy flux measurement height), and e_a is vapour pressure (hPa) at the reference height. $e^{\ast }\left(T_{\mathrm{a}}\right)-e_{\mathrm{a}}$ is known as atmospheric vapour pressure deficit (VPD, expressed in units of hPa). r_a is aerodynamic resistance to heat and water vapour transfer (s m⁻¹), and r_s is surface resistance to water vapour transfer (s m⁻¹) representing drying soil and/or plant stomatal closure.

In principle, high Q and VPD at the reference height increase the diabatic and the adiabatic terms respectively in the PM equation, and as such, Q and VPD are the two primary drivers of evaporation (Monteith and Unsworth, 2013). Yet, this “high VPD leads to high LE” interpretation cannot be generalized because r_s increases with VPD due to stomatal closure by vegetation under high-VPD conditions (Tan et al., 1978; Novick et al., 2016; Massmann et al., 2019). While the PM equation is useful to explore biological control of LE through r_s (Jarvis and McNaughton, 1986; Peng et al., 2019), physical mechanisms corresponding to each term in Eq. (1) are less intuitive due to the sensitivity of r_s to VPD. As a result, how the atmospheric state affects evaporation and vice versa remains ambiguous in the PM equation.

Is there a way to mathematically express the physical mechanisms of LE without requiring r_s? In this paper, we present a pair of equations expressing actual LE as a combination of diabatic and adiabatic processes without requiring r_s. Similar to the PM equation, our new equations are derived by combining the energy balance equation with the flux gradient equations, but crucially ours do not include r_s. The novel equations are applied empirically to eddy-covariance observation sites and a global LE dataset to explore land–atmosphere coupling processes at various spatiotemporal scales. To do this, we decomposed observed LE into adiabatic and diabatic components and discuss how these patterns can help to understand land–atmosphere interactions and potential responses under future climatic conditions.

2 Theory

2.1 A pair of evaporation equations for an unsaturated surface

In this section, we derive a pair of evaporation equations for an unsaturated surface. In this derivation, we assume a horizontally homogenous landscape where the sources of water vapour and heat are identical. Under this idealized condition, aerodynamic resistances (r_a) to heat and water vapour transfers are identical. Here, r_a is a parameterization of turbulent mixing due to mechanical turbulence and buoyancy driven by surface heating.

We first express LE using a flux gradient equation as $\mathit{LE}=\frac{\mathit{\rho }{c}_p}{\mathit{\gamma }}\frac{e_{\mathrm{s}}-e_{\mathrm{a}}}{r_{\mathrm{a}}}$, where e_s is the surface vapour pressure. Here, the subscript s indicates the land surface which is defined as an idealized plane specified as the sum of displacement height and roughness length for heat (Knauer et al., 2018a; Novick and Katul, 2020). If the land surface is saturated, e_s becomes equivalent to the saturation vapour pressure (i.e., $e_{\mathrm{s}}=e^{\ast }\left(T_{\mathrm{s}}\right)$). For an unsaturated land surface, however, relative humidity should be introduced as $e_{\mathrm{s}}={\mathrm{RH}}_{\mathrm{s}}{e}^{\ast }\left(T_{\mathrm{s}}\right)$, where RH_s is surface relative humidity, i.e., the ratio of e_s to e^∗(T_s). For a vegetated surface, RH_s as defined in this study represents relative humidity of the foliage surface and is conceptually equivalent to surface water availability in Li and Wang (2019). For a bare soil land surface, RH_s represents soil surface relative humidity, which can be found using the “alpha” method that is parameterized using soil moisture content or soil water potential (Lee and Pielke, 1992; Wu et al., 2000; Cuxart and Boone, 2020). Using RH_s, LE can be written as $\mathit{LE}=\frac{\mathit{\rho }{c}_p}{\mathit{\gamma }}\frac{{\mathrm{RH}}_{\mathrm{s}}{e}^{\ast }\left(T_{\mathrm{s}}\right)-e_{\mathrm{a}}}{r_{\mathrm{a}}}$ for an unsaturated surface condition.

In order to decompose LE into two individual fluxes related to temperature and relative humidity gradients, we add $-{\mathrm{RH}}_{\mathrm{s}}{e}^{\ast }\left(T_{\mathrm{a}}\right)+{\mathrm{RH}}_{\mathrm{s}}{e}^{\ast }\left(T_{\mathrm{a}}\right)$ (or $-{\mathrm{RH}}_{\mathrm{a}}{e}^{\ast }\left(T_{\mathrm{s}}\right)+{\mathrm{RH}}_{\mathrm{a}}{e}^{\ast }\left(T_{\mathrm{s}}\right)$) to the numerator of the flux gradient equation and rewrite LE as follows.

$$\begin{array}{}\text{(2)}& {\displaystyle }\begin{array}{rl}\mathit{LE}& ={\displaystyle \frac{\mathit{\rho }{c}_p}{\mathit{\gamma }}}{\mathrm{RH}}_{\mathrm{s}}{\displaystyle \frac{e^{\ast }\left(T_{\mathrm{s}}\right)-e^{\ast }\left(T_{\mathrm{a}}\right)}{r_{\mathrm{a}}}}\\ & +{\displaystyle \frac{\mathit{\rho }{c}_p}{\mathit{\gamma }}}{e}^{\ast }\left(T_{\mathrm{a}}\right){\displaystyle \frac{{\mathrm{RH}}_{\mathrm{s}}-{\mathrm{RH}}_{\mathrm{a}}}{r_{\mathrm{a}}}}\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}\mathit{LE}& ={\displaystyle \frac{\mathit{\rho }{c}_p}{\mathit{\gamma }}}{\mathrm{RH}}_{\mathrm{a}}{\displaystyle \frac{e^{\ast }\left(T_{\mathrm{s}}\right)-e^{\ast }\left(T_{\mathrm{a}}\right)}{r_{\mathrm{a}}}}\\ & +{\displaystyle \frac{\mathit{\rho }{c}_p}{\mathit{\gamma }}}{e}^{\ast }\left(T_{\mathrm{s}}\right){\displaystyle \frac{{\mathrm{RH}}_{\mathrm{s}}-{\mathrm{RH}}_{\mathrm{a}}}{r_{\mathrm{a}}}}\end{array}\end{array}$$

We then approximate $e^{\ast }\left(T_{\mathrm{s}}\right)-e^{\ast }\left(T_{\mathrm{a}}\right)=S(T_{\mathrm{s}}-T_{\mathrm{a}})$ using the saturation vapour pressure slope at the air temperature (S), and we introduce a flux gradient equation for sensible heat flux (i.e., $H=\mathit{\rho }{c}_p\frac{T_{\mathrm{s}}-T_{\mathrm{a}}}{r_{\mathrm{a}}}$) into Eqs. (2) and (3). Then, the energy balance equation is combined to substitute H with Q−LE. As results, LE can be expressed as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{rl}\mathit{LE}& =\underset{\text{Diabatic process}:{\mathit{LE}}_Q}{\underbrace{{\displaystyle \frac{{\mathrm{RH}}_{\mathrm{s}}S}{{\mathrm{RH}}_{\mathrm{s}}S+\mathit{\gamma }}}\cdot Q}}+\underset{\text{Adiabatic process}:{\mathit{LE}}_G}{\underbrace{{\displaystyle \frac{\mathit{\rho }{c}_p{e}^{\ast }\left(T_{\mathrm{a}}\right)}{{\mathrm{RH}}_{\mathrm{s}}S+\mathit{\gamma }}}\cdot {\displaystyle \frac{{\mathrm{RH}}_{\mathrm{s}}-{\mathrm{RH}}_{\mathrm{a}}}{r_{\mathrm{a}}}}}},\\ & ={\mathit{LE}}_Q+{\mathit{LE}}_G\end{array}\text{(5)}& {\displaystyle }\begin{array}{rl}\mathit{LE}& =\underset{\text{Diabatic process}:{\mathit{LE}}_Q^{\prime }}{\underbrace{{\displaystyle \frac{{\mathrm{RH}}_{\mathrm{a}}S}{{\mathrm{RH}}_{\mathrm{a}}S+\mathit{\gamma }}}\cdot Q}}+\underset{\text{Adiabatic process}:{\mathit{LE}}_G^{\prime }}{\underbrace{{\displaystyle \frac{\mathit{\rho }{c}_p{e}^{\ast }\left(T_{\mathrm{s}}\right)}{{\mathrm{RH}}_{\mathrm{a}}S+\mathit{\gamma }}}\cdot {\displaystyle \frac{{\mathrm{RH}}_{\mathrm{s}}-{\mathrm{RH}}_{\mathrm{a}}}{r_{\mathrm{a}}}}}}\\ & ={\mathit{LE}}_Q^{\prime }+L{E}_G^{\prime },\end{array}\end{array}$$

where LE_Q (and LE${}_Q^{\prime }$) is a diabatic component, and LE_G (and LE${}_G^{\prime }$) is an adiabatic component of latent heat flux. While the diabatic component is mainly determined by available energy (Q), the adiabatic component is driven by turbulent mixing and vertical gradient of RH. Monteith (1981) originally suggested an equation equivalent to Eq. (4) for the case when the surface does not reach saturation. To our knowledge, Eq. (5) is derived for the first time here. Equations (4) and (5) include RH_s to compensate for eliminating r_s from the original PM equation.

Since the adiabatic processes in Eqs. (4) and (5) are controlled by the vertical difference of RH, we refer to Eqs. (4) and (5) as the proposed PM_RH model (Penman–Monteith equation expressed using RH) to distinguish it from the original PM model. The two Eqs. (4) and (5) are complementary to each other in that they represent distinct thermodynamic paths, each of which will be discussed in the next section. Arguably, applying PM_RH can provide new insights into the fundamental mechanisms of LE, particularly when it is decomposed into its diabatic component (LE_Q or LE${}_Q^{\prime }$) and its adiabatic component (LE_G or LE${}_G^{\prime }$). In the following sections, we will discuss theoretical meanings of Eqs. (4) and (5) in depth.

2.2 Generalized Penman equation

Before discussing PM_RH in depth, we revisit the Penman equation (Penman, 1948) to help with the physical reasoning behind our proposed framework. The widely recognized form of the Penman equation, which was developed as an LE model for a saturated surface, is as follows:

$$\begin{array}{}\text{(6)}& \mathit{LE}=\underset{\text{Diabatic process}}{\underbrace{{\displaystyle \frac{S}{S+\mathit{\gamma }}}\cdot Q}}+\underset{\text{Adiabatic process}}{\underbrace{{\displaystyle \frac{\mathit{\rho }{c}_p\left[e^{\ast }\left(T_{\mathrm{a}}\right)e_{\mathrm{a}}\right]}{\left[S+\mathit{\gamma }\right]r_{\mathrm{a}}}}}}.\end{array}$$

We rearrange this formulation to derive Eq. (7) by factoring out e^∗(T_a) and introducing ${\mathrm{RH}}_{\mathrm{a}}=\frac{e_{\mathrm{a}}}{e^{\ast }\left(T_{\mathrm{a}}\right)}$ into the second term.

$$\begin{array}{}\text{(7)}& \mathit{LE}=\underset{\text{Diabatic process}}{\underbrace{{\displaystyle \frac{S}{S+\mathit{\gamma }}}\cdot Q}}+\underset{\text{Adiabatic process}}{\underbrace{{\displaystyle \frac{\mathit{\rho }{c}_p{e}^{\ast }\left(T_{\mathrm{a}}\right)}{S+\mathit{\gamma }}}\cdot {\displaystyle \frac{\mathrm{1}-{\mathrm{RH}}_{\mathrm{a}}}{r_{\mathrm{a}}}}}}\end{array}$$

Equations (6) and (7) are mathematically equivalent, but their interpretations are quite different. In Eq. (6), the adiabatic process is controlled by VPD at the reference height. However, in Eq. (7), the adiabatic process acts over the vertical RH gradient, i.e., the difference in RH from the surface to the reference height (RH_a). Since the Penman equation is a model for saturated surfaces, 1−RH_a in Eq. (7) indicates the difference in RH over the vertical distance between the ground surface and the reference height. Arguably, Eq. (7) is more thermodynamically sound compared to Eq. (6) since RH is an ideal-gas approximation to the water activity (Lovell-Smith et al., 2015) which represents the chemical potential of water (μ_w) (Monteith and Unsworth, 2013; Kleidon and Schymanski, 2008). When the vertical gradient of RH dissipates owing to well-developed turbulence, the land surface and the atmosphere are in thermodynamic equilibrium (Kleidon et al., 2009). Therefore, taking Eq. (7) instead of Eq. (6) allows us to view the adiabatic process of the Penman model as an equilibration process driving land–atmosphere equilibrium by bringing the surface μ_w to that of the atmosphere.

As with our interpretation of the Penman model, we can view Eqs. (4) and (5) as a generalized form of the Penman model. Here, the LE_G (or LE${}_G^{\prime }$) term is an equilibration process between the land and the atmosphere when the land surface is not saturated. It is worth noting that LE_G can be negative when RH_s is less than RH_a. Thus, the LE_G term operated by turbulent mixing acts to reduce the vertical RH gradient. This physical interpretation is consistent with recent findings that the variance of the RH gradient tends to be minimized over the course of the day, implying that the difference between RH_s and RH_a is reduced (Salvucci and Gentine, 2013; Rigden and Salvucci, 2015). The diabatic LE_Q (or LE${}_Q^{\prime }$) term can be understood as equilibrium LE for an unsaturated surface, which we discuss later in Sect. 2.4.

2.3 Thermodynamic paths

How can we interpret the two formulas of PM_RH in Eqs. (4) and (5)? To explain the two forms, the psychrometric relationship is applied to a parcel of air near an unsaturated land surface that is under constant pressure and steadily receiving radiation energy. The psychrometric diagram in Fig. 1 describes the magnitude of turbulent flux (where the length of the arrow corresponds to the magnitude) from the view point of a parcel of air located at a reference height (an approach based on work by Monteith, 1981). Since the parcel of air receives heat and water vapour from the land surface, the final state is represented by the surface condition, while the initial state is represented by the atmospheric conditions at the reference height. Therefore, the initial thermodynamic state of the air parcel can be represented by its temperature and water vapour pressure such as point A in Fig. 1. The initial state is changed by two processes as follows: (1) equilibrating between the land surface (RH_s) and the air parcel (RH_a) and (2) increasing enthalpy forced by the incoming energy. It should be noted that the changing process (i.e., thermodynamic path) from the initial to the final states in this discussion should be understood as the magnitude of the turbulent heat fluxes.

Figure 1 Schematic conceptualization of the PM_RH model and psychrometric relationship of PM_RH. The example psychrometric chart is modified from Marsh (2018). Path 1 represents Eq. (4) divided by $\frac{\mathit{\rho }{c}_p}{\mathit{\gamma }{r}_{\mathrm{a}}}$ while Path 2 represents Eq. (5) divided by $\frac{\mathit{\rho }{c}_p}{\mathit{\gamma }{r}_{\mathrm{a}}}$. Here, the enthalpy change of the air parcel is defined as $\frac{Q\cdot r_{\mathrm{a}}}{\mathit{\rho }}$ (kJ kg⁻¹). It should be noted that the difference between the initial and the final states represents the magnitude of the turbulent heat fluxes instead of changes in atmospheric state.

In the RH equilibrating process, the air parcel is adiabatically cooled (or heated when RH_s < RH_a) due to turbulent mixing, while the enthalpy of the parcel is not changed. Therefore, the increase (decrease) in latent heat content in the parcel is exactly balanced by a decrease (increase) in sensible heat (A→B in Fig. 1: trajectory along constant enthalpy line). This process is equivalent to the LE_G term in Eq. (4). Now, the air parcel is in thermodynamic equilibrium with the land surface (point B in Fig. 1). Then, the air parcel receives energy while the equilibrium is sustained (i.e., RH_s is steady), which increases both the temperature and absolute water vapour content of the air parcel (B→C in Fig. 1). This process can be expressed as LE_Q of Eq. (4). Consequently, the thermodynamic state of the air parcel approaches point C in Fig. 1.

However, we should recognize that temperature and vapour pressure are “state” variables, meaning that they do not depend on the thermodynamic path by which the system arrived at its final state (Iribarne and Godson, 1981). In the above example, we conceptually followed the adiabatic process first and then the diabatic process (Path 1 in Fig. 1), but one can imagine the opposite order. If we choose Path 2 in Fig. 1, the diabatic process comes first, and thus RH_a instead of RH_s is preserved while enthalpy increases (i.e., LE${}_Q^{\prime }$), and the adiabatic process is followed at temperature of T_S (i.e., LE${}_G^{\prime }$). Path 2 is described by Eq. (5).

Therefore, one can interpret the two forms of PM_RH in Eqs. (4) and (5) as two thermodynamic paths where the diabatic and adiabatic processes occur simultaneously. It should be noted that the diabatic and adiabatic processes in PM_RH are “path” functions and thus they vary by path. For instance, LE_Q is slightly higher than LE${}_Q^{\prime }$ when RH_s > RH_a. Also, the absolute magnitude of LE${}_G^{\prime }$ is always bigger than that of LE_G when Q>0 (i.e., vector $B^{\prime }\to C$ is longer than vector A→B in Fig. 1).

2.4 Equilibrium LE for an unsaturated surface

Another distinct characteristic of the PM_RH model is the way it defines equilibrium at the land–atmosphere interface. Unlike many previous studies which focused on the steady state of VPD (McNaughton and Jarvis, 1983; Priestley and Taylor, 1972; Raupach, 2001), land–atmosphere equilibrium is achieved in the PM_RH model when the vertical RH gradient (i.e., the μ_w gradient) dissipates. That is, if RH_s≈RH_a, then it follows that LE_G (or LE${}_G^{\prime }$) is zero and thus LE becomes

$$\begin{array}{}\text{(8)}& \mathit{LE}\approx {\displaystyle \frac{{\mathrm{RH}}_{\mathrm{a}}S}{{\mathrm{RH}}_{\mathrm{a}}S+\mathit{\gamma }}}Q.\end{array}$$

We note that Eq. (8) is identical to the SFE theory recently introduced by McColl et al. (2019). They hypothesized that in many continental regions, the near-surface atmosphere is in a state of equilibrium, where RH is steady with time in an idealized atmospheric boundary layer at longer than daily timescales. Equation (8) successfully predicted observed LE at daily and monthly timescales for inland regions (McColl and Rigden, 2020; Chen et al., 2021), which implies the vertical RH gradient tends to evolve toward zero at longer timescales than the sub-daily scale. This is logical in that LE_G itself diminishes the vertical RH gradient over the course of a day.

From a different standpoint, if an observed LE is bigger or smaller than Eq. (8) at a longer timescale such as monthly, it may indicate that the land surface conditions are not completely embedded in the near-surface atmospheric state due to highly wet or dry land conditions. Therefore, LE_G (or LE${}_G^{\prime }$) value and sign at monthly timescales could be a useful indicator reflecting land surface dryness relative to the atmosphere.

When both land surface and atmosphere are saturated (i.e., ${\mathrm{RH}}_{\mathrm{s}}\approx {\mathrm{RH}}_{\mathrm{a}}\approx \mathrm{1}$), Eq. (8) becomes classical equilibrium LE (i.e., $\mathit{LE}\approx \frac{S}{S+\mathit{\gamma }}Q$). This is consistent with one of the classical definitions of equilibrium LE that defines equilibrium LE as evaporation from a saturated surface into saturated air (Schmidt, 1915; Eichinger et al., 1996; Raupach, 2001; McColl, 2020). Therefore, we can regard Eq. (8) as a generalized equilibrium LE for an unsaturated surface.

3 Materials and methods

In the following sections, we present a novel physical decomposition of LE from PM_RH into LE_Q and LE_G components to aid in understanding the governing physics of LE. Also, the proportion of net available energy consumed in evapotranspiration, known as the evaporative fraction (EF) $\left(\mathrm{EF}=\frac{\mathit{LE}}{Q}\right)$ is decomposed into $\frac{{\mathit{LE}}_Q}{Q}$ and $\frac{{\mathit{LE}}_G}{Q}$. We conducted a detailed diagnostic analysis of the PM_RH model using the multi-year record of an eddy covariance (EC) flux observation site located in a wet–dry tropical climate. We also applied the PM_RH model to the 212 EC sites represented in the FLUXNET2015 dataset (Pastorello et al., 2020) and to the FLUXCOM global LE product (Jung et al., 2019). We describe the local and global datasets and analysis methods here before presenting the results.

3.1 In situ EC flux observation

In situ half-hourly EC observations used in this study were made from 2015 to 2018 on a ratoon sugarcane farm in the province of Guanacaste, Costa Rica (10^∘25^′07.60^′′ N, 85^∘28^′22.22^′′ W). The site has a wet–dry tropical climate with a dry season from December to March and a median monthly air temperature ranging from 27 to 30 ^∘C. The study site experienced a significant drought in 2015 as the lowest precipitation rate in Fig. 2b (Hund et al., 2018; Morillas et al., 2019). The site was irrigated occasionally during dry seasons via furrow irrigation events, except for 2016 when there was no irrigation due to crop replanting. Due to the ratooning practice (i.e., sugarcane cutting each year followed by resprouting without replanting, detailed explanation in the Supplement), the sugarcane growing seasons varied by year, which provided an opportunity to explore distinct and varied combinations of land surface vs. atmospheric aridity conditions.

Figure 2 Time series for the sugarcane EC tower site in Costa Rica. Panel (a) is daily heat fluxes and panel (b) is monthly precipitation. Panel (c) is mean annual LE and its two components, and (d) is time series of LE, LE_Q, and LE with a background colour of RH_a. Dashed lines with “h” in panel (d) indicate sugarcane harvest. Panels (e) and (f) are wavelet coherence of LE with LE_Q and LE with LE_G. Red and blue colours indicate high and low correlation, respectively. Arrows (pointing right: in-phase; left: antiphase) only appear when the coherence is significant (p < 0.01).

The measured LE and sensible heat flux (H) were quality controlled following Morillas et al. (2019) (details in the Supplement). For the study period, the surface energy balance closure (i.e., $\frac{\mathit{LE}+H}{R_{\mathrm{n}}-G}$) of 30 min data was 86 %, which is typical of high-quality eddy-covariance datasets (Wilson et al., 2002). When canopy height was less than 1 m, the surface energy balance was almost closed (97 %), whereas the closure was 83 % when canopy height was higher than 1 m. It is expected that unmeasured canopy and soil heat storages in this site are significant because the sugarcane canopy grew up to 3.6 m tall with a dense canopy. For instance, Meyers and Hollinger (2004) showed that storage term comprised 14 % of net radiation for a maize field with a 3 m canopy height and 8 % of net radiation for a soybean field with a 0.9 m canopy height, implying larger heat storage capacities for taller crop canopies. Also, since our study site is located within a homogenous landscape (Fig. S1 in the Supplement), horizontal and vertical advective flux divergence and the influence of secondary circulations on the energy balance closure may be marginal (Mauder et al., 2020; Leuning et al., 2012). Therefore, considering the homogenous landscape of the study site as well as a possible significant role of unmeasured canopy and soil heat storages, we did not force the energy closure. Consequently, we defined Q as the sum of LE and H instead of R_n−G. In doing so, we in effect attribute the cause of the surface energy imbalance to unmeasured heat storage terms following Moon et al. (2020).

In order to decompose LE into LE_Q and LE_G, we first estimated half-hourly aerodynamic resistance (r_a) by considering aerodynamic resistance to momentum transfer and the additional boundary layer resistance for heat and mass transfer (or excess resistance) (Thom, 1972; Knauer et al., 2018a).

$$\begin{array}{}\text{(9)}& r_{\mathrm{a}}={\displaystyle \frac{\ln \left[\frac{z_{\mathrm{r}}-d}{z_{\mathrm{0}\mathrm{m}}}\right]-{\mathit{\psi }}_{\mathrm{h}}}{k{u}_{\ast }}}+\mathrm{6.2}{u}_{\ast }^{-\mathrm{0.67}}\end{array}$$

The first term on the right-hand side of Eq. (9) is the aerodynamic component, and the second term is the boundary layer component. Here, u_∗ is friction velocity, k is the von Kármán constant (0.41), d is the zero-plane displacement height (d=0.7z_h), z_0 m is the roughness length for momentum (z_0 m=0.1z_h), and ψ_h is the integrated form of the stability correction function. z_h is canopy height based on manual measurements taken during regular maintenance visits. r_a was estimated using the bigleaf R package (Knauer et al., 2018a).

By rearranging Eq. (2), RH_s can be calculated using

$$\begin{array}{}\text{(10)}& {\mathrm{RH}}_{\mathrm{s}}={\displaystyle \frac{\mathit{\gamma }\mathit{LE}{r}_{\mathrm{a}}/\mathit{\rho }{c}_p+e_{\mathrm{a}}}{SH{r}_{\mathrm{a}}/\mathit{\rho }{c}_p+e^{\ast }\left(T_{\mathrm{a}}\right)}}.\end{array}$$

Negative H and inaccurate r_a modelling sometimes yielded negative RH_s or values greater than 1, especially at nighttime. In these cases, RH_s was assigned the value of 1 following the approach described in the bigleaf R package (Knauer et al., 2018a). We then estimated LE_Q and LE_G from Eq. (4).

In order to explore the timescale of the covariances for LE ∼ LE_Q and LE ∼ LE_G in the frequency domain, we applied wavelet coherence analysis using the WaveletComp R package (Roesch and Schmidbauer, 2014). The package is designed to apply the continuous wavelet transform using the Morlet wavelet, which is a popular approach to analyze hydrological and micrometeorological datasets (Hatala et al., 2012; Johnson et al., 2013). A total time series of half-hourly decomposed LE for the 4-year measurement period was used to estimate localized coherence and phase angle. The wavelet coherence can be interpreted as the local correlation between two variables in the frequency–time domain (where red indicates high correlation). A 0^∘ phase angle (arrow pointing right) indicates periods of positive correlation, while a 180^∘ phase angle (arrow pointing left) indicates periods of negative correlation.

3.2 FLUXNET2015

The daily-scale FLUXNET2015 dataset, which includes 212 empirical eddy-covariance flux tower sites around the globe (Pastorello et al., 2020), was used in this study. The turbulent heat fluxes, net radiation, soil heat flux, air temperature, relative humidity, wind speed, friction velocity, and barometric pressure were obtained from the dataset. For this analysis, we only included daily data for periods for which the quality control flag indicated more than 80 % half-hourly data were present (i.e., measured data in general, or good-quality gap-filled data in cases of partially missing data).

In order to decompose daily LE into LE_Q and LE_G, we estimated daily aerodynamic resistance (r_a) by Eq. (11) instead of Eq. (9) since canopy and measurement heights are unknown (Thom, 1972; Knauer et al., 2018a).

$$\begin{array}{}\text{(11)}& r_{\mathrm{a}}={\displaystyle \frac{u_{\ast }^{\mathrm{2}}}{u\left(z_{\mathrm{r}}\right)}}+\mathrm{6.2}{u}_{\ast }^{-\mathrm{0.67}},\end{array}$$

where u(z_r) is reference height wind speed. r_a was estimated using the bigleaf R package (Knauer et al., 2018a), and RH_s was calculated from Eq. (8).

LE_Q and LE${}_Q^{\prime }$ were calculated using RH_a and RH_s following Eqs. (4) and (5), and then LE_G and LE${}_G^{\prime }$ were calculated by subtracting LE_Q and LE${}_Q^{\prime }$ from LE. To calculate LE_Q and LE${}_Q^{\prime }$, we define Q as LE + H, but it should be noted that this approach can include systematic uncertainty since the sum of LE and H measured by eddy covariance is typically lower than R_n−G (i.e., conditions referred to as the energy balance closure problem; Wilson et al., 2002). To investigate the effect of a lack of energy balance closure on resulting LE terms, we provide Fig. S2 that was generated by (1) defining Q as R_n−G and (2) correcting LE and H based on the assumption that the Bowen ratio ($B=H/\mathit{LE}$) is correct (Pastorello et al., 2020).

3.3 FLUXCOM

The FLUXCOM dataset (Jung et al., 2019) is a global-scale machine learning ensemble product which upscales FLUXNET observations (Baldocchi et al., 2001) using Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data and reanalysis meteorological data. In this study we used the monthly LE FLUXCOM dataset (0.5^∘ resolution) modelled using MODIS and ECMWF ERA5 reanalysis data (Hersbach et al., 2020).

We obtained Q and LE from the FLUXCOM output, and air temperature and dew point temperature were retrieved from ERA5 monthly averaged data (2 m height). RH, S, and γ were calculated from ERA5 data, and then LE${}_Q^{\prime }$ was calculated. LE${}_G^{\prime }$ was then estimated by subtracting LE${}_Q^{\prime }$ from LE.

4 Results

4.1 Decomposition analysis of in situ EC flux observation

Application results of the PM_RH model to the observed LE at an irrigated sugarcane farm in Costa Rica are depicted in Fig. 2. The decomposition analysis of observed LE shows that while LE_Q is the major component of LE, LE_G variability plays a non-negligible role in seasonal and interannual behaviour of LE. In terms of absolute magnitude, the LE_Q term can closely approximate LE, and LE_G only represents 15 % of total evaporation (Fig. 2c). Also, positive coherence between LE and LE_Q was strong over the entire period of observation, particularly at diurnal to multiday timescales (0.5–32 d), implying variability of LE is largely determined by LE_Q variability (i.e., red coloured regions in Fig. 2e).

Although absolute magnitude of LE_G was much smaller than that of LE_Q, the interannual variability of LE_G was larger than the interannual variability of LE_Q (Fig. 2c). Furthermore, LE and LE_G also had a strong positive correlation on longer timescales (32–365 d) (i.e., red coloured regions in Fig. 2f). Unexpectedly, a negative correlation between LE and LE_G at the diurnal timescale was observed in the wavelet analysis only when the land surface was dry and there was little vegetation (i.e., after harvest) or during a year in which there was no dry season irrigation applied. Also, we found that EF variability is mostly determined by LE_G variability since the diurnal and seasonal signals of Q are removed from LE in EF. Interestingly, the annual mean LE_G was the highest in 2015, a drought year in which RH_a and precipitation were generally lower than for the other years, while the annual mean LE_G was close to zero in 2016 when there was no application of dry season irrigation due to crop replanting.

Figure 3 Half-hourly time series of heat fluxes and the two components indicated in Fig. 2d. The background colour represents RH_s. Here, RH_a is mean atmospheric relative humidity, SM is volumetric soil water content, normalized difference vegetation index (NDVI) represents vegetation status, and z_h is canopy height. Panel (e) presents the long-term mean diurnal cycle of decomposed LE (dots) and EF (dashed line).

To explore the diurnal behaviour of decomposed LE, we selected different surficial and atmospheric conditions when LE_G was zero, positive, or negative in Fig. 3. In the 2016 dry season, LE_G was close to zero as a daily average value, as a result of negative daytime and positive nighttime LE_G values due to dry air and dry soil conditions (no irrigation) and an undeveloped vegetation canopy (Fig. 3a). Daily LE_G was also close to zero during wet season conditions (e.g., Fig. 3b). In this case, LE_G was near zero during both daytime and nighttime periods due to near-saturated atmospheric and land surface conditions. These two cases show that “dry land–dry air” or “wet land–wet air” conditions can each lead to daily scale land–atmosphere equilibrium, although the diurnal pattern of LE_G is starkly different for dry land–dry air vs. wet land–wet air conditions.

Meanwhile, when RH_a was low and the canopy was well-developed, LE_G was found to be positive during both daytime and nighttime periods (Fig. 3c). On the other hand, during post-harvest conditions when vegetative canopy cover was minimal and air and soil moisture levels were low, daily LE_G was found to be negative as a result of negative daytime and positive nighttime LE_G (Fig. 3d). Diurnal variation in RH_s was maximized when daily LE_G was negative, implying a large diurnal fluctuation in surface temperature under the drier land surface conditions. Regarding the overall diurnal pattern, LE_G generally declined during the morning and increased in the afternoon, which is consistent with the well-known diurnal pattern of EF (Gentine et al., 2011, 2007) (Fig. 3e).

4.2 Decomposition analysis of the FLUXNET2015 dataset

Decomposition analysis of the daily FLUXNET2015 dataset is illustrated in Fig. 4. In terms of absolute magnitude of each term, the majority of LE_G (and LE${}_G^{\prime }$) values ranged from −50 to 50 W m⁻², with some exceptional values approaching ±100 W m⁻². On the other hand, LE_Q (and LE${}_Q^{\prime }$) values ranged from 0 to 150 W m⁻².

Figure 4 FLUXNET2015 daily-scale decomposed LE for 212 sites and 1532 site years. Panels (a1) and (a2) are linear regressions of LE${}_Q^{\prime }$ on LE_Q and LE${}_G^{\prime }$ on LE_G. Panels (b1) and (b2) are linear regressions of LE on LE_Q and LE on LE_G. Panels (c1) and (c2) are linear regressions of EF on LE${}_Q/Q$ and EF on LE${}_G/Q$. In these panels, daily EF data within a range from −1 to 1.5 are only shown. Here, dashed lines are one-to-one lines, blue lines are regression lines, and colour represents RH_s − RH_a. Panels (d1) and (d2) are histograms of decomposed LE and EF with mean values (dotted lines). To correct for lack of energy balance closure, Q was set equal to LE + H in all calculations.

One of the interesting findings from the decomposition analysis of the FLUXNET2015 dataset was that differences between LE_Q and LE${}_Q^{\prime }$ are marginal at a daily timescale (the slope is close to 1 in Fig. 4a1). This result implies that although the diabatic processes expressed by Eqs. (4) and (5) are different in magnitude due to the difference between RH_s and RH_a (see Sect. 2.3), these differences are practically negligible. This is an important point since LE${}_Q^{\prime }$ can be determined simply and directly using reference height meteorological measurements, while LE_Q is required to know RH_s.

As for the adiabatic terms, LE${}_G^{\prime }$ is roughly 1.1 times LE_G at a daily timescale (Fig. 4a2), which is consistent with the theory regarding their respective thermodynamic paths. As we discussed in Sect. 2.3, the absolute magnitude of LE${}_G^{\prime }$ must be bigger than that of LE_G when available energy is positive. Therefore, the empirical relationship between LE_G and LE${}_G^{\prime }$ in Fig. 4a2 is a consequence of a physical principle, and this result may provide the following empirical relationship.

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& \mathrm{1.1}{\displaystyle \frac{\mathit{\rho }{c}_p{e}^{\ast }\left(T_{\mathrm{a}}\right)}{{\mathrm{RH}}_{\mathrm{s}}S+\mathit{\gamma }}}\cdot {\displaystyle \frac{{\mathrm{RH}}_{\mathrm{s}}-{\mathrm{RH}}_{\mathrm{a}}}{r_{\mathrm{a}}}}\approx {\displaystyle \frac{\mathit{\rho }{c}_p{e}^{\ast }\left(T_{\mathrm{s}}\right)}{{\mathrm{RH}}_{\mathrm{a}}S+\mathit{\gamma }}}\cdot {\displaystyle \frac{{\mathrm{RH}}_{\mathrm{s}}-{\mathrm{RH}}_{\mathrm{a}}}{r_{\mathrm{a}}}}\\ & \mathrm{1.1}{\displaystyle \frac{e^{\ast }\left(T_{\mathrm{a}}\right)}{{\mathrm{RH}}_{\mathrm{s}}S+\mathit{\gamma }}}\approx {\displaystyle \frac{e^{\ast }\left(T_{\mathrm{s}}\right)}{{\mathrm{RH}}_{\mathrm{a}}S+\mathit{\gamma }}}\end{array}\end{array}$$

Equation (12) reveals an emergent daily timescale relationship between temperature and relative humidity which has the potential to be used as a supplementary equation in future research.

Another important finding of the decomposition analysis is the global-scale land–atmosphere equilibrium. Our analysis in Fig. 4d1 indicates that the mean value of daily LE_G of all FLUXNET2015 sites is close to zero, implying the global mean RH gradient is near zero at a daily timescale. Importantly, LE is primarily determined by LE_Q (R²=0.65) instead of LE_G (R²=0.18) as depicted in Fig. 4b1 and b2. Nevertheless, FLUXNET2015 data also suggest that LE_G is the main driver of local-scale variability of EF at the daily timescale (Fig. 4c1 and c2). Although the mean value of daily EF is close to the mean value of $\frac{{\mathit{LE}}_Q}{Q}$, the variation in daily EF depends more on the variation in $\frac{{\mathit{LE}}_G}{Q}$ (Fig. 4d2). It should be noted that Figs. 4 and S2 are almost identical (Fig. S2 repeats the presentation shown in Fig. 4 using the value computed when forcing energy balance closure), implying that the lack of surface energy balance closure for EC observations does not significantly impact our analyses and interpretations.

4.3 Decomposition analysis of FLUXCOM dataset

We then applied the PM_RH model to the FLUXCOM dataset, a benchmark global LE data product (Jung et al., 2019). As shown in Fig. 5a and c, the spatial patterns of the annual mean LE and LE${}_Q^{\prime }$ were similar. For instance, both LE and LE${}_Q^{\prime }$ show the highest values around the Equator at an annual timescale, which is mainly due to the energy available in this region. Also, spatial variability of LE is mostly determined by LE${}_Q^{\prime }$ (R²=0.85 and slope = 1) rather than by LE${}_G^{\prime }$ (R²=0.18) (Fig. 5f1 and f2). This result is consistent with Eq. (8) and the SFE theory. In other words, the land surface is generally under thermodynamic equilibrium with the atmosphere at the global–annual scale (i.e., RH_s≈RH_a). Furthermore, the monthly time series of global LE and its two components in Fig. 5e1 show that (i) LE${}_G^{\prime }$ is consistently close to zero at the global scale and (ii) the seasonal variability of global LE is primarily determined by LE${}_Q^{\prime }$.

Figure 5 Mean annual LE, EF, LE${}_Q^{\prime }$, and LE${}_G^{\prime }$ from 2001 to 2005 (panels a, b, c, and d, respectively). Panel (e1) is a time series of monthly global average LE and the two components LE${}_G^{\prime }$ and LE${}_Q^{\prime }$. Panels (e2), (e3), and (e4) are time series at specific locations highlighted in panel (d). Panels (f1), (f2), (f3), and (f4) are spatial linear regressions of LE on LE${}_Q^{\prime }$, LE on LE${}_G^{\prime }$, EF on LE${}_Q^{\prime }/Q$, and EF on LE${}_G^{\prime }/Q$, respectively.

However, while mean annual LE${}_G^{\prime }$ was close to zero in broad areas (particularly in high-latitude regions) as exemplified in Fig. 5e2, it was distinctly positive or negative at the annual scale for many regions (Fig. 5d). In humid tropical regions like the Amazon basin where moisture convergence is large, LE${}_G^{\prime }$ was generally positive, whereas arid regions such as Australia were characterized by negative LE${}_G^{\prime }$ (Fig. 5e3 and e4). Here, positive LE${}_G^{\prime }$ (i.e., RH_s >  RH_a) indicates the land surface is wetter than the near-surface atmosphere while negative LE${}_G^{\prime }$ (i.e., RH_s < RH_a) implies a drier land surface than the atmosphere. Therefore, the sign of LE${}_G^{\prime }$ in Fig. 5d can be interpreted as representing land surface dryness relative to the atmosphere at an annual timescale.

The spatial pattern of LE${}_G^{\prime }$ is similar to the spatial pattern of EF, but differs from the spatial pattern of LE${}_Q^{\prime }$ (Fig. 5b). For example, EF was high in not only tropical regions but also temperate climates such as Mediterranean regions, and this spatial pattern is well matched with the spatial pattern of LE${}_G^{\prime }$ but not LE${}_Q^{\prime }$. The finding that the spatial variation in EF is primarily controlled by LE${}_G^{\prime }$ instead of LE${}_Q^{\prime }$ was supported by correlation analyses in Fig. 5f3 and f4 (R²=0.60 for EF ∼ $\frac{{\mathit{LE}}_G^{\prime }}{Q}$ and R² = 0.28 for EF ∼ $\frac{{\mathit{LE}}_Q^{\prime }}{Q}$). This is understandable in that EF is a reflection of the land surface dryness (Gentine et al., 2011), and LE${}_G^{\prime }$ reveals the land surface dryness relative to the atmosphere.

5 Discussion

5.1 LE_G at sub-daily scale

Salvucci and Gentine (2013) found that the variance of the RH gradient tends to be minimized over the course of the day. Based on this empirical finding, they developed an approach to predict LE only using standard meteorological measurements, and this approach accurately predicted actual LE (Rigden and Salvucci, 2015, 2017). Our PM_RH model provides theoretical support for their approach in that LE_G acts to reduce the RH gradient. Indeed, the U-shape diurnal cycles of LE_G in Fig. 3 (positive nighttime and negative daytime) show the direction and the magnitude of the equilibration process of RH gradient at a sub-daily scale resulting in a small gradient of RH on daily average.

We found positive nighttime LE_G values regardless of wet or dry conditions, except when the atmosphere was fully saturated (Fig. 3). This result is a natural consequence since the land surface is close to saturation at night. This finding suggests that LE_G is a dominant contributor to nighttime evaporation since available energy is close to zero at night. It is important since nighttime evaporation is a non-negligible component of total ET (Padrón et al., 2020).

Unlike nighttime LE_G values, the direction and the magnitude of the daytime LE_G are highly dependent on the land surface dryness. For example, the U-shape diurnal cycles of LE_G are apparent only when the land surface is dry, which is confirmed by the negative wavelet coherence between LE and LE_G at the diurnal timescale in Fig. 2f. When the land surface is wetter than the atmosphere, LE_G values were positive even in the daytime, and thus the U-shape diurnal cycles of LE_G did not appear (Fig. 3c). The positive LE_G during daytime periods may be explained as a consequence of warm and dry air entrainment at the top of the atmospheric boundary layer and/or horizontal advection of sensible heat which may reduce atmospheric relative humidity (Baldocchi et al., 2016; de Bruin et al., 2016). Indeed, a strong entrainment effect and/or local advection of sensible heat are common phenomena for irrigated agriculture (Baldocchi et al., 2016; de Bruin and Trigo, 2019), and the irrigated sugarcane site shows that the annual mean LE_G was always positive except for in 2016 when there was no application of dry season irrigation.

5.2 LE_G at daily to annual timescales

As described in the theory section, land–atmosphere equilibrium is achieved when LE_G approaches zero and thus LE reduces to Eq. (8) at a timescale longer than sub-daily. The decomposed terms derived from both the empirical FLUXNET2015 and model-based FLUXCOM datasets show that the global mean for LE_G is near zero, implying global-scale land–atmosphere equilibrium (Figs. 4 and 5). This result extends the SFE theory of McColl et al. (2019). Although LE_G is not always near zero for timescales longer than sub-daily, moisture convergence and divergence at the global scale tend to balance each other out, resulting in global-scale land–atmosphere equilibrium on longer timescales (Fig. 5e1).

From a different perspective, non-zero LE_G value and its sign (+ vs. −) at the local scale can be understood as an indicator reflecting land surface dryness relative to the atmosphere. We found that LE_G clearly distinguishes spatially wet and dry regions around the world (Fig. 5 d). We also found that the spatiotemporal variability of EF was largely explained by LE_G instead of LE_Q (Figs. 4c2 and 5f4). These results demonstrate the usefulness of LE_G to quantify land surface dryness. Some previous studies introduced evaporative stress index or evaporation deficit index based on the ratio (or difference) between potential evaporation and actual evaporation (Anderson et al., 2015; Kim and Rhee, 2016; Fisher et al., 2020; Baldocchi et al., 2021), but these methods are highly dependent on the way one calculates potential evaporation. Unlike potential evaporation (which is a theoretical value), LE_G is a true physical quantity, and negative LE_G values directly reflect water-limited land surface conditions. Therefore, we suggest applying our decomposition method to better quantify evaporative stress.

5.3 Future applications

In this study, we present the PM_RH model and demonstrate its utility for exploratory and diagnostic purposes. However, the model has potential applications for other purposes. One possible application is to use PM_RH to predict actual ET. Although the original PM equation is widely used to predict evapotranspiration (e.g., Leuning et al., 2008; Mu et al., 2011; Mallick et al., 2014), its accuracy often relies on parameterized surface resistance models (Polhamus et al., 2013). Since our PM_RH formulation does not include surface resistance, it could represent a good alternative. As shown in the results section, LE${}_Q^{\prime }$ can be calculated using typical meteorological data without additional surface parameters. Also, we found that LE${}_Q^{\prime }$ alone can be used effectively to approximate LE. If the remotely sensed land surface conditions are known (e.g., soil moisture and/or land surface temperature), actual LE may be more accurately predicted by incorporating the LE_G term. To estimate the LE_G term, RH_s may be estimated based on soil moisture and/or land surface temperature data. For example, Eq. (12), which is well supported by observations (Fig. 4a2) and on a physical basis, may be used to calculate RH_s.

Another possible application is to study impacts of climate change and land use and land cover change on surface energy partitioning and evaporation. Changes in the atmospheric state such as temperature and relative humidity, as well as changes in the land surface characteristics such as albedo and aerodynamic roughness, can alter evaporation (Lee et al., 2011; Wang et al., 2018). However, how these changes affect terrestrial energy partitioning and ET is still unclear. Indeed, there is a large discrepancy in long-term LE trends among current land surface models (Pan et al., 2020). Since PM_RH makes it possible to physically decompose LE into adiabatic and diabatic thermodynamic components, the PM_RH approach can be useful to understand how environmental changes affect surface energy partitioning. For instance, trend analysis for the decomposed LE using PM_RH could contribute to improve our understanding of earth's climate system and water cycle.

5.4 Potential caveats

Despite the insights it can offer, the PM_RH model shares several limitations with the traditional PM model. First, PM-style equations linearize the exponential relationship between saturation vapour pressure and temperature (Clausius–Clapeyron relation), but the linearization can cause bias when the temperature difference between surface and atmosphere is substantial (Paw U and Gao, 1988; McColl, 2020). Second, the PM_RH model assumes that aerodynamic resistance for heat and water vapour is identical, which implicitly relies on the assumption that the ratio of the turbulent Schmidt to Prandtl numbers is unity (Knauer et al., 2018a). This similarity assumption cannot be held in some cases, especially when advective flux divergence is significant (Lee et al., 2004). The third potential caveat is that we define the land surface as an idealized single plane that is equivalent to the “bigleaf” representation of the traditional PM equation, but this approach ignores profiles of temperature and humidity inside the canopy (Bonan et al., 2021).

Another potential caveat concerns the surface energy balance. The surface energy balance is a governing equation of the PM-style models, but it is not satisfied in typical EC observations, which is referred to as the “surface energy balance closure problem” (Wilson et al., 2002). This closure problem could cause a systematic uncertainty in estimating r_s when using the PM equation (Knauer et al., 2018b; Wohlfahrt et al., 2009), and this issue may affect the diagnostic analyses using the PM_RH. In this study, we did not force the energy balance closure and attributed the cause of observed surface energy imbalances to unmeasured heat storage terms for the Costa Rica EC site due to the possibly significant role of the heat storage term (details in the Sect. 3.1). Wehr and Saleska (2021) recently demonstrated that regardless of whether the lack of energy balance closure of EC observations is due to LE+H or is due to R_n−G, applying the flux gradient equation to observed LE and H without energy balance correction is the best way to determine r_s. This is because applying the flux gradient equation to observed LE and H can dispense with the unnecessary assumption of energy balance closure. They showed that bias introduced by underestimated LE and H is smaller than the bias introduced by the energy balance closure assumption. This finding may be applied to our analysis in calculating RH_s instead of r_s.

As for the FLUXNET2015 dataset, we provide an alternate analysis using energy-balance-corrected LE and H (Bowen ratio preserving method in Pastorello et al., 2020) in the Supplement. We found that the results for corrected and uncorrected versions were almost identical, which can be viewed as a natural consequence since in Eq. (10) LE and H are included in the numerator and denominator respectively. Multiplying the same ratio to LE and H in Eq. (10) to correct LE and H based on the Bowen ratio method does not significantly change the resulting RH_s. Therefore, the lack of surface energy balance closure does not significantly impact our analyses and interpretations unless the lack of energy balance is dominated by LE only or H only. If the lack of energy balance is dominated by LE only or H only, our results and interpretation may include systematic bias.

Finally, there are several ways to calculate aerodynamic resistance, and the chosen form for r_a may affect the results. However, the influence of this choice is expected to be marginal compared to the energy balance problem. Knauer et al. (2018b) showed that uncertainty caused by different r_a values on surface conductance is low compared to the energy balance closure problem. This finding can be applied to our analysis. Specifically, in Eq. (10), r_a is multiplied by both denominator and numerator, and thus a small difference in r_a should not significantly affect the resulting RH_s.

6 Conclusions

We have shown that our novel PM_RH model provides a new opportunity to understand the governing physics of the terrestrial energy budget. Specifically, the PM_RH model helps to illustrate how the land surface conditions become encoded to the atmospheric state by partitioning LE into two thermodynamic processes. “Dry land–dry air” or “wet land–wet air” conditions can each lead to daily scale land–atmosphere equilibrium, although the diurnal pattern of the equilibration process (i.e., LE_G) is starkly different. Our findings suggest that while LE_G is a primary component determining EF, spatiotemporal variability of LE_Q alone can adequately represent the variability of LE. We found global-scale land–atmosphere equilibrium at daily to annual scales, which implies that global LE can be simply determined by the atmospheric state and radiative energy without any surface constraint required to represent spatial heterogeneity and physiological influences. From a different perspective, the non-zero LE_G value at a local scale can be understood as an indicator revealing land surface dryness. Questions remain regarding how LE_Q and LE_G will be influenced in relation to changing climatic and land surface conditions and how these changes might affect the climate system at differing spatial and temporal scales through positive or negative feedbacks.