Satellite‐Based Assessment of Land Surface Energy Partitioning–Soil Moisture Relationships and Effects of Confounding Variables

Land surface energetic partitioning between latent, sensible, and ground heat fluxes determines climate and influences the terrestrial segment of land‐atmosphere coupling. Soil moisture, among other variables, has a direct influence on this partitioning. Dry surfaces characterize a water‐limited regime where evapotranspiration and soil moisture are coupled. This coupling is subdued for wet surfaces, or an energy‐limited regime. This framework is commonly evaluated using the evaporative fraction–‐soil moisture relationship. However, this relationship is explicitly or implicitly prescribed in land surface models. These impositions, in turn, confound model‐based evaluations of energetic partitioning‐–soil moisture relationships. In this study, we use satellite‐based observations of surface temperature diurnal amplitude (directly related to available energy partitioning) and soil moisture, free of model impositions, to estimate characteristics of surface energetic partitioning–‐soil moisture relationships during 10–‐20‐day surface drying periods across Africa. We specifically estimate the spatial patterns of water‐limited energy flux sensitivity to soil moisture (m) and the soil moisture threshold separating water and energy‐limited regimes (θ*). We also assess how time evolution of other factors (e.g., solar radiation, vapor pressure deficit, surface albedo, and wind speed) can confound the energetic partitioning–‐soil moisture relationship. We find higher m in drier regions and interestingly similar spatial θ* distributions across biomes. Vapor pressure deficit and insolation increases during drying tend to increase m. Only vapor pressure deficit increases in the Sahelian grasslands systematically decrease θ*. Ultimately, soil and atmospheric moisture availability together play the largest role in land surface energy partitioning with minimal consistent influences of time evolution of other forcings.


Introduction
The water, carbon, and energy cycles are coupled at the land surface principally through the soil moisture and soil temperature state variables. These state variables can also restore perturbations or prolong them. For example, soil moisture modulated land-atmosphere coupling has implications for feedbacks that drive rainfall or the lack thereof and can lead to persistence of droughts and heatwaves (Miralles et al., 2019;Rodriguez-Iturbe et al., 1991;Santanello et al., 2018). Under water-limited conditions, soil moisture (θ) directly controls evapotranspiration (ET) and surface temperature, driving a positive land-atmosphere feedback loop (Seneviratne et al., 2010). In this regime, higher θ increases ET which is coupled with greater near surface humidity and a reduction in surface temperature, cloud base, and dry air entrainment into the boundary layer (Betts, 2004;Betts et al., 1996;Entekhabi et al., 1996;Seneviratne et al., 2010). There is general agreement that these conditions can increase the probability of rainfall keeping the surface wet (Findell et al., 1997;Findell & Eltahir, 2003;Koster et al., 2004;Tuttle & Salvucci, 2016), though caveats of convective triggering under dry surface anomalies are noted (Findell & Eltahir, 2003;Taylor et al., 2012Taylor et al., , 2011. Conversely, a drying surface without moisture replenishment leads to less evaporative cooling. Corresponding increases in surface albedo further inhibit vertical motion and rainfall (Charney, 1975), and the lower atmosphere ultimately dries and warms. These feedbacks lead to persistence of drought and heatwaves (Hirschi et al., 2011;Miralles et al., 2019;Roundy et al., 2012), until the system is sufficiently shocked by exogenous (generally advective) restoring forces (Brubaker & Entekhabi, 1996). At higher mean θ availability, the system is considered energy-limited as more moisture does not necessarily lead to greater ET. Here, the water, carbon, and energy cycle couplings are subdued. ET is at or near its potential value where net radiation and atmospheric resistance, for example, are instead limiting.
These coupled processes are typically evaluated using the relationship between evaporative fraction (EF) and θ where EF is (1) where LE is latent heat flux, H is sensible heat flux, G is ground heat flux, and R n is net radiation. The available net radiative energy at the surface can be dissipated via LE, H, and G in equation (2). All flux terms have positive signs when away from the surface. EF is therefore the fraction of available surface energy partitioned into evaporating water (LE) instead of heating the overlying air (H). The canonical EF-θ relationship is characterized by a transition point in soil moisture (θ*; some other studies label θ Crit ; Seneviratne et al., 2010) separating the water (Stage II ET) and energy-limited (Stage I ET) regimes (Budyko, 1974;Eagleson, 1978). A conceptual drawing of this function is shown in Figure 1. Below θ*, EF and θ are positively related Figure 1. Schematic of the typical relationship between evaporative fraction and soil moisture. A third dimension indicates that aspects of this relationship could change with other factors. ET=evapotranspiration; Rs= insolation; VPD=vapor pressure deficit; ∂=surface albedo; θ=soil moisture. (Teuling et al., 2006) where reduced θ lowers EF. The relationship may be nonlinear, but the first-order linear term of the function across this segment below θ* is indicative of the degree of coupling between the surface water and energy balances. Above θ*, EF is decoupled from soil moisture, and hence, its functional relationship is a flat line with an approximately zero slope. Here, EF is an influenced by other variables such as net radiation (Gentine et al., 2007).
θ* itself is largely defined by soil texture and plant physiology. Specifically, θ* defines the θ value below which soil water controls the combined rates of bare soil evaporation and water loss through stomata (Eagleson, 1978;Seneviratne et al., 2010). In other words, ET drops below its potential or energy-limited value. At a given location, the probability density function of θ together with θ* estimate describes the time spent in each hydrologic regime, key for determining the existence of land-atmosphere feedbacks. While θ* divides the system into water and energy-limited states of the system, ∂EF ∂θ (the sensitivity of EF to θ below θ*) determines the magnitude of θ control on the boundary layer (Dirmeyer, 2011). The local slope is only a potential metric of influence, however, as θ variations must also be of sufficient magnitude to lead to variability in boundary layer processes (Dirmeyer, 2011;Guo et al., 2006). Together, both bare soil and vegetation components contribute to the functional behavior below θ*. For bare soils, greater θ leads to less soil capillary pressure retention and more efficient atmospheric extraction of soil vapor. For vegetated surfaces, greater θ is associated with plants keeping stomatal apertures open to transpire and assimilate carbon. It is thus no surprise that the general relationship between EF and θ in Figure 1 is also observed between primary productivity and θ, evidence for coupling to the carbon cycle (Short Gianotti et al., 2019). Ultimately, accurate representation of surface energy partitioning in land surface models is essential for weather forecasting and climate projections of ET and gross primary productivity.

Surface Energy Partitioning-Soil Moisture Relationship: Model Assessments
Energetic partitioning-θ relationships remain largely unobserved due to the challenge of directly measuring surface energy fluxes and θ at climate model grid scales across different climates and biomes. These relationships have primarily been evaluated and visualized via EF-θ within model frameworks (Dirmeyer et al., 2000;Koster et al., 2009;Schwingshackl et al., 2017). These studies revealed much about vegetation dependencies, inter-model differences, and land-atmosphere interactions. However, any model-based evaluation of EF-θ is subject to a combination of explicit and/or implicit impositions that ultimately formulate their relationship (Dirmeyer et al., 2007;Koster et al., 2006). Any parameterization that translates even a purely observed θ to an ET estimate is based on land surface energy partitioning assumptions that impose θ* and ∂EF ∂θ , for instance. Model-free observations of this relationship are required to fully evaluate climate and vegetation dependencies and are critical for diagnosing land-atmosphere feedbacks in models (Dirmeyer et al., 2007;Guo et al., 2006). Observed temperature-precipitation relationships (with measurements of each gridded at large spatial scales) were used to compare with or serve as a proxy for EF-θ relationships (Gevaert et al., 2018;Koster et al., 2009). However, differences between EF and temperature and precipitation and θ hinder one-to-one comparison between these relationships.

Surface Energy Partitioning-Soil Moisture Relationship: Observations
The most direct observations of energetic partitioning can be obtained from flux towers (e.g., using eddy covariance measurements of latent and sensible heat fluxes), though with known errors associated with scaling measurements up to a land surface model grid (Crow & Wood, 2002;Teuling et al., 2006). Also, the flux towers do not adequately cover the full range of combinations of climates and biomes. The newly available θ satellite remote sensing products, especially based on lower microwave frequencies (L-band) with reduced vegetation canopy opacity and greater moisture sensing depth (Entekhabi et al., 2010;Kerr et al., 2010), allow more direct observations of surface energy partitioning-θ relationships at spatial scales of land surface models. These remote sensing θ observations, when used alongside surface temperature metrics, have been shown to provide robust climate coupling assessments (Dong & Crow, 2018;Dong & Crow, 2019). Time series of satellite-derived θ during inter-storm periods  indicate that surface drying holds information about ET and hydrologic regimes comparable to those from models (Koster et al., 2009). Satellite θ combined with station-based observations of EF revealed that these EF-θ relationships hold temporally at a given location over seasonal timescales as well as spatially using mean fields (Short Gianotti et al., 2019). Climate and vegetation constraints on surface energy partitioning are apparent from these studies. However, fundamental characteristics of the surface energy partitioning-θ relationships outlined in Figure 1 remain to be quantified from observations. We do not use EF here as no direct, model-free estimations of ET exist at the land surface model grid scale. Instead, information from the diurnal evolution of land surface temperature is used here as it bears the signature of surface energy partitioning and can be observed using satellite remote sensing (Bateni & Entekhabi, 2012;Panwar et al., 2019). Specifically, increasing surface temperature diurnal amplitude (dT) in the water-limited regime can largely reflect the surface sensible heat flux becoming more efficient than latent heat flux at dissipating surface energy (Panwar et al., 2019). This manifests as less evaporative cooling from morning to afternoon at the land surface. dT is also expected to exhibit a lower magnitude of coupling with θ at high θ (under energy limitation). The θ-dependent division of behavior between high and low θ, if estimated robustly, is θ*. Positive coupling between dT and sensible heating and their negative coupling with EF and θ have been noted previously (Amano & Salvucci, 1999;Betts et al., 2014).

Surface Energy Partitioning-Soil Moisture Relationship: Potential Confounding Effects
The surface energetic partitioning-θ relationship is affected by other factors. The EF-θ relationship, for example, is a multidimensional function with dependence on factors, such as insolation (Rs), albedo (α), wind speed, and vapor pressure deficit (VPD), as confirmed by flux tower observations (Ford et al., 2014;Gu et al., 2006;Haghighi et al., 2018). The scatter in EF-θ plots is partly due to ignoring these dependencies and projecting EF only onto θ space. This can manifest as shifts in θ* and ∂EF ∂θ as shown by changes from Curves 1 to 2 in Figure 1. This framework was investigated for bare surfaces with findings that θ* can change significantly with temperature and humidity . The θ* threshold defines the θ value at which ET equals or sufficiently approaches potential ET. When considering widely used ET theory, this breakpoint depends on when the ratio of surface to atmospheric resistance approaches a small value (Peng et al., 2019). This condition can be seen when setting Penman potential ET and Penman-Monteith actual ET equal, for instance. This ratio is a function of all aforementioned factors: VPD, wind speed, Rs, α, and θ. While theoretical formulations like Penman-Monteith do not always match observations, it is conceivable that systematic changes in these variables to varying degrees can translate to different θ*. In the pursuit of detecting θ* and Stage II energetic partitioning sensitivity to θ, it is important to quantify and isolate the effects of the aforementioned potential confounding factors.

Study Goals
In this study, we use satellite observations of surface θ and dT with near daily sampling to evaluate the relationship between land surface energy partitioning and θ across a range of natural biomes in Africa. The θ fields are based on a low-Earth orbit satellite with a low-frequency microwave radiometer and global coverage. The dT range fields require resolving the diurnal cycle of the variable and hence require geostationary satellite records. For this, we use a multichannel optical to infrared sensor suite in geostationary orbit with full coverage of Africa (nadir of disc is centered on the Equator in Africa). We ask: (1) across climates and biomes, what are model-free estimates of (a) θ-dependent critical points between water and energy-limited regimes (θ*) and (b) surface energetic partitioning sensitivity to θ in the water-limited regime.
(2) Where are frequent transitions between hydrologic regimes occurring? (3) What is the temporal evolution of potential confounding factors during transitions from energy to water-limited regimes, and how might their dynamics confound surface energy partitioning-θ relationships at large spatial scales? Our objectives are thus to observe surface energy partitioning-θ relationships at land surface model scales, identify locations of hydrologic regime transitions, and assess whether variability of these relationships as found at tower sites can be consistently explained by other factors.

Soil Moisture Drydown Framework
The surface energetic partitioning-θ relationship has been shown to have strong climate/vegetation controls with spatial mean patterns and averaged temporal patterns (monthly) at a given location (Koster et al., 2009;Short Gianotti et al., 2019). In this study, we instead investigate these relationships during drying periods following rainfall where θ abruptly rises then decays until the next rain event, referred to henceforth as drydowns. A canonical drydown is shown in Figure 2a where the surface begins in an energy-limited state and decays below θ* into a water-limited state. We condition on these components of the time series for several reasons. (1) We can potentially isolate transitions from energy to water-limited regimes allowing more direct estimation of θ* and water-limited energetic partitioning sensitivity to θ. (2) These shorter time periods allow evaluation of potential confounding factor dynamics that might be otherwise suppressed with temporal aggregation (i.e., monthly and annually). (3) We are removing confounding effects of ET physics due to rainfall (McColl et al., 2019), frozen periods, and completely dry periods with no θ variations (dry season). (4) We take advantage of evaluating short-term (near daily) time evolution during drying segments to remove seasonal behavior and its many confounding factors associated with the time of year (vegetation phenology, mesoscale processes, etc.) rather than direct variable interactions (Tuttle & Salvucci, 2017).
We initially investigate how dT evolves in time along with surface energy fluxes during drying using eddy covariance measurements in a semiarid savanna ( Figure 2b) (Scholes et al., 2001). dT and EF evolve as expected during drying. However, G, a function of both soil thermal conductivity and soil temperature gradient (Castelli et al., 1999), tends to increase along with H (Santanello & Friedl, 2003). Despite soil drying leading to lower thermal conductivity, surface soil heating during drying likely results in a significant thermal gradient driving an increasing energy flux into the surface (greater G). G is commonly not discussed as a factor in climate coupling assessments using EF-soil moisture relationships despite being potentially a significant energy flux (Figure 2b) closely tied to the surface energy balance (Gentine et al., 2007). We argue here that dT is measuring a sensible heat flux signature that has a greater magnitude than G (Figure 2b) that might be subdued by G increases (less heating of atmosphere due to available energy used for soil thermal heating). Similarly, G was previously found to be less effective at dissipating incoming energy than H and LE (Bateni & Entekhabi, 2012). Ultimately, dT is a surface energetic partitioning metric and is expected to capture the surface energy balance's transition from energy to water-limited regimes during drying regardless of which energy flux component is dominant. In other words, relative changes in LE, H, and G are all components of the dT signal, not confounding factors.

Study Region
Africa is chosen as the study region due to collocation of θ from NASA's Soil Moisture Active Passive (SMAP) satellite (Entekhabi et al., 2010) and observations from the Spinning Enhanced Visible and Infrared Imager (SEVIRI) on-board European and EUMETSAT space agency's Meteosat Second Generation (MSG-2) geostationary satellite series . This synergy of datasets allows subweekly observations of energetic partitioning-θ relationships that can resolve post rainfall period dynamics. Africa additionally includes strong climate/vegetation cover gradients for evaluating this relationships' dependence on climate and vegetation.

Datasets 2.3.1. Soil Moisture
Surface θ (approximately top 5 cm) is estimated from NASA's SMAP satellite measurements of Earth's natural upwelling emission (Entekhabi et al., 2010). This includes 4 years (1 April 2015 to 31 March 2019) of the L1C brightness temperature product at an enhanced 9-km resolution on Equal-Area Scalable Earth 2 grid using Backus Gilbert optimal interpolation (Chaubell, Chan, et al., 2016;. The multitemporal dual channel algorithm approach is used to estimate θ which is not subject to ancillary vegetation information (Konings et al., 2016). Comparison statistics to in situ measurements are similar to those from the validated SMAP baseline θ product . See Feldman et al. (2018) methods for more details. The SMAP satellite samples a given location every 2 to 3 days. All further discussed datasets are linearly resampled to this grid and temporal resolution.

Meteorological Data
Meteorological observations are obtained from the MSG-2 geostationary satellite . It measures at multiple bands within the visible and infrared spectra with a maximum spatial sampling of 3 km at the subsatellite point (nadir), which is located at 0°longitude and 0°latitude. The daily downward surface solar radiation (Carrer et al., 2011;Ceamanos et al., 2014; and daily surface albedo (Carrer et al., 2018; products are used for confounding factor evaluations. Additionally, land surface temperatures (LST) (Freitas et al., 2010) at 6:00 AM and 1:30 PM are obtained by averaging the 15-min measurements at a 1-hr window around these times. All MSG-2 products are provided by the EUMETSAT Satellite Applications Facility on Land Surface Analysis . The 1:30 PM and 6:00 AM temperatures are then subtracted to obtain dT with only days with both values available retained. These times are representative of the temperature amplitude for this dataset in Africa (Holmes et al., 2015(Holmes et al., , 2013. Since the land surface temperature product is rigorously filtered for cloud contamination, days during rainfall are removed, replacing the need to otherwise filter out drydowns using relatively sparse precipitation datasets over Africa. This results in a large reduction in sample size (~60%). Therefore, for a given day with a SMAP overpass, any available dT observations since the previous SMAP overpass are averaged for a given pixel. This typically doubles available (θ and dT) pairs with the cost of introducing a slight mismatch between exact timing of θ and dT for a given SMAP overpass time. However, the temporal mismatch between dT and θ samples is small or 1-3 days due to SMAP θ sampling interval. Also, drydown detection (see below) intrinsically filters out cases of non-negligible storms between these 1-3-day SMAP θ sampling intervals that would introduce representativeness errors in the aggregated dT value. Alternatively, one could gap fill by linearly interpolating θ or combining θ products from other satellites (Kerr et al., 2010), though these methods are not attempted here.
Satellite LST is a directional variable since it corresponds to the radiometric temperature of the surface as seen by the sensor (Cao et al., 2019). A bias can therefore be introduced at high incidence angles where sensor view of vegetation, for example, becomes elongated. However, the LST measurement deviation from a reference view measurement (e.g., nadir view) is less than 1 to 2K for incidence angles of less than 40° ( Ermida et al., 2018), the case over the majority of the this study's domain. Furthermore, the use of a geostationary platform is an advantage because the viewing geometry is time invariant (i.e., no change between drydowns), removing confounding effects from changing viewing geometry. Additionally, the geostationary orbit captures the entire LST diurnal cycle providing the most reliable dT estimates.
Wind speed and VPD, also used for confounding factor evaluations, are obtained from MERRA-2 reanalysis as no known gridded observational dataset over this region and time period exist (GMAO, 2015(GMAO, , 2008. Daily VPD is computed using MERRA-2 surface pressure and 2-m air temperature and specific humidity fields at 1:30 PM local time (no temporal interpolation needed as fields are reported every half hour).

Ancillary Information
International Geosphere-Biosphere Programme land cover classifications are used to segment the results across ecohydroclimatic regions in Africa (Kim, 2013). Mean annual precipitation between 2015 and 2017 from GPM IMERG V5 is used to compare to spatial variations in metric estimates (Huffman, 2015). This ancillary information is not used in the estimation of thresholds and sensitivities.

Energetic Partitioning-θ Relationship Analysis 2.4.1. θ Drydown Identification
A similar procedure to previous studies is used to identify SMAP θ drydowns Feldman et al., 2018;McColl et al., 2017;Shellito et al., 2018). θ drydowns are identified for each 9 × 9-km 2 pixel as periods of consecutively decreasing θ for at least four SMAP overpasses (10-12 days) and continuing for no more than 20 days. The choice of these time periods is to filter out noisy θ periods and long drydowns that incorporate seasonal behavior. All drydowns must have at least a 0.01 m 3 m −3 θ increase from the previous overpass to begin the drydown, suggesting occurrence of rainfall. All aforementioned meteorological variables coinciding with these times are retained for the remainder of the analysis.
Seasonal trends in each variable during drydowns are removed with the same procedure as Feldman et al. (2018). Seasonal cycles are obtained by creating a mean climatology from averaging all 4 years of the data and smoothing the series with a 30-day moving average (Dong & Crow, 2018). A tutorial of this process is shown in Figure S1. This seasonality removal approach retains the absolute magnitude of all variables as opposed to converting into anomaly metrics. 2.4.2. θ* and m Estimation θ* and energetic partitioning sensitivity to θ in the water-limited regime (m) in the context of Figure 1 are estimated by fitting a piecewise linear model to the dT-θ relationship during drydowns. To ensure adequate (θ and dT) pairs during drying periods, the analysis is conducted at a 0.5°resolution containing

Water Resources Research
approximately 36 SMAP pixels. A representative example relationship is shown in Figure 3a. Visual inspections of numerous dT-θ relationships reveal approximately piecewise linear relationships with water limitation (with high dT-θ slope) at lower θ and energy limitation (with negligible dT-θ slope) at higher θ, as in Figure 3a. Our purposes here are to estimate an identifiable breakpoint of when dT increases at lower values (θ*) and the magnitude of dT increases per θ decrease of 0.1 m 3 /m 3 below this breakpoint (m). A simple piecewise linear model was deemed adequate to estimate these metrics, even if slight curvilinearities in the data exist. These simple models have been previously used in similar frameworks Schwingshackl et al., 2017). A shortcoming is that m as computed from dT may not have a one-to-one equivalence with the magnitudes of EF estimated m. This is due to physical differences between EF and dT where, for example, warm air advection or perhaps ground heat flux will differentially affect dT more than EF. The estimated θ* is considered a mean state value as it was established that this value is not necessarily constant . We evaluate in later sections how anomalous conditions can systematically shift θ*.
In some regions, a θ* and thus a transitional hydrologic regime may not be present due to θ varying only within a water or energy-limited regime. Therefore, two simpler models reflecting these respective conditions are also fit to the data in each pixel as shown in Figure 3a. Modeled dT is compared to the dT measurements with residual sum of squares (RSS) for each model. These models are equivalent to classes A to C in Akbar et al. (2018). Note that while θ* is only estimated for the transitional model case, m is estimated for both water limited and transitional model cases.
To select the "best" model and avoid overfitting the data in adding parameters progressively from energylimited (two parameters), to water-limited (three parameters), to transitional regime models (four parameters), the Akaike Information Criterion (AIC) (Akaike, 1974) is used to minimize RSS and number of model parameters (k): where N is number of sample (θ and dT) pairs. The lowest AIC determines which model is selected.
Fitting a model only to the phase space of dT and θ ignores the temporal evolution of dT during drydowns. This therefore leaves the possibility open that the phase space as in Figure 3a shows a trend of higher dT at lower θ that is only formulated from mean, seasonal conditions with actual negligible changes to dT during drying. To avoid this scenario and explicitly account for temporal evolution of each individual drydown in the phase space, two constraints on this piecewise linear model are included.
(1) The median of ΔdT/Δθ increments for θ < θ* must be negative. (2) The median of ΔdT/Δθ increments for θ < θ* must be less than the median of ΔdT/Δθ increments for θ > θ*. Note that since we are evaluating drying periods, Δθ is always negative and incorporates time with consecutive decreases. These conditions state that (1) dT must generally increase in time during drying in the water-limited regime and (2) that these dT increases are typically greater in the water-limited than energy-limited regime, respectively. Effectively, θ* is constrained to θ values that meet these conditions. If these two ΔdT/Δθ increment conditions are not met for the transitional hydrologic regime, the two simpler models (water or energy limitation only) are considered. These criteria removed 22% of θ* estimates where the simpler models were subsequently considered instead. Further, Condition (1) is also imposed for the water-limited only regime where the median of all ΔdT/Δθ increments must be negative. If this is not met, only the energy-limited regime is considered. Without imposing ancillary data, these constraints set a definition for θ* based on what is broadly expected of dT behavior during drying in water and energy-limited regimes. This ultimately ensures that when θ* is detected, the short-term temporal evolution of dT during drydowns and not seasonal relationships formulates the relationship shown in Figure 3a.
Model fitting is validated with a 10-fold cross-validation procedure (James et al., 2014) as used in previous hydrologic regime classification studies Schwingshackl et al., 2017). This provides both a test of robustness for the AIC method and bootstrapped confidence intervals on m and θ*. The approach is as follows. At each location, all observed (θ and dT) pairs are randomly divided into 10 equal-sized groups. Next, each model is fit to data using 9 of the 10 bins as training data to estimate their respective model parameters. Then, the fit models and parameters are compared against the data in the 10th group via a mean squared error (MSE) for each model. Analogously to AIC, this cross-validation procedure tests whether

10.1029/2019WR025874
Water Resources Research additional model parameters result in overfitting of the data which would manifest in high MSE in comparing model to validation data in the 10th group. Each group is used once as the validation data resulting in 10 iterations. Thirty Monte Carlo replicates of this process are obtained by randomly binning the (θ and dT) pairs into the 10 groups for each of the 30 iterations. Therefore, for each of the three models, an average of 300 MSE's is computed. The "best" model is that with the lowest average MSE. Note that this procedure was only used as a validation technique.

Potential Confounding Variable Impacts on θ* and m
We finally assess how potential confounding factors can influence θ* and m estimates (see section 1.3 discussion). We use two independent approaches to detect shifts in these metrics: (1) multiple regression and (2) conditioning techniques.
The multiple regression technique includes fitting a multiple linear regression model to dT during drydowns: where β 0 is the intercept, all other β are regression coefficients of respective dependent variable, and ε is residual error. Then, the variations of the potential confounding factors are subtracted from dT to obtain a dT time series that does not include the time variations of these factors (dT C ): Using the estimation procedure in 2.4.2., θ* and m are estimated from both (θ and dT) and (θ and dT C ) pairs and are compared. An example is shown in Figure 3b. dT C thus represents dT during drying periods with variability and coupling due to other factors removed. Only effects of θ principally remain. We remove all factors' variations together as in equation (5) as well as remove each individual factor's variations, where equation (5) is modified to subtract out only one factor's variations. Caveats exist in this framework including assumption of linear interactions, high covariation between regressors, and inability to account for autoregressive terms due to irreconcilable time gaps between drydowns (Gu et al., 2006). Finally, while we remove seasonality of the temporal evolution during each drydown, some seasonality may still exist within the drydown initial conditions. We repeat the analysis on the anomalies of each variable (not shown), though the interpretability of surface energy partitioning relationship with θ becomes nebulous. Overall, barring some caveats, this regression approach is attractive because it allows quantification of changes in θ* and m due to all variables together and each variable individually.
We also employ a conditioning technique that is independent of regression technique assumptions. In this conditioning approach, θ* and m are estimated from dT-θ relationships given only drydowns with a factor (1) significantly changing or (2) not changing. Specifically, in a given pixel, all drydowns are ranked from highest to lowest median change of a selected single factor (e.g., VPD). In the first case, θ* and m are estimated for (θ and dT) pairs during drydowns with above the 60th percentile of changes in the factor (increases in variable). For the second case, θ* and m are estimated for (θ and dT) pairs during drydowns with below the 40th percentile of changes (inactivity/decreases in variable). An example is shown for VPD in Figure 3c. Like in the multiple regression technique, the θ* and m from both methods are compared. This evaluates how drying periods with the presence of large increases in other factors (with significant magnitude of change as will be shown later in Figure 6) propagate into amplified/subdued water-limited couplings with dT (impacting m) and shifting a mean state θ*. This approach has its own shortcomings including foremost a reduced sample size from conditioning only on 40% of the data for the θ* and m estimates. It also only includes conditioning on single factors at a time ignoring interactions with other variables.
Ultimately, while both methods require caveats, they provide independent, data-driven tests of how systematic variations in other factors influence θ* and m estimates. Note that magnitudes of changes in θ* or m are not expected to be equivalent between methods. Also, some nonlocal factors (i.e., mesoscale advection of warm air) could influence dT along with the local factors considered here, but their impacts would require use of a model-based assessment which we do not attempt in this study.

Energetic Partitioning-Soil Moisture Relationship Observations
We discretize ecohydroclimatic regions by including pixels with the same International Geosphere-Biosphere Programme land cover classification (>75% homogeneity within each pixel) and collocation within similar climate regime as shown in Figure 4. Thus, each region roughly has the same vegetation type and climatic conditions. dT behavior during all drydowns within each regime is plotted together in Figure 5. General behavior is as expected from Figures 1 and 2, with dT negatively correlated with EF from previous studies (Dirmeyer et al., 2000;Koster et al., 2009;Short Gianotti et al., 2019). Biomes with predominantly herbaceous vegetation (e.g., grasslands, savannas, and shrublands) show stronger dT-θ coupling at lower soil moisture values than those with woody vegetation (e.g., woody savannas and broadleaf forests). The available energy is increasingly partitioned towards sensible heat flux with decreasing θ during the Stage II or water-limited evaporation regime. Slopes tend to be slightly lower at higher initial θ presumably due to initial well-watered conditions where most energy is partitioned into latent energy which keeps the surface temperature steadier throughout the day. Results from land surface schemes seem to generally reflect sensitivity of energy partitioning to θ at a large range of θ values (e.g., higher θ*) for drier regions as shown here (Dirmeyer et al., 2000). Further, increased surface energy partitioning sensitivity to θ (e.g., more negative slopes) as shown in these drier regions is currently captured in reanalysis (Schwingshackl et al., 2017). Interestingly, evergreen broadleaf forests (Congo Basin) show no dT sensitivity to θ even at low θ, despite some land surface models parameterizing a high sensitivity and low θ* (Dirmeyer et al., 2000). Further

Water Resources Research
investigations reveal that low θ satellite observations in the Congo Basin are likely the result of noisy retrievals. The behavior in evergreen broadleaf forests should therefore be treated with caution as low θ (water limitation) may not have been actually observed in this region during the study period.
The behavior in Figure 5 is not due to seasonal mean states. Evaluating individual drydowns on this phase space shows that the dT-θ relationships are formulated from the time evolution of dT increases during the 10-20 days of drying. Note that progression in time is from high to low θ. Again, seasonal trends are removed from these relationships indicating that the behavior is not due to phase differences associated with seasonal cycles (Williams & Torn, 2015).

Potential Effects of Confounding Factors
Directly addressing how other relevant factors evolve following a rain event (see section 1.3 and Figure 1), Figure 6 shows short-term behavior of Rs, α, wind speed, and VPD during drying with progression in time from high to low θ. Observed Rs in Figure 6a shows large increases during drying following rain events. Tests with SEVIRI cloud reflectance and MERRA-2 reanalysis cloud cover information revealed that this is due to progressively decreasing cloud cover after precipitation ceases. SMAP's low spatial resolution is likely sampling mesoscale events rather than smaller scale, convective cells. The temporal persistence of cloud cover is generally greater for these mesoscale events as observed in the Sahel (Mathon & Laurent, 2001). Some studies noted the effects of solar radiation on temperature (Teuling et al., 2006), while others attempt to remove its effects (Miralles et al., 2012). Regardless, one can hypothesize that the Rs increases could amplify partitioning into sensible heating (and ground heat flux) during drying. Note, the Rs decrease at high initial θ in open shrublands is likely an artifact of averaging on the phase space as all individual Rs trajectories increase under all θ (not shown).
Observed α in Figure 6b, at least in the driest biomes, appears to show a~0.03 increase during drying periods consistent with previous in situ investigations in the US Great Plains (Amano & Salvucci, 1999). α is expected to increase from drying soil and desiccating vegetation. Wind speed exhibits a tendency to increase during drying under drier surfaces in Figure 6c. Bare soil modeling results showed that increasing wind speed could strongly shift θ* to wetter values . Finally, VPD increases significantly due to both increasing temperature and decreasing humidity as shown in Figure 6d. These effects are likely coupled to amplifying sensible heating (as well as ground heat flux) during surface drying. Overall, the Figure 6 behavior indicates that factors other than θ may systematically explain the dT-θ trends and scatter shown in Figure 5.

θ*, m, and Hydrologic Regime Estimation
Estimated θ* from observations (θ and dT) are shown in Figure 7a from those pixels robustly identified as having transitional regimes in Figure 3a. Areas without an identifiable θ* are characterized by either a persistent water-limited or energy-limited evaporation regimes. Spatial patterns of θ* are coherent and show a tendency to decrease with sand fraction (ρ = −0.43; p < 0.01) similarly to Akbar et al. (2018). While the time it takes to reach θ* may be quicker in coarse soils during drying (Aminzadeh & Or, 2014), moisture extraction via ET can occur at lower θ due to less matric suction than in finer soils. This is consistent with findings of lower θ thresholds on plant water uptake for coarser soils (Feldman et al., 2018). θ* and mean annual precipitation are weakly, positively correlated (ρ = 0.26; p < 0.01) with the sign of dependence consistent with previous studies Schwingshackl et al., 2017). Reanalysis and model estimates show a much stronger climate dependence with higher θ* generally closer to the equator, especially in Africa (Schwingshackl et al., 2017). In contrast, these estimates here interestingly indicate that all biomes have similar θ* distributions (see Table 1), despite a wide range of mean annual precipitation values. Pairwise Kolmogorov-Smirnov tests indicate that grasslands and woody savannas have θ* spatial distributions typically significantly different from the other biomes' distributions (p < 0.05) (Massey, 1951). However, the

10.1029/2019WR025874
Water Resources Research θ* medians and standard deviations of each biome (as well as their probability density functions; not shown) are all broadly similar (see Table 1). Therefore, θ* itself is somewhat independent of the mean moisture state. To a first order, one would expect lower θ* in humid regions that have generally greater vegetation and less bare soil cover. Root systems are typically accessing deeper soil water pools than the soil surface interface (as represented by SMAP θ) where bare soil evaporation occurs. These deeper water pools have slower changing dynamics and may allow transpiration to continue at its potential even after evaporation from the soil surface interface becomes water limited. However, a study contrasting semiarid savanna and grassland appear to show θ* of near 0.15 m 3 /m 3 for both biomes, supporting some independence of vegetation cover on θ* (Baldocchi et al., 2004). Thus, while a spatially coherent θ* pattern (evidence for some spatial autocorrelation in Figure 7a instead of complete randomness) emerges here, coarse climate considerations do not explain it.
The slope, m, is greatest in transitional regions between humid and arid climates as shown in Figure 7b. This appears to be captured in some reanalysis outputs (Schwingshackl et al., 2017) but is inconsistent with some land surface schemes that are parameterized with greater m in wetter, forested regions than for drier, less vegetated regions (Dirmeyer et al., 2000). Generally, higher m is expected in drier regions with herbaceous vegetation cover because of (1) vegetation with more responsiveness to individual rainfall events (and less stomatal control) and (2) lower vegetation cover. In tropical dry regions, vegetation hydraulic systems largely respond to individual rainfall events (Feldman et al., 2018). Thus, a transpiration response is expected given a perturbation of θ. In contrast, woody vegetation in more humid environments may be exhibiting greater stomatal control under water limitation (Konings & Gentine, 2017;Lin et al., 2015), accessing deeper water stores not observed by the satellite and/or showing generally less responsiveness to changes in surface soil water availability. These scenarios could reduce m. Also, the drier regions generally have a significant contribution from bare soil due to patchy or light vegetation cover. Bare soil is always responding to changes in θ under water limitation. These factors could in part explain the spatial pattern in Figure 7b. Note that for simplicity, m is reported as the increase in dT given a unit soil moisture decrease of 0.1 m 3 /m 3 during a drydown. The patterns of m are interpretable as a proxy for coupling between the surface water and energy balances without impositions of land-atmosphere model assumptions. Given physical differences between dT and EF, further work is required to evaluate whether m magnitudes uniquely relate to ∂EF ∂θ magnitudes. The 10-fold cross-validation procedure indicates that the θ* and m estimates with AIC are robust ( Figures S2  and S3). The spatial median θ* bootstrapped standard deviation is 0.004 m 3 /m 3 . Additionally, AIC estimated θ* is within 0.01 m 3 /m 3 of the mean θ* from the 10-fold cross validation replicates 86% of the time. Similarly, the median m standard deviation is 0.3 K. AIC estimated m is within 1 K of the mean m from the 10-fold cross validation replicates 83% of the time. A slight bias between AIC and cross-validation m estimation exists near the Kalahari Desert in Southern Africa. The 10-fold cross-validation procedure selected the transitional model regime 59% of the time, while the AIC method selected it 57% of the time. Thus, AIC is at least marginally conservative in selecting transitional regimes, the most complex model. Therefore, the AIC approach is deemed sufficient for use in this study, especially for the purposes of evaluating θ* and m parameters rather than details of the hydrologic regime classification. This justifies the sole use of the AIC Note. All values are reported in spatial medians with the spatial standard deviation in parenthesis. m is reported as dT change [K] per θ decrease of 0.1 m 3 /m approach in quantifying effects of other factors on θ* and m. Model RSS from AIC are also relatively uniform throughout the study region (see Figure S4).
One utility of observation-based θ* estimates is to assess where the local environment varies across this threshold and transitions from uncoupled to coupled water and energy balance regimes. We estimate the fraction of time in the water-limited regime or percentage of time the observed θ is below θ* (see Figure 7 c). The dominant regime classification is shown in the inset of this figure. The dominant regimes observed broadly corroborate previous findings from land surface models (Koster et al., 2009;Schwingshackl et al., 2017). However, regime dominance only reflects the spatiotemporal scales of the datasets. Any given location from rainforest to desert will experience energy and water-limited regimes at some point given sufficient years of data at a frequent sampling interval (i.e., hourly). Instead, the time spent in a given regime holds more information as a continuous metric with less dependence on the dataset's spatiotemporal scales (scales may still impact this metric's robustness). Note that this uses the full θ time series to avoid biases from inconsistent drydown sampling between pixels. A spatial gradient is observed with drier regions generally spending more time in the water-limited regime as expected (Table 1). It is worth noting that biomes in the Sahel tend to exhibit collocation of frequent transitioning between hydrologic regimes and large m (compare Figures 7b and 7c). The Sahelian grasslands spend on average 84% of the time in the water-limited regime. However, when considering only drydowns that occur primarily in the wet season, 66% of the time is spent in the water-limited regime reflecting rainfall regimes that promote frequent transitions between Stage I and II ET. This indicates that ET may be often close to potential ET following rainfall when in the water-limited regime. This suggests a potential for land surface control on the lower atmosphere as there is both strong energetic partitioning sensitivity to θ and likely sufficiently high ET (Koster et al., 2004). Many previous studies labeled this region as a hotspot showing strong soil moisture-climate coupling (Dirmeyer, 2011;Shinoda & Yamaguchi, 2003) and probable soil moisture-precipitation feedbacks (Koster et al., 2004). In contrast, woody savannas south of the Congo Basin also exhibit frequent transitioning between energy and water limitation, but m is lower on average. Spatial variability of both metrics is greater south of the equator with potentially isolated coupling hotspots.
θ retrieval uncertainty is higher in dense broadleaf forests and woody savannas as shown in Table 1. Here, lower mean difference and standard deviation of each brightness temperature (TB) polarization (wave orientations used to estimate θ) indicate lower signal to noise ratio. Both of these metrics need to be sufficiently high to robustly retrieve θ. Specifically, high mean difference between TB polarizations allows partitioning between soil and vegetation surface emission, and high TB standard deviation (larger dynamic range) results in less noise in θ estimation. Land surface temperatures used to derive dT from SEVIRI have biases of less than 1 K throughout vegetated biomes (Freitas et al., 2010;Göttsche et al., 2016;Trigo et al., 2008). Instead, variations in humidity and aerosols optically thicken the atmosphere and result in significant uncertainty (Freitas et al., 2010;Göttsche et al., 2016). Also, the retrieved temperature characterizes the surface skin which may represent the canopy crown in forests as opposed to the soil surface. Soil-vegetation temperature disequilibrium can indeed occur during the day. Ultimately, broadleaf forest and less so woody savanna results should be considered with caution.

Effects of Other Factors on θ* and m Estimates
Regression approach results are discussed first, quantifying the relative effects of VPD, α, Rs, and wind speed (both individually and altogether) on m and θ*. The conditioning approach is used to corroborate individual factor effects on m and θ* as determined in the regression approach. The metrics are interpreted as the shift in m and θ* due to VPD, wind speed, α, and/or Rs variations.
As shown in Figure 8a, the variations of all factors together as determined from the regression approach tended to reduce θ* in general with a median 0.01 m 3 /m 3 decrease. VPD variations have the greatest influence on this θ* decrease as shown in Figures 8b and 8d (see also Table S1). This VPD effect occurred in semiarid regions such as the Sahel. Elsewhere, the relative effects of factors are random showing that wind speed, α, and Rs show no consistent, significant impact on θ*.
Analogously, Figure 9 shows these factors' influence on m. Their variations tend to increase dT sensitivity to θ by about 2 K per 0.1 m 3 /m 3 of soil drying compared to not including their variations. In this case, the effects of all variables on m show more coherent spatial patterns than for θ*. VPD again dominates contributions to the increase in m in subhumid regions. However, solar radiation and α do show coherent, dominant impacts in the arid regions in Southern Africa and boundaries of the Sahara Desert.
To corroborate the regression approach findings of shifts in θ* and m due to individual factors, we compare to conditioning approach results. We evaluate whether an individual factor's variations changed θ* or m (1) significantly and (2) in the same direction for both independent methods (regression and conditioning approaches). This is evaluated for all pixels in Africa as well as pixels in each ecohydroclimatic region. Evergreen broadleaf forests are excluded due to predominant detection of energy limitation. For each biome and potential confounding factor pair, a Mann-Whitney U test (Mann & Whitney, 1947) is performed between spatial θ* distributions with the factor's variations included and removed. This is repeated for m spatial distributions. This nonparametric test is performed as Jarque-Bera normality tests (Jarque & Bera, 1980) reveal all θ* and m spatial distributions are non-Gaussian. Results are shown in Table S1 and summarized here. Both regression and conditioning approaches independently show that θ* decreased (~0.01 m 3 / m 3 ) in the Sahelian grasslands when VPD increased during drydowns. In other words, increasing atmospheric demand during the drydown shifts the water-limited regime transition to a lower θ. However, a robust shift in θ* is not detected in other biomes due to any other factor. m increases due to VPD are on the order of 1 K and are more widespread across generally drier biomes. Rs increases also are associated with the same magnitude of increase m but only in shrublands of Southern Africa. No consistent, significant effects of α and wind speed are found. This is likely due to minimal α dynamics and inconsistent wind speed influences on energy partitioning during drying.

Water Resources Research
These results suggest that (1) there is a spatial dependence on where these factors tend to impact surface energy partitioning as coupled to θ on weather timescales. (2) VPD and Rs are the only variables that are robustly identified as confounding the surface energy partitioning-θ relationship.
(3) The actual impacts of these factors on θ* and m are typically minor. Ultimately, despite clear variations in these variables during transitions from energy to water limitation as in Figure 6, moisture availability remains to be the dominant factor in partitioning surface energy. Specifically, VPD appears to influence θ* and m more than Rs, α, and wind speed, but these VPD variations are highly coupled to θ. VPD, the difference between saturation vapor pressure and actual vapor pressure, increases with warmer and drier air conditions. Soil moisture contributes to both of these factors: decreasing latent heat flux favors sensible and ground heat flux for dissipating available energy. This warms the overlying air and soil and reduces ET, decreasing humidity. Considering the multiple regression coefficients multiplied by their respective factor's standard deviation, VPD and θ variations are coupled to the largest changes in dT (~2 K) with typically little influence from the other factors (~0.1 K). External forcings, such as Rs and wind speed, therefore do not appear to have widespread, consistent influences on surface energy partitioning, despite some coupling to θ as shown in Figure 6. Further work is required to separate effects of surface and atmospheric moisture availability and their relative contributions to surface energy partitioning. We caution that these results do not suggest that time evolution of external forcings can be ignored in coupling analyses, for example. As seen in Figures 8 and 9, there are cases where confounding effects on surface energy partitioning are potentially large with θ* and m changing in time at a given location. Further work should also investigate whether these results are sensitive to spatiotemporal scales. It is suspected that anomalous forcings become more apparent on energetic partitioning-θ relationships at finer spatiotemporal scales. Indeed, pore-scale bare soil modeling studies showed vast variations in θ * especially with temperature and wind speed .

Water Resources Research
It is interesting that VPD anomalous increases tend to decrease θ*. One might expect the opposite where water limitation would be reached at a wetter condition during drying. Assuming similar relationships between θ* and potential confounding variables in both ET and dT space, θ* therefore relies on the interplay between actual and potential ET. Specifically, θ* should be reached during drying when actual ET initially decreases below potential ET. All potential confounding factors generally increase during drying (see Figure 6), and all should contribute to an increasing potential ET. Thus, if actual ET hypothetically remained constant or decreased while potential ET increased during drying, θ* would increase. However, a decrease in θ* suggests that VPD anomalously increasing is due to actual ET continuing near its potential under drier surface conditions. In other words, VPD increasing might cause both actual and potential ET to increase together during drying such that θ* decreases.
For m, we postulate that the increases in m due to anomalous increases in VPD and Rs are likely due to amplified sensible heating. In comparing Figures 5 and 6, increasing VPD and insolation both increase along with dT under water limitation which would suggest that they could amplify dT increases during drying.
It is important to note that this analysis is conducted on the time evolution of dT changes with drying soil. This analysis (Figures 8 and 9 and Tables S1) quantifies how short-term (1-3 day) changes, rather than mean state, of other factors influence dT-θ relationships. Conditioning dT-θ relationships on the mean state of net radiation (high or low), for example, would only suggest how dT-θ relationships are different between seasons. Since net radiation is coupled to nearly every other factor between land surface and boundary layer (Betts, 2004), this would say little about net radiation's direct influence on surface energy partitioning. Thus, conducting the analysis on the variables' short-term time evolution should be relatively isolated from mean climate spurious relationships.

Discussion
With VPD and Rs increases during drying on subweekly scales as in Figure 6, it is plausible that positive feedbacks between these factors can be detrimental to vegetation health. At low θ, vegetation water content at large scales appears to decrease following rainfall along with θ (Feldman et al., 2018). Thus, the sequence of increasing incoming shortwave radiation, increasing atmospheric demand, increasing sensible heating, and/or drying soil may exacerbate vegetation function over days to weeks given insufficient rainfall. This has implications for θ limited transpiration which also has a positively feedback that further limits rainfall (Green et al., 2017). There is thus cause to evaluate vegetation parameter time evolution within this multivariable framework in addition to its influence on land-atmosphere coupling (Zscheischler et al., 2015).
With advances in raw and assimilated remote sensing soil moisture, soil moisture-temperature coupling estimates are becoming more robust (Dong & Crow, 2019, replacing previously used precipitationtemperature metrics. The dT-θ relationships found here present caution to soil moisture-temperature coupling metrics. Above θ*, θ controls on afternoon surface temperature, like dT, are subdued (not shown). Thus, coupling metrics may need to carefully account for differences in θ controls on temperature considering water and energy-limited regimes instead of all regimes together.
The "triangle method" has been used to estimate θ from satellite optical vegetation information and surface temperature (Carlson et al., 1994). The changing sensitivity of dT to θ in in the water-limited regime due to changes in insolation and VPD can confound this θ estimation technique.
Implications for land-atmosphere feedbacks are only for the terrestrial segment, with soil moisture influencing the surface latent and sensible heat fluxes (Dirmeyer, 2011). The time evolution during 10-20-day drying periods following rainfall shows a widespread tightly coupled system of soil water loss and increasing atmospheric demand, sensible heat, and incoming radiation. Despite this coupling, θ* and m appear to sustain only minor systematic shifts at spatiotemporal scales of land surface models. Further, under coarse scale considerations, it is noteworthy that the θ* distribution is relatively narrow with similar values in each ecohydroclimatic regime. Therefore, moderate moisture regimes varying just under θ* with high m are expected to show strong land influence on the lower atmosphere. This appears to occur broadly across the Sahel region. In contrast, minimal land control on the lower atmosphere is expected in humid regions that have θ variability above θ* (e.g., energy limited).

Conclusions
We use two remote sensing datasets to examine key aspects of the relationship between surface water and energy balances over land surfaces. We use satellite-based dT range from a sensor in geostationary orbit and soil moisture (θ) from a low-frequency sensor in low-Earth orbit with shared coverage of diverse climates and biomes across the African continent. These two datasets are independent of any models or parameterizations linking the water and energy balances. They are used to quantify key characteristics of this coupling using the time evolution of dT and θ observations during drying periods: (1) the soil moisture threshold (θ*) that marks the transition from energy-limited to water-limited evaporation regimes and (2) the strength of the coupling between dT and θ in the water-limited regime. On average, ecohydroclimatic regions with drier conditions and herbaceous vegetation cover tend to exhibit stronger coupling of water and energy balances during water-limited ET regimes, while the θ* values and patterns are similar under nearly all of the diverse ecohydroclimate regions of Africa. The regions with the highest land surface flux partitioning sensitivity to θ are interestingly collocated with regions that frequently transition between energy and water limitation. This specifically occurs in the Sahelian grasslands, a region previously receiving much attention as a land-atmosphere feedback hotspot.
The results shown here are purposely free of land surface model and parameterization influences. Estimates of these key characteristics based on land surface model outputs may only reflect explicit and/or implicit model formulations. These estimates provide observations of θ control on energetic partitioning (as commonly visualized in the EF-θ relationship) at the spatial scales of land surface models, noting that dT is used here as proxy for sensible heat flux.
We also assess the possible influence of other environmental variables on dT during drydowns. Vapor pressure deficit, α, wind speed, and insolation generally significantly increase across Africa during drying periods. In order to test how these variations may confound surface energy partitioning-θ relationships, we use regression and conditioning approaches to test whether effects of each factor shift the value of our study metrics. Ultimately, only vapor pressure deficit and less so insolation tend to have any influence on these relationships. Specifically, anomalous increases in both variables tend to increase dT by approximately 1 K per 0.1 (m 3 m −3 ) of soil drying. Vapor pressure deficit increases tend to reduce θ* by 0.01 (m 3 m −3 ) though only in Sahelian grasslands. Vapor pressure deficit and insolation thus show greater potential as confounding factors than α and wind speed. The median 0.01 m 3 m −3 and 1 K shifts in θ* and m are, however, arguably small, though isolated cases show larger changes. Additionally, vapor pressure deficit is strongly coupled to θ which confounds effects of direct control of humidity deficits on surface energy partitioning as separated from soil water deficits. Ultimately, both surface and atmospheric moisture availability together strongly constrain surface energy partitioning despite significant time evolution of other surface and atmospheric variables during these drying periods. We suspect some dependence of results on spatiotemporal scales. More advanced sensors with finer spatiotemporal resolutions will likely be able to detect stronger, systematic effects of confounding factors on energetic partitioning-θ relationships.