**Contribution of GRACE Satellite Gravimetry in Global and Regional Hydrology, and in Ice Sheets Mass Balance**

Frappart Frédéric and Ramillien Guillaume *Université de Toulouse, OMP-GET, CNRS, IRD, France* 

#### **1. Introduction**

190 Water Resources Management and Modeling

Wiener, J.G., & Spry, D.J. (1996). Toxicological significance of mercury in freshwater fish.

128p.

Publishers.

(1998). Technical protocol for evaluating natural attenuation of chlorinated solvents in groundwater. EPA/600/R-98/128. USEPA, Washington DC, USA.

In *Environmental Contaminants in Wildlife: Interpreting Tissue Concentrations*, ed. W. N. Beyer, G. H. Heinz & A. W. Redmon, 297-339. Boca Raton, FL: Lewis

> Continental water storage is a key component of global hydrological cycle which play a major role in the Earth's climate system via controls over water, energy and biogeochemical fluxes. In spite of its importance, the total continental water storage is not well-known at regional and global scales because of the lack of *in situ* observations and systematic monitoring of the terrestrial water reservoirs, especially the groundwaters component (Alsdorf & Lettenmaier, 2003). Although some local hydrological monitoring networks exist, the current description of water storage comes primarily from global hydrology models (*e.g.*, WaterGAP Global Hydrology Model - WGHM (Döll et al., 2003) or Global Land Data Assimilation System - GLDAS (Rodell et al., 2004a)), which still suffer from important limits, such as the absence of different compartments of the terrestrial water compartments (*e.g.*, groundwater and floodplains), the errors introduced by external forcings (mostly the precipitation), the lack of reliability of external parameters as the soil type or the vegetation cover, the parameterization. Unfortunately, direct hydrological measurements are fairly limited, while imaging satellite techniques and satellite altimetry only give access to surface water and superficial soil moisture variations, that represent just the more accessible components of total water storage. Since its launch in 2002, the Gravity Recovery and Climate Experiment (GRACE) gravimetry from space mission offers a unique alternative to classical remote sensing technique for measuring changes in total water storage (ice, snow, surface waters, soil moisture, groundwater) over continental areas, representing a new source of information for hydrologists and the global hydrological modelling community.

> We propose a review of the major techniques and results in extracting continental hydrology signals from GRACE satellite gravimetry, aiming at studying large-scale processes in hydrology and glaciology and the impact of the climatic variations on the global hydrological cycle. The two first parts describe the principle of the GRACE mission and the nature of the data measured. The third part is devoted to the post-processings applied to the GRACE-based products to improve their quality and accuracy. The last part presents the most significant results obtained for land hydrology and ice sheets mass balance using GRACE.

#### **2. Principle of GRACE space gravimetry mission**

GRACE, a joint US-German (NASA/DLR) mission, was launched in March 2002 and placed at low altitude on a quasi-polar orbit (inclination of 89.5°), to map the time variations of the Earth's gravity field.

#### **2.1 Description of the "Low-Low" time-variable gravity measurements**

The mission consists of two identical satellites moving along the same orbit at an altitude of ~ 450 km, one following the other with a mean inter-satellite distance of approximately 220 km in a satellite-to-satellite tracking configuration in the low-low mode initially proposed by Wolf (1969). The satellites use K-Band microwave ranging (KBR) system (Tapley et al., 2004a) to monitor continuously the relative motion of the spacecraft (*i.e.*, the range and the range rate), which varies proportionally to the integrated differences of the gravity accelerations felt by the satellites (Schmidt et al., 2008). The KBR system provides measurements of the distance between the two satellites with an accuracy better than 1 micrometre using carrier phase measurements in the K (26 GHz) and Ka (32 GHz) frequencies. Each GRACE satellite contains a 3-axis accelerometer which gives access to the dynamical effects of the non-conservative forces such as solar pressure and atmospheric drag. Once removed the non-gravitational effects from the inter-satellite range, the corrected variations of distance are used to estimate the geopotential with an unprecedented accuracy.

#### **2.2 Time-variations of mass in the Earth system**

Models of the Earth's geopotential are classically obtained by inverting GRACE data at monthly timescales and a spatial resolution of a few hundred kilometres (Tapley et al., 2004b). For the very first time, GRACE has enabled the monitoring of the spatiotemporal variations of the mass in the Earth system. The tiny variations of the Earth gravity field are due to the redistribution of fluid mass inside the surface fluid envelopes of the planet (atmosphere, oceans, continental water storage) (Dickey et al., 1999).

#### **2.3 Expected accuracy and resolution for continental hydrology**

The main application of GRACE is quantifying the terrestrial hydrological cycle through measurements of vertically-integrated water mass changes inside aquifers, soil, surface reservoirs and snow pack, with a precision of a few millimeters in terms of water height and a spatial resolution of ~ 400 km (Wahr et al., 1998; Rodell & Famiglietti, 1999). Pre-launch studies focused on the ability of GRACE to provide realistic hydrological signals, especially on continents. Wahr et al., (1998) highlighted the need to consider short-wavelength noise and leakage errors, as well as the size of the river basin. This early work largely inspired subsequent studies addressing the accuracy of recovering continental water mass signals from synthetic GRACE geoids (Rodell and Famiglietti 1999, 2001, 2002; Swenson and Wahr, 2002; Ramillien et al., 2004). The results of these simulations suggest that final accuracy would increase with spatial and temporal scales. The comparison of a large a number of modelled outputs of Terrestrial Water Storage (TWS) and the expected GRACE measurements accuracy showed that water storage changes would be detectable at spatial scales greater than 200,000 km2, at monthly and longer timescales, and with monthly accuracies of roughly 1.5 cm (Rodell & Famiglietti, 1999). Similar conclusions were obtained using a network of hydrological observations of snow, surface water, soil moisture (SM) and groundwater (GW) in Illinois (Rodell & Famiglietti, 2001).

At basin-scale, the accuracy of GRACE measurements is expected to be 0.7 cm equivalent water height (EWH) for a basin with an area of 0.4 x 106 km², and about 0.3 cm EWT for a basin with an area of 3.9 x 106 km² (Swenson et al., 2003).

Combined with external information, GRACE offers the potentiality to monitor large-scale groundwater changes. Using a dense network of soil moisture data for the High Plains aquifer in the Central USA, Rodell and Famiglietti (2002) demonstrated the feasibility to dectect changes in aquifers since the uncertainties of mean variations of GRACE-derived GW for the ~ 450,000 km² aquifer was approximately 8.7 mm, whereas the amplitudes of the groundwater storage change signals were 20 and 45 mm for annual and 4-year periods respectively.

### **3. The GRACE data**

192 Water Resources Management and Modeling

GRACE, a joint US-German (NASA/DLR) mission, was launched in March 2002 and placed at low altitude on a quasi-polar orbit (inclination of 89.5°), to map the time variations of the

The mission consists of two identical satellites moving along the same orbit at an altitude of ~ 450 km, one following the other with a mean inter-satellite distance of approximately 220 km in a satellite-to-satellite tracking configuration in the low-low mode initially proposed by Wolf (1969). The satellites use K-Band microwave ranging (KBR) system (Tapley et al., 2004a) to monitor continuously the relative motion of the spacecraft (*i.e.*, the range and the range rate), which varies proportionally to the integrated differences of the gravity accelerations felt by the satellites (Schmidt et al., 2008). The KBR system provides measurements of the distance between the two satellites with an accuracy better than 1 micrometre using carrier phase measurements in the K (26 GHz) and Ka (32 GHz) frequencies. Each GRACE satellite contains a 3-axis accelerometer which gives access to the dynamical effects of the non-conservative forces such as solar pressure and atmospheric drag. Once removed the non-gravitational effects from the inter-satellite range, the corrected variations of distance are used to estimate the geopotential with an unprecedented accuracy.

Models of the Earth's geopotential are classically obtained by inverting GRACE data at monthly timescales and a spatial resolution of a few hundred kilometres (Tapley et al., 2004b). For the very first time, GRACE has enabled the monitoring of the spatiotemporal variations of the mass in the Earth system. The tiny variations of the Earth gravity field are due to the redistribution of fluid mass inside the surface fluid envelopes of the planet

The main application of GRACE is quantifying the terrestrial hydrological cycle through measurements of vertically-integrated water mass changes inside aquifers, soil, surface reservoirs and snow pack, with a precision of a few millimeters in terms of water height and a spatial resolution of ~ 400 km (Wahr et al., 1998; Rodell & Famiglietti, 1999). Pre-launch studies focused on the ability of GRACE to provide realistic hydrological signals, especially on continents. Wahr et al., (1998) highlighted the need to consider short-wavelength noise and leakage errors, as well as the size of the river basin. This early work largely inspired subsequent studies addressing the accuracy of recovering continental water mass signals from synthetic GRACE geoids (Rodell and Famiglietti 1999, 2001, 2002; Swenson and Wahr, 2002; Ramillien et al., 2004). The results of these simulations suggest that final accuracy would increase with spatial and temporal scales. The comparison of a large a number of modelled outputs of Terrestrial Water Storage (TWS) and the expected GRACE measurements accuracy showed that water storage changes would be detectable at spatial scales greater than 200,000 km2, at monthly and longer timescales, and with monthly accuracies of roughly 1.5 cm (Rodell & Famiglietti, 1999). Similar conclusions were obtained

**2.1 Description of the "Low-Low" time-variable gravity measurements** 

**2. Principle of GRACE space gravimetry mission** 

**2.2 Time-variations of mass in the Earth system** 

(atmosphere, oceans, continental water storage) (Dickey et al., 1999).

**2.3 Expected accuracy and resolution for continental hydrology** 

Earth's gravity field.

Three processing centers including the Centre for Space Research (CSR), Austin, Texas, USA, the GeoForschungs Zentrum (GFZ), Potsdam, Germany, and Jet Propulsion Laboratory (JPL), Pasadena, California, USA), and forming the Science Data Center (SDS) are in charge of the processing of the GRACE data and the production of Level-1 and Level-2 GRACE products. They are distribruted by the GFZ's Integrated System Data Center (ISDC – http://isdc.gfzpotsdam.de) and the JPL's Physical Oceanography Distributive Active Data Center (PODAAC – http://podaac-wwww.jpl.nasa.gov). Preprocessing of Level-1 GRACE data (*i.e.*, positions and velocities measured by GPS, accelerometer data, and KBR inter-satellite measurements) is routinely made by the SDS, as well as monthly GRACE global gravity solutions (Level-2).

#### **3.1 Level-2 GRACE products**

The Level-2 GRACE gravity solutions consist of monthy time-series of Stokes coefficients (*i.e.*, dimensionless spherical harmonics coefficients of the geopotential) developed up to a degree between 50 and 120, that are adjusted for each monthly time-period from raw alongtrack GRACE measurements. Formal errors associated with the estimated Stokes coefficients are provided to the users. A dynamic approach, based on the Newtonian formulation of the satellite's equation of motion in an inertial frame centered at the Earth's center combined with a dedicated modelling of gravitational and non-conservative forces acting on the spacecraft, is used to compute the monthly GRACE geoids (Schmidt et al., 2008). During the process, atmospheric and ocean barotropic redistribution of masses are removed from the GRACE coefficients using ECMWF and NCEP reanalyses for atmospheric mass variations and ocean tides respectively, and global circulation oceanic models (Bettadpur, 2007; Fletchner 2007). The GRACE coefficients are hence residual values that should represent continental water storage, errors from the correcting models, and noise over land. The monthly solutions differ from one group to another due to differences in the data processing, the choice of the correcting models, and the data selection during the monthly time-period. Four different completely reprocessed releases were distributed yet presenting a continuous improvement in the quality of the GRACE data. Some other research groups also propose alternate processing of the GRACE data at monthly or sub-monthly time-scales using global (harmonics coefficients) and regional approaches. They include the Groupe de Recherche en Géodésie Spatiale (GRGS), Toulouse, France, the Institute of Theoretical Geodesy (ITG) at the University of Bonn, Germany, the Delft Institute of Earth Observations and Space Systems (DEOS), Delft, Nederlands, the Goddard Space Flight Center (NASA GSFC), Greenbelt, Maryland, USA.

#### **3.2 Conversion into coefficients of water mass anomalies**

Variations in continental water storage cause changes of mass distribution in the Earth system, and thus spatiotemporal changes of the geoid, defined as the gravitational equipotential surface that best coincides with the mean sea surface. In global representation, it can be expressed as (Heiskanen & Moritz, 1967):

$$
\begin{Bmatrix}
\delta \mathbf{C}\_{nm}(t) \\
\delta \mathbf{S}\_{nm}(t)
\end{Bmatrix} = \frac{1+k\_n^{\cdot}}{(2n+1)} \frac{R\_\varepsilon}{M} \begin{Bmatrix}
\mathbb{I}\_\mathbb{S} \ \delta \mathbf{q}(\theta, \mathcal{X}, t)
\end{Bmatrix} \begin{Bmatrix}
\cos(m\mathcal{X}) \\
\sin(m\mathcal{X})
\end{Bmatrix} P\_{nm}(\cos\theta) \delta \mathbf{S} \tag{1}
$$

for a given time *t*, where *q* is the surface distribution of mass, *M* and *S* are the total mass and the surface of the Earth, respectively, and are co-latitude and longitude. Cnm and Snm are the (dimensionless) fully-normalized Stokes coefficients and *Pnm* are associated Legendre functions, *n* and *m* are harmonic degree and order respectively. *Re* is the mean Earth's radius (~6,371 km). *N* is the maximum degree of the development, ideally *N* . In practice, the Stokes coefficients are estimated from satellite data with a finite degree N *<* and this maximum value defines the spatial resolution of the geoid ~ *Re/N*. The load Love number coefficients *k'n* account for elastic compensation of Earth's surface in response to mass load variations. A list of values of the first Love number coefficients can be found in Wahr et al. (1998).

By removing harmonic coefficients of a reference ''static'' gravity field from GRACE solutions, time-variations of the Stokes coefficients *Cnm* and *Snm* are computed for each monthly or 10-day period *t*.

For a monthly period *t*, the surface water storage anomaly within a region of angular area (*i.e.*, surface of the studied basin divided by *Re²*) can be obtained as the scalar product between time variations of Stokes coefficients *Cnm*(t) and *Snm*(t), and averaging kernel coefficients *Anm* and *Bnm* corresponding to the geographical region to be extracted (Swenson & Wahr, 2002):

$$
\Delta\sigma\_{\Omega}(\Delta t) = \frac{\rho\_e R\_e}{3\Omega} \sum\_{n=1}^{N} \sum\_{m=1}^{n} \frac{2n+1}{1+k\_n^{'}} \{\Delta C\_{nm}(\Delta t) A\_{nm} + \Delta S\_{nm}(\Delta t) B\_{nm}\} \tag{2}
$$

where *<sup>e</sup>* is the mean Earth's density (~5,517 kg/m3). In case of a perfect kernel and error free data, *Anm* and *Bnm* are the harmonic coefficients of the basin function or mask, which is equal to 1 inside the basin and zero outside.

#### **3.3 The problem of aliasing**

The monthly or sub-monthy GRACE solutions have a sufficient temporal resolution to monitor the long wavelengths variations of the gravity over land. The phenomenon of short term mass variability with period of hours to day, ocean tides and atmosphere are removed using de-aliasing techniques consisting of removing model outputs of these contributions. If the atmospheric pressure fields from ECMWF allow a reasonable de-aliasing of the higher frequencies caused by non-tidal atmospheric mass changes, errors due to tide models appear in the GRACE solutions, especially the aliasing from diurnal (S1) and semi-diurnal (S2) tides (Han et al., 2004; Ray & Lutcke, 2006). The aliasing and the modelling errors are responsible for the spurious meridional undulations present in the GRACE monthly grids and known as north-south striping. They are associated to spherical harmonics order 15 and its multiples (Swenson & Wahr, 2006a; Schrama et al., 2007; Seo et al., 2008). Thompson et al. (2004) showed that the degree error increased by factors ~ 20 due to atmospheric aliasing, ~ 10 due to the ocean model, and ~ 3 due to continental hydrology model. The S2 aliasing errors have also a strong impact on the C20 spherical harmonics coefficient of the geopotential and (Seo et al., 2008; Chen et al., 2009a) and are significant over the Amazon basin but dismish using longer time-series (Schrama et al., 2007; Chen et al., 2009a).

#### **3.4 Accuracy of GRACE measurements: formal, omission and leakage errors**

#### **3.4.1 Formal error**

194 Water Resources Management and Modeling

Geodesy (ITG) at the University of Bonn, Germany, the Delft Institute of Earth Observations and Space Systems (DEOS), Delft, Nederlands, the Goddard Space Flight Center (NASA

Variations in continental water storage cause changes of mass distribution in the Earth system, and thus spatiotemporal changes of the geoid, defined as the gravitational equipotential surface that best coincides with the mean sea surface. In global representation,

> ' 2 ( ) <sup>1</sup> cos( ) ( , ,) (cos ) ( ) (2 1) sin( ) *nm n e nm*

 

 

*q* is the surface distribution of mass, *M* and *S* are the total mass

 

are co-latitude and longitude.

(1)

*Re/N*. The load Love

*Snm* are computed for each

*Snm*(t), and averaging kernel

Cnm and

*qt P S <sup>s</sup> S t n M <sup>m</sup>*

Snm are the (dimensionless) fully-normalized Stokes coefficients and *Pnm* are associated Legendre functions, *n* and *m* are harmonic degree and order respectively. *Re* is the mean Earth's radius (~6,371 km). *N* is the maximum degree of the development, ideally *N* . In practice, the Stokes coefficients are estimated from satellite data with a finite degree N *<*

number coefficients *k'n* account for elastic compensation of Earth's surface in response to mass load variations. A list of values of the first Love number coefficients can be found in

By removing harmonic coefficients of a reference ''static'' gravity field from GRACE

(*i.e.*, surface of the studied basin divided by *Re²*) can be obtained as the scalar product

coefficients *Anm* and *Bnm* corresponding to the geographical region to be extracted (Swenson

*k*

2 1 ( ) () () <sup>3</sup> <sup>1</sup>

free data, *Anm* and *Bnm* are the harmonic coefficients of the basin function or mask, which is

The monthly or sub-monthy GRACE solutions have a sufficient temporal resolution to monitor the long wavelengths variations of the gravity over land. The phenomenon of short term mass variability with period of hours to day, ocean tides and atmosphere are removed using de-aliasing techniques consisting of removing model outputs of these contributions. If

*<sup>R</sup> <sup>n</sup> <sup>t</sup> C tA S tB*

1 1

*n m n*

*N n e e nm nm nm nm*

*<sup>e</sup>* is the mean Earth's density (~5,517 kg/m3). In case of a perfect kernel and error

*Cnm* and

*t*, the surface water storage anomaly within a region of angular area

'

(2)

*Cnm*(t) and

 

 and 

*C t k R m*

and this maximum value defines the spatial resolution of the geoid ~

GSFC), Greenbelt, Maryland, USA.

**3.2 Conversion into coefficients of water mass anomalies** 

it can be expressed as (Heiskanen & Moritz, 1967):

solutions, time-variations of the Stokes coefficients

between time variations of Stokes coefficients

equal to 1 inside the basin and zero outside.

**3.3 The problem of aliasing** 

*t*.

*nm*

and the surface of the Earth, respectively,

for a given time *t*, where

Wahr et al. (1998).

monthly or 10-day period

For a monthly period

& Wahr, 2002):

where 

The degree amplitude of the GRACE error is defined as:

$$\sigma\_N(l) = R\_e \sqrt{\sum\_{m=0}^{l} \left( \sigma\_{C\_{lm}}^2 + \sigma\_{S\_{lm}}^2 \right)} \tag{3}$$

where Clm σ and lm σ*<sup>S</sup>* are the errors on the gravity potential coefficients.

It can be seen as the square-root of the total variance from all terms of a given spatial scale, as the degree *l* is the measure of the spatial scale of a spherical harmonics (*i.e.*, a half wavelength of 20,000/*l* km). These errors increase at degrees 20 to 30 and become dominent at degrees 40 to 50. As a consequence, GRACE monthly solutions are low-pass filtered at degree 50 or 60 to remove the noise contained in the high frequency domain.

#### **3.4.2 Omission or cut-off frequency error**

Error in frequency cut-off represents the loss of energy in the short spatial wavelength due to the low-pass harmonic decomposition of the signals that is stopped at the maximum degree *N1*. For the GRACE land water solutions; *N1*=60, thus the spatial resolution is limited and stopped at ~330 km by construction. This error is simply evaluated by considering the difference of reconstructing the remaining spectrum between two cutting harmonic degrees *N1* and *N2*, where *N2* > *N1* and *N2* should be large enough compared to *N1* (*e.g.*, *N2*=300):

$$
\sigma\_{\text{truncation}} = \sum\_{n=0}^{N\_2} \xi\_n - \sum\_{n=0}^{N\_1} \xi\_n = \sum\_{n=N\_1}^{N\_2} \xi\_n \tag{4}
$$

using the scalar product:

$$\xi\_n = \sum\_{m=0}^{n} \left( \Delta \mathbf{C}\_{nm} \text{(At)} A\_{nm} + \Delta \mathbf{S}\_{nm} \text{(At)} B\_{nm} \right) \tag{5}$$

These errors are generally lower than 1% of the amplitude of the signal.

#### **3.4.3 Leakage error**

Due to the limited space resolution of the GRACE solutions, the estimate of mass variations inside a region of interest such as a drainage basin is affected by the leakage effect in two different ways. First, the mass changes inside a region spread out the whole globe. On the contrary, signals from other regions of the world pollute the estimate inside the considered region of interest. These two opposite effects have to be accounted for to estimate accurate mass variations at the surface of the Earth. Several techniques were proposed to quantify the leakage error. Swenson & Wahr (2002) proposed a geometric formulation of the leakage departure of the shape of the smoothed averaging kernel from that of the exact averaging kernel, Ramillien et al. (2006a) used an « inverse » mask, which is 0 and 1 in and out of the region respectively, developed in spherical harmonics and then truncated at degree 60, Baur et al. (2009) proposed an iterative approach based on forward modelling, Longuevergne et al. (2010) used a method optimizing the basin shape description. The leakage effects are generally estimated using global models as they can not be determined directly from the GRACE observations. Leakage effects are smaller close to the ocean. Over large areas, the signals leaking in and out tend to compensate each other. For the Amazon basin, leakage was found to have a small impact on TWS estimates (Chen et al., 2009b, Xavier et al., 2010). For small basins, such as the High Plain Aquifers (with an area ~ 450,000 km²), an error due to the leakage is of 25 mm in terms of EWH and remains important, compared with annual variations of TWS varying from 100 to 200 mm (Longuevergne et al., 2010).

### **4. Post-processing techniques**

To filter out the spurious north-south stripes present in the GRACE TWS products different techniques based on empirical, statistical, inverse, or regional approaches were applied.

#### **4.1 Empirical methods**

#### **4.1.1 Isotropic and non-isotropic filters**

The simple isotropic Gaussian filter proposed earlier by Jekeli (1981) was commonly used in pre-launched studies (Wahr et al., 1998), and then applied to real data (Wahr et al., 2004; Tapley et al., 2004b) to demonstrate the abilities of GRACE to retrieve time variations of water mass anomalies over land. The choice of the smoothing radius and the inhability of this type of filter to distinguish between geophysical signals and noise are its major drawbacks. More elaborate filters were then developed to reduce the errors on the GRACE signals. Several averaging kernels based on the separate minimization of measurement and leakage error to signal were proposed, by using Lagrange multipliers and a leakage function defined as the difference in shape between exact and smoothed averaging kernels, with no *a priori* (Swenson & Wahr, 2002), or *a priori* (*i.e.*, the shape of the covariance function assuming it is azimuthally symmetric) information and then then minimizing the sum of variance of the satellite and leakage errors (Swenson & Wahr, 2003; Seo & Wilson, 2005).

To take into account the meridian structures of the noise polluting the GRACE signals, several non-isotropic smoothing kernels were applied, using spectral Legendre coefficients dependent on both degree and order (Han et al., 2005), or Tikhonov-type regularization (Kusche, 2007).

#### **4.1.2 Destriping method**

196 Water Resources Management and Modeling

Due to the limited space resolution of the GRACE solutions, the estimate of mass variations inside a region of interest such as a drainage basin is affected by the leakage effect in two different ways. First, the mass changes inside a region spread out the whole globe. On the contrary, signals from other regions of the world pollute the estimate inside the considered region of interest. These two opposite effects have to be accounted for to estimate accurate mass variations at the surface of the Earth. Several techniques were proposed to quantify the leakage error. Swenson & Wahr (2002) proposed a geometric formulation of the leakage departure of the shape of the smoothed averaging kernel from that of the exact averaging kernel, Ramillien et al. (2006a) used an « inverse » mask, which is 0 and 1 in and out of the region respectively, developed in spherical harmonics and then truncated at degree 60, Baur et al. (2009) proposed an iterative approach based on forward modelling, Longuevergne et al. (2010) used a method optimizing the basin shape description. The leakage effects are generally estimated using global models as they can not be determined directly from the GRACE observations. Leakage effects are smaller close to the ocean. Over large areas, the signals leaking in and out tend to compensate each other. For the Amazon basin, leakage was found to have a small impact on TWS estimates (Chen et al., 2009b, Xavier et al., 2010). For small basins, such as the High Plain Aquifers (with an area ~ 450,000 km²), an error due to the leakage is of 25 mm in terms of EWH and remains important, compared with annual

variations of TWS varying from 100 to 200 mm (Longuevergne et al., 2010).

the satellite and leakage errors (Swenson & Wahr, 2003; Seo & Wilson, 2005).

To filter out the spurious north-south stripes present in the GRACE TWS products different techniques based on empirical, statistical, inverse, or regional approaches were applied.

The simple isotropic Gaussian filter proposed earlier by Jekeli (1981) was commonly used in pre-launched studies (Wahr et al., 1998), and then applied to real data (Wahr et al., 2004; Tapley et al., 2004b) to demonstrate the abilities of GRACE to retrieve time variations of water mass anomalies over land. The choice of the smoothing radius and the inhability of this type of filter to distinguish between geophysical signals and noise are its major drawbacks. More elaborate filters were then developed to reduce the errors on the GRACE signals. Several averaging kernels based on the separate minimization of measurement and leakage error to signal were proposed, by using Lagrange multipliers and a leakage function defined as the difference in shape between exact and smoothed averaging kernels, with no *a priori* (Swenson & Wahr, 2002), or *a priori* (*i.e.*, the shape of the covariance function assuming it is azimuthally symmetric) information and then then minimizing the sum of variance of

To take into account the meridian structures of the noise polluting the GRACE signals, several non-isotropic smoothing kernels were applied, using spectral Legendre coefficients dependent on both degree and order (Han et al., 2005), or Tikhonov-type regularization

**4. Post-processing techniques** 

**4.1.1 Isotropic and non-isotropic filters** 

**4.1 Empirical methods** 

(Kusche, 2007).

**3.4.3 Leakage error** 

This filtering technique, known as destriping or correlated-error filter, was designed to remove the correlated errors in the Stokes coefficients responsible for the north-south stripes in the GRACE solutions. For a particular order *m*, the Stokes coefficients are smoothed with a quadratic polynomial in a moving window of width *w* centered about degree *l*. The sum only concerns the terms of the same parity as *l* (Swenson and Wahr, 2006a). This method provides better results at high latitudes than close to the equator as residual errors due to the near-sectorial coefficients (*l*~*m*) are not completely removed (Swenson and Wahr, 2006a; Klees et al., 2008a). Besides, some short-wavelength features are removed, mainly at high latitudes. Monthly CSR, GFZ and JPL destriped and smoothed EWH grids of 1-degree spatial resolution are available over 2002-2010 at: http://grace.jpl.nasa.gov.

#### **4.1.3 Stabilization criteria**

The Level-2 GRGS-EIGEN-GL04 models are derived from Level-1 GRACE measurements including KBRR, and from LAGEOS-1/2 SLR data for enhancement of lower harmonic degrees (Lemoine et al., 2007a; Bruinsma et al., 2010). These gravity fields are expressed in terms of normalized spherical harmonic coefficients from degree 2 up to degree 50-60 using an empirical stabilization approach without any smoothing or filtering. This stabilization approach consists in adding empirically determined degree and order dependent coefficients to minimize the time variations of the signal measured by GRACE over ocean and desert without significantly affecting the amplitude of the signal over large drainage basins. Monthly (Release-1) and 10-day (Release-2) TWS grids of 1-degree spatial resolution are available over 2002-2010 at: http://grgs.obs-mip.fr.

#### **4.2 Statistical methods**

#### **4.2.1 Wiener filtering**

The spatial averaging of monthly GRACE using the Wiener filter data is based on the leastsquare minimization of the difference between modelled and filtered signals. This filter is isotropic, depending only on the degree power of the signal and noise models. Using hydrology and ocean models outputs for simulating the degree power spectrum of the GRACE signal, the formal error of the GRACE coefficients, and GRACE solutions, Sasgen et al. (2006) were able to determine the optimal coefficients of the Wiener filter, and to demonstrate that this method is robust approach for low-pass filtering the GRACE data which does not need the specification of any averaging radius.

#### **4.2.2 Principal Components Analysis**

Principal Components Analysises (PCA) were performed on time-series of Gaussian-filtered EWH grids (Schrama et al., 2007) and Stokes coefficients (Wouters & Schrama, 2007). PCA was able to efficiently separate modes of geophysical signals from modes corresponding to the North-South striping in the EWH grids. The three first modes explain 73.5% of the continental hydrology with a Gaussian prefiltering of 6.25° of radius. The residual modes contain S2 aliasing errors and a semiannual continental hydrology signal contained in the Global Land Data Assimilation Systems (GLDAS) model (Schrama et al., 2007). Directly applied to the Stokes coefficients, PCA also successfully removes the meridian undulations present in the GRACE solutions. The results of this technique are expected to be better as the time-series of the GRACE solutions become longer (Wouters & Schrama, 2007).

#### **4.2.3 Independent Components Analysis**

The Independent Component Analysis (ICA) approach is used to extract hydrological signals from the noise in the Level-2 GRACE solutions by considering completely objective constraints, so that the gravity component of the observed signals is forced to be uncorrelated numerically. This blind separation method is based on the assumption of statistical independence of the elementary signals that compose the observations, *i.e.*, geophysical signals and spurious noise, and does not require other *a priori* information (Frappart et al., 2010a). Comparisons at a global scale showed that the ICA-based solutions present less north–south stripes than Gaussian and destriped solutions on the land, and more realistic hydrological structures than the destriped solutions in the tropics. ICA filtering seems to allow the separation of the GIA from the TWS as negative trends were found over the Laurentides and Scandinavia. Unfortunately, this important geophysical parameter does not appear clearly in an independent mode yet. At basin-scale, the ICAbased solutions allowed us to filter out the unrealistic peaks present in the time series of TWS obtained using classical filtering for basins with areas lower than 1 million km² (Frappart et al., 2011a). The major drawback of this approach is that it cannot directly be applied to the GRACE Level-2 raw data, as a first step of (Gaussian) prefiltering is required. Monthly CSR, GFZ and JPL ICA-filtered TWS grids of 1-degree spatial resolution are available over 2002-2010 at: http://grgs.obs-mip.fr.

#### **4.2.4 Kalman smoothing**

Daily water masses solutions have been derived from GRACE observations using a Kalman filter approach. They are estimated using a Gauss Markov model, the solution at day *t+1* slightly differs from solution at day *t* from the noise prediction (first-order Markov process) estimated using a Kalman filter. A *priori* information on the hydrological patterns from the WGHM model was introduced in the covariance matrix to compensate the small number of available observations over a daily time-span. This approach enables to produce temporarily enhanced solutions without loss in spatial resolution and resolution (Kurtenbach et al., 2009). Daily TWS Stokes coefficients are available over 2002-2010 at: http://www.igg.unibonn.de/apmg.

#### **4.3 Inverse methods**

#### **4.3.1 Iterative least-squares approach**

This method is based on the matrix formalism of the generalized least-squares criteria (Tarantola, 1987) and consists in estimating separately the spherical harmonics coefficients, in terms of EWH, of four different fluid reservoirs (atmosphere, oceans, soil water, and snow) from the monthly GRACE geoids. For each water reservoir, the anomaly of mass is estimated after the convergence of an inverse approach combining the GRACE observations with stochastic properties of the unknown hydrological signal (Ramillien et al., 2004; 2005).

#### **4.3.2 Optimal filter**

198 Water Resources Management and Modeling

present in the GRACE solutions. The results of this technique are expected to be better as the

The Independent Component Analysis (ICA) approach is used to extract hydrological signals from the noise in the Level-2 GRACE solutions by considering completely objective constraints, so that the gravity component of the observed signals is forced to be uncorrelated numerically. This blind separation method is based on the assumption of statistical independence of the elementary signals that compose the observations, *i.e.*, geophysical signals and spurious noise, and does not require other *a priori* information (Frappart et al., 2010a). Comparisons at a global scale showed that the ICA-based solutions present less north–south stripes than Gaussian and destriped solutions on the land, and more realistic hydrological structures than the destriped solutions in the tropics. ICA filtering seems to allow the separation of the GIA from the TWS as negative trends were found over the Laurentides and Scandinavia. Unfortunately, this important geophysical parameter does not appear clearly in an independent mode yet. At basin-scale, the ICAbased solutions allowed us to filter out the unrealistic peaks present in the time series of TWS obtained using classical filtering for basins with areas lower than 1 million km² (Frappart et al., 2011a). The major drawback of this approach is that it cannot directly be applied to the GRACE Level-2 raw data, as a first step of (Gaussian) prefiltering is required. Monthly CSR, GFZ and JPL ICA-filtered TWS grids of 1-degree spatial resolution are

Daily water masses solutions have been derived from GRACE observations using a Kalman filter approach. They are estimated using a Gauss Markov model, the solution at day *t+1* slightly differs from solution at day *t* from the noise prediction (first-order Markov process) estimated using a Kalman filter. A *priori* information on the hydrological patterns from the WGHM model was introduced in the covariance matrix to compensate the small number of available observations over a daily time-span. This approach enables to produce temporarily enhanced solutions without loss in spatial resolution and resolution (Kurtenbach et al., 2009). Daily TWS Stokes coefficients are available over 2002-2010 at: http://www.igg.uni-

This method is based on the matrix formalism of the generalized least-squares criteria (Tarantola, 1987) and consists in estimating separately the spherical harmonics coefficients, in terms of EWH, of four different fluid reservoirs (atmosphere, oceans, soil water, and snow) from the monthly GRACE geoids. For each water reservoir, the anomaly of mass is estimated after the convergence of an inverse approach combining the GRACE observations with stochastic properties of the unknown hydrological signal (Ramillien et

time-series of the GRACE solutions become longer (Wouters & Schrama, 2007).

**4.2.3 Independent Components Analysis** 

available over 2002-2010 at: http://grgs.obs-mip.fr.

**4.2.4 Kalman smoothing** 

bonn.de/apmg.

al., 2004; 2005).

**4.3 Inverse methods** 

**4.3.1 Iterative least-squares approach** 

This method takes into account the statistical information present in each GRACE monthly solution. The noise and full signal variance-covariance matrix is used to tailor the filter to the error characteristics of a particular monthly solution. The resulting filter was found to be both anisotropic and non-symmetric to accommodate noise of an arbitrary shape as the north-south stripes present in the GRACE monthly solutions. It is optimal as it minimizes the difference between the signal and the filtered GRACE estimate in the least-squares sense. This filter was found to perform better than any isotropic, or anisotropic and symmetric filters preserving better the signal amplitudes better (Klees et al., 2008a). Monthly GRACE-derived TWS grids of 1-degree spatial resolution computed using the optimal filter are available over 2002-2010 at: http://lr.tudelft.nl/en/organisation/departments-andchairs/remote-sensing/physical-and-space-geodesy/data-and-models/dmt-1/

#### **4.4 Regional approaches**

Since the begining of the low-low satellite gravimetry missions dedicated to Earth's gravity field measurement (*i.e.*, CHAMP, GRACE, GOCE), determining locally the time variations of the surface water storage has been of great interest for region where very localized strong mass variations occur, such as flood and glaciers fields. Moreover, it represents an alternative to the problem of numerical singularity at poles by choosing a suitable geometry of surface tiles. So far, rectangular surface elements have been considered but other geometries have successfully tested (Eicker, 2008; Ramillien et al., 2011).

#### **4.4.1 Mascons**

In this local approach, the mass of water in surface 4-by-4 degree blocks has been explicitly solved using the GRACE inter-satellite KBR Rate (KBRR) data for continental hydrology and collected over the region of interest. This size of the surface elements has been chosen as the limit of *a priori* spatial resolution of the GRACE data (*i.e.*, 400 km). This method uses inherently empirical spatial and temporal constraints among the coefficients that are dependent on geographical locations. The local representation of gravity minimizes the leakage error from other areas due to aliasing or mis-modelling (Rowlands et al., 2005; Lemoine et al., 2007b). It has revealed more information about interannual variability in the subregions than any Stokes coefficients-based approach. Global and mascons methods based on spherical harmonics have provided similar results for the mass trend in the drainage basins, especially when no spatial constraints are taken into account (Rowlands et al., 2010). But, inherent problem of spectral truncation such as leakage are not avoided while using spherical harmonics. 10-day and 4° Mascons solutions and corresponding GLDAS model time series for each continental surface element are available at: http://grace.sgtinc.com/V2/Global.html.

#### **4.4.2 Regional method**

Another type of regional approach has been recently proposed by adjusting the surface mass density from the Level-1 GRACE data, in particular the accurate satellite-to-satellite velocity variations (Ramillien et al., 2011). Once these observations are corrected from known accelerations (atmosphere and ocean mass changes, tides, static gravity field…), along track residuals of KBR-rate are used to compute kinetic energy (or equivently potential anomaly) variations between the twin GRACE satellites, corresponding mainly to continental water redistributions. The linear system of equations to solve between observations (*i.e.*, potential anomalies) and parameters (*i.e.*, surface elements) is built according to the first Newton's law of attraction of masses. Unlike the mascons solutions, spherical harmonics are not used as a basis of regional orthogonal functions, consequently this regional method does not suffer from the drawbacks related to any spectral truncation. Singular Value Decomposition (SVD) and L-curve analysis are considered to compute regularized 10-day 2-by-2 degree solutions of the ill-posed of gravimetry. Lately, another version of the inversion scheme takes spatial correlations versus the geographical distance between surface elements into account (Ramillien et al., in press). Time series of regional solutions over South America have been produced and compared to global solutions and mascons, they reveal very comparable amplitudes of seasonal cycle over the Amazon basin, as well as sudden flooding events.

### **5. Applications to continental hydrology: results and discussion**

#### **5.1 Application in continental hydrology and validation**

Once sufficiently long series of GRACE solutions were made available by the former official providers (CSR, JPL, GFZ), various regional studies for validating these products were made. Examples of GRACE-based TWS maps are given in Figure 1 showing the ability of the gravimetry from space mission to restitue both seasonal and interannual variability.

Comparisons with global hydrology model outputs and surface measurements revealed acceptable agreements between GRACE-derived and model changes of continental water storage versus time, especially at seasonal time scale. GRACE provides directly the term of TWS variation of the net balance equation. In general, GRACE-based estimates of TWS compare favourably with those based on land surface models as well as atmospheric and terrestrial water balances (Rodell et al., 2004a; Ramillien et al., 2005; Syed et al., 2005; Klees et al., 2008b). The structures seen on global maps of seasonal amplitudes are comparable to those described by the WaterGap model developed by Döll et al. (2003) and GLDAS (Rodell et al., 2004a). Over the entire Northern Hemisphere, GLDAS water storage simulations with a resolution of 1,300 km and accuracy of 9 mm in terms of EWH were found to have a spatial correlation of 0.65 with GRACE data, suggesting that gravity field changes may be related to TWS variations (Andersen & Hinderer, 2005). However, TWS variations tend to be slightly over-estimated by GRACE (Schmidt et al., 2006; Syed et al., 2008). GRACE-based TWS changes were also validated by accurate observations of superconducting gravimeters during the 2003 heat wave that occured in Central Europe (Andersen et al., 2005). They were used to detect the exceptional drought (Chen et al., 2009b) and flooding (Chen et al., 2010) affected the Amazon basin in 2005 and 2009 respectively (Figure 2a). Correlations of 0.7-0.8 between GRACEbased TWS and *in situ* measurements of water level along the Amazon River were found by Xavier et al. (2010). Combined with other satellite techniques such as imagery, GRACE geoid data were used to study the mechanisms of seasonal flooding in large inundation areas like the Mekong delta (Frappart et al., 2006a) and large river basins as the Amazon River (Frappart et al., 2008; Papa et al., 2008).

velocity variations (Ramillien et al., 2011). Once these observations are corrected from known accelerations (atmosphere and ocean mass changes, tides, static gravity field…), along track residuals of KBR-rate are used to compute kinetic energy (or equivently potential anomaly) variations between the twin GRACE satellites, corresponding mainly to continental water redistributions. The linear system of equations to solve between observations (*i.e.*, potential anomalies) and parameters (*i.e.*, surface elements) is built according to the first Newton's law of attraction of masses. Unlike the mascons solutions, spherical harmonics are not used as a basis of regional orthogonal functions, consequently this regional method does not suffer from the drawbacks related to any spectral truncation. Singular Value Decomposition (SVD) and L-curve analysis are considered to compute regularized 10-day 2-by-2 degree solutions of the ill-posed of gravimetry. Lately, another version of the inversion scheme takes spatial correlations versus the geographical distance between surface elements into account (Ramillien et al., in press). Time series of regional solutions over South America have been produced and compared to global solutions and mascons, they reveal very comparable amplitudes of seasonal cycle over the Amazon basin,

**5. Applications to continental hydrology: results and discussion** 

Once sufficiently long series of GRACE solutions were made available by the former official providers (CSR, JPL, GFZ), various regional studies for validating these products were made. Examples of GRACE-based TWS maps are given in Figure 1 showing the ability of the gravimetry from space mission to restitue both seasonal and interannual variability.

Comparisons with global hydrology model outputs and surface measurements revealed acceptable agreements between GRACE-derived and model changes of continental water storage versus time, especially at seasonal time scale. GRACE provides directly the term of TWS variation of the net balance equation. In general, GRACE-based estimates of TWS compare favourably with those based on land surface models as well as atmospheric and terrestrial water balances (Rodell et al., 2004a; Ramillien et al., 2005; Syed et al., 2005; Klees et al., 2008b). The structures seen on global maps of seasonal amplitudes are comparable to those described by the WaterGap model developed by Döll et al. (2003) and GLDAS (Rodell et al., 2004a). Over the entire Northern Hemisphere, GLDAS water storage simulations with a resolution of 1,300 km and accuracy of 9 mm in terms of EWH were found to have a spatial correlation of 0.65 with GRACE data, suggesting that gravity field changes may be related to TWS variations (Andersen & Hinderer, 2005). However, TWS variations tend to be slightly over-estimated by GRACE (Schmidt et al., 2006; Syed et al., 2008). GRACE-based TWS changes were also validated by accurate observations of superconducting gravimeters during the 2003 heat wave that occured in Central Europe (Andersen et al., 2005). They were used to detect the exceptional drought (Chen et al., 2009b) and flooding (Chen et al., 2010) affected the Amazon basin in 2005 and 2009 respectively (Figure 2a). Correlations of 0.7-0.8 between GRACEbased TWS and *in situ* measurements of water level along the Amazon River were found by Xavier et al. (2010). Combined with other satellite techniques such as imagery, GRACE geoid data were used to study the mechanisms of seasonal flooding in large inundation areas like the Mekong delta (Frappart et al., 2006a) and large river basins as the Amazon

**5.1 Application in continental hydrology and validation** 

River (Frappart et al., 2008; Papa et al., 2008).

as well as sudden flooding events.

Fig. 1. Examples of GRACE-derived TWS maps processed using a Gaussian filtering of radius 400 km and ICA (Level-2 Release 4 GFZ solutions) for a) March 2005, b) September 2005, c) March 2006, d) September 2006, e) March 2007, f) September 2007, g) March 2008, h) September 2008, i) March 2009, j) September 2009.

#### **5.2 Isolating groundwater changes**

Variations in the groundwater storage can be extracted from *TWS* measured by GRACE using external information on the other hydrological reservoirs. The time variations in *TWS* are the sum of the contributions of the different reservoirs present in a drainage basin:

$$
\Delta TV \text{US} = \Delta SV + \Delta SN + \Delta TSS \tag{6}
$$

with:

$$
\Delta TSS = \Delta RZ + \Delta GW + \Delta P \tag{7}
$$

where *SW* represents the total surface water storage including lakes, reservoirs, in-channel and floodplains water, *SN* is the snow storage, *TSS* is the total soil storage including *RZ* the water contained in the root zone of the soil (generally representing a depth of 1 or 2 m), *GW* the groundwater storage in the aquifers, and *P* the permafrost storage at boreal latitudes.

Post-launch studies of GRACE-based groundwater remote sensing have clearly demonstrated that when combined with ancillary measurements of surface waters and soil moisture, either modeled or observed, GRACE is capable of monitoring changes in groundwater storage with reasonable accuracy. Important seasonal correlations of 0.8-0.9

Fig. 2. Time series of TWS anomalies (mm of EWH) from Level-2 Release 4 GFZ GRACE solutions (Gaussian filtering of radius 400 km and ICA) for a) Amazon, b) Ganges.

were found by confronting GRACE to recordings of well network in Illinois (Yeh et al., 2006), Oklahoma (Swenson et al., 2008), the High Plain aquifer (Strassberg et al., 2007) and the Mississippi basin (Rodell et al. 2007). This method was successfully applied to monitor large-scale climatic and anthropogenic impacts on water resources, as the detection of the recent severe drought in the Murray-Darling basin in Southern Australia (Leblanc et al., 2009) or the depletion of the aquifers in India (Rodell et al., 2009) (Figure 2b), but also to estimate aquifer storage parameters as the specific yield or storativity (Sun et al., 2010). For large drainage basins covered with extensive floodplains, changes in water stored in the aquifer is isolated from the TWS measured by GRACE by removing contributions of both the surface reservoir, derived from satellite imagery and radar altimetry, and the root zone reservoir simulated by hydrological models (Frappart et al., 2010b; 2011b). In the Negro River basin, the groundwater anomalies show a realistic spatial pattern compared with the hydrogeological map of the basin (DNPM, 1983), and similar temporal variations to local in situ groundwater observations (Figure 3).

#### **5.3 Estimates of snow mass changes at high latitudes**

After a separation of water and snow contributions from observed gravity field using an inverse iterative approach based on the least-squares criteria, Frappart et al. (2006b) used satellite microwave data and global hydrological model outputs to validate GRACE-derived snow mass changes at high-latitudes. In particular, rms errors of 10-20 mm were found for Yenisey, Ob, McKenzie and Yukon basins. Comparison at basin scales between GRACEderived snow mass variations, GRACE-based TWS variations and river discharge temporal evolution shows that the snow component has a more significant impact on river discharge at high latitudes than TWS, and that snow and streamflow present similar interannual variabilities (Frappart et al., 2011c).

#### **5.4 Estimates of regional geo-hydrological parameters**

#### **5.4.1 River discharges**

202 Water Resources Management and Modeling

Variations in the groundwater storage can be extracted from *TWS* measured by GRACE using external information on the other hydrological reservoirs. The time variations in *TWS* are the sum of the contributions of the different reservoirs present in a drainage basin:

where *SW* represents the total surface water storage including lakes, reservoirs, in-channel and floodplains water, *SN* is the snow storage, *TSS* is the total soil storage including *RZ* the water contained in the root zone of the soil (generally representing a depth of 1 or 2 m), *GW* the groundwater storage in the aquifers, and *P* the permafrost storage at boreal latitudes.

Post-launch studies of GRACE-based groundwater remote sensing have clearly demonstrated that when combined with ancillary measurements of surface waters and soil moisture, either modeled or observed, GRACE is capable of monitoring changes in groundwater storage with reasonable accuracy. Important seasonal correlations of 0.8-0.9

a) Amazon

2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 **Time**

b) Ganges

Fig. 2. Time series of TWS anomalies (mm of EWH) from Level-2 Release 4 GFZ GRACE solutions (Gaussian filtering of radius 400 km and ICA) for a) Amazon, b) Ganges.

2002 2004 2006 2008 2010 2012 **Time**

*TWS SW SN TSS* (6)

*TSS RZ GW P* (7)

Trend : -16 mm.yr-1

**5.2 Isolating groundwater changes** 



**TWS anomaly (mm of EWH)**

**TWS anomaly (mm of EWH)**

with:

GRACE-based estimates of river discharges are based on the resolution of the water balance equation at basin-scale:

$$\frac{\partial \mathbf{\widetilde{S}}}{\partial \mathbf{t}} = P \text{ - } \mathbf{E} \text{ - } \mathbf{R} \tag{8}$$

where *S* is land water storage (TWS), *P* and *ET* the basin-wide totals of precipitation and evapotranspiration, and *R* the total basin discharge, or the net surface and groundwater outflow.

As *ET* presents large uncertainties, the quantity *P-ET* is eliminated using the water balance of the atmospheric branch of the water cycle:

$$\frac{\partial \vec{\text{MV}}}{\partial t} = ET \cdot P \text{ - } div \vec{\text{Q}} \tag{9}$$

where *W* is the vertically-integrated precipitable water, *divQ* the divergence of the verticallyintegrated average atmospheric moisture flux vector, computed using the following equations:


Fig. 3. When combined with other data, GRACE-derived TWS anomalies provide valuable information on changes affecting the aquifers. Above a comparison between a) Map of annual amplitude of GW in the Negro River Basin. b) Hydrogeological map of Brazil from DNPM (1983), available at: http://eusoils.jrc.ec.europa.eu/esdb\_archive/EuDASM/EUDASM.htm, in the framework of the European Commission — Joint Research Centre through the European Digital Archive of Soil Maps (EuDASM) (Panagos et al., in press). Adapted from Frappart et al. (2011b).

$$\mathbf{W} = \int\_{\mathbf{P}\_{\perp}}^{\mathbf{P}\_{\ast}} \frac{\mathbf{q}}{\mathbf{g}} \, \mathrm{d}\mathbf{p} \tag{10}$$

$$
\vec{\mathbf{Q}} = \int\_{\mathbf{p}\_r}^{\mathbf{p}\_s} \frac{\mathbf{q}}{\mathbf{g}} \vec{V} \mathbf{d} \mathbf{p} \tag{11}
$$

where *ps* and *pT* are the pressure at the surface and the top of the atmosphere respectively, *q* is the specific humidity, *g* is gravitational acceleration and *V* is the horizontal wind vector.

This approach was firstly applied to the Amazon and the Mississippi River basins (Syed et al., 2005), and extended to obtain a global volume of freshwater discharge estimated here is 30,354 ± 1,212 km3/y (Syed et al., 2009).

#### **5.4.2 Vertical water fluxes**

204 Water Resources Management and Modeling

Fig. 3. When combined with other data, GRACE-derived TWS anomalies provide valuable information on changes affecting the aquifers. Above a comparison between a) Map of annual amplitude of GW in the Negro River Basin. b) Hydrogeological map of Brazil from DNPM (1983), available at: http://eusoils.jrc.ec.europa.eu/esdb\_archive/EuDASM/EUDASM.htm, in the framework of the European Commission — Joint Research Centre through the European Digital Archive of Soil Maps (EuDASM) (Panagos et al., in press). Adapted from Frappart et al.

Local aquifers restricted to fractured zones. Small hydrogeological reservoir.

Local aquifers or continuous aquifers of limited extension. Two levels of water:

Continuous aquifers of regional extension, free or confined. Medium

Almost no aquifer. Very small hydrogeological reservoir.

free and/or confined. Small hydrogeological reservoir.

hydrogeological reservoir.

(2011b).

GRACE-based TWS can be also used to estimate changes in vertical water fluxes solving the water balance equation (8). Changes of regional evapotranspiration (ET) rate over Mississippi basin were estimated by combining 600-km low-pass filtered GRACE TWS data with observed precipitation and streamflow in the water balance equation (Rodell et al., 2004b). Similarly, Ramillien et al. (2006b) solved the water balance equation for ET rate. Time variations of ET were evaluated over 16 large drainage basins solving the water balance equation and using precipitation from GPCC and runoff from the WGHM model (Döll et al., 2003), and they revealed that GRACE-based ET variations were comparable to the ET value simulated by four differend global hydrology models. Swenson and Wahr (2006b) estimated the difference precipitation minus evapotranspiration using the water balance framework and comparing to surface parameters variations from global analysis.

#### **5.5 Contribution of TWS to sea level variations**

Seasonal and interannual changes in land water storage derived from GRACE were used to estimate the contribution of land hydrology to the sea level variations for 27 of the largest drainage basins in the world. Estimation of 2002-2006 sea level contribution of GRACEderived TWS of the largest basins was made using the GRGS GRACE solutions, and corresponds to a water loss of ~ 0.5 mm per year ESL (Ramillien et al., 2008). Wouters et al. (2011) recently estimated the global mean eustatic cycle, i.e., the total amount of water exchanged between continents and ocean, to be 9.4 ± 0.6 mm equivalent water level over the period 2003-2010. In the Pan Arctic regions, snow volume derived from the Special Sensor Microwave/Imager (SSM/I) radiometric measurements over 1989-2006 and GRACE over 2002-2007 is a key driver of the sea level seasonal cycle, but snow volume trend indicates a negligible and not statistically significant contribution to sea level variations (Biancamaria et al., 2011).

#### **5.6 Optimization and assimilation into hydrological models**

The GRACE data have been widely confronted to global hydrological ouputs for comparison and validation, and then to improve the description and parameterization of the hydrological processes in the models. The GRACE geoid data were used as a proxy to test and improve the efficiency of surface waters schemes in the ORCHIDEE land surface hydrology model (Ngo-Duc et al., 2007). Unfortunately, a lack of consistency in the seasonal cycle and a time shift between GRACE TWS and global hydrology models may persist. The comparison between six GRACE-derived TWS and 9 land surface models forced with the same forcings over West Africa permits to identify the processes needing improvement in the land surface models, especially the correct simulation of slow water reservoirs as well as evapotranspiration during the dry season for accurate soil moisture modeling over West Africa (Grippa et al., 2011).

Güntner (2008) pointed out the importance of integrating the GRACE data into hydrological models for improving the reliability of their prediction through advanced methods of multiobjective calibration and data assimilation. A multicalibration approach to constrain model predictions by both measured river discharge and TWS anomalies from GRACE was apllied to the WGHM model, improving simulation results with regard to both objectives. Using only monthly TWS variations, the RMSE was reduced of about 25 mm for the Amazon, 6 mm for the Mississippi and 1 mm for the Congo river basins (Werth et al., 2009). GRACEderived monthly TWS anomalies were also assimilated into one of the GLDAS models using a Kalman smoother approach. Compared with open-loop simulations, assimilated ones exhibited better performance thanks to improvements in the surface and groundwater estimates (Zaitchik et al., 2008).

#### **5.7 Detecting long-term variations over ice sheets**

GRACE-derived mass balance estimates of Antarctica (Velicogna & Wahr, 2006a), Greenland (Velicogna & Wahr, 2006b; Chen et al., 2006) or both ice sheets (Ramillien et al., 2006a) have been determined from the Level-2 solutions of the GRACE project. Velicogna and Wahr (2006b) found a lost rate of -82 ± 28 Gt/y for Greenland with an acceleration of melting in Spring 2004. Mass changes of this ice sheet resolved by drainage basin were estimated using the NASA mascons solutions as -101 ± 16 Gt/y after correction for postglacial rebound (Luthcke et al., 2006). This latter result is consistent with the one found previously by Ramillien et al. (2006a) using 10-day GRGS GRACE solutions. Velicogna & Wahr (2006a) estimated an extreme decrease of ice mass of 152 ± 80 Gt/y for Antarctica. The problem is that this main value represents the Post-Glacial Rebound (PGR) correction itself that the authors removed from the GRACE data using the IJ05 (Ivins and James, 2005) and ICE-5G (Peltier, 2004). If GRACE data are not corrected from PGR, the mass balance of Antarctica appears close to zero, with large uncertainties due to North-South striping and leakage effects. This result shows it is not clear whether this continent actually looses or gains ice mass during the recent period (Ramillien et al., 2006a). PGR phenomena are not well-modelled, especially over the whole Antarctica where available long-term observational constraints remain rare. Consequently, removing PGR that still cannot be modelled accurately represent a source of important errors. Velicogna & Wahr (2002) have shown earlier that the detection of PGR remains possible by combination of different satellite techniques, using 5 years of simulated GRACE and Laser GLAS data over the whole Antarctica. However, an effective extraction of PGR signals using GRACE products still needs to be demonstrated.

As mentioned by Baur et al. (2009), even if the PGR correction is less than 10 Gt/y over Greenland, mass balance for this continental area using GRACE solutions can vary strongly from simple to twice, according to which provider (CSR, JPL, GFZ) is considered, and which processing is applied. Velicogna and Wahr (2005) found a constant decrease of -84 ±-28 Gt/y for the period 2002-2004. Later, Wouters et al. (2008) adjusted a rate of -171 Gt/y for 2003- 2008, with an acceleration during summer. Interestingly, Ramillien et al. (2006a) and Luthcke et al. (2006) found comparable loss rates of about 110 Gt/y considering 10-day GRGS and mascons solutions, respectively, for 2003-2005. Recent re-evaluation of Greenland ice mass loss provides -65 Gt/y using ICA-400 solutions and around 100 Gt/y using classical Gaussian filtered solution for 2003-2010 (Bergmann et al., 2012). It is clear that an interannual variations of Greenland ice sheet exists, as an acceleration of the ice melting in Spring 2004 and 2007-2009 (Velicogna & Wahr, 2006b; Velicogna, 2009). A recent study pointed a clear de-acceleration of the ice melting over Greenland for the very last years (Bergmann et al., 2012), which was confirmed by a recent analysis of Rignot et al. (2011).

Level-2 GRACE solutions have also been used to quantify ice mass loss of coastal glaciers in Western Antarctica, Patagonia (Chen et al., 2007) and in the Gulf of Alaska (Luthcke et al., 2008) as the response of the recent global warming. Chen et al. (2007) found a depletion of mass of -27.9 ± 11 Gt/y for the Patagonian glaciers between April 2002 and December 2007. For glaciers in Alaska, negative trends of -84 Gt/y during 2003-2007 and -102 Gt/y for 2003- 2006 were deduced from mascon solutions. Once again, trend estimates depend strongly upon either using a model of PGR, or estimating it together with the current mass loss. Another complication is the coarse footprint of GRACE which mixes signals from neighbouring basins (Horwath and Dietrich, 2009).

#### **6. Conclusion**

206 Water Resources Management and Modeling

hydrology model (Ngo-Duc et al., 2007). Unfortunately, a lack of consistency in the seasonal cycle and a time shift between GRACE TWS and global hydrology models may persist. The comparison between six GRACE-derived TWS and 9 land surface models forced with the same forcings over West Africa permits to identify the processes needing improvement in the land surface models, especially the correct simulation of slow water reservoirs as well as evapotranspiration during the dry season for accurate soil moisture modeling over West

Güntner (2008) pointed out the importance of integrating the GRACE data into hydrological models for improving the reliability of their prediction through advanced methods of multiobjective calibration and data assimilation. A multicalibration approach to constrain model predictions by both measured river discharge and TWS anomalies from GRACE was apllied to the WGHM model, improving simulation results with regard to both objectives. Using only monthly TWS variations, the RMSE was reduced of about 25 mm for the Amazon, 6 mm for the Mississippi and 1 mm for the Congo river basins (Werth et al., 2009). GRACEderived monthly TWS anomalies were also assimilated into one of the GLDAS models using a Kalman smoother approach. Compared with open-loop simulations, assimilated ones exhibited better performance thanks to improvements in the surface and groundwater

GRACE-derived mass balance estimates of Antarctica (Velicogna & Wahr, 2006a), Greenland (Velicogna & Wahr, 2006b; Chen et al., 2006) or both ice sheets (Ramillien et al., 2006a) have been determined from the Level-2 solutions of the GRACE project. Velicogna and Wahr (2006b) found a lost rate of -82 ± 28 Gt/y for Greenland with an acceleration of melting in Spring 2004. Mass changes of this ice sheet resolved by drainage basin were estimated using the NASA mascons solutions as -101 ± 16 Gt/y after correction for postglacial rebound (Luthcke et al., 2006). This latter result is consistent with the one found previously by Ramillien et al. (2006a) using 10-day GRGS GRACE solutions. Velicogna & Wahr (2006a) estimated an extreme decrease of ice mass of 152 ± 80 Gt/y for Antarctica. The problem is that this main value represents the Post-Glacial Rebound (PGR) correction itself that the authors removed from the GRACE data using the IJ05 (Ivins and James, 2005) and ICE-5G (Peltier, 2004). If GRACE data are not corrected from PGR, the mass balance of Antarctica appears close to zero, with large uncertainties due to North-South striping and leakage effects. This result shows it is not clear whether this continent actually looses or gains ice mass during the recent period (Ramillien et al., 2006a). PGR phenomena are not well-modelled, especially over the whole Antarctica where available long-term observational constraints remain rare. Consequently, removing PGR that still cannot be modelled accurately represent a source of important errors. Velicogna & Wahr (2002) have shown earlier that the detection of PGR remains possible by combination of different satellite techniques, using 5 years of simulated GRACE and Laser GLAS data over the whole Antarctica. However, an effective extraction of PGR signals using GRACE products still

As mentioned by Baur et al. (2009), even if the PGR correction is less than 10 Gt/y over Greenland, mass balance for this continental area using GRACE solutions can vary strongly from simple to twice, according to which provider (CSR, JPL, GFZ) is considered, and which

Africa (Grippa et al., 2011).

estimates (Zaitchik et al., 2008).

needs to be demonstrated.

**5.7 Detecting long-term variations over ice sheets** 

Many studies have demonstrated the possibility of satellite gravimetry to detect and monitor spatial versus time, at a precision of only tens of mm of EWH and spatial resolution of 400 km. Accuracy of the results still depends upon the level of noise in the GRACE data and the post-processing used. In addition to measurement errors and loss of signals by lowpass filtering, uncertainty on model-predicted GIA represents a major source of error in the mass balance estimates of the ice sheets melting, in particular for Antarctica. Besides, a resolution of 400 km still represents an important limitation for small-scale studies, especially in field hydrology, which often requires surface resolution of tens of kilometers. Improvements in pre- and post-processing techniques should increase the quality of the GRACE products, such as the "regional" or local approaches. Steady improvement of precision and resolution of GRACE products remains encouraging and it opens a wider class of applications than previously possible.

Even if the dying GRACE mission may not be able to provide data of the same accuracy to be exploitable in continental hydrology in 2012, other follow-on projects have already been proposed by a motivated community to ensure the continuity of such a satellite system. Satellite gravimetry remains the only remote sensing technique to map directly large scale mass variations. Any Follow-On GRACE mission in preparation by NASA and GFZ would be of first interest for studying the impacts on surface redistribution forced by climate change.

#### **7. Acknowledgment**

This review on the applications of GRACE to land water hydrology and ice sheets mass balance was funded by Université de Toulouse.

#### **8. References**


Alsdorf, D.E. & Lettenmaier, D.P. (2003). Tracking fresh water from space. *Science*, Vol.301,

Andersen, O. B. & Hinderer, J. (2005). Global inter-annual gravity changes from GRACE:

Andersen, O.B.; Seneviratne, S.I.; Hinderer, J. & Viterbo P. (2005). GRACE-derived terrestrial

Bergmann, I.; Ramillien, G. & Frappart, F. (2012). Climate-driven interrannual variations of

Bettadpur, S. (2007). *CSR Level-2 processing standards document for level-2 product release 04,* 

Biancamaria, S.; Cazenave, A.; Mognard, N.M.; Llovel, W. & Frappart, F. (2011). Satellite-

Bruinsma, S.; Lemoine, J.M.; Biancale, R. & Valès, N. (2010). CNES/GRGS 10-day gravity

Chen, J.L.; Wilson, C.R.; Blankenship, D.D. & Tapley, B. D. (2006). Antarctic mass rates from

Chen, J.L.; Wilson, C.R.; Tapley, B.D.; Blankenship, D.D. & Ivins, E.R. (2007). Patagonia Icefield

*Research Letters*, Vol.34, No.22, (November 2007), pp. L22501, ISSN 0094-8276. Chen, J.L.; Wilson, C.R. & Seo, K.W. (2009a). S2 tide aliasing in GRACE time-variable gravity solutions*. Journal of Geodesy*, Vol.83, No.7, (July 2009), pp. 679-687, ISSN 0949-7714. Chen, J.L.; Wilson, C.R.; Tapley B.D.; Zang, Z.L. & Nyu, G.Y. (2009b). 2005 drought event in the

*of Geophysical Research*, Vol.114, No.B5, (May 2009), pp. B05404, ISSN 0148-0227. Chen, J.L.; Wilson, C.R. & Tapley B.D. (2010). The 2009 exceptional Amazon flood and

Departamento Nacional da Produção Mineral/DNPM (1983). Mapa hidrogeológico do

*Research*, Vol.46, No.12, (December 2010), pp. W12526, ISSN 0043-1397. Dickey, J.O.; Bentley, C.R.; Bilham, R.; Carton, J.A.; Eanes J.A.; Herring, T.A.; Kaula, W.M.;

Early results. *Geophysical Research Letters*, Vol.32, No.1, (January 2005), pp. L01402,

water storage depletion associated with the European heat wave. *Geophysical Research Letters*, Vol.32, No.18, (September 2005), pp. L18405, ISSN 0094-8276. Baur, O.; Kuhn, M. & Featherstone, W.E. (2009). GRACE-derived ice mass variations over

Greenland by accounting for leakage effects. *Journal of Geophysical Research*, Vol.114,

the mass balance of Greenland, *Global and Planetary Change*, Vol.82-83, (February

*GRACE*. The GRACE project, Center for Space Research, University of Texas at

based high latitudes snow volume trend, variability and contribution to sea level over 1989-2006. *Global and Planetary Change*, Vol.75, No.3-4, (February 2011), pp. 99-

field models (release 2) and their evaluation. *Advances in Space Research*, Vol.45,

GRACE. *Geophysical Research Letters*, Vol.33, No.11, (June 2006), pp. L11502, ISSN

melting observed by Gravity Recovery and Climate Experiment (GRACE). *Geophysical* 

Amazon River basin as measured by GRACE and estimated by climate models*. Journal* 

interannual terrestrial water storage change observed by GRACE*. Water Resources* 

Lagerloef, G.S.E.; Rojstaczern S.; Smith, W.H.F.; van den Dool, H.; Wahr J.M. & Zuber, M.T. (1999). Gravity and the hydrosphere: new frontier. *Journal des Sciences Hydrologiques – Hydrological Sciences Journal*, Vol.44, No.3, (June 1999), pp. 407-415,

No.5639, (September 2003), pp. 1491−1494, ISSN 0036-8075.

No.B13, (June 2009), pp. B06407, ISSN 0148-0227.

No.4, (February 2010), pp. 587–601, ISSN 0273-1177.

2012), pp. 1-11, ISSN 0921-8181

Austin, pp. 327–742.

107, ISSN 0921-8181.

0094-8276.

ISSN 0262-6667.

Brasil, escala 1:5,000,000.

**8. References** 

ISSN 0094-8276.


Han, S.-C.; Shum, C.; Jekeli, C.; Kuo, C.; Wilson, C. & Seo, K.W. (2005). Non-isotropic Filtering

Horwath, M. & Dietrich, R. (2009). Signal and error in mass change inferences from GRACE:

Ivins*,* E.R. & James, T.S. (2005). Antarctic glacial isostatic adjustment: A new assessment. *Antarctic Science*, Vol.17, No.4, (December 2005), pp. 541-553, ISSN 0954-1020. Jekeli C. (1981). *Alternative methods to smooth the Earth's gravity field*. Report 327, Department of Geodesy Sciences and Survey, Ohio State University, Columbus, USA. Klees, R.; Revtova, E.A.; Gunter, B. C.; Ditmar, P.; Oudman, E.; Winsemius, H. C. & Savenije,

Klees, R.; Liu, X.; Wittwer, T.; Gunter, B. C.; Revtova, E.A.; Tenzer, R.; Ditmar, P.;

Kurtenbach, E.; Mayer-Gürr, T. & Eicker, A. (2009). Deriving daily snapshots of the Earth's

Leblanc, M. J.; Tregoning, P.; Ramillien, G.; Tweed, S.O. & Fakes, A. (2009). Basin-scale integrated

*Resources Research*, Vol.45, No.4, (April 2009), pp. W04408, ISSN 0043-1397. Lemoine, J.M.; Bruinsma, S.; Loyer, S.; Biancale, R.; Marty, J.C.; Perosanz, F. & Balmino, G.

*Space Research*, Vol.39, No.10, (May 2007), pp. 1620–1629, ISSN 0273-1177. Lemoine, F.G.; Luthcke, S.B., Rowlands, D.D.; Chinn, D.S.; Klosko, S.M. & Cox, C.M.

Longuevergne, L.; Scanlon, B.R. & Wilson, C.R. (2010). GRACE hydrological estimates for

Luthcke, S.B.; Arendt, A.A.; Rowlands, D.D.; McCarthy, J.J. & Larsen, C.F. (2008). Recent

*Resources*, Vol.46, No.11, (November 2010), pp. W11517, ISSN 0043-1397. Luthcke, S.B.; Zwally, H.J.; Abdalati, W.; Rowlands, D.D.; Ray, R.D.; Nerem, R.S.; Lemoine,

(November 2006), pp. 1286–1289, ISSN 0036-8075.

*Letters*, Vol.36, No.17, (September 2009), pp. L17102, ISSN 0094-8276. Kusche, J. (2007). Approximate decorrelation and non-isotropic smoothing of timevariable

071-6702-23-37, San Francisco, USA.

2008), pp. 335−359, ISSN 1573-0956.

2007), pp. 733−749, ISSN 0949-7714.

540-49349-5, Berlin, Germany.

pp. 849–864, ISSN 0956-540X.

0956-540X.

of GRACE Temporal Gravity for Geophysical Signal Enhancement. *Geophysical Journal International*, Vol.163, No.1, (September 2005), pp. 18-25, ISSN 0956-540X. Heiskanen, W.H. & Moritz H. (1967). *Physical Geodesy*. W.H. Freeman and Co Ltd, ISBN 978-

the case of Antarctica. *Geophysical Journal International*, Vol.177, No.4, (June 2009),

H.H.G. (2008a). The design of an optimal filter for monthly GRACE gravity models. *Geophysical Journal International*, Vol.175, No.2, (November 2008), pp. 417−432, ISSN

Winsemius, H. C. & Savenije, H.H.G. (2008b). A comparison of global and regional GRACE models for land hydrology. *Surveys in Geophysics*, Vol.29, No.4-5, (October

gravity field from GRACE L1B data using Kalman filtering. *Geophysical Research* 

GRACE-type gravity field models. *Journal of Geodesy*, Vol.81, No.11, (November

observations of the early 21st century multiyear drought in southeast Australia, *Water* 

(2007a). Temporal gravity field models inferred from GRACE data. *Advances in* 

(2007b). The use of mascons to resolve time-variable gravity GRACE, In: *Dynamic Planet: Monitoring and Understanding A Dynamic Planet with Geodetic and Oceanographic Tools*, P. Tregoning & C. Rizos, (Eds.), 231-236, Springer, ISBN 978-3-

small basins : evaluating approaches on the High Plains aquifer, USA. *Water* 

F.G.; McCarthy, J.J. & Chinn, D.S. (2006). Recent Greenland Ice Mass Loss by Drainage System from Satellite Gravity Observations. *Science*, Vol.314, No.5803,

glacier mass changes in the Gulf of Alaska region from GRACE mascon solutions. *Journal of Glaciology*, Vol.54, No.188, (December 2008), pp.767-777, ISSN 1727-5652.


Rodell, M. & Famiglietti, J.S. (2001). An analysis of terrestrial water storage variations in

Rodell, M. & Famiglietti, J.S. (2002). The potential for satellite-based monitoring of

2001), pp. 1327– 1339, ISSN 0043-1397.

L20504, ISSN 0094-8276.

0148-0227.

Illinois with implications for GRACE. *Water Resources Research*, Vol.37, No.5, (May

groundwater storage changes using GRACE: The High Plains aquifer, central US. *Journal of Hydrology*, Vol.263, No.1-4, (June 2002), pp. 245–256, ISSN 0022 1694. Rodell, M.; Houser, P.R.; Jambor, U.; Gottschalck, J.; Mitchell, K.; Meng, C.-J.; Arsenault, K.;

Cosgrove, B.; Radakovich, J.; Bosilovich, M.; Entin, J.K.; Walker, J.P.; Lohmann, D. & Toll, D. (2004a). The Global Land Data Assimilation System. *Bulletin of the American Meteorological Society*, Vol.85, No.3, (March 2004), pp. 381−394, ISSN 1520-0477. Rodell, M.; Famiglietti, J.S.; Chen, J.; Seneviratne, S.I.; Viterbo, P.; Holl, S. & Wilson, C.R.

(2004b). Basin scale estimates of evapotranspiration using GRACE and other observations. *Geophysical Research Letters*, Vol.31, No. 20, (October 2004), pp.

groundwater changes in the Mississippi River basin using GRACE. Hydrogeology

depletion in India. *Nature*, Vol.460, No.7257, (August 2009), 999-1003, ISSN 0028-0836.

Cox, C.M. & Andersen, O.B. (2005). Resolving mass flux at high spatial and temporal resolution using GRACE intersatellite measurements. *Geophysical Research* 

J.P. & Sabaka, T.J. (2010). Global mass flux solutions from GRACE: a comparison of parameter estimation strategies – Mass concentrations versus Stokes coefficients. *Journal of Geophysical Research*, Vol.115, No.B1, (January 2010), pp. B01403, ISSN

Cazenave, A.; Petrovic, S.; Jochmann, H. & Wünsch, J. (2006). GRACE observations of changes in continental water storage. *Global and Planetary Change*, Vol.50, No.1-2,

(2008). Hydrological signals observed by the GRACE satellites. *Surveys in* 

and Climate Experiment (GRACE) observed surface mass variations. *Journal of* 

storage variations from GRACE with groundwater-level measurements from the

Rodell, M.; Chen, J.; Kato, H.; Famiglietti, J.S.; Nigro, J. & Wilson C.R. (2007). Estimating

Rowlands, D.D.; Luthcke, S.B.; Klosko, S.M.; Lemoine, F.G.R.; Chinn, D.S.; McCarthy, J.J.;

Rowlands, D.D.; Luthcke, S.B.; McCarthy, J.J.; Klosko, S.M.; Chinn, D.S.; Lemoine, F.G.; Boy,

Sasgen, I.; Martinec, Z. & Fleming, K. (2006). Wiener optimal filtering of GRACE data. *Studia Geophysica et Geodaetica*, Vol.50, No.4, (October 2006), 499−508, ISSN 1573-1628. Schmidt,R.; Schwintzer, P.; Flechtner, F.; Reigber, C.; Güntner, A.; Döll, P.; Ramillien, G.;

Schmidt, R.; Flechtner, F.; Meyer, U.; Neumayer, K. -H.; Dahle, Ch.; Koenig, R. & Kusche, J.

*Geophysical Research*, Vol.112, No.B8, (August 2007), B08407, ISSN 0148-0227. Seo, K.W. & Wilson, C.R. (2005). Simulated estimation of hydrological loads from GRACE. *Journal of Geodesy*, Vol.78, No.7–8, (March 2005), pp. 442–456, ISSN 0949-7714. Seo, K.W.; Wilson, C.R.; Chen, J. & Waliser D. (2008). GRACE's spatial errors. *Geophysical Journal International*, Vol.172, No.3, (January 2008), pp.41-48, ISSN 0956-540X. Strassberg, G.; Scanlon, B.R. & Rodell, M. (2007). Comparison of seasonal terrestrial water

*Geophysics*, Vol.29, No.4-5, (October 2008), pp. 319−334, ISSN 1573-0956. Schrama, E. J. O.; Wouters, B. & Lavallée, D.A. (2007). Signal and noise in Gravity Recovery

Journal, Vol.15, No.1, (February 2007), pp. 159-166, ISSN 1431-2174. Roddell, M.; Velicogna, I. & Famiglietti, J. (2009). Satellite-based estimates of groundwater

*Letters*, Vol.32, No.4, (February 2005), pp. L04310, ISSN 0094-8276.

(February 2006), pp. 112–126, ISSN 0921-8181.

High Plains Aquifer (USA). *Geophysical Research Letters*, Vol.34, No. 14, (July 2007), pp. L14402, ISSN 0094-8276.


Recovery and Climate Experiment, and GPS satellite data. *Journal of Geophysical Research*, Vol.107, No.B10, (October 2002), pp. ETG 20-1, ISSN 0148-0227.


**Part 2** 

**Groundwater Modeling** 

214 Water Resources Management and Modeling

*Research*, Vol.107, No.B10, (October 2002), pp. ETG 20-1, ISSN 0148-0227. Velicogna, I. & Wahr, J. (2005). Greenland mass balance from GRACE. *Geophysical Research Letters*, Vol.32, No.18, (September 2005), pp. L18505, ISSN 0094-8276. Velicogna, I. & Wahr, J. (2006a). Measurements of Time-Variable Gravity Show Mass Loss in

Velicogna, I. & Wahr, J. (2006b). Acceleration of Greenland Ice Mass Loss in Spring 2004. *Nature*, Vol.443, No.7109, (September 2006), pp. 329–331, ISSN 0028-0836. Velicogna, I. (2009). Increasing rates of ice mass loss from the Greenland and Antarctic ice

Wahr, J.; Molenaar, M. & Bryan F. (1998). Time variability of the Earth's gravity field:

Wahr, J.; Swenson, S.; Zlotnicki, V. & Velicogna, I. (2004). Time variable gravity from

Werth, S.; Güntner, A.; Schmidt, R. & Kusche, J. (2009). Evaluation of GRACE filter tools

Wolf, M. (1969). Direct measurement of the Earth's gravitational potential using a satellite

Wouters, B. & Schrama, E.O.J. (2007). Improved accuracy of GRACE gravity solutions

Wouters, B.; Riva, R. E. M.; Lavallée, D. A. & Bamber J.L. (2011). Seasonal variations in sea

*Research Letters*, Vol.38, No.3, (February 2011), pp. L03303, ISSN 0094-8276. Xavier, L.; Becker, M.; Cazenave, A.; Longuevergne, L.; Llovel, W. & Rotunno Filho, O. C. (2010).

*Environment*, Vol.114, No.8, (August 2010), pp. 1629–1637, ISSN 0034-4257. Yeh, P.J.F.; Swenson, S.C.; Famiglietti, J.S. & Rodell, M. (2006). Remote sensing of

Zaitchik, B.F.; Rodell, M. & Reichle, R.H. (2008). Assimilation of GRACE Terrestrial Water

2009), pp. L19503, ISSN 0094-8276.

(December 2009), pp. 1499-1515, ISSN 0956-540X.

0148-0227.

L11501, ISSN 0094-8276.

ISSN 0148-0227.

ISSN 0094-8276.

W12203, ISSN 0043-1397.

Recovery and Climate Experiment, and GPS satellite data. *Journal of Geophysical* 

Antarctica. *Science*, Vol.311, No.5768, (March 2006), pp. 1754–1756, ISSN 0036-8075.

sheets revealed by GRACE. *Geophysical Research Letters*, Vol.36, No.19, (October

hydrological and oceanic effects and their possible detection using GRACE. *Journal of Geophysical Research*, Vol.103, No.B12, (December 1998), pp. 30205– 30229, ISSN

GRACE: First results. *Geophysical Research Letters*, Vol.31, No.11, (June 2004), pp.

from a hydrological perspective. *Geophysical Journal International*, Vol.179, No.3,

pair. *Journal of Geophysical Research*, Vol.74, No.22, (November 1969), pp. 5295-5300,

through empirical orthogonal function filtering of spherical harmonics, *Geophysical Research Letters*, Vol.34, No.23, (December 2007), pp. L23711, ISSN 0094-8276. Wouters, B.; Chambers, D. & Schrama, E.J.O. (2008). GRACE observes small-scale mass loss

in Greenland. *Geophysical Research Letters*, Vol.35, No.20, (October 2008), pp. L20501,

level induced by continental water mass: First results from GRACE. *Geophysical* 

Interannual variability in water storage over 2003–2008 in the Amazon Basin from GRACE space gravimetry, in situ river level and precipitation data. *Remote Sensing of* 

groundwater storage changes in Illinois using the Gravity Recovery and Climate Experiment (GRACE). *Water Resources Research*, Vol.42, No.12, (December 2006), pp.

Storage Data into a Land Surface Model: Results for the Mississippi River Basin. *Journal of Hydrometeorology*, Vol.9*,* No.3, (June 2008), pp. 535-548, ISSN 1525-755X.

## **Simplified Conceptual Structures and Analytical Solutions for Groundwater Discharge Using Reservoir Equations**

Alon Rimmer1 and Andreas Hartmann2

*1Israel Oceanographic and Limnological Research Ltd., The Yigal Alon Kinneret Limnological Laboratory, 2Institute of Hydrology, Freiburg University, 1Israel 2Germany* 

#### **1. Introduction**

The approaches to the study of hydrological issues are generally divided into two very different groups: (1) the physical approach; and (2) the system approach (Singh 1988). The physical approach is motivated primarily by scientific study and understanding of the physical phenomena, whereas the practical application of this knowledge to engineering and water resources management is recognized but not always fully required. Unlike detailed physical studies of each hydrological problem, the system approach is driven by the need to establish working relationships between measured parameters for solving practical hydrological problems. This approach simplifies the issue because it is unfeasible to consider the entire physical system. Therefore, a logical approach consists of measuring those variables in the hydrologic cycle, which appear significant to the problem, and establish explicit mathematical relationships between them.

An initial step and a well-recognized part of groundwater flow analysis is the definition of a conceptual model. It is usually a simplified perception of the dominant physical components of the studied groundwater system. The main purpose for constructing a conceptual model is concentrating on the parts relevant for solving the hydrological problem.

Common ways to convert a conceptual model of a groundwater system into mathematical formulations are reservoir (or 'tank') type models (Dooge, 1973; Sugawara, 1995). These model types are often used as a theoretical tool in surface and subsurface hydrology, for water management, control of inflows and outflows in lakes, rivers, reservoirs, and aquifers. The linear reservoir concept is an important component of many widely used hydrological models like the TOPMODEL (Beven & Kirby, 1979), HBV (Lindström et al., 1997) or WaSiM-ETH (Schulla & Jasper, 2007). Reservoir type models are especially useful in karst environments, because the essential information for physical approaches is usually not available (Jukic & Denic-Jukic, 2009). The lack of information and the necessity to use simplified reservoir models become evident in the high number of recently published studies on karst hydrology (Fleury et al., 2007; Geyer et al., 2008; Hartmann et al., 2011; Jukic & Denic-Jukic, 2009; Kessler & Kafri, 2007; Le Moine et al., 2008; Rimmer & Salingar, 2006; Tritz et al., 2011).

In this chapter, a set of typical groundwater modeling problems is described, exemplifying the use of simple reservoir structures to model spring discharge and/or groundwater level during time. In each example, we will explicate the use of the proposed reservoir type system. Moreover, in each case, we will examine an analytical solution associated with the proposed system using simple domain geometries. The advantage of analytical solutions is that their equations offer quick answers to the proposed mechanism based on a few basic parameters. These solutions therefore allow an immediate system understanding and provide a meaning value for each parameter or group of parameters. Given the differential equation that describes the groundwater system, most of the presented analytical solutions can be found using the 'symbolic mathematical toolbox' of MATLAB (http://www.mathworks.com).

#### **2. Examples**

Our set of example models include: (1) the classic formation of the linear reservoir problem for an aquifer drained by a single spring; (2) spring discharge potentially fed by two parallel aquifers; (3) spring discharge potentially fed by two serial groundwater aquifers; (4) two parallel aquifers with linear exchange and linear discharges; (5) the discharge from an aquifer with two outlets; (6) the discharge from an aquifer into a lake (submerged springs); and (7) the cases of long-term change of groundwater level and annual spring discharge. Although in most cases the models will be applied with a given set of measured data, it is important to clarify that these types of models are not location-specific, and can be used to model various groundwater flow systems.

#### **2.1 The formation of linear reservoir problem for a single spring discharge**

In a traditional hydrology, a spring discharge is often conceptually described and modelled using simple linear reservoirs. We can start the simplification of a system by examining the spring discharge Q (L3 T-1) according to Darcy's Law:

$$Q(t) = -k\_i G \frac{h(t) - H\_0}{\Delta x} \tag{1}$$

where h (with units of length, L) represents an equivalent unknown hydraulic head in the aquifer, H0 (L) is the head at the spring outlet (if an exact head can be evaluated) so that h-H0 represents the equivalent hydraulic head difference between two points, located at x (L) distance one from the other. The ki (units of length over time, L T-1) is the saturated hydraulic conductivity, and G (L2) represents an "equivalent" cross section of the flow. For practical purposes, it is assumed that ki, G and x are constant for a given natural aquifer, and therefore Eq. 1 can be simplified to:

$$Q(t) = a \cdot \left[ h(t) - H\_0 \right] \quad ; \quad a = -\frac{k\_i G}{\Delta x} \tag{2}$$

considering H0=0 in Eq. 2, further simplification can be conducted by conceptualising the drained aquifer as a reservoir (0) with storage V (L3) varying in time; constant recharge area A (L2) and a given effective porosity n (-):

In this chapter, a set of typical groundwater modeling problems is described, exemplifying the use of simple reservoir structures to model spring discharge and/or groundwater level during time. In each example, we will explicate the use of the proposed reservoir type system. Moreover, in each case, we will examine an analytical solution associated with the proposed system using simple domain geometries. The advantage of analytical solutions is that their equations offer quick answers to the proposed mechanism based on a few basic parameters. These solutions therefore allow an immediate system understanding and provide a meaning value for each parameter or group of parameters. Given the differential equation that describes the groundwater system, most of the presented analytical solutions can be found using the 'symbolic mathematical toolbox' of

Our set of example models include: (1) the classic formation of the linear reservoir problem for an aquifer drained by a single spring; (2) spring discharge potentially fed by two parallel aquifers; (3) spring discharge potentially fed by two serial groundwater aquifers; (4) two parallel aquifers with linear exchange and linear discharges; (5) the discharge from an aquifer with two outlets; (6) the discharge from an aquifer into a lake (submerged springs); and (7) the cases of long-term change of groundwater level and annual spring discharge. Although in most cases the models will be applied with a given set of measured data, it is important to clarify that these types of models are not location-specific, and can be used to

In a traditional hydrology, a spring discharge is often conceptually described and modelled using simple linear reservoirs. We can start the simplification of a system by examining the

> <sup>0</sup> ( ) ( ) *<sup>i</sup> ht H Qt kG*

where h (with units of length, L) represents an equivalent unknown hydraulic head in the aquifer, H0 (L) is the head at the spring outlet (if an exact head can be evaluated) so that h-H0 represents the equivalent hydraulic head difference between two points, located at x (L) distance one from the other. The ki (units of length over time, L T-1) is the saturated hydraulic conductivity, and G (L2) represents an "equivalent" cross section of the flow. For practical purposes, it is assumed that ki, G and x are constant for a given natural aquifer,

<sup>0</sup> ; *<sup>i</sup> k G Qt ht H*

considering H0=0 in Eq. 2, further simplification can be conducted by conceptualising the drained aquifer as a reservoir (0) with storage V (L3) varying in time; constant recharge area

*x*

(1)

*x*

(2)

 

**2.1 The formation of linear reservoir problem for a single spring discharge** 

MATLAB (http://www.mathworks.com).

model various groundwater flow systems.

and therefore Eq. 1 can be simplified to:

A (L2) and a given effective porosity n (-):

spring discharge Q (L3 T-1) according to Darcy's Law:

**2. Examples** 

$$V(t) = A \cdot n \cdot h(t) \tag{3}$$

according to Eqs. 2 and 3, in such a reservoir model, spring discharge through the outlet, Qout, is proportional to storage.

$$V(t) = KQ\_{out}(t) \quad \text{ :} \quad K = \frac{\left(A \cdot n\right)}{a} \tag{4}$$

where K (given in units of time, T) is known as the reservoir constant or storage, representing the recharge area, the porosity, the saturated hydraulic conductivity, and the equivalent path and cross section of the flow within the aquifer. Usually, changes of K in time or from one season to another are not physically justified, and it should be independent of both the selected period of modeling, and the boundary conditions (amount of precipitation).

The equation for the continuous water balance in this kind of reservoir is:

Fig. 1. (a) Schematic description of groundwater system; (b) linear reservoir model Incorporating Eq. 4 into Eq. 5 results in the linear reservoir differential equation:

$$K\frac{dQ\_{out}(t)}{dt} = Q\_{in}(t) - Q\_{out}(t) \quad ; \quad t\_0 \ge t \ge 0 \tag{6}$$

A well-known application in hydrology is the determination of K. This task becomes significantly easier in the dry period that follows the rainy season, since the flow is then a smooth, physical and unidirectional process, with no random processes (such as rainstorms) to be taken into account. At this time Qin(t)=0 and the mathematical description of linear reservoir model (Eq. 6) reduces to the homogeneous equation:

$$\frac{dQ\_{out}(t)}{dt} + \frac{Q\_{out}(t)}{K} = 0 \qquad ; \qquad Q\_{out}\left(t=0\right) = Q\_0 \tag{7}$$

Eq. 7 is solved analytically by:

$$\begin{aligned} a. \quad &Q\_{out}\left(t\right) = Q\_0 \exp\left(-\frac{t}{K}\right) \\ b. \quad &V\left(t\right) = K Q\_{out}\left(t\right) = K Q\_0 \exp\left(-\frac{t}{K}\right) \end{aligned} \tag{8}$$

In Eq. 8, V is the volume (assume 103 m3), t is the time (day), Qout is the outflow (103 m3 day-1), Q0 is the outflow (103 m3 day-1) at the day when Qin vanished, and K is the reservoir constant with units identical to the units of t (day).

Analyzing spring recession in this way is known as Maillet's approach (Maillet, 1905). An application of this fundamental method is presented in Fig. 2, with measured discharge flow from the Carcara Springs in the Western Galilee, Israel, during the dry period starting in March 1981. The springs emerge from the aquifer of the Upper Judea Group formation, which appears to be connected to the aquifer of the Lower Judea Group formations. In this time of the year, the regional groundwater level is usually high. Data from 1950-1985 indicated that the spring had never dried, a situation that changed significantly since the beginning of pumping in 1985 (These changes are discussed in section 2.5).

Fig. 2. The discharge of Carcara Spring during the dry period starting in March1981. K was calibrated to 117 days.

#### **2.2 Parallel linear reservoirs**

During a dry season that follows a rainy season, the discharge of a spring reduces in time. The shape of the graph discharge vs. time corresponds to the sum of several exponential functions (Bonacci, 1993; Grasso & Jeannin, 1994). Often, such spring discharge is represented as a combination of two parallel linear reservoirs (Fig. 3), mathematically represented by:

$$\begin{aligned} a. \quad Q\_{out}(t) &= Q\_{out1}\left(t\right) + Q\_{out2}\left(t\right) = \\ Q\_{01} \exp\left(-\frac{t}{K\_1}\right) &+ Q\_{02} \exp\left(-\frac{t}{K\_2}\right) \\ b. \quad V(t) &= K\_1 Q\_{out1}\left(t\right) + K\_2 Q\_{out2}\left(t\right) \end{aligned} \tag{9}$$

A simple optimization algorithm can be applied to identify the K1 and K2 constants, as well as the initial flows Q01 and Q02 during a recession period. When K is small, the recession of

0

Model Measured K=117 days (8)

(9)

 

beginning of pumping in 1985 (These changes are discussed in section 2.5).

.

.

Qout

calibrated to 117 days.

represented by:

**2.2 Parallel linear reservoirs** 

(1000 m3day-1)

*out*

with units identical to the units of t (day).

. exp

*<sup>t</sup> aQ t Q <sup>K</sup>*

0

. exp

*<sup>t</sup> b V t KQ t KQ <sup>K</sup>*

In Eq. 8, V is the volume (assume 103 m3), t is the time (day), Qout is the outflow (103 m3 day-1), Q0 is the outflow (103 m3 day-1) at the day when Qin vanished, and K is the reservoir constant

Analyzing spring recession in this way is known as Maillet's approach (Maillet, 1905). An application of this fundamental method is presented in Fig. 2, with measured discharge flow from the Carcara Springs in the Western Galilee, Israel, during the dry period starting in March 1981. The springs emerge from the aquifer of the Upper Judea Group formation, which appears to be connected to the aquifer of the Lower Judea Group formations. In this time of the year, the regional groundwater level is usually high. Data from 1950-1985 indicated that the spring had never dried, a situation that changed significantly since the

> 0 60 120 180 240 days

Fig. 2. The discharge of Carcara Spring during the dry period starting in March1981. K was

During a dry season that follows a rainy season, the discharge of a spring reduces in time. The shape of the graph discharge vs. time corresponds to the sum of several exponential functions (Bonacci, 1993; Grasso & Jeannin, 1994). Often, such spring discharge is represented as a combination of two parallel linear reservoirs (Fig. 3), mathematically

*out out out*

*b V t KQ t KQ t*

*aQ t Q t Q t*

1 2

1 2

 

A simple optimization algorithm can be applied to identify the K1 and K2 constants, as well as the initial flows Q01 and Q02 during a recession period. When K is small, the recession of

11 22

*out out*

exp exp

*t t Q Q K K*

01 02

*out*

the reservoir will be fast, and its discharge and volume will reach zero within a short time. When K is large, the recession will be slow, and the reservoir outflow will last for a long time. If K1 >> K2 the discharge from reservoir 2 (second component in the right hand side of Eq. 9) decreases much faster than the discharge from reservoir 1 (first component in the right hand side of Eq.9), which will still be active much longer.

Fig. 3. Two parallel reservoirs. The reservoirs are fed by groundwater recharge originating from the surface and drain simultaneously to the stream.

The parallel groundwater reservoir structure is incorporated in many hydrological models, e.g., the Vensim model (Fleury et al., 2007). Here it is exemplified by applying it to the recession discharge of the Hermon Stream (Israel) during the year of 1996 (Fig. 4). The stream is one of the three main tributaries of the Upper Jordan River (Rimmer and Salingar, 2006). It is fed mainly by the Banias Spring, located at the edge of the karst exposures on the lower parts of the Hermon Mountain, at an altitude of 359 m a.s.l. The Banias annual average discharge is ~67 M m³ (~2.15 m³/s). The spring exhibits behaviour of pluvio-nival regime, where discharge is mainly due to precipitation, but also slightly influenced by snowmelt (Gilad & Bonne, 1990; Samuels et al., 2010).

Fig. 4. Stream discharge of the Hermon stream during 1996. The stream is fed by two parallel reservoirs during the recession period: Since K2>>K1 the reservoir 1 discharge represents most of the sharp changes following the rainy season, while reservoir 2 represents the more stable component.

The optimization algorithm revealed that K1 = 56 days and Q01 = 420,000 m3day-1, representing the immediate aquifer that contributes to the spring, while K2 = 300 days and Q01 = 145,000 m3day-1, representing the discharge from a large stable aquifer, which also drains into the Dan Spring located nearby. During the recession period, the discharge of reservoir 1 ceases after ~170 days, while the memory of reservoir 2 remains for ~2.5 years (Rimmer & Salingar, 2006).

#### **2.3 Two serial linear reservoirs**

Usually, the discharge recession of a karst spring is fast at the beginning of a dry season and slows at its end (see Eqs. 8 and/or 9). However, there are cases in which the recession is rather slow at the beginning, and increases towards the end of the dry season (Rimmer & Salingar, 2006). Moreover, the recession is faster following a low precipitation season, than after a high precipitation one. One reason for such behaviour can be explained by the interplay of two systems in series, e.g., a large vadose zone on top of the phreatic zone or two groundwater systems of which one is recharged by leakage from the other (Fig. 5a). The pattern analysis of such measured spring discharge requires a different setup. The proposed mechanism for examining this type of observed curve is a system with two serial linear reservoirs (Fig. 5b).

Fig. 5. (a) schematic description of the proposed groundwater system; (b) the system represented by two serial reservoirs. Excess saturation flow from the earth surface feed the upper reservoir (1), which recharges the lower reservoir (2).

In this example, the simplified system is described by an upper linear reservoir, contributing to a lower reservoir, draining through a spring outlet. Similarly to the previous case, we are particularly interested in determining the system storage coefficients K1 and K2, and the initial conditions (flow at the beginning of the dry season) Q01 and Q02. By defining the input to the upper reservoir (1) during the dry season as zero, the input to the lower reservoir (2) is an exponential recession with time, typical to a linear reservoir system (Section 2.1; see also Nash 1957; Huggins & Burney, 1982). We therefore write the differential equation for the lower reservoir (2) for a single dry season as follows:

$$\begin{aligned} \frac{dQ\_{out2}(t)}{dt} + \frac{Q\_{out2}(t)}{K\_2} &= \frac{Q\_{out1}(t)}{K\_2} \quad ; \quad t \ge 0\\ \text{with} \quad 1. \ Q\_{out1} \left(t = 0\right) &= Q\_{01} \quad & 2. \ Q\_{out2} \left(t = 0\right) = Q\_{02} \; ; \end{aligned} \tag{10}$$

The optimization algorithm revealed that K1 = 56 days and Q01 = 420,000 m3day-1, representing the immediate aquifer that contributes to the spring, while K2 = 300 days and Q01 = 145,000 m3day-1, representing the discharge from a large stable aquifer, which also drains into the Dan Spring located nearby. During the recession period, the discharge of reservoir 1 ceases after ~170 days, while the memory of reservoir 2 remains for ~2.5 years

Usually, the discharge recession of a karst spring is fast at the beginning of a dry season and slows at its end (see Eqs. 8 and/or 9). However, there are cases in which the recession is rather slow at the beginning, and increases towards the end of the dry season (Rimmer & Salingar, 2006). Moreover, the recession is faster following a low precipitation season, than after a high precipitation one. One reason for such behaviour can be explained by the interplay of two systems in series, e.g., a large vadose zone on top of the phreatic zone or two groundwater systems of which one is recharged by leakage from the other (Fig. 5a). The pattern analysis of such measured spring discharge requires a different setup. The proposed mechanism for examining this type of observed curve is a system with two serial linear

Fig. 5. (a) schematic description of the proposed groundwater system; (b) the system represented by two serial reservoirs. Excess saturation flow from the earth surface feed the

In this example, the simplified system is described by an upper linear reservoir, contributing to a lower reservoir, draining through a spring outlet. Similarly to the previous case, we are particularly interested in determining the system storage coefficients K1 and K2, and the initial conditions (flow at the beginning of the dry season) Q01 and Q02. By defining the input to the upper reservoir (1) during the dry season as zero, the input to the lower reservoir (2) is an exponential recession with time, typical to a linear reservoir system (Section 2.1; see also Nash 1957; Huggins & Burney, 1982). We therefore write the differential equation for

1. 0 2. 0 ;

1 01 2 02

; 0

(10)

upper reservoir (1), which recharges the lower reservoir (2).

the lower reservoir (2) for a single dry season as follows:

*dQ t Q t Q t <sup>t</sup> dt K K*

2 2

*out out*

*with Q t Q Q t Q*

221

*out out out*

(Rimmer & Salingar, 2006).

reservoirs (Fig. 5b).

**2.3 Two serial linear reservoirs** 

where Q01 and Q02 are the initial conditions, yet to be determined from the measured data of each season. In Eq. 10 the contribution from the upper to the lower groundwater reservoir and the upper reservoir volume are determined by:

$$\begin{aligned} Q\_{out1}\left(t\right) &= Q\_{01} \exp\left(-\frac{t}{K\_1}\right) \\ V\_1\left(t\right) &= K\_1 Q\_{out1} \end{aligned} \tag{11}$$

and Eq. 10 can be solved analytically so that the discharge from the lower groundwater reservoir and its volume are determined by:

$$\begin{aligned} Q\_{out2}\left(t\right) &= Q\_{out2/1} + Q\_{out2/2} \\ Q\_{out2/1} &= \frac{Q\_{01}K\_1}{\left(K\_1 - K\_2\right)} \Bigg[ \exp\left(\frac{-t}{K\_1}\right) - \exp\left(\frac{-t}{K\_2}\right) \Bigg] \\ Q\_{out2/2} &= Q\_{02} \exp\left(\frac{-t}{K\_2}\right) \\ V\_2\left(t\right) &= K\_2 Q\_{out2} \end{aligned} \tag{12}$$

Here, the outflow from the lower reservoir is combined of the contribution from the upper reservoir Qout2/1 and the self-discharge of the lower reservoir Qout2/2. With an optimization algorithm, Eq. 12 may be used to evaluate K1, K2, Q01 and Q02 for each season, so that it

Fig. 6. Illustration of the terms in Eq. 12- the upper reservoir contribution Qout2/1 and the self-discharge of the lower reservoir Qout2/2 combine the total outflow Qout2. The K1 = 70 days and K2 =300 days are identical for both rainy (1993) and dry (1990) years, while the two initial conditions Q01 and Q02 are different. (a, b): 1993; (c, d): 1990. (a, c): spring discharge; (b, d): Aquifer volume.

would match the measured spring discharge. Two restrictions should be imposed on the calibration procedure in order to take into account the physical conditions of the entire system. First, the same K1 and K2 must be used for all seasons, and second, there should be a good correlation between Q01 and Q02 and the annual precipitation during the years, since the entire system is driven by the same precipitation.

Fig. 6 illustrates the curve fitting of Eq. 12 to the discharge of the Dan Spring, Israel, during the dry season that followed two different rainy seasons. The K1 = 70 days and K2 = 300 days were evaluated as the best fit. The initial conditions of Q01=1900 m3 day-1 and Q02 =800 m3 day-1 were valid following a very rainy season (1992-1993), while Q01=16 m3 day-1 and Q02 =580 m3 day-1 were valid for extremely dry season (1989-1990). Following the rainy season, both reservoirs were partly filled according to the amount of precipitation. However, while the recession of the lower reservoir follows the same rate (exp(-t/K2)) under any initial condition, the additional recharge from the upper reservoir changes significantly the Qout2(t) curve during the dry season. Consequently, the flow rate of the spring may increase first, following a very rainy year (1992-1993) or reduce immediately following a very dry year (1989-1990). Similar applications can be found in Kiraly, (2003) or Rimmer & Salingar, (2006).

#### **2.4 Two reservoirs with linear exchange**

The karst environment is often described as a system with dual porosity (Goldscheider & Drew, 2007), including the fast flow component within the preferential flow paths (karstic conduits), and the slower Darcian groundwater flow within the fissure matrix (Fig. 7a). This process can be conceptualized by dividing the groundwater system in two reservoirs, one representing the conduits and the other representing the fissure matrix (Fig. 7b). Similar to section 2.2, the water exchange between the reservoirs is controlled by the difference in their levels and with similar considerations as in Eqs. 1-4 the spring discharge Q1 (L³T-1) (or the conduit outflow) may be derived by:

$$Q\_1(t) = \frac{V\_1(t)}{K\_1} \tag{13}$$

Applying the same procedure to the exchange flow between the fissure matrix and the conduit reservoir (reservoirs 2 and 1 in in Fig. 7), and aggregating all constant parameters in a single exchange constant KE(T), the exchange flow QE (L3T-1) can be described as a simple linear relation between water level differences (Fig. 7),

$$Q\_E\left(t\right) = k\_{12}G\_{12}\frac{h\_2\left(t\right) - h\_1\left(t\right)}{\Delta x\_{12}} = \frac{V\_2\left(t\right) - f\_P V\_1\left(t\right)}{K\_E} \tag{14}$$

Where:

$$h\_2(t) = \frac{V\_2(t)}{A n\_2} \qquad \frac{1}{K\_E} = \frac{k\_{12} G\_{12}}{A n\_2 \Delta x\_{12}} \qquad f\_P = \frac{n\_2}{n\_1} \tag{15}$$

Similar to Eqs. 1-4, h2 (L) is the water table elevation of the fissure matrix, k12 (L T-1) is a representative saturated hydraulic conductivity, G12 (L2) is an equivalent cross-section, and

would match the measured spring discharge. Two restrictions should be imposed on the calibration procedure in order to take into account the physical conditions of the entire system. First, the same K1 and K2 must be used for all seasons, and second, there should be a good correlation between Q01 and Q02 and the annual precipitation during the years, since

Fig. 6 illustrates the curve fitting of Eq. 12 to the discharge of the Dan Spring, Israel, during the dry season that followed two different rainy seasons. The K1 = 70 days and K2 = 300 days were evaluated as the best fit. The initial conditions of Q01=1900 m3 day-1 and Q02 =800 m3 day-1 were valid following a very rainy season (1992-1993), while Q01=16 m3 day-1 and Q02 =580 m3 day-1 were valid for extremely dry season (1989-1990). Following the rainy season, both reservoirs were partly filled according to the amount of precipitation. However, while the recession of the lower reservoir follows the same rate (exp(-t/K2)) under any initial condition, the additional recharge from the upper reservoir changes significantly the Qout2(t) curve during the dry season. Consequently, the flow rate of the spring may increase first, following a very rainy year (1992-1993) or reduce immediately following a very dry year (1989-1990). Similar applications can be found in Kiraly, (2003) or Rimmer & Salingar,

The karst environment is often described as a system with dual porosity (Goldscheider & Drew, 2007), including the fast flow component within the preferential flow paths (karstic conduits), and the slower Darcian groundwater flow within the fissure matrix (Fig. 7a). This process can be conceptualized by dividing the groundwater system in two reservoirs, one representing the conduits and the other representing the fissure matrix (Fig. 7b). Similar to section 2.2, the water exchange between the reservoirs is controlled by the difference in their levels and with similar considerations as in Eqs. 1-4 the spring discharge Q1 (L³T-1) (or the

<sup>1</sup>

Applying the same procedure to the exchange flow between the fissure matrix and the conduit reservoir (reservoirs 2 and 1 in in Fig. 7), and aggregating all constant parameters in a single exchange constant KE(T), the exchange flow QE (L3T-1) can be described as a simple

21 2 1

<sup>2</sup> 12 12 <sup>2</sup>

Similar to Eqs. 1-4, h2 (L) is the water table elevation of the fissure matrix, k12 (L T-1) is a representative saturated hydraulic conductivity, G12 (L2) is an equivalent cross-section, and

1

*E*

2 2 12 1

12

*h t h t V t fV t Q t kG x K*

1 *V t*

*Q t <sup>K</sup>* (13)

*P*

*P*

*V t kG n h t <sup>f</sup> An K An x n* (15)

(14)

*E*

1

the entire system is driven by the same precipitation.

**2.4 Two reservoirs with linear exchange** 

conduit outflow) may be derived by:

linear relation between water level differences (Fig. 7),

*E*

2

12 12

(2006).

Where:

x12 (L) an average flow distance; all are parameters representing the interface between conduits and fissure matrix.

Fig. 7. (a) Schematic description of the groundwater system with karstic conduits and fissure matrix; (b) the system represented by a reservoir combination with a fissure matrix reservoir (left, 2) and conduit reservoir (right, 1).

With the effective porosity of the fissure matrix n2, and the area A2, the relation between water level and stored water volume V2 (L3) can be established. Note that in Eqs. 13-15 the same area A is used because the simplification approach assumes that the conduits section is embedded within the fissure matrix (double porosity approach) and that the porosity differences between the conduits and fissure matrix were taken into account by the porosity factor fP. With it A, A1 and A2 in Fig. 7 are related to each other as follows:

$$A = A\_1 + A\_2 = A\_1(1 + f\_P) \tag{16}$$

Having defined the flow processes of the conduit and the fissure matrix, water balance for both reservoirs can be established:

$$\begin{aligned} \frac{dV\_1(t)}{dt} &= Q\_{IN1} + \frac{V\_2(t) - f\_p V\_1(t)}{K\_E} - \frac{V\_1(t)}{K\_1} \\ \frac{dV\_2(t)}{dt} &= Q\_{IN2} - \frac{V\_2(t) - f\_p V\_1(t)}{K\_E} \end{aligned} \tag{17}$$

Rearranging Eq. (17) results in a linear system of inhomogeneous differential equations:

$$
\begin{pmatrix} D + \frac{1}{K\_1} + \frac{f\_P}{K\_E} & -\frac{1}{K\_E} \\\\ -\frac{f\_P}{K\_E} & D + \frac{1}{K\_E} \end{pmatrix} \begin{pmatrix} V\_1(t) \\\\ V\_2(t) \end{pmatrix} = \begin{pmatrix} Q\_{IN1} \\\\ Q\_{IN2} \end{pmatrix} \tag{18}
$$

Hereby, D is the differential operator d/dt. Assuming constant inflows QIN1 and QIN2, Eq. (18) can be solved analytically with standard methods (Kramer's rule, variation of constants; e.g. Boyce and DiPrima, 2000) and yield:

$$\begin{aligned} V\_1(t) &= B\_1 \exp(A\_1 t) + B\_2 \exp(A\_2 t) + C\_1 \\ V\_2(t) &= B\_3 \exp(A\_1 t) + B\_4 \exp(A\_2 t) + C\_2 \end{aligned} \tag{19}$$

With the constants

$$A\_{1,2} = -\frac{1}{2} \left( \frac{1}{K\_1} + \frac{1 + f\_P}{K\_E} \right) \pm \sqrt{\frac{1}{4} \left( \frac{1}{K\_1} + \frac{1 + f\_P}{K\_E} \right)^2} - \frac{1}{K\_1 K\_E} \tag{20}$$

$$\begin{aligned} B\_1 &= \frac{V\_{10}' - A\_2 V\_{10} + A\_2 C\_1}{A\_1 - A\_2} \\ B\_3 &= \frac{V\_{20}' - A\_2 V\_{20} + A\_2 C\_2}{A\_1 - A\_2} \end{aligned} \tag{21}$$

$$B\_4 = \frac{V\_{20}' - A\_2 V\_{20} + A\_2 C\_2}{A\_1 - A\_2} \quad \qquad \qquad \qquad B\_4 = V\_{20} - B\_3 - C\_2$$

$$\begin{aligned} \mathbf{C}\_1 &= K\_1 \left( Q\_{IN1} + Q\_{IN2} \right) \\ \mathbf{C}\_2 &= K\_E Q\_{IN2} + K\_1 f\_P \left( Q\_{IN1} + Q\_{IN2} \right) \end{aligned} \tag{22}$$

Where V10 and V20 are the reservoir volumes and *V*10 and *V*20 the storage change at t= 0. *V*10 and *V*20 can be obtained by Eq. 18):

$$\begin{split} V\_{10}' &= \frac{dV\_1 \left(t=0\right)}{dt} = Q\_{IN1} + \frac{V\_{20} - f\_p V\_{10}}{K\_E} - \frac{V\_{10}}{K\_1} \\ V\_{20}' &= \frac{dV\_2 \left(t=0\right)}{dt} = Q\_{IN2} - \frac{V\_{20} - f\_p V\_{10}}{K\_E} \end{split} \tag{23}$$

Except for A1,2, the constants, as well as the initial conditions refer to a single time step, and have to be calculated each time step again. For instance at time step t V10 would be equal to V1(t-1) and V20 equal to V2(t-2), respectively.

Methods that consider the exchange between fissure matrix and conduits can be found in Cornaton & Perrochet (2002) and Sauter (1992). In Fig. 8, the exchange reservoirs solution was applied to the last recharge event and the dry season recession in 1998 at the Banias Spring (see section 2.2). The exchange between the conduit and the fissure matrix reservoir resulted in a buffering of the recharge signal and a slow increase in fissure matrix storage. Exchange flow was negative, indicating flow towards the fissure matrix. Around the end of

Fig. 8. Left: stored water in the conduit reservoir V1, the fissure matrix reservoir V2 and total storage V1+V2; Middle: total recharge QIN, spring discharge Q1 and exchange flow QE vs. observed spring flow; Right: conduit and matrix water SO4 concentrations, CConduits and CMatrix, discharge concentrations C1 vs. observed SO4 concentrations

April, it changed its direction, which means that parts of the stored water in the fissure matrix were released again to the springs. This switch of direction of the exchange flow was nearly insignificant in terms of flow rates but had immense impact on the water quality. This is exemplified by the SO4 variations observed at the same spring during the same time: by simply attributing constant SO4 concentrations to the conduit and matrix flows their mixing at the spring outlet resulted in an acceptable agreement with the observations.

#### **2.5 The linear reservoir with two outlets**

226 Water Resources Management and Modeling

 

1 2 10 1 1

*V AV AC <sup>B</sup> BV BC*

3 4 20 3 2

*V AV AC <sup>B</sup> BV BC*

*C KQ Kf Q Q*

11 1 2

*C KQ Q*

0

0

*dV t V fV V Q dt <sup>K</sup>*

10 1

20 2

2 21 1 2 *IN IN*

*E IN P IN IN*

1 20 10 10

2 20 10

*IN*

*dV t V fV <sup>V</sup> V Q dt K K*

*IN*

Except for A1,2, the constants, as well as the initial conditions refer to a single time step, and have to be calculated each time step again. For instance at time step t V10 would be equal to

Methods that consider the exchange between fissure matrix and conduits can be found in Cornaton & Perrochet (2002) and Sauter (1992). In Fig. 8, the exchange reservoirs solution was applied to the last recharge event and the dry season recession in 1998 at the Banias Spring (see section 2.2). The exchange between the conduit and the fissure matrix reservoir resulted in a buffering of the recharge signal and a slow increase in fissure matrix storage. Exchange flow was negative, indicating flow towards the fissure matrix. Around the end of

Fig. 8. Left: stored water in the conduit reservoir V1, the fissure matrix reservoir V2 and total storage V1+V2; Middle: total recharge QIN, spring discharge Q1 and exchange flow QE vs. observed spring flow; Right: conduit and matrix water SO4 concentrations, CConduits and

CMatrix, discharge concentrations C1 vs. observed SO4 concentrations

2 4

10 2 10 2 1

1 2 20 2 20 2 2

*A A*

1 2

*A A*

Where V10 and V20 are the reservoir volumes and *V*<sup>10</sup>

can be obtained by Eq. 18):

V1(t-1) and V20 equal to V2(t-2), respectively.

*f f <sup>A</sup>*

1 11 1 1 1 1 1 1 1

*K K K K KK*

*E EE*

and *V*<sup>20</sup>

*P*

*E P*

*E*

1

*P P*

2

(20)

(21)

(22)

the storage change at t= 0.

(23)

With the constants

*V*<sup>10</sup>

and *V*<sup>20</sup>

1,2

In this section, the case of the effect of additional outlet is discussed. Consider the case of a spring discharge, which differs from the basic case (section 2.1) in two important elements: (1) the spring may dry out completely, so that the exponential recession (Eq. 8) is not valid for low flow rates; and (2) From the water mass balance calculations, it is assumed that the groundwater recharge is larger than the spring discharge, and therefore part of the water continues to flow downstream the aquifer to deeper layers. When these two conditions are valid, the linear reservoir with upper and lower outlets (Fig. 9) may represent the system rather well.

With similar considerations, we can handle the problem in Eq. 5 with two outlets, and no recharge, described as follows:

$$\frac{dV\left(t\right)}{dt} = \begin{cases} Q\_{out0} + Q\_{out1} & ; \quad 0 \le t \le t\_1 \\ Q\_{out0} & ; \quad t \ge t\_1 \end{cases} \tag{24}$$

where t1 is the time in which the upper outlet (Qout1) is drying, leaving only the flow in the lower outlet (Qout0).

Fig. 9. Linear reservoir with two outlets at different levels

If it is assumed that outlet (1) changes the pressure field only locally, we can consider each outflow separately as a linear function of the head above it, so that:

$$\begin{aligned} Q\_{out1}(t) &= -\alpha\_1 \cdot \left( h(t) - H\_1 \right) \\ Q\_{out0}(t) &= -\alpha\_2 \cdot \left( h(t) - H\_0 \right) \end{aligned} \tag{25}$$

With Eq. 25 incorporated into Eq. 24, assuming no inflow and H0=0 the problem is defined as follows:

$$\frac{dh\,\mathrm{h}\,(t)}{dt} = \begin{cases} -\left(\frac{1}{K\_1} + \frac{1}{K\_2}\right)h(t) + \left(\frac{1}{K\_1}H\_1\right) & ; \quad 0 \le t \le t\_1\\ -\frac{1}{K\_2}h(t) & ; \quad t\_1 \le t \end{cases}; \quad \text{with} \quad h\,(t=0) = h\_0 \tag{26}$$

The analytical solution to the problem in Eq. 26 is:

$$\begin{aligned} h(t) &= \frac{\mathbb{C}}{q} + \left(h\_0 - \frac{\mathbb{C}}{q}\right) \exp\left(-qt\right) \quad \quad ; \quad 0 \le t \le t\_1 \\ &\text{where} \\ \mathbb{C} = \frac{H\_1}{K\_1} \quad ; \quad q = \frac{1}{K\_1} + \frac{1}{K\_2} \\ &\text{and} \\ h(t) &= \mathbb{C} \exp\left(-\frac{1}{K\_2}(t - t\_1)\right) \quad ; \quad \qquad t \ge t\_1 \end{aligned} \tag{27}$$

In order to keep continuous recession curve, the head at time t1 should be equal to H1, and identical for the two problems, therefore:

$$h(t\_1) = \frac{C}{q} + \left(h\_0 - \frac{C}{q}\right) \exp\left(-qt\_1\right) = H\_1 \tag{28}$$

From Eq. 28 we can define the time t1 in which the flow from the upper outlet vanishes. It is a function of 0 11 2 *h H K and K* , , .

$$\left.t\_1 = t\right|\_{h=H\_1} = -\ln\left(\frac{H\_1 - \mathcal{C}/q}{h\_0 - \mathcal{C}/q}\right)\Big/q\tag{29}$$

That type of groundwater reservoir is also included in the HBV model (Lindström et al., 1997). An application of the proposed mechanism is presented in Fig. 10, with measured discharge flow from the Carcara Springs in the Western Galilee, Israel, during the dry period starting in March 2002. Note that this spring is similar to the one presented in section 2.1 and therefore K1 was calibrated to 117 days. However, there is a major difference between section 2.1 and 2.5; since the beginning of pumping in 1985, water levels have been dropping significantly, so that the spring has been drying completely almost every dry season since 1995. The drying requires analysing the spring discharge with Eq. 26 rather than with Eq. 7, and additional calibration of K2=100 days, which, as expected, turned to be nearly similar to K1. On the regular scale (Fig. 10a) the difference between a drying and not drying spring is not easily perceived, but it becomes clear, and the value of t1 = 161 days is obvious when plotted on a logarithmic scale (Fig. 10b).

Fig. 10. The outflows through the upper outlet in a linear reservoir with two outlets in different level Qout1 and, for comparison, the outflow from a regular linear reservoir without lower outlet Qout1,exp . (a) regular scale; (b) logarithmic scale.

#### **2.6 Aquifer drainage to submerged springs**

228 Water Resources Management and Modeling

With Eq. 25 incorporated into Eq. 24, assuming no inflow and H0=0 the problem is defined

1 1

0 1

exp ; 0

1 1

1

1 1

*C C ht h qt tt*

;

*ht H t t dh t KK K with h t h*

; 0

0

(27)

(28)

(29)

(26)

as follows:

2

The analytical solution to the problem in Eq. 26 is:

1

*K*

*dt*

12 1

1

*where*

*and*

identical for the two problems, therefore:

a function of 0 11 2 *h H K and K* , , .

11 1 ; 0

*h t t t*

1 1 2

2

*q q* 

1

ln *h H H Cq t t <sup>q</sup> h Cq* 

1

obvious when plotted on a logarithmic scale (Fig. 10b).

*ht C t t t t K*

In order to keep continuous recession curve, the head at time t1 should be equal to H1, and

 1 0 1 1 exp *C C ht h qt H*

From Eq. 28 we can define the time t1 in which the flow from the upper outlet vanishes. It is

That type of groundwater reservoir is also included in the HBV model (Lindström et al., 1997). An application of the proposed mechanism is presented in Fig. 10, with measured discharge flow from the Carcara Springs in the Western Galilee, Israel, during the dry period starting in March 2002. Note that this spring is similar to the one presented in section 2.1 and therefore K1 was calibrated to 117 days. However, there is a major difference between section 2.1 and 2.5; since the beginning of pumping in 1985, water levels have been dropping significantly, so that the spring has been drying completely almost every dry season since 1995. The drying requires analysing the spring discharge with Eq. 26 rather than with Eq. 7, and additional calibration of K2=100 days, which, as expected, turned to be nearly similar to K1. On the regular scale (Fig. 10a) the difference between a drying and not drying spring is not easily perceived, but it becomes clear, and the value of t1 = 161 days is

1

0

 

1 exp ;

;

*<sup>H</sup> C q K KK*

*q q*

The same physical factors were considered in modelling the process of groundwater discharge into springs onshore and offshore a lake, or a river (Fig. 11). Unlike in previous cases, here, the analytical solution was applied to the entire annual cycle, in order to exemplify the case where the spring outflows are dictated by the downstream head at the lake or river.

Fig. 11. (a) Schematic description of groundwater system that discharges simultaneously from confined aquifer into springs onshore and offshore a lake. (b) A model where the reservoir drains through a constant level (z1) onshore spring Qout1(t), and a time varying boundary HL(t) representing offshore spring Qout0(t). The hydraulic head within a short distance (up to several hundred meters) from the lake is h(t).

The proposed simplified model aims to link the time-dependent spring discharge to the hydraulic heads in the contributing aquifer under the fluctuations of the lake level. These fluctuations are independent of the aquifer system, and affect the spring flow as a close boundary condition. Here we assume constant recharge Qin, and time dependent discharge from the aquifer to the onshore Qout1(t) and offshore Qout0(t) springs:

$$\begin{aligned} Q\_{out1}(t) &= \alpha\_1 \cdot \left[ h(t) - z\_1 \right] \\ Q\_{out0}(t) &= \alpha\_0 \cdot \left[ h(t) - H\_L(t) \right] \end{aligned} \tag{30}$$

We regard the elevation of the onshore spring z1=0, and lake HL level fluctuating as a sin function around an average level. Under these conditions:

$$\begin{aligned} Q\_{out1}(t) &= \alpha\_1 \cdot h(t) \\ Q\_{out0}(t) &= \alpha\_0 \cdot \left( h(t) - \left\lceil b\_0 + b\_1 \sin\left(\beta\_1 t + \gamma\right) \right\rceil \right) \end{aligned} \tag{31}$$

Here b0 (m) is the average lake level below the spring outlet (b0 is negative); b1 (m) is the lake fluctuations amplitude; is the angular frequency in radians which for yearly rotations has a set value of 1 2 365.25 ; and is the phase shift (radians). With similar considerations as in earlier problems, we can handle the reservoir mass balance with two outlets, and constant recharge:

$$\frac{dV\{t\}}{dt} = Q\_{in} - Q\_{out0} - Q\_{out1} \tag{32}$$

Incorporating Eq. 31 into Eq. 32 and rearranging using Eqs. 2, 3 and 4 results in:

$$\begin{aligned} \frac{dh(t)}{dt} + \left(\frac{1}{K\_1} + \frac{1}{K\_0}\right)h(t) &= \frac{Q\_{in}}{A \cdot n} + \frac{b\_0}{K\_0} + \frac{b\_1}{K\_0} \Big[\sin\left(\beta\_1 t + \gamma\right)\Big] \\\text{with } : \qquad h\left(t = 0\right) = H\_0 \end{aligned} \tag{33}$$

This equation is solved analytically by:

$$\begin{aligned} h(t) &= H\_0 \exp(-qt) + \frac{b\_0}{K\_0 q} + \frac{Q\_{in}}{A \cdot n \cdot q} - b\_1 \frac{\beta \cos(\beta\_1 t + \gamma) + q \sin(\beta\_1 t + \gamma)}{K\_0 \left(\beta\_1^2 + q^2\right)} \\\ q &= \frac{1}{K\_0} + \frac{1}{K\_1} \end{aligned} \tag{34}$$

An initial test of this solution reveals that if *t* , the lake level assumed to be steady with no fluctuations (b1=0), and the inflow Qin=0, then:

$$h(t) = \frac{b\_0}{K\_0 q} = \frac{b\_0}{\left(1 + \frac{K\_0}{K\_1}\right)}\tag{35}$$

From Eq. 35 it can be concluded that if the connection between the aquifer and the lake is significantly stronger than the connection to the onshore spring (K0<<K1), then the aquifer hydraulic head assumes the level of the lake <sup>0</sup> *ht b* , but if K0>>K1 the aquifer hydraulic head adapts to the level of the spring outlet *h t* 0 . If Qin>0 then h(t) increases by Qin resulting higher discharge through the spring outlets. Discharge of an onshore spring Qout1(t) is straight forward to measure. Therefore, we can evaluate it according to Eq. 31 and calibrate 1. However the offshore spring discharge Qout2(t) is usually difficult to measure resulting in infinite possibilities to evaluate it since 0 is also unknown.

We regard the elevation of the onshore spring z1=0, and lake HL level fluctuating as a sin

*in out out* 0 1

0 1

0 1 1

*Kq Anq K q*

  

; and is the phase shift (radians). With similar

*QQ Q dt* (32)

1

 

(35)

 (31)

(33)

(34)

0 0 01 1 sin

Here b0 (m) is the average lake level below the spring outlet (b0 is negative); b1 (m) is the lake fluctuations amplitude; is the angular frequency in radians which for yearly rotations

considerations as in earlier problems, we can handle the reservoir mass balance with two

cos sin exp( )

An initial test of this solution reveals that if *t* , the lake level assumed to be steady with

0 0

*b b h t*

resulting in infinite possibilities to evaluate it since 0 is also unknown.

0 0

*K q K*

From Eq. 35 it can be concluded that if the connection between the aquifer and the lake is significantly stronger than the connection to the onshore spring (K0<<K1), then the aquifer hydraulic head assumes the level of the lake <sup>0</sup> *ht b* , but if K0>>K1 the aquifer hydraulic head adapts to the level of the spring outlet *h t* 0 . If Qin>0 then h(t) increases by Qin resulting higher discharge through the spring outlets. Discharge of an onshore spring Qout1(t) is straight forward to measure. Therefore, we can evaluate it according to Eq. 31 and calibrate 1. However the offshore spring discharge Qout2(t) is usually difficult to measure

1

1

*K*

 

0 1 2 2 0 0 1

*in b Q tqt h t H qt <sup>b</sup>*

1 1 sin

1 0 0 0 0

*in dh t Q b <sup>b</sup> h t <sup>t</sup>*

*Q t ht b b t*

function around an average level. Under these conditions:

*out out*

 2 365.25 

This equation is solved analytically by:

0 1

no fluctuations (b1=0), and the inflow Qin=0, then:

1 1

*<sup>q</sup> K K*

has a set value of 1

outlets, and constant recharge:

Incorporating Eq. 31 into Eq. 32 and rearranging using Eqs. 2, 3 and 4 results in:

*dt K K A n K K*

*dV t*

: 0

*with h t H*

1 1

*Q t ht*

As an example, the analytical solution is applied to the Fuliya Springs (Fig. 12) onshore and offshore lake Kinneret Israel. The case of the Fuliya saline springs was classified as confined carbonate aquifer, interacting with the lake through fractures and faults (Goldshmidt et al., 1967; Gvirtzman et al., 1997; Bergelson et al., 1998). The carbonate aquifer system of these springs overlays deep-seated brine, from which saline flux mixes with the fresh groundwater. Diluted saline water drains through fracture springs to both onshore and offshore springs (Rimmer et al., 1999, Abbo et al., 2003). Hydrogeological studies of this natural group of springs, as well as their intensive monitoring (Rimmer et al., 1999) allow us to analyze the simultaneous discharge processes of both onshore and offshore springs in more detail. The observations show that the measured hydraulic head of the aquifer and the discharge to the onshore springs follows the fluctuations (increase or decrease) of the measured lake level (Fig. 12). Discharge to offshore springs could not be measured directly. There is however clear evidence (Simon & Eizik, 1991) that it behaves as a "mirror" picture of the lake level. These results were previously verified by a partial analytical solution proposed by Rimmer et al., (1999) and later by a detailed numerical model (Abbo et al., 2003). With the current analytical solution in Eq. 34 we can test the offshore and the total discharge in time by changing 0 (Fig. 13). The 'real' value of Qout2(t) remain however unknown.

Fig. 12. Application of the analytical solution (Eq. 34) to a. the measured Lake Kinneret level and the measured hydraulic head in the aquifer ~100 m from the lake, and to b. Fuliya Spring discharge through onshore spring. Discharge to the onshore spring vanishes when the hydraulic head drops below the level of the spring outlet.

Fig. 13. Application of the analytical solution (Eq. 34 and 31) to Fuliya Spring discharge through both onshore and offshore springs, with three different values of 0 (a: 0=-2, b: 0=- 10, c: 0=-20).

#### **2.7 Long term reduction of groundwater level and spring discharge**

The same physical considerations may be used to examine the process of long-term changes of groundwater level and annual spring discharge (Fig. 14). Unlike previous cases, the time scale here is much larger than a daily scale. The analytical solution is applied here for multi annual changes of hydrological variables such as groundwater level and spring discharge, to test whether the aquifer storage is affected by the initiation of large changes upstream. Such changes are for example the initiation of pumping wells, or construction of large water storage reservoirs, which started at a certain point in time.

We consider an average constant annual inflow Qin0 to the reservoir that represents the aquifer storage. The outflow is similar to Eq. 2, where elevation of the spring outlet is set to H0=0. A constant flow Qp represents an outflow from the reservoir in addition to the spring outlet, such as pumping wells or reduction of inflows due to significant land use changes. Under these definitions:

$$\begin{aligned} Q\_{in}\left(t\right) &= Q\_{in0} \\ Q\_{out}\left(t\right) &= \alpha\_0 \cdot \left[ h\left(t\right) - H\_0 \right] \\ Q\_p\left(t\right) &= Q\_p \end{aligned} \tag{36}$$

Fig. 14. (a) Schematic description of groundwater system; (b) linear reservoir model -the water flux through the outlet is proportional with storage.

With similar considerations as in the problems described above, the reservoir mass balance is controlled by two outflows (Fig. 14) – one is constant in time Qp, whereas the other is time-dependent spring discharge Qout(t). On the annual time scale, the natural recharge Qin0 is considered as constant. The time when the change occurred (pump, land use change) is considered as t=0. The reservoir equation is therefore:

$$\begin{aligned} \frac{dV(t)}{dt} &= Q\_{in0} - a \left[ h(t) - H\_0 \right] - Q\_p \\ t < 0 \to Q\_p &= 0 \qquad ; \quad t \ge 0 \to Q\_p > 0 \end{aligned} \tag{37}$$

Rearranging the problem results in

$$\begin{cases} \frac{dh(t)}{dt} = -\frac{1}{K}h(t) + \frac{Q\_{in0} - Q\_p}{A \cdot n} \\ t < 0 \to Q\_p = 0 \end{cases}; \quad \text{with} \quad h(t = 0) = h\_0 \tag{38}$$

The same physical considerations may be used to examine the process of long-term changes of groundwater level and annual spring discharge (Fig. 14). Unlike previous cases, the time scale here is much larger than a daily scale. The analytical solution is applied here for multi annual changes of hydrological variables such as groundwater level and spring discharge, to test whether the aquifer storage is affected by the initiation of large changes upstream. Such changes are for example the initiation of pumping wells, or construction of large water

We consider an average constant annual inflow Qin0 to the reservoir that represents the aquifer storage. The outflow is similar to Eq. 2, where elevation of the spring outlet is set to H0=0. A constant flow Qp represents an outflow from the reservoir in addition to the spring outlet, such as pumping wells or reduction of inflows due to significant land use changes.

> 0

*Q t ht H*

Fig. 14. (a) Schematic description of groundwater system; (b) linear reservoir model -the

0 0

*dV t Q ht H Q dt tQ tQ*

0 0; 0 0 *in p*

 

<sup>0</sup>

*dh t Q Q h t with h t h*

*in p*

0 0; 0 0

*dt K A n tQ tQ*

*p p*

*p p*

<sup>1</sup> ; 0

With similar considerations as in the problems described above, the reservoir mass balance is controlled by two outflows (Fig. 14) – one is constant in time Qp, whereas the other is time-dependent spring discharge Qout(t). On the annual time scale, the natural recharge Qin0 is considered as constant. The time when the change occurred (pump, land use change) is

0 0

**h**

**<sup>H</sup> Qout <sup>0</sup>**

**A**

**Qp**

**K**

**Qin**

(36)

(b)

(37)

(38)

0

*Qt Q*

*in in*

( )

*Qt Q*

*p p*

*out*

**H0 Qout**

**A**

**<sup>h</sup>** (a)

impermeable

water flux through the outlet is proportional with storage.

considered as t=0. The reservoir equation is therefore:

Rearranging the problem results in

**Qin**

**2.7 Long term reduction of groundwater level and spring discharge** 

storage reservoirs, which started at a certain point in time.

Under these definitions:

**Qp**

It should be emphasized that K in this case represents a timescale by far larger than the seasonal timescale. Eq. (38) is solved analytically with

$$\begin{split} h(t) &= \left[ h\_0 - \frac{K}{A \cdot n} (Q\_{in0} - Q\_p) \right] \exp\left( \frac{-t}{K} \right) + \frac{K}{A \cdot n} (Q\_{in0} - Q\_p) \\ Q\_{out}(t) &= \frac{A \cdot n}{K} h(t) \end{split} \tag{39}$$

It is assumed that prior to t=0, a steady state had been reached with Qp=0 and therefore h(t)=h0 =(K/An)Qin0. At t the system is set on a new steady state h(t)= (K/An)(Qin0- Qp). The expression [(h0- K/An)(Qin0-Qp)] is the aquifer system long-term full response to the change Qp in water inflows and outflows. If this expression is zero, aquifer level and spring discharge will remain unchanged in time. If the expression is positive, hydraulic head decrease from one steady state to another, and vice versa.

As an example, this analytical solution is applied to the groundwater level in the Lower Judea Group Aquifer near the Uja spring, located in the Eastern Basin of the Judea-Samaria Mountains, ~12 km north-west of the town of Jericho. According to water level measurements and the stratigraphic analysis in this region (Guttman, 2007; Laronne Ben-Itzhak & Gvirtzman, 2005), the Judea Group aquifer, with a thickness of about 800 to 850 m, is comprised of two sub-aquifers: the upper and the lower aquifers. The upper and lower sub-aquifers are separated by relatively low permeability formations, causing groundwater levels in the upper aquifer to be significantly higher than those in lower aquifer do.

Near the Uja Spring there are four wells. (Mekorot Uja-Na'aran wells 1,2,3,4) drilled into the lower aquifer (Guttman, 2007). The first well (Uja 1) was drilled in 1964 by the Jordanian authority to a depth of 288 m and later was deepened by the Israeli authorities to 536 m. This well pumped from the upper part of the lower aquifer. In 1974, a new well (Uja 2) was drilled to a depth of 615 m in order to replace the Uja 1 well. At the beginning of the 1980s, two more wells were drilled (Uja 3, to a depth of 738 m and Uja 4 to 650.5 m) a few kilometres south of the other two wells. The three wells (Uja 2,3,4) currently pump ~3×106 m3 annually from the lower aquifer of the Judea Group.

It is assumed that the steady state of groundwater levels in the wells stood at 100 m below sea level (bsl.) prior to the significant pumping in 1974, whereas currently, the new steady state is ~280 m bsl. The long-term measured reduction of groundwater level is nearly exponential from 1974 to 1991 (Fig. 15). Following the extremely rainy season of 1991-1992 the levels increased to ~220 bsl, but since the year 2000 it returned to the steady state of ~280 m bsl. The proposed solution for this case was reached assuming K=1980 days; Qin=8200 m3 day-1 (3×106 m3 annually); Qp=8200 m3 day-1; A×n=90,000 m2; and reduction of level (t=0) initiated in 1974.

The physical interpretation of these results is that prior to the year 1974 a flux of ~3×106 m3 passed through the local Lower Judea Group aquifer annually (both Qin and Qout were ~ 8200 m3 day-1). The continuous pumping caused a significant reduction of groundwater level, and brought the system to a new steady state in which the natural flow of groundwater is reduced. The artificial deployment replaced the natural groundwater outflow, which originally travelled downstream following the hydraulic gradient.

Fig. 15. Measured and modeled multiannual ground water level in the Lower Judea Group aquifer near Uja Spring from 1974 to 2007.

Fig. 16. Measured and modeled multiannual trends of (a) groundwater level (monthly values from 1988 to1999), and (b) average annual discharge of the Masrefot Spring from 1952 to 2008.

Another simplified solution can be obtained for the special case in which long-term regional groundwater level and spring discharge is being constantly reduced. One possible explanation for such reduction is the increasing deployment of the aquifer. From the modelling point of view this is a case where Qp in Eqs. 36-39 follows a linear change in time (Qp = at+ b). We consider an average constant annual inflow Qin and outflow Qout similar to Eq. (2), but the additional outflow Qp, representing local pumping, is evolving and increasing linearly in time. When Qp = at+ b is implemented in equations (36-39) the analytical solution is:

$$\begin{aligned} h(t) &= \left( h\_0 + \frac{\left( Kb - KQ\_{in} \right)}{A \cdot n} \right) \cdot \exp\left( -\frac{t}{K} \right) + \frac{K}{A \cdot n} \{Q\_{in} - (at + b)\} \\ Q\_{out}(t) &= \frac{A \cdot n}{K} h(t) \end{aligned} \tag{40}$$

In this case, the exponential term (first term in the right hand side of Eq. 40) may approach zero very quickly, while most of the reduction of groundwater level and spring discharge depends on the decrease of recharge expressed in the second term in Eqs. 40 by Qin–(at+b). At large t the system continuously reduces in time as expected.

As an example, this analytical solution is applied to the Masrefot Spring (Fig. 16), which is affected by the hydraulic heads of the Lower Judea Group aquifer in the Western Galilee, Israel. This spring was selected for this case since on the one hand, according to Kafri (1970), its seasonal changes in discharge are hardly noticed due to the large storage of water that feeds them. On the other hand, the long-term history of measured spring discharge (~60 years, Fig. 16) may reflect the reduction of regional groundwater level. The regional water supply system includes dozens of pumping wells. Analysis of the actual annual pumping in these wells revealed a nearly linear increase of pumped water (*r*2=0.87) at least between 1960 and 1990, with an average increase of ~350,000 m3 year-1. This is probably the reason for the systematic linear decrease of groundwater levels and Masrefot Spring discharge.

#### **3. Summary**

234 Water Resources Management and Modeling

01/1974 01/1978 01/1982 01/1986 01/1990 01/1994 01/1998 01/2002 01/2006 **date**

Fig. 15. Measured and modeled multiannual ground water level in the Lower Judea Group

Jan-52 Sep-65 May-79 Jan-93 Oct-06

Fig. 16. Measured and modeled multiannual trends of (a) groundwater level (monthly values from 1988 to1999), and (b) average annual discharge of the Masrefot Spring from

**date**

Measured discharge Modeled discharge

Modeled level

Measured groundwater level Trend of groundwater level

**a.**

**b.**

Exp function

Mek Uja-Na'aran 1 Mek Uja-Na'aran 2 Mek Uja-Na'aran 3 Mek Uja-Na'aran 4


**Annual discharge (m**

1952 to 2008.

**3day-1)**

**Groundwater level (m)**

aquifer near Uja Spring from 1974 to 2007.




**Groundwater level (m)**



The steps towards modeling groundwater usually include (1) definition of the modeling domain; (2) definition of the hydrogeological structure; and (3) evaluation of initial and boundary conditions. The objective of this paper is to suggest an additional step (4), estimating the dominant parts that define the timely response of the hydrological system. This step is particularly important in developing conceptual models that simplify the hydrological problem to its relevant processes. We showed that using an analytical solution with this methodology could result in some important understanding of the system in question. Although the analytical solution can sometimes be the entire required modeling, usually the usage of analytical solution is only the first idea that we have on the timedependent system. The cases described in section 2.1 and 2.2 are indeed well known, and often found in the literature. However, the cases described in sections 2.3-2.7 are less familiar, but can be used for creating other new types of models. In section 2.3, we showed that a system of two serial reservoirs might be used to characterize the flow instead of parallel reservoirs. In section 2.4, the possibility of exchanging parallel reservoirs was discussed. Section 2.5 proposed that the recession curve might not fall into the well-known exponential shape due to downward flow to lower outlet, while section 2.6 showed how spring flow could change in time due to nearby dominant boundary condition (such as lake). Finally, section 2.7 suggested simple modeling for multiannual groundwater levels and spring discharge under reduction in water availability, a specific environmental problem that is occurring now, and expected in the future. Altogether, our set of examples can help in developing new process-based models for better system understanding and forward prediction.

#### **4. References**


parallel reservoirs. In section 2.4, the possibility of exchanging parallel reservoirs was discussed. Section 2.5 proposed that the recession curve might not fall into the well-known exponential shape due to downward flow to lower outlet, while section 2.6 showed how spring flow could change in time due to nearby dominant boundary condition (such as lake). Finally, section 2.7 suggested simple modeling for multiannual groundwater levels and spring discharge under reduction in water availability, a specific environmental problem that is occurring now, and expected in the future. Altogether, our set of examples can help in developing new process-based models for better system understanding and

Abbo, H., U. Shavit, D. Markel, and A. Rimmer, (2003). A Numerical Study on the Influence

Bakalowicz, M., 2005. Karst groundwater: a challenge for new resources. Hydrogeology

Bergelson, G., Nativ, R. & Bein, A. 1998. Assessment of hydraulic parameters in the aquifers sorrounding and underlying Sea of Galilee. GroundWater 36:409-417. Beven, K.J., Kirby, M.J., 1979. A physically based, variable contributing area model of basin

Bonacci, O., 1993. Karst springs hydrographs as indicators of karst aquifers / Les

Boyce, W.E., DiPrima, R.C., 2000. Elementary Differential Equations. Wiley, New York, 700

Cornaton, F., Perrochet, P., 2002. Analytical 1D dual-porosity equivalent solutions to 3D

Dooge, J.C.I. 1973. Linear theory of hydrologic systems. US Dept. Agric. Tech. Bull. No.

Fleury, P., Plagnes, V., Bakalowicz, M., 2007. Modelling of the functioning of karst aquifers

Geyer, T., Birk., S., Liedl, R., Sauter, M., 2008. Quantification of temporal distribution of

Goldscheider, N., Drew, D., 2007. Methods in Karst Hydrogeology. Taylor & Francis Group,

Goldshmidt, M.J., Arad A., Neev, D., 1967. The mechanism of the saline springs in the Lake

Guttman, J. 2007. The Karstic Flow System in Uja Area – West Bank: An Example of two

[Eds.] Water Resources in the Middle East. Springer. ISBN 9783540695097.

hydrogrammes des sources karstiques en tant qu'indicateurs des aquifères

discrete-continuum models. Application to karstic spring hydrograph modelling.

with a reservoir model: Application to Fontaine de Vaucluse (South of France).

recharge in karst systems from spring hydrographs. Journal of Hydrology, 348:

Tiberias depression. Min. Dev. Geol. Surv., Jerusalem, Hydrol. Pap. #11, Bull. 45. 19

Separated Flow Systems in the Same Area. Chapter6 in: Shuval H. and Dweik, H.

Journal of Hydrology, 283/1-4 pp. 225-243.

hydrology. Hydrological Sciences Bulletin, 24: 43-69.

karstiques. Hydrological Sciences Journal, 38(1): 51 - 62.

Journal, 13: 148-160.

pp, ISBN 9780470039403

1468, pp. 267–293.

452-463.

pp.

Journal of Hydrology, 262: 165-176.

Journal of Hydrology, 345: 38-49.

264 p., ISBN 9780415428736

of Fractured Regions on Lake/Groundwater Interaction; the Lake Kinneret Case.

forward prediction.

**4. References** 


Singh, V.P., 1988. Hydrologic systems, rainfall-runoff modeling. Prentice Hall, NJ.

