Analyses of uncertainties and scaling of groundwater level fluctuations

Analytical solutions were derived for the variance, covariance, and spectrum of groundwater level, h(x, t), in an unconfined aquifer described by a linearized Boussinesq equation, with random source/sink and initial and boundary conditions. It was found that in a typical aquifer, the error in h(x, t) at an early point in time is mainly caused by the random initial condition, and the error reduces as time progresses to reach a constant error at a later time. The duration for which the effect of the random initial condition is significant may be a few hundred days in most aquifers. The constant error in h(x, t) at a later time is due to the combined effects of the uncertainties in the source/sink and flux boundary: the closer to the flux boundary, the larger the error. The error caused by the uncertain head boundary is limited to a narrow zone near the boundary and remains more or less constant over time. The aquifer system behaves as a low-pass filter which filters out high-frequency noises and keeps lowfrequency variations. Temporal scaling of groundwater level fluctuations exists in most parts of a low permeable aquifer whose horizontal length is much larger than its thickness, caused by the temporal fluctuations of areal source/sink.


Introduction
Groundwater level or hydraulic head (h) is the main driving force for water flow and advective contaminant transport in aquifers and thus the most important variable studied in groundwater hydrology and its applications.Knowledge about h is critical in dealing with groundwater-related environmental problems, such as over-pumping, subsidence, sea water intrusion, and contamination.One often finds that data about groundwater level are limited or unavailable in a hydrogeological investigation.In such cases the groundwater level distribution and its temporal variation are usually obtained with an analytical or numerical solution for a groundwater flow model.
It is obvious that errors always exist in the groundwater levels calculated or simulated with analytical or numerical solutions.The main sources of errors include the simplification or approximation in a conceptual model and uncertainties in the model parameters.Problems in conceptualization or model structure have been dealt with by many researchers (Neuman, 2003;Rojas et al., 2008Rojas et al., , 2010;;Ye et al., 2008;Refsgaard et al., 2007;Zeng et al., 2013).Uncertainties in the model parameters (e.g., hydraulic conductivity, recharge rate, evapotranspiration, and river conductance) have been investigated, based on generalized likelihood uncertainty estimation and Bayesian methods (Nowak et al., 2010;Neuman et al., 2012;Rojas et al., 2008Rojas et al., , 2010)).The uncertainty in groundwater level has been one of the main research topics in stochastic subsurface hydrology for more than 3 decades.Most of these studies were focused on the spatial variability of groundwater level due to aquifers' heterogeneity (Dagan, 1989;Gelhar, 1993;Zhang, 2002).Little attention has been given to the uncertainties in groundwater level due to temporal variations in hydrological processes, e.g., recharge, evapotranspiration, discharge to a river, and river stage (Bloomfield and Little, 2010; Zhang and Schilling, 2004;Schilling and Zhang, 2012;Liang and Zhang, 2013a;Zhu et al., 2012).
Uncertainties in groundwater level fluctuations have been studied by Zhang andLi (2005, 2006) and most recently X. Y. Liang and Y.-K.Zhang: Analyses of uncertainties and scaling of groundwater level fluctuations by Liang and Zhang (2013a).Based on a linear reservoir model with a white noise or temporally correlated recharge process, Zhang andLi (2005, 2006) derived the variance and covariance of h(t) by considering only a random source or sink process, assuming deterministic initial and boundary conditions.Liang and Zhang (2013a) extended the studies of Zhang andLi (2005, 2006) and carried out nonstationary spectral analysis and Monte Carlo simulations using a linearized Boussinesq equation, and investigated the temporospatial variations in groundwater level.However, the only random process considered by Liang and Zhang (2013a) is the source/sink.Temporal scaling of groundwater levels, discovered first by Zhang and Schilling (2004), was verified in several studies (Zhang andLi, 2005, 2006;Bloomfield and Little, 2010;Zhang and Yang, 2010;Zhu et al., 2012;Schilling and Zhang, 2012).However, we do not know the effect of random boundary conditions on temporal scaling of groundwater levels.
In this study we extended the above-mentioned work by considering the groundwater flow in a bounded aquifer, described by a linearized Boussinesq equation, with a random source/sink as well as random initial and boundary conditions, since the latter processes are known to give uncertainties.The objectives of this study are (1) to derive analytical solutions for the covariance, variance, and spectrum of groundwater level, and (2) to investigate the individual and combined effects of these random processes on uncertainties and scaling of h(x, t).In the following, we will first present the formulation and analytical solutions, then discuss the results, and finally draw some conclusions.

Formulation and solutions
Under the Dupuit assumption, the one-dimensional transient groundwater flow in an unconfined aquifer near a river (Fig. 1) can be approximated with the linearized Boussinesq equation (Bear, 1972) with the initial and boundary conditions, i.e., where T [L T −1 ] is the transmissivity, h [L] is the hydraulic head or groundwater level above the bottom of the aquifer, which is assumed to be horizontal, W (t) [L T −1 ] is the timedependent source/sink term, representing areal recharge or evapotranspiration, S Y is the specific yield,  spatially random variable.The source/sink, W (t), the flux to the left boundary, Q(t), and the head at the right boundary, H (t), are all taken to be temporally random processes and spatially deterministic.The parameters T and S Y are taken to be constant.The groundwater level, h(x, t), the three random processes, W (t), Q(t), and H (t), and the random variable, H 0 (x), are expressed in terms of their respective ensemble means plus small perturbations, where stands for ensemble average and for perturbation.The initial condition H 0 (x) in Eq. (1) can be any function.For the conceptualization of the groundwater flow presented in Fig. 1, the steady-state condition can be reached in this aquifer after a rainfall or during a wet season.Thus the steady-state solution to this model was often adopted as the initial condition in previous research (Liang andZhang, 2012, 2013a, b).Thus, in this study, we set initial condition H 0 (x) to be the steady-state solution to the one-dimensional groundwater flow equation, i.e., H 0 (x) = h 0 + 0.5W 0 (L 2 − x 2 )/T , whereh 0 [L] is the constant groundwater level at the right boundary and W 0 [L T −1 ] is the spatially constant recharge rate (Liang and Zhang, 2012).Since h 0 is taken to be constant, the source of the uncertainty in the initial head H 0 (x) is due to random W 0 only.Thus, the mean and perturbation of H 0 (x) can be written as H 0 (x) = h 0 + 0.5 W 0 (x) (L 2 − x 2 )/T and H 0 (x) = 0.5W 0 (L 2 − x 2 )/T , respectively.By substituting Eq. ( 2), H 0 (x) , and H 0 (x) into Eq.( 1 Subtracting Eq. (3) from Eq. (1) leads to the following perturbation equation with the initial and boundary conditions The analytical solution to Eq. ( 4) can be derived with integral-transform methods (Özisik, 1968) given by where β = T /S Y , b n = (2n + 1)π /2L.Using Eq. ( 5), the temporal covariance of the groundwater level fluctuations can be derived as , and H (t), respectively.We assume that W (t), Q(t), and H (t) are uncorrelated in order to simplify our analyses.It is shown in Eq. ( 6) that the head covariance depends on the variance of W 0 and the covariances of W (t), Q(t), and H (t) and this equation can be evaluated for any random W (t), Q(t), and H (t). We assume that these processes are white noise, as employed in previous studies (Gelhar, 1993;Hantush and Marino, 1994;Liang and Zhang, 2013a).More realistic randomness of these processes will be considered in future studies.
Following Gelhar (1993, p. 34), we express the spectra of W (t), Q(t), and and S H H = σ 2 H λ H /π , respectively, where σ 2 W , σ 2 Q , and σ 2 H are the variances and λ W , λ Q , and λ H are the correlation time intervals of these three processes, respectively.The corresponding covariances of W (t), Q(t) and . Substituting these covariances into Eq.( 6) and taking integration, one obtains an analytical solution of head covariance where τ = t 2 − t 1 and t = (t 2 + t 1 )/2.The analytical solution for the head variance can be obtain by setting τ = 0 where (Gelhar, 1993) where the transmissivity (T ) is replaced by the product of the hydraulic conductivity (K) and the average saturated thickness (M) of the aquifer.The characteristic timescale (t c ) is an important parameter and its value for most shallow aquifers is usually larger than 100 days, since the horizontal extent of a shallow aquifer is usually much larger than its thickness.For instance, the value of t c is 250 days for a sandy aquifer with L = 100 m, M = 10 m, K = 1 m day −1 , and S Y = 0.25.The spectral density of h(x, t) cannot be derived by ordinary Fourier transform, since the head covariance and variance depend on time t , and thus h(x, t) are temporally non-stationary as shown in Eqs. ( 7) and (8).Priestley (1981) defined the spectral density of non-stationary processes (Wigner spectrum) as the Fourier transform of time-dependent auto-covariance with fixed reference time t and derived time-dependent spectral density.In order to obtain the spectrum of h(x, t), we applied Priestley's method and obtained the time-dependent spectral density (Priestley, 1981;Zhang and Li, 2005;Liang and Zhang, 2013a), i.e., where ω is angular frequency and ω = 2π f , f is frequency, and i = √ −1.It is seen in Eq. ( 9) that the spectrum S hh is dependent on not only frequency and locations but also time t.The time-dependent term (i.e., first term) in Eq. ( 9) is caused by the random initial condition and is proportional to e −β(b 2 m +b 2 m )t which decays quickly with t.We evaluated the first term in the Eq. ( 9) by setting t = 0 and found that it is much smaller than the second term in Eq. ( 9).We thus ignored the first term and evaluated the spectrum using the approximation, 3 Results and discussion

Variance of groundwater levels
The general expression of the head variance in Eq. ( 8) depends on the variances of the four random processes, Q , and σ 2 H .In the following, we will study their individual and combined effects on the head variation and focus our attention only on the variance of h(x, t).The dimensionless standard deviation of h(x, t), σ h , and the square root of the dimensionless variance, (σ 2 h ), as a function of the dimensionless time (t ), are evaluated and presented in the left column of Fig. 2 at fixed dimensionless locations (x).The σ h as a function of x was evaluated and is presented in the right column of Fig. 2 at fixed t .
We first evaluate the effect of the random initial condition due to the random term, W 0 , by setting σ 2 W = σ 2 Q = σ 2 H = 0.In this case, the dimensionless variance in Eq. ( 8) reduces to    12) where σ 2 and (f) are based on Eq. ( 13) where σ 2 (g) and (h) are based on Eq. ( 14) where σ 2 and (i) and (j) are based on Eq. ( 15) where σ 2 where σ 2 h = σ 2 h T 2 /(4L 4 σ 2 W 0 ).The changes of the σ h with x and t are presented in Fig. 2a and b 2a that for a fixed location the σ h is at its maximum at t = 0 and it decreases gradually with time to a negligible number at t = 1.0.This means that the error in h(x, t), predicted by an analytical or numerical solution due to the uncertain initial condition, is significant at an early point in time, especially near a flux boundary.The duration for which the effect of the uncertain initial condition is significant depends on the value of the characteristic timescale (t c ), since t = t/t c .In most aquifers this duration may be many days.
In the typical aquifer studied, the effect of the uncertainty in the initial condition on h(x, t) is significant during the first 250 days (t = 1.0).This duration should be relatively short, however, in a more permeable aquifer whose horizontal extent (L) is relatively smaller than its thickness (M).It is seen in Fig. 2b that for a fixed time the σ h is the largest at the left flux boundary (x = 0.0) and becomes zero at the right constant head boundary (x = 1.0), since the right boundary is deterministic.This means that the error in h(x, t) predicted by an analytical or numerical solution due to the uncertain initial condition is significant almost everywhere in the aquifer: the further away from a constant head boundary, the larger the error.
We then consider the uncertainty in the areal source/sink term (W ) by setting σ 2 In this case the dimensionless variance in Eq. ( 8) reduces to where σ 2 h = σ 2 h T S Y /(4L 2 σ 2 w λ W ). The changes of the σ h with x and t are presented in Fig. 2c and d, respectively.It is noticed in Fig. 2c that at a fixed location, the σ h is zero initially, gradually increases as time goes, and approaches a constant limit at later time.This means that the error in h(x, t) due to an source/sink is at its minimum at early time and increases with time to approach a constant limit at later time: the closer to the left flux boundary, the larger the limit.For a fixed time the σ h decreases smoothly from the left to the right boundary (Fig. 2d).The error in h(x, t) due to the uncertainty in the source/sink is significant almost everywhere in the aquifer: the further away from the constant head boundary, the larger the error, similar to the previous case with the random initial condition (Fig. 2b).
Thirdly, we investigate the effect of the left random flux boundary by setting σ 2 W 0 = σ 2 W = σ 2 H = 0 in Eq. ( 8).In this case the dimensionless head variance is given by where The changes of the σ h with x and t are presented in Fig. 2e and f, respectively.At any location, the σ h in Fig. 2e or the error in h(x, t), due to an uncertain flux boundary, is at its minimum at an early point in time, and it increases quickly with time to approach a constant limit: the closer to the left flux boundary, the larger the limit.At any time, the σ h in Fig. 2f or the error in the water head due to the uncertain flux boundary is at its maximum at the left boundary but decreases quickly away from the boundary to become insignificant for x > 0.8.
Fourthly, we investigated the effect of the random head boundary by setting σ 2 W 0 = σ 2 W = σ 2 Q = 0 in Eq. ( 8).The dimensionless head variance in this case is given by where The changes of this σ h with x and t are presented in Fig. 2g and h, respectively.It is seen in Fig. 2g that at any location, the σ h or the error in h(x, t), due to the random head boundary, increases quickly with time to approach a constant limit: the closer to the uncertain head boundary, the larger the error.The spatial variation in σ h can be clearly observed in Fig. 2h for fixed t .At any time, σ h is at its maximum at the right boundary (x = 1) where the head is uncertain, and it decreases quickly away from the boundary.The error in h(x, t) due to the uncertain head boundary is limited to a narrow zone near the boundary (x > 0.8) (Fig. 2h).
Finally, we consider the combined effects of the uncertainties from all four sources, i.e., the initial condition, sources, and flux and head boundaries.The head variance in Eq. ( 8) is written in the dimensionless form as where The dimensionless variances, σ 2 W 0 , σ 2 Q , and σ 2 H , need to be specified in order to evaluate the dimensionless σ 2 h (x , t ) in Eq. ( 15).For the typical aquifer mentioned above W λ W ) = 10 3 , and σ 2 H λ H /(σ 2 W λ W ) = 10 4 , and obtain σ 2 W 0 = 25, σ 2 Q = 0.1 and σ 2 H = 0.01.The changes of this σ h with x and t are presented in Fig. 2i and j, respectively.It is observed in Fig. 2i that at any location, the σ h is at its maximum due to the uncertainty in the initial condition, it gradually decreases with time, and approaches a constant limit at a later time (t > 0.6), which is due to the combined effects of the uncertain source/sink and flux and head boundaries.This means that the error in the water head at an early time is significant if the initial condition is uncertain and reduces with time to reach a constant limit.The error in the water head at a later time is determined by the uncertainties in the source/sink, and flux and head boundaries.It can be observed in Fig. 2j that σ h is relatively larger near both boundaries.The values of σ h at the two boundaries are equivalent (∼ 1.3) at an early time, say t = 0.01 (the top curve in Fig. 2j) and it reduces slowly away from the flux boundary, but quickly away from the head boundary.As time progresses, the σ h near the head boundary stays more or less the same but reduces significantly in most parts of the aquifer.This means that early on, the error in h(x, t) in most parts of the aquifer is mainly caused by the initial condition and at a later time it is due to the combined effects of the uncertain areal source/sink and flux boundary.The effect of the uncertain head boundary on h(x, t) does not significantly change with time, but it is limited to a narrow zone near the boundary.

Spectrum of groundwater levels
We first evaluated S hh in Eq. ( 10) due to the effect of the white noise flux boundary only by setting S QQ = 0, S W W = 0, and S H H = 0.The dimensionless spectrum S hh /S QQ as a function of the frequency (f ) was evaluated and presented in the log-log plot (Fig. 3a-c) for three values of t c (40, 400, and 4000 days), since the value of t c is 250 days for a sandy aquifer, as we mentioned above, and also at the six locations (x = 0.0, 0.2, 0.4, 0.6, 0.8, and 0.9).The spectrum S hh /S QQ in Fig. 3a is more or less horizontal (i.e., white noise) at low frequencies and it decreases gradually as f increases, indicating that an aquifer acts as a lowbass filter that filters signals at high frequencies and keeps signals at low frequencies.The aquifer significantly dampened the fluctuations of the groundwater level.The spectrum varies with the location x : the smaller the value of x or the closer to the left flux boundary (x = 0), the larger the spectrum (Fig. 3a-c).All spectra in Fig. 3a are not a straight line in the log-log plot, meaning that the temporal scaling of h(x, t) does not exist in the range of f = 10 −3 -10 0 when t c = 40 days.As t c increases to 400 and 4000 days, however, the spectrum at x = 0 becomes a straight line (the top curve in Fig. 3b and c) or has a power-law relation with f , i.e., S hh /S QQ ∝ 1/f , since its slope is approximately 1.The fluctuations of h(0, t) are pink noise due to the white noise fluctuations flux boundary when the characteristic timescale (t c ) is large which means that the aquifer is relatively less permeable and/or has a much larger horizontal length than its thickness.
Secondly, the spectrum S hh /S H H due to the sole effect of the random head boundary was evaluated by setting S H H = 0, S W W = 0, and S QQ = 0 in Eq. ( 10) for the same three values of t c and six locations and presented in Fig. 3d-f as a function of f .It is shown that similar to Fig. 3a-c, the spectrum decreases as f increases but different from Fig. 3ac, the spectrum is larger at x = 0.9 near the right boundary (the top curves in Fig. 3d-f) than at x = 0.0 (the bottom curves).Furthermore, none of the spectra are a straight line in the log-log plot, indicating that the temporal scaling of groundwater level fluctuations does not exist in the case of the white noise head boundary.
Thirdly, the spectrum S hh /S W W under the white noise recharge was evaluated by setting S W W = 0, S QQ = 0, and S H H = 0 in Eq. ( 10) for the same values of t c and x and presented in Fig. 3g-i as a function of f .It is shown that when t c = 40 days, the spectrum in Fig. 3g is horizontal at low frequencies and becomes a straight line at high frequencies: the closer to the right head boundary, the later it approaches a straight line (Fig. 3h).As t c increases to 400 and 4000 days, the slope of the spectrum at all locations, except at x = 0.9, approaches a straight line with a slope of 2 (Fig. 3h and i), indicating a temporal scaling of h(x, t).The fluctuations of groundwater level reflect a Brownian motion, i.e., S ∝ 1/f 2 , when t c ≥ 4000 days or in a relatively less permeable and/or has a much larger horizontal length than its thickness.
Finally, the head spectrum due to the combined effect of all three random sources (the white noise recharge, and flux and head boundaries) was evaluated, i.e., S W W = 0, S QQ = 0, and S H H = 0 in Eq. ( 10).The spectrum of S hh /S W W as a function of f is presented in Fig. 3j-l for the same values of t c and x , where S QQ /S W W = 1000 and S H H /S W W = 1000, which are the same as the values used in the previous section.It is noticed that the general patterns of S hh /S W W in the combined case are similar to the case of the random source/sink only (Fig. 3g-i), except at x = 0.0 and 0.9 (the dashed and dotted curves in Fig. 3j, respectively) due to the strong effects of the boundary conditions at these two locations.At t c = 4000 days, the spectra at all locations except x = 0.0 (Fig. 3l) are similar to those in Fig. 3i, indicating the dominating effect of the random areal source/sink.The spectrum at x = 0 in this case is also a straight line (the dashed curve in Fig. 3l) but with a different slope due to the effect of the random flux boundary which is similar to the top straight line in Fig. 3c.The above results provide a theoretical explanation as to why temporal scaling exists in the observed groundwater level fluctuations (Zhang and Schilling, 2004;Bloomfield and Little, 2010;Zhu et al., 2012).We thus conclude that x'=0.0 x'=0.9 x'=0.0 x'=0.9 x'=0.0 x'=0.9 c c t =40 day t =4000 day c t =400 day j.
Figure 3.The dimensionless power spectrum versus frequency (f ) at the dimensionless locations x = 0.0, 0.2, 0.4, 0.6, 0.8, and 0.9.The graphs on the left column show t c = 40 days, the graphs on the middle column show t c = 400 days, and the graphs on the right column show t c = 4000 days.The graphs on the first row show the dimensionless spectrum S hh /S QQ when S W W = 0, S H H = 0, and S QQ = 0 in Eq. ( 10); the graphs on the second row show S hh /S H H when S W W = 0, S QQ = 0, and S H H = 0; the graphs on the third row show S hh /S W W when S QQ = 0, S H H = 0, and S W W = 0; and the graphs on the bottom row show S hh /S W W when S QQ = 0, S H H = 0, and S W W = 0.
temporal scaling of h(x, t) may indeed exist in real aquifers due to the strong effect of the areal source/sink.

Conclusions
In this study the effects of random source/sink, and initial and boundary conditions on the uncertainty and temporal scaling of the groundwater level, h(x, t) were investigated.Analytical solutions were derived for the variance, covariance, and spectrum of h(x, t) in an unconfined aquifer, described by a linearized Boussinesq equation with white noise source/sink, and initial and boundary conditions.The standard deviations of h(x, t) for various cases were evaluated.Based on the results, the following conclusions can be drawn.
1.The error in h(x, t), due to a random initial condition, is significant at an early time, especially near a flux boundary.The duration for which the effect is significant may be a few hundred days in most aquifers.
2. The error in h(x, t) due to a random areal source/sink is significant in most parts of an aquifer: the closer to a flux boundary, the larger the error.
3. The errors in h(x, t) due to random flux and head boundaries are significant near the boundaries: the closer to the boundaries, the larger the errors.The random flux boundary may affect the head over a larger region than the random head boundary.
4. In the typical sandy aquifer studied (with the length of aquifer at the direction of water flow L = 100 m, the average saturated thickness M = 10 m, hydraulic conductivity K = 1 m day −1 , and specific yield S Y = 0.25) the error in h(x, t) at an early point in time is mainly caused by an uncertain initial condition, and the error reduces with time, reaching a constant error at a later time.The constant error in h(x, t) is mainly due to the combined effects of uncertain source/sink and boundaries.
5. The aquifer system behaves as a low-pass filter which filters the short-term (high frequencies) fluctuations and keeps the long-term (low frequencies) fluctuations.
6. Temporal scaling of groundwater level fluctuations may indeed exist in most parts of a low permeable aquifer whose horizontal length is much larger than its thickness, caused by the temporal fluctuations of areal source/sink.
Finally, it is pointed out that the analyses carried out in this study are under the assumption that the processes W (t), Q(t), and H (t) are uncorrelated white noise.In reality, they may be correlated and spatially varied.We plan to relax those constraints and study more realistic cases in the near future.
It is also noted that the analytical solutions for head variances derived in this study provide a way to identify and quantify the uncertainty.The spectrum relationship obtained among the head, recharge, and boundary conditions can help one to improve spectrum analysis for a groundwater level time series and remove the effects of the boundary conditions.

Figure 1 .
Figure 1.A schematic of the unconfined aquifer studied, where W (t) is the random time-dependent source/sink, H 0 (x) is the random initial condition, Q(t) is the random time-dependent flux at the left boundary, H (t) is the random time-dependent water level at the right boundary, L is the distance from the left to the right boundary, and h(x, t) is the random groundwater level in the aquifer.