Stochastic analysis of field-scale heat advection in heterogeneous aquifers

Owing to the analogy between the solute and heat transport processes, it can be expected that the rate of growth of the spatial second moments of the heat flux in a heterogeneous aquifer over relatively large space scales is greater than that predicted by applying the classical heat transport model. The motivation of stochastic analysis of heat transport at the field scale is therefore to quantify the enhanced growth of the field-scale second moments caused by the spatially varying specific discharge field. Within the framework of stochastic theory, an effective advection-dispersion equation containing effective parameters (namely, the macrodispersion coefficients) is developed to model the mean temperature field. The rate of growth of the field-scale spatial second moments of the mean temperature field in the principal coordinate directions is described by the macrodispersion coefficient. The variance of the temperature field is also developed to characterize the reliability to be anticipated in applying the mean heat transport model. It is found that the heterogeneity of the medium and the correlation length of the log hydraulic conductivity are important in enhancing the field-scale heat advection, while the effective thermal conductivity plays the role in reducing the field-scale heat advection.


Introduction
The temperature of the land surface is influenced by seasonal heating and cooling.Water seepage near the land atmosphere interface results in a heat transport that modifies the temperature profile and, in turn, affects most reactions occurring in the aquifers.In addition, the information on the heat trans-fer provides better-constrained groundwater flow and permeability estimates.It is therefore of great importance in characterizing and predicting the heat transport processes in the aquifers.Comprehensive overviews of selected work on heat are given by Anderson (2005) and Saar (2010).
The spatially varied velocity field creates the degree of spreading of a solute plume in a heterogeneous aquifer that is greater than what would occur by local dispersion alone in the uniform velocity field.Motivated by that, a stochastic methodology is devoted to relating this enhanced spreading to the characteristics of the velocity field and thus to the statistical properties of hydraulic conductivity field based on the representation of natural heterogeneity as a spatial random variable characterized by a limited number of statistical parameters.This leads to a solution in terms of an effective dispersion coefficient (macrodispersion coefficient) for describing the rate of growth of the second moments of the ensemble averaged concentration field.The stochastic methodology has successfully provided a basis framework for quantifying and understanding the effect of the natural heterogeneity on the field-scale spreading process.
The stochastic methodology is generally built around either the Eulerian or the Lagrangian framework for analyzing the solute transport in heterogeneous media.More details on the construction of the Eulerian and Lagrangian approaches and their application to the analysis of the solute transport in heterogeneous media are provided in Rubin (2003).The Eulerian approach develops an effective advection-dispersion equation (mean transport equation) containing effective parameters and seeks a quantitative measure of the uncertainty (the variance) anticipated in applying the effective transport Published by Copernicus Publications on behalf of the European Geosciences Union.
Similar to Taylor's (1921) classics analysis of turbulent diffusion, the Lagrangian analysis of field-scale solute transport is focused on the statistical properties of displacements of solute particles through a random velocity field.It offers an alternative and allows the development of preasymptotic coefficients, travel time statistics of solute particles and solute fluxes.Note that the effective dispersion coefficient is determined by half the rate of change of the particle displacement variance (or the spatial second moment of a concentration distribution).This approach has been applied to analyze the nonreactive solute transport in heterogeneous media in a number of papers (e.g.Dagan, 1984Dagan, , 1987;;Neuman and Zhang, 1990;Dagan et al., 1992;Rubin and Seong, 1994;Indelman and Rubin, 1996;Cvetkovic et al., 1996;Dagan and Fiori, 1997;Andričević, 1998;Fiori and Dagan, 1999;Destoumi et al., 2001;Riva et al., 2001;Fiorotto and Caroni, 2002;Guadagnini et al., 2003;Caroni and Fiorotto, 2005;Bellin and Tonina, 2007).
Due to the analogy between the contaminant and heat transports, it is expected that the heterogeneity of natural formations also plays an important role in influencing the heat advection at field scale.In other words, predictions from the classical heat transport equation in the uniform velocity field subject to a great deal of uncertainty.However, the application of stochastic methods, used to predict the field-scale solute transport, to the analysis of heat transport by groundwater in heterogeneous aquifers has so far not been attempted, and this is the task undertaken here.This task is performed using a spectral approach (e.g.Gelhar and Axness, 1983) based on Fourier-Stieltjes representations for the perturbed quantities under the assumption of local statistical homogeneity.We seek a mean heat transport containing macrodispersion coefficients for describing the fieldscale heat advection and the variance for quantifying the uncertainty anticipated in applying the mean heat transport equation.It is hoped that our findings will provide a basic framework for understanding and quantifying field-scale heat transport processes and be useful in stimulating further research in this area.

Mathematical formulation of the problem
Owing to the analogy between the solute and heat transport processes, the governing equations for transport in the aquifers can be represented by similar advection-dispersion equations.Following de Marsily (1986), a temperaturebased advection-dispersion equation at the local scale is of the form (1) In Eq. ( 1), T is the temperature, K T is the effective thermal conductivity, D m is mechanical dispersion coefficient for heat transport, ρ w and C w are density and specific heat of fluid, respectively, q i is the specific discharge in the principal coordinate directions, and ρ and C are density and specific heat of rock-fluid matrix, respectively.The effect of thermal dispersion is very small and negligible when compared with that of conduction (Bear, 1972;Hopmans et al., 2002).This simplifies Eq. (1) to where µ = K T /ρC and ν = (ρ w C w )/(ρC).Note that K T , ρ w , C w , ρ, and C in Eq. ( 2) are treated as constants since their variabilities are usually smaller than the variability in hydraulic conductivity (Anderson, 2005).It has been concluded from the analytical model of cold water front movement in a geothermal reservoir by Stopa and Wojnarowski (2006) that the velocity of the thermal front found from the weak model solution differs from the velocity obtained under the assumption of constant thermal properties by about 1 to 14 %, depending on the temperatures used for evaluation of the thermal properties.In addition, Lo Russo and Taddia (2010) mentioned that the density and viscosity variations with temperature can be considered negligible for systems in an unconfined aquifer with temperature changes below 10-15 K. On the other hand, the heat transport simulation should take account of the physical temperature dependencies of the thermal parameters for systems with higher temperature changes ( 10 K).
In the analysis that follows, the log-hydraulic conductivity field (lnK) is assumed second-order stationary and it is characterized by its variance and correlation scale.Note that through the Darcy's law the specific discharge and the hydraulic conductivity are directly related.Thus, spatially correlated random heterogeneity in the lnK field is the cause of the stochastic specific discharge, which in turn results in spatially correlated random perturbations in the temperature field.
The random space fields of specific discharge and temperature are typically represented, respectively, by the sum of a mean and a small zero-mean perturbation as: Hydrol.Earth Syst.Sci., 16, 641-648, 2012 www.hydrol-earth-syst-sci.net/16/641/2012/ The equation for mean temperature is found by substituting perturbation expansions in specific discharge and temperature, Eqs. ( 3) and ( 4), into Eq.( 2), assuming incompressibility of the fluid and taking the ensemble average, while the equation for the temperature perturbation is obtained by subtracting the resulting mean equation from Eq. ( 2) where the mean fluid flow is parallel to the X 1 coordinate axis so that q 1 = q and q 2 = q 3 = 0, and < > stands for the ensemble average.The last term on the left-hand side of Eq. ( 5) is referred to as macrodispersive flux in the work of Gelhar and Axness (1983) for the case of the solute transport in a saturated heterogeneous aquifer.It reflects the additional (and indeed dominant) heat advection produced as the result of the correlation between specific discharge and temperature fluctuations.This term may be fully characterized by solving Eq. ( 6), which describes the temperature perturbation due to the variation of specific discharge.The evaluation of the macrodispersive flux is the focus of this paper.As will be seen below, this term introduces the dispersive effect of the variability in the specific discharge on the temperature field.
In this study, the flow domain under consideration is assumed to be of a sufficiently large extent.Note that to fully characterize the variation in flow field, one must know the spatial behavior of the mean hydraulic head.The spectral representation theorem of random head perturbations in Fourier space will not be feasible if the mean hydraulic head has the pattern of spatial variability.In other words, the critical condition needed in the solution of the flow perturbations is that the mean hydraulic gradient must be approximately constant (or the mean fluid flow is uniform).The uniform mean flow condition corresponds to the case where the size of the flow domain becomes infinite.In practice, the validity of the uniform mean flow assumption (or infinite-domain assumption) requires that the correlation length of the random fields is much smaller than the domain size (e.g.Ababou et al., 1988;Dagan, 1989).

Spectral solutions of macrdispersion coefficients
In analyzing changes in the temperature field in time, it is convenient to introduce a moving coordinate system, ξ 1 = X 1νqt, ξ 2 = X 2 , and ξ 3 = X 3 , that follows the mean ad-vective movement of the heat distribution.As such, Eqs. ( 5) and ( 6) take the following forms Equation ( 8) can be solved using the spectral representation (e.g.Gelhar and Axness, 1983;Rehfeldt and Gelhar, 1992) based on Fourier-Stieltjes representation for the perturbed quantities in wave number domain.By using this approach, the random perturbations in Eq. ( 8) are represented by the following 3-D wave number integrals: where R = (R 1 , R 2 , R 3 ) is the wave number vector, and dZ T (R, t) and dZ qi (R) are the complex Fourier-Stieltjes increments.The transient-state spectral relation follows from Eq. ( 8) through the application of Eqs. ( 9) and ( 10) and the use of uniqueness of the representations: 3 and ε = µ/q.Note that Eq. ( 11) has been obtained under the assumption of negligible perturbation-boundary effects.Naff and Vecchia (1986) studied the impervious boundary effects on the head covariances for a steady three-dimensional flow in a formation of infinite horizontal extent, bounded above and below by impervious horizontal boundaries.They demonstrated that the boundary effect is largely limited to a zone near the medium boundary.Similar results were also obtained by Rubin andDagan (1988, 1989), who analyzed the effects of constant head and impervious boundary conditions on the head variation in semi-infinite aquifers.In addition, it is of interest to note how the type of boundary conditions affects the head covariance function in heterogeneous media.Bonilla and Cushman (2000) found that for a flow of constant mean head gradient, the effects of the Dirichlet boundary conditions (prescribed head) on the head covariance function is restricted to a boundary layer from three to four integral scales.However, under the same circumstances, effects of the Neumann boundary conditions (prescribed flux) may persist as far as three to eight integral scales from the boundary.Therefore, it may conclude that the assumption of negligible boundary effects is applicable, at least far enough from the boundary.The general solution to Eq. ( 11) can be found as in which dZ 0 is a constant of integration that depends on the initial temperature distribution.In general, the mean temperature field is a smooth function of space and time and the perturbations fluctuate on a much smaller scale than that associated with variations in the mean.The mean temperature gradient in Eq. ( 12) may therefore be approximated as a constant and it may be taken outside of the time integration.It is expected that the assumption of a constant mean temperature gradient will not be valid near the heat source where a large temperature gradient and sharp curvature occur.
It is assumed that the initial temperature distribution in the aquifer is known.Thus, there is no temperature perturbation at t = 0.That is, with dZ T = 0 at t = 0.The solution of Eq. ( 12) with dZ where G i = −∂ T/∂ξ i .The random temperature perturbation results from Eqs. ( 9) and ( 13) as follows: The macrodispersive heat flux term in Eq. ( 7) (or Eq. 5) is then found by multiplying Eq. ( 14) by the complex conjugate of dZ qi and taking the ensemble average where the macrodispersivity for heat transport β ij is given by the integral The spectrum of the specific discharge in Eq. ( 16) can be determined from the Fourier-Stieltjes representation for the first-order perturbed form of the Darcy equation and the spectral solution for the head perturbations (e.g.Gelhar and Axness, 1983): where S ff (R) is the spectrum of lnK perturbations.The macrodispersivities can thus be written from Eq. ( 16) with the spectrum of the specific discharge given by Eq. ( 17) as Generally, the expression for the effective parameter (Eq.18) is a second-order tensor, however, the off-diagonal components are zero due to an odd term in R 2 or R 3 involved in the integration over the wave number domain.By substituting Eqs. ( 15) and ( 18) into Eq.( 7), we then obtain the following mean heat transport equation: where (a) for i = j , β ii is defined by Eq. ( 18), and (b) for i = j , β ij = 0.It is clear from Eq. ( 18) that the macrodispersivities, while being dependent on travel time (or mean travel distance), do not depend on the type of coordinates used.Thus, these parameters are identical in the mean heat transport equation that uses the fixed coordinates (X 1 , X 2 , X 3 ), and we can rewrite Eq. ( 19) as

Spectral solution of temperature variance
We have made use of the spectral representation and a perturbation approximation to develop an effective advectiondispersion Eq. ( 20) containing effective parameters (namely, the macrodispersion coefficients, Eq. 18) in quantifying the field-scale heat transport processes.The macrodispersion coefficients in Eqs. ( 20) or ( 18) are used to determine the fieldscale rate of growth of the spatial second moments of the heat flux in the principal coordinate directions in a heterogeneous aquifer.Since the enhanced growth of the second moments is caused by the spatially varying specific discharge field, the outcome of the spatial variation in hydraulic conductivity field, it has been related to the statistical properties of the hydraulic conductivity.This implies large uncertainty to be anticipated in applying the mean transport Eq. ( 20).Therefore, there is a need to provide a basis for judging reliability of the field-scale mean transport model.The theoretical result (Eq.14) may also be used to determine the variance of temperature fluctuations which can be developed within the spectral framework as follows (e.g.Vomvoris and Gelhar, 1990) where the asterisk denotes the operation of complex conjugation.The temperature variance is a measure of the variation of the temperature about the mean temperature.Thus, it turns out that Eq. ( 21) will give us a quantitative measure of the reliability to be anticipated in applying the mean transport model (Eq.20).
In general, the stochastic analysis leading to the analytical results (Eqs.18 and 21) relies on some kind of small parameter expansions and the assumption of stationarity for the distribution of flow properties, where the small parameter corresponds to the variability of the underlying random log hydraulic conductivity field.The validity of those assumptions requires that the standard deviation of the random log hydraulic conductivity fluctuations, σ f , should be smaller than unity (Gutjahr and Gelhar, 1981).However, from a Monte Carlo simulation study, Zhang and Winter (1999) found it to be accurate for the head moment solutions for the value of variance of log hydraulic conductivity (σ 2 f ) as high as 4.38.A similar finding was reported in Guadagnini and Neuman (1999b) even for strongly heterogeneous media with σ 2 f up to 4 from the comparison of the moments of hydraulic head with the results of numerical Monte Carlo simulations.
Based on the application of the Continuous Time Random Walk transport formalism to the fracture networks, Geiger and Emmanuel (2010) concluded that in systems with a low degree of heterogeneity, heat transfer is Fickian-like (or Fourier-like).The implication is clear that in the flow fields with a low degree of heterogeneity an effective advectiondispersion equation containing the macrodispersion coefficients can be used to model the heat transport at the field scale, as adopted in the paper.However, in systems with a high degree of spatial variability in the flow fields, such as often occurs in fractured systems, the thermal transport is expected to become anomalous like (non-Fickian heat transport) (Geiger and Emmanuel, 2010).

Closed-form expressions for macrdispersion coefficients
The field-scale dispersivity in Eq. ( 18) and temperature variance in Eq. ( 21) are determined once the form of the spectrum of log hydraulic conductivity perturbations is specified.The evaluation of Eqs. ( 18) and ( 21) cannot be performed analytically for the general case of statistically anisotropic lnK distribution.However, to take the advantage of closed-form expressions, which provide a clear insight of the impact of heterogeneity on the behavior of heat transport, we assume statistical isotropy of the lnK field.In addition, the random lnK perturbation field under consideration is characterized by a Gaussian covariance function (Dagan, 1994;Zhang and Di Federico, 1998) which has the following spectral density function where S is the separation distance or lag, σ 2 f is the variance of lnK and λ is the correlation scale of lnK.
In the limit of τ → ∞, the longitudinal and transverse macrodispersivities converge, respectively, to )} (26) The linear relationship between the longitudinal and transverse macrodispersivities and σ 2 f in Eqs. ( 24) and ( 25) suggests that the dispersive flux of the heat at the field scale increases linearly with the heterogeneity of the medium.and transverse macrodispersivities, respectively, as a function of normalized time.Similar to the case of nonreactive solute transport, the longitudinal macrodispersivity for heat transport increases monotonically with time to its asymptotic value, while the transverse macrodispersivity ?rst increases to its maximum and then decreases to its asymptotic value.These figures also indicate that macrodispersivities are inversely dependent upon the effective thermal conductivity K T at fixed travel time.The larger the effective thermal conductivity, the higher the capacity of the medium to attenuate the advection of heat front in the longitudinal and transverse directions.
The result of Eq. ( 26) is presented graphically in terms of a function of P (or the correlation scale of lnK, λ) in Fig. 2. It is clear that the correlation scale of lnK has a positive effect on the asymptotic longitudinal macrodispersivity (or the field-scale heat advection).An increase in λ produces more persistence of the specific discharge fluctuation and in turn leads to larger deviations of specific discharge from the mean specific discharge.It is clear from Eq. ( 14) that the fluctuations in the temperature field are positively correlated to those in the specific discharge.Therefore, the field-scale heat advection increases with the correlation scale of lnK.

Conclusions
Within the framework of stochastic theory, this paper has presented the analysis of field-scale heat transport in heterogeneous aquifers.Making use of the spectral representation and a perturbation approximation leads to closed-form solutions for the field-scale dispersive heat flux in the principal coordinate directions in terms of macrodispersion coefficients.These solutions, expressed in terms of the statistical properties of log hydraulic conductivity and the effective thermal conductivity, are allowed to investigate the influence of the effective thermal conductivity and correlation scale of log hydraulic conductivity on the field-scale heat advection.The general expression for the variance of the temperature field is also developed to characterize the reliability to be anticipated in applying the mean heat transport model.
Our results indicate that the heterogeneity of the medium has a positive influence on the dispersive flux of the heat at the field scale.Larger effective thermal conductivity results in reduced the field-scale dispersive heat flux and produces less heat advection.The correlation scale of log hydraulic conductivity is important in enhancing the variability of the specific discharge and in turn the field-scale heat advection.It is hoped that our findings will stimulate further research in this area.

Fig. 2 .
Fig. 2. Dimensionless longitudinal macroscopic dispersivity as a function of P where = K T /(ρ w C w q).