Uncertainty in applying the temperature time-series method to the field under heterogeneous flow conditions

Uncertainty in applying the temperature time-series method

to the field under heterogeneous flow conditions

Ching-Min Chang, Hund-Der Yeh

⇑

Institute of Environmental Engineering, National Chiao Tung University, Hsinchu, Taiwan

a r t i c l e

i n f o

Article history: Received 30 May 2014

Received in revised form 29 July 2014 Accepted 11 August 2014

Available online 19 August 2014 This manuscript was handled by Peter K. Kitanidis, Editor-in-Chief, with the assistance of Jean-Raynald de Dreuzy, Associate Editor

Keywords:

Field-scale heat transport Variance of temperature Heterogeneous aquifers

s u m m a r y

Due to the irregular distributions of aquifer hydraulic properties, the detail on the characterization of flow field cannot be anticipated. There can be a great degree of uncertainty in the prediction of heat trans-port processes anticipated in applying the traditional deterministic transtrans-port equation to field situations. This article is therefore devoted to quantification of uncertainty involving predictions over larger scales in terms of the temperature variance. A stochastic frame of reference is adopted to account for the spatial variability in hydraulic conductivity and specific discharge. Within this framework, the use of the first-order perturbation approximation and spectral representation leads to stochastic differential equations governing the mean behavior and perturbation of the temperature field in heterogeneous aquifers. It turns out that the mean equation developed in this sense is equivalent to the traditional deterministic transport equation and the temperature variance gives a measure of the prediction uncertainty from the traditional transport equation. The closed-form expression for the temperature variance developed here indicates that the controlling parameters such as the correlation scale of specific discharge, which measures the spatial persistence of the flow field, and the periodicity of the source term tend to increase the variability in temperature field in heterogeneous aquifers. The uncertainty of the traditional heat transport model increases as the penetration depth of thermal front through the aquifer increases. This suggests that prediction of temperature distribution using the traditional heat transport model in heter-ogeneous aquifers is expected to be subject to large uncertainty at a large depth (in the downstream region). For the management purpose, the variance of temperature could serve as a calibration target when applying the traditional model to field situations. It may be more reasonable to make conclusions from, say, the mean temperature with one or two standard deviations rather than only the mean temper-ature drawn from the traditional heat transport equation.

Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction

It is well known that the transport of heat in aquifers is partly driven by the flowing groundwater. Especially vertical water fluxes are prone to propagate temperature differences. The fluctuations in aquifer properties are often viewed as random processes as a result of the details of which cannot be described precisely. The spatial variations in hydraulic conductivity cause a non-uniform velocity field. Many practical problems of heat transport involve predic-tions over much larger scales than these at which direct measure-ments are possible. It can thus be expected that there can be large uncertainty in predictions of heat transport in the field based on

the traditional deterministic heat transport equation for a homog-enous porous medium. Therefore, it is useful to provide a quantita-tive measure of uncertainty, such as the variance of the predicted temperature, as a calibration target when applying the determinis-tic model to field situations. This could be performed using a sto-chastic approach.

Stochastic modeling of subsurface flow and transport recognizes hydrological properties of the porous medium to be affected by uncertainty and regards these as random. This randomness leads to predictions of the flow or transport process in terms of a relatively small number of statistical properties, such as the first and second moments of hydraulic head or concentration (namely, the mean and variance, respectively). With the introduction of statistical inference, a field-scale equation containing effective coefficients such as effective hydraulic conductivities or macrodispersivities is developed to model the ensemble mean behavior of the dependent variable. In the case of natural formations, the mean stochastic

http://dx.doi.org/10.1016/j.jhydrol.2014.08.029

0022-1694/Ó 2014 Elsevier B.V. All rights reserved.

⇑Corresponding author. Address: Institute of Environmental Engineering, National Chiao Tung University, 1001 University Road, Hsinchu 300, Taiwan. Tel.: +886 3 5731910; fax: +886 3 5725958.

E-mail address:hdyeh@mail.nctu.edu.tw(H.-D. Yeh).

Contents lists available atScienceDirect

Journal of Hydrology

solution is useful to make decisions (e.g.,Andricevic and Cvetkovic, 1996; Maxwell et al., 1999) in real life transport events, but there will be variations around the mean. Therefore, for a successful diction a quantification of the degree of variability around the pre-dicted mean behavior (the variance) should be established.

Determination of ground water flux using the analytical solu-tion to the one-dimensional heat transport model has been dem-onstrated and applied to situations of stream–aquifer interactions (e.g.,Stallman, 1965; Silliman et al., 1995; Hopmans et al., 2002; Hatch et al., 2006; Keery et al., 2007; Rau et al., 2010; Jensen and Engesgaard, 2011) and groundwater recharge (e.g., Suzuki, 1960; Taniguchi, 1993; Taniguchi and Sharma, 1993; Tabbagh et al., 1999; Bendjoudi et al., 2005; Cheviron et al., 2005). Interpretation of field observations using one-dimen-sional analytical results appropriate for a homogenous system may lead to significant errors in the predicted vertical flux in situations where the flow field is non-uniform (e.g.,Shanafield et al., 2010; Schornberg et al., 2010; Jensen and Engesgaard, 2011; Ferguson and Bense, 2011; Rau et al., 2012b; Roshan et al., 2012; Cuthbert and Mackay, 2013). In other words, the prediction can be subject to high levels of uncertainty.

As will be seen in the next section given below, the mean heat transport equation is identical to the traditional equation except that the mean specific discharge is replaced by the local specific discharge. The traditional analytical result describing the tempera-ture distribution may be interpreted as the mean of temperatempera-ture distribution, while the temperature variance may then be viewed as the uncertainty anticipated in applying the deterministic analyt-ical result. For the prediction of an actual temperature distribution in the field, it may be more reasonable to draw conclusions from the mean value (the analytical result) and the variance rather than only the mean temperature. This research is primarily concerned with the development of a quantification of deviation around the mean temperature field in a non-uniform flow field and the analy-sis of the influence of controlling parameters on that. The analyanaly-sis

we perform is relevant mainly to shallow subsurface situations that receive and transfer cyclic temperature fluctuations (i.e., daily or seasonal) over depth. The temperature fluctuations are damped with depth depending on their periodicity, so the solution gener-ally applies to the surficial zone (Anderson, 2005). We hope that the findings provided here will be useful for interpretation of field data.

2. Mathematical statement of the problem

The heat transport equation for three-dimensional saturated flow in a porous medium at the local level can be written as (e.g.,de Marsily, 1986;Demenico and Schwartz, 1998)

Ke

q

C @2T @X2 i 

q

WCW

q

C @ @Xi ðqiTÞ ¼ @T @t i ¼ 1; 2; 3 ð1Þ

where T is the temperature, Keis the effective thermal conductivity,

C and

q

are specific heat capacity and density of the fluid–rock matrix, respectively, Cwand

q

ware specific heat capacity and

den-sity of the fluid, respectively, and qiis the ith component of the

spe-cific discharge vector q = (q1, q2, q3). The effective thermal

conductivity takes into account the effects of thermal dispersion and conduction through the rock–fluid matrix. It is worth mention-ing that the effect of thermal dispersion is very small and negligible (Bear, 1972; Hopmans et al., 2002; Rau et al., 2012a). The parame-ters in Eq. (1), such as Ke, Cw, C,

q

w and

q

, are considered fixed

parameters for their variations in space and time may be assumed to be negligible (e.g., Demenico and Schwartz, 1998; Anderson, 2005).

To account for the natural heterogeneity of geological forma-tions, the log hydraulic conductivity ln K is regarded as the spatially correlated random function. Spatially correlated random heteroge-neity in ln K field results in spatial perturbations in specific dis-charge in Eq. (1)and in turn in the modeled temperature field. Nomenclature

A amplitude of temperature variations

C specific heat capacity of the fluid–rock matrix Cw specific heat capacity of the fluid

G Eq.(11)

K hydraulic conductivity

Ke effective thermal conductivity

L length of the domain

P period of temperature variations

R wave number

Sqq specific discharge spectrum

T temperature T mean temperature T0 _{fluctuation in temperature} T0* complex conjugate of T0 U =

c

q

Z vertical space coordinate

dZqZ complex random amplitude of specific discharge pro-cess

qi ith component of the specific discharge vector



qi ith component of the mean specific discharge vector

q0

i fluctuation in ith component of the specific discharge

vector q ¼ qZ t time C1 Eq.(22) C2 Eq.(23) HTq transfer function k1 Eq.(14) k2 Eq.(15) N =(

r

2_T/A2₎0.5 U1 Eq.(24) U2 Eq.(25) U3 Eq.(26) W ¼ T=A

a

e =Ke/(

q

C) b =

pa

e/(UL)

c

=

q

wCw/(

q

C)

e

Eq.(16)

g

=PU/L k correlation scale of ln K

l

1 =4

p

2

t

2+ 1

l

2 =

p

2

t

2+ 1 f =Z/L

q

density of the fluid–rock matrix

q

w density of the fluid

r

2

f variance of ln K

r

2

q variance of the specific discharge

r

2 T variance of temperature

s

=

p

2

_a

et/L2

m

=k/L

-

=exp(1/

t

)

Spatial flux variability has been discussed recently based on small-scale experimental observations byRau et al. (2012b). On a larger scale, the propagation of the temperature signal over depth in a heterogeneous streambed environment and its implications on flux estimates have been investigated numerically (e.g.,Ferguson, 2007; Schornberg et al., 2010; Ferguson and Bense, 2011).

In this study, the flow field we are concerned with is under the steady-state condition, i.e., oqi/oXi= 0. This simplifies(1)to

Ke

q

C @2T @X2i 

q

WCW

q

C qi @T @Xi¼ @T @t ð2Þ

Consider a decomposition of variables qiand T in space into a

mean and a fluctuation about the mean represented, respectively, by

T ¼ T þ T0

ð3aÞ

and

qi¼ qiþ q0i ð3bÞ

InEq. (3)the bar represents the mean value while the prime denotes the small perturbation about the mean. The perturbation is considered to be a zero-mean, spatial stochastic process. In gen-eral, it is preferable to work with perturbations which are small such that the products of perturbations are small and negligible.

Following the approach ofGelhar and Axness (1983), we substi-tute (3) into (2) and subsequently take the expectation of the resulting equation to yield the equation governing the mean temperature: Ke

q

C @2T @X2 i 

q

WCW

q

C qi @T @Xi ¼@T @t ð4Þ

In the development of Eq.(4), terms involving products of the perturbations are disregarded. The differential equation governing the perturbations of temperature, T0_{, is then obtained by}

subtract-ing the mean Eq.(4)from(2), after using(3)into(2):

Ke

q

C @2T0 @X2i 

q

WCW

q

C qi @T0 @Xi

q

WCW

q

C q 0 i @T @Xi¼ @T0 @t ð5Þ

In the present study we are interested in the case where only the mean vertical heat transport is preponderant (e.g.,Reiter, 2001). i.e., @ T=@X3 @T=@X1 and @T=@X2, and @2T=@X23 @

2_T=@X2 1 and

@2T=@X22. We also consider here the steady-state flow assumption,

where the uniform mean flow is in the vertical direction (Z-direction or X3-axis), q1¼ q2¼ 0 and q3¼ q, but perturbations to the flow are

in three dimensions. As such,(4) and (5)reduce, respectively, to

a

e @2T @Z2

c

q @T @Z¼ @T @t ð6Þ

a

e @2T0 @Z2

c

q @T0 @Z

c

q 0 Z @T @Z¼ @T0 @t ð7Þ

where

a

e= Ke/(

q

C),

c

=

q

wCw/(

q

C), and q0Z is the perturbation to the

flow in the Z-direction. Note that in the development of Eq.(7), the contribution of conduction and thermal dispersion in the transverse heat transport process is disregarded.

The mean transport Eq.(6)is identical in form to the traditional one-dimensional heat transport equation for a deterministic sys-tem if the mean specific discharge parameter in Eq.(6)is replaced with the local specific discharge parameter. The third term on the left-hand side of Eq.(7)is the sink term and reflects the dissipation produced by the mean temperature gradient interacting with the fluctuations in specific discharge. Therefore, the solution to Eq. (7), providing the relationship between the temperature and

specific discharge perturbations, forms the basis for characterizing the variability (or uncertainty) of the mean (or traditional) heat transport model. Determination of the variation of temperature field from the use of the representation theorem is the line of the research pursued here.

Note that the representation theorem applied by this work is referred to the expectation of the product of the Fourier-Stieltjes integral representation for T0 _{and its complex conjugate together}

with the orthogonality property of random Fourier increments of q0

Z. The representation theorem has been widely applied to

com-pute the variances of hydraulic head and concentration fields in the stochastic subsurface hydrology literature (e.g.,Gelhar, 1993; Zhang, 2002; Rubin, 2003).

To provide a complete description of the heat transport pro-cesses given by Eqs.(6) and (7), it is necessary to specify the initial and boundary conditions. The conditions we are concerned with are deterministic and similar to those imposed byHatch et al. (2006): TðZ; 0Þ ¼ 0 ð8aÞ Tð0; tÞ ¼ A cos 2

p

P t   ð8bÞ TðL; tÞ ¼ 0 ð8cÞ and T0 ðZ; 0Þ ¼ 0 ð9aÞ T0ð0; tÞ ¼ 0 ð9bÞ T0ðL; tÞ ¼ 0 ð9cÞ

where A and P are the amplitude and the period of temperature variations at the upper boundary, respectively (Stallman, 1965), and L denotes the maximal depth so that Z 2 [0, L]. Note that Hatch et al. (2006) reformulated Stallman’s solution (1965) to reveal the amplitude and phase features.

In the next section, we proceed to develop the analytical solu-tion of Eq.(7), which requires(6)to be solved first in order to know the mean temperature gradient.

3. Solution to the stochastic perturbation equation

The analytical solution to Eq.(6)with boundary conditions(8) can be found by using the method of eigenfunction expansions (e.g.,Farlow, 1993; Haberman, 1998) as:

TðZ;tÞ ¼ 2AL2X 1 n¼1 1 n

p

sin n

p

LZ   exp UZ 2

a

e  ₁ G exp  n2

_p

2

_a

e L2 þ U2 4

a

e ! t " # ( ð64

a

2 eL 2

_p

2_þ4

_a

2 en 2_P2_U2

_p

2_þL2_P2_U4_Þ ½L2_P2_U4_þ4

_a

2 e

p

2_ð16L2_þn2_P2_U2_Þcos 2

p

Pt   þ32

a

3 en 2_P

_p

3_sin 2

p

Pt   2AX 1 n¼1 1 n

p

sin n

p

LZ   exp UZ 2

a

e   exp  n 2

_p

2

_a

e L2 þ U2 4

a

e ! t " # þALZ L exp UZ 2

a

e   cos 2

p

Pt   ð10Þ where U =

c

q and G ¼ 16

a

4 en 4_P2

_p

4_{þ L}4 P2U4þ 64

a

2 eL 4

_p

2_{þ 8}

_a

2 eL 2 n2_P2 U2

p

2 _ð11Þ

In the large-time limit

a

et/L2 1/

p

2, we arrive at a simplified

TðZ; tÞ ¼ 2

p

A exp

p

2 n b   sinð

p

nÞ exp  1 þ1 4 1 b2  

s

  ð

K

1 1Þ þ

K

2sin 2

s

b

g

  

K

1cos 2

s

b

g

  þ Að1  nÞ exp

p

2 n b   cos 2

s

b

g

  ð12Þ

and its spatial gradient takes the form

@T @Z¼ A Lexp

p

2 n b   exp  1 þ1 4 1 b2  

s

  ð

K

1 1Þ þ

K

2sin 2

s

b

g

   

K

1cos 2

s

b

g

  ₁ bsinð

p

nÞ þ 2 cosð

p

nÞ   þ

p

2 1 bð1  nÞ  1   cos 2

s

b

g

  ð13Þ

where f = Z/L, b =

pa

e/(UL),

s

=

p

2

a

et/L2,

g

= PU/L, and

K

1¼ 64ðb2₌

_g

2_{Þ þ 4b}2_{þ 1}

e

ð14Þ

K

2¼ 32 b3

ge

ð15Þ

e

¼ 16b4_{þ 8b}2_{þ 64ðb}2₌

_g

2_{Þ þ 1 ð16Þ}

Eq.(13)puts us in a position to develop the analytical solution to Eqs.(7)and(9).

To determine the variance of temperature using the representa-tion theorem, we need to construct a wave domain solurepresenta-tion to Eqs. (7)and(9). When the specific discharge, the input parameter in Eq. (1), is defined as a second-order stationary process, its perturbed quantities can then be represented by a Fourier-Stieltjes integral representation q0 Z¼ Z1 1 exp½iRZdZqZðRÞ ð17Þ

where dZq1ðRÞ is a complex random amplitude of the process and R is the wave number. It is important to know that although the heat transport process analysis is carried out herein within a one-dimensional framework, the perturbation flow field is not a strictly one-dimensional flow. However, to simplify the analysis,Gelhar (1993) pointed out that the variability in longitudinal specific discharge can be determined within the context of such a quasi-one-dimensional treatment by using that obtained from a three-dimensional analysis.Duffy (1982)andGelhar and Gutjahr (1982) used that conceptual framework in analysis of the one-dimensional transport problem. That is why the one-dimensional representation for the longitudinal specific discharge perturbation in Eq.(17)is used.

In addition, the non-stationary Fourier-Stieltjes integral repre-sentation (e.g.,Li and McLaughlin, 1991) for the perturbed quanti-ties allows us to relate the output (T0_{) to input ðq}0

ZÞ perturbations in

the following way

T0_¼Z 1 1

H

TqdZqzðRÞ ð18Þ

where HTq represents the transfer function which describes the

relationship between spectral fluctuations in specific discharge and spatial fluctuations in temperature. Introducing(17) and (18) into(7)yields

a

e @2

H

Tq @Z2  U @

H

Tq @Z 

c

expðiRZÞ @T @Z¼ @

H

Tq @t ð19Þ

where @T=@Z is defined in Eq.(13). Transformation of the initial and boundary conditions leads(9)to

H

TqðZ; 0Þ ¼ 0 ð20aÞ

H

Tqð0; tÞ ¼ 0 ð20bÞ

H

TqðL; tÞ ¼ 0 ð20cÞ

For the case of

a

et/L2 1/

p

2, the transfer function can be

expressed as:

H

TqðZ; t; RÞ ¼ 4

p

c

A LU2

_e

exp

p

2 n b   sinð

p

nÞ  4

a

e

C

1ð

s

Þ½1  expðiRLÞ R þ iðU=

a

eÞ RðL2R2 4

p

2_Þ ( þ

C

2ð

s

Þ LU  2

a

eð1 þ exp½iRLÞ L2_R2 

p

2  i2L 2_{Uð1 þ exp½iRLÞ}   K ðL2R2

p

2_Þ2 #) ð21Þ where

C

1ð

s

Þ ¼

U

1exp  1 þ 1 4 1 b2  

s

  þ

U

2sin 2

s

g

b   

U

3cos 2

s

g

b   ð22Þ

C

2ð

s

Þ ¼ ð4b2þ 1Þ cos 2

s

g

b    exp 1 þ1 4 1 b2  

s

   þ 8b

g

sin 2

s

g

b   ð23Þ

U

1¼ ð4b2þ 1Þ

K

1þ 8 b

g

K

2þ

e

ð

K

1 1Þ

s

4b2 ð24Þ

U

2¼ ð4b2þ 1Þ

K

2 8 b

g

K

1 ð25Þ

U

3¼ ð4b2þ 1Þ

K

1þ 8 b

g

K

2 ð26Þ

Combining(21)with(18)gives

T0_{¼ 4}

_p

c

A LU2

_e

exp

p

2 n b   sinð

p

nÞ  Z1 1 4

a

eC1ð

s

Þ½1  expðiRLÞ R þ iðU=

a

eÞ RðL2_R2  4

p

2_Þ ( þC2ð

s

Þ LU  2

a

eð1 þ exp½iRLÞ L2_R2 

p

2  i2L 2_{Uð1 þ exp½iRLÞ} K ðL2_R2_

_p

2_Þ2 " #) dZqzðRÞ ð27Þ 4. Variance of temperature

It follows from the use of the representation theorem for T0_that

r

2 T¼ hT 0_T0_{i ¼ 16}

_p

2

c

2A 2 L2_U4

_e

2exp

p

n b   sin2ð

p

nÞ  Z 1 1 32

a

2 e

C

2 1 ½R2þ ðU=

a

eÞ2½1  cosðRLÞ R2ðL2R2 4

p

2_Þ2 ( þ8

a

e

C

1

C

2 LU½1  cosðRLÞ ðL2R2  4

p

2_ÞðL2_R2 

p

2_Þþ

ðLU2=

a

e 4UÞ sinðRLÞ

RðL2_R2  4

p

2_ÞðL2_R2 

p

2_Þ " 4 L 2_{UR sinðRLÞ} ðL2R2  4

p

2_ÞðL2_R2 

p

2_Þ2 # þ

C

22 L2_U2 þ 8

a

2 e 4LU

a

e ðL2R2 

p

2_Þ2 " þ4ð2

a

2 e

a

eLUÞ cosðRLÞ ðL2R2 

p

2_Þ2 þ 4 L3_U2_{R sinðRLÞ} ðL2R2 

p

2_Þ3 þ8U 2_L4_R2 ½1 þ cosðRLÞ ðL2R2

p

2_Þ4 #) SqqðRÞdR ð28Þ

where

r

2_{T is the variance of temperature, the angle bracket denotes}

the expected value operator, T0*_{is the complex conjugate of T}0_{, and}

Sqq(R) is the specific discharge spectrum in wave number domain.

Eq.(28)provides a means of quantifying the temperature variability for the mean heat transport process or the uncertainty in applying the traditional (deterministic) heat transport model.

Before evaluation of Eq.(28)can be completed, it is necessary to select the spectrum of the specific discharge process. We consider a particular form for Sqq(R) (Bakr et al., 1978; Duffy, 1982)

SqqðRÞ ¼ 2

p

K

3R2 ð1 þ

K

2R2Þ2

r

2 q ð29Þ

which is widely applicable to modeling of natural phenomena. In Eq.(29), where k and

r

2

qrepresent the correlation scale and the

var-iance of the specific discharge process.

Substituting(29)into(28)and performing integration yields-where

t

= k/L,

-

= exp(1/

t

),

l

1= 4

p

2

t

2+ 1 and

l

2=

p

2

t

2+ 1. From

a three-dimensional analysis of first-order fluctuations in flow field, the variance of the specific discharge can be expressed in the form (Gelhar and Axness, 1983; Dagan, 1987; Chang and Yeh, 2007)

r

2 q q2¼ 8 15

r

2 f ð31Þ where

r

2

f is the variance of ln K. With(31), the final result is now

given by

5. Discussion

The analytical results above are developed based on the key assumptions of smallness of the perturbations of specific dis-charge and temperature (the first-order perturbation approxima-tion), second-order stationarity of the specific discharge perturbations, and nonstationary representation for the tempera-ture perturbation. At this point it is appropriate to review those assumptions. In terms of the variability of ln K, the first-order per-turbation approximation leading to the analytical results is valid only if the variance of ln K  1 (Gutjahr and Gelhar, 1981). That is, the variance of temperature developed here based on the first-order approximation is restricted to the case of mildly heter-ogeneous media. However, the study of Monte Carlo simulations of flow through heterogeneous formations shows agreement with

the small perturbation approximation for the moments of hydrau-lic head with variance up to 4 (Zhang and Winter, 1999; Guadagnini and Neuman, 1999).

The assumption of stationarity of the specific discharge field is valid when the mean hydraulic head field is uniform (or relatively smooth). In other words, the only source of variability in specific discharge is the hydraulic conductivity perturbation field. The log-arithm of hydraulic conductivity in this work is modeled as a real-ization of a stationary random field and, in turn, stationarity of the specific discharge field is presumed. On the other hand, the space-dependent mean temperature gradient (see Eq. (13)) produces

r

2 T¼ 64A 2b2

e

2

r

2 q q2exp

p

n b   sin2ð

p

nÞ 32

C

21

t

3 1 4

t

l

2 1 

t

l

3 1 þ

p

2 b2 1 16

p

2 1

l

2 1 þ1 4 1

l

2 1 ½1 þ ð1 

t

Þ

-



t

2_{ð1 þ}

_t

_Þ

_-

_þ

_t

3_{ þ}

t

ð1 

-

Þð

t

2 1Þ

l

3 1     þ8

C

1

C

2

t

1 3

p

b 1 4 1

l

₁ 4 3

p

2ð1=

t

 2Þ

-

þ 4

t

2

-

t

 ð1 

t

Þ

-  þ1 2 1

l

2 1 4 3

p

2

-

þ 16 3

t

2 

t

ð1 

-

Þð4

t

 1Þ     þ1 4 1

l

2

t

þ ð1 

t

Þ

-

 4

t-

þ 4 3

p

2ð2  1=

t

Þ

-  þ1 2 1

l

2 2 ð

-

 1Þ

t

 2

t

2 ð2  1=

t

Þ

-

 4

t

2

-þ4 3

t

2  4 3

p

2

-  þ4

t

2

l

3 2 ½

p

2

_t

2 

-

  þ 1 6b2  1 2

t-l

1 þ

t

2_{ð1 }

_-

_Þ

l

2 1 þ1 2

t-l

2 þ

t

2_{ð1 þ}

_-

_Þ

l

2 2   þ

C

2 2

t

3 ₂1 þ

t

 ð1 

t

Þ

-l

2 2  8

t

ð1 þ

-

Þ

l

3 2 þ

p

b ð1 

t

Þ

-



t

 1

l

2 2 þ 4

t

ð1 þ

-

Þ

l

3 2    þ

p

2 b2 1 12

p

2 3

p

2

_t

_{þ 2}

_p

2_{þ 3}

l

2 2 

t

l

3 2 ½

t

þ ð2

t

 1Þ

-

 1  2

t

2

l

4 2 ½3

t

þ ð3

t

 4Þ

-

 þ 16

t

31 þ

-l

5 2   ð30Þ

r

2 T¼ 512 15 A 2

r

2 f b2

e

2 exp

p

n b   sin2ð

p

nÞ 32

C

21

t

3 1 4

t

l

2 1 

t

l

3 1 þ

p

2 b2 1 16

p

2 1

l

2 1 þ1 4 1

l

2 1 ½1 þ ð1 

t

Þ

-



t

2_{ð1 þ}

_t

_Þ

_-

_þ

_t

3_{ þ}

t

ð1 

-

Þð

t

2 1Þ

l

3 1     þ8

C

1

C

2

t

1 3

p

b 1 4 1

l

₁ 4 3

p

2ð1=

t

 2Þ

-

þ 4

t

2

_-

_

_t

_{ ð1 }

_t

_Þ

-  þ1 2 1

l

2 1 4 3

p

2

-

þ 16 3

t

2_

_t

_{ð1 }

_-

_Þð4

_t

_{ 1Þ}     þ1 4 1

l

₂

t

þ ð1 

t

Þ

-

 4

t-

þ 4 3

p

2ð2  1=

t

Þ

-  þ1 2 1

l

2 2 ð

-

 1Þ

t

 2

t

2_{ð2  1=}

_t

_Þ

_-

_{ 4}

_t

2

_-

_þ4 3

t

2_ 4 3

p

2

-  þ4

t

2

l

3 2 ½

p

2

_t

2_

_-

_  þ 1 6b2  1 2

t-l

₁þ

t

2_{ð1 }

_-

_Þ

l

2 1 þ1 2

t-l

₂þ

t

2_{ð1 þ}

_-

_Þ

l

2 2   þ

C

2 2

t

3 ₂1 þ

t

 ð1 

t

Þ

-l

2 2  8

t

ð1 þ

-

Þ

l

3 2 þ

p

b ð1 

t

Þ

-



t

 1

l

2 2 þ 4

t

ð1 þ

-

Þ

l

3 2    þ

p

2 b2 1 12

p

2 3

p

2

_t

_{þ 2}

_p

2_{þ 3}

l

2 2 

t

l

3 2 ½

t

þ ð2

t

 1Þ

-

 1  2

t

2

l

4 2 ½3

t

þ ð3

t

 4Þ

-

 þ 16

t

31 þ

-l

5 2   ð32Þ

nonstationarity of temperature perturbation process, which excludes the direct applicability of the stationary spectral repre-sentation. The nonstationary Fourier-Stieltjes integral representa-tion (Li and McLaughlin, 1991) is then used to represent the temperature perturbation process instead.

The result in Eq.(32)shows that the amplitude of the temper-ature variance is linearly proportional to the variance of ln K. The textural variations exhibited in natural porous media give rise to spatial variability of their constitutive properties. This implies that the temperature variability increases linearly with the heterogene-ity of the aquifer for mildly heterogeneous formations.

The variance of temperature as a function of correlation scale of specific discharge for various values of

g

is presented graphically in Fig. 1. Similar to the field-scale solute transport process, the corre-lation scale k has a positive influence on the temperature variance. As k increases, the persistence of correlations increases and the fluctuations spend less time around the mean. This results in a large variability in temperature field. The figure also indicates that increasing period of temperature variations P tends to increase the temperature variance with k held constant. As pointed out byGoto et al. (2005)and Wörman et al. (2012), the damping of thermal front is related to the periodicity of the source term. The thermal front with small period of temperature variations penetrates more rapidly into the aquifer than the large one does, but is dampened more abruptly with depth, which leads to a less variability in tem-perature field.

Fig. 2shows how the variability in temperature field varies with the depth. As the thermal front penetrates large regions of the aquifer, the transported heat responds to larger and larger hetero-geneities. There is a change in the size of those heterogeneities with the depth in the associated flow field that affects the move-ment of heat transport. This is why the temperature variability

increases with the depth. The increase of variability in temperature field with the depth reveals that the prediction of temperature dis-tribution is subject to large uncertainty in the far-source region (downstream region) in heterogeneous aquifers.

The analytical solution(12)to the mean temperature equation described in this work is equivalent to that to the traditional one-dimensional deterministic heat transport equation (e.g., Stallman, 1965). We can anticipate irregular variations in temper-ature around the mean in natural porous media. Therefore, the var-iance(32) gives us a quantitative measure of the uncertainty in applying the traditional transport model to field situations. The most challenging types of heat transport problems involve predic-tions over much larger scales where direct measurements are not possible. Under such conditions, the mean profile along with stan-dard deviations provides a useful way of evaluating the model pre-diction. For practical applications of heat transport modeling in the field, for example for management purposes, it may be more rea-sonable to consider, say, the mean temperature with one standard deviation (square root of Eq.(32)) rather than only the mean tem-perature drawn from the traditional heat transport equation.Fig. 3 indicates that the level of uncertainty grows with the depth and is largest in the downstream. Presented in solid line is the predicted mean temperature field, while the dashed lines present the tem-perature field corresponding to the ± one standard deviation.

6. Conclusions

The perturbation-based nonstationary spectral techniques have been applied to quantify the variability in temperature field in a heterogeneous aquifer. The closed-form expressions developed here apply to the region of shallow subsurface. We conclude from the analysis of the closed-form expression for the temperature var-iance that the correlation length scale of the specific discharge and the period of temperature variations have a strong influence on increasing the variability of temperature field. In addition, there can be large uncertainty in the prediction of temperature distribu-tion at a large depth in heterogeneous aquifers. From the practical application viewpoint, a result such as(32)could serve as a calibra-tion target when applying the tradicalibra-tional deterministic transport equation to the field situations.

Acknowledgements

This research work is supported by the Ministry of Science Technology under the Grants NSC 101-2221-E-009-105-MY2, 102-2221-E-009-072-MY2 and NSC 102-2218-E-009-013-MY3. We are grateful to the anonymous reviewers for constructive com-ments that improved the quality of the work.

Fig. 1. Dimensionless temperature variance as a function of dimensionless corre-lation length of specific discharge.

Fig. 2. Dimensionless temperature variance as a function of dimensionless depth.

Fig. 3. Dimensionless temperature profile with mean ± one standard deviation envelopes.

References

Anderson, M.P., 2005. Heat as a ground water tracer. Ground water 43 (6), 951–968.

Andricevic, R., Cvetkovic, V., 1996. Evaluation of risk from contaminants migrating by groundwater. Water Resour. Res. 32 (3), 611–621.

Bakr, A.A., Gelhar, L.W., Gutjahr, A.L., MacMillan, J.R., 1978. Stochastic analysis of spatial variability in subsurface flows: 1. Comparison of one- and three-dimensional flows. Water Resour. Res. 14 (2), 263–271.

Bear, J., 1972. Dynamics of Fluids in Porous Media. American Elsevier Publishing Company Inc, New York.

Bendjoudi, H., Cheviron, B., Guérin, R., Tabbagh, A., 2005. Determination of upward/ downward groundwater fluxes using transient variations of soil profile temperature: test of the method with experimental data. Hydrol. Process. 19 (18), 3735–3745.

Chang, C.-M., Yeh, H.-D., 2007. Large-time behavior of macrodispersion in heterogeneous trending aquifers. Water Resour. Res. 43 (11), W11501.http:// dx.doi.org/10.1029/2007WR006017.

Cheviron, B., Guérin, R., Tabbagh, A., Bendjoudi, H., 2005. Determining long-term effective groundwater recharge by analyzing vertical soil temperature profiles at meteorological stations. Water Resour. Res. 41 (9), W09501. http:// dx.doi.org/10.1029/2005WR004174.

Cuthbert, M.O., Mackay, R., 2013. Impacts of nonuniform flow on estimates of vertical streambed flux. Water Resour. Res. 49, 19–28. http://dx.doi.org/ 10.1029/2011WR011587.

Dagan, G., 1987. Theory of solute transport by groundwater. Ann. Rev. Fluid Mech. 19, 183–215.

de Marsily, G., 1986. Quantitative Hydrogeology. Academic Press, San Diego, California.

Domenico, P.A., Schwartz, F.W., 1998. Physical and Chemical Hydrogeology. John Wiley & Sons, New York.

Duffy, C.J., 1982. Stochastic modeling of spatial and temporal water quality variations in groundwater. Ph. D. thesis, New Mexico Institute of Mining and Technology, Socorro.

Farlow, S.J., 1993. Partial Differential Equations for Scientists and Engineers. Dover, New York, NY.

Ferguson, G., 2007. Heterogeneity and thermal modeling of ground water. Ground Water 45 (4), 485–490.

Ferguson, G., Bense, V., 2011. Uncertainty in 1D heat-flow analysis to estimate groundwater discharge to a stream. Ground Water 49 (3), 336–347.

Gelhar, L.W., 1993. Stochastic Subsurface Hydrology. Prentice Hall, Englewood Cliffs, New Jersey.

Gelhar, L.W., Gutjahr, A.L., 1982. Stochastic solutions of the one-dimensional convective dispersion equation, Progress report to Sandia National Laboratories on the project: A stochastic approach to risk assessment for geologic disposal of radioactive waste, Contract 20-P1-M8-0-000, Hydrology Research Program Report H-11, Geophysical Research Center, New Mexico Institute of Mining and Technology, Socorro.

Gelhar, L.W., Axness, C.L., 1983. Three-dimensional stochastic analysis of macrodispersion in aquifers. Water Resour. Res. 19 (1), 161–180.

Goto, S., Yamano, M., Kinoshita, M., 2005. Thermal response of sediment with vertical fluid flow to periodic temperature variation at the surface. J. Geophys. Res. (Solid Earth) 110, B01106.http://dx.doi.org/10.1029/2004JB003419.

Guadagnini, A., Neuman, S.P., 1999. Nonlocal and localized analyses of conditional mean steady state flow in bounded, randomly nonuniform domains: 2. Computational examples. Water Resour. Res. 35 (10), 3019–3039.

Gutjahr, A.L., Gelhar, L.W., 1981. Stochastic models of subsurface flow: infinite versus finite domains and stationarity. Water Resour. Res. 17 (2), 337–350.

Haberman, R., 1998. Elementary Applied Partial Differential Equations: With Fourier Series and Boundary Value. Prentice Hall, Upper Saddle River, New Jersey.

Hatch, C.E., Fisher, A.T., Reveaugh, J.S., Constantz, J., Ruehl, C., 2006. Quantifying surface water-groundwater interactions using time series analysis of streambed thermal records: method development. Water Resour. Res. 42, W10410.http:// dx.doi.org/10.1029/2005WR004787.

Hopmans, J.W., Simunek, J., Bristow, K.L., 2002. Indirect estimation of soil thermal properties and water flux using heat pulse probe measurements: geometry and dispersion effects. Water Resour. Res. 38 (1), 1006.http://dx.doi.org/10.1029/ 2000WR000071.

Jensen, J.K., Engesgaard, P., 2011. Nonuniform groundwater discharge across a streambed: heat as a tracer. Vadose Zone J. 10 (1), 98–109.

Keery, J., Binley, A., Crook, N., Smith, J.W., 2007. Temporal and spatial variability of groundwater–surface water fluxes: development and application of an analytical method using temperature time series. J. Hydrol. 336 (1–2), 1–16.

Li, S.-G., McLaughlin, D., 1991. A nonstationary spectral method for solving stochastic groundwater problems: unconditional analysis. Water Resour. Res. 27 (7), 1589–1605.

Maxwell, R.M., Kastenberg, W.E., Rubin, Y., 1999. A methodology to integrate site characterization information into groundwater-driven health risk assessment. Water Resour. Res. 35 (9), 1855–2841.

Rau, G.C., Andersen, M.S., Acworth, R.I., 2012a. Experimental investigation of the thermal dispersivity term and its significance in the heat transport equation for flow in sediments. Water Resour. Res. 48, W03511.http://dx.doi.org/10.1029/ 2011WR011038.

Rau, G.C., Andersen, M.S., Acworth, R.I., 2012b. Experimental investigation of the thermal time-series method for surface water–groundwater interactions. Water Resour. Res. 48, W03530.http://dx.doi.org/10.1029/2011WR011560.

Rau, G.C., Andersen, M.S., McCallum, A.M., Acworth, R.L., 2010. Analytical methods that use natural heat as a tracer to quantify surface water–groundwater exchange, evaluated using field temperature records. Hydrogeol. J. 18 (5), 1093– 1110.

Reiter, M., 2001. Using precision temperature logs to estimate horizontal and vertical groundwater flow components. Water Resour. Res. 37 (3), 663–674. Roshan, H., Rau, G.C., Andersen, M.S., Acworth, I.R., 2012. Use of heat as tracer to

quantify vertical streambed flow in a two-dimensional flow field. Water Resour. Res. 48, W10508.http://dx.doi.org/10.1029/2012WR 011918.

Rubin, Y., 2003. Applied Stochastic Hydrogeology. Oxford University Press, New York.

Schornberg, C., Schmidt, C., Kalbus, E., Fleckenstein, J.H., 2010. Simulating the effects of geologic heterogeneity and transient boundary conditions on streambed temperatures – implications for temperature-based water flux calculations. Adv. Water Resour. 33 (11), 1309–1319.

Silliman, S.E., Ramirez, J., McCabe, R.L., 1995. Quantifying downflow through creek sediments using temperature time series: one-dimensional solution incorporating measured surface temperature. J. Hydrol. 167 (1–4), 99–119. Shanafield, M., Pohll, G., Susfalk, R., 2010. Use of heat-based vertical fluxes to

approximate total flux in simple channels. Water Resour. Res. 46, W03508.

http://dx.doi.org/10.1029/2009WR007956.

Stallman, R.W., 1965. Steady 1-dimensional fluid flow in a semi-infinite porous medium with sinusoidal surface temperature. J. Geophys. Res. 70 (12), 2821– 2827.

Suzuki, S., 1960. Percolation measurements based on heat flow through soil with special reference to paddy fields. J. Geophys. Res. 65 (9), 2883–2885.

Tabbagh, A., Bendjoudi, H., Benderitter, Y., 1999. Determination of recharge in unsaturated soils using temperature monitoring. Water Resour. Res. 35 (8), 2439–2446.

Taniguchi, M., 1993. Evaluation of vertical groundwater fluxes and thermal properties of aquifers based on transient temperature-depth profiles. Water Resour. Res. 29 (7), 2021–2026.

Taniguchi, M., Sharma, M.L., 1993. Determination of groundwater recharge using the change in soil temperature. J. Hydrol. 148 (1–4), 219–229.

Wörman, A., Riml, J., Schmadel, N., Neilson, B.T., Bottacin-Busolin, A., Heavilin, J.E., 2012. Spectral scaling of heat fluxes in streambed sediments. Geophys. Res. Lett. 39 (23), L23402.http://dx.doi.org/10.1029/2012GL053922.

Zhang, D., 2002. Stochastic Methods for Flow in Porous Media: Coping with Uncertainties. Academic Press, San Diego, Calif.

Zhang, D., Winter, C.L., 1999. Moment equation approach to single phase fluid flow in heterogeneous reservoirs. Soc. Petrol. Eng. J. 4 (2), 118–127.