Stochastic analysis of the variability of groundwater flow fields in heterogeneous confined aquifers of variable thickness

The problem of flow through heterogeneous confined aquifers of variable thickness is analyzed from a stochastic point of view. The analysis is carried out on the basis of the integrated equations of the depth-averaged hydraulic head and integrated specific discharge, which are developed by integrating the continuity equation and equation for the specific discharge over the thickness, respectively. A spectrally based perturbation approach is used to arrive at the general results for the statistics of the flow fields in the Fourier domain, such as the variance of the depth-averaged head, and the mean and variance of integrated discharge. However, the closed-form expressions are obtained under the condition of steady unidirectional mean flow in the horizontal plane. In developing stochastic solutions, the input hydraulic conductivity parameter is viewed as a spatial random field characterized by the theoretical spatial covariance function. The evaluation of the closed-form solutions focuses on the influence of the controlling parameters, namely as a geometrical parameter defining the variation of the aquifer thickness and the correlation scale of log hydraulic conductivity, on the variability of the fluid fields. The application of the present stochastic theory to predict the total specific discharge under uncertainty using the field data is also provided.


Introduction
There is an increasing demand for science-based approaches to the management and use of groundwater resources at regional scale (e. g., Zhou and Li 2011; Barthel and Banzhaf 2015;Sishodia et al. 2017;Pétré et al. 2019).
Hydrologists are therefore often interested in the solutions of regional groundwater flow models. The solutions may be required for purpose of prediction and design. As such, it is important in constructing a regional flow model that may accurately represent the effects of hydrological changes on the behavior of the aquifer.
Groundwater flow at regional scale is concerned with the movement of flow through an aquifer whose lateral extent of the formation is much larger than the formation thickness (e. g., Nastev et al. 2008;Vassena et al. 2012;Rotiroti et al. 2019). When dealing with regional groundwater flow problems, it is more practical to treat the flow essentially horizontal and neglect the vertical component of the flow. That is, the flow behavior at regional scale can be appropriately described by a two-dimensional continuity equation, whereby the parameters such as the specific storage coefficient and hydraulic conductivity are replaced by the storativity and transmissivity, respectively. The output of the two-dimensional equation is then interpreted as the depth-averaged hydraulic head. The assumption of an essentially horizontal flow helps by reducing the effective dimension of the problem by one, which simplifies the problem, and the two-dimensional problem needs less detailed knowledge about the medium properties. Note that the analysis of regional groundwater flow based on the two-dimensional horizontal flow approximation is referred to as the hydraulic approach to flow in aquifers (Bear 1979;Bear and Cheng 2010).
Most research on essentially horizontal flow problems in confined aquifers considers the thickness of the aquifer to be uniform throughout the entire flow domain. However, natural aquifers often exhibit nonuniformity in the thickness of confined aquifers at a regional scale (Refsgaard 1997;Leray et al. 2012;Rotzoll et al. 2013;Andreasen et al. 2013). Hence, there arise pertinent questions as follows: (1) How can the existing model, which is limited to the flow through the confined aquifer of constant thickness, be generalized to the case of flow through the aquifer of variable thickness? (2) How does the thickness of the aquifer affect the flow field? These questions define the theme of this paper. The influence of the nonuniform thickness of confined aquifers on the flow behavior will be analyzed within a stochastic framework.
It is generally recognized that the movement of groundwater is strongly influenced by the complex heterogeneity of natural aquifer materials, the details of which cannot be precisely predicted. Therefore, for accurate predictions of the groundwater movement, a large number of measurements are required to represent the actual heterogeneity of the aquifer. This situation is, however, far from that encountered in practice. The amount of information on all coefficients and parameters appearing in the groundwater prediction model for accurate predictions is never sufficient, so that a certain degree of uncertainty (variability) about the average predictions must be accepted.
For the problems of flow in heterogeneous aquifers, a logical approach such as a stochastic methodology may be to ignore the details of the fine-scale variations and to develop a theoretical basis for the description of the largescale flow behavior in terms of effective parameters (e.g., Gelhar 1986;Dagan and Neuman 1997;Christakos 2003). The next step is to quantify the degree of variability (uncertainty) about the predicted mean. This approach will be appropriate if the overall scale of the problem is large compared to the scale of the heterogeneity. The detailed theoretical development of the modeling of the mean flow behavior and the variation about the predicted mean involved in the study of subsurface-flow in heterogeneous aquifers may be referred to the textbooks such as Dagan (1989), Gelhar (1993), Zhang (2002) and Rubin (2003). A stochastic methodology will be used to arrive at the results of this work. This paper begins with the development of integrated equation of the depth-averaged hydraulic head, for an essentially horizontal flow in confined aquifers of variable thickness. These provides a basis for a stochastic description of the statistics of the flow fields. General results such as the variance of the depth-averaged head, and the mean and variance of integrated discharge are obtained for the temporally and spatially varied flow fields, but closed-form expressions are limited to the case of steady unidirectional mean flow in the horizontal plane. The stochastic methodology employed to derive the closed-form expressions is based on the representation theorem and Fourier-Stieltjes representation. Focus of analysis is placed on the influence of the controlling parameters on the variability of fluid field. The application of the present stochastic theory is demonstrated by predicting the total specific discharge under uncertainty from the field data. To our knowledge, the influence of aquifer thickness on the variability of flow field has not been investigated in the literature so far because of the difficulty of using the traditional approach. The stochastic approach presented in this study provides an efficient way to analyze flow fields in confined aquifers with non-uniform thickness. It is our hope that the present study will stimulate further research in this area.

Mathematical framework
The present research is concerned with a confined aquifer of variable thickness bounded by two confining beds. The governing equation of transient flow in confined aquifers in the absence of recharge can be written in the form (e.g., Bear 1979;de Marsily 1986) as where S s is the specific storage coefficient of the aquifer, K is the hydraulic conductivity (tensor), h = h(x,t) is the hydraulic head, and x (= (x 1 ,x 2 ,x 3 )) is the spatial coordinate vector. Equation (1) is obtained from coupling the mass conservation equation with Darcy's law. For an isotropic, but inhomogeneous medium, Eq. (1) becomes For later applications, it is convenient to choose one of the principal coordinates, x 3 , perpendicular to the confining beds and the other two principal coordinates, x 1 and x 2 , parallel to the plane of the confining beds. For flow in the x 1 -x 2 plane (horizontal flow), integration of Eq.
(2) along the x 3 -axis and use of Leibniz' rule leads to the following general structure: where b 1 (x 1 ,x 2 ) and b 2 (x 1 ,x 2 ) describe the surfaces of the fixed bottom and ceiling from the horizontal plane of reference, and The last two terms in Eq. (3a) represent the leakage fluxes exchanged between the confined aquifer and its confining beds in the direction of x 3 -axis. In the development of Eq. (3a), all the terms involved in the fluxes at the boundaries are removed from Eq. (3a) on the basis of the no-slip condition at the boundary. Equation (3a) is a depthintegrated continuity equation for an essentially horizontal flow in a confined aquifer of variable thickness. By introducing the depth-averaged hydraulic head operator, and the assumption of K = K(x 1 ,x 2 ), Eq. (3a) may be written as where S(x 1 ,x 2 ) is the storage coefficient (or storativity) of the aquifer (= S s B(x 1 ,x 2 )), B(x 1 ,x 2 ) = b 2 (x 1 ,x 2 )b 1 (x 1 ,x 2 ) (an aquifer's thickness) and T(x 1 ,x 2 ) is the transmissivity of the aquifer (= K(x 1 ,x 2 )B(x 1 ,x 2 )), interpreted as the depthintegrated hydraulic conductivity. Furthermore, introducing the assumption of vertical equipotentials (i.e., h x 1 ; x 2 ; ð b 2 ; tÞ %hðx 1 ; x 2 ; tÞ % h x 1 ; x 2 ; b 1 ; t ð Þ ; Bear 1979; Bear and Cheng 2010) into Eq. (4b) yields Equation (5) indicates that the essentially horizontal flow is suitably described in terms of the average hydraulic head, taken over the thickness of the aquifer. The two aquifer properties characterizing the medium are the depth-integrated hydraulic conductivity and specific storage coefficient, namely transmissivity and storativity, respectively. Alternatively, in the absence of the leakage fluxes Eq. (5) may be rewritten as which is convenient for the stochastic analysis of the effect of the aquifer's thickness. It is worth mentioning that the assumption of vertical equipotentials or, equivalently, the assumption of essential two-dimensional flow parallel to the confining beds, is valid when the variation of the aquifer's thickness is much smaller than the average thickness (Bear 1979;Bear and Cheng 2010). The error introduced by this assumption is very small in most cases of practical interest, which greatly simplifies the analysis of flow in confined aquifers. Note that using the essentially horizontal flow approximation introduced above, the diffusion equation in phreatic aquifers can be obtained by integrating (or averaging) the threedimensional flow model along the vertical coordinate axis (e.g., Bear, 1979;de Marsily, 1986).
Equation (6) serves as the starting point for the present study, where the spatial variation of hydraulic conductivities in the continuum sense is treated as a random field characterized by the theoretical spatial covariance function. The random heterogeneity in hydraulic conductivity as an input parameter in Eq. (6) leads to the random heterogeneity in depth-averaged hydraulic head, the output of Eq. (6). Consequently, Eq. (6) is a stochastic differential equation with a stochastic outputhðx 1 ; x 2 ; tÞ. In order to obtain a complete description of model output, one would have to specify its joint probability density function over the entire flow field. A more realistic approach, however, is to represent the resulting predictions in terms of certain statistical moments of the random output field and this is the task to be performed here. The mean or expected value (the first-moment) represents the optimum unbiased predictor of the output. The variance-covariance (the secondmoment) of the output field serves to quantify predictive uncertainty (variability).

General derivation of variances of the head and integrated discharge
Representing each of random variables in Eq. (6), namely lnK andh, in terms of a mean and a zero-mean perturbation about the mean, results in the first-order stochastic partial differential equations for the mean depth-averaged head and the depthaveraged head perturbation, respectively, as In the derivation of Eqs. (9) and (10), F (the mean lnK) is treated as constant. Note that since small changes in lnK may correspond to large variations in K, there is an obvious advantage to work in terms of lnK rather than K. Unlike the classic stochastic equations for the confined aquifer flow, the variable thickness of confined aquifer introduces the second terms on the right-hand sides of Eqs. (9) and (10), respectively. Equations (9) and (10) form a system of coupled partial differential equations in the sense that one must know the behavior of the gradient of the mean depth-averaged hydraulic head in order to solve Eq. (10). Upon the solution of Eq. (9), Eq. (10) may be solved using Fourier-Stieltjes representations for the perturbed quantities in Eq. (10). From the functional spectral relationship between the fluctuations in lnK andh (the solution of Eq. (10)), it is then possible to determine the variance of the depth-averaged hydraulic head within the following framework for known autocorrelation function of the log-hydraulic conductivity: Many applications require information such as the fluxes as a prelude to more complex calculations. In addition, since solute transport is controlled by the underlying velocity field, the statistical moments of the velocity field are essential for studying solute transport process at the field scale. Therefore, the characteristics of specific discharge variability are important in the analysis of groundwater related problems. Similarly, applying Leibnitz' rule to the equation for the specific discharge through the confined aquifer of variable thickness under the aforementioned assumption of vertical equipotentials yields the vertically integrated specific discharge (total specific discharge) in the x i direction as The reader is referred to Bear and Cheng (2010) for a detailed derivation of Eq. (12). Due to the random hydraulic conductivity as an input parameter, Eq. (12) is the stochastic differential equation governing the random vertically integrated specific discharge through the confined aquifer of variable thickness. Following the approach of Gelhar (1993), the mean flow characteristics of a heterogeneous confined aquifer can be evaluated by substituting Eqs. (7) and (8) into Eq. (12), taking the expected value of it, treating F as a constant, and disregarding the third-order perturbation terms where Q x i is the mean integrated discharge, r f 2 is the variance of log-conductivity, J i = -qH/qx i , and The variations of the integrated discharge about the mean can be found simply by subtracting the mean expression of Eq. (13) from Eq. (12) where terms involving the products of perturbed quantities have been neglected in the final expression. The variance of the integrated discharge may be determined by employing the representation theorem for the discharge perturbation as follows Natural aquifer materials are heterogeneous in terms of their flow properties such as the hydraulic conductivity so that the local water flow rate can vary in different regions of the aquifer. The effect of heterogeneity can be accounted for by formulating the specific discharge problem in a stochastic framework. The mean specific discharge in Eq. (13) gives us an unbiased estimate of a system state in a heterogeneous medium, and the corresponding variance in Eq. (16) characterizes the uncertainty associated with the estimation. These two moments construct confidence intervals for the discharge field. The following section is devoted to derive closed-form analytic expressions for the variance of the depth-averaged hydraulic head, the mean and variance of integrated specific discharge in Eqs. (11), (13) and (16), respectively, based on the Fourier-Stieltjes representation approach.

Results and discussion
To gain a clear insight into the impact of the aquifer's thickness, the following theoretical development restricts attention to the case of steady unidirectional mean flow (e. g., Naff and Vecchia 1986;Gelhar 1993;Rubin and Bellin 1994;Ni and Li 2006) in the x 1 direction through a confined aquifer whose thickness varies as B(x 1 ) = b 0 exp(-ax 1 ) (Hantush 1962b;Marino and Luthin 1982) where b 0 and a are positive geometrical parameters. Note that three types of expressions have been used in the literature to describe the thickness of the aquifer, such as linear (Hantush 1962a;Zamrsky et al., 2018), quadratic (Cuello et al. 2017;Zamrsky et al. 2018), and exponential (Hantush, 1962b;Marino and Luthin 1982) varying thickness expressions. As such, the two-dimensional aspect of the flow is entirely due to local deviations of the gradient from its mean. In practice, the flow is generally unsteady, but after a prolonged period of time, the head changes very slowly in time, and a steady state is practically reached.
Under the above assumptions, Eqs. (9) and (10) can be simplified, respectively, as Equation (17) admits a representation of the mean depthaveraged head gradient of the form where J 0 is a known value of J(x 1 ) at the reference location x 1 = 0. Furthermore, this work treats lnK as a weakly stationary (second-order stationary) random process and considers the autocovariance of fluctuations in lnK field to be (Whittle 1954;Rodriguez-Iturbe and Mejia 1974;Gelhar 1977) Cf where n is the magnitude of the separation, b = p/(2k), k is the correlation scale of lnK field, and K 1 (-) is the first-order modified Bessel function of second kind. Whittle (1954) stated that the covariance function in Eq. (20) is more suitable than the exponential form for the analysis of twodimensional random processes. The spectrum associated with Eq. (20) is where the R i are the components of the wave number vector R (= (R 1 , R 2 )). For every second-order stationary process f, there is a complex-valued, random distribution Z f (R) in the wavenumber domain such that (e.g., Lumley and Panofsky 1964;Christakos 1992) where dZ f (R) represents an orthogonal increment. Note that when the process f is non-stationary, the presentation (22) is no longer valid. On the other hand, with no restriction on stationarity in a random process, the process h 0 may be represented in the form (e.g., Priestley 1965;Ni et al. 2010Ni et al. , 2011 where K(-) is an oscillatory function in the sense of Priestley's definition (Priestley, 1965). In fact, Eq. (22) is simply a special case of Eq. (23) with the oscillatory function chosen as K(x 1 ,x 2 ;R 1 ,R 2 ) = exp(x 1 R 1 ? x 2 R 2 ). The use of the representation (23) for the process h 0 is due to the effect of variable thickness of the confined aquifer, causing the nonstationarity in the depth-averaged hydraulic head field. Note that nonuniformity in the mean head gradient caused by the effects of the physical flow boundaries (e.g., Oliver and Christakos 1995), uniformly distributed recharge (e.g., Rubin and Bellin 1994) and trends in mean log hydraulic conductivity fields (e.g., Ni and Li 2005) can also lead to nonstationarity in the statistics of flow fields in heterogeneous aquifer. Substituting Eqs. (19), (22) and (23) which admits the solution given by The solution applies only to a region of interest sufficiently far from the boundary so that its effect on head fluctuations is not felt. Combining Eq. (25) with Eq. (23) gives

Variance of depth-averaged hydraulic head
Application of the spectral representation theorem (Eq. (11)) for h 0 in Eq. (26) yields Finally, the integral in Eq. (27) with respect to R is evaluated along with Eq. (21) to yield the closed-form expression for the variance of depth-averaged hydraulic head as where g = p/(4ak). It is apparent from Eqs. (19) and (28) that the statistics of random depth-averaged head fields is nonstationary. That is, the nonuniform thickness of the confined aquifer introduces nonuniformity in the mean depth-averaged head gradient and, in term, nonstationarity in the statistics of the head fields. Figure 1 displays dimensionless results of the variance of depth-averaged hydraulic head in Eq. (28) versus the aparameter for various values of the positions. It can be clearly seen that the larger the a parameter, the greater the variability of depth-averaged hydraulic head. Notice from Eq. (18) that the variation of fluctuations in lnK appears as a forcing term that derives the variation of fluctuations in depth-averaged head. The behavior of an increase in the variability of head with the a-parameter may therefore be explained by the analysis of the input-output spectrum relation. From Eqs. (11) and (23), an evolutionary power spectral density of a non-stationary process, S hh (-), can be defined as, according to (Priestley 1965) which reveals from Eq. (25) that Shhð x 1 ; R1; R2Þ ¼ Kðx 1 ; x 2 ; R1; R2Þ j j 2 Sff ðR 1 ; R2Þ Since the head variance may be interpreted as a measure of the total energy of the nonstationary process at the position x 1 , Eq. (29) gives a decomposition of total energy in which the contribution from wave numbers R 1 , R 2 is S hh (x 1 ;R 1 ,-R 2 )dR 1 dR 2 . The result in Fig. 2, a plot of the dimensionless oscillatory function in Eq. (30) versus the dimensionless wave number, shows that the contribution of the local spectrum of lnK processes to that of the head processes is enhanced by a larger a-parameter, which, in turn, leads to a higher depth-averaged head variability. It is worth noting that Kðx 1 ; x 2 ; R1; R2Þ j j 2 acts as a low-pass filter in the wave number domain which attenuates the high wave number portion of the spectrum. The filtering associated with the flow process smooths-out much of the small-scale variations caused by the input lnK parameter. Physically, this feature implies that the head field is much smoother than the lnK field. The spatial correlation of the head perturbations is larger than that of the conductivity perturbations.
In addition, using Eq. (22) along with the fluctuations in depth-averaged head given by Eq. (26), the cross-correlation between the lnK perturbation and the perturbation in depth-averaged head can be obtained as The result of Eq. (31) is presented graphically in Fig. 3, which indicates an increase in the f-h 0 cross-correlation with the a-parameter. The f-h 0 cross-correlation may be interpreted as the projection of amplitude of log-conductivity fluctuations onto that of head fluctuations. Since the spatial variation of hydraulic conductivities is the source of the depth-averaged head variation, a larger f-h 0 cross-correlation increases the correlation of head perturbations. This provides an alternative way to explain the behavior shown in Fig. 1. Another aspect of interest is the impact of the correlation scale of lnK field on the variability of depth-averaged hydraulic head. An evaluation of the effect of the correlation scale of lnK field is shown in Fig. 4. As expected, a larger correlation scale causes greater variability in the depth-averaged head. Stochastic processes with a larger correlation scale tend to show an increase in values after previous increases, and a decrease after previous decreases. The data profile with the larger correlation scale exhibits quite clear trends with relatively little noise, i.e. a fewer number of mean value crossings. As such, a larger lnK correlation scale means larger inclusions, which, in

Mean and variance of integrated specific discharge
In the case of steady unidirectional mean flow, the mean and variance of integrated specific discharge (total discharge) in Eqs. (13) and (16), respectively, take the form The terms, E fðx 1 ; (32) and (33) (35), are expressed, respectively, by making use of the representation theorem as follows Equations (36) where Figures 5 through 7 show the effect of the a-parameter on the mean integrated specific discharge and variability of integrated discharge in the directions of x 1 and x 2 , respectively. These figures clearly indicate that introduction of a larger a-parameter leads to greater mean and variability of integrated specific discharge. As shown in previous sections, nonuniformity in the mean flow field introduced by the spatial variable thickness of the confined aquifer leads to nonstationarity in statistics of the head fields. Consequently, the statistics of the integrated specific discharge fields is nonstationary (position-dependent). Due to this dependence, an effective transmissivity throughout the flow domain does not exist following the rigorous definition of effective properties (e.g., Dagan 1989), representing the aquifer's material properties. That is, the heterogeneous material cannot be replaced by a homogeneous one of an effective property. When the integrated discharge divided by the mean head gradient varies spatially, it is sometimes referred to as the pseudo-effective parameter (e. g., Sanchez-Vila 1997), apparent transmissivity (e. g., Dagan 2001) or equivalent transmissivity (e. g., Rubin 2003).
By the representation theorem for the perturbation of integrated discharge Q 0 x i , the variance of integrated discharge in the i-direction can be written in the following simple form as where S QiQi ( -) is the evolutionary power spectral density of the integrated discharge in the i-direction. Combining this with Eqs. (16), (22) and (23) yields SQ 2 Q 2 ðx 1 ; R1; R2Þ ¼ K Q 2 ðx 1 ; R1; R2Þ j j 2 Sff ðR 1 ; R2Þ where R 2 = R 1 2 ? R 2 2 . The incremental S QiQi (x 1 ;R 1 ,R 2 )-dR 1 dR 2 in Eq. (44) describes a wave number decomposition of the total energy of the process Q 0 x i . It can be demonstrated from the plots of Eqs. (45) and (46) in the wave number domain that the terms K Q 1 ðx 1 ; R1; R2Þ j j 2 and K Q 2 ðx 1 ; R1; R2Þ j j 2 in Eqs. (45) and (46), respectively, grow with the a-parameter. This means that the contribution to the spectrum of the integrated discharge from the space variation of the local spectrum of lnK processes grow with the a-parameter, resulting in a higher variability of the integrated discharge for a larger a-parameter (see Figs. 6 and 7). In addition, similar to the filtering role of Kðx 1 ; x 2 ; R1; R2Þ j j 2 in Eq. (30), K Q 1 ðx 1 ; R1; R2Þ j j 2 and K Q 2 ðx 1 ; R1; R2Þ j j 2 attenuate the high wave number portions of the spectrum in the x 1 and x 2 directions, respectively. The integrated discharge perturbations are correlated over a larger distance than the conductivity perturbations.

The case of increasing thickness of the confined aquifer
The analysis leading to the closed-form solutions is limited to the case of flow in a confined aquifer whose thickness decreases in the x 1 -direction. In the case of increasing thickness in the x 1 -direction (B(x 1 ) = b 0 exp(ax 1 ), where b 0 and a are positive constants), the resultant expressions for the depth-averaged head, and mean and variance of integrated specific discharge are similar to those obtained when the aquifer's thickness decreases in the x 1 -direction, with the exception that the term exp(ax 1 ) is replaced by exp(ax 1 ) in Eqs. (28), (40), (42) and (43), namely ð49Þ These equations suggest that the a-parameter plays a role in reducing the variance of the depth-averaged head, and mean and variance of the integrated specific discharge in the case of an increasing aquifer's thickness in the x 1direction.

Differences between traditional and current approaches
It is important to recognize some significant differences between traditional and current approaches to describing flow in two-dimensional confined aquifers. The traditional approach to regional groundwater flow problems assumes h(x 1 ,x 2 ,x 3 ) & h(x 1 ,x 2 ) and introduces the depth-integrated hydraulic conductivity operator to map groundwater flow equation from three dimensions to two. This means that the effects of the variation of K in x 3 -direction and the aquifer thickness are implicit in the term T(x 1 ,x 2 ). As such, the diffusion equation can be approximated as Since Eq. (52), along with Eq. (51), is originally designed for the analysis of flow field in a confined aquifer with uniform thickness, it is very difficult to assess the influence of thickness on the flow field with them. An example that outlines the evaluation of the head variance using Eq. (52) is given below. Spectral analysis of two-dimensional steady-state groundwater flow in a confined aquifer was performed by Mizell et al. (1982). Under the assumption of a unidirectional mean flow, the head variance is given by (Mizell et al. 1982) where J 1 = -dE[h]/dx 1 , S yy (-) is the spectrum of fluctuations in y and y = lnT. When the hydraulic conductivity field is a second-order stationary process, the covariance function for the transmissivity field defined by Eq. (51) is of the form (Vanmarcke 1983) CT ðn 1 ; n 2 Þ ¼ 2 where n (= (n 1 , n 2 , n 3 )) is the separation vector, B is the depth of the aquifer and C K (-) is the covariance function of K. Following the approach of Gutjahr et al. (1978), the relationship between the covariance function of lnT and that of T can be determined as where Ç = E[lnT], and r T 2 is the variance of T. The spectrum S yy of fluctuations in lnT can be obtained by taking the inverse Fourier transform of C y in Eq. (55a), i.e., Finally, the head variance considering the variable aquifer thickness can be determined by substituting Eqs. (54)-(55) into Eq. (56) for the Fourier transform inversion and after the inversion, substituting the result into Eq. (53) and then integrating over the wavenumber domain.
On the other hand, this work uses the hydraulic approach (Bear 1979;Bear and Cheng 2010), namely the assumption of K(x 1 ,x 2 ,x 3 ) & K(x 1 ,x 2 ) and introduction of a depth-averaged head operator (Eq. (4a)), to reduce the three-dimensional diffusion equation to two-dimensional one: which is the rewrite of Eq. (5) in the absence of the leakage fluxes. The effect of variation of K in x 3 -direction is implicitly reflected in the depth-averaged head term hðx 1 ; x 2 ; tÞ. Equation (57) provides an efficient way for the analysis of flow fields in confined aquifers of non-uniform thickness.
In addition, the common observations of flow in porous media are hydraulic head measurements from wells screened over extended sections of the medium. The measurement at a given location represents approximately a depth-averaged actual hydraulic head, resulting from flow through a three-dimensional hydraulic conductivity field, over the thickness of the medium. This means that the depth-averaged head representation used in this work is consistent with what is observed from the fields. The problem solved in this work is practical.

Application in the prediction of total discharge considering uncertainty
Many practical applications of groundwater modeling involve predictions of total flow discharge over a relative large space scale, where only limited field data are available. There will be a great deal of uncertainty in the application of the solution of the mean model (or the deterministic model) as an estimate in field situations. An important feature of the stochastic approach described here is its essentially predictive character. The variance Eq. (42) (or Eq. (49)) can be regarded as a quantitative measure of the uncertainty associated with the prediction in field situations using the mean model. As such, Eq. (40) (or Eq. (48)) ± twice the square root of Eq. (42) (or Eq. (49)) (standard deviations) offers a useful basis for predicting the total discharge under uncertainty. Below is an example to illustrate how this can be done using the field data.
To apply the present stochastic model to a particular area or region, such as the Patuxent aquifer system in Anne Arundel county, Maryland, the specific parameters such as the geometrical parameters describing the thickness of the confined aquifer (b 0 anda) and the correlation scale of lnK field (k) must be estimated. For the analyses that followed, the estimation of these parameters is based on the field data presented in Andreasen et al. (2013) and Andreasen (2017). Figure 8a, b show the variation in aquifer thickness and the distribution of hydraulic conductivity as a function of distance from the aquifer outlet in the x 1 direction, respectively. The data points in Fig. 8a are interpolated from Fig. 4 of Andreasen et al. (2013). The curve with the best fit to the data points is given by The data points in Fig. 8b are obtained from the interpolated transmissivity from Fig. 8 of Andreasen (2017) divided by the best-fit B (Eq. (58)). If the number of samples for the hydraulic conductivity measurements is not large enough, the variogram, an alternative description of the structure of the random hydraulic conductivity field, can be used to estimate the variance and correlation scale of the lnK field. For the case where the lnK process is second-order stationary, the covariance function and the variogram are linearly related by (Mathero 1971;de Marsily 1986) where c f (-) is the variogram of lnK field. For the following analyses, the traditional experimental variogram estimator is used, defined as follows (de Marsily 1986; Kovitz and Christakos 2004): where lnK i and lnK j are two measurements separated by a distance (separation) n and N is the number of pairs of data points in a selected separation interval. The points shown in Fig. 9 represent the experimental variogram of lnK in the x 1 -direction calculated from Eq. (60) and the hydraulic conductivity measurements described in Fig. 8b. From the fit of the theoretical variogram (i.e., the combination of Eq. (59) and Eq. (20)) to the experimental variogram Eq. (60), the variance and correlation scale of lnK field are determined as r f 2 = 0.000173 and k = 128 m, respectively. As shown in Fig. 9, the line that best fits the experimental variogram is used to estimate the variance and correlation scale of the log-conductivity field.
Finally, based on the specific parameters identified above, the predicted mean longitudinal total discharge (Eq. (48)) and associated standard deviation (square root of Eq. (49)) are presented in Fig. 10, where the mean and standard deviation of total discharge are normalized by b 0 e F J 0 . The general stochastic framework outlined here has potential applications in predictions over relatively large spatial scales where direct observations of such a dependent variable are not feasible.
It is now generally accepted that the transport of solutes in the subsurface is controlled by spatial variability in flow fields (e.g., Cvetkovic et al., 1991;Fiori et al., 2003;Darvini, 2013). Stochastic environmental risk assessment considers the effects of uncertainty (variability) in flow Fig. 10 Predicted normalized mean longitudinal total discharge profiles along with two standard deviation intervals as a function of dimensionless distance fields on contaminated groundwater (e.g., Benekos et al., 2006;Mumford et al., 2015;Henri et al., 2020). The results of this work for flows in random, nonstationary flow fields should provide useful information for stochastic environmental risk assessment.

Conclusions
Taking into account the influence of the thickness, the integrated equations of the depth-averaged hydraulic head and integrated specific discharge are generated for flow through confined aquifers of nonuniform thickness. They are based on the integration of the continuity equation and equation for the specific discharge over the thickness, respectively. The application of stochastic methodology for solving the integrated equations leads to the general results for the statistics of the flow field in the Fourier domain, such as the variance of the depth-averaged head, and the mean and variance of integrated discharge. The closedform expressions are obtained for the case of steady unidirectional mean flow in the horizontal plane, which, however, allow to evaluate the influence of the aquifer thickness on the variability of the flow field. The application of the present stochastic theory in predicting the total specific discharge under uncertainty in the field is also demonstrated.
Equation (6) developed in this work provides an efficient way to evaluate the effects of nonuniform aquifer thickness on flow fields, whereas it is very difficult to do so with a conventional diffusion equation for regional groundwater flow through confined aquifers. This paper also proposes an analytical methodology for quantifying the variability of nonstationary flow fields using the Fourier-Stieltjes representation approach. The variance can be viewed as a quantitative measure of the uncertainty that can be expected when the mean model is applied. This work provides a practical way to calculate the mean discharge ± twice the square root of the variance to solve the problem of regional groundwater discharge prediction under uncertainty on relatively large spatial scales where direct measurements are not possible.
The analysis of our closed-form solutions reveals that in the case of aquifers with decreasing thickness the a-parameter has a positive impact on the variability of the flow field. The a-parameter acts as a low-pass filter which smooths out much of the small-scale variations caused by the input lnK parameter. Therefore, the larger the a-parameter, the greater the variability of the depth-averaged hydraulic head and integrated specific discharge. A larger lnK correlation scale represents greater consistency of head variation around the mean, which means fewer mean crossings for the depth-averaged head process around the mean head value, resulting in a larger variability of the depth-averaged head field. In addition, an introduction of the a-parameter leads to nonstationarity in the statistics of the head field and thus to a position-dependent mean integrated specific discharge. On the other hand, in aquifers with increasing thickness the a-parameter plays a role in reducing the variability of the flow field.