MONITORING GROUNDWATER STORAGE IN NORTHERN CHILE BASED ON SATELLITE OBSERVATIONS AND DATA SIMULATION

Groundwater is one of the most valuable sources of fresh water in many places worldwide, especially in regions with low pluviometric indices such as northern Chile. Thus, it is mandatory to monitor this precious resource in space and time domains. Currently, groundwater in Chile is monitored using sparse stations of water table observations. Although other indirect alternatives such as space-borne observations can contribute to regional understanding of groundwater variations, they have been poorly studied in Chile. In this study, groundwater monitoring is carried out based on 104 monthly solutions of the Gravity Recovery and Climate Experiment (GRACE) mission between 2004 and 2013. The extraction of the groundwater storage (GWS) signal obtained from GRACE was recovered once the effects of soil moisture and snow storage, retrieved from the Global Land Data Assimilation System (GLDAS), were removed. Analysis of the data was performed point-wise (six stations) and at regional scale (Northern Chile). Overall, the results are correlated with wells observations obtained by the General Directorate of Water Resources (DGA) of the Ministry of Public Works of Chile. Point-wise comparison shows root mean square error (RMSE) large than 30.0 mm while regional scale validation shows RMSE of 21.5 mm. Furthermore, regional groundwater variations obtained from GRACE/GLDAS are highly consistent in terms of trend with results obtained from well observations in the DGA network. The Empirical Orthogonal Function (EOF) analysis revealed higher annual groundwater variability in the metropolitan region and a higher inter-annual variability in the north. The methodology used may contribute to the regional study of spatial-temporal variations of groundwater in regions with sparse hydrometric network. Montecino,H.C. et al. 2 Bol. Ciênc. Geod., sec. Artigos, Curitiba, v. 22, n1, p.01-15, jan-mar, 2016.


Introduction
It has been reported in the early 1990s that Chile's sustainable development would depend on the strategic use of groundwater resources (Peña et al., 2004).The main advantage identified in the Chilean context however is that the groundwater generally has lower variability compared to surface water resources (e.g., rivers and lakes).It also presents less risk of contamination, appropriate chemical and bacteriological properties, and exploitation near consumption centers is generally feasible.In the area from the Metropolitan Region of Santiago to the north, groundwater resources have vital importance, since precipitation decreases (see, e.g., Ribera and Lucero, 2006;Squeo et al., 2006).It has to be noted that periods of 20 years have occurred without rainfall in the city of Iquique, Region of Tarapaca (Ribera and Lucero, 2006).For this reason, knowledge about the aquifer formations and disputes arising from their usage is greater at north of Santiago than to the south.Groundwater assessment, modeling and management are difficult due to the lack of in-situ observations, especially in arid and semi-arid zones where only sparse surveillance infrastructure exists (see, e.g., Brunner et al., 2007).In these environments, generally only a limited number of in-situ/well measurements are available, whereas groundwater models need temporal and spatial distribution for input data and calibration.If such Monitoring ground water... Bol.Ciênc.Geod., sec.Artigos, Curitiba, v. 22, n o 1, p. 01-15, jan-mar, 2016.data are not available, the models cannot play a relevant role in decision support, because under these circumstances they are inaccurate (see, e.g., Brunner et al., 2007;Becker, 2006).Methods for monitoring groundwater are traditionally based on in-situ data (e.g., water wells and lysimeters), which involve high costs (Tuinhof et al., 2004).Alternatively, some researches related to the study of groundwater using remotely sensed data have been reported (see, e.g., Meijerink, 1996;Travaglia and Ammar, 1998;Jha et al., 2007;Strassberg et al., 2009;Kolokoussis et al. 2011;Frappart et al., 2011;Muskett and Romanovsky, 2011;Boronina and Ramillien, 2008;Chinnasamy et al., 2013).Thereby, these investigations rely essentially on the qualitative characterization of hydro-geological mapping units and detection of characteristic features.Most applications refer to the crystalline subsoil, limestone and quaternary volcanic terrain.Recent developments are related to the occurrence of groundwater discharge in areas of groundwater flow systems using thermal and multispectral images (see, e.g., Bobba, et al., 1992;Meijerink, 1996;Banks et al., 1996;Jackson 2002;Mutiti et al., 2010).Nevertheless, satellite gravimetry is the space-borne geodetic technique most recently used to assess changes in the total water storage and its changes (Cazenave and Chen, 2010).As it is well known, the composition and structure of the Earth, including the distribution of the atmosphere, surface and under groundwater masses, is reflected in the Earth's gravity field (Tapley et al., 2004).The launch of GRACE in March 2002 opened a wide range of possibilities for studies in the geosciences and hydrology, due to the recovery of temporal variations of the gravity field.Several studies of groundwater using GRACE observations have already been reported (see, e.g., Rodell et al., 2009;Henry et al., 2011;Muskett and Romanovsky, 2011;Scanlon et al., 2012;Jin and Feng 2013;Awange et al., 2014;Joodaki et al., 2014;Forootan, et al., 2014;Xiao et al., 2015).However, as far as we know none of them refers to the Chilean situation, although the lack of research on these precious resources was pointed out by McPhee et al. (2012).Thus, the main aim of this study is to assess the groundwater storage over Chile using GRACE satellite mission.However, the content of the recovered signal by GRACE satellites contains changes associated with gravity changes coupled to changes in water storage in the atmosphere, surface water, soil moisture and groundwater.Therefore, in the case of groundwater studies it is necessary to isolate the other water components.One possibility to filter out the effects of changes in water storage in the atmosphere, surface water and soil moisture is through hydrological surface models such as the Global Land Data Assimilation System (GLDAS) (Rodell, et al., 2004).From GRACE and GLDAS data, we estimated and characterized the spatial and temporal variations of groundwater in northern Chile since there are no significant surface water reservoirs.Additionally, we performed the validation of GRACE-based groundwater changes using in-situ water table level measurements provided by the General Directorate of Water Resources of Chile (Dirección General de Aguas -DGA, Chile).

GRACE Level 2 Products
Monthly solutions acquired by GRACE from January 2004 to December 2012, corresponding to 104 monthly spherical harmonic coefficients (Level 2 products) release 05 (RL05), up to degree and order (d/o) 60 were applied.The set of spherical harmonic coefficients were obtained from GeoForschungsZentrum (GFZ) (Dahle et al. 2013) University of Texas at Austin (Bettadpur 2012), and Jet Propulsion Laboratory (JPL) center (Watkins 2012).Since 2011, data gaps occur in the GRACE K-Band data due to the aging of the batteries.Specifically for the periods of 01/2011, 06/2011, 05/2012, and 10/2012.For the months without data, linear interpolation was performed by considering the previous GRACE derived terrestrial water storage anomalies.The coefficients , and were replaced from those obtained from Satellite Laser Ranging (SLR).Here, we use as a solution the average of the solutions of the three processing centers GFZ, JPL, and CSR.Average deviations between the different solutions reach less than 15 mm.

Global Land Data Assimilation System -GLDAS
GLDAS drives multiple, offline (not coupled to the atmosphere) land surface models, integrates a huge quantity of observation based data, and executes globally at 2.5° to 1 km resolutions, enabled by the Land Information System (LIS) (see Kumar et al., 2006).Currently, GLDAS drives four land surface models: Mosaic, Noah, the Community Land Model (CLM), and the Variable Infiltration Capacity (VIC).GLDAS products include land surface state (e.g., soil moisture and surface temperature) and flux (e.g., evaporation and sensible heat flux) parameters.The temporal resolution for the GLDAS products is 3 hours.Monthly products are also generated through temporal averaging of the 3-hourly products (Fang et al., 2009).GLDAS water content (soils moisture and snow storages) correspondent to the solutions from the Noah (Ek, et al., 2003) LSM2.7.1 were used in this study.The GLDAS fields cover the same time span as the GRACE observations.The parameters contained in the GLDAS data refer to all soil moisture layers and the snow storage component (canopy water is negligible over the study region).

DGA Groundwater stations
A small network of six hydrometric stations extended over a period from 2004 to 2013 controlled by the DGA were used (Fig. 1) in order to assess the quality of the analysis based on satellite and model datasets.The selection of these stations is justified by their particular characteristics and availability of data in the study period.The Arica station has particular characteristics in terms of tectonic and hydrologic streams.The Alana station is located in the basin of the Salar de Atacama.The Atacama station is located in the Copiapó River basin, which anthropogenic activities such as industrial, mining, fishing, and agriculture are taking place.The location of the Coquimbo station is near the Elqui River basin, which is essential for the development of agricultural activities in the region.Additionally, Valpo and R.M. stations are located near to important watersheds: the basin of the rivers Petorca and La Ligua, in the case of Valpo station; and associated watersheds Chacabuco systems, Colina, Rio Mapocho, Melipilla and Puange in the case of the R.M. station.The time series of in-situ groundwater observations were examined in order to remove gross errors were removed, for example, variations due to change of limnimetric ruler.Additionally, in order to harmonize the in-situ data with GRACE results, the raw well depth observation was converted into changes of water table heights in relation to the 2004 to 2013 mean for that well location.The groundwater storage (GWS) changes as derived from GRACE are parameterized in units of equivalent water heights.Under the assumption that the well measurements represent unconfined, static water table condition, the GWS for each well can be calculated as (Telford et al., 1990): where is the change in water table elevation, is the specific yield, and is the thickness of an equivalent layer of water that is required to produce the measured change in gravity.It has to be mentioned that, according to the soil type in each study region, the specific yield values were applied and not empirical values.The specific yield values associated to each well were obtained from values published in Heath (1983).The regional averaged groundwater storage anomalies (in terms of Equivalent Water Heights -EWH) can be calculated by (e.g., Xiao et al., 2015): where refers to the number of subareas or zones divided in the study region; are the specific yield values of the unconfined aquifers; are the sizes of subareas; and refer to the mean values of the well or GRACE/GLDAS water-level variations in each subarea.The areas of influence related to each station ( ) were obtained by the Thiessen polygons algorithm (e.g., Raghunath, 2006).

Empirical Orthogonal Functions (EOF)
In order to identify and quantify the predominant spatio-temporal signal patterns in groundwater, an analysis by Empirical Ortogonal Functions (EOF) was performed in the study region.Other studies (see, e.g., Winter et al., 2000) had indicated that EOF analysis would be particularly useful for long-term groundwater monitoring, which would greatly reduce the cost of monitoring programs.In short, EOF analysis uses a set of orthogonal functions (i.e., EOFs) for representing a time series as: where is the original time series as a function of time ( ) and space ( ) of the major factors that account for the temporal variations of while the principal components PC shows how the amplitudes of each EOF varies with time.To address some limitations of classical EOF analysis, such as the difficulty of physical interpretability, it is recommended rotating the EOFs (Hannachi, 2004).Thus, in this study we have applied the rotated EOFs (i.e., REOFs).

Groundwater anomalies from GRACE and GLDAS
Mass changes detected by GRACE (ΔMGRACE) are related to changes of hydrological mass (ΔMHydrologic) and mass changes provided by other causes (ΔMother).Considering that ΔMHydrologic is composed by the contribution of changes in surface mass water (ΔMSW), soil moisture (ΔMSM), and groundwater (ΔMGW), it can be defined as (see, e.g., Longuevergne et al., 2012;Tapley et al., 2004): where is the precipitation, is the evaporation (describing all process of vaporization), is the lateral transport of water (generally equivalent to discharge), and are the first and last days of the month, respectively.In case that the region under investigation does not contain any major surface water bodies, as in Northern Chile, (4) can be simplified as follow: The term ΔMGRACE is the residual variation of the hydrological mass, obtained as the difference of mean spherical harmonic coefficients with respect to the monthly fields at each epoch (Ramillien et al., 2005).Based on the methodology proposed by Wahr, et al. (1998), it can be expressed as: where ΔClm and ΔSlm are the residual of the spherical harmonic coefficients (i.e., the mean value over the time period has been subtracted from the monthly values).ρEarth is the average density of the Earth, ρwater is the density of fresh water, a is the equatorial radius of the Earth, are the associated Legendre completely normalized functions of degree l and order m, kl is the Love number of degree l, θ is the spherical co-latitude, and λ is the longitude.
Terrestrial water storage anomalies derived from GRACE measurements were filtered by a destriping filter proposed by Swenson and Wahr (2006), which is based on a polynomial spline for consecutive coefficients, which can be odd or even, for a particular degree.The correlation between the coefficients is represented by a polynomial fit, and is removed from every coefficient.In addition to this filter, a Gaussian filter of 350 km was used to remove the remaining errors at short wavelengths (Swenson and Wahr, 2006).The filters previously applied to GRACE data reduce noise and remove some correlated errors, but also cause loss of signal strength.In order to restore the geophysical signal, a scaling coefficient was applied on the equivalent water storage according to Landerer and Swenson (2012).

Comparison of groundwater estimations
In order to investigate the relative difference between the groundwater time series, the results of groundwater derived from GRACE were analyzed for each DGA station (i.e., Arica, Alana, Atacama, Coquimbo, Valpo, and R.M). Figure 2 shows the behavior of the variations of groundwater for each station.It is possible to see that the time series of Coquimbo (d), Valpo (e), and R.M. (f) reflect the impact of drought episodes over the periods of 2007, 2008, 2010, 2011 and 2012 on groundwater variations.The Region of Coquimbo, for example, is characterized by a strong seasonal semi-arid climate, it is strongly affected by inter-annual variability, and high evaporative demands and variable rains make droughts a recurrent phenomenon (Meza, 2013).Alana (Fig. 2-b) and Atacama (Fig. 2-c) stations had the lowest and highest, respectively, dispersion among the solutions obtained from the CSR, GFZ and JPL processing centers as shown by their respectively error bars.
In order to have an overview and to contrast with GRACE/GLDAS estimations of GWS, from the metropolitan region to the north of Chile, a time series was constructed from 2004-2013 (Fig. 3) based on the six well stations data averaged through Eq. ( 2).It is possible to see the increase tendency of both time series.Moreover, in the period 2008.5-2013.0,both data sources are highly consistent, present trends of 1.5 mm/yr for GRACE and 2.9 mm/yr for well-based observation.The relative comparison between GRACE-based GWS and well observation time series at the specific sites (Fig. 1) shows root mean square error (RMSE) of 44.25 mm for Arica, 30.53 mm for Alana, 58.10 mm for Atacama, 32.24 mm for Coquimbo, 60.79 mm for Valpo, and 76.42 mm for R.M.For instance, Valpo station show a large RMSE while the other five stations are within the results presented by other authors.Nevertheless, the finite number of spherical harmonic coefficients, as provided by GRACE, prohibits the estimation of TWS anomalies at a point, but regional averages should be made (see Swenson and Wahr, 2002).Thus, considering the average of all wells stations, as shown in Fig. 3 computed by Eq. ( 2), the RMSE is 21.46 mm.In Xiao et al. (2015) it has been reported an RMSE of 26 mm over the Mid-Atlantic Region of the United States.The large RMSEs at specific sites and regional scales comparisons could be explained by the uncertainties of GRACE-derived TWS, especially at higher degrees and the uncertainties related to the hydrologic models as well (i.e., GLDAS) relatively to the other six stations.Additionally, this could also be due to the sparse temporal sampling of the ground measurements, and difference in spatial scales represented by the GRACE and GLDAS data.The trend estimates obtained from GRACE are quite consistent with well observations (Table 1).However, in some cases there are significant differences in terms of magnitude.These differences might be due to local effects not perceived by GRACE or errors such as sensor replacements at the DGA wells and gross errors as well (see Figure 2-c, the possible outliers between 2008 and 2011).Overall, the results of groundwater temporal variations from GRACEand well-based observations in the study period reveal a mass gain of groundwater in Arica and Alana stations, and a loss in the Atacama, Coquimbo, and R.M. stations (Table 1).For instance, at Valpo station, the results obtained from GRACE and well observations show a positive and negative trend of 2.3 and -8.4 mm/yr for in-situ and GRACE-derived GWS, respectively.However, during the period from September 2008 to December 2012 a increase in the mass gain reaching ~11 mm/year is observed (Fig. 2-b).This effect reflected in this period could be due to a lack of rainfall and a high evaporation rate.Moreover, our estimation reveals a different groundwater behavior (Table 1) in terms of seasonal variations.Estimates of GWS obtained from GRACE in the stations and regional analysis reveal a moderate consistency with the wells observations in terms of trend.Therefore, it should be noted some of the differences between the time series of GWS obtained from GRACE and wells, may be associated with particular characteristics (e.g.basins, sub-basins and rivers) located very close to each other, and are not recovered by the limited resolution of satellite data.Additionally, at regional scale, the Northern Chile ground water storage shows an overall mass gain during early 2004 to middle 2012 of 0.6 mm/yr for GRACE and 5.7 mm/yr for in-situ data.Despite this difference in terms of linear trend, it is possible to see the overall agreement in terms of patterns between the two time series, for example, the periods 2004-2007 and 2008-2012 are clearly depicted in both.Considering the period of 2007 afterwards, the mass gain is approximately 2 mm/yr for GRACE-based GWS and 6 mm/yr for in-situ measurements.This differences call for future investigations, which will be the next step of the current study.
Table 1: Linear trend and its standard deviation (mm/yr)) of the time six stations' series and we as for the averaged values as per Eq. ( 2) for GRACE-derived groundwater storage anomalies and in-situ data.

EOF analysis of groundwater anomalies
In order to reveal patterns at regional scale and for a better physical interpretation of the spatiotemporal of GWS, Rotated EOF (REOF) were used (Fig. 4).A set of 104 monthly GWS grids with a spacing of 0.5° were used for the EOF analysis.Most of the signal variances of GWS are explained by the first three modes, which show annual and semi-annual patterns, these modes account for 60.7% ( ), 13.9% ( ), and 11.2% ( ).The Principal Components Rotated 1 (PCR1, see Fig. 4 -a) is related to the annual cycle, and exhibits low variability during the period of early 2004 to latter 2006, along with a negative trend.Starting on early 2007, PRC1 of GWS presents a positive trend with high variability (Fig. 4-a).The PCR2 (Fig. 4-b) is related to inter-annual variability effects with more variability than the PCR1.The PCR3 (Fig. 4-c) also shows annual fluctuations, but with overlapping inter-annual effects, in addition it presents a negative slope until early 2008, and a positive trend from 2008 to 2013.In relation to spatial variations (Fig. 4-d), the REOF1 (Fig. 4-a) corresponding to annual variations, these have the greatest negative variability in the Santiago Metropolitan Region, a slight positive variability in Arica, and slightly between Alana and Coquimbo.It is also possible to see that the center is located out of Chile (located at Bolivia and Argentina) and it seems that the patterns are influenced by the local factors such as orography.The REOF2 (Fig. 4-e) shows a significant positive variation from the north of R.M. On the other hand, REOF3 shows a negative variability in Valpo and Coquimbo regions (Fig. 4-f).However, with the center of EOF2 out of Chile while the center of EOF3 over Santiago Metropolitan Region.

Conclusions and Outlook
Groundwater storage (GWS) changes based on GRACE and GLDAS for northern Chile is presented for the first time.In this sense, the GWS were isolated from the GRACE-derived terrestrial water storage (TWS) by considering the groundwater, soil moisture and snow as the only significant contributors to the regional water resources.Hence, the groundwater variations were isolated by subtracting the GLDAS-Noah simulated soil moisture and snow storages from GRACE-derived TWS.Additionally, water level records were employed for validation and evaluation purposes.The time series of the groundwater stations show particular tendencies depending on the period, and the specific characteristics (e.g., river basins) of each location.The behavior of the GWS obtained from GRACE/GLDAS in relation to well observations is moderately consistent in terms of trends in the stations and regional scale.The EOF analysis indicates that the largest annual spatial variations over the time window of 2004-2013 are experienced in the metropolitan region, and decreases to the north.However, inter annual variations are seen from the metropolitan region to the north.Differences between GRACE solutions and wells observations might be the result of deficiencies in GRACE and/or the GLDAS hydrological model that does not allow considering short term events.Nevertheless, our study shows that GRACE-based GWS estimations depicted the droughts episodes in the years of 2007, 2008, 2010, 2011, and 2012.Thus, the methodology used may contribute to regional study of temporal variations of groundwater in regions with sparse hydrometric network.As a future work, specific models and/or data that include soil moisture (e.g., Soil Moisture and Ocean Salinity -SMOSsatellite mission) in order to study the spatial distribution of recharge and discharge of major aquifers in northern Chile will be considered.Additionally, optical satellite imagery (e.g., Landsat 7, Landsat 8) will be analyzed with the aim to study vegetation growth and health as well as their correlations with groundwater changes.

Figure 1 :
Figure 1: Study region and well stations.
Additionally, the amplitudes and patterns of both signals are quite different.GRACE-based GWS shows much larger water increases in the austral spring since 2004 until 2008, while it starts to show a water level decreasing until middle of 2012.This period of fall (2008-2012) is reported as period of drought that has blighted Chile since 2007.The drought in Chile, which began in 2007, has hampered copper production, exacerbated forest fires, and affected agriculture activities.During spring 2007 until winter of 2011, the well-based estimates and those from GRACE GWS show almost constant water level, which seems to be due to the drought episodes.

Figure 2 :
Figure 2: Time series of in-situ and GRACE-based groundwater storage changes obtained for Arica (a), Alana (b), Atacama (c), Coquimbo (d), Valpo (e), and R.M. (f).The error bars represent the standard deviation among GRACE processing centers.

Figure 3 :
Figure 3: Monthly groundwater storage anomalies for Northern Chile based on six wells and GRACE/GLDAS data.