1 Introduction
Warming associated with anthropogenic climate change has led to a doubling in the frequency of extreme hydroclimate events in the Colorado River Basin since 2010, including droughts, heatwaves, and floods (Bennett et al., 2021). Since 2000, the Colorado River Basin has experienced a historic drought (Meko et al., 2022; Williams et al., 2022). Over this time period, streamflow in the Colorado River has declined by 19% (3.6E9 m3) relative to the 1906-1999 average (18.7E10 m3); (Hogan and Lundquist, 2024; Udall and Overpeck, 2017). Rapid population growth in the Southwestern US, including Arizona, is increasing the demands on already strained water supplies in the State. Reductions in streamflow have increased reliance on groundwater pumping, resulting in groundwater declines across Arizona (Tadych et al., 2024). Analyses of regional gravity data have suggested that the rate of groundwater loss in the Colorado River Basin may far exceed the depletion rate of Lake Powell and Lake Mead and that groundwater may account for a more significant portion of water use than previously anticipated (Castle et al., 2014).
Warming temperatures have tripled the frequency and quadrupled the size of wildfires since 2000 (Iglesias et al., 2022). The risk of catastrophic wildfires is increasing in Western forests–an emerging driver of runoff change that will increase the impact on the water supply (Williams et al., 2022). Forest structure in Northern Arizona has changed significantly post-Euro-American settlement. Many forests are overstocked relative to pre-settlement conditions due to grazing, logging, and wildfire exclusion (Covington and Moore, 1994; Friederici, 2013). These changes have increased the risk of catastrophic wildfires (Allen et al., 2002) post-fire flooding, and debris flows (Youberg et al., 2019). Rising temperatures and droughts have contributed to extensive tree mortality from disease and insect infestation, making forested areas more vulnerable to wildfires and altering natural fire regimes, leading to hotter and more destructive wildfires (Berner et al., 2017; Wasserman and Mueller, 2023).
Forest treatments such as mechanical thinning and prescribed burning can reduce the risk of catastrophic wildfires and alter forests’ hydrologic cycles (Allen et al., 2002; Del Campo et al., 2022). For example, forest thinning in Arizona has been associated with increased snow cover days (Belmonte et al., 2021; Donager et al., 2021; Sankey et al., 2015), greater soil moisture (Belmonte et al., 2022; Sankey and Tatum, 2022), and greater forest canopy moisture (Sankey et al., 2021).
Numerous studies have linked forest treatments to increased water yields in semi-arid forests and have emphasized the role of forest thinning in improving hydrologic services and increasing water availability (Baker, 1986; Bosch and Hewlett, 1982; Gottfried, 1991; Hibbert, 1979; Moreno et al., 2015; O’Donnell et al., 2018; Simonit et al., 2015; Smerdon et al., 2009; Wyatt et al., 2015; Wyatt, 2013a; Zou et al., 2010). While the connection between forest treatment and water yield is well documented, the response of forests to treatments is complex and non-linear and differs across forest types, with treatment level, and along aspect and elevational gradients (Biederman et al., 2022; Del Campo et al., 2022; Hibbert, 1979; Moore and Wondzell, 2005; Zou et al., 2010).
Regardless of the potential for increased water yield, the enhancement of groundwater recharge rarely, if ever, ranks among the primary motivations for forest treatment, even among projects with the stated goal of improving watershed health (Allen et al., 2002; Filoso et al., 2017; Friederici, 2013; O’Donnell, 2016; Stanturf et al., 2014). Responses from surveys of public willingness to pay for forest restoration to enhance groundwater recharge were the lowest among the environmental services discussed in the study, behind critical habitat protections, surface water quality improvements, and the preservation of cultural values(Mueller et al., 2019; Soder et al., 2022). This suggests that groundwater recharge may be considered the least valued among other hydrologic services when considered as stand-alone services and may indicate that the general public is unaware of the complex interdependence between groundwater and surface water (Mueller et al., 2019). Groundwater is difficult to measure and visualize, so it is often undervalued and misunderstood (Campana, 2014; Megdal et al., 2014).
Forest health and water security are intimately linked. About 66% of the water supply in 11 western states comes from forested lands (Brown et al., 2005). Despite this, land, surface water, and groundwater are still managed separately. The Western Water Network (WWN) identified the need to promote better collaboration between researchers, managers, educators, industry, and stakeholders across the American West (Hansen et al., 2024).
This study uses an ecohydrological lens to link insights from multiple fields to examine forest restoration to groundwater recharge enhancement, and aims to identify potential recharge zones. Specifically, we demonstrate a suitability mapping approach using multi-criteria decision analysis (MCDA) to identify areas where forest thinning may enhance recharge as an additional benefit. Suitability maps like these may complement (or supplement) existing frameworks for jointly prioritizing landscape-scale forest management and managing lands and water.
MCDA has long been used in forest management for strategic (long-term, organizational) and tactical (annual) planning (Blagojević et al., 2019; Mendoza and Martins, 2006). MCDA in combination with GIS, hereafter GIS-MCDA, is a methodology that combines geographical data and value judgments to support decision-making (Malczewski, 2006). Suitability mapping, particularly GIS-MCDA, has been widely used to map potential recharge zones and areas suitable for Managed Aquifer Recharge (MAR)(Fathi et al., 2021; Rahman et al., 2012; Raja Shekar and Mathew, 2023; Shaban et al., 2006). However, this study is the first to examine where forest thinning may be sited to enhance recharge at a landscape scale.
1.1 Study Area
Our study area includes all of the HUC8 (Hydrologic Unit Code 8) watersheds of Arizona, extending into 5 other U.S. states. The portions of HUC8 watersheds in Mexico were excluded due to a lack of data. However, within that study area, we focus on the ponderosa pine (PIPO) woodland type, which covers approximately 1.8 million ha (1,835,530.65) of this area, comprising 28% of the forested area within its bounds. Most PIPO forest is found along the Mogollon Rim, around the San Francisco Peaks, and on the Kaibab Plateau–all within the Colorado Plateau physiographic region at elevations between 2000 m and 2400 m. Smaller isolated patches of PIPO forest can be found on mountain ranges and plateaus throughout Arizona at similar elevations (Figure 1 ).
The Mogollon Rim is a topographic feature forming the southern edge of the Colorado Plateau. It contains most of the state’s PIPO forests and has been identified as an essential groundwater recharge area for regional aquifers (Parker et al., 2005). Of the estimated 2.1 billion m3 (174,000 acre-feet) of precipitation that falls on the Mogollon Rim annually, about 8% is estimated to recharge the regional groundwater aquifers(Parker et al., 2005) Similarly, The Kaibab Plateau, an uplifted region north of the Grand Canyon with large stands of PIPO, is a critically important recharge area for springs in the Grand Canyon which provide the sole source of drinking water for over 6 million annual visitors and the residents of Grand Canyon National Park(Tobin et al., 2018). The Kaibab Plateau lacks perennial surface streams and is drained rapidly by sinkholes, faults, and fractures or more slowly through diffuse infiltration with little streamflow leaving the plateau (Huntoon, 1974; Jones et al., 2018).
Most of the PIPO forests within Arizona are within the Kaibab, Coconino, Apache-Sitgreaves, Tonto National Forests, The Fort Apache Reservation, and Navajo Nation. Precipitation in Arizona is bi-modal, with wet winters and a late-summer monsoon season. As much as 60% of precipitation at elevations above 2000 m comes in the form of snow, and snow is responsible for 80 - 95% of streamflow and much of the recharge (Baker, 2013; Earman et al., 2006; Eastoe, 2023). Therefore, this research focuses on high-elevation (>= 2000 m) PIPO forests, which receive a high proportion of annual precipitation as snow. Of the 1.8 million ha of PIPO forests in our study area, about 43% or 790,000 ha are underlain by lithologies that may include karst or pseudo-karst features, including carbonates (24%), evaporates (6%), and fractured volcanics (13%). These areas merit extra attention due to enhanced infiltration potential through fractures, faults, lava tubes, or sinkholes.
1.2 Literature Review and Criteria
This research aims to identify areas where mechanical thinning is most likely to enhance groundwater recharge. We based our suitability criteria on literature primarily from regional studies, as they likely provide usful information about hydrologic response to thinning in Arizona’s forests (Wyatt, 2013b). Landscape-scale forest restoration efforts have been planned or implemented across much of Arizona. For example, the Four Forest Restoration Initiative (4FRI) includes a goal of restoring over 1 million hectares of Arizona’s forests (Schultz et al., 2012). A synthesis study of 4FRI treatments found that thinned and burned forests have significantly greater total ecosystem moisture and are thus more resilient to drought and wildfire (Sankey et al., 2021; Sankey and Tatum, 2022). Thinned forests also have more snow and soil moisture (Belmonte et al., 2022; O’Donnell et al., 2021).
Snow and its persistence appear particularly important for recharge in Arizona’s semi-arid and high-elevation conifer forests. An analysis of groundwater coming from springs along the Mogollon rim found a strong relationship between snow persistence and the duration of spring discharge (Donovan et al., 2022). Isotopic groundwater analyses have revealed that Northern Arizona’s recharge is dominated by winter precipitation from altitudes above 1500 m(Blasch et al., 2006; Earman et al., 2006; Eastoe, 2023; Springer et al., 2017). Studies of plant water use have found that larger PIPO trees primarily utilize deeper soil moisture from winter precipitation, while understory vegetation and smaller trees utilize shallow soil water from monsoonal storms(Kerhoulas et al., 2023; Kerhoulas et al., 2013).
Thinning in PIPO forests is expected to increase available water in two ways: 1) reducing transpiration of deep soil water by overstory vegetation and 2) reducing sublimation of intercepted snowfall, which could increase below canopy snow depth and persistence. Thinning in semi-arid forested watersheds can significantly alter snowmelt timing (Dwivedi et al., 2024a). Reduced forest cover can delay snowmelt at colder sites, higher elevations, and northern aspects or advance snowmelt through increasing sub-canopy solar radiation and wind, particularly at lower elevations, southern aspects, and warmer sites (Biederman et al., 2015; Dwivedi et al., 2024a).
Total annual precipitation also appears to influence whether thinning reduces or enhances water yield at particular sites. In the Colorado River Basin, there appears to be a precipitation threshold for hydrologic responses. Below ~500 mm of annual precipitation, precipitation and runoff may become decoupled (Biederman et al., 2022), likely because below ~500 mm, most precipitation is evaporated regardless of forest condition. Semi-arid forests with high inter-annual precipitation variability may show effects in wet years when precipitation is greater than ~500 mm or in snow-dominant watersheds (Adams et al., 2012; Carroll et al., 2016).
Reductions in canopy cover can reduce transpiration losses and decrease soil moisture stress (Frank et al., 2019; Knowles et al., 2023; O’Donnell et al., 2021; Sankey and Tatum, 2022) or increase surface evaporation and sublimation depending on the elevation, aspect, and the amount of canopy removed (Biederman et al., 2014; Dwivedi et al., 2024b; Harpold et al., 2013). Canopy cover between 25% and 40% appears optimal for net snow accumulation at continental mid-latitude sites (Veatch et al., 2009), though this is dependent on elevation and aspect, with lower forest densities being optimal areas with shorter snowpack duration (such as low elevations) and higher forest densities being optimal for areas with longer snowpack duration (Broxton et al., 2020). UAV-based studies within the study area have shown that 24 - 35% canopy cover was optimal for snow persistence. (Belmonte et al., 2021; Donager et al., 2021; Sankey et al., 2015). However, water yield enhancement from thinning likely requires at least a 20% reduction in canopy cover (Adams et al., 2012). Hydrologic models in northern Arizona watersheds suggest that reductions in basal area of at least 30% can increase groundwater recharge by up to 15% (Wyatt et al., 2015).
Research on the effect of thinning to below-ground hydrological processes in semi-arid Mediterranean forests found that sites with high antecedent soil moisture had the highest response, with drainage to deeper soil layers increasing by 50 mm/year relative to control sites(Del Campo et al., 2019). Soil type (texture, composition, field capacity, and permeability) and topography affect soil moisture conditions. Slope, aspect, landscape convexity or concavity, and geomorphons (or landform types) such as pits, flats, hills, valleys, ridges, and hollows strongly affect the accumulation and residence time of overland flow and shallow soil water (Jasiewicz and Stepinski, 2013; Parker, 1982). In dry soils, capillary forces generally dominate over gravitational forces, and areas with higher flow accumulation and residence time are more likely to reach field capacity, leading to more percolation (Woessner and Poeter, 2020).
Another important aspect affecting recharge enhancement potential is the underlying lithology, which controls the subsurface flow of groundwater. Arizona currently lacks detailed hydrogeologic mapping. However, representative ranges of primary (matrix) permeability and porosity values for rock types are available based on 1:1,000,000 and 1:500,000 geologic mapping data(Freeze and Cherry, 1979; Huscroft et al., 2018; Woessner and Poeter, 2020). Secondary permeability and porosity, which involve flow in faults and fractures, can be estimated by measuring the density of linear surface features, or ‘lineaments,’ which often correspond to underlying geologic structures (O’leary et al., 1976; Sander, 2007). Lineaments have been used to identify potential groundwater supplies, and studies have found higher well yields along lineaments (Sander, 2007). In some crystalline bedrock, groundwater flow occurs exclusively in faults and fractures (Meijerink et al., 2007; Nyborg et al., 2007). Orientations of lineaments can help determine the anisotropy of the fracture environment(Meijerink et al., 2007). Caves and sinkholes (or dolines) often form along and at the intersection of faults and fractures in karst terrains (Ford and Williams, 2007). Karst and pseudo-karst terrains can have enhanced infiltration potential far beyond the hydraulic conductivity of the rock matrix (Krešić, 2013). For example, the Edwards Aquifer in Texas was found to have permeability values that vary by up to nine orders of magnitude (Halihan and Wicks, 1998). The presence or absence of karst or pseudo-karst in concert with lineament density is a critical criteria in our suitability analysis.
In summary, suitable sites where thinning enhances recharge would most likely occur in areas with significant snowfall, maximum annual precipitation >= 500 mm, higher antecedent soil moisture, NE aspects, in valleys or flat areas where a 20% reduction in canopy cover and a 30% reduction in basal area would result in a thinned canopy cover of between 25 and 35%. Suitability would be lowest in sites with minimal snowfall, SW aspects, lower elevations, ridge tops, or steep slopes, where thinning might reduce canopy cover to below 24%. Soil hydraulic conductivity, flow accumulation, and topographic metrics can reveal the areas likely to have higher antecedent soil moisture. Lithology, primary permeability and porosity, and lineament density can provide insight into where rapid subsurface infiltration is possible.
2 Methods
The MCDA-GIS suitability mapping process generally involves defining the objective, screening suitable areas, classifying thematic layers or criteria, standardizing, weighing, and sensitivity analysis (Malczewski, 2006; Mendoza and Martins, 2006). In this case, the objective is to determine the relative suitability of areas to achieve enhanced recharge through forest treatment. The criteria were selected based on a review of the literature in Section 1.2 . Standardization was achieved by scaling variables and final suitability on a scale of 1 - 10, from lowest to highest suitability, to facilitate straightforward interpretation. Weighting the relative contribution of 12 different variables, some likely exhibiting strong collinearity, through pairwise comparison would be complex. To address this challenge, similar variables were consildated into indices, which were subsequently weighted against each other, streamlining analysis and reducing redundancy.
2.1 Suitability Criteria
We used ArcGIS Pro 3.4.0 to create a weighted suitability model consisting of the 12 thematic data layers (Table 1). These data layers were re-scaled, weighted, and combined into four final thematic layers: (1) Snowfall fraction (SF), (2) Vegetation Density Index (VDI), (3) Soil Moisture and Infiltration Index (SMII), and (4) Subsurface Infiltration Index (SbII : Table 1) weights within and between indices were assigned using pairwise comparisons, and applied using weighted linear sums.
| Thematic Layer | Intermediate Index | Scaled Variable | Data Layers | Source |
|---|---|---|---|---|
| SF | - | sSF | Snowfall Fraction | UA SnowData (800m) |
| VDI | - | sBA | Basal Area | TreeMap 2016 (30m) |
| VDI | - | sCC | Canopy Cover | NLCD 2021 (30m) |
| SMII | TRMI | sApt | Aspect | Derived from 1-Arcsecond DEM |
| SMII | TRMI | sGeo | Geomorphon | Derived from 1-Arcsecond DEM |
| SMII | TRMI | sSlop | Slope | Derived from 1-Arcsecond DEM |
| SMII | TRMI | sTPI | Topographic Position Index | Derived from 1-Arcsecond DEM |
| SMII | - | sSoilK | Soil Saturated Hydraulic Conductivity | ggNATSGO (30m) |
| SbII | - | sPo | Matrix Porosity | GLHYMPS v2 |
| SbII | - | sPm | Matrix Permeability | GLHYMPS v2 |
| SbII | - | sLD | Lineament Density | Derived from 1-Arcsecond DEM |
| SbII | - | sPK | Potential Karst or Pseudo-Karst | Karst in the United States (USGS) |
2.1.1 Snowfall Fraction (SF)
Snowfall Fraction (SF) is the percent of annual precipitation that comes in the form of snow. The SF metric was created using accumulated snowfall calculated from the University of Arizona snowpack (UA Snow) data (Zeng et al., 2018) by summing positive increments of snow water equivalent (SWE; to estimate snowfall) and dividing by accumulated precipitation from PRISM (Daly et al., 2008). The UA Snow data, itself is based on the interpolation of SWE and snow depthh measurements from ground stations against a background field of snowpack estimates created using PRISM precipitation and temperature data, and importantly, it contains an algorithm that separates rainfall and snowfall based on a spatiotemporally variable rain-snow transition temperature that varies based on station whether stations record snow accumulation (Broxton et al., 2016). Note that in this study, we use the 800 m version of the UA data (Broxton et al., 2024).
Snowfall is important for recharge at elevations over 2000 m in Arizona (Baker, 2013; Earman et al., 2006; Eastoe, 2023). Snow dominance and snowmelt duration are critical factors in whether thinning enhances recharge (Adams et al., 2012; Carroll et al., 2016; Saksa et al., 2017). Therefore, our suitability mapping prioritizes areas with more snowfall relative to rainfall (higher SF). SF was scaled linearly from 1 to 10 to create a scaled snowfall fraction (sSF) thematic layer.
2.1.2 Vegetation Density Index (VDI)
The Vegetation Density Index (VDI) is created from two data layers, canopy cover (CC) and basal area (BA).
2.1.2.1 Canopy Cover (CC)
Canopy Cover was obtained from the 2021 National Land Cover Database (NLCD), which estimates total canopy cover at a 30m resolution. CC values ranged from 0% - 86%. The variable NLCD TCC (total canopy cover) was reclassified from 1 to 10 based on suitability for thinning to enhance recharge. Forest canopy cover<= 12.5% was reclassified as equal to 1, scaling suitability linearly up to CC >= 28.75% equal to 10 (see appendix: Figure 12;Figure 11 ). These threshold values were determined by considering a minimum 20% reduction in CC without dropping below the pre-settlement minimum CC of 10% (Battaglia and Shepperd, 2007). Studies examining the relationship between CC and snow retention in PIPO forests found that 23% to 35% is the ideal CC range for retaining snow(Belmonte et al., 2021; Donager et al., 2021; Sankey et al., 2015; Veatch et al., 2009). Areas with CC currently exceeding 28.75%, were considered ideal candidates for thinning. The scaled canopy cover (sCC) thematic layer resulting from this reclassification function was used as an input layer for further analysis.
2.1.2.2 Basal Area (BA)
We extracted BA estimates from the TreeMap 2016 (Riley et al., 2022) CONUS dataset with a resolution of 30 meters. BA was rescaled similarly to canopy cover assuming a 30% reduction required to increase water yield, and a pre-settlement estimate of between 9.2 and 18 m2ha-1 for PIPO (Battaglia and Shepperd, 2007). BA was therefore rescaled linearly from <= 13.14 m2ha-1 equal to 1, up to >= 25.71 m2ha-1 equal to 10, see appendix (Figure 11;Figure 12). The scaled basal area (sBA) layer resulting from this reclassification function was used as an input layer for further analysis.
2.1.2.3 Creation of VDI
BA and CC are often highly correlated but represent different aspects of forest structure. Therefore, they were combined into a single index. Vegetation density index (VDI) is the weighed sum of sBA, and sCC was combined with weights of 0.6667 and 0.3333, respectively, to create the vegetation density index (VDI)(Equation 1). We used pairwise comparisons to compare the importance of the two variables and ultimately assigned a higher weight to the basal area. This reflects the fact that reductions in basal area are more likely correlated with reductions in transpiration, while canopy cover primarily affects precipitation partitioning and sub-canopy solar radiation.
\[ VDI =sBA* 0.6667 + sCC * 0.3333 \tag{1}\]
2.1.3 Soil Moisture and Infiltration Index (SMII)
Soil Moisture and Infiltration Index (SMII) consists of two data layers, Soil Hydraulic Conductivity (SoilK) and Topographi Relative Moisture Index (TRMI), which includes several topographic paramters derived from a digital elevation model (DEM).
2.1.3.1 Soil Hydraulic Conductivity (SoilK)
SoilK is used to estimate the max infiltration rate of the soil profile. SoilK comes from the gridded National Soil Geographic Database produced by the USDA-NRCS. The gNATSGO Soil hydraulic conductivity layer estimates the rate at which water moves through the pores of saturated soil based on point field measurements of soil texture, structure, and porosity. Data was extracted from the gNATSGO database for the AZHUC8 study area. The depth-weighted average method was used to calculate values for soil layers from 0 - 200 cm for each 30 by 30 m pixel, using the Soil Data Development Toolbox in ArcMap 10.8.3. Values are reported in micrometers per second and broken into the standard 6 classes from very low (0.00 - 0.01) to Very High (100 - 705). These standard classes were assigned suitability values between 1 and 10 to create a scaled SoilK layer (sSoilK)(Figure 2). See Table 2 in the appendix for more details.
2.1.3.2 Topographic Relative Moisture Index (TRMI)
TRMI incorporates several topographic parameters that influence moisture dynamics, including slope gradient, aspect, relative elevation (or topographic position), and landscape convexity or concavity (Parker, 1982) (see appendix: Figure 15). TRMI was calculated as per Parker (1982) by reclassifying and then summing each of the input layers; slope (degrees) (1-10), slope configuration or topographic position index (TPI)(1-10) (De Reu et al., 2013), geomorphon (ArcGIS Pro 3.4.0 Spatial Analyst Toolbox–Geomorphon Landform Tool) (1-20), and aspect (degrees azimuth) (1-20). The resulting TRMI values were between 4 and 60, subsequently reclassified into 10 classes with equal interval breaks to create a scaled TRMI layer (sTRMI). Additional details regarding the geomorphon layer and specific reclassification functions used for each input layer can be found in the appendix (Table 3;Table 4;Table 5).
2.1.3.3 Creation of SMII
The Soil Moisture and Infiltration Index (SMII) represents where water accumulates on the landscape and where it is likely to infiltrate through the top 2 meters of the soil profile. This index is a weighted sum of the TRMI and SoilK (Figure 2). We prioritized moisture supply (TRMI) over infiltration (SoilK), reasoning that infiltration matters less in areas with no accumulated water to infiltrate. Hence, we applied a weight of 0.6667 to sTRMI and 0.33333 to sSoilK (Equation 2;see appendix:Figure 13).
\[ SMII = sTRMI * 0.6667 + sSoilK * 0.3333 \tag{2}\]
2.1.4 Subsurface Infiltration Index (SbII)
The Subsurface Infiltration Index (SbII) was created by combining layers of Matrix Pereabolity (Pm) and Porosity (Po), Lineament Density (LD), and Potential karst of Pseudo-karst (PK).
2.1.4.1 Matrix Permeability (Pm) and Porosity (Po)
Pm and Po values were extracted from the GLobal Hydrogeology MaPS 2.0 (GLHYMPS v2) dataset (Huscroft et al., 2018). Our study area contained 11,057 polygons with values for saturated log-permeability ranging from -16 to -10 and porosity values ranging from 0.01 - 0.28 (Figure 10, see appendix). The attribute values were extracted, the polygons were rasterized at a 30m resolution and then re-scaled to values (1-10). See the appendix for more details on the rescaling of Po and Pm (Table 6).
2.1.4.2 Lineament Density (LD)
Lineaments are linear or curvilinear surface features that may reflect subsurface geologic structures(O’leary et al., 1976). Lineament density(LD) is commonly used in studies evaluating potential recharge zones (Al-Adamat, 2012; Chenini et al., 2010; Chowdhury et al., 2010). Lineaments were extracted from U.S. Geologic Survey, 2019 3D Elevation Program 1-arc-second tiles using the LINE algoithm in Catalyst 3.0.2. the process is described in detail in the appendix. The ouput of the LINE tool was exported as a shapefile and analysed using the Line Density Tool in ArcGIS Pro to calculate the density of lineaments within a circular neighborhood with a radius of 1km, providing lineament density values between 0.001 and 5.55 km/km2. While the presence of lineaments likely enhances the secondary Pm and Po, the absence of lineaments does not necessarily reduce overall recharge suitability. Therefore, we scaled lineament density between 5 and 10 for the thematic layer scaled lineament density (sLD) (see appendix: Figure 16); this is consistent with other GIS-MCDA studies evaluating recharge potential (Shaban et al., 2006).
2.1.4.3 Potential Karst or Pseudo-Karst (PK)
Shapefiles of potential karst or pseudo-karst lithology (PK) were sourced from the Karst in the United States digital map compilation and database (Weary and Doctor, 2014). The polygons were clipped to the AZHU8 study area and converted to a 30m raster, and assigned values based on their potential for enhanced recharge through faults, fractures, sinkholes, or caves. The classes were non-karst = 5 (neutral), volcanics with potential karst or pseudo-karst (such as lava tubes) = 7, evaporates at or near the surface = 8, and carbonates at or near the surface = 10, resulting in a scaled potential karst layer (sPK).(Figure 3)
2.1.4.4 Creation of SbII
SbII is a weighted linear sum of sPm, sPo, sLD, and sPK, with weights 0.4, 0.25, 0.15, and 0.2, respectively; the equation is shown below (Equation 3)(Figure 3; see appendix: Figure 14). The SbII can be seen as an estimate of the infiltration potential of the subsurface lithology excluding soils (which are included in the SMII index). Due to the resolution of the geologic mapping that produced the GLHYMPS v2 dataset and general lack of detailed geologic mapping in much of Arizona, focused recharge in stream channels is not addressed by this index, and it can be seen as a conservative estimate of recharge potential state-wide.
\[ SbII = sPm * 0.4 + sPo * 0.25 + sLD * 0.15 + sPK * 0.2 \tag{3}\]
2.2 Overall Suitability Weighting
The four final thematic layers: sSF, VDI, SMII, SbII were ranked in order of importance and assigned weights of 0.2, 0.2, 0.4, and 0.2 respectively and combined using a weighted sum resulting in the final suitability values for thinning to enhance recharge (Equation 4). The relative proportion, or percent contribution of each thematic layer to overall suitability is shown in Figure 4
\[ Final Suitability = 0.2 * sSF + 0.2 * VDI + 0.4 * SMII + 0.2 * SbII \tag{4}\]
2.3 Sensitivity Analysis
A sensitivity analysis was conducted on a subset of the data using a one-at-a-time (OAT) analysis, where each parameter was systematically removed from the VDI, SMII, SbII and final suitability. Also, for each of the indices and final suitability, the applied weights were shuffled to determine the effect of differing the weights on mean pixel values and the percentage of pixels with values over 5, 6, and 7, representing areas of moderate, high, and very high suitability, respectively, for thinning to enhance recharge. Altogether, 32 different weighting simulations and 14 removal simulations were conducted on the various indices and final suitability(see appendix tables 8-14). For each index and final suitability, the resulting indices from weight shuffling and variable removal were stacked. The Coefficient of Variation (CoV), mean, range and standard deviation were calculated for all pixels (x,y) and plotted to indicate where variability in suitability numbers manifested spatially and to understand the uncertainty in index values (see appendix: Section 8.3).
3 Results
3.0.1 Overall Suitability
The study area contained 1.8 million ha of PIPO forest. Approximately 44% (807,000 ha) of the PIPO forest was deemed suitable (> 6) and about 10% (182,000 ha) were deemed highly suitable (>7) for thinning to enhance recharge (Figure 5; Figure 9 in appendix). The remainder of the PIPO forest is moderate suitability, low suitability, or not suitable for thinning to enhance recharge 41% (762,000 ha). About 48% of the area was marginally suitable or marginally unsuitable (4 - 6), and about 8% (141,000 ha) was found to be unsuitable for thinning to enhance recharge.
4FRI has planned (through Environmental Impact Statements) or implemented thinning on about 180,000 ha of PIPO forest (Stewart and Williams, 2025). Within those areas, about 16% or 30,000 ha are highly suitable for thinning to enhance groundwater recharge (or have suitability values >= 7). Within the 4FRI planning area, there are still about 540,000 ha, which are not part of areas planned for thinning in one of the two currently existing environmental impact surveys. About 12% or 65,000 ha of that area are highly suitable for thinning to enhance groundwater recharge (or have suitability values >= 7) (Figure 6). This information can help forest managers decide which areas to include in the next phase of planning and provide information on the co-benefits of thinning for recharge enhancement as well as forest health and fire risk mitigation.
Large areas of high-suitability PIPO forest can also be found on the Fort Apache Reservation along the Mogollon Rim, and in Navajo Nation on the Defiance Plateau and in the Chuska Mountains, which are located to the north of the Mogollon Rim near the border between Arizona and New Mexico(Figure 7). Groundwater security is particularly important in native communities. Over 70,000 members of the Navajo Nation lack running water and haul water, often from wells (Credo et al., 2019; Jones and Ingram, 2022).
3.0.2 Sensitivity Analysis Results
3.0.2.1 VDI Sensitivity
Shuffling the weights of sCC and sBA resulted in a relatively small (0.3) change in the mean value, of VDI, and a relatively small change in the percent of pixels which were highly suitable (>7). Suggesting that either approach likely would not change the results of the final suitability very much. There were some marginal changes though the percent of moderately suitable (>5) and moderately high suitable pixels (>6). Some differences are also apparent when looking at the pixel-wisee variability, where large patches of forest have relatively high CoV values. After taking a closer look, these are areas with high canopy cover but low basal area, this is reassuring because it shows that these two variables are non-redundant (see appendix:Table 8; Figure 17). The OAT removal of each variable in VDI resulting in slightly larger changes to overall suitability than the weight swapping, removing canopy cover increased the mean value, suggesting it might be the limiting factor(see appendix:Figure 18; Table 9).
3.0.2.2 SMII Sensitivity
SMII values increased when sTRMI had a lower weight relative to sSoilK and when sTRMI was removed. This suggests that sTRMI is limiting suitability in areas where there is high soil permeability. The inclusion of TRMI is an important aspect of recharge potential, and it is acceptable that its addition reduces suitability in some places (See appendix: Table 10; Figure 19; Table 11; Figure 20).
3.0.2.3 SbII Sensitivity
There were 24 possible weighting combinations for SbII, because sLD and sPK were scaled 5-10 nearly all of the area had values >5. However, lower mean values and a lower percentage of high suitability pixels (>7) were observed when sPo was weighted heavily. The same pattern is found when sPo was removed–higher mean and high suitability values. Suggesting the sPo is the limiting factor. sLD also appears to be a factor limiting suitability as shown by the removal of that variable and resulting lower mean suitability values. High SbII then appears to be driven permeability and Karst lithology, which is what we would expect. The variability between runs shows that these variables are non-redundant (see appendix: Table 12; Figure 21; Table 13; Figure 22).
3.0.2.4 Sensitivity of Overall Suitability
Variability, or uncertainty in final suitability values was systematically examined by removing variables OAT, and shuffling weights (see appendix: Table 14; Figure 23; Table 15; Figure 24). The max pixel-wise CoV when swapping weights was about 0.5 or a 50% change in final suitability values. The high CoV values are mostly found in square blocks, which correspond to the SF raster’s 800m pixels before resampling, suggesting that high variability is being driven by the sSF raster primarily; their scattered distribution suggests that they may be related to noise in the precipitation data. Otherwise, variability appears high along major river channels, likely due to edge effects and having no data for some valley bottoms (Figure 8).
4 Discussion and Conclusion
The sensitivity analysis showed that results after weighting changes and variable removal were relatively consistent and interpretable. Changes in suitability when removing or re-weighting variables are expected if variables are non-redundant. Most of the changes in suitability values for all indices could be explained and are not at odds with the reviewed literature. The overall mean pixel-wise range for suitability values was 1.4; therefore, the suitability values for any given spot are likely uncertain by +/—1.4. Modeling or field studies could help reduce the overall uncertainty.
This research demonstrates a novel application of GIS-MCDA for mapping areas where forest thinning could enhance groundwater recharge. The methodology here is adaptable and could be implemented in other semi-arid forested regions. More work is needed to validate suitability maps using process-based hydrologic models.
Forest thinning is already planned in over 1 million ha of Arizona’s PIPO forest. Thinning is intended to reduce the risk of catastrophic wildfires, protect forest resources, enhance surface water provisioning, and improve habitat. However, we have also demonstrated that thinning will likely enhance groundwater recharge in much of the PIPO forest. This important co-benefit may increase stakeholder buy-in and strengthen justification for management actions.
As water managers across the state struggle to secure new water supplies to meet the demands of population growth and maintain the agricultural economy, insecurity in supplies from the Colorado River means that more and more water users are turning to groundwater to fill the gap. Increased interest in thinning to enhance recharge could encourage land managers to speed-up the pace of restoration and demonstrate the value of thinning to water managers looking to improve water supply resiliency, potentially providing an alternative funding source.
The Forest Service is mandated to manage for multiple uses; however, groundwater recharge is not currently an explicitly managed use. In 2014, the US Forest Service introduced a proposed directive that would mandate the consideration of groundwater resources within its activities. However, push-back from stakeholders and state water management agencies led to the withdrawal of this directive, with the Forest Service vowing to start over and develop ways to include groundwater resources within their planning(Lexology, 2014; U.S. Forest Service, 2014; U.S. House Committee on Agriculture, 2014). In 2024, the President’s Council of Advisers on Science and Technology, as part of the Bipartisan Infrastructure Law and Inflation Reduction Act explicitly recommended acceleration in the development of collaborative databases and tools to address groundwater storage, withdraw and recharge at spatial scales useful to water managers and users(President’s Council of Advisors on Science and Technology et al., 2024). This analysis offers an a new tool to assess potential groundwater recharge at landscape scales through forest management, allowing for groundwater recharge to be managed as a co-benefit of thinning for forest health and wildfire risk reduction.
One of the main challenges in mapping suitability at landscape scales is the lack of high-resolution data. Arizona lacks detailed hydrogeologic mapping throughout much of the state. The State Geologic map within the State Geologic Mapping Compilation contains geologic units and faults mapped to the 1:500,000 scale. However, digitized lineaments, particularly faults, are limited to very large faults state-wide and are almost completely absent in the northeastern part of the state on the Navajo Nation and Colorado Plateau.
Groundwater is and always has been difficult to monitor and measure. However, springs are convenient locations to monitor the effects of land management on groundwater recharge. Though tracer studies are needed to connect specific springs to particular recharge areas, continuous monitoring of springs and the base flow component of perennial spring-fed rivers can inform land managers of the effectiveness of thinning in enhancing recharge and help validate suitability analyses.
5 Acknowledgments
Thank you to the Arizona Board of Regents for funding this research through the Technology and Research Initiative Fund (TRIF) through the Arizona Tri-University Recharge and Water Reliability Project (ATUR). Thank you to all the members of ATUR for your support, guidance, and suggestions. Special thanks to Abdul Moiz for generating hydroclimate datasets used in this analysis.
6 Open research
The data, methods, and code base for this research can be found on github: https://github.com/Ryan3Lima/ATUR-Thinning-to-enhance-recharge , in the Open Sciences Framework Project associated with this work. The data will also be available on hydro-share.
7 Declaration of generative AI and AI-assisted technologies in the writing process
During the preparation of this work, the authors used ChatGPT4o and Research Rabbit to assist in our literature review. ChatGPT4o and Github Copilot were used in developing the Python code for the analysis tools created to combine and weight layers and conduct the sensitivity analyses. The authors reviewed and edited the content as needed and take full responsibility for the content of the published article.
References
8 Appendix
8.1 Additional Figures
8.2 Addtional Methods
8.2.1 VDI
8.2.2 SMII
| Ksat | Standard Class | Assigned suitability score | Present in Study Area |
|---|---|---|---|
| 0.0 - 0.01 | Very Low | 1 | No |
| 0.01 - 0.1 | Low | 2 | No |
| 0.1 - 1 | Moderately Low | 4 | Yes |
| 1 - 10 | Moderately High | 6 | Yes |
| 10 - 100 | High | 8 | Yes |
| 100 - 705 | Very High | 10 | Yes |
| Geomorphon landform | Wetness value (1-20) |
|---|---|
| Pit | 20 |
| Valley | 18 |
| Hollow | 14 |
| Footslope | 12 |
| Flat | 10 |
| Slope | 8 |
| Spur | 6 |
| Shoulder | 4 |
| Ridge | 2 |
| Peak | 1 |
| Steepness (degrees) | Wetness Value (1-10) |
|---|---|
| <3.0 | 10 |
| 3.0 - 5.9 | 9 |
| 6.0 - 8.9 | 8 |
| 9.0 - 11.9 | 7 |
| 12.0 - 14.9 | 6 |
| 15.0 - 17.9 | 5 |
| 18.0 - 20.9 | 4 |
| 21.0 - 23.9 | 3 |
| 24.0 - 26.9 | 2 |
| >27.0 | 1 |
| Aspect ° | Description | Wetness Value |
|---|---|---|
| -1 | Neutral (flat areas) | 10 |
| 0 ° (N) | Wettest (Shaded, least evaporation) | 20 |
| 45 ° (NE) | Moist, receives moderate sun | 16 |
| 90 ° (E) | Neutral (morning sun, less drying) | 10 |
| 135 ° (SE) | Drying increases | 6 |
| 180 ° (S) | Driest (maximum sun exposure) | 1 |
| 225 ° (SW) | Quite dry | 4 |
| 270 ° (West) | Neutral (afternoon sun, retains some moisture) | 10 |
| 315 °(NW) | Moist, cooler | 16 |
| 360 ° (N) | Wettest | 20 |
Topographic Relative Moisture Index (TRMI) was calculated by summing the wetness values for aspect, slope steepness, geomophon, and TPI by summing, resulting in values between 3 and 60 (Equation 5). The TRMI values were then reclassified to a 1-10 scale using (Equation 6).
\[ TRMI = (Aspect + Slope + Geomorphon + TPI) \tag{5}\]
\[ TRMI scaled = \frac{(TRMI - 3)}{((60-3)*9)} \tag{6}\]
8.2.3 SbII
| LogK x 100 (Permeability) | Scaled Permeability (sPm) | Porosity x 100 (%) | Scaled Porosity (sPo) |
|---|---|---|---|
| -1650 to -1595 | 1 | 0 - 2 | 1 |
| -1595 to -1540 | 2 | 2 - 5 | 2 |
| -1540 to -1485 | 3 | 5 - 8 | 3 |
| -1485 to -1430 | 4 | 8 - 10 | 4 |
| -1430 to -1375 | 5 | 10 - 13 | 5 |
| -1375 to -1320 | 6 | 13 - 17 | 6 |
| -1320 to -1265 | 7 | 17 - 21 | 7 |
| -1265 to -1210 | 8 | 21 - 24 | 8 |
| -1210 to -1155 | 9 | 24 - 26 | 9 |
| -1155 to -1052 | 10 | 26 - 28 | 10 |
8.2.3.1 LD
Lineaments were extracted from U.S. Geologic Survey, 2019 3D Elevation Program 1-arc-second tiles using a combination of methods. First, topographic position index (TPI) which estimates neighborhood convexity and concavity was calculated using a circular neighborhood with a radius of 5m The TPI layer was then rendered as a multi-directional hillshade and exported from ArcGIS to Catalyst 3.0.2, where the LINE module was used to detect edges and link lines (Salui, 2018). The default parameters for the LINE function were used (Table 7), and the resulting lineaments were brought back into ArcGIS and cleaned to remove any edge effects. A second lineament analysis was conducted as with the TPI, using a standard deviation focal statistic of elevation with a 5 m radius. This raster was also rendered, exported, and processed as before. These two rasters provided complementary sets of combined lineaments. Due to the resolution and parameters chose, man-made lineaments (roads, railways, and power-lines) were not primarily delineated by the algorithm in either analysis; however buffers around roads were randomly sampled and spot-checked to ensure man-made lineaments were not included.
| Parameters | Description | Units | Default Values | Distance (m) |
|---|---|---|---|---|
| RADI | filter radius | Pixels (30m) | 10 | 300 |
| GTHR | edge gradient threshold | Luminance (0-255) | 100 | NA |
| FTHR | line fitting threshold | Pixels (30m) | 3 | 90 |
| LTHR | curve length threshold | Pixels (30m) | 30 | 900 |
| ATHR | angular difference threshold | degrees | 30 | NA |
| DTHR | linking distance threshold | Pixels (30m) | 20 | 600 |
8.3 Sensitivity Analysis Results
8.3.1 VDI Sensitivity Analysis
| Run | sBA weight | sCC weight | %>5 | %>6 | %>7 | Mean value |
|---|---|---|---|---|---|---|
| VDI_0 | 0.6667 | 0.3333 | 72.45 | 67.17 | 63.62 | 7.25 |
| VDI_1 | 0.3333 | 0.6667 | 78.52 | 73.52 | 63.70 | 7.54 |
| Removed | % > 5 | % > 6 | % > 7 | Mean Value |
|---|---|---|---|---|
| None | 37.46 | 33.35 | 30.86 | 7.40 |
| sBA | 39.98 | 38.28 | 36.22 | 7.83 |
| sCC | 37.13 | 32.67 | 31.72 | 7.94 |
8.3.2 SMII Sensitivity Analysis
| Run | sSoilK weight | sTRMI weight | %>5 | %>6 | %>7 | Mean value |
|---|---|---|---|---|---|---|
| SMII_0 | 0.3333 | 0.6667 | 54.32 | 20.84 | 6.90 | 5.01 |
| SMII_1 | 0.6667 | 0.3333 | 66.61 | 33.14 | 15.33 | 5.38 |
| Removed | % > 5 % > 6 % > 7 Mean | |||
|---|---|---|---|---|
| None | 27.00 | 12.01 | 2.58 | 5.20 |
| sSoilK | 17.55 | 8.89 | 2.67 | 5.06 |
| sTRMI | 36.22 | 14.55 | 14.55 | 6.09 |
8.3.3 SbII Sensitivity Analysis
| File | sLD | sPK | sPm | sPo | % > 5 | % > 6 | % > 7 | Mean |
|---|---|---|---|---|---|---|---|---|
| SbII_0 | 0.15 | 0.2 | 0.4 | 0.25 | 99.92 | 73.78 | 53.95 | 6.65 |
| SbII_1 | 0.15 | 0.2 | 0.25 | 0.4 | 99.95 | 58.80 | 0.72 | 5.97 |
| SbII_2 | 0.15 | 0.4 | 0.2 | 0.25 | 99.97 | 59.41 | 55.62 | 6.58 |
| SbII_3 | 0.15 | 0.4 | 0.25 | 0.2 | 99.97 | 59.41 | 54.42 | 6.81 |
| SbII_4 | 0.15 | 0.25 | 0.2 | 0.4 | 99.92 | 59.37 | 0.72 | 5.95 |
| SbII_5 | 0.15 | 0.25 | 0.4 | 0.2 | 99.92 | 74.82 | 54.05 | 6.86 |
| SbII_6 | 0.2 | 0.15 | 0.4 | 0.25 | 99.93 | 73.78 | 36.96 | 6.55 |
| SbII_7 | 0.2 | 0.15 | 0.25 | 0.4 | 99.95 | 40.03 | 0.72 | 5.86 |
| SbII_8 | 0.2 | 0.4 | 0.15 | 0.25 | 99.96 | 59.41 | 55.62 | 6.45 |
| SbII_9 | 0.2 | 0.4 | 0.25 | 0.15 | 99.94 | 58.97 | 53.41 | 6.91 |
| SbII_10 | 0.2 | 0.25 | 0.15 | 0.4 | 99.89 | 59.36 | 0.78 | 5.82 |
| SbII_11 | 0.2 | 0.25 | 0.4 | 0.15 | 99.93 | 94.63 | 53.64 | 6.97 |
| SbII_12 | 0.4 | 0.15 | 0.2 | 0.25 | 99.97 | 58.58 | 0.54 | 6.04 |
| SbII_13 | 0.4 | 0.15 | 0.25 | 0.2 | 99.96 | 58.60 | 7.97 | 6.27 |
| SbII_14 | 0.4 | 0.2 | 0.15 | 0.25 | 99.97 | 58.95 | 0.53 | 6.02 |
| SbII_15 | 0.4 | 0.2 | 0.25 | 0.15 | 99.97 | 58.63 | 36.33 | 6.48 |
| SbII_16 | 0.4 | 0.25 | 0.15 | 0.2 | 99.97 | 61.61 | 7.89 | 6.23 |
| SbII_17 | 0.4 | 0.25 | 0.2 | 0.15 | 99.97 | 61.63 | 36.38 | 6.46 |
| SbII_18 | 0.25 | 0.15 | 0.2 | 0.4 | 99.92 | 13.23 | 0.72 | 5.73 |
| SbII_19 | 0.25 | 0.15 | 0.4 | 0.2 | 99.93 | 73.80 | 36.97 | 6.65 |
| SbII_20 | 0.25 | 0.2 | 0.15 | 0.4 | 99.92 | 42.42 | 0.72 | 5.71 |
| SbII_21 | 0.25 | 0.2 | 0.4 | 0.15 | 99.92 | 94.56 | 53.63 | 6.86 |
| SbII_22 | 0.25 | 0.4 | 0.15 | 0.2 | 99.97 | 59.41 | 55.62 | 6.55 |
| SbII_23 | 0.25 | 0.4 | 0.2 | 0.15 | 99.97 | 59.53 | 54.35 | 6.78 |
| Removed | % > 5 | % > 6 | % > 7 | Mean |
|---|---|---|---|---|
| None | 51.84 | 29.63 | 4.14 | 6.38 |
| sLD | 51.82 | 29.20 | 28.00 | 6.60 |
| sPK5 | 50.91 | 14.58 | 0.37 | 5.88 |
| sPm | 32.28 | 22.01 | 2.53 | 5.75 |
| sPo | 51.83 | 38.63 | 27.67 | 7.28 |
8.3.4 Final Suitability Sensitivity Analysis
| File | Sbii | sSF | SMII | VDI | % > 5 | % > 6 | % > 7 | Mean |
|---|---|---|---|---|---|---|---|---|
| Final_0 | 0.2 | 0.4 | 0.2 | 0.2 | 87.41 | 63.28 | 25.75 | 6.23 |
| Final_1 | 0.2 | 0.2 | 0.4 | 0.2 | 82.07 | 55.72 | 16.32 | 6.00 |
| Final_2 | 0.2 | 0.2 | 0.2 | 0.4 | 78.67 | 66.14 | 49.18 | 6.45 |
| Final_3 | 0.4 | 0.2 | 0.2 | 0.2 | 90.33 | 67.17 | 27.60 | 6.33 |
| Removed | % > 5 | % > 6 | % > 7 | Mean |
|---|---|---|---|---|
| None | 43.75 | 33.03 | 15.90 | 6.25 |
| Sbii | 39.67 | 31.15 | 17.56 | 6.12 |
| sSF | 41.29 | 33.28 | 19.74 | 6.29 |
| SMII | 45.17 | 35.77 | 27.06 | 6.68 |
| VDI | 46.75 | 26.13 | 1.58 | 5.92 |
Citation
@online{e_lima,
author = {E Lima, Ryan and Gupta, Neha and Zalesky, Travis and Tsagaan
Sankey, Temuulen and E Springer, Abraham and D. Broxton, Patrick and
Jacobs, Katharine},
title = {Mapping Landscape Suitability for Thinning to Reduce
Evapotranspiration and Enhance Groundwater Recharge in Semi-Arid
Ponderosa Pine Forests},
langid = {en},
abstract = {Literature on the relationship between forest thinning and
water yield was used to develop suitability criteria to map where
forest treatment is most likely to enhance groundwater recharge
across the ponderosa pine (*Pinus ponderosa*) forests of Arizona.
Recharge in ponderosa pine forests is ephemeral and focused in
periods of snowmelt and locations of enhanced permeability when soil
moisture exceeds threshold levels. Our approach combines thematic
maps of criteria such as fraction of precipitation as snow, slope,
aspect, landscape morphology, forest basal area, canopy cover,
lineament density, lithology, and hydrologic soil type into a
GIS-Multi-Criteria Decision Analysis. Thematic layers were grouped
into indices and suitability was determined through a weighted sum
model of variables and indices. Approximately 10\% (182,000
hectares) of the 1.8 million *ha* of ponderosa pine forests across
western and northern Arizona were found to be highly suitable for
groundwater recharge enhancement from forest thinning. This study
finds that much of the area where thinning is already planned for
the purpose of restoring other ecosystem services may also improve
groundwater management in a state facing serious groundwater
declines. This GIS-based approach enables hydrologically informed
and spatially targeted forest management in semi-arid landscapes
from catchment to regional scales.}
}