Managed aquifer recharge as a drought mitigation strategy in heavily-stressed aquifers

Increasing meteorological drought frequency and rising water demand drive groundwater exploitation beyond sustainable limits. In heavily-stressed aquifers mitigation strategies, such as Managed Aquifer Recharge (MAR), are needed to restore depleted groundwater storage. MAR is also designed to overcome short dry periods. However, wider impacts of MAR as a drought mitigation strategy remain to be quantified. The objective of this study is to assess impacts of MAR in heavily-stressed aquifers using a case study of the Central Valley in California (USA). The novelty of this study lies in its analytical approach based on long-term observational data of precipitation, groundwater levels, and MAR operations. The impact of MAR operations is assessed regionally and for different temporal scales. Results show spatially-coherent clusters of groundwater level time series in the Central Valley representing three main patterns that manifest themselves in different groundwater drought characteristics and long-term trends. The first regional pattern shows lengthened groundwater droughts and declining groundwater levels over time, indicating effects of over abstraction in aquifer sections without MAR. The second regional pattern shows reduced groundwater drought duration and magnitude related to periodically rising groundwater levels, showing short-term MAR impacts. The third regional pattern shows alleviated groundwater droughts and groundwater levels show a long-term rise, representing long-term MAR impacts. Mitigated groundwater droughts and long-term rise in groundwater levels reveal the value of long-term MAR operations and their contribution toward sustainable groundwater management. Increased institutional support is recommended to ensure longevity of MAR and thereby amplify its success as regional drought mitigation strategy in heavily-stressed aquifers.


Introduction
Groundwater resources are increasingly pressured due to growing domestic, industrial, and agricultural water demand (Siebert et al 2010, Döll et al 2012. Overuse of groundwater resources can lead to heavily-stressed aquifers and severe groundwater depletion (Custodio 2002, Gleeson et al 2012, Famiglietti 2014) that can be aggravated following increased water use during meteorological drought events (Taylor et al 2013, Russo andLall 2017), which are more likely to occur in near future (Swain et al 2018). Groundwater depletion rates have already increased globally (Konikow 2011). Only few aquifers show decreasing groundwater depletion rates that are realised due to reduced withdrawals (mandatory or voluntarily), improved groundwater management, conjunctive use of surface water and groundwater, and enhanced groundwater recharge (Konikow 2011, Dillon et al 2019. These exceptional examples show the potential of mitigation strategies to ensure sustainable groundwater use and avoid maladaptation to meteorological droughts in heavily-stressed aquifers. Enhanced groundwater recharge can contribute to sustainable groundwater management and mitigate groundwater droughts, defined as a below-normal groundwater levels (Mishra and Singh 2010). However, empirical evidence to support this hypothesis is difficult to obtain, because few regions globally have a long-term practice of reported Managed Aquifer Recharge (MAR; Dillon et al 2019). One of these regions is the Central Valley in California (USA) that has a long history of MAR and is classified as a heavily-stressed aquifer (DWR 2016).
In the Central Valley, groundwater resources are heavily-stressed due to the extensive groundwater use to meet the high agricultural demand (Faunt 2009, Scanlon et al 2012b, which is likely to rise with an increased frequency of extreme meteorological droughts in near future (Swain et al 2018, Alam et al 2019. Current surface water imports from North California are insufficient to meet the agricultural water demand during meteorological droughts (Faunt 2009, Li et al 2018 exacerbating groundwater depletion with irreversible consequences , Ojha et al 2017. Without a change in its water management strategies, current groundwater depletion is likely to increase in the near future (Alam et al 2019), because existing drought mitigation strategies focus on maintaining agricultural production (Medellín-Azuara et al 2015) resulting in a direct increase in groundwater use to complement the deficit in surface water (Christian-Smith et al 2015) and an indirect increase due to a change in cropping pattern (Xiao et al 2017, Li et al 2018, Mall and Herman 2019. Initiatives to increase drought resilience, such as MAR, receive increasing attention with the implementation of Sustainable Groundwater Management Act (SGMA; SGMA 2014), although more evidence-based research is required to evaluate current impacts of MAR operations.
MAR, as a management strategy, is increasingly practised globally enhancing unmanaged surface water infiltration to benefit water quality and quantity (Dillon et al 2019). Reported global MAR capacity is greatest in India (31%) and USA (26%) aiming to meet the extensive agricultural groundwater demand (Dillon et al 2019). Most of India's MAR structures are, however, unreported and millions of structures are still in planning (Shah 2009, Stefan andAnsems 2018). In the USA, and in particular in Arizona and California, large-scale surface water basins are in use since the 1960s to increase drought resilience of agricultural water users (Scanlon et al 2012a, Megdal et al 2014. On a smaller scale, MAR contributes significantly to sustainable groundwater use in many European countries and Australia. For example, Germany's MAR, representing 9% of global MAR capacity, has increased the resilience of urban water supply systems since 1870 (Sprenger et al 2017). In Australia (4% of global MAR capacity), urban water supply was adapted to include MAR recharging groundwater in the aftermath of the Millennium Drought (Grant et al 2013, Radcliffe 2015. Moreover, national guidelines for applying MAR have been developed to aid safe development of MAR embedded in national water policies (Dillon et al 2020). In contrast to the growth of MAR globally, few areas are monitored that limits the assessment of MAR impacts.
The long-term MAR practice in the USA is therefore an exceptional example, illustrated by the well-monitored and studied MAR facilities in California. Previous MAR studies focused on individual (Scanlon et al 2012a) or multiple MAR facilities , temporal variability in groundwater levels (Thomas and Famiglietti 2015), and alternative strategies for enhanced recharge using agricultural fields for MAR (Ag-MAR; Dahlke et al 2018 (Alam et al 2020). A limiting factor is, however, the required additional infrastructure for conveyance, storage, and recharging of high magnitude flows (Yang and Scanlon 2019). Increased interest and highlighted potential of MAR contrast sharply with limited evidence-based research on the spatio-temporal impact of existing MAR operations. Assessing existing MAR impacts would show how MAR contributes to sustainable groundwater use and that would be of interest for further expansion of MAR within SGMA.
The aim of this study is to assess the impact of long-term MAR practices on groundwater droughts using a case study from the Central Valley Aquifer of California. To meet this aim, we focus on (1) spatial patterns in groundwater level time series, (2) shortterm and long-term patterns in precipitation and groundwater level time series, and (3) groundwater drought characteristics. The novelty of this study lies in its analytical approach that infers the long-term impact of MAR on groundwater droughts using longterm observational data of groundwater levels, precipitation, and MAR operations on a regional scale. This analytical approach differs from previous studies in the Central Valley, which used either groundwater models or water budgets to estimate impacts of actual or potential MAR sites , Xiao et al 2017, Ghasemizade et al 2019, Kourakos et al 2019, Maples et al 2019, Alam et al 2020. Previously, this method has been applied to study groundwater droughts in near-natural settings (Bloomfield and Marchant 2013, Kumar et al 2016, Haas and Birk 2017, characterising of groundwater dynamics (Haaf and Barthel 2018, Heudorfer et al 2019, and it was a starting point for modelling groundwater time series (Marchant and Bloomfield 2018) and investigating the asymmetric impact of groundwater use on groundwater droughts (Wendt et al 2020). This study shows how a similar analytical approach is useful to infer MAR impacts.

Study area
The study area in the Central Valley in California is the Tulare Basin, which extends across five counties (Kings, Kaweah, Tulare Lake, Tule, and Kern; figure 1). The Tulare Basin is currently overdrafted and represents critical groundwater conditions in Southern California (DWR 2016). The area has a hot Mediterranean climate, which is dry with winter precipitation (including rain and snow in the mountains) that is the main source of groundwater recharge. Average annual precipitation is 152-254 mm, which is exceeded by reference evaporation of 1295-1422 mm. In addition to incoming precipitation, the landscape and streams also recharge the alluvial aquifer representing 36% of the modelled long-term water budget (Faunt 2009). Outgoing components of the water budget are groundwater pumping (50%) and losses due to (in)elastic matrix storage, compression of specific yield, and groundwater flow leaving the Tulare Basin (14%; Faunt 2009). The long-term modelled groundwater balance is thus strongly negative due to large-scale groundwater abstractions (Faunt 2009, Alam et al 2020. These abstractions are currently not monitored or restricted, nor is the depth of abstractions reported, even though this might change with the implementation of SGMA (Thomas 2019). The impact of MAR is currently not included in the groundwater model.
The Tulare Basin is set between two mountain ranges and is topographically confined by the Tehachapi Mountains in the south and the Sierra Nevada Mountains in the east (Faunt 2009). From these mountain ranges, small streams fed by snow melt and run-off develop into rivers that recharge the layered alluvial fan aquifer in the valley. Discharge used to accumulate in the topographicallyclosed aquifer resulting in the second largest lake (Tulare Lake) in the U.S. until 1900, when it dried up due to intensification of groundwater use (Bertoldi et al 1991). Regional groundwater flow and recharge in the area are still affected by lake bed sediments that consist of fine-grained material and (Corcoran) clays forming a confining layer in the alluvial aquifer. Unconfined aquifer layers define the first 900 m in the Tulare Basin and groundwater levels, when unaffected by groundwater abstractions, drain internally to Tulare Lake county (Faunt 2009).
In addition to natural recharge, artificial recharge is provided by MAR facilities operating in conjunction with the State Water Project and Central Valley Project (USBR 2019, California Department of Water Resources 2020). These projects convey water from North California to the southern part of Central Valley since the 1960s (Kletzing 1987, Scanlon et al 2012a. To date, approximately 15 MAR facilities are in operation aiming to sustain agricultural water demand during dry years (KCWA 2010). These facilities use infiltration basins to recharge a mixture of available surface water, storm run-off, and imported surface water. All facilities comply with SGMA and water is treated prior to infiltration to maintain high standards for water quality (SGMA 2014, Kern County (CA) 2020).

Data and methods
The regional groundwater analysis is based on spatially-distributed precipitation, groundwater level, and operational MAR data in the Tulare Basin that cover a 35-year time period. The period of investigation started in 1980, when eight MAR facilities were (starting to be) operational, and ended in 2015 due to reduced availability of groundwater level observations after 2015. This 35-year period covers both extreme wet and dry periods including the exceptional drought in 2012-2015 (DWR 2016, Griffin and Anchukaitis 2014, Robeson 2015).

Data
Climate information was derived from two datasets with different spatial and temporal resolutions. The first dataset, the North American Drought Monitor (NADM), shows drought intensity and occurrence based on monthly standardised precipitation index (SPI) . NADM data are divided into climatological regions and San Joaquin basin that includes the Tulare Basin was used representing basin-wide monthly meteorological drought conditions (Lawrimore et al 2002). These meteorological drought conditions are defined using the NADM threshold for moderate droughts (SPI < − 0.8; Svoboda et al 2002). For consistency, we applied the same (opposite) threshold (SPI > 0.8) for wet conditions. The second used dataset has a higher temporal and spatial resolution and consisted of gridded (1 km 2 ) daily precipitation estimates from the DAYMET database by Thornton et al (2017), who spatially interpolated station measurements. This location-specific data were used to compare groundwater at a site with extracted grid cells of precipitation estimates. Extracted daily precipitation amounts were summed to 6-month totals to meet the biannual of groundwater time series and allow standardisation into SPI (Mckee et al 1993, Wu et al 2007. Groundwater data were obtained from the CAS-GEM dataset managed by the Department of Water Resources (CASGEM 2017) that contains about 12,100 groundwater monitoring sites (monitored biannually) in the Tulare Basin. Despite the high number of monitoring sites, few sites are monitored regularly, which is a key requirement to standardise  (1983-2011), West Kern (1988, Rosedale (1989-2011), Pioneer (1995-2011), Semitropic (2005, and Waldron (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007). Selected groundwater monitoring sites (149) are shown in red. groundwater levels in a regional comparison (Bloomfield and Marchant 2013). Therefore, time series were initially selected on the starting year (1980) and further screened for missing data, which were linearly interpolated in case of minor gaps (⩽3 observations; Thomas et al 2017, Tallaksen andvan Lanen 2004). This reduced the dataset from 12.128 sites to 3030 sites that starting monitoring in 1980, and finally to 149 sites that are regularly monitored. Regular groundwater monitoring was discontinued in 2015 for 39% of the selected wells, hence we set the time period to 1980-2015 to optimise the spatial representation of groundwater data. Final selection of sites consisted thus of 149 sites, of which 20 were located inside MAR facilities (figure 1). These sites were flagged as 'MAR observation wells' , as MAR impacts would first be observed in these sites .
The operational MAR data obtained from KCWA (2010) and Scanlon et al (2016) contained detailed information regarding recharged (infiltrated) and retrieved (extracted) water volumes at MAR facilities. Reported annual volumes were obtained for 11 out of 15 MAR facilities. Annual volumes were converted to millimetres per year using the reported size of recharge basins at the time of infiltration (KCWA 2010). Two MAR facilities out of eleven were in operation during the entire 35-year time period, other facilities partially covered the time period (see caption figure 1) resulting in a variable temporal and spatial extend of the dataset. The total reported volume of MAR is 106.4 mm yr −1 that is equivalent to 42-70% of the long-term precipitation (152-254 mm yr −1 ) in the Tulare Basin.

Methods
Regional patterns in groundwater level time series were investigated by clustered standardised time series. First, the standardised groundwater index (SGI) was calculated for all 149 groundwater level time series using the non-parametric approach of Bloomfield and Marchant (2013). Two clustering techniques (K-means and Ward's minimum) were applied to the SGI time series, applying both numerical and hierarchical reduction of similarity to reduce bias in the choice of clustering method (Haaf and Barthel 2018). Both methods resulted in similar regional patterns for 85% of the cluster composition, as shown in figure S1 (stacks.iop.org/ERL/16/014046/mmedia) modified from Wendt et al (2019). Only in Kern county, SGI clusters were slightly more spatially-coherent for Ward's minimum clustering (Wendt et al 2019). Hence, these SGI clusters were used for further analysis.
Temporal patterns were analysed using the 35year time series of groundwater levels and precipitation estimates for each site. The occurrence and strength of anomalies in time series were analysed for both short-term (in decades) and long-term (whole time series) patterns, following the same approach as Thomas and Famiglietti (2015). Short-term analysis consisted of SPI and SGI time series subdivided in three decades starting from 1980 (1980-1989, 1990-1999, and 2000-2009) and a remainder of 5 years using the average of 10 years to evaluate climatic controls to groundwater variation (Thomas and Famiglietti 2015). Long-term temporal patterns were analysed using a monotonic Mann-Kendall trend test (Mann 1945, Kendall 1949) that was modified for serial correlated groundwater data (Yue et al 2002, Hamed 2008. Biannual groundwater level observations were averaged for each calendar year to meet the requirements of the modified trend test. Trends in precipitation were tested using a standard Mann-Kendall trend test for annual precipitation totals. Trends were considered significant when the trend Z value was either < − 2.56 or > 2.56 (α = 0.01) and its strength was measured by the Sen slope of the (modified) Mann-Kendal trend test.
Both meteorological and groundwater droughts were quantified in the drought analysis by applying the NADM threshold to SPI and SGI time series (SPI or SGI < − 0.8). This threshold approach was used to calculate drought duration, length of time in which SPI or SGI time series were below the drought threshold, and magnitude, which is the accumulated SPI or SGI during a drought (Tallaksen andvan Lanen 2004, Mishra andSingh 2010).

Results
The results consist of three subsections. First, spatial patterns in groundwater level time series are presented showing all SGI clusters and main regional patterns identified in the Tulare Basin. Second, shortterm (decades) and long-term (35-year) patterns in precipitation and groundwater time series are shown, and lastly, groundwater drought characteristics are presented for the 35-year period and the 2012-2015 meteorological drought event.

Spatial patterns in groundwater level time series
Three regional spatial patterns were found in the clustered SGI time series (figures 2(a) and (b)). Six out of seven SGI clusters represent a spatially coherent group of groundwater monitoring sites (figure 2(a)) that can be summarised in three main regional patterns showing a declining, variable, and rising SGI in the 35-year period ( figure 2(b)). The first regional pattern (RP1) shows a decline in SGI over time ( figure 2(b)). RP1 is represented by cluster 1 (CL1) that consists of 31% of all monitoring sites in the northern counties of the Tulare Basin (Kings, Kaweah, and Tulare Lake counties). The cluster mean of CL1 shows a strong decline with below-average SGI (SGI <0) since 2002. The second regional pattern (RP2) shows periodic variation in SGI and is distinguished from the first by peaks in SGI in 2000-2001, 2005-2006, and 2011-2012 ( figure 2(b)). During these periods, SGI rises to average or above-average conditions (SGI > 0). RP2 consists of clusters 2, 3, 4 (CL2, CL3, CL4) and represents 55% of all monitoring sites located mainly in Tule and Kern counties. RP2 also includes 16 of the 20 'MAR observation wells' grouped in CL3, which cluster mean is sharply rising and falling over time.
The third regional pattern (RP3) shows a consistent rise in SGI for most of the 35 years ( figure 2(b)). Clusters 5 and 6 (CL5 and CL6) represent RP3, located in Tule and Kern counties (14% of monitoring sites). The consistent rise in SGI contrasts with the declining SGI in RP1 for the same period. Not surprisingly, CL5 and CL6 are located in different aquifer sections. CL5 is a smaller cluster of six wells located along a line between Tule and northern Kern counties. CL6 is located in southern Kern county, which is the most southern and the lowest section of the Central Valley Aquifer.

Temporal patterns in precipitation and groundwater level time series
Contrasting temporal patterns are also found in short-term (decadal) SPI and SGI averages in the Tulare Basin ( figure 2(c)). Short-term SPI is aboveaverage in the first decade (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989) and belowaverage in the second and third decade (1990-1999 and 2000-2009) for all clusters. The last 5 years (2010-2015) include an extremely wet period and an extreme meteorological drought resulting in, on average, normal conditions during the 5 years. SGI differs from this decadal SPI pattern (see figure 2(c). In CL1, moderately negative to strongly negative SGI are observed in 2000-2015 exceeding the SPI in this period. The SGI in CL2, CL3, and CL4 follow the SPI pattern relatively well, although the SGI declines in the last 5 years resulting in strongly negative SGI for both clusters. In CL5 and CL6, a rise in SGI is found that exceeds the SPI in the last 15 years (CL5) and in the second and third decade (CL6), but declines to below-average conditions in 2010-2015. This decline is thus reflected in all but one cluster and shows average groundwater conditions during extreme meteorological drought in 2012-15, impacting most groundwater monitoring sites in the Tulare Basin.
Long-term (35-year) temporal patterns in groundwater level time series show both negative and positive trends that are distributed unevenly in the Tulare Basin. Significant negative trends are primarily found in northern counties (Kings, Kaweah, and Tulare Lake) compared to moderately negative, neutral, and positive trends in the southern counties (Tule and Kern; figure 3(a)). Significant positive trends are detected for 12 sites suggesting groundwater levels have risen consistently. This is remarkable considering negative precipitation trends ( figure 3(b)) and the . Six out of seven clusters represent a regional pattern and are further analysed. CL7 represents a local pattern and is therefore excluded. The centre figure (b) shows clustered SGI time series of these six regional clusters. For each cluster, the cluster mean (thick line) and the range (minimum and maximum) in SGI are shown in the six panels. Colours are matching the map 2(a). On the right (c), the average cluster SPI (grey) and SGI (coloured) is plotted for each decade (horizontal bars and individual points for groundwater sites). Decades are divided into 1980-1989, 1990-1999, 2000-2009, and the 5 remaining years (2010-2015).
considerable groundwater abstraction in the Tulare Basin.
The strength of identified trends reflect the long-term regional spatial patterns (RP1-3; figure 4). Significant negative trends have a linear decrease in groundwater level exceeding of, on average, 0.79 m yr −1 with extreme outliers in RP1 (CL1). Sites in RP2 (CL2-4) decrease less per year compared to RP1, but average trend Z values remain significant for CL2 and CL4. CL3 marks a transition from negative to positive trends. This cluster includes most 'MAR observation wells' and represents sites located in the vicinity of four large MAR facilities (Kern Water Bank, City of Bakersfield, Kern River Channel, and Pioneer) suggesting that long-term MAR impact contributes to these moderately and significant positive trends. RP3 (CL5 and CL6) consists of mainly positive trends, increasing on average 0.2 m yr −1 , which are found close to the Arvin Edinson MAR facility (CL6).

Groundwater drought characteristics
Groundwater drought characteristics are also summarised by RP1-3 and show three different patterns in the Tulare Basin (RP1-3; figure 5). In RP1 (CL1), groundwater drought ( In RP2 (CL2-4), groundwater droughts occurred more often and recovered during above-normal precipitation conditions. SGI values in RP2 recovered more quickly compared to CL1 resulting in a (brief) periods of above-normal SGI between drought events. Drought recovery is highest for sites in CL3 that includes most 'MAR monitoring wells' . The sharp rise in SGI could be due to a combination of abovenormal precipitation and additional recharge supplied by MAR facilities, which recharge most water during above-normal precipitation conditions (see bottom panel in figure 5). The synchronised rising and falling SGI in CL2-4 suggests that most groundwater monitoring sites reflect this combined effect in Tule and Kern counties. The last regional pattern (RP3; based on CL5 and CL6) shows a rise in SGI resulting in entirely different groundwater drought characteristics. Groundwater droughts are observed at the start of the investigation period in 1991-1996(CL5) and 1980-1985. During meteorological droughts in 2007-2009 and 2012-2015, SGI declined gradually for both clusters, but the NADM threshold was not exceeded resulting in alleviated groundwater droughts compared to the other regional patterns.
The increment of MAR volumes since 1993 coincides with short-term increases in SGI in RP2 and long-term increasing SGI in RP3 resulting in different groundwater droughts. This is remarkable, as precipitation trends are decreasing in Tule and  Kern counties ( figure 3(b)). Most MAR volumes were recharged during periods of above-normal precipitation resulting in large MAR contributions in addition to the natural recharge. However, actual recharged volumes remain uncertain as documented volumes did not cover the complete 35-year period and might be higher than shown in figure 5.
In general, the groundwater drought duration in RP2 and RP3 halved compared to RP1 and drought magnitude reduced (figure 6). Average drought durations in RP2 and RP3 were both around a year, compared to two years for RP1. Average drought magnitude was close to −1.2 compared to −2.1 for RP1. The maximum drought duration was slightly shorter for RP3 and a lot shorter for most events in RP2. Maximum drought magnitude also decreased, although high outliers are still observed in RP2 and RP3.
Most severe groundwater droughts events (measured in drought magnitude) occurred at different times in the three regional patterns. In RP1, most severe droughts occurred before and during the extreme meteorological drought in 2012-2015. This contrasts with RP2 and RP3, where most severe drought events occurred in the period 1980-00. This is remarkable, as driving meteorological droughts in 1980-00 were less severe than the 2012-2015 drought (Robeson 2015).
The extreme meteorological drought in 2012-2015 had also a mixed impact in the Tulare Basin (figure 7). In RP1, groundwater droughts were severe and lasted 3.3 years on average. Drought duration in RP1 covered the entire meteorological drought and possibly longer, as the drought continues beyond the analysis period. Groundwater droughts in RP2 and RP3 started later and lasted 2.4 years on average, which is significantly (p = 4.6 × 10 −6 ) shorter. The later start resulted in lower drought magnitude (−0.8 less on average), which is a significantly (p = 7.3 × 10 −3 ) different compared to RP1. Largest reductions and absent droughts are observed in the vicinity of MAR facilities or within MAR facility boundaries, where above-average antecedent conditions prevented SGI values from crossing below the drought threshold. Figure 5. Groundwater droughts are shown for the identified three regional patterns in the cluster analysis. In the first panel, standardised precipitation (SPI12) is shown based on Drought Monitor data for San Joaquin basin (includes the Tulare Basin) . Meteorological droughts (below-average precipitation SPI12 < −0.8) are shaded and shown as light red surfaces in the other panels. Similarly, above-average precipitation (SPI12 > 0.8) are marked by light blue surfaces. The three regional patterns are shown in panels 2-4 with cluster means in matching colours according to figure 2. Groundwater drought events are shaded. The fifth panel shows recharged (blue) and extracted (red) MAR volumes of 11 MAR facilities. The stacked bar plot visualises recharged or extracted volumes (in (m) in each year, reported by individual facilities. Most reports (9 out of 11) did not contain most up-to-date MAR volumes and were updated until 2010 (see caption of figure 1). It is therefore plausible that more water was recharged during 2010-15 than shown here.  35-year period (1980-2015), as observed in the three regional patterns. Groundwater drought duration is measured in years. Groundwater drought magnitude is measured in accumulated SGI over the drought period.

Discussion
In the Tulare Basin, groundwater drought occurrence, duration, and magnitude change from north to south according to three regional patterns in longterm groundwater level variations. The first regional pattern (North Tulare Basin: Kings, Kaweah, and Tulare Lake counties) shows a long-term decline in groundwater levels, which resulted in extended groundwater droughts, as deficits in groundwater storage were not replenished despite above-average precipitation. The second regional pattern (South Tulare Basin: Tule and Kern counties) shows rising groundwater levels during periods of above-normal precipitation resulting in shorter droughts and rapid drought recovery. Long-term trends are moderately negative, neutral or even positive. The third regional pattern is found in a smaller southern section of South Tulare Basin. Here, groundwater levels rose consistently from 1995 onwards. Significant positive trends suggest an increase in groundwater storage over the past 35 years that alleviated droughts in 2012-2015.

Regional patterns in groundwater level variations
The long-term decline in groundwater levels in the first regional pattern has been related to the continuous overuse of groundwater (Faunt 2009, Famiglietti et al 2011, Thomas et al 2017, Scanlon et al 2012b. Faunt (2009) found an additional non-linear increase of groundwater abstractions during dry years that explains the discrepancy between precipitation and groundwater anomalies in the short-term (decadal) SPI and SGI comparison ( figure 2(c)). The natural drought propagation was altered, as groundwater conditions remained below-normal despite periods of above-normal precipitation resulting in an extended groundwater drought in 2007-2015 presumably driven by long-term overuse of groundwater.
Regional short-term MAR impacts in the Tulare Basin are seen in the second regional pattern, representing a gradual change from declining to rising groundwater levels. Since the 1960s, groundwater has been recharged in MAR facilities aiming to overcome short dry periods showing temporary rising groundwater levels for single sites due to MAR recharge (KCWA 2010. This is confirmed by clusters in Tule and Kern counties showing rising SGI following above-normal precipitation and MAR recharge. The amplified periodic rise in groundwater levels was also noted by Xiao et al (2017). However, our results suggest that regional MAR impact is larger than previously assumed. Short-term MAR impacts are observed for the majority of sites in Tule and Kern counties and were not, or less strongly, observed in other counties. Moreover, groundwater deficits were quicker replenished in Tule and Kern counties. This rapid drought recovery is largest in the vicinity of MAR facilities and synchronises with recharged MAR volumes. As a result of regional short-term MAR impact, groundwater droughts reduced significantly in duration and magnitude.
The long-term rise in groundwater levels, found in the third regional pattern, shows that groundwater storage is (slowly) increasing despite a negative precipitation trend. Observed rising groundwater levels are probably due to a combination of long-term MAR practice and regional hydrogeological conditions. In this region, groundwater storage (natural and artificial) accumulates resulting from the dominant North-South regional groundwater flow and topographic confinement of the Central Valley Aquifer (Faunt 2009, p.49). Previous studies indicate a steady increase in groundwater as a consequence of MAR practices in the Central Valley (Faunt 2009, Scanlon et al 2012a, Thomas 2019. A similar, local increase in groundwater storage due to MAR impact was found in Coachella Valley in California (Thomas and Famiglietti 2015). However, findings of this study show a larger extent of the long-term rising groundwater levels in particular in aquifer sections where groundwater storage naturally accumulates. This accumulation of groundwater storage results in alleviated groundwater droughts in Kern and Tule counties during the last extreme drought in 2012-2015. This illustrates the potential of MAR as a measure to enhance drought resilience  that is also effective on a regional scale. The regional increase in drought resilience was also found by Thomas (2019), who analysed observed and remotely-sensed groundwater anomalies.

Implications for management
The short-term and long-term MAR impacts highlight the contribution of MAR in a heavily-stressed aquifer to sustainable groundwater management. However, increasing groundwater recharge using MAR is only a partial solution (Dillon et al 2012). Sustainable use of groundwater implies that groundwater use is in balance with (natural and artificial) recharge. Monitoring groundwater levels, MAR volumes, and groundwater abstractions would enable water managers to evaluate sustainability of abstractions, as measurable objectives are essential for sustainable groundwater management (Gleeson et al 2012, Thomas 2019) that would also inform water managers whether groundwater is abstracted from deeper confined layers (non-renewable) or shallower unconfined (renewable) sections of the aquifer. This information is crucial to assess contribution of MAR to sustainable groundwater management, as MAR in (semi-)confined aquifers requires a different technique compared to unconfined aquifers (Bouwer 2002). Unconfined aquifers can be recharged with enhanced surface water infiltration, such as Flood-MAR and Ag-MAR that can result in a short-term increase (Kocis and Dahlke 2017, Dahlke et al 2018, Ghasemizade et al 2019 and a long-term rise in unconfined groundwater storage (Niswonger et al 2017, Gailey et al 2019. MAR in semi-confined aquifer sections would only be impacted if lateral spread of additional recharge allows seepage (i.e. preferential path ways) to deeper sections of the aquifer (Faunt 2009). The value of MAR to sustainable groundwater management depends thus partly on the type of groundwater use and MAR contribution for which regular monitoring can be used to ensure sustainability objectives are met.
Encouraging MAR in heavily-stressed aquifers requires not only available water to infiltrate and potential to store water safely in the aquifer, but also careful implementation of MAR practices (Dillon et al 2019). The growth in MAR in the past 60 years suggests that MAR is going to play an important role in groundwater management (Dillon et al 2019). Although this should also be accompanied by a policy framework to ensure its correct implementation and safe development of MAR. Australia (only representing 4% of global MAR capacity) is the first by having a risked-based MAR guidelines in place since 2009 (Dillon et al 2020). Less strict guidelines are found in India (30% of global MAR capacity; Dillon et al 2019). In Europe, EU member states are encouraged to develop their own policies resulting in varying practices and applications (Capone andBonfanti 2015, Sprenger et al 2017). In the US, legislation is in place to secure the water quality of infiltrated water and guidelines on safe implementation, but MAR should be further included in water policies before encouraging further expansions (i.e. Flood-MAR, Ag-MAR) (Kiparsky et al 2017). Potential MAR capacity to store water may exceed surface reservoir capacity in many US states , Maples et al 2019 and additional funding could, for example, facilitate additional infrastructure to capture high magnitude flows (Kocis and Dahlke 2017) increasing the water availability for MAR (Dahlke et al 2018, Alam et al 2020. Despite some state funding being available, water demand still exceeds the planned capacity (Rohde et al 2014). Using high magnitude flows would, however, also require additional water treatment to avoid deterioration of groundwater quality (Dillon 2005, Yang and, which highlights the importance of careful implementation of MAR practices.

Conclusions
The impact of long-term MAR on groundwater droughts has been identified using an analytical regional groundwater analysis applied to groundwater observations in the Central Valley of California. Presented results show that regional MAR impact is larger than previously estimated confirming the importance of MAR in heavily-stressed aquifers. Groundwater droughts are reduced and even alleviated in aquifer sections as a result of short-term and long-term impacts of MAR. Short-term MAR impacts result in rapid groundwater drought recovery and thereby reduced groundwater drought duration and magnitude. Long-term MAR impacts are reflected in neutral and positive groundwater trends in the Central Valley. Despite a negative trend in longterm precipitation, groundwater trends are neutral and even significantly positive in the vicinity of MAR facilities located in the most southern, topographically confined, section of the Central Valley Aquifer. The consistent increase in groundwater levels resulted in alleviated groundwater droughts during the extreme meteorological drought in 2012-2015 showing the potential of MAR as regional drought mitigation strategy.
Neutral and positive trends in groundwater level data show that groundwater levels were maintained thanks to long-term (35 years) MAR practices despite a long-term reduction in precipitation and continuous use of groundwater. The transition from negative trends in the North to neutral and positive trends in the South stresses the significant contribution of long-term MAR practices to sustainable groundwater use in heavily-stressed aquifers. However, this success highly depends on water availability, capacity and infrastructure to store water, and careful MAR implementation. Institutional support and guidance for a safe implementation are required to ensure longevity of MAR facilities and their success as a drought mitigation strategy.
Further research on MAR impacts and sustainable groundwater management in California could address the largely unknown groundwater use, which could not be included in this study due to unknown groundwater abstractions. Enhanced monitoring of the Californian groundwater resource might be encouraged with the implementation of SGMA and wider inclusion of MAR within the new water policies is recommended. Advancing techniques to fill in missing groundwater level observations in a human-modified are also recommended, as this limited the spatial coverage of the current analysis.
Future applications of the presented analytical method could aid to assess the impact of individual MAR facilities or to identify MAR facilities of unknown performance. For example, groundwater monitoring sites located within MAR facility boundaries were primarily found in one cluster (CL3) showing that short-term MAR impact can be identified regionally without documented MAR recharge, which opens the door for advanced analytical methods focusing on quantification of groundwater level dynamics (Heudorfer et al 2019) that could isolate short-term MAR impacts. In conclusion, we presented a versatile method that enables researchers to evaluate short-term and long-term MAR impacts on groundwater droughts and thereby, we have shown that MAR can be used as a regional drought mitigation strategy in heavily-stressed aquifers.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https://water.ca.gov/Programs/Groundwater-Management/Groundwater-Elevation-Monitoring-CASGEM.

Acknowledgments
In the process of developing this research, we would like to thank a number of colleagues and collaborators for their support. DW would like to thank the Bureau of Economic Geology at University of Texas and Amir Aghagouchak and his research group for hosting a research visit to CA organised within the IAHS Panta Rhei network, working group 'Drought in the Anthropocene' . DW also thanks Mike Jones (Thames Water) for his support and thoughts on the potential of MAR, Claudia Faunt (US Geological Survey) for her thoughts on initial findings and sharing useful resources and model outputs, and lastly, the Female Geography writing group for supporting the writing process. In addition, DW would like to acknowledge her funding by CENTA NERC (NE/lL002493/1) and CASE studentship of British Geological Survey (GA/16S/023). AVL would like to acknowledge funding from the NWO Rubicon project 'Adding the Human Dimension to Drought' (2004/08338/ALW). This work also contributes to the objectives of the NERC-funded 'Groundwater Drought Initiative' (NE/R004994/1).