Lessons to Be Learned: Groundwater Depletion in Chile’s Ligua and Petorca Watersheds through an Interdisciplinary Approach

: Groundwater (GW) is the primary source of unfrozen freshwater on the planet and in many semi-arid areas, it is the only source of water available during low-water periods. In north-central Chile, there has been GW depletion as a result of semi-arid conditions and high water demand, which has unleashed major social conﬂicts, some due to drought and others due to agribusiness practices against the backdrop of a private water management model. The Ligua and Petorca watersheds in the Valpara í so Region were studied in order to analyze the inﬂuence of climatic and anthropogenic factors on aquifer depletion using an interdisciplinary approach that integrates hydroclimatic variables, remote sensing data techniques, and GW rights data to promote sustainable GW management. The Standardized Precipitation Index (SPI) and Normalized Di ﬀ erence Vegetation Index (NDVI) were calculated and the 2002–2017 land-use change was analyzed. It was shown that GW decreased signiﬁcantly (in 75% of the wells) and that the hydrological drought was moderate and prolonged (longest drought in the last 36 years). The avocado-growing area in Ligua increased signiﬁcantly—by 2623 ha—with respect to other agricultural areas (higher GW decrease), while in Petorca, it decreased by 128 ha. In addition, GW-rainfall correlations were low and GW rights were granted continuously despite the drought. The results conﬁrmed that aquifer depletion was mostly inﬂuenced by human factors due to overexploitation by agriculture and a lack of water management.


Introduction
Groundwater (GW) is a very valuable resource as it accounts for 96% of the planet's unfrozen freshwater [1]. As a part of the hydrological cycle, GW is an essential source of water for human use and agricultural and industrial activities [2], as well as the main source of water in arid and semi-arid regions [3,4]. It is estimated to supply 36% of drinking water, 42% of irrigation water, and 24% of industrial water worldwide [1].

Figure 1.
Conceptual diagram of the relationship among groundwater (GW), agriculture, and rainfall in the study area. ET-evapotranspiration, R-rainfall, and SNW-snow.

Study Area
The Ligua and Petorca river watersheds (32° S, 1980 km 2 and 1986 km 2 , respectively) are located in the Valparaíso Region (north-central Chile) (Figure 2b,c). The two rivers meet at the Salinas de Pullally coastal wetland and their watersheds have similar physical-geographical characteristics [17,20,26]. They are characterized by semi-arid conditions; a predominantly pluvial (pluvio-nival), mixed regime; and an absence of glaciers in their Andean sections, which makes them vulnerable in terms of water availability [16,17,23]. Annual rainfall varies between 253 and 319 mm, increasing in the upper parts and concentrated in the austral winter (April-October) [17,22,26]. The primary occurrence of GW is around the main channels in pre-cordillera valleys, where there is permeable fill made up of granular sediment, ranging from boulders to fine sand. In each of the valleys, there is an unconfined aquifer, both of which are highly permeable, resulting in a high degree of interaction between surface water and GW [17,26]. According to the study of Ayala et al. [19] and reaffirmed by an official report of the Chilean Water Directorate (DGA, for its initials in Spanish) [27], there are 12 hydrogeological sectors of common GW use (Figure 2c).

Study Area
The Ligua and Petorca river watersheds (32 • S, 1980 km 2 and 1986 km 2 , respectively) are located in the Valparaíso Region (north-central Chile) (Figure 2b,c). The two rivers meet at the Salinas de Pullally coastal wetland and their watersheds have similar physical-geographical characteristics [17,20,26]. They are characterized by semi-arid conditions; a predominantly pluvial (pluvio-nival), mixed regime; and an absence of glaciers in their Andean sections, which makes them vulnerable in terms of water availability [16,17,23]. Annual rainfall varies between 253 and 319 mm, increasing in the upper parts and concentrated in the austral winter (April-October) [17,22,26]. The primary occurrence of GW is around the main channels in pre-cordillera valleys, where there is permeable fill made up of granular sediment, ranging from boulders to fine sand. In each of the valleys, there is an unconfined aquifer, both of which are highly permeable, resulting in a high degree of interaction between surface water and GW [17,26]. According to the study of Ayala et al. [19] and reaffirmed by an official report of the Chilean Water Directorate (DGA, for its initials in Spanish) [27], there are 12 hydrogeological sectors of common GW use (Figure 2c).

GW and Rainfall Data
The 20 observation wells were obtained from the official site of the DGA, https://dga.mop.gob.cl/Paginas/default.aspx; those with the best data quality, temporal range, and hydrogeological characteristics were selected, constituting a period of 16 years (2002-2017) ( Figure  2c). The rainfall data of the 1980-2017 period were acquired from the Center for Climate and Resilience Research (CR2) and are available at http://www.cr2.cl/. Twenty weather stations were used and the gridded monthly rainfall dataset CR2MET was also obtained from the same web site. CR2MET is a spatially distributed dataset developed for Chile that is currently being used by the DGA to update the national water balance; it has a resolution of 0.05° and is available from 1976 to 2016 [28]. This gridded dataset is based on networks of land stations with the most complete data, a statistical downscaling of ERA-Interim reanalysis data and local topography [28,29].

Satellite Data
The MODIS/006/MOD13Q1 product (version 6) of the Land Processes Distributed Active Archive Center (LP DAAC, http://lpdaac.usgs.gov) was used to obtain MODIS (Moderate Resolution Imaging Spectroradiometer) NDVI time series data. The Terra images of this product have a spatial and temporal resolution of 250 m and 16 days, respectively [30,31]. To generate the cover and landuse results, Landsat images with a spatial resolution of 30 m from 2002 and 2017 were used (Table 1), which are available on the official site of the United States Geological Survey (USGS) https://landsatlook.usgs.gov/.

GW and Rainfall Data
The 20 observation wells were obtained from the official site of the DGA, https://dga.mop. gob.cl/Paginas/default.aspx; those with the best data quality, temporal range, and hydrogeological characteristics were selected, constituting a period of 16 years (2002-2017) (Figure 2c). The rainfall data of the 1980-2017 period were acquired from the Center for Climate and Resilience Research (CR2) and are available at http://www.cr2.cl/. Twenty weather stations were used and the gridded monthly rainfall dataset CR2MET was also obtained from the same web site. CR2MET is a spatially distributed dataset developed for Chile that is currently being used by the DGA to update the national water balance; it has a resolution of 0.05 • and is available from 1976 to 2016 [28]. This gridded dataset is based on networks of land stations with the most complete data, a statistical downscaling of ERA-Interim reanalysis data and local topography [28,29].

Satellite Data
The MODIS/006/MOD13Q1 product (version 6) of the Land Processes Distributed Active Archive Center (LP DAAC, http://lpdaac.usgs.gov) was used to obtain MODIS (Moderate Resolution Imaging Spectroradiometer) NDVI time series data. The Terra images of this product have a spatial and temporal resolution of 250 m and 16 days, respectively [30,31]. To generate the cover and land-use results, Landsat images with a spatial resolution of 30 m from 2002 and 2017 were used (Table 1), which are available on the official site of the United States Geological Survey (USGS) https://landsatlook.usgs.gov/. Water rights are granted to private users in perpetuity and at no cost. Current Chilean legislation (1981 Water Code) establishes that water is a national asset for public use and that individuals are granted the right to use it. This right of use is a real right and consists of the use and enjoyment of water [15,16]. There are two types of rights registered: surface and GW rights [32]. The GW rights data used were recorded by the DGA, which, according to Article 122 of the Water Code, collects information on original rights and requests associated with groundwater. Rights are granted to users, in accord with established regulations, for a given use and are granted by GW quantity and annual average streamflow [14,32]. In the DGA's Manual of Regulations and Procedures for Water Resources Administration [32], the nomenclature of each type of right is established, as shown below.
Use rights (ND, for the abbreviation in Spanish) are new rights that are granted in wells, pumping stations, culverts, or driven point wells. Intake point changes (VPC, for the abbreviation in Spanish) are rights that are authorized by the DGA to be transferred from an intake point within a hydrogeological sector. Right regularization (NR, for the abbreviation in Spanish) is requested by a user and consists of registration of a right in his or her name, be it an unregistered right or one that belongs to someone else, while old users (UA, for the abbreviation in Spanish) are those who have had GW rights for many years, be it through a water grant or a transaction [32].

Methods
This section presents the data processing steps, the analysis methods employed, and their integration through an interdisciplinary approach. The spatial-temporal and statistical analysis are carried out between watersheds and by the hydrogeological sector. Additionally, satellite image processing, drought and vegetation index calculation, land-use modeling using new programming technology, trend test statistical analysis, and general statistical analysis are described below. The general diagram is shown in Figure 3.

GW, Rainfall, and GW Rights Processing
The GW values of each well from the 2002-2017 period were analyzed on a monthly and annual basis, as were the in situ rainfall data. The annual GW averages were obtained by watershed and wells and represented on a map. In addition, the annual GW values were standardized in accordance with Hu et al. [33]. The standardization, Z, is carried out by subtracting the long-term mean value of the distribution from the annual value and dividing the result by the standard deviation, as shown in the following equation: where X is the value that is being standardized, µ the mean of the distribution, and σ the standard deviation of the distribution. The gridded rainfall data were processed and correlated with the in situ data to evaluate the representativeness of this source. The results presented Pearson coefficients between 0.98 and 0.99 (pvalue < 0.0001) at all the weather stations. The statistical analyses were carried out in XLSTAT 2019 v3.2 (Addinsoft, Paris, France).
The granted GW rights were processed at a commune or municipality level within the studied watersheds, with only those in the area selected. The data were processed at monthly and annual levels, in terms of both quantity of rights and average streamflow, resulting in a time series, and then new rights were correlated with GW levels (Table 1).

Drought Index
The Standardized Precipitation Index (SPI) [34] is an index that is widely used to identify droughts and wet events based on temporal rainfall data [35][36][37]. It is among the most used indices worldwide for drought monitoring and the World Meteorological Organization (WMO) considers it the universal drought index [35]. It is calculated by taking the difference between precipitation and the mean of a given time period and then dividing it by the standard deviation [36], as shown in Equation (2). Subsequently, it is fit to a gamma distribution and is transformed into a normal distribution [36][37][38].

GW, Rainfall, and GW Rights Processing
The GW values of each well from the 2002-2017 period were analyzed on a monthly and annual basis, as were the in situ rainfall data. The annual GW averages were obtained by watershed and wells and represented on a map. In addition, the annual GW values were standardized in accordance with Hu et al. [33]. The standardization, Z, is carried out by subtracting the long-term mean value of the distribution from the annual value and dividing the result by the standard deviation, as shown in the following equation: where X is the value that is being standardized, µ the mean of the distribution, and σ the standard deviation of the distribution. The gridded rainfall data were processed and correlated with the in situ data to evaluate the representativeness of this source. The results presented Pearson coefficients between 0.98 and 0.99 (p-value < 0.0001) at all the weather stations. The statistical analyses were carried out in XLSTAT 2019 v3.2 (Addinsoft, Paris, France).
The granted GW rights were processed at a commune or municipality level within the studied watersheds, with only those in the area selected. The data were processed at monthly and annual levels, in terms of both quantity of rights and average streamflow, resulting in a time series, and then new rights were correlated with GW levels ( Table 1).

Drought Index
The Standardized Precipitation Index (SPI) [34] is an index that is widely used to identify droughts and wet events based on temporal rainfall data [35][36][37]. It is among the most used indices worldwide for drought monitoring and the World Meteorological Organization (WMO) considers it the universal drought index [35]. It is calculated by taking the difference between precipitation and the mean of a given time period and then dividing it by the standard deviation [36], as shown in Equation (2). Subsequently, it is fit to a gamma distribution and is transformed into a normal distribution [36][37][38].
where Xi is monthly rainfall, Xj the mean value of the time period, and σ the standard deviation. The SPI has been used to establish relationships between drought and GW levels in some research [2,6,9,35,36]. This affirmation is based on the fact that a rainfall deficit could lead to a lack of GW recharge [2,6]. Therefore, the SPI between 1982 and 2017 was calculated at 19 weather stations and 191 grid points of the CR2MET dataset at different time scales (1,3,6,9,12, and 24 months) using the SPEI package in R Student software. The advantage of these time scales is their use, in the case of three months or less, in drought monitoring tasks; in the case of six months or less, in monitoring effects on agriculture; and, in the case of 9 months or more, in monitoring hydrological effects [39]. SPI 12 and 24 from 2002 to 2016 in various seasons of the year in Chile-summer (December, January, February), autumn, (March, April, May), winter (June, July, August), and spring (September, October, November) [40]-were represented spatially using the ordinary Kriging interpolation method in ArcGis 10.7.1, as it presented lower semivariogram errors compared to other methods. This method has proven to be superior to others and has been used in other investigations to interpolate rainfall, the SPI, and even GW [33,41,42]. Then, the annual SPI values were calculated to analyze them with the annual standardized GW values.

Normalized Difference Vegetation Index NDVI
The so-called spectral vegetation greenness indices derived from satellite data such as the NDVI are tools for monitoring vegetation conditions [44]. This index consists of the spectral reflectances of the normalized ratio of the near infrared and red bands [44] (see Equation (3)); it has been widely used in remote sensing of vegetation for decades [42]. It shows variations in sensitivity to global climate changes and sensitivity to precipitation; this sensitivity is particularly high in arid and semi-arid areas [42]. It is important to consider that studies with different data such as rainfall and NDVI, among others, provide key information on regional vegetation changes resulting from climatic and human factors [44].
where NIR is the near infrared band and R the red band. The NDVI data were processed using the programming technology offered by Google Earth Engine (GEE). Google Engine is a cloud platform for rapid analysis using the information technology infrastructure of Google [42,45,46]. This NDVI was calculated based on atmospherically corrected bidirectional surface reflectances that have been masked for water, clouds, heavy aerosols, and cloud shadows [30]. A monthly time series from 2002 to 2017 throughout the study area was obtained after processing 368 images. Then, these values and the filtered NDVI were analyzed for the area below 800 m as this is the maximum elevation that avocado crops have reached (previously identified in Google Earth software and extracted through the Digital Elevation Model (DEM) in ArcGIS 3D analysis tools); therefore, it was used as an estimate of the total agricultural area.

Land Use/Land Cover Change LULCC
The composite preprocessing began with geometric correction (reprojection to UTM-WGS84 coordinate system, zone 19 S) [11]. Subsequently, a topographic correction was applied to eliminate the illumination effects of the relief by applying the normalization method described by Richter et al. [47]. Then, for the atmospheric correction (TOA), the methodologies of Chander et al. [48] and the USGS (https: //www.usgs.gov/land-resources/nli/landsat/using-usgs-landsat-level-1-data-product) were followed. The TOA correction method consists of a transformation of the digital values (DN) of the pixels into reflectance values. Next, to obtain a better atmospheric correction, the Dark Object Subtraction technique was applied [49]. Both corrections were applied in ArcGIS, ENVI, and TerrSet software.
The land-use classification proposed in this study was based on the classification (or nomenclature) of Zhao et al. [50], which was modified according to the interests of the investigation ( Table 2). To generate the classification, approximately 600 field monitoring points were used, which served as the basis for taking 150 points for each class (10 classifications), for a total of 1500 points. Based on these points and other new points generated using band combinations, 1500 points were also established for the 2002 classification. Subsequently, a non-parametric classifier called Random Forest (RF) [51][52][53] was selected to generate the classification using the GEE platform. RF is a classifier that uses multiple decision trees to classify pixels and then performs voting to give each pixel a final category. This classifier uses two hyperparameters: the number of trees and √ n, where n is the number of bands. These hyperparameters were configured with 500 decision trees, in accordance with Zhang et al. [54].
To increase the discrimination and separability of the proposed classes, a series of additional bands was calculated and incorporated into the model. The bands generated to improve the classification of each image from 2002 and 2017 were: NDVI, NDWI (Normalized Difference Water Index), SAVI (Soil Adjusted Vegetation Index), MSAVI (Adjusted SAVI), NBR (Normalized Burn Ratio), DEM (Digital Elevation Model), and Slope. Once the classifications were carried out, the precision parameters of each image were calculated; to this end, the Kappa index and overall accuracy were calculated for each of the classifications.

Statistical Tests and Correlation
To analyze the trends of the monthly and annual GW levels, rainfall and NDVI data, non-parametric Mann Kendall tests [55,56], and Sen's slope [57] were used as they do not require normal distribution and homoscedasticity of data and present the advantage of accommodating missing measurements and atypical values in a time series [2,58]. The level of significance used in all the tests was a p-value < 0.05. The Pettitt test was used for the significant change points in the time series or abrupt changes [59]. Pettitt is another non-parametric test and has been widely used to detect change points in environmental and hydrological data [60]. When the p-value is less than 0.05, there is a significant change point and the series is divided into two parts [7]. The correlation and cross-correlation analyses were carried out by the hydrogeological sector ( Figure 2c) and not at a pixel scale as vegetation or crops can be affected both in the vicinity of the extraction well and in a large area very distant from it, as the water is extracted from one site and used for irrigation in others. A Pearson correlation was applied to the average of GW, gridded rainfall, SPI (scales 1, 3, 6, 9, 12, and 24 months), and NDVI by the hydrogeological sector up to 800 m. The aim of this analysis was to investigate spatial and temporal relationships among GW and rainfall (effect of rainfall over GW recharge), drought, and vegetation (likely effect of irrigation). In addition, a cross correlation (time lag of 1 to 6 months) was made between GW-rainfall and GW-NDVI. Cumulative GW rights, both number of rights and annual average streamflow granted, were also correlated with GW levels.
The k-means non-hierarchical cluster used by Chu [9] was applied to the GW values of each well to establish similarities or differences between the behavior of each well and thus obtain groupings of wells with similar characteristics through statistical analyses and then to find links between aquifer depletion and anthropogenic factors.

Temporal and Spatial Variation of GW and Rainfall
In Figure 4, a water level decrease is observed in the 2002-2016 period at all the analyzed wells, which is steeper at some wells in the Ligua watershed located in the south of the study area ( Figure 4b). In general, a water level decrease and many wells with values below the average, mostly from 2009 and 2010 (75% respect to annual average), are observed (Figure 4a The k-means non-hierarchical cluster used by Chu [9] was applied to the GW values of each well to establish similarities or differences between the behavior of each well and thus obtain groupings of wells with similar characteristics through statistical analyses and then to find links between aquifer depletion and anthropogenic factors.

Temporal and Spatial Variation of GW and Rainfall
In Figure 4, a water level decrease is observed in the 2002-2016 period at all the analyzed wells, which is steeper at some wells in the Ligua watershed located in the south of the study area ( Figure  4b). In general, a water level decrease and many wells with values below the average, mostly from 2009 and 2010 (75% respect to annual average), are observed (Figure 4a   Annual rainfall presented a decrease that exceeded 350 mm starting in 2003 (67.7% of the annual average), but, as with GW, the decrease was more critical from 2009 and 2010, as rainfall was below 200 mm for several years (Figure 4a,b). However, no differences in rainfall behavior between watersheds were observed, unlike with GW, as the decrease is homogenous throughout the area. Annual rainfall presented a decrease that exceeded 350 mm starting in 2003 (67.7% of the annual average), but, as with GW, the decrease was more critical from 2009 and 2010, as rainfall was below 200 mm for several years (Figure 4a,b). However, no differences in rainfall behavior between watersheds were observed, unlike with GW, as the decrease is homogenous throughout the area. In general, the SPI shows a moderate but very prolonged hydrological drought (2010-2015) that coincides with the previously identified period of greatest GW decrease, as well as the period with the longest drought of the last 36 years. In general, the SPI shows a moderate but very prolonged hydrological drought (2010-2015) that coincides with the previously identified period of greatest GW decrease, as well as the period with the longest drought of the last 36 years. Water 2020, 12, x FOR PEER REVIEW 11 of 25

GW and Drought
As previously verified, the most critical drought period was between 2010 and 2015, mainly in SPI 9, 12, and 24, which is reaffirmed in Figure 8, where each SPI scale is compared with the standardized GW values. In this case, the rainfall deficits as seen through SPI values present a relationship with the decrease in standardized GW values, which is strongest between 2010 and 2015. While the SPI values do not show different behaviors between the watersheds, the standardized GW values are lower in the Ligua watershed. This demonstrates that other factors such as GW rights, GW overexploitation for irrigation, and land-use change could have a greater influence on GW depletion than decreased rainfall. This idea or assumption will be addressed in the following sections.

GW and Drought
As previously verified, the most critical drought period was between 2010 and 2015, mainly in SPI 9, 12, and 24, which is reaffirmed in Figure 8, where each SPI scale is compared with the standardized GW values. In this case, the rainfall deficits as seen through SPI values present a relationship with the decrease in standardized GW values, which is strongest between 2010 and 2015. While the SPI values do not show different behaviors between the watersheds, the standardized GW values are lower in the Ligua watershed. This demonstrates that other factors such as GW rights, GW overexploitation for irrigation, and land-use change could have a greater influence on GW depletion than decreased rainfall. This idea or assumption will be addressed in the following sections.

GW and Drought
As previously verified, the most critical drought period was between 2010 and 2015, mainly in SPI 9, 12, and 24, which is reaffirmed in Figure 8, where each SPI scale is compared with the standardized GW values. In this case, the rainfall deficits as seen through SPI values present a relationship with the decrease in standardized GW values, which is strongest between 2010 and 2015. While the SPI values do not show different behaviors between the watersheds, the standardized GW values are lower in the Ligua watershed. This demonstrates that other factors such as GW rights, GW overexploitation for irrigation, and land-use change could have a greater influence on GW depletion than decreased rainfall. This idea or assumption will be addressed in the following sections. Water 2020, 12, x FOR PEER REVIEW 12 of 25

NDVI Time Series Data Analyses
The NDVI time series ( Figure 9) present a very slight decrease between 2009 and 2015; thus, vegetation is also affected in this period of aquifer depletion, which is also found in previous analysis.

LULCC 2002-2017. Natural and Human Implications
The land-use classifications ( Figure 10) present good precision (>90%) ( Table 3), with a general Kappa index of 0.91. The most difficult-to-classify uses were plantations of fruit trees such as avocado, due mainly to the similarity of their spectral signatures. Meanwhile, the uses that allowed the best discrimination were bare soil and water as their signatures are not particularly difficult to separate and classify. The LULCC (Land Use/Land Cover Change) up to 800 m is shown as it is the most important zone for the study area since the agricultural area is concentrated there. Regarding the distribution of the land-use categories, the predominance of native covers in the upper parts of the watershed, mainly shrubs and native forest, associated with sclerophyll forest, can be observed, while the valley regions were found to be dominated mainly by avocado, fruit tree, agricultural, and urbanization uses.

NDVI Time Series Data Analyses
The NDVI time series ( Figure 9) present a very slight decrease between 2009 and 2015; thus, vegetation is also affected in this period of aquifer depletion, which is also found in previous analysis. The NDVI is greater starting at 800 m in comparison with the rest of the area, which is a result of farmland and the effect of irrigation. The values in the watersheds vary between 0.24 and 0.49, decreasing to 0.24 in 2014, while up to 800 m, they range between 0.30 and 0.67, with lower values also found in 2014.

LULCC 2002-2017. Natural and Human Implications
The land-use classifications ( Figure 10) present good precision (>90%) ( Table 3), with a general Kappa index of 0.91. The most difficult-to-classify uses were plantations of fruit trees such as avocado, due mainly to the similarity of their spectral signatures. Meanwhile, the uses that allowed the best discrimination were bare soil and water as their signatures are not particularly difficult to separate and classify. The LULCC (Land Use/Land Cover Change) up to 800 m is shown as it is the most important zone for the study area since the agricultural area is concentrated there. Regarding the distribution of the land-use categories, the predominance of native covers in the upper parts of the watershed, mainly shrubs and native forest, associated with sclerophyll forest, can be observed, while the valley regions were found to be dominated mainly by avocado, fruit tree, agricultural, and urbanization uses.

LULCC 2002-2017. Natural and Human Implications
The land-use classifications ( Figure 10) present good precision (>90%) ( Table 3), with a general Kappa index of 0.91. The most difficult-to-classify uses were plantations of fruit trees such as avocado, due mainly to the similarity of their spectral signatures. Meanwhile, the uses that allowed the best discrimination were bare soil and water as their signatures are not particularly difficult to separate and classify. The LULCC (Land Use/Land Cover Change) up to 800 m is shown as it is the most important zone for the study area since the agricultural area is concentrated there. Regarding the distribution of the land-use categories, the predominance of native covers in the upper parts of the watershed, mainly shrubs and native forest, associated with sclerophyll forest, can be observed, while the valley regions were found to be dominated mainly by avocado, fruit tree, agricultural, and urbanization uses. Water 2020, 12, x FOR PEER REVIEW 13 of 25  Scrubland was the largest land-use, with 83,974 ha in 2002 and 90,208 ha in 2017. Native forest also occupied a large area. The scrubland is made up of sclerophyll arborescent species and is found near prairies or areas of prairie crop rotation, while native forest such as Mediterranean sclerophyll forest and Mediterranean spiny forest is a vegetation type that is denser [61]. Both uses decreased in the analyzed period. Reservoirs and bodies of water are not identified at the scale of the map, but were considered in this study as they are an important cover to analyze in the context of the investigation. Reservoir area increased from 102 ha a 121 ha, meaning that they continued to be built during the analyzed years despite the drought. The total agricultural area increased, but this growth was small and the individual behavior in each watershed was very different. The avocado (Persea americana Mill.) use was identified as a separate category due to the importance of the crop for the study area in economic and political terms, as well as its water demand. It accounted for 5023 ha in 2002 and 7518 ha in 2017. Table 4 shows the net change in land uses. Fruit trees, avocado, bare soil, scrubland, exotic plantations, reservoirs, and urban zones presented increases in area, while agricultural land, native forest, and water bodies decreased. Bare soil was the use with the biggest gain, which means that its area increased, mainly at the expense of native forest, thereby decreasing the moisture retention of the soil and infiltration. Reservoirs had a net change of 19 ha as the gains were greater than the losses. Although fruit tree plantations decreased in some areas, their gains far exceeded their losses; however, avocado had a large positive net change (2495 ha) compared to the other agricultural areas.
There were differences in these land-use changes between watersheds. Reservoirs increased and agricultural area decreased in Petorca, while only farmland dedicated to traditional crops and prairie crop rotation decreased in Ligua. In the latter watershed, fruit tree plantations increased by 631.8 ha and avocado plantations by 2623.6 ha. Most significant was the increase in the bare soil area in both watersheds and the large avocado increase in Ligua compared to other farmland, unlike in Petorca, where the area of this crop decreased by 128 ha (Figure 11).  Scrubland was the largest land-use, with 83,974 ha in 2002 and 90,208 ha in 2017. Native forest also occupied a large area. The scrubland is made up of sclerophyll arborescent species and is found near prairies or areas of prairie crop rotation, while native forest such as Mediterranean sclerophyll forest and Mediterranean spiny forest is a vegetation type that is denser [61]. Both uses decreased in the analyzed period. Reservoirs and bodies of water are not identified at the scale of the map, but were considered in this study as they are an important cover to analyze in the context of the investigation. Reservoir area increased from 102 ha a 121 ha, meaning that they continued to be built during the analyzed years despite the drought. The total agricultural area increased, but this growth was small and the individual behavior in each watershed was very different. The avocado (Persea americana Mill.) use was identified as a separate category due to the importance of the crop for the study area in economic and political terms, as well as its water demand. It accounted for 5023 ha in 2002 and 7518 ha in 2017. Table 4 shows the net change in land uses. Fruit trees, avocado, bare soil, scrubland, exotic plantations, reservoirs, and urban zones presented increases in area, while agricultural land, native forest, and water bodies decreased. Bare soil was the use with the biggest gain, which means that its area increased, mainly at the expense of native forest, thereby decreasing the moisture retention of the soil and infiltration. Reservoirs had a net change of 19 ha as the gains were greater than the losses. Although fruit tree plantations decreased in some areas, their gains far exceeded their losses; however, avocado had a large positive net change (2495 ha) compared to the other agricultural areas.
There were differences in these land-use changes between watersheds. Reservoirs increased and agricultural area decreased in Petorca, while only farmland dedicated to traditional crops and prairie crop rotation decreased in Ligua. In the latter watershed, fruit tree plantations increased by 631.8 ha and avocado plantations by 2623.6 ha. Most significant was the increase in the bare soil area in both watersheds and the large avocado increase in Ligua compared to other farmland, unlike in Petorca, where the area of this crop decreased by 128 ha (Figure 11).   Figure 11. Gains, persistence, and losses of avocado land.

GW Rights
The GW rights granted by the DGA (Figure 12) have contributed to the decrease in GW levels. Rights continued to be granted during the years analyzed in this study and the number of rights granted was greatest in 2006 and 2008. Use rights (ND) account for the great increase of rights in this analyzed time period because they are the new rights granted.
Although in recent years-since 2015-the number of new rights granted has decreased, possibly due to measures taken by the authorities, and these adjustments were made somewhat late considering the effect that they could have had on GW depletion. According to the data, most of the rights are consumptive for agricultural irrigation. Although this analysis does not consider values by drainage basin, as there is not enough information to place each right geographically and the data are available at a municipality level, these results show-at the scale of the entire study area-that the number of rights and associated streamflow are related to the decrease in GW. The number of new rights increased in all municipalities and it was also verified that most are concentrated among families dedicated to avocado exports.

GW Rights
The GW rights granted by the DGA (Figure 12) have contributed to the decrease in GW levels. Rights continued to be granted during the years analyzed in this study and the number of rights granted was greatest in 2006 and 2008. Use rights (ND) account for the great increase of rights in this analyzed time period because they are the new rights granted.
Although in recent years-since 2015-the number of new rights granted has decreased, possibly due to measures taken by the authorities, and these adjustments were made somewhat late considering the effect that they could have had on GW depletion. According to the data, most of the rights are consumptive for agricultural irrigation. Although this analysis does not consider values by drainage basin, as there is not enough information to place each right geographically and the data are available at a municipality level, these results show-at the scale of the entire study area-that the number of rights and associated streamflow are related to the decrease in GW. The number of new rights increased in all municipalities and it was also verified that most are concentrated among families dedicated to avocado exports.

GW, Rainfall, and Vegetation Trend Tests
The trend analysis of the non-parametric tests is shown in Table 5; the results indicate a significant decrease in GW (in 75% of the wells) and most of the Ligua watershed. GW decreases by season were also significant. The Pettitt test showed the greatest number of significant points in 2009, but each abrupt change in the data series may be determined by local anthropogenic factors. Rainfall was also analyzed during 1980-2017, according to a period established by the World Meteorological Organization [62]. The results showed that it also decreased, but not to a statistically significant extent; therefore, it cannot be stated that there was a trend in the series. Nor were there any abrupt changes during these 36 years; rather, rainfall decreased slightly and the decrease in 2009 and 2010 was also too small to generate abrupt changes. In addition, with the tests applied to the NDVI, it cannot be stated that there is a trend in the series.

GW, Rainfall, and Vegetation Trend Tests
The trend analysis of the non-parametric tests is shown in Table 5; the results indicate a significant decrease in GW (in 75% of the wells) and most of the Ligua watershed. GW decreases by season were also significant. The Pettitt test showed the greatest number of significant points in 2009, but each abrupt change in the data series may be determined by local anthropogenic factors. Rainfall was also analyzed during 1980-2017, according to a period established by the World Meteorological Organization [62]. The results showed that it also decreased, but not to a statistically significant extent; therefore, it cannot be stated that there was a trend in the series. Nor were there any abrupt changes during these 36 years; rather, rainfall decreased slightly and the decrease in 2009 and 2010 was also too small to generate abrupt changes. In addition, with the tests applied to the NDVI, it cannot be stated that there is a trend in the series.  Figure 13 shows the Pearson correlation coefficient values between the studied variables and GW by hydrogeological sector. The analysis was not performed in sectors 2, 9, and 11 due to the absence of wells there. The highest coefficients found were with SPI 12 and 24, oscillating between 0.46 and 0.68 at the 12-month scale and 0.52 and 0.80 at the 24-month scale, which reflects the relationship between GW and hydrological drought. In sector 7 (Ligua central valley), the correlations with all variables were low, which could be due to local human factors such as excessive GW extraction. The values with rainfall were very low (<0.20), meaning a low correlation with GW. There was a low, negative correlation in sector 11, but it is located near the coast and its rainfall behavior is determined by oceanic effects such as sea breezes, which promote the formation of fog and coastal clouds, leading to a cooler, wetter coastal climate [61]. The NDVI-GW correlation values were higher than the GW-rainfall correlation values, even reaching 0.49 in sector 3 (upper eastern of Petorca), and significant in most cases ( Figure 13). These GW-NDVI values are high considering that the study area is a semi-arid zone, where a low correlation between vegetation and GW is expected; this result reflects the presence of farmland. In general, the coefficients with the SPI and NDVI were statistically significant, unlike with rainfall.
Monthly cumulative new GW rights were also correlated with average GW levels (only at a municipal level). The value was −0.70 for number of rights and −0.71 for annual average streamflow; thus, a high, significant inverse correlation was found, which reaffirms that human activity affects GW levels.

Cross Correlation
To investigate if the delayed effect of rainfall and vegetation on GW levels was important, they were cross correlated for up to six subsequent months (time lag of 1 to 6 months). In some cases, the GW-rainfall correlation values increased very little in the first months of time lag, sometimes even decreasing ( Figure 14). The GW-NDVI correlation coefficient values continued to be greater than those of GW-rainfall and were also greater than in the correlation without time lag. In summary, the correlations were slightly more influenced by the time lag, mainly GW-NDVI. The NDVI-GW correlation values were higher than the GW-rainfall correlation values, even reaching 0.49 in sector 3 (upper eastern of Petorca), and significant in most cases ( Figure 13). These GW-NDVI values are high considering that the study area is a semi-arid zone, where a low correlation between vegetation and GW is expected; this result reflects the presence of farmland. In general, the coefficients with the SPI and NDVI were statistically significant, unlike with rainfall.
Monthly cumulative new GW rights were also correlated with average GW levels (only at a municipal level). The value was −0.70 for number of rights and −0.71 for annual average streamflow; thus, a high, significant inverse correlation was found, which reaffirms that human activity affects GW levels.

Cross Correlation
To investigate if the delayed effect of rainfall and vegetation on GW levels was important, they were cross correlated for up to six subsequent months (time lag of 1 to 6 months). In some cases, the GW-rainfall correlation values increased very little in the first months of time lag, sometimes even decreasing ( Figure 14). The GW-NDVI correlation coefficient values continued to be greater than those of GW-rainfall and were also greater than in the correlation without time lag. In summary, the correlations were slightly more influenced by the time lag, mainly GW-NDVI.

Cluster Analysis
The cluster analysis grouped the wells into three groups or classes ( Table 6). The first class of five wells, which is located in areas with a large increase in avocado crop and urban areas, presented high GW decreases. The second group (12 wells) included the wells with less of a decrease in GW and their location did not show an association with land use. The third group consists of the O. Chalaco, Aconcagua, and S. Lorenzo wells ( Figure 5), which were those that presented the highest levels of GW decrease as of 2010 and which are located in large agricultural areas where avocado land increased considerably between 2002 and 2017.

Assessment of GW Depletion and its Implicants
GW depletion in the 2002-2017 period was evident when its behavior by well in the two watersheds was analyzed. A greater decrease in some wells in the Ligua central valley, a reflection of local factors, was observed; however, rainfall and the SPI results do not indicate spatial and temporal differences between the watersheds. The 2009-2015 period was the most critical and a relationship between GW-SPI was found, mainly at the 9-24-month scales, indicating a moderate hydrological drought, but the longest in the last 36 years. This drought is consistent with that addressed in the study carried out by Garreaud et al. [63] in south-central Chile, a megadrought lasting from 2010 to 2015. Other investigations have delved into the climate dynamics of this drought [16,64] and analyzed the negative impacts on the hydrological regime of Chile [16,65,66]. The most critical years were 2012 to 2014 in terms of GW, rainfall, and the NDVI and a review of the fruit

Cluster Analysis
The cluster analysis grouped the wells into three groups or classes ( Table 6). The first class of five wells, which is located in areas with a large increase in avocado crop and urban areas, presented high GW decreases. The second group (12 wells) included the wells with less of a decrease in GW and their location did not show an association with land use. The third group consists of the O. Chalaco, Aconcagua, and S. Lorenzo wells ( Figure 5), which were those that presented the highest levels of GW decrease as of 2010 and which are located in large agricultural areas where avocado land increased considerably between 2002 and 2017.

Assessment of GW Depletion and its Implicants
GW depletion in the 2002-2017 period was evident when its behavior by well in the two watersheds was analyzed. A greater decrease in some wells in the Ligua central valley, a reflection of local factors, was observed; however, rainfall and the SPI results do not indicate spatial and temporal differences between the watersheds. The 2009-2015 period was the most critical and a relationship between GW-SPI was found, mainly at the 9-24-month scales, indicating a moderate hydrological drought, but the longest in the last 36 years. This drought is consistent with that addressed in the study carried out by Garreaud et al. [63] in south-central Chile, a megadrought lasting from 2010 to 2015. Other investigations have delved into the climate dynamics of this drought [16,64] and analyzed the negative impacts on the hydrological regime of Chile [16,65,66]. The most critical years were 2012 to 2014 in terms of GW, rainfall, and the NDVI and a review of the fruit registry of the Office of Agriculture Studies and Policies (ODEPA, for its acronym in Spanish) between 2008 and 2014 confirmed that in 2014, the fruit-growing area in both watersheds had decreased, possibly due to the drought or more likely due to increases in avocado crops at the expense of other fruits; however, GW for irrigation continued to be extracted.
The results demonstrate that there is a relationship between the GW decrease and the drought, although it cannot be stated that there is a trend in rainfall behavior related to the significant GW decrease; however, the statistical correlation between GW and rainfall was very low. The standardized GW values showed spatial difference between the two watersheds, while the SPI was spatially homogeneous. In Ayala et al. [19], there was also discussion of the decrease in well levels and the difficulty of quantifying the actual amount of water level variation using historical pumping data on the area as it was undergoing irregular extraction of water resources. Celedón [20] stated that 2010 was the year with the greatest GW use for irrigation as it was the year with the greatest amount of planted farmland and lowest surface water availability. These ideas show that a large quantity of water was being extracted in the studied valleys and confirm the results.
It was possible to confirm that the agricultural factor played the most important role in aquifer depletion. The drought affected both valleys with the same intensity, but Ligua crops were maintained and avocado trees continued to be planted and were expanded into the hills (Figure 15a-c). In this, watershed irrigation and GW extraction continued despite recharge by rain decreasing, without regard for aquifer sustainability. Petorca was slightly more affected by the drought and the effects of agricultural areas on GW were lower, as can be observed in Figure 15f. This photograph was taken by the authors in the Petorca valley, where it was confirmed that many agricultural areas had been abandoned. In Figure 15d, avocado growing on a hill in Petorca is seen, showing that although some areas were lost, many were maintained. Nonetheless, water rights continued to be granted and reservoirs increased (Figure 15e). These reservoirs are located near the cultivated areas to guarantee irrigation. Many of these reservoirs were built by the Chilean government to "mitigate" drought effects, according to Bolados et.al. [14] and Panez-Pinto et al. [17], however based on the results of this investigation, this measure taken by the government could have had the opposite effect by contributing to greater groundwater extraction to keep the reservoirs supplied, considering that the Ligua and Petorca watersheds are characterized by semi-arid conditions and low annual rain, in addition to being affected by a prolonged drought process.
In addition, it was confirmed that the wells with the greatest decreases were those located in zones where avocado growing areas increased or were maintained. Avocado is a fruit with a tropical origin that requires a large water supply throughout the year for adequate vegetative and productive growth [67][68][69]. The water requirements for achieving maximal yields of this crop in a mediterranean region are over~7500 m 3 ha −1 year −1 , according to [70]; therefore, the increase or maintenance of this crop has repercussions for GW use and extraction.
Other studies have investigated the situation in the Petorca Province and Ligua watershed regarding avocado growing and water rights. In [14,15,22,23], the conflicts that were taking pace over the human right to water, agribusiness, and the usurpation of GW in Petorca were analyzed. Bolados et al. [14] carried out a study on the beginnings of avocado growing in Petorca and how the province became an exporter to large international markets such as China and Europe, as well as an analysis on the great inequalities regarding water rights stemming from the Water Code, as a result of legislation imposed during the military dictatorship. The results of these studies are reaffirmed in this work, which found a large area of avocado compared to other crops, with its expansion in valleys and hills and the continuous granting-in perpetuity and at no cost-of water rights, which were highly correlated to aquifer depletion. Although the number of rights granted decreased in recent years, probably due to measures taken by the government and the DGA, the effect is still noticeable and cumulative, with repercussions on GW depletion.
Through the analysis of data records of GW rights, it was confirmed that a large part of the rights granted belonged to entrepreneurs dedicated to the cultivation of avocado and lemon, as mentioned in Bolados et al. [14].
was taken by the authors in the Petorca valley, where it was confirmed that many agricultural areas had been abandoned. In Figure 15d, avocado growing on a hill in Petorca is seen, showing that although some areas were lost, many were maintained. Nonetheless, water rights continued to be granted and reservoirs increased (Figure 15e). These reservoirs are located near the cultivated areas to guarantee irrigation. Many of these reservoirs were built by the Chilean government to "mitigate" drought effects, according to Bolados et.al. [14] and Panez-Pinto et al. [17], however based on the results of this investigation, this measure taken by the government could have had the opposite effect by contributing to greater groundwater extraction to keep the reservoirs supplied, considering that the Ligua and Petorca watersheds are characterized by semi-arid conditions and low annual rain, in addition to being affected by a prolonged drought process.

GW Resources Management in the Ligua and Petorca Watersheds
GW is a resource for current and future generations and its sustainable use is possible only through sustainable GW resources management, as explained in [1,13]. The results showed GW overexploitation as there is an imbalance between availability and extraction. Therefore, GW management in the Ligua and Petorca watersheds is not sustainable. For several decades, with the beginning of agribusiness and the privatization of water in the area [14], water resources management has been absent despite great efforts made by inhabitants and social organizations linked to the human right to water such as the Movimiento de Defensa del Agua, la Tierra y el Medioambiente (Water, Land and Environmental Defense Movement; MODATIMA, for its acronym in Spanish), as well as the DGA. This agency has used various regulatory mechanisms such as the creation of Groundwater User Communities and declarations of water scarcity and restriction zones [16,20,21], but they have not been effective.
It is important to establish scientific and political agreements between academic and governmental agencies inside and outside the watersheds. A new reform of the Water Code must include environmental streamflows and new, effective regulations that guarantee that the supply capacity of aquifers does not exceed their sustainable volume. New water management mechanisms that ensure greater power for regions and local authorities as regulating agents and improvements in agricultural practices amid drought scenarios are necessary. The Chilean government and the DGA must implement effective actions to equitably guarantee the human right to water. There also must be more inspections of water use and sanctions in the event of violations of groundwater rights. Additionally, GW communities must be supported to guarantee increased user interest and greater functioning of these already established organizations, since they were created to achieve GW management, according to [71] and [23]. All these proposed measures can contribute to GW resources management.
To continue delving into this issue of GW loss, more interdisciplinary studies that include citizen perceptions are recommended to contribute to decision-making. Providing these results to local, regional, and national governments could improve the water management. In addition, further in-depth study of land-use dynamics before, during, and after the studied period is needed.

Conclusions
The GW level presented a statistically significant decreasing trend, which was greater in the Ligua watershed wells, while the decrease in rainfall was not statistically significant. Starting with SPI 6, a drought period between 2010 and 2015 was observed, which was more evident in SPI 12 and SPI 24, which showed a hydrological drought classified as moderate but very long. A relationship between rainfall and GW was demonstrated, but the SPI and rainfall analysis did not show spatial differences between the two watersheds. The statistical analysis offered an accurate view of the relationship between variables. The GW-SPI (12 and 24) correlations were high, which showed a relationship between GW depletion and hydrological drought. However, GW-rainfall correlations were low; therefore, although rainfall is the main recharge source of GW in the study area and aquifer depletion was mostly influenced by human factors. In addition, the GW-NDVI correlation was greater than that of GW-rainfall, reflecting the influence of irrigation, and the correlation between GW and GW rights was very high. The land-use change analysis allowed the conclusion to be drawn that the areas with the greatest GW decreases were those in which avocado growing increased. The Ligua watershed was more affected, with an increase of 2623.6 ha and a greater GW decrease, while in Petorca, the avocado area decreased very little (128 ha); however, it also presented areas in which the area dedicated to this crop increased or was maintained, along with an increase in reservoirs.
The increase in high-water demand or the overgranting of GW rights, associated with overexploitation by agriculture due to increases in the avocado-growing area, irrigation, and reservoirs, along with the lack of sustainable water resources management, confirmed that groundwater depletion was mostly influenced by human factors regarding climatic factors.
The management measures proposed and the results obtained through an interdisciplinary approach, which reflect water systems that cannot support the high demand of agriculture in semi-arid conditions, the impacts of drought, and the need for more sustainable GW resources management, all against the backdrop of a private water management model, are the main lessons to be learned from this research.