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

. Earth’s climate and water cycle are highly dependent on terrestrial evapotranspiration and the associated ﬂux 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 ﬂux 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.


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;Wil-lett et al., 2014;Byrne and O'Gorman, 2018), but landatmosphere 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 masstransfer 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).
where S is the linearized slope of saturation vapour pressure versus temperature (hPa K −1 ), γ is the psychrometric constant (hPa K −1 ), ρ is the air density (kg m −3 ), c p is the specific heat capacity of air at constant pressure (MJ kg −1 K −1 ), 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 −2 ). 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 * (T a ) − e 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 −1 ), and r s is surface resistance to water vapour transfer (s m −1 ) 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 landatmosphere 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.

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 LE = ρc p γ e s −e a r 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 s = e * (T s )). For an unsaturated land surface, however, relative humidity should be introduced as e s = RH s e * (T s ), 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 LE = ρc p γ RH s e * (T s )−e a r a for an unsaturated surface condition.
In order to decompose LE into two individual fluxes related to temperature and relative humidity gradients, we add −RH s e * (T a ) + RH s e * (T a ) (or −RH a e * (T s ) + RH a e * (T s )) to the numerator of the flux gradient equation and rewrite LE as follows.
We then approximate e * (T s ) − e * (T a ) = S(T s − T 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 = ρc p T s −T a r 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: Adiabatic process: Adiabatic process: where LE Q (and LE Q ) is a diabatic component, and LE G (and LE G ) 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 ) and its adiabatic component (LE G or LE G ). In the following sections, we will discuss theoretical meanings of Eqs. (4) and (5) in depth.

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: We rearrange this formulation to derive Eq. (7) by factoring out e * (T a ) and introducing RH a = e a e * (T a ) into the second term.
Adiabatic process (7) 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 ) 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 ) term can be understood as equilibrium LE for an unsaturated surface, which we discuss later in Sect. 2.4.

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.
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 ), and the adiabatic process is followed at temperature of T S (i.e., LE G ). 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 when RH s > RH a . Also, the absolute magnitude of LE G is always bigger than that of LE G when Q > 0 (i.e., vector B → C is longer than vector A → B in Fig. 1).

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 ) is zero and thus LE becomes We note that Eq. (8) 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 ) 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., RH s ≈ RH a ≈ 1), Eq. (8) becomes classical equilibrium LE (i.e., LE ≈ S S+γ 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;Mc-Coll, 2020). Therefore, we can regard Eq. (8) as a generalized equilibrium LE for an unsaturated surface. Here, the enthalpy change of the air parcel is defined as Q·r a ρ (kJ kg −1 ). 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.

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) 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.

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;. The site was irrigated occasionally during dry sea-sons 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.
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., LE+H R 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).
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 RH s = γ LEr a /ρc p + e a SH r a /ρc p + e * (T a ) .
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 follow-ing 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.

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 gapfilled 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).
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 were calculated using RH a and RH s following Eqs. (4) and (5), and then LE G and LE G were calculated by subtracting LE Q and LE Q from LE. To calculate LE Q and LE Q , 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 /LE) is correct (Pastorello et al., 2020).

FLUXCOM
The FLUXCOM dataset (Jung et al., 2019) is a globalscale 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 was calculated. LE G was then estimated by subtracting LE Q from LE.

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.
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  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). 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 welldeveloped, 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(Gentine et al., , 2007 (Fig. 3e).

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 ) values ranged from −50 to 50 W m −2 , with some exceptional values approaching ±100 W m −2 . On the other hand, LE Q (and LE Q ) values ranged from 0 to 150 W m −2 . One of the interesting findings from the decomposition analysis of the FLUXNET2015 dataset was that differences between LE Q and LE Q 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 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 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 must be bigger than that of LE G when available energy is positive. Therefore, the empirical relationship between LE G and LE G in Fig. 4a2 is a consequence of a physical principle, and this result may provide the following empirical relation-  (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. ship.
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 2 = 0.65) in-stead of LE G (R 2 = 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 LE Q Q , the variation in daily EF depends more on the variation in 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.

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 were similar. For instance, both LE and LE Q 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 (R 2 = 0.85 and slope = 1) rather than by LE G (R 2 = 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 is consistently close to zero at the global scale and (ii) the seasonal variability of global LE is primarily determined by LE Q .
However, while mean annual LE G 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 was generally positive, whereas arid regions such as Australia were characterized by negative LE G (Fig. 5e3  and e4). Here, positive LE G (i.e., RH s > RH a ) indicates the land surface is wetter than the near-surface atmosphere while negative LE G (i.e., RH s < RH a ) implies a drier land surface than the atmosphere. Therefore, the sign of LE G 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 is similar to the spatial pattern of EF, but differs from the spatial pattern of LE Q (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 but not LE Q . The finding that the spatial variation in EF is primarily controlled by LE G instead of LE Q was supported by correlation analyses in Fig. 5f3 and f4 (R 2 = 0.60 for EF ∼ LE G Q and R 2 = 0.28 for EF ∼ LE Q Q ). This is understandable in that EF is a reflection of the land surface dryness (Gentine et al., 2011), and LE G reveals the land surface dryness relative to the atmosphere. 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 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 equilibra-tion 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).

LE G at sub-daily scale
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.

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 subdaily. 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 evapora- tive 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.

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 can be calculated using typical meteorological data without additional surface parameters. Also, we found that LE Q 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.

Potential caveats
Despite the insights it can offer, the PM RH model shares several limitations with the traditional PM model. First, PMstyle 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 .

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 landatmosphere 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.
Author contributions. YK and MSJ designed research; UW provided FLUXCOM data; YK, LM, and MSJ performed research; YK analyzed data; YK, MG, TAB, LM, and MSJ wrote the paper.
Competing interests. The authors declare that they have no conflict of interest.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.