Technical note: Discharge response of a confined aquifer with variable thickness to temporal, nonstationary, random recharge processes

This work develops a transfer function to describe the variation in the integrated specific discharge in response to the temporal variation in the rainfall event in the frequency domain. It is assumed that the rainfall–discharge process takes place in a confined aquifer with variable thickness, and it is treated as nonstationary in time to represent the stochastic nature of the hydrological process. The presented transfer function can be used to quantify the variability in the integrated discharge field induced by the variation in rainfall field or to simulate the discharge response of the system to any varying rainfall input at any time resolution using the convolution model. It is shown that, with the Fourier–Stieltjes representation approach, a closed-form expression for the transfer function in the frequency domain can be obtained, which provides a basis for the analysis of the influence of controlling parameters occurring in the rainfall rate and integrated discharge models on the transfer function.

1 Introduction

Quantifying the variability in the specific discharge response of an aquifer system to fluctuations in inflow recharge is essential for efficient groundwater resources management. However, this requires extensive and continuous hydrological time series data, and these data are very often not available in practice. A possible approach (namely convolution or transfer function approach) to this problem is to simulate the discharge response by convoluting the time-varying recharge input with the corresponding impulse response. In convolution models, the aquifer is regarded as a filter that converts recharge signals into fluctuations of the aquifer head or discharge. Lumped conceptual–convolution models have been shown to be an efficient means for the simulation of the time series of groundwater levels (e.g., Gelhar, 1974; Molénat et al., 1999; Olsthoorn, 2007; Long and Mahler, 2013; Pedretti et al., 2016).

Since the impulse response function in the convolution model contains all information of the system necessary to relate its input to its output, it may be determined from the analytical solution of the linear system equation governing the input–output process (e.g., Cooper and Rorabaugh, 1963). Once a suitable impulse response function can be specified, it allows the simulation of the linear system response to any varying input at any time resolution.

In this work, a regional-scale flow in a confined aquifer with variable thickness, which is recharged by rainfall through an outcrop, is analyzed by deriving transfer functions to characterize the rainfall–discharge process in the frequency domain. The stochastic analysis of groundwater flow is traditionally based on the assumption of the stationarity of the recharge and discharge processes. However, the hydrologic process in nature is nonstationary stochastic (e.g., Christensen and Lettenmaier, 2007; Milly et al., 2008; Sang et al., 2018). In order to improve the quantification of the natural recharge–discharge process, the nonstationary rainfall–discharge process is assumed in this study. The Fourier–Stieltjes representation approach is used to achieve the goal of this work. The analysis of the results is focused on the influence of controlling parameters in the rainfall–discharge models on the transfer function.

2 Problem formulation

In certain areas, aquifer recharge can vary greatly over time, so determining the discharge of the aquifer at the outlet for regional groundwater problems, which involves transferring recharge at the aquifer outcrop over a relatively large space scale, can be quite difficult. However, it is very important for the planning and management of regional groundwater resources that require knowledge of discharge at the aquifer outlet over a long period of time. This study is, therefore, devoted to quantifying the discharge response of the confined aquifer at the outlet to the temporal variation in aquifer recharge.

In this study, a confined aquifer with variable thickness is considered as a linear block box system, with a stochastic rainfall recharge input and, therefore, a stochastic runoff output. Both inputs and outputs are variable in time. In a linear system, the output of the system can be represented as a linear combination of the responses to each of the basic inputs through the convolution integral on a continuous timescale as (e.g., Rugh, 1981; Rinaldo and Marani, 1987) as follows:

$$\begin{array}{}\text{(1)}& Q\left(t\right)=\underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\phi }(t,\mathit{\tau })R\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau },\end{array}$$

where Q and R denote the output flow (discharge) rate and the input flow (recharge) rate of the system, respectively, and φ is the impulse response function of the system. As shown in Fig. 1, once an appropriate impulse response function can be specified at the scale of the aquifer, it is possible to evaluate the system response from records of the input without the need to specify smaller-scale heterogeneity. As will be shown below, the transfer function of the system can be used to characterize the uncertainty (variability) expected in applying the convolution integral Eq. (1) to the regional groundwater flow problems.

Figure 1 Schematic representation of a linear block box system.

When using the nonstationary Fourier–Stieltjes representations for the perturbed quantities of random recharge and outflow discharge processes, namely (e.g., Priestley, 1965) in the following:

$$\begin{array}{}\text{(2)}& {\displaystyle }r\left(t\right)=R\left(t\right)-E\left[R\right(t\left)\right]=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathrm{\Lambda }}_r(t;\mathit{\omega })\mathrm{d}{Z}_{\mathit{\xi }}\left(\mathit{\omega }\right),\text{(3)}& {\displaystyle }q\left(t\right)=Q\left(t\right)-E\left[Q\right(t\left)\right]=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathrm{\Lambda }}_q(t;\mathit{\omega })\mathrm{d}{Z}_{\mathit{\xi }}\left(\mathit{\omega }\right),\end{array}$$

the power spectrum of the mean removed convolution (1) can be written in the following form:

$$\begin{array}{}\text{(4)}& S_{qq}(t;\mathit{\omega })={\left|{\mathrm{\Lambda }}_q(t;\mathit{\omega })\right|}^{\mathrm{2}}{S}_{\mathit{\xi }\mathit{\xi }}\left(\mathit{\omega }\right),\end{array}$$

where, in the following:

$$\begin{array}{}\text{(5)}& {\mathrm{\Lambda }}_q(t;\mathit{\omega })=\underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\phi }(t,\mathit{\tau }){\mathrm{\Lambda }}_r(\mathit{\tau };\mathit{\omega })\mathrm{d}\mathit{\tau }.\end{array}$$

In Eqs. (2) and (3), Λ_r and Λ_q are the oscillatory functions (Priestley, 1965) of the recharge and outflow processes, respectively, ω is the frequency, ξ is a zero-mean random stationary forcing process, which generates the variations in the recharge and, thus, the output flow processes, with an orthogonal increment dZ_ξ. In Eq. (4), S_qq and S_ξξ represent the power spectra of the processes q and ξ, respectively, and $|{\mathrm{\Lambda }}_q{|}^{\mathrm{2}}$ is termed the transfer function.

In practice, the interest in many cases resides in evaluating the influence of the variation in the recharge on the variation in the outflow discharge. Equation (4) provides an efficient way of quantifying the variability in the outflow induced by the fluctuations in the inflow process in the frequency domain, since it relates the fluctuations in an output time series to those of an input series.

It is worthwhile mentioning that, for the case of second-order stationary rainfall processes, the representations of the forms (2) and (3) are reduced, respectively, to the following:

$$\begin{array}{}\text{(6)}& {\displaystyle }r\left(t\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathit{\omega }t}\mathrm{d}{Z}_r\left(\mathit{\omega }\right),\text{(7)}& {\displaystyle }q\left(t\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathrm{\Lambda }}_q(t;\mathit{\omega })\mathrm{d}{Z}_r\left(\mathit{\omega }\right),\end{array}$$

and, correspondingly, to the following:

$$\begin{array}{}\text{(8)}& S_{qq}(t;\mathit{\omega })={\left|{\mathrm{\Lambda }}_q(t:\mathit{\omega })\right|}^{\mathrm{2}}{S}_{\mathrm{rr}}\left(\mathit{\omega }\right),\end{array}$$

where, in the following:

$$\begin{array}{}\text{(9)}& {\mathrm{\Lambda }}_q(t;\mathit{\omega })=\underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\phi }(t,\mathit{\tau })e^{i\mathit{\omega }t}\mathrm{d}\mathit{\tau }.\end{array}$$

Equations (1) and (4) reveal that, once the transfer function for the linear lumped system is identified, the first two moments of the temporal random discharge fields can be determined. That is, the transfer function approach provides a basic framework for the characterization of large-scale flow processes, which may serve as a basis for an efficient management of groundwater resources. Furthermore, Eq. (4) provides another possible way to identify the aquifer parameters, as it relates the observed fluctuations of an output discharge process to those of a recharge process in the frequency domain.

In the following, the focus is on the development of a closed-form expression for the transfer function for a linear lumped confined flow model in which the regional confined aquifer is directly recharged by rainfall in the area corresponding to the high-elevation outcrop.

3 Theoretical development

The differential equation describing the transient flow of groundwater in inhomogeneous isotropic confined aquifers is of the following form (e.g., Bear, 1979; de Marsily, 1986):

$$\begin{array}{}\text{(10)}& S_{\mathrm{s}}{\displaystyle \frac{\partial }{\partial t}}h(\mathit{x},t)={\displaystyle \frac{\partial }{\partial x_i}}\left[K\right(\mathit{x}\left){\displaystyle \frac{\partial }{\partial x_i}}h\right(\mathit{x},t\left)\right]i=\mathrm{1},\mathrm{2},\mathrm{3},\end{array}$$

in which S_s represents the specific storage coefficient of the aquifer, $h=h(x,t$) is the hydraulic head, K(x) is the hydraulic conductivity, and x (= (x₁, x₂, x₃)) is the spatial coordinate vector. Many problems in groundwater flow are regional in nature, with the horizontal extent of the formation being much larger than the vertical extent. It is more practical to regard the flow as essentially horizontal. The regional-scale flow equations can be derived by integrating Eq. (10) along the thickness of the confined aquifer using the assumption of vertical equipotential surfaces (e.g., Bear, 1979; Bear and Cheng, 2010).

Integrating Eq. (10) along the x₃ axis perpendicular to the confining beds and using Leibniz's rule results in the following:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}S& (x_{\mathrm{1}},x_{\mathrm{2}}){\displaystyle \frac{\partial }{\partial x_i}}\stackrel{\mathrm{̃}}{h}(x_{\mathrm{1}},x_{\mathrm{2}},t)={\displaystyle \frac{\partial }{\partial x_i}}\left[T\right(x_{\mathrm{1}},x_{\mathrm{2}}\left){\displaystyle \frac{\partial }{\partial x_i}}\stackrel{\mathrm{̃}}{h}\right(x_{\mathrm{1}},x_{\mathrm{2}},t\left)\right]\\ & +T(x_{\mathrm{1}},x_{\mathrm{2}}){\displaystyle \frac{\partial }{\partial x_i}}\ln B(x_{\mathrm{1}},x_{\mathrm{2}}){\displaystyle \frac{\partial }{\partial x_i}}\stackrel{\mathrm{̃}}{h}(x_{\mathrm{1}},x_{\mathrm{2}},t),i=\mathrm{1},\mathrm{2},\end{array}\end{array}$$

where S(x₁, x₂) is the storage coefficient (or storativity) of the aquifer (=S_sB(x₁, x₂)), B(x₁, x₂)= b₂(x₁, x₂)−b₁(x₁, x₂) (an aquifer's thickness), b₁(x₁, x₂) and b₂(x₁, x₂) are the elevations of the fixed bottom and ceiling of the confined aquifer, respectively, T(x₁, x₂) is the transmissivity of the aquifer (=K(x₁, x₂)B(x₁, x₂)) interpreted as the depth-integrated hydraulic conductivity, and $\stackrel{\mathrm{̃}}{h}(x_{\mathrm{1}},x_{\mathrm{2}},t$) is the depth-averaged hydraulic head defined as follows:

$$\begin{array}{}\text{(12)}& \stackrel{\mathrm{̃}}{h}(x_{\mathrm{1}},x_{\mathrm{2}})={\displaystyle \frac{\mathrm{1}}{b_{\mathrm{2}}(x_{\mathrm{1}},x_{\mathrm{2}})-b_{\mathrm{1}}(x_{\mathrm{1}},x_{\mathrm{2}})}}\underset{b_{\mathrm{1}}(x_{\mathrm{1}},x_{\mathrm{2}})}{\overset{b_{\mathrm{2}}(x_{\mathrm{1}},x_{\mathrm{2}})}{\int }}h(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}},t)\mathrm{d}{x}_{\mathrm{3}}.\end{array}$$

Equation (11) is derived under the following assumptions: (1) there is no exchange of leakage fluxes between the confined aquifer and its confining beds in the direction of x₃ axis, (2) h(x₁, x₂, b₂, t) $\approx \stackrel{\mathrm{̃}}{h}(x_{\mathrm{1}},x_{\mathrm{2}},t)\approx h(x_{\mathrm{1}}$, x₂, b₁, t) (vertical equipotentials; Bear, 1979; Bear and Cheng, 2010), and (3) all terms involved in the fluxes in the directions of x₁ and x₂ at the boundaries are removed due to the no-slip condition at the boundaries.

The use of the depth-averaged hydraulic head operator for modeling regional groundwater flow is valid when the variation in aquifer thickness is much smaller than the average thickness (Bear, 1979; Bear and Cheng, 2010). The error introduced by the use of this operator is very small in most cases of practical interest, greatly simplifying the analysis of flow in confined aquifers.

Similarly, when applying Leibniz's rule to the Darcy equation, the vertically integrated specific discharge in the x_i direction is given by the following:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{Q}_{x_i}(x_{\mathrm{1}},x_{\mathrm{2}},t)& =-K(x_{\mathrm{1}},x_{\mathrm{2}})B(x_{\mathrm{1}},x_{\mathrm{2}}){\displaystyle \frac{\partial }{\partial x_i}}\stackrel{\mathrm{̃}}{h}(x_{\mathrm{1}},x_{\mathrm{2}},t)\\ & =-T(x_{\mathrm{1}},x_{\mathrm{2}}){\displaystyle \frac{\partial }{\partial x_i}}\stackrel{\mathrm{̃}}{h}(x_{\mathrm{1}},x_{\mathrm{2}},t),i=\mathrm{1},\mathrm{2}.\end{array}\end{array}$$

In this study, the regional confined aquifer is considered with a nonuniform, unidirectional mean flow in the x₁-axis direction but with small flow variations in the x₁-axis and x₂-axis directions and time-varying recharge at the aquifer outcrop (x₁=0). Since the regional flow domain considered in the x₁ direction is much larger than that in the x₂ direction, Eqs. (11) and (13) can be approximated as one-dimensional by the following:

$$\begin{array}{}\text{(14)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{S\left(x\right)}{\stackrel{\mathrm{‾}}{T}}}{\displaystyle \frac{\partial }{\partial t}}\stackrel{\mathrm{̃}}{h}(x,t)& ={\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}\stackrel{\mathrm{̃}}{h}(x,t)+{\displaystyle \frac{\partial }{\partial x}}\ln \stackrel{\mathrm{‾}}{T}\left(x\right){\displaystyle \frac{\partial }{\partial x}}\stackrel{\mathrm{̃}}{h}(x,t)\\ & +{\displaystyle \frac{\partial }{\partial x}}\ln B\left(x\right){\displaystyle \frac{\partial }{\partial x}}\stackrel{\mathrm{̃}}{h}(x,t)+{\displaystyle \frac{R\left(t\right)}{\stackrel{\mathrm{‾}}{T}}},\end{array}\text{(15)}& {\displaystyle }{Q}_x(x,t)=-\stackrel{\mathrm{‾}}{T}\left(x\right){\displaystyle \frac{\partial }{\partial x}}\stackrel{\mathrm{̃}}{h}(x,t),\end{array}$$

where $\stackrel{\mathrm{‾}}{T}=\stackrel{\mathrm{‾}}{K}B$, $\stackrel{\mathrm{‾}}{K}$ represents the spatial average of the hydraulic conductivity, and R is the recharge rate. It is worth noting that a one-dimensional flow equation with the transmissivity parameter has been widely used to predict the regional groundwater flow fields in the downstream region of the aquifer in field applications (e.g., Gelhar, 1974; Onder, 1998; Molénat et al., 1999; Russian et al., 2013). Equation (14) can be expressed alternatively as follows:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\displaystyle \frac{S_{\mathrm{s}}}{\stackrel{\mathrm{‾}}{K}}}{\displaystyle \frac{\partial }{\partial t}}\stackrel{\mathrm{̃}}{h}(x,t)& ={\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}\stackrel{\mathrm{̃}}{h}(x,t)+\mathrm{2}{\displaystyle \frac{\partial }{\partial x}}\ln B\left(x\right){\displaystyle \frac{\partial }{\partial x}}\stackrel{\mathrm{̃}}{h}(x,t)\\ & +{\displaystyle \frac{R\left(t\right)}{\stackrel{\mathrm{‾}}{K}B\left(x\right)}},\end{array}\end{array}$$

for the convenient analysis of the effect of the thickness of the aquifer.

In the following analysis, the recharge rate is considered a random function of time. Equation (15) is then regarded as a stochastic differential equation with a stochastic input in time and, therefore, a stochastic output in time. Introduction of the decomposition of the depth-averaged hydraulic head into a mean and zero-mean perturbation into Eq. (16), after subtracting the mean of the resulting equation from Eq. (16), means that the result is the following equation describing the depth-averaged head perturbation:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\displaystyle \frac{S_{\mathrm{s}}}{\stackrel{\mathrm{‾}}{K}}}{\displaystyle \frac{\partial }{\partial t}}{h}^{\prime }(x,t)& ={\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}{h}^{\prime }(x,t)+\mathrm{2}{\displaystyle \frac{\partial }{\partial x}}\ln B\left(x\right){\displaystyle \frac{\partial }{\partial x}}{h}^{\prime }(x,t)\\ & +{\displaystyle \frac{r\left(t\right)}{B\left(x\right)\stackrel{\mathrm{‾}}{K}}},\end{array}\end{array}$$

where $h^{\prime }(x,t)$ is the fluctuations in depth-averaged head.

If it is assumed that the thickness of the confined aquifer increases exponentially in the x direction, in accordance with Hantush (1962) and Marino and Luthin (1982), as follows:

$$\begin{array}{}\text{(18)}& B\left(x\right)=\mathit{\beta }{e}^{\mathit{\alpha }x},\end{array}$$

then Eq. (17) becomes the following:

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{S_{\mathrm{s}}}{\stackrel{\mathrm{‾}}{K}}}{\displaystyle \frac{\partial }{\partial t}}{h}^{\prime }(x,t)={\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}{h}^{\prime }(x,t)+\mathrm{2}\mathit{\alpha }{\displaystyle \frac{\partial }{\partial x}}{h}^{\prime }(x,t)+{\displaystyle \frac{e^{-\mathit{\alpha }x}}{\mathit{\beta }\stackrel{\mathrm{‾}}{K}}}r\left(t\right).\end{array}$$

In Eq. (18), β and α are positive geometrical parameters. Furthermore, the outcrop (x=0) and outlet (x=L) of the confined aquifer are considered as constant head boundaries. Since Eq. (19) only quantifies the response of the depth-averaged head to changes in the recharge rate, the initial and boundary conditions for Eq. (19) may be represented as follows:

$$\begin{array}{}\text{(20a)}& {\displaystyle }{h}^{\prime }(x,\mathrm{0};\mathit{\omega })=\mathrm{0},\text{(20b)}& {\displaystyle }{h}^{\prime }(\mathrm{0},t;\mathit{\omega })=\mathrm{0},\text{(20c)}& {\displaystyle }{h}^{\prime }(L,t;\mathit{\omega })=\mathrm{0}.\end{array}$$

The following Fourier–Stieltjes integral representation of a depth-averaged head process is used to solve Eqs. (19) and (20) for the fluctuations h′ in terms of r:

$$\begin{array}{}\text{(21)}& h^{\prime }(x,t)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathrm{\Lambda }}_{\mathrm{h}}(x,t;\mathit{\omega })\mathrm{d}{Z}_{\mathit{\xi }}\left(\mathit{\omega }\right),\end{array}$$

where Λ_h is the oscillatory function of depth-averaged head process. The resulting differential equation for the oscillatory functions is found from using Eqs. (2) and (21) in Eqs. (19) and (20) as follows:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{\displaystyle \frac{S_{\mathrm{s}}}{\stackrel{\mathrm{‾}}{K}}}{\displaystyle \frac{\partial }{\partial t}}{\mathrm{\Lambda }}_{\mathrm{h}}(x,t;\mathit{\omega })& ={\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}{\mathrm{\Lambda }}_{\mathrm{h}}(x,t;\mathit{\omega })+\mathrm{2}\mathit{\alpha }{\displaystyle \frac{\partial }{\partial x}}{\mathrm{\Lambda }}_{\mathrm{h}}(x,t;\mathit{\omega })\\ & +{\displaystyle \frac{e^{-\mathit{\alpha }x}}{\mathit{\beta }\stackrel{\mathrm{‾}}{K}}}{\mathrm{\Lambda }}_r(t;\mathit{\omega }),\end{array}\end{array}$$

with the following conditions:

$$\begin{array}{}\text{(23a)}& {\displaystyle }{\mathrm{\Lambda }}_{\mathrm{h}}(x,\mathrm{0};\mathit{\omega })=\mathrm{0},\text{(23b)}& {\displaystyle }{\mathrm{\Lambda }}_{\mathrm{h}}(\mathrm{0},t;\mathit{\omega })=\mathrm{0},\text{(23c)}& {\displaystyle }{\mathrm{\Lambda }}_{\mathrm{h}}(L,t;\mathit{\omega })=\mathrm{0}.\end{array}$$

By solving the above boundary value problem, the oscillatory function of depth-averaged head process is found to be the following (see Appendix A):

$$\begin{array}{}\text{(24)}& \begin{array}{rl}{\mathrm{\Lambda }}_{\mathrm{h}}(x,t;\mathit{\omega })& ={\displaystyle \frac{\mathrm{2}}{S_{\mathrm{s}}\mathit{\beta }}}\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(n\mathit{\pi }\right)}{n\mathit{\pi }}}\exp \left(-\mathit{\mu }{\displaystyle \frac{x}{L}}\right)\\ & \sin \left(n\mathit{\pi }{\displaystyle \frac{x}{L}}\right)\underset{\mathrm{0}}{\overset{t}{\int }}\exp [-{\mathit{\theta }}_n(t-\mathit{\tau }\left)\right]{\mathrm{\Lambda }}_r(\mathit{\tau };\mathit{\omega })\mathrm{d}\mathit{\tau },\end{array}\end{array}$$

where μ=αL and ${\mathit{\theta }}_n=\stackrel{\mathrm{‾}}{K}(n^{\mathrm{2}}{\mathit{\pi }}^{\mathrm{2}}+{\mathit{\mu }}^{\mathrm{2}})/(S_{\mathrm{s}}{L}^{\mathrm{2}}$). It implies from Eqs. (3), (15) and (24) that, at the arbitrary location x=x_ε, one finds the following:

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{\mathrm{\Lambda }}_q(t;\mathit{\omega })& ={\mathrm{\Lambda }}_{q_x}(x_{\mathit{\epsilon }},t;\mathit{\omega })\\ & =-\mathrm{2}{\displaystyle \frac{\stackrel{\mathrm{‾}}{K}}{S_{\mathrm{s}}L}}\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(n\mathit{\pi }\right)}{n\mathit{\pi }}}\left[n\mathit{\pi }\cos \right(n\mathit{\pi }\mathrm{{\rm Y}})\\ & -\mathit{\mu }\sin \left(n\mathit{\pi }\mathrm{{\rm Y}}\right)\left]\underset{\mathrm{0}}{\overset{t}{\int }}\exp \right[-{\mathit{\theta }}_n(t-\mathit{\tau })\left]{\mathrm{\Lambda }}_r\right(\mathit{\tau };\mathit{\omega })\mathrm{d}\mathit{\tau },\end{array}\end{array}$$

where $\mathrm{{\rm Y}}=x_{\mathit{\epsilon }}/L$. This means that the impulse response function of the system φ in Eqs. (1) or (5) is taken in the following form:

$$\begin{array}{}\text{(26)}& \begin{array}{rl}\mathit{\phi }(t,\mathit{\tau })& =-\mathrm{2}{\displaystyle \frac{\stackrel{\mathrm{‾}}{K}}{S_{\mathrm{s}}L}}\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(n\mathit{\pi }\right)}{n\mathit{\pi }}}\\ & \left[\cos \right(n\mathit{\pi }\mathrm{{\rm Y}})-\mathit{\mu }\sin (n\mathit{\pi }\mathrm{{\rm Y}}\left)\right]\exp [-{\mathit{\theta }}_n(t-\mathit{\tau }\left)\right].\end{array}\end{array}$$

4 Results and discussion

Equation (25) implies that the transfer function $|{\mathrm{\Lambda }}_q{|}^{\mathrm{2}}$ depends on the oscillatory function of the temporal random rainfall process; consequently, to complete the analysis of the transfer function, the oscillatory function of the temporal random rainfall process must be specified. It is assumed that the generated temporal random perturbations of the rainfall field are governed by the noise-forced diffusive rainfall model (North et al., 1993) as follows:

$$\begin{array}{}\text{(27)}& {\mathit{\tau }}_{\mathrm{0}}{\displaystyle \frac{\partial }{\partial t}}\mathit{\rho }(x,t)={\mathit{\lambda }}_{\mathrm{0}}^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}\mathit{\rho }(x,t)-\mathit{\rho }(x,t)+\mathit{\xi }\left(t\right),\end{array}$$

where ρ is a zero-mean rainfall rate perturbation, τ₀ and λ₀ are the characteristic timescales and length scales, respectively, which are inherent to the rainfall field, and ξ is a zero-mean random stationary forcing process which has a spectral representation of the following form (e.g., Lumley and Panofsky, 1964):

$$\begin{array}{}\text{(28)}& \mathit{\xi }\left(t\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i\mathit{\omega }t}\mathrm{d}{Z}_{\mathit{\xi }}\left(\mathit{\omega }\right).\end{array}$$

In Eq. (27), the rainfall rate field is represented as a first-order continuous autoregressive process in time and an isotropic second-order autoregressive process in space.

Furthermore, the rest of this study takes into account that rain falls within a defined period of time over a certain area of horizontal extension from $x=-\mathrm{\ell }$ to x=ℓ. As such, the initial and boundary conditions for rainfall rate perturbations may be represented by the following:

$$\begin{array}{}\text{(29a)}& {\displaystyle }\mathit{\rho }(x,\mathrm{0})=\mathrm{0},\text{(29b)}& {\displaystyle }\mathit{\rho }(-\mathrm{\ell },t)=\mathrm{0},\text{(29c)}& {\displaystyle }\mathit{\rho }(\mathrm{\ell },t)=\mathrm{0}.\end{array}$$

4.1 Nonstationary random rainfall fields in time

Using the Fourier–Stieltjes integral representation for the perturbation ρ, as follows:

$$\begin{array}{}\text{(30)}& \mathit{\rho }(x,t)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathrm{\Lambda }}_{\mathit{\rho }}(x,t;\mathit{\omega })\mathrm{d}{Z}_{\mathit{\xi }}\left(\mathit{\omega }\right),\end{array}$$

and Eq. (28) in Eq. (27), it follows that:

$$\begin{array}{}\text{(31)}& \begin{array}{rl}{\mathit{\tau }}_{\mathrm{0}}{\displaystyle \frac{\partial }{\partial t}}{\mathrm{\Lambda }}_{\mathit{\rho }}(x,t;\mathit{\omega })& ={\mathit{\lambda }}_{\mathrm{0}}^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}{\mathrm{\Lambda }}_{\mathit{\rho }}(x,t;\mathit{\omega })-{\mathrm{\Lambda }}_{\mathit{\rho }}(x,t;\mathit{\omega })\\ & +e^{i\mathit{\omega }t},\end{array}\end{array}$$

where Λ_ρ is the oscillatory function of the rainfall rate processes. With the application of the initial and boundary conditions, as follows:

$$\begin{array}{}\text{(32a)}& {\displaystyle }{\mathrm{\Lambda }}_{\mathit{\rho }}(x,\mathrm{0};\mathit{\omega })=\mathrm{0},\text{(32b)}& {\displaystyle }{\mathrm{\Lambda }}_{\mathit{\rho }}(-\mathrm{\ell },t;\mathit{\omega })=\mathrm{0},\text{(32c)}& {\displaystyle }{\mathrm{\Lambda }}_{\mathit{\rho }}(\mathrm{\ell },t;\mathit{\omega })=\mathrm{0},\end{array}$$

the solution of Eqs. (31) and (32) is given by the following (see Appendix B):

$$\begin{array}{}\text{(33)}& \begin{array}{rl}{\mathrm{\Lambda }}_{\mathit{\rho }}(x,t;\mathit{\omega })& =\mathrm{2}\sum _{m=\mathrm{1}}^{m=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(m\mathit{\pi }\right)}{m\mathit{\pi }}}\sin \left(m\mathit{\pi }{\displaystyle \frac{x+\mathrm{\ell }}{\mathrm{2}\mathrm{\ell }}}\right)\\ & {\displaystyle \frac{\exp \left(i{\mathrm{\Omega }}_t\right)-\exp (-{\mathrm{\Theta }}_{m}t/{\mathit{\tau }}_{\mathrm{0}})}{{\mathrm{\Theta }}_m+i\mathrm{\Gamma }}},\end{array}\end{array}$$

where ${\mathrm{\Theta }}_m=\mathrm{1}+m^{\mathrm{2}}{\mathit{\pi }}^{\mathrm{2}}{\mathit{\eta }}^{\mathrm{2}}$, $\mathit{\eta }={\mathit{\lambda }}_{\mathrm{0}}/\left(\mathrm{2}\mathrm{\ell }\right)$, Ω_t=ωt, and Γ=ωτ₀.

In the case where the regional confined aquifer is directly recharged by rainfall at the aquifer outcrop (x=0), the oscillatory function is reduced to the following:

$$\begin{array}{}\text{(34)}& \begin{array}{rl}{\mathrm{\Lambda }}_r(t;\mathit{\omega })& ={\mathrm{\Lambda }}_{\mathit{\rho }}(\mathrm{0},t;\mathit{\omega })=\mathrm{2}\sum _{m=\mathrm{1}}^{m=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(m\mathit{\pi }\right)}{m\mathit{\pi }}}\sin \left(m{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right)\\ & {\displaystyle \frac{\exp \left(i{\mathrm{\Omega }}_t\right)-\exp (-{\mathrm{\Theta }}_{m}t/{\mathit{\tau }}_{\mathrm{0}})}{{\mathrm{\Theta }}_m+i\mathrm{\Gamma }}}.\end{array}\end{array}$$

Correspondingly, the power spectrum of the rainfall rate, S_rr(t, ω), can be expressed by the following:

$$\begin{array}{}\text{(35)}& \begin{array}{rl}{S}_{\mathrm{rr}}(t;\mathit{\omega })& ={\left|{\mathrm{\Lambda }}_r(t;\mathit{\omega })\right|}^{\mathrm{2}}{S}_{\mathit{\xi }\mathit{\xi }}\left(\mathit{\omega }\right)\\ & =\mathrm{4}\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}\sum _{m=\mathrm{1}}^{m=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(m\mathit{\pi }\right)}{m\mathit{\pi }}}{\displaystyle \frac{\mathrm{1}-\cos \left(n\mathit{\pi }\right)}{n\mathit{\pi }}}\\ & \sin \left(m{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right)\sin \left(n{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right){\displaystyle \frac{\mathrm{1}}{{\mathrm{\Theta }}_m^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{1}}{{\mathrm{\Theta }}_n^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}}}}\\ & \mathit{\left\{}\right({\mathrm{\Theta }}_m{\mathrm{\Theta }}_n+{\mathrm{\Gamma }}^{\mathrm{2}}\left)\right[\mathrm{1}+T_{\mathrm{1}}-T_{\mathrm{2}}\cos \left({\mathrm{\Omega }}_t\right)]\\ & -T_{\mathrm{3}}\mathrm{\Gamma }({\mathrm{\Theta }}_m-{\mathrm{\Theta }}_n)\sin \left({\mathrm{\Omega }}_t\right)\mathit{\left\}}{S}_{\mathit{\xi }\mathit{\xi }}\right(\mathit{\omega }),\end{array}\end{array}$$

where T₁= exp[$-({\mathrm{\Theta }}_m+{\mathrm{\Theta }}_n)t/{\mathit{\tau }}_{\mathrm{0}}$], T₂= exp($-{\mathrm{\Theta }}_{m}t/{\mathit{\tau }}_{\mathrm{0}})+\exp (-{\mathrm{\Theta }}_{n}t/{\mathit{\tau }}_{\mathrm{0}}$), and T₃ = exp($-{\mathrm{\Theta }}_{m}t/{\mathit{\tau }}_{\mathrm{0}})-\exp (-{\mathrm{\Theta }}_{n}t/{\mathit{\tau }}_{\mathrm{0}}$).

The transfer function of the rainfall processes in Eq. (35) behaves like a filter, attenuating the high-frequency part of the rainfall spectrum. The graph of the transfer function, which is characterized by the characteristic timescale τ₀ for different characteristic length scales, is shown in Fig. 2. It clearly shows a reduction in the transfer function with increasing τ₀, implying a reduction in the variability in the rainfall field with the characteristic timescale of the rainfall field. A larger τ₀ decreases the temporal persistence of the rainfall fluctuations, resulting in a smaller transfer function. It is also seen that, for a fixed value of the timescale, the transfer function of the rainfall processes tends to decrease as the length scale of the rainfall field increases. The influence of the length scale plays a similar role to the influence of the timescale in reducing the temporal persistence of the rainfall fluctuations and, thus, the variability in the rainfall field.

Figure 2 Graphical representation of the transfer function of the rainfall processes in Eq. (35) characterized by the timescale for different length scales, where the series calculation is truncated up to $M=N=$ 100.

Through the use of Eqs. (25) and (34), the oscillatory function of the integrated discharge process could be represented as follows:

$$\begin{array}{}\text{(36)}& \begin{array}{rl}{\mathrm{\Lambda }}_q(t;\mathit{\omega })& =-\mathrm{4}{\displaystyle \frac{\stackrel{\mathrm{‾}}{K}}{S_{\mathrm{s}}L}}\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(n\mathit{\pi }\right)}{n\mathit{\pi }}}\left[n\mathit{\pi }\cos \right(n\mathit{\pi }\mathrm{{\rm Y}})\\ & -\mathit{\mu }\sin \left(n\mathit{\pi }\mathrm{{\rm Y}}\right)]\\ & \times \sum _{m=\mathrm{1}}^{m=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(m\mathit{\pi }\right)}{m\mathit{\pi }}}{\displaystyle \frac{\sin \left(m\frac{\mathit{\pi }}{\mathrm{2}}\right)}{{\mathrm{\Theta }}_m+i\mathrm{\Gamma }}}[{\displaystyle \frac{\exp \left(i{\mathrm{\Omega }}_t\right)-\exp (-{\mathit{\theta }}_{n}t)}{{\mathit{\theta }}_n+i\mathit{\omega }}}\\ & -{\displaystyle \frac{\exp (-{\mathrm{\Theta }}_{m}t/{\mathit{\tau }}_{\mathrm{0}})-\exp (-{\mathit{\theta }}_{n}t)}{{\mathit{\theta }}_n-{\mathrm{\Theta }}_m/{\mathit{\tau }}_{\mathrm{0}}}}].\end{array}\end{array}$$

Thus, the transfer function of the integrated discharge flux takes the following form:

$$\begin{array}{}\text{(37)}& \begin{array}{rl}{\displaystyle \frac{S_{qq}(t;\mathit{\omega })}{S_{\mathit{\xi }\mathit{\xi }}\left(\mathit{\omega }\right)}}& ={\left|{\mathrm{\Lambda }}_q(t;\mathit{\omega })\right|}^{\mathrm{2}}=\mathrm{16}{L}^{\mathrm{2}}{\mathit{\vartheta }}^{\mathrm{2}}\mathit{\left\{}\right[\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}\sum _{m=\mathrm{1}}^{m=\mathrm{\infty }}{\mathrm{\Psi }}_{\mathrm{1}}{\mathrm{\Psi }}_{\mathrm{2}}\\ & \left({\displaystyle \frac{{\mathrm{\Theta }}_m{\mathrm{\Psi }}_{\mathrm{3}}+\mathrm{\Gamma }{\mathrm{\Psi }}_{\mathrm{4}}}{{\mathit{\theta }}_n^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{0}}^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}}}}+{\displaystyle \frac{{\mathrm{\Theta }}_m{\mathrm{\Psi }}_{\mathrm{5}}}{{\mathrm{\Theta }}_m-{\mathit{\theta }}_n{\mathit{\tau }}_{\mathrm{0}}}}\right){]}^{\mathrm{2}}\\ & +\left[\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}\sum _{m=\mathrm{1}}^{m=\mathrm{\infty }}{\mathrm{\Psi }}_{\mathrm{1}}{\mathrm{\Psi }}_{\mathrm{2}}\right({\displaystyle \frac{{\mathrm{\Theta }}_m{\mathrm{\Psi }}_{\mathrm{4}}-\mathrm{\Gamma }{\mathrm{\Psi }}_{\mathrm{3}}}{{\mathit{\theta }}_n^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{0}}^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}}}}\\ & -{\displaystyle \frac{\mathrm{\Gamma }{\mathrm{\Psi }}_{\mathrm{5}}}{{\mathrm{\Theta }}_m-{\mathit{\theta }}_n{\mathit{\tau }}_{\mathrm{0}}}}\left){]}^{\mathrm{2}}\mathit{\right\}},\end{array}\end{array}$$

where $\mathit{\vartheta }=\stackrel{\mathrm{‾}}{K}{\mathit{\tau }}_{\mathrm{0}}/\left(S_{\mathrm{s}}{L}^{\mathrm{2}}\right)$, and in the following:

$$\begin{array}{}\text{(38a)}& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{{\mathrm{\Theta }}_m^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{1}-\cos \left(m\mathit{\pi }\right)}{m\mathit{\pi }}}\sin \left(m{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right),\text{(38b)}& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{2}}={\displaystyle \frac{\mathrm{1}-\cos \left(n\mathit{\pi }\right)}{n\mathit{\pi }}}\left[\cos \right(n\mathit{\pi }\mathrm{{\rm Y}})-\mathit{\mu }\sin (n\mathit{\pi }\mathrm{{\rm Y}}\left)\right],\text{(38c)}& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{3}}=\mathrm{\Gamma }\sin \left({\mathrm{\Omega }}_t\right)+{\mathit{\theta }}_n{\mathit{\tau }}_{\mathrm{0}}\left[\cos \right({\mathrm{\Omega }}_t)-\exp (-{\mathit{\theta }}_{n}t\left)\right],\text{(38d)}& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{4}}={\mathit{\theta }}_n{\mathit{\tau }}_{\mathrm{0}}\sin \left({\mathrm{\Omega }}_t\right)-\mathrm{\Gamma }\left[\cos \right({\mathrm{\Omega }}_t)-\exp (-{\mathit{\theta }}_{n}t\left)\right],\text{(38e)}& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{5}}=\exp (-{\mathrm{\Theta }}_{m}t/{\mathit{\tau }}_{\mathrm{0}})-\exp (-{\mathit{\theta }}_{n}t).\end{array}$$

Note that the linearity in modeling the recharge–discharge response of a catchment in Eq. (1), which was originally developed for large catchments, increases with catchment area (e.g., Chow et al., 1988). This implies that the impulse responses and transfer functions derived here are valid in large, confined aquifers.

An essential feature of the transfer function of the integrated discharge flux in Eq. (37) is the resulting filtering associated with the flow process, as shown in Fig. 3. The attenuation of the high-frequency part of the flow discharge spectrum means that the flow process smooths out much of the small-scale variations caused by the rainfall field. Physically, this feature implies that the flow field is much smoother than the rainfall field. The figure also shows that the transfer function at fixed values for frequency and time increases with the increasing thickness of the confined aquifer. An increase in the thickness of the aquifer leads to an increased temporal persistence of the flow discharge fluctuations caused by the variation in the rainfall field and, thus, to an increase in the variability in integrated discharge field. As shown in Fig. 4, the ratio of the mean hydraulic conductivity to the storage coefficient (often referred to as the aquifer diffusivity) plays a similar role in influencing the variation in the transfer function to the thickness of the confined aquifer. The introduction of a larger aquifer diffusivity leads to a larger transfer function of integrated discharge and, thus, to a larger variability in the discharge field. Since the variability in the discharge field is positively correlated with that of rainfall field, the variability in the integrated discharge field will decrease with the increasing characteristic timescale or length scale of the rainfall field (see Fig. 2).

Figure 3 Influence of the thickness of the confined aquifer on the transfer function of the discharge flux, where the series calculation is truncated up to $M=N=$ 100.

Figure 4 Influence of the aquifer diffusivity on the transfer function of the discharge flux, where the series calculation is truncated up to $M=N=$ 100.

From Eqs. (4) or (8), the transfer function can be defined as the ratio of the fluctuations of an observation of output time series to those of input time series in a frequency domain. Equations (35) and (37) indicate that the transfer functions are related to the properties of the rainfall field and the aquifer, such as the characteristic scales of the time and length of the rainfall field and the diffusivity and thickness parameters of the aquifer. Therefore, the transfer function derived here has the potential to perform a parameter estimation based on the observations of input and output time series using the inverse modeling approach.

The traditional approach to regional groundwater flow problems introduces the transmissivity term, namely the depth-integrated hydraulic conductivity operator, as follows:

$$\begin{array}{}\text{(39)}& T(x_{\mathrm{1}},x_{\mathrm{2}})=\underset{b_{\mathrm{1}}(x_{\mathrm{1}},x_{\mathrm{2}})}{\overset{b_{\mathrm{2}}(x_{\mathrm{1}},x_{\mathrm{2}})}{\int }}K(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}})\mathrm{d}{x}_{\mathrm{3}},\end{array}$$

into to the groundwater flow equation (diffusion equation) to reduce the three-dimensional equation to a two-dimensional one as follows:

$$\begin{array}{}\text{(40)}& \begin{array}{rl}S(x_{\mathrm{1}},x_{\mathrm{2}}){\displaystyle \frac{\partial }{\partial t}}h(x_{\mathrm{1}},x_{\mathrm{2}},t)& ={\displaystyle \frac{\partial }{\partial x_i}}\left[T\right(x_{\mathrm{1}},x_{\mathrm{2}}\left){\displaystyle \frac{\partial }{\partial x_i}}h\right(x_{\mathrm{1}},x_{\mathrm{2}},t\left)\right]\\ & i=\mathrm{1},\mathrm{2}.\end{array}\end{array}$$

This means that the effects of both the variation in K in x₃ direction and the aquifer thickness are implicitly reflected in the term T(x₁, x₂). This leads to great difficulties in assessing the influence of aquifer thickness on the flow field with Eq. (40).

The proposed diffusion equation of this work is the following:

$$\begin{array}{}\text{(41)}& \begin{array}{rl}{S}_{\mathrm{s}}& (x_{\mathrm{1}},x_{\mathrm{2}}){\displaystyle \frac{\partial }{\partial x_i}}\stackrel{\mathrm{̃}}{h}(x_{\mathrm{1}},x_{\mathrm{2}},t)={\displaystyle \frac{\mathrm{1}}{B(x_{\mathrm{1}},x_{\mathrm{2}})}}{\displaystyle \frac{\partial }{\partial x_i}}\\ & \left[K\right(x_{\mathrm{1}},x_{\mathrm{2}}\left)B\right(x_{\mathrm{1}},x_{\mathrm{2}}\left){\displaystyle \frac{\partial }{\partial x_i}}\stackrel{\mathrm{̃}}{h}\right(x_{\mathrm{1}},x_{\mathrm{2}},t\left)\right]\\ & +K(x_{\mathrm{1}},x_{\mathrm{2}}){\displaystyle \frac{\partial }{\partial x_i}}\ln B(x_{\mathrm{1}},x_{\mathrm{2}}){\displaystyle \frac{\partial }{\partial x_i}}\stackrel{\mathrm{̃}}{h}(x_{\mathrm{1}},x_{\mathrm{2}},t)i=\mathrm{1},\mathrm{2},\end{array}\end{array}$$

which, when derived by the hydraulic approach (Bear, 1979; Bear and Cheng, 2010), provides an efficient way of analyzing flow fields in confined aquifers of nonuniform thickness. Note that Eq. (41) is the reformulation of Eq. (11). In addition, the usual observations of flow in porous media are measurements of hydraulic head from wells screened over extended sections of the medium. The measurement at a given location approximately represents a depth-averaged actual hydraulic head resulting from flow through a three-dimensional hydraulic conductivity field across the thickness of the medium. This means that the depth-averaged head representation used in Eq. (41) is consistent with what is observed in the fields.

Climate changes have a direct influence on the rainfall event (e.g., Trenberth, 2011; Pendergrass et al., 2014; Eekhout et al., 2018). The nonstationarity in the statistical properties of a rainfall field is a representation of climate change (e.g., Razavi et al., 2015; López and Francés, 2013; Benoit et al., 2020). The nonstationary effect of climatic change over time on variability in groundwater-specific discharge has not yet been well characterized. The transfer function in Eq. (37), which relates the nonstationary spectra of the rainfall fluctuations to those of integrated discharge variation, generalizes existing studies that considered stationary recharge and discharge fields. To our knowledge, it has not been previously presented in the literature and has the potential for analyzing the effects of climate change on temporal groundwater-specific discharge variability.

4.2 Application in the prediction of outflow discharge

The usefulness of the stochastic theory presented here lies in its essentially predictive nature. The variance can be used as a quantification of the uncertainty associated with the prediction in field situations using the linear system model. In this sense, the solution of Eq. (1) ±, i.e., two times the square root of the variance, provides a rational framework for predicting discharge over a relatively large spatial scale where direct observations of such a dependent variable are not possible.

For large times, the first term in Eq. (37) dominates the sum of the other terms, and therefore, the transfer function can be approximated by the following:

$$\begin{array}{}\text{(42)}& \begin{array}{rl}{\left|{\mathrm{\Lambda }}_q(t;\mathit{\omega })\right|}^{\mathrm{2}}& ={\displaystyle \frac{\mathrm{256}}{{\mathit{\pi }}^{\mathrm{2}}}}{L}^{\mathrm{2}}{\mathit{\vartheta }}^{\mathrm{2}}\left[\mathit{\pi }\cos \right(\mathit{\pi }\mathrm{{\rm Y}})-\mathit{\mu }\sin (\mathit{\pi }\mathrm{{\rm Y}}){]}^{\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{{\mathrm{\Theta }}_{\mathrm{1}}^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}}}}\\ & \mathit{\{}{\mathrm{\Xi }}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{{\forall }^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}}}}[\mathrm{1}+\mathrm{2}\forall \mathrm{\Xi }{T}_A+T_A^{\mathrm{2}}-\mathrm{2}(\forall \mathrm{\Xi }+T_A)\\ & \cos \left({\mathrm{\Omega }}_t\right)-\mathrm{2}\mathrm{\Xi }\mathrm{\Gamma }\sin \left({\mathrm{\Omega }}_t\right)\left]\mathit{\right\}},\end{array}\end{array}$$

where Θ₁= 1+π²η², $\forall =\stackrel{\mathrm{‾}}{K}{\mathit{\tau }}_{\mathrm{0}}({\mathit{\pi }}^{\mathrm{2}}+{\mathit{\mu }}^{\mathrm{2}})/(S_{\mathrm{s}}{L}^{\mathrm{2}}$), Ξ= ($T_R-T_A)/(\forall -{\mathrm{\Theta }}_{\mathrm{1}}$), $T_R=\exp (-{\mathrm{\Theta }}_{\mathrm{1}}t/{\mathit{\tau }}_{\mathrm{0}}$), and T_A= exp($-\forall t/{\mathit{\tau }}_{\mathrm{0}}$). If the variation in the rainfall event is generated by a random white noise forcing, the variance in the outflow discharge at large timescales can then be calculated using Eq. (42) as follows:

$$\begin{array}{}\text{(43)}& \begin{array}{rl}{\mathit{\sigma }}_q^{\mathrm{2}}\left(t\right)& =\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{S}_{qq}(t;\mathit{\omega })\mathrm{d}\mathit{\omega }=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\left|{\mathrm{\Lambda }}_q(t;\mathit{\omega })\right|}^{\mathrm{2}}{S}_{\mathit{\xi }\mathit{\xi }}\left(\mathit{\omega }\right)\mathrm{d}\mathit{\omega }\\ & ={\displaystyle \frac{\mathrm{256}}{\mathit{\pi }}}{\displaystyle \frac{G_{\mathrm{0}}{L}^{\mathrm{2}}}{{\mathit{\tau }}_{\mathrm{0}}}}{\mathit{\vartheta }}^{\mathrm{2}}\left[\mathit{\pi }\cos \right(\mathit{\pi }\mathrm{{\rm Y}})-\mathit{\mu }\sin (\mathit{\pi }\mathrm{{\rm Y}}\left){]}^{\mathrm{2}}\mathit{\right\{}{\displaystyle \frac{{\mathrm{\Xi }}^{\mathrm{2}}}{{\mathrm{\Theta }}_{\mathrm{1}}}}\\ & +{\displaystyle \frac{\mathrm{1}+\mathrm{2}\forall \mathrm{\Xi }{T}_A+T_A^{\mathrm{2}}}{\forall {\mathrm{\Theta }}_{\mathrm{1}}\left(\forall +{\mathrm{\Theta }}_{\mathrm{1}}\right)}}\\ & -\mathrm{2}(\forall \mathrm{\Xi }+T_A){\displaystyle \frac{{\mathrm{\Theta }}_{\mathrm{1}}{T}_A-\forall T_R}{\forall {\mathrm{\Theta }}_{\mathrm{1}}({\mathrm{\Theta }}_{\mathrm{1}}^{\mathrm{2}}-{\forall }^{\mathrm{2}})}}-\mathrm{2}\mathrm{\Xi }{\displaystyle \frac{T_R-T_A}{{\mathrm{\Theta }}_{\mathrm{1}}^{\mathrm{2}}-{\forall }^{\mathrm{2}}}}\mathit{\}},\end{array}\end{array}$$

where G₀ represents a constant spectral density of a white noise process. Note that white noise is a signal that contains all frequencies in equal proportions; that is, a signal for which the spectrum is flat.

After observing the recharge rate R(t) over time at the outcrop of the aquifer and identifying input parameters, such as the specific storage coefficient, mean hydraulic conductivity and geometrical parameters of the aquifer and the characteristic timescales and length scales of the rainfall event for a given area or region, the discharge can be determined under uncertainty in the far downstream aquifer area, i.e., Eq. (1) together with Eq. (26) ± being two times the square root of Eq. (43). It provides an important basis for the rational management of regional groundwater resources in complex geologic settings under uncertainty.

4.3 A note on stationary random rainfall fields in time

If the temporal random rainfall fields are stationary, there exists a representation of the rainfall perturbation process in terms of a Fourier–Stieltjes integral in the form of Eq. (6). Substituting Eqs. (6) and (21) into Eq. (19) gives the following:

$$\begin{array}{}\text{(44)}& \begin{array}{rl}{\displaystyle \frac{S_{\mathrm{s}}}{\stackrel{\mathrm{‾}}{K}}}{\displaystyle \frac{\partial }{\partial t}}{\mathrm{\Lambda }}_{\mathrm{h}}(x,t;\mathit{\omega })& ={\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}{\mathrm{\Lambda }}_{\mathrm{h}}(x,t;\mathit{\omega })+\mathrm{2}\mathit{\alpha }{\displaystyle \frac{\partial }{\partial x}}{\mathrm{\Lambda }}_{\mathrm{h}}(x,t;\mathit{\omega })\\ & +{\displaystyle \frac{e^{-\mathit{\alpha }x}}{\mathit{\beta }\stackrel{\mathrm{‾}}{K}}}{e}^{i\mathit{\omega }t}.\end{array}\end{array}$$

The solution of Eq. (44) with conditions Eq. (23) is as follows:

$$\begin{array}{}\text{(45)}& \begin{array}{rl}{\mathrm{\Lambda }}_{\mathrm{h}}(x,t;\mathit{\omega })& ={\displaystyle \frac{\mathrm{2}}{S_{\mathrm{s}}\mathit{\beta }}}\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(n\mathit{\pi }\right)}{n\mathit{\pi }}}\exp \left(-\mathit{\mu }{\displaystyle \frac{x}{L}}\right)\\ & \sin \left(n\mathit{\pi }{\displaystyle \frac{x}{L}}\right){\displaystyle \frac{\exp \left(i{\mathrm{\Omega }}_t\right)-\exp (-{\mathit{\theta }}_{n}t)}{{\mathit{\theta }}_n+i\mathit{\omega }}},\end{array}\end{array}$$

so that, in the following:

$$\begin{array}{}\text{(46)}& \begin{array}{rl}{\mathrm{\Lambda }}_q(t;\mathit{\omega })& =-\mathrm{2}{\displaystyle \frac{\stackrel{\mathrm{‾}}{K}}{S_{\mathrm{s}}L}}\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(n\mathit{\pi }\right)}{n\mathit{\pi }}}\left[n\mathit{\pi }\cos \right(n\mathit{\pi }\mathrm{{\rm Y}})\\ & -\mathit{\mu }\sin \left(n\mathit{\pi }\mathrm{{\rm Y}}\right)]{\displaystyle \frac{\exp \left(i{\mathrm{\Omega }}_t\right)-\exp (-{\mathit{\theta }}_{n}t)}{{\mathit{\theta }}_n+i\mathit{\omega }}},\end{array}\end{array}$$

and, thus, in the following:

$$\begin{array}{}\text{(47)}& \begin{array}{rl}{\left|{\mathrm{\Lambda }}_q(t;\mathit{\omega })\right|}^{\mathrm{2}}& =\mathrm{4}{L}^{\mathrm{2}}{\mathit{\vartheta }}^{\mathrm{2}}\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}\sum _{m=\mathrm{1}}^{m=\mathrm{\infty }}{\displaystyle \frac{\mathrm{\Phi }\left(m\right)\mathrm{\Phi }\left(n\right)}{({\mathit{\theta }}_n^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{0}}^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}})({\mathit{\theta }}_m^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{0}}^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}})}}\\ & \left[\right({\mathit{\theta }}_m{\mathit{\theta }}_n{\mathit{\tau }}_{\mathrm{0}}^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}}\left)\right(\mathrm{1}+{\mathrm{\Delta }}_{\mathrm{1}}-\cos \left({\mathrm{\Omega }}_t\right){\mathrm{\Delta }}_{\mathrm{2}})\\ & -\mathrm{\Gamma }\sin \left({\mathrm{\Omega }}_t\right)({\mathit{\theta }}_m-{\mathit{\theta }}_n){\mathit{\tau }}_{\mathrm{0}}{\mathrm{\Delta }}_{\mathrm{3}}],\end{array}\end{array}$$

where, in the following:

$$\begin{array}{}\text{(48)}& \mathrm{\Phi }\left(y\right)={\displaystyle \frac{\mathrm{1}-\cos \left(y\mathit{\pi }\right)}{y\mathit{\pi }}}\left[y\mathit{\pi }\cos \right(y\mathit{\pi }\mathrm{{\rm Y}})-\mathit{\mu }\sin (y\mathit{\pi }\mathrm{{\rm Y}}\left)\right],\end{array}$$

Δ₁= exp[-(θ_m+θ_n)t], Δ₂= exp(−θ_mt) + exp(−θ_nt), and Δ₃= exp(−θ_mt)-exp(−θ_nt).

At large timescales, Eq. (35) approaches a finite value as follows:

$$\begin{array}{}\text{(49)}& \begin{array}{rl}{S}_{\mathrm{rr}}\left(\mathit{\omega }\right)& =\mathrm{4}\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}\sum _{m=\mathrm{1}}^{m=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(m\mathit{\pi }\right)}{m\mathit{\pi }}}{\displaystyle \frac{\mathrm{1}-\cos \left(n\mathit{\pi }\right)}{n\mathit{\pi }}}\\ & \sin \left(m{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right)\sin \left(n{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right){\displaystyle \frac{{\mathrm{\Theta }}_m{\mathrm{\Theta }}_n+{\mathrm{\Gamma }}^{\mathrm{2}}}{({\mathrm{\Theta }}_m^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}})({\mathrm{\Theta }}_n^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}})}}{S}_{\mathit{\xi }\mathit{\xi }}\left(\mathit{\omega }\right),\end{array}\end{array}$$

and the corresponding rainfall process is stationary. Combining Eqs. (47) with (49) gives the following:

$$\begin{array}{}\text{(50)}& \begin{array}{rl}{\displaystyle \frac{S_{qq}\left(\mathit{\omega }\right)}{S_{\mathit{\xi }\mathit{\xi }}\left(\mathit{\omega }\right)}}& =\mathrm{16}{L}^{\mathrm{2}}{\mathit{\vartheta }}^{\mathrm{2}}\mathit{\{}\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}\sum _{m=\mathrm{1}}^{m=\mathrm{\infty }}\\ & {\displaystyle \frac{\mathrm{\Phi }\left(m\right)\mathrm{\Phi }\left(n\right)}{({\mathit{\theta }}_n^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{0}}^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}})({\mathit{\theta }}_m^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{0}}^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}})}}\left[\right({\mathit{\theta }}_m{\mathit{\theta }}_n{\mathit{\tau }}_{\mathrm{0}}^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}})\\ & (\mathrm{1}+{\mathrm{\Delta }}_{\mathrm{1}}-\cos ({\mathrm{\Omega }}_t\left){\mathrm{\Delta }}_{\mathrm{2}}\right)-\mathrm{\Gamma }\sin \left({\mathrm{\Omega }}_t\right)({\mathit{\theta }}_m-{\mathit{\theta }}_n){\mathit{\tau }}_{\mathrm{0}}{\mathrm{\Delta }}_{\mathrm{3}}\left]\mathit{\right\}}\\ & \times [\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}\sum _{m=\mathrm{1}}^{m=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(m\mathit{\pi }\right)}{m\mathit{\pi }}}{\displaystyle \frac{\mathrm{1}-\cos \left(n\mathit{\pi }\right)}{n\mathit{\pi }}}\sin \left(m{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right)\\ & \sin \left(n{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right){\displaystyle \frac{{\mathrm{\Theta }}_m{\mathrm{\Theta }}_n+{\mathrm{\Gamma }}^{\mathrm{2}}}{({\mathrm{\Theta }}_m^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}})({\mathrm{\Theta }}_n^{\mathrm{2}}+{\mathrm{\Gamma }}^{\mathrm{2}})}}].\end{array}\end{array}$$

Note that the nonstationarity in the hydraulic head or integrated discharge is introduced by a nonuniform thickness of the confined aquifer, even if the recharge field is stationary. Nonuniformity in the mean flow, for example, can also cause the nonstationarity in the statistics of random flow fields in heterogeneous aquifers (e.g., Rubin and Bellin, 1994; Ni and Li, 2006; Ni et al., 2010).

5 Conclusions

An analytical transfer function is developed to describe the spectral response characteristics of confined aquifers with variable thickness to the variation in the rainfall field, where the aquifer is directly recharged by rainfall at the outcrop of the aquifer. The rainfall–discharge process is treated as nonstationary in time, as it reflects the stochastic nature of the hydrological process. Any varying rainfall input at any time resolution can be convolved with the transfer function (or impulse response function) to simulate any discharge output of a linear model. The transfer function derived here, which relates the nonstationary spectra of the rainfall fluctuations to those of integrated discharge variation, has the potential to analyze the influence of climate change on groundwater recharge variability.

The closed-form results of this work are developed on the basis of the Fourier–Stieltjes representation approach, which allows us to analyze the effects of the controlling parameters in the models on the transfer function of the integrated discharge. It is found that the persistence in rainfall fluctuations is greater for a smaller value of the characteristic timescale or length scale of the rainfall field, which, in turn, leads to greater variability in the integrated discharge field. The attenuating characteristic of the confined aquifer flow system is observed in the spectral domain. The variability in the integrated discharge in a confined aquifer with variable thickness is increased with the thickness parameter α. The larger the aquifer diffusivity, the greater the spectrum (variability) of the integrated discharge.

Appendix A: Evaluation of Λ_h in Eq. (20)

The boundary value problem describing the depth-averaged head fluctuations induced by the variation in recharge rate in frequency domain is given by Eqs. (22) and (23). Using the transformation, as follows:

$$\begin{array}{}\text{(A1)}& {\mathrm{\Lambda }}_{\mathrm{h}}(x,t;\mathit{\omega })=\exp [-\mathit{\alpha }\left(x+{\displaystyle \frac{\mathit{\alpha }\stackrel{\mathrm{‾}}{K}}{S_{\mathrm{s}}}}t\right)]U(x,t;\mathit{\omega }).\end{array}$$

Equation (22) in Λ_h(x, t; ω), together with Eq. (23), can be converted into a new (easier) one with a new variable U(x, t; ω) as follows:

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}U(x,t;\mathit{\omega })& ={\displaystyle \frac{\stackrel{\mathrm{‾}}{K}}{S_{\mathrm{s}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}U(x,t;\mathit{\omega })\\ & +{\displaystyle \frac{\mathrm{1}}{\mathit{\beta }{S}_{\mathrm{s}}}}\exp \left({\displaystyle \frac{\stackrel{\mathrm{‾}}{K}{\mathit{\alpha }}^{\mathrm{2}}}{S_{\mathrm{s}}}}t\right){\mathrm{\Lambda }}_r(t;\mathit{\omega }),\end{array}\end{array}$$

with the following:

$$\begin{array}{}\text{(A3a)}& {\displaystyle }U(x,\mathrm{0};\mathit{\omega })=\mathrm{0},\text{(A3b)}& {\displaystyle }U(\mathrm{0},t;\mathit{\omega })=\mathrm{0},\text{(A3c)}& {\displaystyle }U(L,t;\mathit{\omega })=\mathrm{0}.\end{array}$$

The solution of Eqs. (A2) and (A3) can be found through the technique of the separation of variables (e.g., Farlow, 1993) as follows:

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}U(x,t;\mathit{\omega })& ={\displaystyle \frac{\mathrm{2}}{S_{\mathrm{s}}\mathit{\beta }}}\sum _{n=\mathrm{1}}^{n=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(n\mathit{\pi }\right)}{n\mathit{\pi }}}\sin \left(n\mathit{\pi }{\displaystyle \frac{x}{L}}\right)\underset{\mathrm{0}}{\overset{t}{\int }}\\ & \exp [-{\mathit{\upsilon }}_n(t-\mathit{\tau }\left)\right]\exp \left({\displaystyle \frac{\stackrel{\mathrm{‾}}{K}}{S_{\mathrm{s}}}}{\mathit{\alpha }}^{\mathrm{2}}\mathit{\tau }\right){\mathrm{\Lambda }}_r(\mathit{\tau };\mathit{\omega })\mathrm{d}\mathit{\tau },\end{array}\end{array}$$

where ${\mathit{\upsilon }}_n=\stackrel{\mathrm{‾}}{K}{n}^{\mathrm{2}}{\mathit{\pi }}^{\mathrm{2}}/(S_{\mathrm{s}}{L}^{\mathrm{2}}$). With reference to Eq. (A1), the solution to Eqs. (22) and (23) is then given by Eq. (24).

Appendix B: Evaluation of Λ_ρ in Eq. (31)

Making use of the following transformation:

$$\begin{array}{}\text{(B1)}& {\mathrm{\Lambda }}_{\mathit{\rho }}(x,t;\mathit{\omega })=\exp \left(-{\displaystyle \frac{t}{{\mathit{\tau }}_{\mathrm{0}}}}\right)u(x,t;\mathit{\omega }),\end{array}$$

leads Eqs. (31) and (32) to the following:

$$\begin{array}{}\text{(B2)}& {\displaystyle \frac{\partial }{\partial t}}u(x,t;\mathit{\omega })={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{\tau }}_{\mathrm{0}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}u(x,t;\mathit{\omega })+{\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{0}}}}\exp \left[\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{0}}}}+i\mathit{\omega }\right)t\right],\end{array}$$

with the following:

$$\begin{array}{}\text{(B3a)}& {\displaystyle }u(x,\mathrm{0};\mathit{\omega })=\mathrm{0},\text{(B3b)}& {\displaystyle }u(-\mathrm{\ell },t;\mathit{\omega })=\mathrm{0},\text{(B3c)}& {\displaystyle }u(\mathrm{\ell },t;\mathit{\omega })=\mathrm{0}.\end{array}$$

In a similar way, based on the technique of the separation of variables, Eqs. (B2) and (B3) arrive at the solution in the following form:

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}u(x,t;\mathit{\omega })& =\mathrm{2}\sum _{m=\mathrm{1}}^{m=\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\cos \left(m\mathit{\pi }\right)}{m\mathit{\pi }}}\sin \left(m\mathit{\pi }{\displaystyle \frac{x+\mathrm{\ell }}{\mathrm{2}\mathrm{\ell }}}\right)\\ & {\displaystyle \frac{\exp \left[\right(\mathrm{1}+i\mathrm{\Gamma })t/{\mathit{\tau }}_{\mathrm{0}}]-\exp (-{\mathit{\varsigma }}_{m}t/{\mathit{\tau }}_{\mathrm{0}})}{{\mathrm{\Theta }}_m+i\mathrm{\Gamma }}},\end{array}\end{array}$$

where ς_m=m²π²η², η= ${\mathit{\lambda }}_{\mathrm{0}}/(\mathrm{2}\mathrm{\ell }$), ${\mathrm{\Theta }}_m=\mathrm{1}+{\mathit{\varsigma }}_m$, and Γ=ωτ₀. The use of Eqs. (B1) and (B4) results in Eq. (33).