Groundwater effects on net primary productivity and soil organic carbon: a global analysis

Groundwater affects ecosystem services (ES) by altering critical zone ecohydrological and biogeochemical processes. Previous research has demonstrated significant and nonlinear impacts of shallow groundwater on ES regionally, but it remains unclear how groundwater affects ES at the global scale and how such effects respond to environmental factors. Here, we investigated global patterns of groundwater relationships with two ES indicators—net primary productivity (NPP) and soil organic carbon (SOC)—and analyzed underlying factors that mediated groundwater influences. We quantitatively compared multiple high-resolution (∼1 km) global datasets to characterize water table depth (WTD), NPP and SOC, and performed spatial simultaneous autoregressive modeling to test how selected predictors altered WTD-NPP and WTD-SOC relationships. Our results show widespread significant WTD-NPP correlations (61.5% of all basins globally) and WTD-SOC correlations (64.7% of basins globally). Negative WTD-NPP correlations, in which NPP decreased with rising groundwater, were more common than positive correlations (62.4% vs. 37.6%). However, positive WTD-SOC relationships, in which SOC increased with rising groundwater, were slightly more common (53.1%) than negative relationships (46.9%). Climate and land use (e.g., vegetation extent) were dominant factors mediating WTD-NPP and WTD-SOC relationships, whereas topography, soil type and irrigation were also significant factors yet with lesser effects. Climate also significantly constrained WTD-NPP and WTD-SOC relationships, suggesting stronger WTD-NPP and WTD-SOC relationships with increasing temperature. Our results highlight that the relationship of groundwater with ES such as NPP and SOC are spatially extensive at the global scale and are likely to be susceptible to ongoing and future climate and land-use changes.


Introduction
Groundwater is one of the most important freshwater resources, supporting agricultural, industrial and domestic water use (Siebert et al 2010, Gleeson et al 2020. However, accelerated withdrawal and depletion of groundwater, especially in arid and semiarid regions (Aeschbach-Hertig and Gleeson 2012), challenges the long-term sustainability of freshwater resources (Wada et al 2012, de Graaf et al 2019, Elshall et al 2020. By altering ecohydrological and biogeochemical processes across the critical zone (i.e., a permeable layer of the Earth that extends from the top of the canopy to the bottom of weathered bedrock), changes in the water table depth (WTD) can cascade to affect a range of groundwater-dependent or -mediated ecosystem services (ES). Past work has revealed that groundwater impacts ES including biomass production, nutrient cycling, streamflow regulation and greenhouse gas mitigation (Danielopol et al 2000, Qiu et al 2019, Zipper et al 2022a. However, the spatial variations in the magnitude and direction of interactions between groundwater and ES at global scales remain poorly understood. Primary productivity, as a fundamental indicator for many provisioning ES (e.g., crop and timber production), is sensitive to precipitation and limited by water availability (Churkina et al 1999, Wu et al 2011. In the sandy, humid forests of Wisconsin, USA, shallow groundwater (i.e., the near-surface portion of groundwater that can influence land surface processes) (Fan 2015, Hare et al 2021 increases tree growth by 63% compared to trees where groundwater was below the rooting depth (RD) (Ciruzzi and Loheide 2021). Nevertheless, groundwater effects on productivity are not always positive and vary in response to WTD and hydroclimatic conditions (Zipper et al 2015, Qiu et al 2019. In the cornfields of south-central Wisconsin, USA, both biophysical modeling and field experiments have revealed that shallow groundwater can subsidize crop water requirements and increase yield in dry years, while exerting a yield penalty in wet years by waterlogging and creating oxygen stress (Soylu et al 2014, Zipper et al 2015. A more recent global groundwatervegetation analysis depicted complex positive or negative groundwater relationships with gross primary productivity (GPP), which were spatially heterogeneous and seasonally dynamic, and varied with vegetation type and regional climate (Koirala et al 2017). However, global-scale patterns of groundwater effects on net primary productivity (NPP) remain unclear and may differ from GPP because (1) soil evaporation and plant transpiration are highly responsive to soil moisture that is susceptible to groundwater influences (Maxwell and Condon 2016) but not reflected in GPP; and (2) the ratio of NPP to GPP is not always constant spatially (Collalti and Prentice 2019) or temporally. Thus, we may expect that groundwater effects on NPP may show a different pattern from WTD-GPP covariation as previously revealed across the globe.
Groundwater also significantly affects regulating ES; for example, by influencing climate regulation through changes in soil carbon storage (Qiu et al 2019). Many regional studies (e.g., in coastal wetlands and riparian zones) have found that WTD was positively correlated with soil organic carbon (SOC) (Lyon et al 2011, Guan et al 2021, indicating greater SOC content with rising groundwater. In tropical peatlands, mainly located in Southeast Asia, soil CO 2 emissions can be amplified (i.e., leading to reduced SOC) with increases in WTD due to enhanced soil respiration and decomposition (Prananto et al 2020). NPP also provides a net carbon input into ecosystems (Chapin et al 2011), suggesting that positive groundwater effects on NPP could lead to more organic carbon input into soils. However, rising groundwater levels do not always improve SOC stocks and CO 2 emissions are not consistently dependent on WTD (Tiemeyer et al 2016). The net effects of groundwater on SOC, especially at continental and global scales, remain elusive and can be highly variable depending on complex interactions between WTD, WTDmediated soil abiotic and biotic conditions, C inputs and outputs, and climatic conditions (Meersmans et al 2011, Chen et al 2022. In tandem, groundwater effects on ES can be positive or negative, highly variable spatially, and likely differ by ES and mediated by environmental factors (Qiu et al 2019, Wen et al 2020. In a global climate change context, multiple environmental factors could interact to alter the dynamics of groundwater as well as its long-term impacts on ES, including climate (Cuthbert et al 2019), soil texture (Qiu et al 2019, Huang et al 2021 and topographic factors such as slope and elevation (Fan et al 2013, Maxwell andCondon 2016). Yet the global interplay between these potential drivers of WTD-ES relationships has not been thoroughly investigated. Hence, in this paper, we ask two main questions: (1) how does groundwater relate to NPP and SOC at the global scale? and (2) how do environmental factors (i.e., climate, topography, soil and land use/cover) mediate groundwater depth's relationships with NPP and SOC? To answer these questions, we first developed a transferrable analytical framework to quantify the relationship between ES and groundwater availability (indicated by WTD) (figures 1(a) and (b)), which can be further mediated by environmental factors at the global scale (figure 1(c)). We then empirically investigated these relationships by quantitatively comparing high-resolution global datasets of WTD and two ES indicators (NPP and SOC) and an array of selected environmental covariates.

Data sources and preprocessing
Our overall analytical procedure is shown in the flowchart (figure S1) and details of data sources are shown in table 1. The WTD dataset and the RD dataset were obtained from steady-state global model simulations (Fan et al 2013(Fan et al , 2017figures S2 and S3). The long-term NPP datasets were derived from the global annual products of Collection 5.5 MOD17A3, from the Numerical Terradynamic Simulation Group, University of Montana. To match the long-term average WTD dataset, we estimated long-term average global NPP as the mean from 16 years (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) of available NPP data for subsequent analyses ( figure  S4). SOC stock datasets were obtained from the SoilGrids250m (Hengl et al 2017). We calculated SOC for one shallow soil layer from 0 to 15 cm depth (i.e., SOC [0-15 cm] = SOC [0-5 cm] + SOC [5-15 cm] [100-200 cm] ), indicative of the total soil carbon pools ( figure S5). All these global datasets were from independent sources quantified with different approaches and/or estimated using different variables, thus avoiding any potential circularity in the analyses. While these datasets represented the bestavailable global products, there were some inherent challenges linking and analyzing them together. First, the underlying data used to generate these globally gridded products were spatially heterogeneous, with a bias towards more measurements and representation in temperate North America and Europe and fewer in Africa, Asia and high latitudes, which could cause inconsistent accuracy of WTD, NPP and SOC datasets. Second, these datasets were from multiple sources with different temporal coverages, making it impossible to match them temporally under a same time period. For example, the WTD and SOC datasets lumped together observations from many years to create a single static dataset, while mean annual NPPs were calculated by averaging all the annual estimates from 2000 to 2015. Given the spatial heterogeneity and temporal asynchrony of datasets, our study aims to depict the fundamental patterns of WTD-NPP and WTD-SOC relationships with typical data that reflected long-term natural conditions, rather than an analysis of inter-or intra-year dynamics.
All datasets were adjusted to a consistent spatial resolution of 30 arcsec (∼1 km). To filter out grid cells that were not expected to be influenced by groundwater, we calculated the gap distance between expected water table and RDs, since prior research has revealed that groundwater exerted the most profound impacts in shallow groundwater environments, but produced negligible influences in deep groundwater areas (Orellana et al 2012, Fan et al 2013, Qiu et al 2019. Here, we restrict our analysis to areas where the gap distance was ⩽5 m as a conservative threshold (more details in the supplementary information (SI)).

Spatial scale of analysis
The WTD-NPP and WTD-SOC relationships were estimated at the basin scale, as defined using level 07 of the HydroBASINS dataset (figure S6) (Linke et al 2019). We used HydroBASINS since it was derived from the same digital elevation model as the WTD dataset (Fan et al 2013). Basins with areas <100 km 2 or with <30 grid cells of WTD-NPP or WTD-SOC data were filtered out to avoid unrobust results (see the SI for details). After this screening, the final WTD-NPP dataset included 45 963 basins with an average basin size of 2581.3 km 2 and the final WTD-SOC dataset included 50 001 basins with an average basin size of 2585.2 km 2 .

Estimating WTD-NPP and WTD-SOC relationships
To address our first question, Spearman's correlation coefficient, ρ (x−y) , was calculated to quantify the basin-scale relationships between WTD (x) and NPP/SOC (y) across the globe. Spearman's correlation is robust to data distribution and has been previously used to estimate WTD-NPP and WTD-SOC relationships (Mitchell et al 2014, Koirala et al 2017. When calculating ρ ( x − y ) , the WTD value was transformed into a negative value (-WTD) so that interpretation of the ecohydrological implications of WTD-NPP and WTD-SOC correlations was more intuitive. Specifically, this meant that basins with ρ > 0 had greater NPP or SOC when the water table was closer to the land surface (i.e., positive NPP/SOC response to shallow groundwater, as depicted in figure 1(a)), and basins with ρ < 0 had greater NPP or SOC when the water table was deeper from the land surface (i.e., negative NPP/SOC response to shallow groundwater, as depicted figure 1(b)). We analyzed the distributions of Spearman's correlation, ρ (x−y) , within 20 • -latitude bins and across eight geographic regions as defined based on continents from HydroBASINS (excluding Greenland and Antarctica). As noted above, the confidence of estimated ρ (x−y) for each bin or geographic region was subject to confidence in the underlying data. Varying densities of WTD observations were used to generate the global WTD input dataset (see details in section S5 of the SI). Specifically, the confidence of simulated WTD data in North America, Australia and Europe could be higher than the confidence in Asia or Africa (table S1), indicating that the confidence of our derived WTD-NPP/SOC correlations may be similarly higher in North America, Australia and Europe than in other regions. In addition, the terrain and landscape factors affected the confidence of the simulated global WTD dataset (Fan et al 2013), and thus the confidence of our WTD-NPP/SOC correlations may be greater in areas with flatter topography or with more natural vegetation (see the SI for more details).

Analyzing effects of environmental factors on WTD-NPP and WTD-SOC relationships
In addressing our second question, to avoid mixing correlations with potentially different mechanistic reasons, we first separated all basins into positive and negative subsets based on the sign of their WTD-NPP and WTD-SOC correlations. We then extracted from HydroATLAS (Linke et al 2019) a total of 13 environmental factors that were hypothesized to mediate groundwater influences on NPP and SOC (table 2). All environmental factors were standardized before statistical modeling and their correlation matrix was constructed to identify and filter out highly collinear factors. We then built a fully linear model and used Moran's I to test for potential spatial autocorrelation. Since our analyses revealed significant spatial autocorrelations in residuals, we used the spatial simultaneous autoregressive (SAR) model for spatial regression, which was conducted with the spatial dependence (i.e., spdep) package in R (Bivand 2022). Depending on the volume of datasets and number of predictors, we either used the 'dredge' function in the Multi-Model Inference (i.e., MuMIn) package in R (Barton 2009) or a consistent model selection process (see the SI for details) to determine the optimal SAR model for each dataset that had the lowest Akaike's information criteria values (Burnham and Anderson 2004).
Based on the results of optimal SAR models, we then used scatterplots to visualize the responses of WTD-NPP and WTD-SOC correlations to the strongest mediating environmental factors. Since the responses could be highly scattered, we used the constraint line analysis to characterize the response curves (Cade and Noon 2003, Qiao et al 2019, Liu and Wu 2021. Specifically, adopted from Mills et al (2006), the segmented quantile regression method was used to generate response curves for the upper and lower boundaries of each scatter plot to characterize and visualize how environmental factors affected the strongest WTD-NPP or WTD-SOC relationships (see the SI for details). All the analyses were conducted in R 4.2.1 (R Core Team 2022).

Global distribution of WTD-NPP and WTD-SOC relationships
Globally, 61.5% of all basins showed significant correlations between simulated WTD and longterm mean NPP, with Spearman's correlation coefficients ranging from −0.94 to 0.83. There were more basins with negative WTD-NPP correlations than positive correlations (62.4% vs. 37.6%, respectively, out of basins with significant correlations; figure 2(a)). Negative WTD-NPP correlations (indicating decreased NPP at shallower WTD) were primarily located in high northern latitudes and aggregated clusters in southern Asia, mid-southern Africa and southern Australia. Positive WTD-NPP correlations (indicating increased NPP at shallower WTD) were mainly found in Central Europe, mid-western North America, western and central Australia, and Southern Africa. Median values and ranges of Spearman's correlation coefficients, ρ (WTD-NPP mean) , slightly decreased from the equator to the poles, with a similar unimodal distribution for each latitude group, but the distribution of ρ (WTD-NPP mean) in high latitudes showed overall smaller ranges ( figure 2(b)). The global pattern of WTD-NPP standard deviation correlations and distribution of ρ (WTD-NPP standard deviation) were broadly similar to those of WTD-NPP mean correlations, though certain hotpots of positive or negative correlations varied slightly ( figure S7).
For SOC, 64.7% of all basins showed significant correlations between WTD and SOC at the 0-15 cm depth, of which 53.1% were positive and 46.9% negative (figure 2(c)). The positive WTD-SOC correlations (indicating increased SOC with shallower WTD) were most common in arid or seasonally dry regions with diverse geography, such as mid-northern Europe, southern India and mid-western North America, and especially strong in the desert and steppe of Africa and Australia (figures 2(c) and 3(b)). Basins with negative WTD-SOC correlations (indicating decreased SOC with shallower WTD) were mainly located in mideastern and southern North America and southeastern Asia. The distributions of Spearman's correlation coefficients, ρ (WTD-SOC 015) were symmetric with long tails within each 20 • -latitude bin, except for 70 • -90 • North and 50 • -70 • South (figure 2(d)). Median values of the correlation coefficients were generally around zero but decreased at lower latitudes in the Southern Hemisphere. Hotspots of WTD-SOC correlations in the deep soil layer (i.e., 0-200 cm) and distribution patterns of ρ (WTD-SOC 0200) were similar to the results of WTD-SOC relationships at the 0-15 cm depth (SI section S9). However, there were more positive WTD-SOC relationships at the 0-200 cm soil depth, where 70.6% of all basins with significant WTD-SOC correlations showed positive relationships at the global scale (figure S7).
At regional to continental scales, most correlation coefficient distributions were symmetric with long tails spanning nearly the full −1.0-+1.0 range, though distributions for the Arctic and Siberia have smaller ranges than in other regions (figures 3 and S8). However, as noted previously, empirical training and validation data for WTD, SOC and NPP datasets were more limited in these regions, which may affect accuracy and reduce confidence in our findings for these geographic settings. For WTD-NPP correlations, coefficient medians of different regions were consistently negative, aligning with the result that negatively correlated basins were more common across the globe ( figure 3(a)). However, medians of WTD-SOC correlation coefficients were essentially zero except for Africa, which was slightly positive ( figure 3(b)).

Environmental factors mediating WTD-NPP and WTD-SOC relationships
Across the subset of basins with positive or negative WTD-NPP and WTD-SOC relationships, the optimal SAR models showed pseudo R 2 values ranging from 0.19 to 0.29 (table S3). The pseudo R 2 values ranged from 0.21 to 0.29 for WTD-NPP mean and WTD-SOC (at 0-15 cm depth) relationships (figure 4), while pseudo R 2 values ranged from 0.19 to 0.29 for WTD-NPP standard deviation and WTD-SOC (at 0-200 cm depth) relationships (figure S10). Overall, climate factors were the dominant predictors for WTD-NPP and WTD-SOC relationships at the basin scale across the globe. By comparing all predictors in each optimal SAR model, the annual average air temperature showed the strongest influence on mediating WTD-NPP/SOC relationships (figures 4 and S10). Specifically, the strength of positive correlations between WTD vs. NPP and SOC responded positively to temperature, indicating that at higher temperatures, greater NPP and SOC were associated with shallower WTD. On the other hand, for basins with negative WTD-NPP and WTD-SOC relationships, the strength of negative correlations responded negatively to temperature, meaning that at higher temperatures, NPP and SOC were more strongly negatively related with WTD.
Segmented quantile regressions showed that optimal constraint lines (i.e., those with the highest R 2 ) were quadratic for both NPP and SOC. Constraint effects of the annual average air temperature on the potential maximum of positive WTD-NPP/SOC correlations (i.e., points on the upper constraint lines) first increased and then flattened (figures 5(a), (d) and S11(a), (d)). Temperature constraint effects were symmetric for positive and negative WTD-NPP/SOC correlations, and the potential minimum negative WTD-NPP and WTD-SOC coefficients (i.e., points on the lower constraint lines) first decrease and then flattened with rising temperature. Collectively, these results showed that higher temperatures tended to either boost positive effects of groundwater or amplify negative groundwater influences.
Most optimal SAR models for NPP and SOC also contained annual average precipitation as a significant predictor. Precipitation was generally a weaker predictor and showed opposite effects from temperature; for instance, positive WTD-NPP/SOC relationships decreased in magnitude as precipitation increases (figures 4 and S10). Such relationships were consistent with the effects of the global aridity index on WTD-NPP and WTD-SOC correlations (figures 4 and S10), with data points more clustered in the lower range of the global aridity index (figures 5(c), (f) and S11(c), (f)). Segmented quantile regressions showed that the potential maximum of positive WTD-NPP/SOC correlations decreased, and potential minimum of negative WTD-NPP/SOC correlations increased with rising precipitation (figures 5(b), (e) and S11(b), (e)). Collectively, these results showed that stronger precipitation or greater humidity both weakened positive and negative groundwater effects on NPP and SOC.
The correlations with land use differed between WTD-NPP and WTD-SOC, and mainly reflected vegetation extent including natural, semi-natural vegetation and cultivated area (figures 4 and S10). Specifically, the magnitude of positive WTD-NPP correlations was negatively correlated with natural vegetation extent, indicating that potential increases in NPP associated with shallower WTD were lower in settings with more natural vegetation cover. Meanwhile, the magnitude of negative WTD-NPP correlations were positively correlated with natural, semi-natural and cultivated area extent, meaning that negative WTD-NPP relationships were weaker in settings with increasing extent of these cover types. Consistently, both results suggested that groundwater relationships with NPP were weaker in basins with higher vegetation extent. In contrast, models of WTD-SOC correlations contained fewer land use factors with smaller effects. Overall, our results showed that increasing vegetation extent had the potential to enhance positive groundwater effects, while weakening negative groundwater effects on SOC. Topographic factors, irrigated area extent, soil erosion, and soil texture were not retained in the final optimal SAR models, with weak to negligible relationships across all candidate models.

Global pattern of WTD-NPP and WTD-SOC relationships
Groundwater depth variation was significantly correlated with NPP and SOC at the basin scale across the globe, indicating a strong groundwater influence on these ES indicators. However, WTD-NPP and WTD-SOC correlations were spatially heterogeneous and differed between NPP and SOC, indicating that groundwater effects on NPP and SOC are highly context-specific and mediated by environmental factors such as climate or ecohydrological conditions. For instance, WTD can affect plant transpiration and soil evaporation (Maxwell and Condon 2016), which could lead to effects on plant growth by modifying root water availability (Soylu et al 2014). Similar to global WTD-GPP correlations revealed in Koirala et al (2017), we found that positive WTD-NPP correlations were dominant in arid and semiarid regions (e.g., the Great Plains of North America, the Eurasian Pontic-Caspian steppe; figure 2(a)), where shallow groundwater can provide a water subsidy for enhancing productivity during dry periods (Zipper et al 2015, Rizzo et al 2018. In contrast, more negative WTD-NPP correlations were distributed in humid and water abundant regions (e.g., South Asia, northern Eurasia, figure 2(a)), where increasing WTD could lead to oxygen stress from waterlogging that, in turn, negatively affects NPP (Wen et al 2020).
While we generally observed similar WTD-NPP patterns to the WTD-GPP relationships found in Koirala et al (2017), there are key distinctions between WTD, NPP and GPP datasets used in our work as well as our processing workflows for analyzing WTD-ES relationships. To complement our WTD-NPP analysis and better compare findings from Koirala et al (2017), we conducted an additional global analysis on WTD-GPP relationships using our datasets and workflow (details and full results in section S10 in SI). Essentially, when comparing WTD-NPP (figure 2) and WTD-GPP (figure S9) patterns in our study, we found primarily similar patterns globally but some minor regional differences. For example, in eastern North America, WTD-NPP relationships were mostly negative, while WTD-GPP relations were mostly positive. Differences between WTD-NPP and WTD-GPP relationships may be due to large-scale spatial variability of NPP:GPP ratios (Zhang et al 2009). This work found that NPP:GPP ratios were low (from 0.45 to 0.6) in eastern Northern America, which could cause differences in the WTD-NPP and WTD-GPP correlations in this region. The strong WTD-NPP and WTD-GPP coupling as we observed at the global scale suggests that shallow groundwater may be a critical explanatory factor for spatial heterogeneity of NPP: GPP ratios by changing shallow soil moisture that influences autotrophic respiration and CO 2 fluxes (Green et al 2019). However, more work would be needed to investigate the strength of such interactions at the global scale.
Comparing our WTD-GPP patterns with those in Koirala et al (2017), we also found spatially similar patterns with some notable regional differences. For example, we found more basins with significant WTD-GPP correlations (either positive or negative) in our analysis, especially in the middle of North America and some parts of southern Asia (figure S9). These results suggest that the higher-resolution gridded GPP datasets we used (∼1 km resolution, compared to ∼10 km from Koirala et al 2017) can capture fine-scale ecosystem responses to WTD that can be masked in coarser resolution analyses. Further, Koirala et al (2017) did not investigate groundwater effects on GPP above 60 • N, where we observed widespread negative WTD-NPP and WTD-GPP relationships in perennial cold and humid tundra regions at high Northern Hemisphere latitudes (figures 2(a) and (b)), though this finding needs to be cautiously interpreted due to the underrepresentation of empirical observations for developing the WTD and NPP datasets for these high latitude regions. These results suggest potential negative relationships between excessive water availability caused by shallow groundwater and plant productivity extend outside of temperate environments where most research on the topic has occurred (Liu et al 2020). Overall, our work complements and expands on the important past findings of Koirala et al (2017).
As an important predictor of SOC, groundwater can significantly affect greenhouse gas fluxes (Turetsky et al 2014) and influence the dynamics and spatial variations in SOC. Our results revealed more basins with positive WTD-SOC relationships at the global scale, both for the most active soil layer (0-15 cm, figure 2(c)) and deeper SOC pools (0-200 cm, figure S7(c)). Such positive WTD-SOC correlations were most common in arid or seasonally dry regions with diverse geography, especially in the desert and steppe areas in Africa (figures 2(c) and 3(b)). Based on the dynamics observed in regional studies (Lyon et al 2011, Guan et al 2021, the underlying mechanism for positive WTD-SOC may be that when WTD is shallower, pore space oxygen deficiency due to increased soil moisture inhibits decomposition of organic matter, which causes SOC accumulation (Callesen et al 2003). Moreover, groundwater could regulate carbon exchange indirectly by regulating nitrogen availability in the plant-soil system and a shallower WTD may enhance CO 2 inputs from GPP over CO 2 emissions by ecosystem respiration, ultimately affecting SOC (Pohl et al 2015). Variations in deep groundwater can still impact SOC in the top soil layer (Poeplau et al 2020), through regulating heterotrophic respiration (i.e., non-plant root respiration, Prananto et al 2020). However, WTD effects on different soil gas fluxes in previous studies have been inconsistent and nonlinear (Pohl et al 2015), or not even correlated (Tiemeyer et al 2016). For example, in seasonally well-drained wetlands without perennial frozen soil at high latitudes, an elevated WTD potentially promoted soil CH 4 emissions and thus reduced SOC (Turetsky et al 2014), which may explain negative WTD-SOC correlations in some high-latitude Northern Hemisphere basins (figure 2(c)).

Environmental factors mediating WTD-NPP and WTD-SOC relationships
Optimal SAR regression models demonstrated that multiple environmental factors could mediate WTD-NPP and WTD-SOC correlations, of which climate and land use were the most notable factors (figures 4, S10 and table S3).
Specifically, climate effects had a similar mediating role for both WTD-NPP and WTD-SOC relationships (figures 4 and S10). For basins with an annual average air temperature >16 • C (thresholds of constraint effects were shown in table S4), a shallower WTD either promoted positive groundwater effects on ES or amplified negative effects. These results suggested that rising water tables subsidize crop production by enhancing root water uptake to promote NPP or SOC in warmer regions (Lowry andLoheide 2010, Jobbágy et al 2011), while in cooler settings shallow groundwater introduced anaerobic effects to reduce evapotranspiration, NPP and reduced SOC accumulation (Soylu et al 2014). In contrast, for basins with higher annual average precipitation (or higher humidity), both positive and negative WTD-NPP/SOC correlations were weakened. Since vegetation prefers to obtain water from the soil water caused by recent precipitation infiltration or past infiltration stored in deep soil (Miguez-Macho and Fan 2021), increasing precipitation and air moisture may diminish the relative importance of groundwater, but the extent of diminishment of groundwater may differ among climate zones. It is noteworthy that under high-temperature conditions, quadratic constraint regression lines turned from upward to flat (figures 5(a), (d) and S11(a), (d)), indicating that potential maximum groundwater effects on NPP and SOC were limited and groundwater's influences on NPP and SOC may be lessened in future warmer climate conditions.
In contrast to climate, the effects of land use on WTD-NPP and WTD-SOC relationships were different between NPP and SOC, with mediating effects from four types of land use extent (figure 4). Weaker positive and negative WTD-NPP correlations were associated with higher vegetation extent in basins. Hence, vegetation extent showed similar mediating effects on WTD-NPP relationships as precipitation, potentially indicating that the extent of vegetation, especially natural vegetation, may weaken groundwater effects on NPP through increasing forest rainfall (Staal et al 2018) or enhancing rainfall interception (Porada et al 2018). For WTD-SOC correlations, greater vegetation extent appeared to promote positive groundwater effects in basins with a positive WTD-SOC correlation, while weakening negative groundwater effects in basins with a negative WTD-SOC correlation. Similarly, the normalized vegetation index has been proven to have a strong positive correlation with SOC (Bangroo et al 2020), indicating that vegetation extent should be highly positively correlated with SOC.

Implications for groundwater management and research limitations
Our study revealed widespread basin-scale relationships between groundwater and NPP/SOC across multiple ecosystems globally, though the conclusions drawn are necessarily limited by the quality of the input datasets. Spatial patterns of WTD-NPP and WTD-SOC correlations we documented could help guide groundwater management and policy development in different regions (Cord et al 2017).
For example, in regions where positive WTD-NPP and WTD-SOC correlations were concentrated (e.g., central North America), with high temperature and dry conditions, long-term attention should be paid to the sustainable use and extraction of groundwater (Gleeson et al 2010). Drawdown of the water table in these regions can lead to unexpected cascading effects on multiple ES beyond freshwater supply, for example, NPP and SOC as revealed in this study, along with streamflow depletion (Zipper et al 2021(Zipper et al , 2022b) and a loss of stream species (Stubbington et al 2020). In contrast, for basins in South Asia, where negative WTD-NPP and WTD-SOC relationships were predominant, measures, such as drainage or alternate irrigation during land management to lower WTD may have benefits (Sarker et al 2020), though drainage could also be associated with ecosystem disservices such as impaired surface water quality (Gramlich et al 2018, Nazari et al 2021. Due to the bias of WTD, NPP and SOC observational data towards temperate latitudes in North America and Europe, the implications and mechanisms of these WTD-NPP/SOC correlations may be more robust for these geographic regions. Our finding that WTD-NPP and WTD-SOC correlations varied in response to climate factors will provide insights for understanding potential changes of groundwater effects under future climate change (e.g., warming and more frequent extreme events, such as drought). Among climatic factors, temperature emerged as the most important driver mediating groundwater correlations with both NPP and SOC, indicating that regions more susceptible to future global warming would experience the most drastic changes in the influences of groundwater. For example, in regions with projected temperature increases, the importance of groundwater on sustaining NPP or SOC will also be likely to increase, but, ironically, these regions may at the same time experience lower groundwater levels from reduced groundwater recharge and increased irrigation requirements (Green et al 2011, Crosbie et al 2013, Rodell et al 2018. Previous work has found that severe droughts caused multi-year declines in terrestrial NPP and could further degrade terrestrial carbon pools (Zhao and Running 2010), but subsidies from sustained groundwater levels may effectively moderate this scenario (Vicente-Serrano et al 2020). In addition, substantial spatial variations in groundwater effects and their response to climate factors highlight that groundwater effects tend not to be static and thus a more dynamic approach to understanding and managing groundwater to enhance the sustainable provision of ES may be required (Elshall et al 2020).
Our research has some limitations. Due to data availability, only groundwater effects on NPP and SOC were considered in our analysis; thus, groundwater effects on a broad array of other ES still require future investigations. Past local-to regional-scale work has demonstrated that changes in WTD can interact with SOC to further influence crop yield and biogeochemical cycling (Pohl et al 2015, Huang et al 2021. Hence, groundwater influences on synergies or tradeoffs among multiple ES needs to be further explored. Additionally, we used static, global WTD and SOC products derived from observational data that are not uniformly distributed (Zhao et al 2005, Fan et al 2013, Hengl et al 2017 and were therefore unable to identify how NPP and SOC respond to interannual or seasonal groundwater level changes that can be caused by natural fluctuations or artificial disturbances. These interannual WTD changes and effects could be a major concern for water use policy or groundwater management (Tiemeyer et al 2016), and merit future investigation (e.g., using coupled dynamic large-scale groundwater-ecosystem models or future datasets with improved spatial and temporal resolutions). The development of such interdisciplinary research on the dynamics of groundwater and its social-ecological effects is necessary in the context of environmental change and socioeconomic development.

Conclusion
Our work revealed that groundwater effects on NPP and SOC were spatially heterogeneous and differed depending on the type of ES across the globe. We identified multiple environmental factors that mediated groundwater relationships with two ES indicators, particularly for climate (temperature, precipitation and aridity) and land use (vegetation extent) factors. Constraint effects on WTD-NPP/SOC correlations imply that climate changes may saturate potential groundwater relationships with NPP and SOC. In tandem, our work suggests that future ES research should address potential conservation and holistic management implications of changes to groundwater flows and levels. The widespread WTD-NPP and WTD-SOC relationships we detected indicate that climate and land-use changes that alter groundwater dynamics may have cascading impacts on the sustainable provision of ES through different biogeochemical and hydrological pathways. Our study highlights the importance of integrating multiple sources and disparate datasets to understand large-scale ES patterns and dynamics, their environmental controls and responses to current and future environmental changes.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files). The R code and data for statistical analysis can be downloaded in: https://doi.org/10.6084/m9. figshare.23694354.