Local and General Patterns of Terrestrial Water‐Carbon Coupling

Terrestrial carbon uptake and water availability have coupled feedbacks; specifically water uptake for plant growth and soil drying via transpiration. While we might expect this coupling over time at arid sites, climatic water availability also widely covaries geographically with biomass variables that control photosynthetic rates. Using eddy covariance data globally, we find convex, positively‐covarying relations between carbon uptake and a turbulent flux metric controlled by land surface moisture (r = 0.73 monthly across sites) at the site level. We estimate a general, empirical relationship based on site‐wise water‐carbon dynamics. Most sites, and the general relationship, show strong power‐law dependence, implicating the role of sub‐seasonal land‐cover dynamics. We also find that long‐term mean carbon/water states follow a similar convex relationship to the site‐specific temporal dynamics. We discuss opportunities and caveats for space‐for‐time frameworks of carbon/water feedback processes globally.


Introduction
The terrestrial water and carbon cycles are coupled from leaf-to global-scale.This means that there are feedbacks in both directions-from water state variables to carbon fluxes, and from carbon state variables to water fluxes.Biosphere-hydrosphere feedbacks drive the heat fluxes that control surface climate variables, and they integrate complexity across timescales, across spatial heterogeneity, and across disciplinary boundaries (Abbott et al., 2019;Green et al., 2017;Humphrey et al., 2018;Piao et al., 2020;Yang et al., 2019).
In the water-to-carbon direction, water is a primary control on the land carbon sink (Humphrey et al., 2021;Trugman et al., 2018).Water availability (soil moisture, surface water resources, accumulated precipitation, etc.) shapes biome type, stability, succession and landscape carrying capacity.On sub-annual time-scales, plants grow and rehydrate in response to discrete precipitation events (Feldman et al., 2019;Konings et al., 2021), and suppress shoot and leaf growth in response to water stress, particularly in arid regions/conditions.This can be framed as a "growth-acclimation," growth-inhibition, or "pulse-response" mechanism, depending on whether viewing vegetation as in dynamic equilibrium or permanently-transient (Chaves et al., 2002;Claeys & Inzé, 2013;Reynolds et al., 2004).Water is not the sole resource determining landscape vegetation; both water and biomass variables co-evolve in relation to each other and to independent drivers.
In the carbon-to-water direction, vegetation redistributes water between land and atmosphere, creates a biotic surface for the surface energy budget, and changes the broad flows of the hydrologic cycle (Lemordant et al., 2018).Plants transport water from soils to the lower atmosphere through transpiration and (less-directly) subsurface water redistribution (Caldwell et al., 1998;Stoy et al., 2019).Further downstream, these processes alter surface hydrology, atmospheric humidity, and precipitation regimes-all of which feed back on vegetation dynamics (Jasechko, 2018;Mankin et al., 2019;Yu et al., 2017).
Myriad other biochemical, physical, and geologic processes link the water and carbon cycles as well, all of which create (a) substantial process complexity and (b) difficulty in properly modeling subcomponents of the coupled water-carbon interactions, even at the point-scale.And while it is well-known that individual plants within-and across-species display large heterogeneity in their water-carbon linkages, it is not well-understood if aggregation across landscapes increases or decreases the signals of these linkages (Harrison et al., 2020;Kattge et al., 2011).

Plant-Driven Controls on Evaporative Fluxes
Water-carbon relationships are often viewed through the lens of water-use efficiency at the leaf-, canopy-, or landscape-averaged-scale.For example, the bulk canopy-scale water-use efficiency, ω c , can be defined as ω c = GPP c /T c , where GPP c (Gross Primary Productivity) is canopy-averaged carbon uptake (mass per canopyarea per time) and T c is canopy-averaged transpiration (volume per area per time, defined over only the vegetated portion of a landscape).The subscript c denotes canopy-averaged quantities.To scale GPP to the landscape level -GPP l , including non-canopy bare soils-we multiply by the vegetation cover fraction f v (canopy area divided by landscape area): Water-use efficiency and transpiration rates fluctuate at diurnal timescales and faster with changes in leaf-level energy balances, hormonal and hydraulic adjustments to stomatal conductance, turbulence generation, and the moisture of plant, soil, and atmospheric pools.Flux tower data suggests that ω c typically decreases in wetter soil conditions (Nie et al., 2021); both landscape transpiration and "T over ET" ratios are strong functions of f v (Nelson et al., 2020;Stoy et al., 2014;Wei et al., 2017).

Water Constraints on Photosynthesis
Another common lens is photosynthetic limitation.A simple model at the scale of a chloroplast or small leaf area gives the carbon uptake rate as where A max is the maximum carbon uptake rate, θ is a water-availability function (from 0 to 1), and λ represents all other limitation processes (also 0 to 1).Total landscape GPP is GPP l is thus a function of (a) the number of landscape chloroplasts n chl -to first-order a function of leaf area index (LAI) or f v , (b) photosynthetic capacity (subject to adaptive dynamics), (c) the amount of available water, and (d) other limitations including light, temperature, and nutrients.
To constrain the global carbon cycle, we are typically interested in the landscape-scale GPP l .We are reminded from the above lenses that (a) GPP l scales with water-supply (T l , θ), energy supply (T l , λ), physiological variables (ω c , A max ), and landcover structural variables ( f v , n chl ).(b) The structural variables are clearly correlated with the water and energy variables at subseasonal timescales, as are water and energy variables with each other.(c) The non-structural variables all covary at sub-daily scales as well.(d) The coupled dynamic relationships between these variables are represented differently across different models of Earth's climate system, ecosystem dynamics and demography, and land-atmosphere interactions (e.g., Hay-Chapman & Dirmeyer, 2023;Ukkola et al., 2018;S. Zhang et al., 2019).
In this study we will analyze observations of coupling between water and carbon variables over time at multiple sites and time-averaged across sites.These analyses will be observationally-driven, rather than following a specific lens above (or similar), to provide benchmarking relationships for conceptualization and models.We will use data from networks of eddy covariance towers, with normalized evaporative fluxes to reduce the effect of variable energy-balance regimes.The analyses will be at monthly scale to capture dominant subseasonal dynamics, averaging over higher-frequency (sub-daily) fluctuations.We will show a relatively unexplored bivariate coupling of water and carbon variables, the monthly co-evolution of landscape carbon uptake (Gross Primary Productivity-GPP l ) and a turbulent surface energy flux partitioning variable (the Evaporative Fraction-EF).
Local relationships are strong at subseasonal scales, across-site patterns of behavior are broadly similar, and simple generalized empirical relationships between these variables have high predictive power.We begin with discussion of EF as a normalized measure of water availability across climates.

Evaporative Fraction
EF is a metric of turbulent energy flux partitioning, EF = LE/(LE + H), where LE is latent heat flux and H is sensible heat flux.This is related to the Bowen ratio B as EF = 1/(1 + B).Unlike the Bowen ratio, EF is bounded.
Most directly, EF is the fraction of available surface energy that goes to evaporating water versus warming the atmosphere (EF = 0 implying only sensible heating, EF = 1 implying only latent heating).
EF is a useful metric for inter-site comparison, as the denominator controls to a large extent for available energy (e.g., radiation) and variable turbulent energy states across weather, climate, and seasonal differences (Haghighi et al., 2018;Short Gianotti et al., 2019).This normalization also keeps EF relatively stable to changes in energy availability (e.g., solar radiation) and wind speed over time, instead responding primarily to the slowly-evolving availability of water (Dirmeyer et al., 2000;Gentine et al., 2007).This is particularly important when comparing latent heat flux observations, due to the non-linear response of saturated humidity to temperature.
Classic agronomy studies estimate carbon assimilation through the empirical relationship GPP c = kET c /E p , where GPP c is canopy GPP, ET c is plant or canopy-scale evapotranspiration, E p is potential evaporation, and k is a species-specific constant (de Wit, 1958).Eddy flux-based observations suggest that the most appropriate landscape E p is the available turbulent energy E p ∼ LE + H, generally scaled by 0.8 (Maes et al., 2018;Milly & Dunne, 2016).This makes the GPP model GPP = k • EF, suggesting a direct GPP-EF relationship in time, or a species-varying relationship across space.
While there have been major advances in our understanding of plant carbon assimilation since these models were developed, there are reasons to revisit the GPP-EF relationship.First, there is evidence of strong temporal GPP-EF correlation at individual eddy covariance sites (Williams & Torn, 2015).Second, there is evidence of acrossspace covariance between GPP and various water metrics that are reinforced to some degree by temporal relationships (Biederman et al., 2016;Huxman et al., 2004;Short Gianotti et al., 2019).
While we will use eddy covariance estimates of latent and sensible heat partitioning in this study, EF is also estimateable from surface meteorological observations due to emergent land-atmosphere relationships with relative humidity (McColl & Rigden, 2020;Salvucci & Gentine, 2013).EF and humidity (hence vapor pressure deficit-VPD) are tightly coupled themselves, although VPD is more directly subject to changing atmospheric variables.For the purposes of this study, EF's observational independence from carbon cycle variables, its relative independence from variable weather conditions, its across-site inter-comparability, and its monotonicity with water availability lead us to select it over other possible water metrics such as unnormalized λE, VPD, water-use efficiency metrics, or estimated surface conductance.

Time Scales
We focus in this study on monthly data, as the most direct water/carbon feedbacks have large dynamics on this timescale.Using available data, we have little capacity to resolve successionary (or even inter-generationally adaptive) processes.Subseasonal data allow us to capture the major variability on seasonal timescales, and they integrate the primary feedback mechanisms (land cover dynamics, transpiration, leaf growth and growth acclimation, senescence, hydration, and surface conductance) on timescales of water availability changes (e.g., root zone soil moisture timescales).

Eddy Covariance Data
It is necessary to use direct observational data to investigate the relationship between GPP and EF, as many estimates of both variables (e.g., from remote sensing) rely on assumed relationships between the water and carbon cycle.Data for this study come from numerous tower-based eddy covariance flux measurement sites (Baldocchi, 2014).We draw from three publicly available data sets representing networks of independentlyoperated sites: (a) the FLUXNET 2015 data set, (b) the Ameriflux FLUXNET SUBSET data set, and (c) the Ameriflux BASE data product (Pastorello et al., 2020).From each data set, we select all sites with records available under the CC-by-4.0 license.For sites in multiple data sets, we select the record with the longest temporal span (321 unique sites).
For each site, we use monthly reported values of GPP using the nighttime respiration partitioning method (GPP_NT_VUT_REF), latent heat flux (LE_F_MDS), and sensible heat flux (H_F_MDS).If monthly data is not available for a given record (Ameriflux BASE product), we aggregate hourly values to monthly; if hourly values are not available we use half-hourly.Only days with at least 22 of 24 hourly or 44 of 48 half-hourly observations are used.Missing hourly or half-hourly data are replaced with the monthly mean for that time-of-day before aggregating to monthly to avoid diurnal sampling bias.
EF is calculated at monthly time-scales as the ratio of mean monthly LE to total available turbulent energy (LE + H).

IGBP Exclusion
Since we are looking at the effect of water availability on carbon fluxes, we exclude those sites denoted as primarily Snow and Ice, Urban, Wetlands, or Water, as well as any sites described as including wetlands in their main vegetation IGBP description.

Sufficient Data Exclusion
From these records, we calculate the maximum monthly GPP for each site.To reduce noise from non-growingseason data, we exclude all individual months at a site with less than 25% of the site's local maximum GPP.We then use only those sites with at least 20 months of useable GPP and EF data.
Following these exclusions, we have 183 sites with at least 20 useable monthly data points.For the list of all sites, their attributes, data period used, and appropriate background material, see Table S1 in Supporting Information S1.

Caveats
Eddy covariance tower data are the current benchmark data sets for stand-and landscape-scale fluxes across land covers and climates.There are noted errors in energy balance closure across most sites-typically a biased underestimation of turbulent fluxes-hypothetically related to sub-mesoscale transport, neglected energy pools, instrumentation, and violations of assumptions about local turbulence and homogeneity (Mauder et al., 2020).To the extent that the errors in sensible and latent heat fluxes are correlated-which is commonly assumed-use of the evaporative fraction removes these biases.Additionally-since the differences in site instrumentation, biome, and climate can affect the magnitude of energy closure biases-removing sources of bias at each site (normalizing by LE + H) helps to reduce across-site (spatial) biases which occur in GPP-LE relationships.

Curve Fits
Because the coupled feedbacks between GPP and EF are so mechanistically-entangled, we fit empirical parameterized curves of the form The power-law relationship is not based on any theory or assumptions, but is selected for its simplicity and empirical flexibility.
The parameters of bivariate relationships can be highly sensitive to which variable is treated as "independent."Because in this case, it is clear that neither GPP = f(EF) nor EF = f(GPP) is physically defensible due to the bidirectional complexity of feedbacks, we instead use a common regression model which treats both variables as equally subject to (observational-or model structural-) errors-total least-squares (TLS).TLS is designed to minimize biases (including omitted variable biases) that can occur from using ordinary least squares methods while assuming a unidirectional independence/dependence relationship, particularly for non-linear covariances.This fitting treats GPP and EF as fully coupled with error structure locally orthogonal to the power law curve fit.Orthogonal errors are minimized through maximum likelihood parameter estimation.
Because the relative magnitudes of orthogonal vector components are unit-dependent, we normalize the GPP data to scale from 0 to 1 locally to match EF.This treats GPP-EF error structure (including both observation and model structural errors) as proportional to observable range.While observation errors can be more directly estimated, we assume that model error structure is much larger than monthly-aggregated measurement error given its purely empirical nature.We then use "orthogonal least squares"-TLS with equal error scaling-to fit the curves, and then transform the parameters back to the standard GPP scale for inter-comparison and plotting.For each (EF, GPP) observation, the nearest point on the curve is required to be within the (rectangular) convex hull of that site's observations, both during parameter optimization and for the final fitted models.
To separate the roles of temporal GPP-EF covariance versus their spatial covariance, we fit two across-site relationships as follows.

General Temporal GPP-EF Relationship
We differentiate each site-specific curve in Figure 1c to determine dGPP/dEF as a function of local EF state.At each of fifteen evenly-spaced EF values, we select all sites spanning this EF value.We calculate dGPP/dEF for each site, and then determine the median of these slopes.These median slopes are shown in Figure 2a, either using only non-evergreen sites (66 sites) or using all sites.Again, recognizing the ambiguity of GPP-EF causality, we highlight that the median (unlike the mean) of a sample dGPP/dEF distribution is the same as the inverse of the median of a dEF/dGPP distribution.
To estimate uncertainty in this median curve, we resample the sites with replacement and calculate the median slopes (both for non-evergreen sites and for all sites).We repeat this 1,000 times to determine a 99% confidence interval for the median estimator at each of the 15 EF values (Figure 2a).
We then integrate the dGPP/dEF slopes with respect to EF for the median curve (no additive integration constants), as well as for the median curve for each bootstrap sample.This gives a median integrated curve and a distribution of re-integrated bootstrapped median curves (Figure 2b).

Spatial GPP-EF Relationship
To calculate long-term site averages, we again use only months with at least 25% of the local maximum GPP.The curves are again fit using total least squares, again scaling by the maximum observable GPP value.This spatial curve is shown in Figure 2c.The estimated confidence interval is again determined by bootstrapping.We randomly select sites with replacement, fit a spatial curve, and repeat 1,000 times.

Results
Figure 1a shows monthly GPP versus monthly EF at three example eddy covariance flux towers across the gradient of climatic water availability: Santa Rita Mesquite (US-SRM) in Arizona, USA (Barron-Gafford et al., 2013); Howard Springs (AU-How) in Northern Territory, Australia (Beringer et al., 2007); and Mead (US-Ne2) in Nebraska, USA (Suyker, 2016).For each site we fit a curve using an empirical power-law fit (Equation 4).Notably, each of the fit curves displays positive correlation and convex curvature.
The curves in Figure 1c show monthly water/carbon relationships site-by-site.At these time scales, correlations are positive (as opposed to the negative correlations that could come from water-use-efficiency vs. EF relationships).If the variation was caused exclusively by solar-driven seasonality, we could think of these curves as tracing out the annual dynamic equilibrium of the coupled water and carbon cycles locally.
We separate IGBP classes containing evergreen species (ENF, EBF, MF) from other classes to see the degree to which leaf phenology drives the relationship.Most curves in non-evergreen ecosystems follow the same general convex, positively-covarying relationship between carbon uptake and water availability seen for non-evergreen sites.Notably, evergreen sites with negative correlations (16 sites) are nearly all on Australia's south coast or in wet taiga regions, suggesting likely contamination of the evaporation signal from nearby surface water.Further tower footprint analyses may be able to clarify these relationships.Parameters a and b in Equation 4are not statistically different in the mean between evergreen and non-evergreen ecosystems (t-test, α = 0.05), but the exponential c-parameter is statistically distinct, with mean values of 4.7 for evergreen ecosystems and 5.9 across all sites.Across all sites, a regression of the estimated c-parameter versus site-mean EF is weakly positivelycorrelated (Pearson r = 0.35).See Figures S1 and S2 in Supporting Information S1 for the same analyses using un-normalized LE in place of EF.

Generalizing
The similarity in curve shapes across ecoclimate regions suggests that there may be utility in a simple, across-site observationally-driven and empirical model of the GPP-EF coupling.While it would be possible to average the curves in Figure 1 directly, this would convolve the temporal dynamics with the (successional-scale) mean-state information.To mitigate this effect, we work with the distribution of local, monthly dGPP/dEF slopes.
Figure 2a shows the median of the local derivatives, evaluated at even EF spacing (solid line for non-evergreens, dashed for all sites), along with a bootstrapped 99% confidence interval in the median.See Figure S3 in Supporting Information S1 for similar analysis for dGPP/dLE in place of dGPP/dEF.
Figure 2b shows the integrals of the median curves in Figure 2a.This provides an across-site relationship based only on monthly temporal dynamics rather than mean state information.As in Figure 2a, the reintegrated curves are quite similar when calculated from only non-evergreen sites and from all sites.Fitting power law curves to the median non-evergreen and all-site relations yields GPP = 7.65 × 10 5 + 29.4EF 5.77 , and GPP = 6.94 × 10 2 + 30.0EF 7.02 respectively, in units of g C /(m 2 • day).When used as models of GPP as a function of EF, they tend to underestimate, suggesting that some integration constant is likely appropriate for direct generalization.
The form of the GPP/EF relationship in Figure 2b is quite similar to those in Figures 1a-1c, presumably because (structural) biomass scaling is still implicitly at play in the differentiated relationships-even when the photosynthetic apparatus of chlorophyll is absent, the carbohydrate stores, root systems, and woody above-ground biomass of ecosystems are still present and can scale the relationship in many ecosystems.It is still noteworthy to see such across-biome similarity (e.g., grasslands, forests), and suggests that, while biome and plant function may covary with water availability on successional timescales, they are not a particularly strong determinant of landscape GPP/EF coupling or of landscape water-use efficiency (closely-related to dGPP/dLE).Biomass and potential vegetation are likely a larger factor in these relationships (of note, most of these sites are not recovering from major disturbance events).
Regardless of fit method (e.g., Equation 4 by total least-squares, Equation 4 by standard non-linear least-squares, linear least-squares of log-transformed data), a global power-law fit using all monthly data has correlations of ∼0.65 and ∼0.73 for GPP = f(EF) and EF = f(GPP), respectively for non-evergreen sites (0.53, and 0.57 when including evergreen sites).The difference in correlation reinforces our choice of total least-squares fitting, and implies that further investigation of appropriate empirical model and/or error structure is worthy of investigation.
Because of the strong relationship between monthly GPP and EF-even across sites with a simple general model, we further investigate the across-site long-term mean relationships.Since biomass is an integrated variable-on successional as well as sub-annual scales-this could be thought of as looking at consequences of decadal covariance using space-for-time substitutions, or as naively modeling landscape GPP-EF covariance as though it were driven by the same mechanisms at all tower sites (e.g., if all sites followed Figure 2b).
Figure 2c shows the long-term mean GPP and EF values for all sites.Power-law fits to the mean values show shallower curvature (c = 3.46 and c = 2.43 for non-evergreen and all-site curves, respectively).This is expected for longer time-aggregation.The average (x,y) of variables x and y with a convex relationship will always fall above the curve, and so averaging or aggregating variables with a nonlinear relationship will always push the relationship f(x,y) towards linearity.This is often stated as "Jensen's inequality" (Denny, 2017;Englund & Leonardsson, 2008) and is a consequence of the site-wise convex, local temporal dynamics.For this reason-even with similar monthly dynamics across sites-space-for-time substitutions will over-(under-) estimate GPP (EF) values.

Discussion
The correlations at a large number of sites show strong and similar covariance relationships between GPP and EF (e.g., Figure 1a).More broadly, a single parameterized relationship between the two variables explains more than half of the variance in monthly EF across climates and biomes.While the rough (positively-correlated) form of the relationship is not entirely surprising, the within and across-site coherence is unexpected since theoretical relationships between GPP and EF involve dozens of other complex variables and representation of plant trait heterogeneity.This emergent relationship suggests that some water/carbon coupling mechanisms may be dominant at sub-annual, landscape-aggregated scale.
EF as a water availability metric points to both water demand and supply, and cannot be used alone to fully disentangle the directionality of water/carbon coupling seen in this study.Rather than a shortcoming of EF as a water availability metric though, this is simply the consequence of the near-total entanglement of water supply and demand due to strong water-cycle coupling.EF has a two-way causal relationship with land surface temperature (Feldman et al., 2019(Feldman et al., , 2022)), surface-level vapor pressure deficits (Seneviratne et al., 2010;Stocker et al., 2018), the humidity profile of the atmospheric column (Rigden & Salvucci, 2015), transpiration (Rigden & Salvucci, 2017;Rigden et al., 2018), and soil moisture (Dong et al., 2022).This coupling is sufficiently strong that EF may be nearly informationally equivalent to relative humidity in theoretical terms (McColl & Tang, 2023).
The carbon cycle is entwined in this coupling as well, in regions with vegetative land cover and organic soils.
Because of the many coupled processes involved in water-carbon coupling, we view the similarity of the GPP-EF relationship across sites as an "emergent" property of these landscapes.Fully disentangling the primacy of each mechanism will require significant process representation and parameterization, likely in dynamic demographic models coupled to the soil and atmosphere.
It is clear that the data used in this analysis suffer from a lack of spatial representativeness, related to the siting of eddy covariance towers globally (see Figure 1b).We recommend extreme caution extrapolating any prognostic GPP-EF relation in Asia, Africa, or the global tropics (see Figure S6 in Supporting Information S1).In general, the relationships in this study are empirical and the product of many causal mechanisms, and are subject to the effects of spatial and temporal scaling.

Structural Role of Vegetation
It is intuitive that GPP l scales with the amount of green vegetation (e.g., LAI, f v )-greenness is the primary basis for remotely sensed GPP estimation.To the degree that there are functional relationships between GPP and greenness, it would not be clear if the GPP-EF relationships in Figures 1 and 2 are better viewed through GPP-EF coupling or f v -EF coupling.Figures S4 and S5 in Supporting Information S1 recreate the analyses of this study using remotely-sensed estimates of vegetation cover fraction and suggest that both temporal and spatial f v -EF relationships are strong and mutually reinforcing.The correspondence between EF and f v has considerably more scatter than between EF and GPP.Moreover, the curvatures (convexity) are not as clearly-linked between spatial and temporal relationships, suggesting more divergence in spatial and temporal mechanisms linking f v and EF relative to GPP.Evergreen vegetation, perhaps unsurprisingly, also diverges substantially in its temporal f v -EF relationships from non-evergreens.
Spatially, the convexity of the GPP-EF relationship can be viewed as the effect of positive covariance between terms on the left side of Equation 1.Since area-averaged transpiration and f v are both positively correlated with both water and energy availability, GPP l should have a steeper-than linear relationship with either variable alone.
A similar argument holds for θ and the number of summed terms in Equation 3.
Importantly, while these water, energy, and carbon variables are all strongly coupled, predictive power from alternative observational data (e.g., EF derived from weather station data) allows for independent GPP estimates, and can help us to more fully understand what remote sensing data sets are telling us about global photosynthetic rates.The similarity of relationships across bioclimate without any representation of either land-cover or plant functional traits suggests that vegetation-water interactions may be more similar in aggregate than between individuals.While there is a wealth of research linking water to vegetation on scales from leaf-to-landscape, more synthesis is needed in this area to fully separate time-scales of acclimation, phenology, plant hydraulics, and demography as they relate to the local water cycle.
As previously mentioned, the most direct mechanisms in the water-to-carbon direction are a positive correlation between EF and local water content (soil moisture and humidity); promoting higher stomatal-controlled carbon flux rates; reducing drought stress and moving towards light-limitation; upregulating photosynthesis; encouraging growth; and increasing leaf area, photochemical concentrations, and carbon conductance/uptake.This is one pathway for the canopy structure-reinforcement of the positive covariance.
In the carbon-to-water direction, the most direct mechanistic pathway for positive subseasonal water/carbon correlation is high carbon uptake leading to higher carbohydrate production and biomass growth (above and below-ground); increasing root and leaf densities; increasing water conductance from otherwise diffusion-limited root-zone soils; moistening the lower atmospheric boundary layer and cooling the land surface through transpirative latent heat flux; increasing evaporative fraction.In this direction, the EF = f(GPP) relationship is concave, suggesting diminishing effects of further carbon uptake on EF.
Both of these mechanistic pathways are certainly true.Both are also dependent on structural biomass-(or leafarea, root area, etc.) scaling and growth feedbacks (i.e., n chl in Equation 3) to create the nonlinear GPP-EF relationship.There are numerous additional feedbacks and mechanisms which contribute as well, but these are sufficient to hypothesize that growth feedbacks-even on subseasonal timescales-drive strong landscape water/ carbon coupling.This coupling on subseasonal scales reinforces similar processes on successionary scales, that is, water-limited biomass (A.Zhang et al., 2021), reproductive/recruitment success (Sapes & Sala, 2021), mortality (Anderegg et al., 2015), and landscape optimality theories (Franklin et al., 2020).The subseasonal coupling, may in fact, be mechanistically sufficient to determine these longer, successional processes, but this remains to be seen.
Additionally, in arid environments, growing seasons are adaptively aligned with wet seasons, which provides an exogenous driver of both EF and GPP.The similarity of responses across hydroclimate (Figures 1a-1c), however, suggests that this is not the sole or even dominant variable determining landscape water/carbon coupling.
Among the sites used in this analysis are 13 managed agricultural sites, many of which are also irrigated (e.g., Mead US-Ne2 in Figure 1a).These sites are not obviously different from typical patterns at other sites (apart from being in the wetter half of the data distribution), suggesting that the some of the across-site patterns we see likely arise due to processes other than succession.

Alternative Fitting Approaches
We fit the GPP-EF relationship using power-law curves for empirical simplicity rather than theory.This is a "process-blind" or data-driven approach so as to not impose our own assumptions onto the data.This is similarly the motivation for using total least-squares rather than simple nonlinear regression.It is worth noting, that residuals will be smaller for predicting one variable from another using normal nonlinear regression, and that power laws-while simple-are not theoretically tied to real processes.EF is nearly always bounded between 0 and 1, so a slightly more physically-constrained approach (with its own assumptions and caveats) might be logistic regression of EF on GPP.If we fit a function EF = β 0 + β 1 [1 + exp β 2 + β 3 GPP)] 1 directly using nonlinear regression to each site, we get slightly different, but qualitatively similar curves to those in Figure 1.A weighted average of these parameters across all sites (weighted by 1/σ 5→95 , where σ 5→95 is the span of the estimated parameter uncertainty from the 0.05 to 0.95 quantiles) gives parameters [β 1 ,β 2 ,β 3 ,β 4 ] = [ 0.032, 0.913, 1.725, 0.436].We present this only because it has the benefit of not exceeding EF = 1, but again, do not claim this to be tied directly to theory.

Conclusions
We present here simple empirical relationships between monthly evaporative fraction and gross primary productivity observations.These empirical relationships are good predictors of EF/GPP behavior at terrestrial eddy covariance towers, both locally and even globally.A single three-parameter model fit to monthly data from all non-evergreen sites has a correlation of 0.73, with no representation of land-cover, plant functional traits, land history, seasonality/phenology, canopy density/structure, photosynthetic pathway, nutrient availability, weather conditions, latitude, irrigation, or other management practices.
The local temporal relationships are convex and lead to an across-space convex relationship as well (Figure 2c, slightly linearized by averaging).Further investigations from scale transition theory (Bauke et al., 2022;Chesson, 2012;Yin & Porporato, 2023) could help to clarify the relationships between monthly and longer timescales.
The similarities in local dynamics and the existence of an across-space relationship-even in areas we might not think of as water-limited-suggests strong coupling of water and carbon across hydroclimates (c.f., Feldman et al., 2018Feldman et al., , 2021)), likely modulated by the spatial density of above-ground biomass.We hypothesize that proper representation of biomass/leaf-area dynamics with water-in particular on phenologic, water stress, and growth acclimation timescales-may be particularly critical in modeling large-scale water and carbon-cycle variables.
Finally, the response of ecosystems to disturbance may be particularly evident within this coupled framework, as water and carbon variables may respond to shocks on different timescales.This could help to pinpoint more theory-driven explanations of this emergent water/carbon coupling pattern.

Figure 1 .
Figure 1.Monthly coupling of the carbon and water cycles at example eddy covariance tower sites.(a) Monthly Gross Primary Productivity (GPP) versus monthly Evaporative Fraction (EF) at three sites.Carbon uptake and water availability covary strongly, with positive, convex local relationships across the range of seasonal and inter-annual variability.(b) Location of 183 tower sites used in this analysis.Blue circles denote non-evergreen ecosystems, green stars denote evergreen ecosystems.(C) Curve fits for all 183 towers as in (b).Most locations show similar positive-sloped convex curves, even in evergreen ecosystems.

Figure 2 .
Figure 2. Local temporal dynamics and site-means for all tower sites.(a) Slopes of GPP-EF relationship as a function of monthly EF.Solid line shows median curve of all slopes shown in Figure 1a (non-evergreen); dashed line shows median of all slopes for all IGBP classes.Gray shading shows bootstrapped 99% confidence interval of the non-evergreen median.(b) The re-integrated general GPP-EF curve, having removed all site-wise mean-state information.Solid line is the integrated slope curve from (a) for all non-evergreen sites.Dashed line shows the integrated slope for all sites.Shaded blue regions shows bootstrapped 99% confidence interval for the nonevergreen median curve.(c) The long-term site-means.Red circles show long-term mean of individual non-evergreen sites; green stars show evergreen sites.Solid line shows TLS fit to non-evergreen long-term means; dashed line shows fit for all sites.Red shading shows bootstrapped 99% confidence interval in the non-evergreen median curve.Note the continued positively-correlated convex relationship.The fit curve can be compared to those in Figures 1 and 2b if GPP-EF dynamics are assumed to be similar across landcover classes, biogeography, and climates, with important averaging/scaling caveats.