A Mathematical View of Water Table
Fluctuations in a Shallow Aquifer in Brazil
by Dagmar C. Neto1 , Hung K. Chang2 , and Martinus Th. van Genuchten3

Abstract
Detailed monitoring of the groundwater table can provide important data about both short- and long-term aquifer processes,
including information useful for estimating recharge and facilitating groundwater modeling and remediation efforts. In this paper,
we presents results of 4 years (2002 to 2005) of monitoring groundwater water levels in the Rio Claro Aquifer using observation wells
drilled at the Rio Claro campus of São Paulo State University in Brazil. The data were used to follow natural periodic fluctuations in
the water table, specifically those resulting from earth tides and seasonal recharge cycles. Statistical analyses included methods oftime-series analysis using Fourier analysis, cross-correlation, and R/S analysis. Relationships could be established between rainfall
and well recovery, as well as the persistence and degree of autocorrelation of the water table variations. We further used numericalsolutions of the Richards equation to obtain estimates of the recharge rate and seasonable groundwater fluctuations. Seasonablesoil moisture transit times through the vadose zone obtained with the numerical solution were very close to those obtained withthe cross-correlation analysis. We also employed a little-used deep drainage boundary condition to obtain estimates of seasonable
water table fluctuations, which were found to be consistent with observed transient groundwater levels during the period of study.

Introduction
Groundwater monitoring constitutes an important
tool for aquifer management by providing information
about the current state and trends of various hydrological processes, as well as recording responses to bothanthropogenic activities and natural phenomena. Detailed
monitoring further provides insight into the nature ofgroundwater table fluctuations in response to such processes as precipitation, evapotranspiration and recharge,
regional flow and pumping, barometric pressure changes,
earth tides, and even earthquakes (Todd 1959; Bredehoeft
1967; Roeloffs 1988; Healy and Cook 2002; Zeumann
et al. 2009; Cutillo and Bredehoeft 2011; Chen et al.
2013; and numerous references therein). Understanding
the processes that cause variations in the water table
and the ability to monitor groundwater levels is critical for
estimating recharge, whether on the scale of an aquifer toassess groundwater resources and management in general,
1
Department of Applied Geology, São Paulo State University,
UNESP, Rio Claro, SP 13506-900, Brazil; dagmarcn@rc.unesp.br2
Department of Applied Geology, São Paulo State University,
UNESP, Rio Claro, SP 13506-900, Brazil; chang@rc.unesp.br3
Corresponding author: Department of Mechanical Engineering, LTTC/COPPE, Federal University of Rio de Janeiro, Riode Janeiro, RJ 21945-970, Brazil; +55-21-3264-7025; rvangenuchten@hotmail.com
Received August 2014, accepted February 2015.
© 2015, National Ground Water Association.
doi: 10.1111/gwat.12329

82

or for local situations such as for assessing groundwatercontamination.
This paper presents the results of detailed monitoringof the groundwater table from 2002 to 2005 in wellsof the Rio Claro Aquifer, drilled on the campus ofSão Paulo State University (UNESP) in the city of Rio
Claro in the Brazilian state of São Paulo. The data
were analyzed statistically using a combination of publicdomain software and software specifically developed forthis study. We used methods of time-series analysis(Fourier analysis, cross-correlation, and R/S analysis) toverify the periodicity of natural cycles such as thosestemming from seasonal precipitation and earth tides,
and to estimate relationships between rainfall and wellrecharge. The persistence and degree of autocorrelationof the water table variations were also evaluated. We
further used one-dimensional (1D) vadose zone numericalflow modeling with a special deep drainage boundary
condition to obtain estimates of the recharge rate, soil
moisture transit times in the unsaturated zone, andresultant seasonable water table fluctuations. Results showthat detailed monitoring provides valuable information
contributing to a more comprehensive understanding of
subsurface flow processes. The information may further
facilitate remediation efforts for situations where verticaltransport of contaminants, especially volatile constituents,
is enhanced by frequent variations in the groundwatertable (e.g., Neeper 2001; Stafford and Rixey 2011; Iglauerand Muggeridge 2013).

Vol. 54, No. 1–Groundwater–January-February 2016 (pages 82–91)

NGWA.org

Study Location
Groundwater levels were monitored in wells of the
Rio Claro Aquifer, drilled for this purpose on the campus
of São Paulo State University in the city of Rio Claro,
Brazil. The Rio Claro Aquifer is a relatively shallow
unconfined aquifer with groundwater tables being quite
variable, generally between depths of about 6 and 10 m
at the study site (Oliva et al. 2005).The aquifer consists
mostly of Cenozoic sedimentary rocks of fluvial origin
in the Rio Claro Formation (DAEE 1981). Sediments
are composed of fine- to medium-grained sands with anabundant clay matrix, with the thickness of the aquifer
varying from a few meters to up to 50 m. Recharge occurs
in most of the aquifer, while discharge areas are locatednear the contact of the formation with Permian clay-richsediments of the Corumbataí Formation of the ParanáBasin, as well as in narrow valleys.
Two wells located 80 m apart were drilled to penetrate
and pass through the Rio Claro Aquifer for detailed
monitoring of the groundwater table. The first well (IGCE3) was drilled to a depth of 17 m at the end of March,
2001. The well was lined with polyvinyl chloride (PVC)
casing, with a section of the screen/mesh extending froma depth of 5.6 m to the bottom of the well. The watertable at the time of drilling was 7.74 m below the soil
surface. This first well was monitored from April 2001
to April 2004, at which time the monitoring devices weretransferred to another well (IGCE-6). That second well
was drilled in September 2003 down to a depth of 20.1 m
to the aquitard of the Corumbataí Formation, with thescreen extending from a depth of 7.1 m to the bottom of
the well. The water table in this second well at the timeof drilling was 9.8 m below the soil surface.

Water Table Fluctuations and Pluviometric
Data
Groundwater levels were monitored using a watercolumn probe connected to a computer recording thewater level every 5 s. Automated measurements wereobtained using a pressure transducer having a range of
5 psi, equivalent to an approximately 3.5 m column ofwater. The transducer was vented to allow barometriccorrections and installed below the water table in thewell. The accuracy of the probe was 0.1% of its fullrange. Electrical signals from the probe were converted
into digital values using a metering device equipped witha standard RS-232 interface for computer data transfer.
The data were registered using a microcomputer running
Wellplex, a software package capable of monitoring up
to eight simultaneous devices, developed by the Basin
Studies Laboratory (LEBAC) at UNESP University in Rio
Claro.
Pluviometric data were obtained from daily total
rainfall records maintained by the meteorological stationof the Center for Environmental Studies and Planning(CEAPLA) at UNESP, located approximately 500 m fromthe observation wells. Figure 1 shows a plot of the
NGWA.org

Figure 1. Tsoft corrected plots of the seasonal cycles ofrecharge and drawdown of the aquifer from 2002 through2005 (top), and of daily rainfall rates (bottom). The intensityand duration of the rainy season influenced the recharge rate
as reflected in the water levels, as well as the interval between
the minimum and maximum water levels. The dashed line in
the top plot indicates the date of probe relocation from well
IGCE-3 to well IGCE-6.

precipitation data, as well as of measured water table data
after time-series preprocessing to allow for corrections
and filtering. Much of the preprocessing was done using
the interactive Tsoft software package of Van Camp andVauterin (2005). This software was found to be very
useful for pre-treatment of our data, including for isolating
unusual factors such as the presence of outliers, and for
identifying very short temporal variations such as theeffects of earthquakes on groundwater levels.
One difficulty during data acquisition was the presence of time-series gaps caused by temporary instrument
failure. However, throughout 4 years of monitoring, thetotal data loss was less than 20%, which were filled bylinear interpolation. This seems acceptable given the factthat the maximum interval of data loss was less than 3 dout of a total of 1460 d (4 years).
A second problem requiring correction involved theabrupt offset in values resulting from probe relocationfrom the first to the second well. This relocation wasnecessary to maintain the probe submerged and within
the specified metering range. The shifts in position werecorrected by adding (or subtracting) a constant value to (orfrom) all points in the time series after the first displacedvalue. This value was easily calculated from the difference
between the value of the last point prior to offset andthe measurement immediately afterwards. After making
these corrections, we adopted an arbitrary reference datum
situated at 633 m above mean sea level (approximately
9.5 m below the soil surface).
The original 5-s sampling interval resulted in asequence of more than 20 million data points. These points
were filtered to reduce the data to generate smaller files
that were easier to manage. A new sample interval of
24 h was subsequently adopted to match the interval of
D.C. Neto et al. Groundwater 54, no. 1: 82–91

83available pluviometric data. The 5-s sampling interval was
initially selected with the intent to record high-frequency
phenomena such as those resulting from earthquakes(Todd 1959; Roeloffs 1988; Zeumann et al. 2009; Chenet al. 2013). We indeed detected several earthquakes with
our monitoring system, including some from as far away
as Alaska (results not further shown here).

Observed Seasonal Variations and Earth Tides
Figure 1 shows the major cycles identified during
4 years of monitoring the groundwater table (from January
2002 to December 2005). Since monitoring involved two
different wells covering different periods of time, data
were concatenated into a single curve. Seasonal variationsin the water table are clearly a result of variations inannual rainfall, with the maximum levels increasing from2002 to 2004, but decreasing in 2005. The response of
the groundwater table is a result of recharge, which inturn is a function not only of the amount of rainfall but
also of its distribution in time. For example, the longer
duration of the summer rainy season (December to March)
in 2003/2004 led to more recharge and a higher peak in
the water table, whereas in early 2005 some of the rains
were heavier but lasted less time, leading to less recharge.
The effects of earth tides on the water table data areshown in Figure 2. Although the effects of tides are more
obvious on the water levels of oceans and other surface
water bodies when compared with those of subsurface
systems, forces generated by gravitational fields docause deformation of rocks, leading to earth tides of
the type shown in Figure 2. The effect of earth tides onwater levels in wells in both unconfined and especiallyconfined aquifers has been well documented (Todd 1959;
Bredehoeft 1967; Roeloffs 1988; Rojstaczer and Agnew
1989; Neeper 2001; Zeumann et al. 2009, among others).
According to Bredehoeft (1967), observations of fluctuations in water levels caused by earth tides have existed
since the end of the 19th century, while Todd (1959)
cited the work of Robinson (1939) describing variationsin the water level of wells in Iowa and New Mexico.
The effects of earth tides on the groundwater tablewere clearly visible in our data (Figure 2), which is not
surprising given the sensitivity of the instruments we used.
Two cycles occurred each day, with the greatest amplitude
occurring at approximately 3 or 4 pm, and the smallest in
the early morning at approximately 3 or 4 am. As notedby Todd (1959), the effect of these tides during new andfull moons is quite regular, with the amplitude of thevariations being greater than those during the first and lastquarters when the cycles are less regular due to opposing
gravitational forces resulting from the positions of the sunand moon relative to Earth.

Analysis of Time Series
The water table data were analyzed statistically using
several methods of time-series analysis: Fourier analysis,
84

D.C. Neto et al. Groundwater 54, no. 1: 82–91

Figure 2. Oscillations in the water table due to earth tides
for several days in late February and early March 2003,
showing two maxima and two minima each day.

cross-correlation, and R/S analysis. The methods andresults are detailed in the following sections.
Fourier Analysis
A Fourier analysis may be used to reveal thefrequencies of periodic and non-periodic functions, withthe frequencies of the periodic functions determined by thecoefficients in the Fourier series. Fourier analyses makeit possible to analyze the components of the time seriesof the water table in the frequency mode, which in termsof our study should confirm the presence of two different
types of periodic variations: one due to the season of
the year involving annual recharge and drawdown of theaquifer, and one caused by earth tides producing two
daily fluctuations. As our study involved data sampled
at discrete time intervals, a discrete Fourier transform(Bracewell 1986; Fleming et al. 2002; Si 2007) was used:
F (n) =

N−1
f (t) e−iωn t

(1)

t=0

where ωn = 2π n/N for n = 0, 1, . . . , N − 1, resulting ina discrete set of Fourier frequencies. When this set of
frequencies is plotted, the resulting graph is a periodogramor spectrogram, which identifies the contribution of thefrequency of each of the periodic phenomena to theformation of the original signal.
For analysis of the seasonal water table variations,
we used the filtered and corrected data time series.
Calculation of the spectrum of frequencies of the periodic
process using the discrete Fourier transform generated thespectrogram shown in Figure 3, portrayed on a log scale
to enhance the spectral peaks. The horizontal coordinateis plotted in terms of cycles per day. Since seasonalvariations occur only once a year, their frequency is
1/365 cycles/d, or approximately 0.0027 cycles/d. Figure 3
shows that the peak of the spectrogram of the entire timeseries indeed occurred at or very close to this value.
For analysis of the earth tides, only a portion of 2003
was used, corresponding to the peak of the seasonal cycle
NGWA.orgbeing investigated (the cause), and a second series, yt (y 1 ,y 2 , . . . , yn ), the output (the result). The cross-correlation
function is then given by the expressions (Padilla and
Pulido-Bosch 1995):
r+k = rxy (k) = 

Cxy (k)r−k = ryx (k) = 

where

(2)

Cx2 (0) Cy2 (0)Cyx (k)(3)

Cx2 (0) Cy2 (0)1(xt − x) (yt+k − y)
n

(4)

1
(yt − y) (xt+k − x)
n

(5)

1
(xt − x)2
n

(6)

1
(yt − y)2
n

(7)

n−k

Figure 3. Spectrogram of the complete water-level timeseries from 2002 to 2005.

Cxy (k) =

t=1

n−k

Cyx (k) =t=1

n

Cx (0) =t=1nCy (0) =t=1

Figure 4. Spectrum of tidal oscillations showing a peak for
two cycles per day related to the tides.

when the effect of tide was clearer. The groundwatertable at that time was very stable and did not reflect
much influence from either recharge or drawdown. The
spectrogram is shown in Figure 4. The larger peak offrequency for a periodicity of 2 cycles/d and the smaller
peak for a single cycle per day reflect the result of thetwice-daily tide (K2) and the daily tide (K1), associatedwith the revolution of the Earth–Moon system around thecenter of mass and the rotation of the Earth around its
own axis, respectively (Pugh 1987).
Cross-Correlation Analysis
Another useful method for analyzing temporal series
in hydrogeological studies is cross-correlation quantifying
the dependency of successive values over a period oftime. Padilla and Pulido-Bosch (1995) and Larocque
et al. (1998) used this method in their studies of thehydrodynamic characteristics of karst aquifers, including
both spatial and temporal variations, as well as the relationof the hydrographs of rivers and creeks to discharge zones,
estimation of specific storage and porosity, and the delay
between input and output signals (such as rainfall andpiezometric levels).
Consider two discrete time series, with one series,
xt (x 1 , x 2 , . . . , xn ), representing the input to the system
NGWA.org

in which x and y are the means of the xt and yt series,
respectively, k = 0, 1, 2, . . . , m, with k being the time
lag from one element to the next one, and where m isthe cutting point that determines the interval over which
the cross-correlation analysis is being carried out. The
cutting point is generally chosen such that it includes
the process whose verification is sought (Larocque et al.
1998). We note that correlation does not always imply
cause as the correlation measures only the degree of theshared increase or decrease of two variables. However,both series (xt and yt ) may be varying as a function of athird variable. For our study, we assumed that xt is thecause of yt due to previous knowledge of how the system
works.
The cross-correlation function generates a graph (thecorrelogram), which shows the correlation (positive or
negative) between the input and output elements for eachtime lag k . The delay between the input and output
series is defined as the time interval between k = 0 andthe maximum value of the correlogram. The shorter thisinterval, the more rapidly the output series yt reacts to asignal of the input series xt . If there is no time lag between
the input and output, the influence of the former on
the latter is instantaneous, and the maximum correlation
is found for k = 0. In this study, the correlogram wascalculated using the time series of rainfall data as input
and the water table data as output (Figure 1). A cutting
point value of 15.552.000 s (180 d) was chosen as only
one major recharge peak occurs each year, a time window
large enough to cover the recharge process. The two timeseries were used to calculate the mean time duration
D.C. Neto et al. Groundwater 54, no. 1: 82–91

85

The difference between the maximum and minimum
values of X is the range R, which for a reservoir is thedifference between the maximum and minimum levelsof water stored over a period of τ years. This range isexpressed as:
R (τ ) = max X (t, τ ) − min X (t, τ ) , (1 ≤ t ≤ τ ) (10)while the standard deviation, S , is given by


τ
 1

2
ξ (t) − ξ τ
S=
τ

(11)

t=1

Figure 5. Correlogram resulting from the analysis of watertable and rainfall data.

t for the rainfall events to move through the vadosezone to the water table. The highest correlation value of
0.278 was obtained for a time period t of 7,516,800 s,
or 87 d (Figure 5). This is hence an estimate of theaverage time for rain water to percolate from the soil
surface to the phreatic surface. We will later test this delay
between rainfall and water table response by carryingout 1D transient simulations of water flow through thevadose zone.
Rescaled Range (R/S) Analysis
Another method known as rescale range analysis,or R/S analysis (Hurst et al. 1965), has been widelyused to analyze time series of natural phenomena,
including soil, hydrologic and meteorological processes(e.g., Mandelbrot and Wallis 1969; Christofoletti 1997;
Breslin and Belwards 1999; Miranda and Andrade 2001;
Peters et al. 2002; Favaretto 2004; Green et al. 2008). The
analysis generates a value known as the Hurst exponent(H ), which is able to distinguish between random andnon-random time series, as well as provides information
about the degree of autocorrelation of the series being
considered. The R/S analysis is best described within
the context of modeling a reservoir, for which it was
originally created. If one considers that the input of water
into a reservoir in year u is represented by ξ (u), thedeviation from the mean ξ (u) − ξ τ for a specific yearτ is calculated on the basis of the average over several
years as follows:
1
ξ (u)
τ
τ

ξ τ =(8)

u=1

The sum of the deviations, X , of the discharge fromthe mean is thus given by:
X (t, τ ) =tξ (u) − ξ τ



u=1

86

D.C. Neto et al. Groundwater 54, no. 1: 82–91(9)

Dividing R by S for various intervals of the timeseries data and calculating the mean (or expected value,E ), Hurst et al. (1965) observed that the values converged
to an exponential relationship:E

R (N )= kN HS (N )

with N → ∞(12)where k is a constant, N is the number of observations,and H is the Hurst exponent.
If H = 0.5, deviations of the values in the time series
from the mean are random and hence independent of thoseof previous and subsequent elements. If 0.5 < H < 1, theseries is persistent, in which case deviations tend toremain stable over time (e.g., 1 year of flooding tends tobe followed by another year of flooding, or a dry year isfollowed by another dry year). The probability of a value
in the series deviating in the same direction (either positiveor negative) as the previous one increases as H increases.
On the other hand, if 0 < H < 0.5, deviations will be antipersistent in that a positive deviation then tends to befollowed by a negative deviation, and vice versa. The
probability of this inversion of the deviation from themean to occur increases as the value of H decreases.
Breslin and Belwards (1999) developed an autocorrelation average, C , related to the Hurst exponent, using
the equation:
C = 22H −1 − 1(13)

If H = 0.5, then C = 0, which implies that for arandom series no autocorrelation exists between themeasurements, and succeeding values are independentof each other. If H > 0.5 or H < 0.5, however, thecorrelation is either positive or negative, meaning thatthe system “remembers” previous events, and hence thatpast events will affect future events. The value of a given
measurement is then affected by previous measurements,
and the current measurement will also affect future results.
The H exponent can be estimated by first establishing
the size of the largest and smallest intervals of thetime series being considered (the various N values inEquation 12). Next the value of R/S is calculated forthe largest interval (a value equal or close to that of
the complete series). This series is then divided into two
NGWA.org

Table 1
Results for Estimating the Hurst Exponent of the
Measured Water Levels
Number of
Readings
8
16
32
64
128
256
512
1024
2048
4096
8192
16,38432,768Mean(RS)

Log 2 of
Number of Readings

Log 2 of
Mean RS

3.491
6.941
13.86127.70855.426110.458
218.893435.093905.375
1719.366
3533.701
5219.362
10118.010

3
4
5
6
7
8
9
10
11
12
13
14
15

1.804
2.795
3.793
4.792
5.793
6.787
7.774
8.765
9.822
10.748
11.78712.35013.305

parts, and the R/S value of each of them is calculated,thus determining the mean. Each of these halves issubsequently halved again, with the process of calculating
the mean of the two halves being repeated until the valueof the original smallest interval is reached.
Converting both sides of Equation 12 to logarithms
gives
log

R (N)
S (N )

= log kN H = H log N + log k

Table 2
Results for Estimating the Hurst Exponent of the
Rainfall Data

(14)

Equation 14 shows that a log–log graph of thecalculated values, or a linear graph of the logarithms ofthe values, will make it possible to estimate the value of
H from the slope of the straight line through simple linear
regression.
The aforementioned procedure for the H exponent
was coded in a computer program and applied to both
the water table data and the rainfall measurements. For
the groundwater table, we used the same time series as
before, but with temporal spacings of 1 h between readings. The shortest interval involved 8 readings and thelongest interval 32,768 readings, the latter number being
the power of 2 closest to the total number of readings(slightly more than 34,000). The points obtained in thismanner are presented in Table 1. A graph of these points,
along with the fitted line, is shown in Figure 6. The valueof the Hurst exponent, which corresponds to the slope ofthe straight line, was 0.971. As this value is greater than0.5, it shows that variations in the groundwater table are
persistent in that the system has a memory of past levels,with strongly autocorrelated values. In comparison,
the value of H obtained by Hurst from an analysis ofminimum water levels of the Nile River from the year
622 to 1469 was 0.91 (Mandelbrot and Wallis 1969).
Using the same process for the rainfall data produced
the results in Table 2 and Figure 6. The H exponent in thiscase was 0.766, again indicating a persistent phenomenon
NGWA.org

Figure 6. Graph of points in Tables 1 and 2 (with marks),
together with the fitted lines, obtained with the R/S analysis
for the water table data (a) and rainfall data (b).

Number of
Readings

Mean
(RS)

Log 2 of
Number of
Readings

Log 2 of
Mean RS

8
16
32
64
128
256
512
1024
2048
4096
8192
16,38432,7684.254
8.868
14.74027.66654.24092.420140.618200.090379.913759.1971463.844
1748.165
2085.897

3
4
5
6
7
8
9
10
11
12
13
14
15

2.089
3.149
3.882
4.790
5.761
6.530
7.136
7.645
8.570
9.568
10.51610.772
11.026and consistent with results obtained for other rainfall timeseries (e.g., Favaretto 2004; Peters and Christensen 2006).
The H exponent for the water table data was found to
be larger than the value for rainfall. This can be explainedby the fact that the aquifer functions as a reservoir which,to a certain extent, integrates the rainfall variations overtime and space. Therefore, for the various time scales, theeffect of occasional dry periods interspersed by periodsof rainfall are somewhat buffered by the aquifer as theinfiltrating rain water is stored and released only slowly,
and then also only in part to the aquifer. This leads toa greater autocorrelation of the water table values andD.C. Neto et al. Groundwater 54, no. 1: 82–91

87suggests that the aquifer system has a relatively long“memory” of the water table fluctuations, thereby showing
the “Joseph effect” of the series. This term, taken from thebiblical story of the Prophet Joseph foreseeing 7 years ofplenty followed by 7 years of famine in Egypt, was usedby Mandelbrot and Wallis (1968) to describe the effect of
persistence in hydrological time series with H > 0.5, in
which an increase or decrease in values tends to follows
the same direction as previous values.

Recharge Calculations
Further insight into recharge and the response of thewater table was obtained by simulating infiltration andvariably saturated flow in the soil profile during the 4-year
experiment. Simulations were carried out using version
4.08 of the HYDRUS-1D software package (http://www.
pc-progress.com/en/Default.aspx?hydrus-1d) of Šimůnek
et al. (2008). This code has been used in several recent
studies to estimate recharge processes (e.g., Lu et al. 2011;
Assefa and Woodbury 2013). The calculations assumed
applicability of the Richards equation to 1D verticalvariably saturated flow subject to root water uptake by thegrass cover at the site. We refer to the HYDRUS technicalmanual (Šimůnek et al. 2008) for a detailed description of
the governing equations, as well as of various initial andboundary conditions that can be implemented.
For the recharge calculations, we assumed a verticalprofile of 10 m containing an observation node at a
depth of 9.5 m, being the approximate location of thepressure transducer. Atmospheric boundary conditions
were assigned to the soil surface assuming runoff to occur
when surface ponding exceeded 2 cm. Daily values ofprecipitation and potential evapotranspiration were used,in combination with the water stress response model ofFeddes et al. (1978) to account for root water uptake
by the grass cover. Daily potential evapotranspiration
rates were calculated using the approach of Hargreaves(Hargreaves 1975; Jensen et al. 1997) based on available
temperature data, the latitude (22.24 S) and altitude(643 m) of the site, the leaf area index of the grasscover (estimated to be 2.0) and the rooting depth (70 cm).
Default parameter values of the water stress response
functions were used based on literature values for grass(Feddes et al. 1978) as implemented in HYDRUS-1D. A
linear decreasing root distribution from the soil surface(maximum) to a depth of 70 cm was assumed.
The van Genuchten–Mualem equations (vanGenuchten 1980) were used to describe the unsaturated
soil hydraulic properties. Soil hydraulic parameters(Table 3) were estimated from available soil texture datausing the pedotransfer functions of Schaap et al. (1998)
as implemented in HYDRUS-1D, but with measured
values of the saturated hydraulic conductivity, K s (Alfaro
Soto et al. 2007).
The initial condition was taken to be simply−200 cm in the upper part of the profile and, in a second
simulation, an equilibrium distribution in the lower partof the profile with the water table at 800 cm. Only an
88

D.C. Neto et al. Groundwater 54, no. 1: 82–91

approximate estimate of the initial condition was needed
as we started the numerical calculations using local
weather data on January 1, 2001 (taken to be t = 0 here),
1 year before the 4-year water table study period covering
the period January 1, 2002 through December 31, 2005.
This ensures that the adopted initial moisture profile att = 365 d was as consistent as possible with the local soil
and hydrologic conditions.
Several alternatives were considered for the lower
boundary condition. Figure 7 shows the calculated
recharge rate for a free-draining profile that neglects thepresence of a water table. The plot indicates excellent
visual correlation with the water table series in Figure 1,
with the highest recharge rates occurring during periods
when groundwater levels increased the fastest (the steepest
rising parts of the graph in Figure 7). An interestingfinding was that the average percolation time of the
moisture front during the rainy season was always very
close to the transit time of 87 d derived using the crosscorrelation analysis. This is best illustrated using the
precipitation data of year 2005, and the recharge data in
Figure 7 of that same year. The rainfall data showed avery high peak at t = 1491 d (the very end of January2005) from the start of the simulation (more than 6 and12 cm on days 1490 and 1491, respectively). The peak of
the recharge plot that year (2005) occurred on day 1580,
which is 89 d after the large rainfall event. This result isvery close indeed to the average transit time of 87 d foundusing the cross-correlation analysis.
Some caution about the aforementioned numericalresults is warranted as the recharge rate depends also on
the transients of the water table and the upward/downward
movement of the capillary fringe during the year. One
alternative boundary condition would be to use a transientcondition specifying the water table vs. time consistent
with the water levels shown in Figure 1, or perhaps evena constant water table. The latter case was considered by
fixing the water table at a depth of 800 cm (150 cm above
the monitoring point). The calculated recharge rates forthis case remained very close to those in Figure 7, exceptthat the peaks occurred slightly earlier. In particular, thepeak recharge rate in 2005 occurred 65 d after the intense
rainfall event at t = 1490 d.
We further implemented a special deep drainage
lower boundary condition to account in an approximate
way for regional flow effects on the recharge rate.
Use of this boundary condition was motivated by anexperimental study by Ernst and Feddes (1979) who
established a relationship between the recharge rate
and the groundwater level at a particular location (seealso Hopmans and Stricker 1989). The relationship was
generalized within the HYDRUS-1D software (Šimůnek
et al. (2008)) to the more general boundary conditionqL (t) = A exp −B hL − hgw (15)

where q L is the imposed lower boundary flux, h L is thetransient pressure head (or water table elevation) at thelower boundary, A and B are adjustable parameters, andNGWA.org

Table 3
van Genuchten–Mualem Soil Hydraulic Parameters Used for the Recharge Calculations (θ r Is the Residual
Water Content, θ s Is the Saturated Water Content, α, n, and l Are Empirical Shape Factors, and K s Is theSaturated Hydraulic Conductivity)
Soil Depth (cm)0–150
150–250
250–350
350–450
450–550550–1000θ r (cm3 /cm3 )θ s (cm3 /cm3 )α (cm−1 )

n (−)

l (−)

K s (cm/d)

0.039
0.043
0.046
0.050
0.056
0.057

0.387
0.386
0.385
0.387
0.386
0.386

0.03340.03110.02940.02570.0266
0.0263

1.42
1.40
1.39
1.39
1.37
1.36

0.50
0.50
0.50
0.50
0.50
0.50

53.0
276
378
436
49.8
99.1

Figure 8. Observed (solid line) and calculated (dashed line)
water table locations obtained for the time period 2002
(t = 365 d) through 2005 (t = 1825 d), assuming a 10-m deepdraining soil profile.
Figure 7. Calculated recharge rates at 900 cm depth forthe time period 2002 (t = 365 d) through 2005 (t = 1825 d),assuming a 10-m free-draining soil profile. The highest
recharge rates occurred at approximately 545, 820, 1230, and1550 d.

h gw is some long-term equilibrium water table position
relative to the lower boundary. In this study, we considered the three parameters (A, B , and h gw ) to be completely
empirical and fitted them to the available water table data.
Figure 8 shows the calculated water table position interms of simulated pressure heads at the monitoring point,
similarly as the measurements in Figure 1. Calculations
used values of −1.8 cm/d for A (the negative sign
indicating a downward flux), −0.007/cm for B , and800 cm for h gw . Calculated water levels in Figure 8 wererelatively close to the observed values in 2004 and 2005,
but underpredicted the water table fluctuations in 2002
and early 2003 (between approximately 365 and 760 d).
Figure 9 shows calculated recharge rates for the deep
drainage boundary condition. As compared with Figure 7,
the recharge rate in Figure 9 is much less variable intime as the net downward flux now is affected by the
upward and downward movement of the capillary fringe.
However, the overall (cumulative) recharge rates during
the 4-year simulation (2002 to 2005) for the two boundary
conditions were relatively close: 137 vs. 121 cm during
4 years (34.3 vs. 30.4 cm/year). Our results indicate thatpredictions of the yearly averaged recharge rate using
root zone modeling using a free-drainage draining lower
boundary condition may give similar results as the more
NGWA.org

Figure 9. Calculated recharge rates at 900 cm depth forthe time period 2002 (t = 365 d) through 2005 (t = 1825 d),assuming a 10-m deep-draining soil profile. The highest
recharge rates occurred at approximately 545, 820, 1230, and1550 d.

realistic deep drainage condition, but not necessarily of
the transient recharge distribution during the year.

Concluding Remarks
Systematic monitoring of the water column made it
possible to obtain a detailed time series involving bothvery small oscillations in amplitude (on the order of
10−3 m) related to natural phenomena such as earth tides
D.C. Neto et al. Groundwater 54, no. 1: 82–91

89and earthquakes, as well as much larger amplitudes (ofapproximately 2 to 3 m) related to seasonal cycles ofrecharge and discharge of the aquifer. Statistical analysesof the data using Fourier analysis, cross-correlation, andR/S analysis (the Hurst exponent) corroborated qualitative
observations, such as the frequency of earth tides (twotimes per day) and seasonal oscillations (one cycle peryear). Oscillations in the water table showed that the fieldsystem has a definite “memory” of previous events, thus
very much reflecting the fact that the effects of futurehydrologic events will be influenced by past results.
Our analysis also identified temporal relationships
between rainfall and aquifer recharge, as well as theaverage time required for seasonal infiltration fronts toreach the aquifer (87 d in this study). The calculated
average percolation time was very close to estimatesobtained using numerical solutions of the Richards
equation for flow in variably saturated vadose zone
systems. Implementation of a special deep drainage
boundary condition in the numerical solution allowed us
to obtain independent estimates of the seasonable watertable fluctuations, with results being very consistent
with observed data. Numerical results indicated thatpredictions of the yearly averaged recharge rate using
root zone modeling using a free-drainage draining lower
boundary condition may give similar results as the more
realistic deep drainage condition, but not necessarily of
the transient recharge distribution during the year.

Acknowledgments
We acknowledge support for this study from theBasin Studies Laboratory (LEBAC) of the Department
of Applied Geology, associated with the Center for
Environmental Studies of UNESP, the Foundation for
Development of UNESP (FUNDUNESP), the NationalCounsel for Technological and Scientific Development(CNPq) (350983/1997-6), and CAPES-Brazil. We thank
Jack Guswa and two anonymous reviewers for theirhelpful comments and suggestions.

References
Alfaro Soto, M.A., D.M. Kumayama, and H.K. Chang. 2007.
Calibration of a reflectometer for water flow studies inunsaturated soils. Geociências 26, no. 4: 357–368 (inPortuguese).
Assefa, K.A., and A.D. Woodbury. 2013. Transient, spatiallyvaried groundwater recharge modeling. Water Resources
Research 49, no. 2: 4593–4606.
Bracewell, R.N. 1986. The Fourier Transform and its Applications. New York: McGraw-Hill.
Bredehoeft, J.D. 1967. Response of well-aquifer systems toearth tides. Journal of Geophysical Research 72, no. 12:
3075–3087.
Breslin, M.C., and J.A. Belwards. 1999. Fractal dimensions
for rainfall time series. Mathematics and Computers in
Simulation 48, no. 4–6: 437–446.
Chen, C.-H., C.-H. Wang, S. Wen, T.-K. Yeh, C.-H. Lin, J.-Y.
Liu, H.-Y. Yen, C. Lin, R.-J. Rau, and T.-W. Lin. 2013.
Anomalous frequency characteristics of groundwater level

90

D.C. Neto et al. Groundwater 54, no. 1: 82–91

before major earthquakes in Taiwan. Hydrology and Earth
System Sciences 17, no. 5: 1693–1703.
Christofoletti, A.L.H. 1997. Fractal and multifractal analysis of
structure of rainy seasons in São Paulo State. Doctoral
dissertation, Instituto de Geociências e Ciências Exatas,
Universidade Estadual Paulista, Rio Claro, SP, Brazil,
263 pp (in Portuguese).
Cutillo, P.A., and J.D. Bredehoeft. 2011. Estimating aquifer
properties from water level response to earth tides. Ground
Water 49, no. 4: 600–610.
DAEE. 1981. Estudo de águas subterrâneas, Região Administrativa 5, Campinas. Departamento de Águas e Energia
Elétrica, São Paulo, SP, Brazil (in Portuguese).
Ernst, L.F., and R.A. Feddes. 1979. Invloed van grondwateronttrekking voor beregening en drinkwater op de grondwaterstand. Nota 116, Instituut voor Cultuurtechniek en
Waterhuishouding, ICW, Wageningen, the Netherlands (inDutch).
Favaretto, A.B. 2004. Estimate of the Hurst exponent of rainfall
time-series in São Paulo State using Fourier transforms,
wavelets and R/S analysis. Master dissertation, Institutode Geociências e Ciências Exatas, Universidade EstadualPaulista, Rio Claro, SP, Brazil, 89 pp (in Portuguese).
Feddes, R.A., P.J. Kowalik, and H. Zaradny. 1978. Simulationof Field Water Use and Crop Yield . Simulation Monograph.
Wageningen, the Netherlands: Pudoc.
Fleming, S.W., A.M. Lavenue, A.H. Aly, and A. Adams. 2002.
Practical applications of spectral analysis to hydrologic timeseries. Hydrological Processes 16, no. 2: 565–574.
Green, T.R., G.H. Dunn, R.H. Erskine, J.D. Salas, and L.R.
Ahuja. 2008. Fractal analyses of steady infiltration andterrain on an undulating agricultural field. Vadose Zone
Journal 8, no. 2: 310–320.
Hargreaves, G.H. 1975. Moisture availability and crop production. Transactions of the American Society of AgriculturalEngineers 18, no. 5: 980–984.
Healy, R.W., and P.G. Cook. 2002. Using groundwater levels
to estimate recharge. Hydrogeology Journal 10, no. 1:
91–109.
Hopmans, J.W., and J.N.M. Stricker. 1989. Stochastic analysis
of soil water regime in a watershed. Journal of Hydrology
105: 57–84.
Hurst, H.E., R.P. Black, and Y.M. Simaika. 1965. Long-Term
Storage: An Experimental Study. London, UK: Constable& Co.
Iglauer, S., and A. Muggeridge. 2013. The impact of tides on
the capillary transition zone. Transport in Porous Media 97,
no. 1: 87–103.
Jensen, D.T., G.H. Hargreaves, B. Temesgen, and R.G. Allen.
1997. Computation of ETo under nonideal conditions.
Journal of the Irrigation and Drainage Division, ASCE 123,
no. 5: 394–400.
Larocque, M., A. Mangin, M. Razack, and O. Banton. 1998.
Contribution of correlation and spectral analyses to theregional study of a large Karst aquifer (Charente, France).
Journal of Hydrology 205, no. 3–4: 217–231.
Lu, X., M. Jin, M.T. van Genuchten, and B. Wang. 2011.
Ground water recharge at five representative sites in theHebei Plain of China: case study. Ground Water 49, no. 2:
286–294.
Mandelbrot, B.B., and J.R. Wallis. 1968. Noah, Joseph andoperational hydrology. Water Resources Research 4, no. 5:
909–918.
Mandelbrot, B.B., and J.R. Wallis. 1969. Some long-run properties of geophysical records. Water Resources Research 5,
no. 2: 321–340.
Miranda, J.G.V., and R.F.S. Andrade. 2001. R/S analysis of pluviometric records: comparison with numerical experiments.
Physica A: Statistical Mechanics and its Applications 295,
no. 1: 38–41.

NGWA.org

Neeper, D.A. 2001. A model of oscillatory transport in
granular soils, with application to barometric pumping
and earth tides. Contaminant Hydrology 48, no. 3–4:
237–252.
Oliva, A., H.K. Chang, and M.R. Caetano-Chang. 2005. Determination of hydraulic conductivity of the Rio Claro Formation: comparative analysis using granulometric analysis andtests with the Guelph permeameter. Águas Subterrâneas 19,
no. 1: 1–18 (in Portuguese).
Padilla, A., and A. Pulido-Bosch. 1995. Study of hydrographs of
Karstic aquifers by means of correlation and cross-spectral
analysis. Journal of Hydrology 168, no. 1: 73–89.
Peters, O., and K. Christensen. 2006. Rain viewed as relaxational
events. Journal of Hydrology 328, no. 1–2: 46–55.
Peters, O., C. Hertlein, and K. Christensen. 2002. A complexity
view of rainfall. Physical Review Letters 88, no. 1:
018701–1.
Pugh, D.T. 1987. Tides, Surges, and Mean Sea-Level: A
Handbook for Engineers and Scientists. Chichester, UK:
John Wiley.
Robinson, T.W. 1939. Earth-tides shown by fluctuations ofwater-levels in wells in New Mexico and Iowa. Transactions of the American Geophysical Union 20, no. 4:
656–666.
Roeloffs, E. 1988. Hydrologic precursors to earthquakes: areview. Pure and Applied Geophysics 126, no. 2–4:
177–209.
Rojstaczer, S., and D.C. Agnew. 1989. The influence of formation material properties on the response of water levels in

NGWA.org

wells to Earth tides and atmospheric loading. Journal ofGeophysical Research B 94, no. B9: 12403–12411.
Schaap, M.G., F.J. Leij, and M.T. van Genuchten. 1998.
Neural network analysis for hierarchical prediction ofsoil hydraulic properties. Soil Science Society of America
Journal 62, no. 4: 847–855.
Si, B.C. 2007. Spatial scaling analyses of soil physicalproperties: a review of spectral and wavelet methods.
Vadose Zone Journal 7, no. 2: 547–562.
Šimůnek, J., M. Šejna, H. Saito, M. Sakai, and M.Th. van
Genuchten. 2008. The HYDRUS-1D software package forsimulating the one-dimensional movement of water, heatand multiple solutes in variably-saturated media. Version
4.08, HYDRUS Software Series 3, Department of Environmental Sciences, University of California, Riverside, CA.
Stafford, B.P., and W.G. Rixey. 2011. Distribution of fuelgrade ethanol near a dynamic water table. Ground Water
Monitoring & Remediation 31, no. 3: 55–60.
Todd, D.K. 1959. Groundwater Hydrology. New York: John
Wiley and Sons.
Van Camp, M., and O. Vauterin. 2005. Tsoft: graphical andinteractive software for the analysis of time series and earth
tides. Computers & Geosciences 31, no. 5: 631–640.
van Genuchten, M.T. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil
Science Society of America Journal 44, no. 5: 892–898.
Zeumann, S., A. Weise, and T. Jahr. 2009. Tidal and nontidal signals in groundwater boreholes in the KTB area,
Germany. Journal of Geodynamics 48, no. 3–5: 115–119.

D.C. Neto et al. Groundwater 54, no. 1: 82–91

91