Overlooked risks and opportunities in groundwatersheds of the world’s protected areas

Protected areas are a key tool for conserving biodiversity, sustaining ecosystem services and improving human well-being. Global initiatives that aim to expand and connect protected areas generally focus on controlling ‘above ground’ impacts such as land use, overlooking the potential for human actions in adjacent areas to affect protected areas through groundwater flow. Here we assess the potential extent of these impacts by mapping the groundwatersheds of the world’s protected areas. We find that 85% of protected areas with groundwater-dependent ecosystems have groundwatersheds that are underprotected, meaning that some portion of the groundwatershed lies outside of the protected area. Half of all protected areas have a groundwatershed with a spatial extent that lies mostly (at least 50%) outside of the protected area’s boundary. These findings highlight a widespread potential risk to protected areas from activities affecting groundwater outside protected areas, underscoring the need for groundwatershed-based conservation and management measures. Delineating groundwatersheds can catalyse needed discussions about protected area connectivity and robustness, and groundwatershed conservation and management can help protect groundwater-dependent ecosystems from external threats. Global initiatives to expand protected areas focus on controlling ‘above ground’ impacts such as land use, overlooking the potential human impacts on protected areas through groundwater flow. This study analyses the potential extent of these impacts by mapping groundwatersheds.

Protected areas are a key tool for conserving biodiversity, sustaining ecosystem services and improving human well-being. Global initiatives that aim to expand and connect protected areas generally focus on controlling 'above ground' impacts such as land use, overlooking the potential for human actions in adjacent areas to affect protected areas through groundwater flow. Here we assess the potential extent of these impacts by mapping the groundwatersheds of the world's protected areas. We find that 85% of protected areas with groundwater-dependent ecosystems have groundwatersheds that are underprotected, meaning that some portion of the groundwatershed lies outside of the protected area. Half of all protected areas have a groundwatershed with a spatial extent that lies mostly (at least 50%) outside of the protected area's boundary. These findings highlight a widespread potential risk to protected areas from activities affecting groundwater outside protected areas, underscoring the need for groundwatershed-based conservation and management measures. Delineating groundwatersheds can catalyse needed discussions about protected area connectivity and robustness, and groundwatershed conservation and management can help protect groundwater-dependent ecosystems from external threats.
Protected and conserved areas are fundamental tools for safeguarding biodiversity and play an important role in improving human well-being and sustaining ecosystem services [1][2][3][4][5] . With the newly agreed Convention on Biological Diversity's Kunming-Montreal Global Biodiversity Framework, Parties to the Convention have committed to ensuring and enabling that by 2030 at least 30 percent of terrestrial, inland water, and of coastal and marine areas are effectively conserved through protected areas and other effective area-based conservation measures (OECMs) 6 . Current and past approaches to inland protection, which have focused principally on land-based objectives and measures, have had clear limitations in conserving freshwater ecosystems and species, which have shown staggering declines 7,8 . One frequently discussed reason for Analysis https://doi.org/10.1038/s41893-023-01086-9

Groundwatersheds are often larger than their protected areas
The combined size of the groundwatersheds of the world's protected areas is 75% larger (22.0 million km 2 , n = 32,490 groundwatersheds) than the combined size of the protected areas that we analysed (12.6 million km 2 ). Most groundwatersheds (85%; 27,651 of 32,490) extend beyond their protected area boundary. Larger protected areas generally have larger groundwatersheds ( Supplementary Fig. 2). Thus, patterns in total groundwatershed size largely reflect patterns in protected area size. Groundwatersheds also span international borders and raise transboundary management concerns. We find that 1,379 groundwatersheds cross international borders, 454 of which do so despite their associated protected area existing entirely inside a single country.
The median RGS is 1.46 (Fig. 3a), with an interquartile range of 1.17-1.94. Overall, RGS tends to be larger in arid regions (Fig. 3b), which means that the size of the area contributing groundwater flow is greater in proportion to the size of the GDEs they are connected to in arid regions compared with more humid regions. These larger RGS values in arid regions are consistent with previous modelling of the impact of aridity on regional groundwater flow 23,24 . Larger RGS values in arid regions (for example, Fig. 3e) suggest that nested and regional flow paths, which are not represented in our water table-based approach, are particularly important in these settings (Supplementary Section 1.1). Lower RGS values, as found in the boreal forest of central North America (for example, Fig. 3c), correspond to groundwatersheds where vegetation is highly connected to shallower water tables. In these humid regions, convergence patterns in the water table are more localized and lead to smaller, shallow flow systems.
Trends in RGS do not differentiate a vulnerability gradient in protected areas but rather provide insights regarding the hydrogeological systems the protected areas depend on, and also provide context to inform protection strategies. That groundwatersheds exist outside of protected areas may appear as an intuitive finding, given that protected areas are rarely established on the basis of hydrological system boundaries or processes. Yet the global prevalence of this misalignment between protected areas and their groundwatersheds necessitates that these ecologically important areas of contributing groundwater flow be considered in conservation priority-setting, and RGS is a metric to help inform these efforts.
For instance, larger RGS generally implies that there is a larger area of contributing groundwater flow to manage. Further, larger the poor performance of protection initiatives for these ecosystems is the lack of consideration given to external (upstream) hydrologically connected freshwater systems 9,10 . The need to manage human activities in connected lands and waters outside protected areas has been absent from effectiveness discussions and indicators 11 , and the designation of new protected areas is rarely based solely on hydrologic boundaries. Connecting protected area initiatives to surface water processes is an important step; and the potential is apparent across the World Database on Protected Areas, which intersect with hundreds of watersheds.
However, consideration must also be given to groundwater systems. Protected areas face impacts from activities occurring outside of the protected area, such as agricultural drainage, mining, contamination and groundwater pumping, which are transmitted to protected areas through connected groundwater systems. Doñana National Park (Spain) 12,13 and Grand Canyon National Park (USA) 14 are two iconic examples where groundwater-mediated impacts have been documented.
The consideration and management of groundwater becomes increasingly important as land and water use intensify around protected areas 2,15 . Yet, no study has systematically investigated the potential for human activities outside protected areas to impact terrestrial and aquatic groundwater-dependent ecosystems (GDEs) in protected areas through groundwater flow (Fig. 1). Lateral groundwater flow supplies a substantial proportion of water used by vegetation 16 , and changes in land use or land cover can impact downgradient terrestrial ecosystems by changing the quantity and distribution of groundwater 17,18 . Groundwater pumping can reduce streamflow and change streamflow from perennial to intermittent or ephemeral 19,20 . Since groundwater has distinct chemical and temperature characteristics and can transmit contaminants such as nutrients 21 , changes in groundwater levels and flow can introduce pollutants or otherwise alter water quality in protected areas 22 .
In this study, we estimate the area from which human impacts outside protected areas may propagate to protected areas through groundwater flow systems. We employ a generic, reproducible workflow to map groundwatersheds for GDEs in protected areas (Box 1 and Fig. 2) at the spatial resolution of 30 arcseconds or ~1 km at the equator (see Methods). We conclude by identifying risks posed to existing protected areas based on levels of human activity and land use modification in the underprotected portions of groundwatersheds and discuss opportunities for improved conservation outcomes.   34,[36][37][38] . b, Area distribution of GDE types. c, Spatial distribution of protected areas 40 . In this study, we map the groundwatersheds of the 'high protection' group of protected areas and not for the 'low protection' group. d, Area distribution of all GDEs across protected area groups. e, Area distribution of protection groups across each GDE type. Note that the vertical axis is the same across b and e (for example, ~6% of all GDEs are lentic with high levels of protection).

Analysis
https://doi.org/10.1038/s41893-023-01086-9 groundwater flow systems generally have longer system response and residence times 25 , meaning that human impacts in larger groundwatersheds may potentially have longer legacy impacts on GDEs than in smaller groundwatersheds. Conversely, protected areas with smaller RGS are typically in regions with a greater density of GDEs. Generally smaller groundwatershed sizes in these humid regions mean that human impacts in these groundwatersheds may more rapidly affect GDEs.
The median UGR is 0.52 (Fig. 3f), with an interquartile range of 0.17-0.80. This means that the median protected area's groundwatershed extent exists 52% outside of the protected area boundary. There are no regional trends visible in UGR; however, we do find that larger protected areas generally have lower UGRs in comparison with smaller protected areas ( Supplementary Fig. 2), meaning that larger protected areas tend to have more-protected groundwatersheds.
UGR does not vary substantially with RGS ( Fig. 3g), nor with national levels of land protection ( Supplementary Fig. 3). Even in countries where international area-based conservation targets have been met, groundwatersheds of protected areas remain underprotected by conventional protected area initiatives. Combined, these findings reveal a global misalignment between protected areas and their connected groundwater flow systems and underscore the challenge of conserving protected area ecosystems above and below ground without consideration of groundwatersheds.

Human activity in groundwatersheds may undermine protection
GDEs inside protected areas could be affected directly by external groundwater pumping and contamination, and indirectly via climate change or land use change through their impact on groundwater recharge [26][27][28] . Activities such as mining, agriculture and urban expansion are key determinants of the potential risk to the quality and quantity of groundwater flow to protected areas. Thus, human activity and land modification in underprotected portions of groundwatersheds represent a potential vulnerability for GDEs in protected areas (Fig. 4a). The timing and severity of these impacts are a function of the type, location and magnitude of the activities in conjunction with the local hydrogeological setting. Assessing this timing and severity of impacts Box 1 What are groundwatersheds?
Note: Key terms are bolded and summarized at the bottom of the box.
A groundwatershed is the contributing area from which a groundwater system flows to a feature or set of features of interest (Fig. 2). In this respect, groundwatersheds are the groundwater analogue of surface watersheds. The groundwatershed concept was first introduced by Haitjema 57 to evaluate groundwater residence times, and similar concepts have been called 'groundwater catchments' 58 , 'groundwater basins' 59 and mapping of 'groundwater divides' 60 . However, the concept has seen limited uptake in water science and management, possibly owing to groundwater being an often-overlooked resource and also possibly due to some characteristic differences between groundwatersheds and surface watersheds, such as what features are used as outlet points, that make their delineation and use more challenging.
In arid environments, flat topographies and regions with complex geologies, groundwatershed divides can be spatially unaligned with surface watershed divides 61 . Groundwatersheds also differ from surface watersheds in their ability to fluctuate with time. Whereas surface watersheds are defined by static topography, groundwatersheds are dynamic and their size and shape can change due to pumping, climate change, land use change or seasonality. Therefore, groundwatersheds can be affected by a multitude of natural and human factors. However, in this analysis, we expect the majority of each groundwatershed's spatial extent to be consistent through time as fluctuations in the water table will only correspond to changes in the groundwatershed extent if the locations of water table divides are altered.
Here we derive groundwatersheds using the water table instead of the land surface in a standard watershed delineation algorithm (see Methods). Using the water table to derive groundwatersheds enables a computationally simple approach to delineate groundwatershed extents. This approach generates groundwatersheds that reflect shallow, local groundwater flow systems (Supplementary Fig. 1) but it does not represent nested regional groundwater flow systems that require particle-tracking simulations that are currently infeasible at the global scale.
The groundwatersheds we derive in this study are for the world's protected areas. Unlike surface watersheds, we do not use river outlets as outlet locations in our groundwatershed delineation approach. Instead, we use the locations of GDEs ( Fig. 1a and described below) that lie inside protected areas (Fig. 1c) as outlet features. Thus, we do not derive contributing areas of groundwater flow for one location per protected area, but for all GDEs inside each protected area.
We identify groundwater-dependent terrestrial and aquatic ecosystems using datasets of groundwater-dependent wetlands, root zone intersections with the water table, and surface water features such as perennial rivers and streams, since there is no existing global dataset of GDEs 62 . The GDEs we identify in this study do not represent a comprehensive and refined global database but rather indicate locations where these ecosystems can potentially occur.
To evaluate the potential importance of groundwatersheds and analyse their relationship with protected areas globally, we defined two metrics (Fig. 2c): relative groundwatershed size (RGS) and the underprotected groundwatershed ratio (UGR). RGS is an ecohydrological index representing size of the groundwatershed relative to the size of GDEs inside the groundwatershed. UGR is primarily a socio-hydrological conservation index that represents the underprotected proportion of each groundwatershed.
Groundwatersheds: Contributing areas of shallow local groundwater flow to a feature.
Protected areas: IUCN protected area categories Ia, Ib, II and III (mapped as the 'high protection' group of protected areas in Fig. 1c).
Underprotected areas: All areas outside of protected areas, as defined above. We use the term underprotected rather than unprotected as other forms of protection and management can exist in these areas, such as the European Union Water Framework Directive or California's Sustainable Groundwater Management Act.
GDEs: Terrestrial or aquatic ecosystems that contain species or habitats that rely on groundwater 62 . Lotic GDEs have running water (for example, rivers and streams), whereas lentic GDEs refer to those with standing waters (for example, lakes and wetlands).
Analysis https://doi.org/10.1038/s41893-023-01086-9 is beyond the scope of this study but could enable improved management as we describe below.
Comparing patterns in human land modification 29 with UGR provides insight on the vulnerability of protected areas (Fig. 4b,c). To assess regional patterns, we summarized the level of land modification and UGR per terrestrial ecoregion. Ecoregions where we find high human modification levels and underprotected groundwatersheds are probably areas where GDEs inside protected areas are most vulnerable to potential impacts through groundwater flow. These ecoregions are scattered across the world and include: the USA's Midwest to east coast, Central America,    coastal Brazil, the majority of Europe, northern Africa, the Sahel and Sudanian savanna in sub-Saharan Africa, Iran, Pakistan and Northern India.
We have contained our analysis to higher levels of protection (International Union for Conservation of Nature (IUCN) management categories I-III). However, we found most groundwatersheds to remain underprotected even when expanding our analysis to include all IUCN management categories (categories I-VI). The median national percentage of underprotected groundwatershed surface area that is in fact protected by lower levels of protection (categories IV-VI) is only 5%. However, Germany, Guatemala, Morocco, Myanmar, South Korea and Uruguay are among the few nations whose lower levels of protected areas cover substantial portions (all over 30%) of their underprotected groundwatershed area identified in this analysis ( Supplementary Fig. 4).

Discussion
Here we have used groundwatersheds to reveal the global potential for distant and long-term subsurface impacts on GDEs inside protected areas. Yet, our mapping of groundwatersheds for protected areas is only one of many possible applications and this work can serve as a proof-of-concept for wider use. Groundwatersheds, such as surface watersheds, can be identified for any feature, such as a protected area, groundwater well, wetland or community. Groundwatersheds have a strong potential to inform a range of decisions and management approaches for sustainability planning and resilience-building, especially when used in conjunction with surface watershed approaches.
As countries work towards expanding and strengthening their systems of protected and conserved areas, in alignment with Target 3 (the 30×30 target) of the Kunming-Montreal Global Biodiversity Framework, there are several pathways that could reduce groundwater-mediated threats to existing and new protected areas (as well as to designated OECMs, which equally contribute to the 30×30 target). First, formal area-based protections could be expanded to cover critical groundwatersheds, especially in conjunction with the addition of new protected areas. Extending the boundaries of existing protected areas to encompass their groundwatersheds may be more challenging, as protected areas are often surrounded by land uses that would limit protected area expansion and there is very little land that can be considered free of trade-offs 30 . The expansion of formal protected areas and OECMs may in fact be unnecessary if activities in groundwatersheds can be managed through means such as groundwater permitting or other restrictions 31 . For any of these pathways, additional information about the timescale, magnitude and distribution of impacts on protected areas and people is necessary.
The variability of potential human impacts and the social, economic and political differences across regions imply that a diverse portfolio of  approaches is necessary to protect groundwatershed water quality and quantity. Enhanced conservation and management of groundwatersheds could be achieved through adoption or expansion of strategies such as groundwater regulation (for example, well permitting), sustainable water policies (for example, the Sustainable Groundwater Management Act in California, USA), source water protection (for example, Edwards aquifer protection in Texas, USA), hydrological protection zones for wetlands 32 , Indigenous-led land and water management and monitoring (for example, guardian programmes such as in northwestern Australia), conservation or regenerative agriculture (for example, practices that reduce groundwater pumping) and nature-based solutions (for example, invasive species removal for the Greater Cape Town Water Fund in South Africa). Management strategies could be borrowed or adapted from these and other conservation and source water protection approaches, rather than developing entirely new policy or management approaches.
Selecting an appropriate strategy depends on the social, economic and political context, as well as on the type of possible impact, from severe (nearby, large magnitude pumping or contamination) to less impactful (distant or minor land use change). While we have mapped groundwatersheds of protected areas to place greater focus on groundwater in conservation initiatives, it is important to note that many protected areas also have surface watersheds extending beyond their boundaries. Directly comparing groundwatersheds with surface watersheds is non-trivial as important differences exist in the conceptualization and analysis of these two different types of watersheds, and a detailed comparison is beyond our scope. In this study, we included groundwater-dependent wetlands and root zone intersections with the water table to derive outlet features for groundwatersheds, but these features are not typical outlets for surface watersheds. Furthermore, as surface watersheds are nested and hierarchical, their delineation also hinges on the spatial scale of study. For example, the surface watershed for Mangroves National Park in the Democratic Republic of Congo (located at the outlet of the Congo River) could range from a localized sub-basin to the entire Congo Basin, depending on the scale of analysis. Yet, it holds that for effective conservation, approaches must consider both contributing areas of groundwater and surface water flow that extend beyond protected area boundaries and the human impacts on these systems. For more discussion on comparing groundwatersheds to surface watersheds for protected areas, see Supplementary Section 1.2 and Fig. 5.
Our results importantly highlight the connection between groundwater and protected areas and reveal the vulnerability of protected areas to potential groundwater impacts. However, our approach has limitations (Supplementary Section 2). For instance, we used a simplified approach to identify potential GDEs, focused on higher levels of protection, and mapped only the spatial extent but not the timing of human impacts acting on protected areas through groundwater flow. This first-order global analysis is not intended to lead to recommendations for specific protected areas but rather identifies regional trends in these relationships and discusses potential strategies. With more detailed information, our water table-based approach can be applied to smaller, specific areas. Alternatively, numerical models including particle-tracking approaches that are computationally feasible at local scales can provide more detailed information about the full hydrogeological system and can produce critical insights when combined with the groundwatershed concept and motivation introduced here. As governments around the world commit to new protected area targets and other actors make their own conservation commitments, our analysis serves as a reminder that protection neither stops at protected area borders, nor at the ground surface.

Methods
We implemented a simple geospatial methodology using best-available, openly accessible global data (Supplementary Table 1) to map the groundwatersheds of the world's protected areas. A flow chart of our methodology is shown in Supplementary Fig. 6. All analyses in this study were performed at the spatial resolution of 30 arcseconds (~1 km at the equator), matching the resolution of the water table data used in this study. See Supplementary Section 2.3 for a discussion on the implications of using this spatial resolution.

A computationally simple approach to groundwatershed mapping
Groundwatersheds were derived by making minor modifications to the D8 surface watershed delineation method 33 . Whereas surface watersheds are derived using an outlet location (or 'pour point') and a digital elevation model of the land surface, groundwatersheds are derived using potentially multiple outlet locations and the water table surface instead of the land surface. Whereas a surface watershed identifies the contributing area of overland flow to a specified outlet, a groundwatershed identifies the contributing area of local groundwater flow to groundwater-connected features. In this study, these features are GDEs inside protected areas.
Using this water table-driven D8 flow direction algorithm to derive groundwatersheds does not enable representation of nested deeper regional groundwater flow systems. For an extended discussion on our approach and its limitations, see Supplementary Sections 1 and 2. In the following sections, we summarize our methods to identify GDEs and protected areas, which are combined to derive groundwatershed outlet locations.

Water table
Our analysis is based on a global depth to water table dataset 34 . This dataset contains mean monthly water table depths averaged over a 2004-2013 model run. To generate a mean annual water table depth layer, we calculated the weighted mean of the monthly mean water table depths using the number of days in each month as weights. As we required water table elevations to derive flow direction, we converted  water table depth to water table elevation by subtracting water table  depth from the land surface elevation. We used mean monthly water table elevations in our derivation of GDEs and in our groundwatershed uncertainty analysis, and we used the mean annual water table elevation in our core groundwatershed delineation.

GDEs
Although we mapped the groundwatersheds of the world's protected areas, we did not map groundwatersheds using the entire extent of protected areas as outlet features. Rather, we identified and used areas inside the protected areas where there are likely GDEs. To identify areas with likely GDEs, we considered ecosystems reliant on surficial expressions (for example, wetlands and rivers) and subsurface expressions (for example, phreatophytes) of groundwater, but not subterranean (for example, hyporheic or karst ecosystems) 35 . GDEs were mapped globally using an inference-based approach based on the following: (1) the interaction between rooting depths and the depth to the water table (terrestrial GDEs), (2) the presence of groundwater-dependent wetlands (lentic aquatic GDEs) and (3) surface water systems interconnected with groundwater (lotic aquatic GDEs) systems. Together, these interactions connect groundwater to terrestrial and aquatic ecosystems and are represented by available global data. Our process to identify these interactions is summarized in the numbered paragraphs below.
(1) To identify likely terrestrial GDEs, we considered the relationship between rooting depth and depth to the water table. We identified grid cells where root systems are probably sourcing groundwater by comparing mean monthly depths to the water table with the depth of the root zone 34 . Any grid cell in which the root zone intersects the water table for at least one month per year is identified as a terrestrial GDE.

Analysis
https://doi.org/10.1038/s41893-023-01086-9 (2) To identify likely aquatic GDEs, we considered multiple forms of groundwater-surface water interactions and classified aquatic GDEs as either lotic or lentic systems. To identify lentic systems, we used existing maps of groundwater-dependent wetlands 36 and lake extents 37 . (3) To identify lotic (riverine) aquatic GDEs, we used a network of perennial rivers 38 . Although not all rivers and surface water bodies depend on groundwater discharge, global data availability did not permit the consideration of hydrologically disconnected surface water bodies. However, our use of only perennial river reaches minimizes this impact. We also did not remove losing river and stream reaches as surrounding water table levels regulate the hydraulic gradient across groundwater-surface water interactions.
Losing stream reaches are reflected in our analysis by lower surrounding water table levels and thus do not receive an associated contributing groundwatershed beyond the GDE grid cell(s). In particular, intermittent rivers with a seasonal connection between the groundwater and surface water system (that is, gaining during the wet season, losing during the dry season) 39 may be sensitive to changes in seasonal groundwater levels, but may have been missed in our analysis that focuses on mean annual conditions.
We then combined these three GDE types (terrestrial, lotic and lentic) into a single GDE map. Among these GDEs, those that are located inside protected areas (see below) were used as outlet features in our groundwatershed delineation.

Protected areas
From the World Database on Protected Areas 40 , we created two groups of IUCN terrestrial protected areas: those with high levels of protection that restrict human activity inside their boundaries, and those with lower levels of protection that are more permissive of human activity. The protected area classes we considered as highly protective are: Ia (Strict Nature Reserve), Ib (Wilderness Area), II (National Park), III (National Monument or Feature), as well as protected areas with 'Not Reported' or 'Not Assigned' categories. We included 'Not Reported' and 'Not Assigned' protected areas in this high protection grouping as we found these categories to be more prevalent in countries with lower levels of development where reporting of protected areas may be less comprehensive. By including these categories, we retained a greater global coverage across the protected areas dataset. The remaining protected area categories, IV (Habitat/Species Management Area), V (Protected Landscape/Seascape), VI (Protected area with sustainable use of natural resources) and 'Not Applicable', were grouped together and were considered as having lower levels of protection.
We rasterized both sets of protected areas to our operating resolution, including all grid cells touching a protected area. As the spatial resolution of our analysis is 30 arcseconds (~1 km), we filtered out any protected areas with a reported surface area less than 1 km 2 before rasterization. As protected areas can overlap or border one another, we subsequently identified all spatially contiguous protected areas when representing the protected areas in a binary map at 30 arcseconds, and when applying an 8-connectedness criterion for spatial contiguity.
This set of spatially contiguous protected areas is the protected area set we used as the basis for all calculated metrics (that is, the RGS and the UGR) and to report summary statistics. Using this spatially contiguous but flattened representation of protected areas enabled a more streamlined approach to handle and report global protected area results. These contiguous protected areas differ in total count from the original protected area dataset from which they are derived (Supplementary Table 2).

Groundwatershed delineation
Our groundwatershed delineation process followed conventional watershed delineation approaches that generate a flow direction raster which is used to derive watersheds for specified features. We did not apply hydrological preconditioning steps to the water table surface,  including the removal of depressions as depressions in the water table  represent local water table gradients which we sought to represent in our study. The flow direction raster was generated using the D8 flow direction algorithm, which can represent 8 possible flow directions to adjacent cells according to the direction of the steepest water table gradient. Although the D8 algorithm has known limitations, such as generating parallel flow paths and poorly depicting watersheds in coastal and endorheic basins for example 41,42 , it remains a common, simple, deterministic and widely used approach to derive flow direction. Improving the sophistication of our flow direction derivation may be unwarranted as our analysis was performed at a spatial resolution (30 arcseconds) that is much coarser than conventional watershed-specific delineation studies.
Once the flow direction raster was generated, groundwatersheds were delineated for each GDE cell inside a protected area. Subsequently, groundwatersheds for individual GDE grid cells were grouped across all GDEs found in each contiguous protected area. To avoid double-counting of groundwatershed area, we assigned a single groundwatershed per protected area even if groundwatershed extents may overlap between protected areas. This is possible when a protected area is found inside the groundwatershed of another protected area ( Supplementary Fig. 7). In these cases, the groundwatershed area for the nested protected area is assigned to this protected area, while the remaining groundwatershed area is assigned as the groundwatershed for the downgradient protected area. For a discussion on how this methodological decision affects our calculated summary metrics (RGS and UGR), see Supplementary Section 1.3. To generate flow direction rasters and delineate groundwatersheds, we used the 'D8Pointer' and 'Watershed' functions in the Hydrological Analysis toolbox of the open source geospatial platform Whitebox Geospatial 43 .

RGS and UGR
Once groundwatersheds were delineated for each contiguous protected area, we subsequently evaluated the two metrics (RGS and UGR) developed to understand regional patterns in groundwatersheds. RGS was calculated by dividing the surface area of each groundwatershed by the surface area of GDEs inside the groundwatershed. Importantly, we also considered GDE surface area that exists outside of the protected area but is inside the groundwatershed as stopping at the protected area boundary introduces a social influence on the ecohydrological metric. The UGR was calculated by dividing the surface area of the groundwatershed that lies outside of the protected area by the total surface area of the groundwatershed. These metrics are summarized in Supplementary Table 3.

Uncertainty analysis
As groundwatersheds are dynamic (that is, fluctuate with the water table), we performed an uncertainty analysis to quantify how groundwatershed extents change throughout a typical year. To accomplish this, we repeated our groundwatershed delineation process for mean monthly water table depths and evaluated the variability in total groundwatershed size throughout a year. The results of this uncertainty analysis are included in Supplementary Section 1.4 and Figs. 8 and 9.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
This study uses data from multiple open-access datasets. Source data are documented in Supplementary Table 1 and can be downloaded from the persistent web-links provided. Data produced in this study, including the GDE map, groundwatershed extents and protected area Analysis https://doi.org/10.1038/s41893-023-01086-9 metrics have been deposited on Borealis, the Canadian Dataverse Repository (https://doi.org/10.5683/SP3/P3OU3A).