Rémi Zamansky

Ivana Vinkovic

Mikhael Gorokhovski
email: mikhael.gorokhovski@ec-lyon.fr

LES approach coupled with stochastic forcing of subgrid acceleration in a high Reynolds number channel flow

Keywords: Turbulent channel flow, intermittency, LES, SGS model

In this paper, the high Reynolds number channel flow is simulated by numerical approach at coarse resolution, in which the instantaneous acceleration is decomposed into filtered and subgrid parts, and then both components are modeled. The filtered acceleration is modeled in the framework of the LES approach. The model for the subgrid acceleration is based on two stochastic processes. The first is for its norm, and is based on statistical universalities in fragmentation under scaling symmetry, providing correlation of subgrid forcing across the channel. The second is for its orientation, and is based on the Brownian motion on a unit sphere in order to represent a stochastic relaxation towards full isotropy away from the wall. Two main parameters of the stochastic process include the Reynolds number based on the friction velocity, and the channel half-width. In order to assess the capability of the model proposed, the paper illustrates its application versus recent high Reynolds number DNS, including DNS performed in this paper.

Introduction

Application of large-eddy simulation (LES) in turbulent channel flow becomes too expensive when the Reynolds number is high, and the resolution of strong gradients in the near-wall region is needed. Starting from the earlier study of Moin and Kim [START_REF] Moin | Numerical investigation of turbulent channel flow[END_REF], different modifications of this approach have been proven successful (Meneveau and Katz [START_REF] Meneveau | Scale-invariance and turbulence models for large-eddy simulation[END_REF], Piomelli and Balaras [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF], Sagaut [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows: An introduction[END_REF], Pantano et al. [START_REF] Pantano | Les approach for high reynolds number wall-bounded flows with application to turbulent channel flow[END_REF], Spalart et al. [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF]). Among them one is to incorporate direct simulation of subgridscale turbulence into LES techniques on the basis of estimation of the unresolved velocity field. Approaches proposed by Domaradzki and Adams [START_REF] Domaradzki | Direct modelling of subgrid scales of turbulence in large eddy simulations[END_REF] Park and Mahesh [START_REF] Park | A velocity-estimation subgrid model constrained by subgrid scale dissipation[END_REF], and Westbury et al. [START_REF] Westbury | Analysis of a stochastic backscatter model for the large-eddy simulation of wall-bonded flow[END_REF] are focused on simulation of the residual-stress tensor. The model proposed by Schmidt et al. [START_REF] Schmidt | Near-wall LES closure based on one-dimensional turbulence modelling[END_REF] is based on the LES/ODT (One Dimensional Turbulence) coupling [START_REF] Kerstein | One-dimensional turbulence: model formulation and application to homogeneous turbulence, shear flows, and buoyant stratified flows[END_REF][START_REF] Kerstein | One-dimensional turbulence: A new approach to high-fidelity subgrid closure of turbulent flow simulations[END_REF] and provides stochastically the unresolved velocity field as the wall is approached. The approach of Kemenov and Menon [START_REF] Kemenov | Explicit small-scale velocity simulation for high-Re turbulent flows. Part II: Non-homogeneous flows[END_REF] is based on the simulation of a turbulent velocity field which is defined by superposition of large-scale and small-scale components. The latter is computed directly on subgrid-scales by a simplified 1D motion equation (in three directions). It is worth noting that the above reported subgrid-scale (SGS) models are invariant to the Reynolds number. Thereby, in these models, the intermittency effects are discarded on subgrid-scales.

In the approach proposed by Sabel'nikov et al. in [START_REF] Sabel'nikov | The coupled LES -sub-grid stochastic acceleration model (LES-SSAM) of a high Reynolds number flows[END_REF], the key element of the SGS model is acceleration in the residual flow. This acceleration has been modeled stochastically in the case of three-dimensional box turbulence, with statistical properties depending on the Reynolds number. The general idea of this approach, referred to as LES-SSAM (stochastic subgrid acceleration model), consists of the following. The total instantaneous acceleration

du i dt = - 1 ρ ∂p ∂x i + ν ∂ 2 u i ∂x k ∂x k ( 1
)
∂u k ∂x k = 0
is composed of two contributions:

du i dt = du i dt + du i dt (2)
where dui dt = ∂ui ∂t + u k ∂ui ∂xk The first component is the filtered total acceleration:

du i dt = du i dt + ∂ (u i u k -u i u k ) ∂x k = - 1 ρ ∂p ∂x i + ν ∂ 2 u i ∂x k ∂x k ( 3
)
∂u k ∂x k = 0
where the overbar denotes the filtering operation, u i denotes the filtered velocity.

The second component represents the total acceleration in the residual field:

du i dt = ∂u i ∂t + ∂ (u k u i -u k u i ) ∂x k = - 1 ρ ∂p ∂x i + ν ∂ 2 u i ∂x k ∂x k ( 4
)
∂u k ∂x k = 0
If both parts are modeled, their sum gives an approximation to [START_REF] Moin | Numerical investigation of turbulent channel flow[END_REF].

The first assumption is to replace (3) by a classical LES-equation with a simple eddy-viscosity model

du i dt ∆ = du i dt - ∂ ∂x k ν t ∂u i ∂x k + ∂u i ∂x k = - 1 ρ ∂p ∂x i + ∂ ∂x k ν ∂u i ∂x k + ∂u i ∂x k (5) ∂u k ∂x k = 0
where ν t denotes the eddy viscosity and subscript ∆ is the length scale associated with the filter size.

The second assumption is to replace (4) by

du i dt S = - 1 ρ ∂p * ∂x i + â i ( 6
)
∂u k ∂x k = 0
where p * is the pressure providing the subgrid velocity model-field to be solenoidal, and âi is attributed to the stochastic acceleration of the residual flow.

The third assumption is to replace the sum of ( 5) and ( 6) by the model-equation for a surrogate field, denoted hereafter by the symbolˆ:

dû i dt = - 1 ρ ∂ p ∂x i + ∂ ∂x k (ν + ν t ) ∂ ûi ∂x k + ∂ ûk ∂x i + â i (7)
In this equation, the new pressure field p is such that

∂ ûk ∂x k = 0 (8)
Equations ( 7) and ( 8) will be applied here to a high Reynolds number channel flow, with a new subgrid acceleration model for â i .

In high Reynolds number channel flows, the acceleration is a highly intermittent variable. It has been reported by direct numerical simulations (DNS) [START_REF] Lee | Intermittent nature of acceleration in near-wall turbulence[END_REF], that this intermittency is mostly related to quasi-coherent vortical structures in the wall-layer. A fluid particle, trapped in such a vortex, is subjected to an intense centripetal acceleration. It was revealed [START_REF] Lee | Intermittent nature of acceleration in near-wall turbulence[END_REF] that those fluid particles preserve the magnitude of their acceleration longer than its direction. This finding is consistent with the experimental studies [START_REF] Mordant | Experimental and numerical study of the Lagrangian dynamics of high Reynolds turbulence[END_REF][START_REF] Voth | Lagrangian acceleration measurements at large reynolds numbers[END_REF][START_REF] Voth | Measurements of particle accelerations in fully developed turbulence[END_REF] of a high Reynolds number "free" turbulence. Here, in wall-bounded turbulent flow, the intermittency effects are manifested by deep cross-channel interaction between energy containing eddies of very different length-scales [START_REF] Jiménez | Recent developments on wall-bounded turbulence[END_REF]. Recent DNS studies identified mechanisms, which may be responsible for such turbulent correlations across the channel. In the near-wall region, there are upward ejections of vortical fluid, as a part of vortices-streaks regeneration cycle [START_REF] Kim | On the structure of wall-bonded flows[END_REF][START_REF] Hamilton | Regeneration mechanisms of ner-wall turbulence structures[END_REF][START_REF] Jiménez | The autonomous cycle of near-wall turbulence[END_REF][START_REF] Jiménez | The large-scale dynamics of near-wall turbulence[END_REF][START_REF] Smith | The characteristics of low-speed streaks in the near-wall region of a turbulent boundary layer[END_REF][START_REF] Swearingen | The growth and breakdown of streamwise vortices in the presence of a wall[END_REF]. Above the buffer layer, where the Reynolds number dependent flow [START_REF] Marusic | Streamwise turbulence intensity formulation for flat-plate boundary layers[END_REF][START_REF] Hutchins | Evidence of very long meandering features in the logarithmic region of turbulent boundary layers[END_REF] has intermittent streaky structure as well [START_REF] Jiménez | The largest structures in turbulent wall flows[END_REF][START_REF] Del Alamo | Linear energy amplification in turbulent channels[END_REF][START_REF] Jiménez | Turbulent fluctuations above the buffer layer of wall-bounded flows[END_REF], the pressure fluctuations acting on the flow in the buffer zone and in the viscous sublayer, may also generate eruptions of small-scale wall structures towards outer flow [START_REF] Del Alamo | Scaling of the energy spectra of turbulent channels[END_REF][START_REF] Toh | Interaction between a large-scale structure and near-wall structures in channel flow[END_REF][START_REF] Del Alamo | Self-similar vortex clusters in the turbulent logarithmic region[END_REF]. Besides, the unsteady bursting events, when the large outer-flow eddies appear closely to the wall, has been reported in [START_REF] Jiménez | Characterization of near-wall turbulence in terms of equilibrium and "bursting" solutions[END_REF]. Concerning LES of such flow, the significant role of the "under-resolved" momentum transfer across the channel was emphasized in [START_REF] Härtel | Subgrid-scale energy transfert in the near-wall region of turbulent flows[END_REF].

The question raised in this study is as follows. Despite the huge number of very different turbulent structures, produced by different random realizations in the wall-bounded flow, and interacting across the scales, is it possible to match recent DNS statistics by a rather simple model in the framework of LES-SSAM approach. The purpose of this paper is to propose such a model with assessment of its capability by comparison with DNS. The subgrid acceleration model proposed here is based on the following assumptions : (i) there are two independent stochastic processes, one for the norm, and another for the orientation; (ii) both processes are correlated over the channel height; (iii) the stochastic process for the norm of acceleration is given by the statistical universalities of fragmentation under scaling symmetry; (iv) the stochastic process for the orientation relaxes towards isotropy in the channel center-line. Two main parameters of the stochastic process include the Reynolds number based on the friction velocity, and the channel half-width. This paper is organized as follows. The stochastic model for the non-resolved acceleration will first be presented. The models for the norm and the orientation of the non-resolved acceleration will be described. Then, details about the performed numerical simulations will be given. Finally, the comparison between the classical LES, DNS and the LES-SSAM approach will be presented and discussed.

Stochastic model for the non-resolved acceleration

The stochastic behavior of the subgrid acceleration is prescribed to its norm and to its orientation. It is assumed that the orientation of subgrid acceleration is mostly controlled by orientation of vortex sheets and tubes, lying mostly in the streamwisespanwise plane near the wall. These structures become randomly distributed with increasing distance from the wall. As to the subgrid acceleration norm, it is constructed in terms of wall-parameters, following DNS observation [START_REF] Sabel'nikov | The coupled LES -sub-grid stochastic acceleration model (LES-SSAM) of a high Reynolds number flows[END_REF], in which the significance of the solenoidal component of acceleration very near the wall was shown.

In this study, the suggestion is to introduce the stochastic processes for this norm with correlation across the channel. The subgrid forcing â in [START_REF] Domaradzki | Direct modelling of subgrid scales of turbulence in large eddy simulations[END_REF] is therefore given by two independent stochastic processes, one for its norm |a| and another for its orientation e i , i.e. :

â i = |a|e i (9)
This model will be different from the model proposed in [START_REF] Sabel'nikov | The coupled LES -sub-grid stochastic acceleration model (LES-SSAM) of a high Reynolds number flows[END_REF], where authors introduced a model for the acceleration in homogeneous and isotropic turbulent flow given also by two independent stochastic processes, one for the norm of the acceleration, and another one for its orientation. In [START_REF] Sabel'nikov | The coupled LES -sub-grid stochastic acceleration model (LES-SSAM) of a high Reynolds number flows[END_REF], the model for the acceleration norm has been derived from the lognormal model for the energy dissipation [START_REF] Kolmogorov | A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number[END_REF][START_REF] Oboukhov | Some specific features of atmospheric tubulence[END_REF][START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF][START_REF] Pope | The velocity-dissipation probability density function model for turbulent flows[END_REF][START_REF] Lamorgese | A conditionally cubic-Gaussian stochastic Lagrangian model for acceleration in isotropic turbulence[END_REF] (see also recent paper advocating the lognormal model, such as [START_REF] Arneodo | Towards log-normal statistics in high Reynolds number turbulence[END_REF] , and its application in [START_REF] Sabel'nikov | The extended iem mixing model in the framework of the composition pdf approach: applications to diesel spray combustion[END_REF]). Simultaneously its orientation was generated randomly in each time step of order of the Kolmogorov time scale.

Model for the norm of the subgrid acceleration

We represent the norm of the non-resolved acceleration by a typical velocity increment at the wall expressed as ∆ u 2 * ν and multiplied by a frequency f , which is a random variable:

|a| = ∆ u 2 * ν f (10)
Here u * is the friction velocity. High values of f correspond to strong velocity gradients. Usually the stochastic process evolves with time. In order to introduce the cross-channel correlation in subgrid forcing by a stochastic process for f , we shall replace time by distance from the wall. Thus the random frequency f will evolve stochastically from the wall to the outer flow. To this end, we introduce the non-dimensional evolution parameter τ , which is defined as follows:

τ = -ln h -y h ( 11
)
where h is the channel half-width, and y is the wall distance. It can be seen from ( 11) that when y = 0 then τ = 0, and when y → h then τ → ∞. Hence f evolves with τ , increasing from zero to infinity, correspondingly to y increasing from zero to the channel half-width. The near-wall region is characterized by strong velocity gradients (high values of f ), which decrease in mean towards the outer flow through a highly intermittent dynamics. We assumed that with increasing distance from the wall (increasing τ ) the frequency f is changing by a random independent multiplier α (0 < α < 1), f → αf . α is governed by distribution q(α), 1 0 q(α)dα = 1. This multiplicative process, often referred to as fragmentation under scaling symmetry, requires knowledge of q(α), which is, in principle, unknown. However if relaxation of f is fast, comparing to advancement from the wall to the outer region, the population balance equation, describing continuous evolution of distribution G(f ; τ ) can be reduced exactly to the Fokker-Planck equation [START_REF] Gorokhovski | Statistical universalities in fragmentation under scaling symmetry with a constant frequency of fragmentation[END_REF]:

∂G(f ; τ ) ∂τ = -ln α ∂ ∂f (f G) + ln 2 α 2 ∂ ∂f f ∂ ∂f (f G) (12)
Here G(f ; τ ) is the normalized distribution function,

∞ 0 G(f )df = 1.
It should be noted that only the first two logarithmic moments of α, denoting here as ln k α =

1 0 q(α) ln k αdα, k = 1, 2, appear in the evoluton of G(f ; τ ) (equation 12) and ln α / ln 2 α = ln f / (ln f -ln f ) 2 , with ln f = ∞ 0 G(f ) ln f df .
The solution of ( 12) and its properties are described in [START_REF] Gorokhovski | Statistical universalities in fragmentation under scaling symmetry with a constant frequency of fragmentation[END_REF]. It is shown that with increasing time (here with increasing of parameter τ , O < τ < ∞ instead of O < t < ∞), the initial distribution G(f ; τ = 0) goes first to the log-normal distribution:

G(f ; τ ) = 1 F 2π ln 2 α τ exp   - ln f F 2 2 ln 2 α τ    f F 1-ln α ln 2 α ( 13
)
where F is the initial scale of frequency, which is introduced bellow. It is worth noting that originally, the lognormal distribution has been identified by Kolmogorov [START_REF] Kolmogorov | On the log-normal distribution of particles sizes during break-up process[END_REF] in his discrete model of cascade breakage events. However the Kolmogorov's lognormality was derived in the context of the central limit theorem. The analytical solution [START_REF] Kemenov | Explicit small-scale velocity simulation for high-Re turbulent flows. Part II: Non-homogeneous flows[END_REF] of the population balance equation, obtained without appealing to the central limit theorem, shows that by further increasing of τ , the "log-normal" multiplier in eq. 13 expln f F 2 /2 ln 2 α τ tends to unity, and the long-time limit distribution is determined by the power law with one universal parameter, ln α / ln 2 α . In our paper, instead of using those universal asymptotics, we use [START_REF] Kerstein | One-dimensional turbulence: A new approach to high-fidelity subgrid closure of turbulent flow simulations[END_REF] in order to derive the stochastic equation. In the Itô interpretation, it reads :

df = lnα + ln 2 α /2 f dτ + ln 2 α /2f dW (τ ) ( 1 4 )
where dW (τ ) is a Wiener process ( dW (τ ) = 0 and dW (τ ) 2 = 2dτ , here and below the brackets denote the ensemble average). Parameters are chosen in the following form:

lnα = -Re 1/3 + ln 2 α = -lnα (15)
where

Re + = u * h/ν.
The main properties of the stochastic process ( 14) are: (i) it providing a subgrid forcing that is correlated across the channel height; (ii) as the wall distance increases, it generates evolution of frequency from high to low mean values through two intermediate universal asymptotics: the earlier is log-normal distribution, while the latter is power law, with high probability of zero-frequency and stretched tail of high-frequency. This stretched tail is determined by Re + , and induces rare events of strong subgrid forcing in the outer-flow, thereby mimicking bursting effects.

The starting condition (first grid cell, τ = 0) is given as follows. First, the mean value of the frequency at the wall is introduced as f + = u * /λ, where λ is a Taylorlike scale. It is estimated by Kolmogorov's scaling in terms of wall parameters. The Reynolds number, based on friction velocity, is:

Re + = u * h/ν = h/y 0 ≈ Re 3/4
h , where y 0 is the thickness of the viscous layer, and Re h is the Reynolds number based on the center-line velocity. One then yields: λ ≈ hRe

-1/2 h ≈ hRe -2/3 +
. Second, the starting value for the random path ( 14) is sampled from the stationary log-normal distribution of f /f + .

The probability density function (PDF) of the frequency f has been computed for different wall distances by means of DNS in a channel flow with Re + = 590. Details about this simulation are given in Table 1. First, the norm of the acceleration was computed and then f was estimated by using Equation [START_REF] Schmidt | Near-wall LES closure based on one-dimensional turbulence modelling[END_REF]. As it can be seen from Fig. 1, surprisingly, the model (Equations 14 and 15) matches well the results of DNS. In agreement with DNS data, the stochastic process provides a broad frequency distribution in the wall region and relaxes toward low frequencies when distance from the wall increases.

Model for the orientation of the acceleration vector

We assumed that with increasing distance from the wall, the orientation of subgrid acceleration relaxes stochastically from the streamwise-spanwise alignment (with random azimuth) to full isotropy.

First, this assumption was assessed by DNS performed in this paper. The orientation vector e i in Equation ( 9) is determined by longitude φ and latitude θ:

e i = â i |a| =    e x = cos(θ) cos(φ) e y = sin(θ) e z = cos(θ) sin(φ) (16)
where -π ≤ φ ≤ π characterizes orientation in the streamwise-spanwise (x, z) plane, and -π/2 ≤ θ ≤ π/2 defines orientation relatively to the normal-to-wall direction (θ = 0 and θ = ±π/2 will correspond to accelerations, which are parallel to the wall and normal to the wall, respectively). The schematic representation is given in Fig. 2. If φ and θ are random, their PDF's corresponding to full isotropy, have the following forms P isotropic (φ) = 1/2π and P isotropic (θ) = | cos(θ)|/2, respectively. The mean value for both distributions is zero, and computation of variance yields:

< φ 2 > isotropic = π -π φ 2 P isotropic (φ)dφ = π 2 3 ( 17
)
< θ 2 > isotropic = π/2 -π/2 θ 2 P isotropic (θ)dθ = π 2 4 -2 ( 1 8 )
By DNS at Re + = 180, 590 and 1000, one can obtain the statistics of φ and θ (see Section 3, and especially, Table 1 for details about the resolution used in the simulation). Then, < φ 2 > and < θ 2 > can be computed and compared with [START_REF] Voth | Lagrangian acceleration measurements at large reynolds numbers[END_REF] and [START_REF] Voth | Measurements of particle accelerations in fully developed turbulence[END_REF], respectively. To this end, large-scale acceleration obtained by our DNS has been removed by using sharp spectral filtering in homogeneous directions and Gaussian filtering in the normal to the wall direction. In such filtering operation, the mean of these two angles is negligibly small. The results of computations are given in Fig. 3 for different distances from the wall. It can be seen that the variance of φ (Fig. 3a) is close to the variance of the isotropic distribution (Equation 17) for almost the entire range of y + , and for all three Reynolds numbers. As to variance of θ (Fig. 3b), it is negligibly small at the wall, and then it grows almost linearly with the distance from the wall up to y + ≈ 30 ∼ 50, after which it tends towards a constant value. This value is relatively close to the variance of the isotropic distribution of θ (Equation 18). From this observation, we conclude that the subgrid acceleration vector relaxes towards isotropy, reaching at a high Reynolds number a nearly isotropic distribution at about y + ≈ 50.

In order to reflect this tendency towards isotropy, the stochastic model for the orientation ( 16) is constructed as a Brownian random walk on the surface of a sphere of unity radius. This random walk takes place with increasing distance from the wall. The evolution of the unit vector e i is defined by the following stochastic process:

γ = 2DdW 0 ≤ β < 2π ( 19
)
where γ is the path length between two successive positions on the sphere, and β is the initial direction from point 1 to point 2 (see Fig. 2). β is chosen randomly from the uniform distribution. D is a diffusion coefficient, and dW is the standard Wiener process defined by: < dW >= 0, < dW 2 >= 2dy + , with dy + representing the cell size (in wall units) in the wall normal direction. The rules of this random walk are given by geodesic calculus:

       θ k+1 = sin -1 (sin θ k cos γ + cos θ k sin γ cos β) φ k+1 = φ k + arg(ξ) (ξ) = sin β sin γ cos θ k (ξ) = cos γ -sin θ k sin θ k+1 (20)
where θ k and φ k are two angles corresponding to the node k, and ξ is a complex number, with real and imaginary part (ξ) and (ξ), respectively. With increasing distance from the wall (y + ), the diffusive equilibrium is attained as a final state, which corresponds to the isotropic PDF's:

P θ (θ, y + → ∞) → P isotropic (θ) P φ (φ, y + → ∞) → P isotropic (φ) (21)
where, at the distance from the wall y + , P θ (θ, y + ) and P φ (φ, y + ) are the PDF's of θ and φ, respectively. An example of one realization of such a Brownian motion on a sphere is given in Fig. 4. Correspondingly to DNS, the process starts on the wall (y + = 0) with:

P θ (θ, y + = 0) = δ(θ) P φ (φ, y + = 0) = P isotropic (φ) ( 22
)
where δ is the Dirac distribution. This implies that the subgrid scale acceleration at y + = 0 is taken to be parallel to the wall with random azimuth. In our case, the diffusion coefficient D controls the relaxation rate towards isotropy. A linear regression of DNS data (Fig. 3b for y + < 20), yields that D = 0.01. Fig. 5 shows statistics of θ computed from the model proposed here (Equations 19, 20 and 22, and D = 0.01) and from our DNS. It can be seen that the model matches well the DNS data. For φ, with the prescribed wall condition [START_REF] Jiménez | The autonomous cycle of near-wall turbulence[END_REF], the model provides exactly < φ 2 >=< φ 2 > isotropic .

Characteristics of the numerical simulations

For assessment of the SGS acceleration model, we ran simulations of a pressure driven turbulent channel flow. A pseudo-spectral computational code is used. The code uses the spectral approximation (Fourier Chebyshev) and a variational projection method on a divergence free space following the original idea of Moser [START_REF] Moser | A spectral numerical method for the Navier-Stokes equations with applications to Taylor-Couette flow[END_REF]. The time integration uses the explicit Adam-Basforth algorithm for non-linear terms, and the semi-implicit algorithm for diffusion terms. A rotational form is used for convective terms in order to ensure energy conservation. Periodic boundary conditions were applied along the streamwise (x) and the spanwise (z) directions, whereas the no-slip condition is imposed on the walls. The computer code was developed in the scientific group of M. Buffat [START_REF] Pascal | Etude d'une turbulence compressée et ou cisaillée entre deux plans parallèles: comparaison entre approche statistique et simulation des équations de Navier-Stokes instantanées[END_REF].

A posteriori tests have been done for three Reynolds numbers: Re + = 590, 1000 et 2000. Results of LES-SSAM approach were compared with the standard LES and the DNS given in the literature [START_REF] Moser | Direct numerical simulation of turbulent channel flow up to Reτ = 590[END_REF][START_REF] Hoyas | Reynolds number effects on the reynolds-stress budgets in turbulent channels[END_REF] and also with the DNS performed specifically in the present work. The conditions for those three approaches were exactly the same, including the same computational mesh for LES and LES-SSAM. The parameters taken for these simulations are summarized in Table 1. [START_REF] Moser | Direct numerical simulation of turbulent channel flow up to Reτ = 590[END_REF] c Hoyas and Jiménez (2008) [START_REF] Hoyas | Reynolds number effects on the reynolds-stress budgets in turbulent channels[END_REF] For LES and LES-SSAM simulations, the classical Smagorinsky model is used with wall damping function for the turbulent viscosity [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows: An introduction[END_REF]:

ν turb = (C s ∆f V D ) 2 |S| |S| = (2S ij S ij ) 1/2 f V D = 1 -e -y/A (23)
where A is specified below. Here C s is the Smagorinsky constant, ∆ = (∆x × ∆y × ∆z) 1/3 is the characteristic grid size,

S ij = 1 2 ∂ ūi ∂xj + ∂ ūj ∂xi
is the resolved strainrate tensor, f V D is the Van Driest function, y is the distance from nearest wall, and A is the constant controlling the damping function f V D . As suggested by [START_REF] Shur | A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities[END_REF], for a detached-eddy simulation (DES), the turbulent length-scale in the wall-layer is defined as = min(y, ∆). In order to fulfill this suggestion, in the Van Driest damping function, we choose A such that ∆ × f V D ∼ min(y, ∆).

It should be noted that in this code, the Reynolds number is given by choosing ν and the pressure gradient-1 ρ ∂Pt ∂x . In order to attain the required Reynolds number we used the Dean's correlation [START_REF] Dean | Reynolds Number dependance on skin Friction and other bulk flow variables in twodimensional rectangular duct flow[END_REF]: ν = 0.110U c hRe -1.1296 + and-

1 ρ ∂Pt ∂x = Re 2 + ν 2 /h 3
, where U c is the center-line velocity. It can be noticed from Table 1 that the LES-SSAM matches better the DNS Reynolds number than the standard LES. Thereby, for a given set of parameters (ν and -1 ρ ∂Pt ∂x ), both the center-line velocity and the mass flow rate are improved by LES-SSAM model.

Results and discussion

For the sake of simplicity, we will only present simulations for the 64 × 65 × 64 grid, for all three Reynolds numbers. For finer resolutions the difference between LES-SSAM and standard LES is less pronounced, but still present. Fig. 6 shows evolution of the mean velocity across the channel. The LES-SSAM reproduces well the logarithmic law of the wall predicted by DNS, which is not the case with standard LES. Fig. 7 shows also the improvement in the prediction of the velocity variance when LES-SSAM is used. This can be seen especially for the streamwise velocity compared to the DNS of [START_REF] Hoyas | Reynolds number effects on the reynolds-stress budgets in turbulent channels[END_REF]. For the other components, the profile is also better predicted by LES-SSAM than by the standard LES. It should be noted that the variance of the streamwise velocity at Re + = 1000 computed by our DNS deviates slightly from the DNS in [START_REF] Hoyas | Reynolds number effects on the reynolds-stress budgets in turbulent channels[END_REF]. This may be due to a smaller computational domain in comparison with [START_REF] Hoyas | Reynolds number effects on the reynolds-stress budgets in turbulent channels[END_REF] and a higher Reynolds number.

Fig. 8 illustrates the vertical profiles of the turbulent and viscous stresses, τ turb = -ρ u v and τ visc = -ρν ∂ u ∂y , respectively. The results are presented as ratios τ turb /(τ turb + τ visc ) and τ visc /(τ turb + τ visc ). Once again, the advantage of the LES-SSAM approach versus the classical LES is explicitly shown. Both the viscous and the turbulent stresses are better predicted by LES-SSAM.

Velocity spectra are illustrated in Fig. 9. For both Reynolds numbers and both y + shown in the figure, it can be clearly seen that the high wave number part of the spectrum is much better resolved by LES-SSAM than by LES. Fig. 10 represents the evolution of the longitudinal autocorrelation coefficient of the velocity components at y + = 5, representing the near-wall region. The improvement of the correlation length due to the better resolution of the small-scale part of the spectrum by LES-SSAM is clearly seen in this figure. Another observation is that for small wall distances, the streamwise longitudinal autocorrelation is more spatially correlated than the remaining two components, because of the alignment of streaky structures with the wall. Fig. 11 illustrates the PDF of the spanwise acceleration. At three very different distances to the wall, the PDF obtained by the DNS displays stretched tails, as a manifestation of intermittency. The same evidence of intermittency is obtained by LES-SSAM, while in the case of the standard LES those tails are much less developed. This effect is also observed for the PDF of the other two components of the acceleration.

Conclusion

In numerical simulation of a high-Reynolds number channel flow at moderate resolution, the instantaneous acceleration has been split into filtered and subgrid parts. Then, both components were modeled, in the framework of LES combined with the stochastic subgrid acceleration model.

The aim of the new model for subgrid acceleration was to account for intermittency in the residual (unfiltered) wall-bounded flow. This model was constructed in accordance with DNS observations, reported in the literature, and DNS performed specifically in the present paper. The model is based on two independent stochastic processes: one for the norm of the subgrid acceleration, the other for its orientation.

The norm is given by a characteristic velocity increment, in terms of wallparameters, multiplied by stochastic frequency. The stochastic process, constructed in lines of statistical universalities in fragmentation under the scaling symmetry, evolves with distance from the wall towards outer flow. To this end the nondimensional evolution parameter was introduced, which increases from zero to infinity, as the distance to the wall increases from the wall to the channel half-width. The multiplicative stochastic process provides correlation of the subgrid forcing across the channel. With increasing wall distance and consequently, decreasing mean frequency, this process has two intermediate universal asymptotics. The earlier is the log-normal distribution, while the later is power distribution. The power law has a high peak of probability for zero-frequency and a stretched tail for high-frequencies, which depends on the Reynolds number. Sampling from that stretched tail induces rare events of strong subgrid forcing in the outer-flow, mimicking thereby bursting events.

It was assumed that orientation of the subgrid acceleration is controlled by orientation of vortices, which become randomly distributed, if distance from the wall increases. In order to fulfill this tendency towards full isotropy, the stochastic model for subgrid acceleration orientation is based on a Brownian random walk on the surface of a sphere of unity radius. Starting close to the wall, from the direction in the streamwise-spanwise plane with random azimuth, the random walk takes place with increasing distance from the wall, and tends to the diffusive equilibrium as a final state, which corresponds to isotropic PDF's for the two geodesic angles, longitude and latitude. The diffusion coefficient was computed by DNS and found to be almost independent of the Reynolds number.

In the case of turbulent channel flow, the results of the LES-SSAM approach were compared with DNS data (Re + = 590, 1000 and 2000) and with standard LES. The LES-SSAM model was proven successful in predictions of DNS results: (i) resolution of the small-scale part of the energy spectrum and simulation of spatial correlations are close to DNS; (ii) stretched tails in the PDF of acceleration are predicted in agreement with DNS; (iii) computation of mean and of variance of the velocity components, as well as of turbulent and viscous stresses are explicitly improved. Along with fundamental interest, the LES-SSAM approach has a practical relevance, when significant physics take place on subgrid scales. In such conditions, simulations at a high Reynolds number warrant advanced SGS models that account for intermittency effects on small spatial scales. Examples include scalar mixing and turbulent combustion, dispersion, vaporization and combustion in twophase flows. This requires the further development of this approach with many associated questions, concerning complex geometry, correlation between norm and orientation of the subgrid acceleration, time-correlations, etc. However, the approach given in [START_REF] Sabel'nikov | The coupled LES -sub-grid stochastic acceleration model (LES-SSAM) of a high Reynolds number flows[END_REF], and in the present paper is very simple, and easily realizable for any LES code without supplementary CPU requirement.

Figure 1 .e 1 = (θ 1 , φ 1 )e 2 = (θ 2 ,

 111122 Figure 1. Distribution of f/f+ from SSAM (cross) and comparison with DNS (line) at Re+ = 590, for several distances from the wall.

Figure 2 .Figure 3 .

 23 Figure 2. Coordinate system, and definition of the angles φ, θ, γ and β.

Figure 4 .Figure 5 .Figure 6 .

 456 Figure 4. Realization of the random walk on the sphere given by Equations 19 and 20.

Table 1 .

 1 Summery of parameters used for numerical simulations. dt + = dt/t * is the time step of the simulation; t * = ν/u 2 * is a viscous time of the order of the Kolmogorov time scale at the wall.

	Name	Re+	Rec	Nx × Ny × Nz	Lx × Ly × Lz	∆x + × ∆y + × ∆z +	dt +a	Cs	A/h
	DNS	180	3280	192 × 193 × 192	3πh × 2h × 4 3 πh	9.0 × (0.02 ∼ 3.0) × 4.0	0.030	-	-
	DNS DNS b LES LES-SSAM	587 587 587 587	12490 12547 14160 12760	384 × 257 × 384 384 × 257 × 384 64 × 65 × 64 64 × 65 × 64	3 2 πh × 2h × 3 4 πh 7.2 × (0.04 ∼ 7.2) × 3.6 2πh × 2h × πh 9.7 × (0.04 ∼ 7.2) × 4.8 3πh × 2h × πh 87 × (0.71 ∼ 29) × 29 3πh × 2h × πh 87 × (0.71 ∼ 29) × 29	0.033 -0.1 0.1	--0.16 0.16	--0.015 0.015
	DNS DNS c LES LES-SSAM 1000 1000 934 1000 LES 1000 LES-SSAM 1000	22250 20960 25430 23380 25500 23700	512 × 385 × 512 3072 × 385 × 2304 96 × 97 × 96 96 × 97 × 96 64 × 65 × 64 64 × 65 × 64	4 3 πh × 2h × 2 3 πh 8.2 × (0.03 ∼ 8.3) × 4.1 8πh × 2h × 3πh 7.6 × (0.06 ∼ 7.6) × 3.8 3πh × 2h × πh 99 × (0.53 ∼ 33) × 33 3πh × 2h × πh 99 × (0.53 ∼ 33) × 33 3πh × 2h × πh 147 × (1.2 ∼ 49) × 49 3πh × 2h × πh 147 × (1.2 ∼ 49) × 49	0.034 -0.1 0.1 0.1 0.1	--0.16 0.16 0.2 0.2	--0.009 0.009 0.015 0.015
	DNS c LES LES-SSAM 2000 2003 2000 LES 2000 LES-SSAM 2000	48680 49350 48950 52640 49050	6144 × 633 × 4608 128 × 129 × 128 128 × 129 × 128 64 × 65 × 64 64 × 65 × 64	8πh × 2h × 3πh 3πh × 2h × πh 3πh × 2h × πh 3πh × 2h × πh 3πh × 2h × πh	8.2 × 8.9 × 4.1 147 × (0.60 ∼ 49) × 49 147 × (0.60 ∼ 49) × 49 295 × (2.4 ∼ 98) × 98 295 × (2.4 ∼ 98) × 98	-0.1 0.1 0.1 0.1	-0.16 0.16 0.2 0.2	-0.006 0.006 0.015 0.015

a b Moser et al. (1999)

Acknowledgments

M. Buffat is acknowledged for the development of the computational code. The authors express their gratitude to F. Laadhari who kindly provided his initial fields for the DNS. It is a pleasure to thank V. Sabel'nikov, P. Moin and J. Jiménez for discussions on this subject with many benefits for us. This work was granted access to the HPC resources of CINES under the allocation 2009-c200902560 made by GENCI (Grand Equipement National de Calcul Intensif)


Franc ¸ois Lavoué
email: francois.lavoue@ens-lyon.org

Bernard Tourancheau

Eléments pour la Modélisation et le Dimensionnement des Sondes Géothermiques de Faible Profondeur : Influence des propriétés thermiques des sols Modeling and Dimensioning Ground Heat Exchangers Principles: Influence of the soil's thermal proprieties

Keywords: Ground Heat Exchangers (GHE), Sensitivity Study, TRNSYS, Ground-Coupled Heat Pump

R ÉSUM É. Dans un contexte de crise énerg étique et de r échauffement climatique, le chauffage des b âtiments gr âce à l' énergie solaire stock ée dans le sol est un principe int éressant. De plus, les sondes g éothermiques utilis ées pour extraire la chaleur du sol permettent aussi de l'y r éinjecter, et donc de refroidir les b âtiments gr âce à l'inertie thermique du sol. Cet article pr ésente nos r éflexions sur la mod élisation et le dimensionnement des échangeurs thermiques en faible profondeur coupl és à des pompes à chaleur (PAC g éothermique). Apr ès une pr ésentation de la physique sous-jacente aux échangeurs thermiques et des m éthodes de r ésolution analytiques et num ériques, nous soulignons l'importance des valeurs des param ètres thermiques du sol utilis ées par les outils num ériques pour simuler les PAC g éothermiques. Nous illustrons cette question à partir de r ésultats de simulation dynamique TRNSYS : ceux-ci sont fortement influenc és par les valeurs de conductivit é thermique du sol. Enfin, nous pr ésentons un état de l'art des mesures de ces param ètres thermiques in situ et sugg érons quelques m éthodes g éophysiques pour les d éterminer indirectement, parmi lesquelles le g éoradar, l'induction électromagn étique et la tomographie électrique. La forte empreinte de la teneur en eau sur la thermique des sols permettrait en effet d'estimer leurs caract éristiques thermiques via des mesures de conductivit é électrique ou de permittivit é di électrique, elles aussi reli ées à la teneur en eau du milieu. MOTS-CL ÉS: Echangeurs thermiques, Etude de sensibilit é param étrique, TRNSYS, Pompe à chaleur (PAC) g éothermique ABSTRACT. In the context of energy crisis and global warming, heating buildings with the solar energy stored in the soil represents a very interesting alternative. Moreover, cooling buildings can also use the soil damping capacity. This paper presents our reflexion about the modeling and dimensioning of the ground heat exchanger part of ground-coupled heat pumps (GCHP). After a physical overview of the ground heat exchanger, we extract from analytical solutions practical consequences of the soil damping behavior and limits that provide guidelines for the dimensioning. We then question the default values of the numerical tools thermal parameters used for the simulation of GCHP. We illustrate this issue through a TRNSYS dynamic simulation of GCHP, demonstrating that soil's thermal parameters have a strong impact on the results. Finally, we give some perspectives for the determination of soil's characteristics in situ, but indirectly, thanks to geophysical prospection methods as ground penetrating radar, electromagnetic induction, or electrical resistivity tomography.

Nomenclature:

Indices: T temperature (K)

INTRODUCTION

In the present energetic context, where fossil energies must face global warming, it appears as a necessity to find new, and renewable, energy sources. In order to heat buildings, an alternative could be the solar energy stored in the soil. The recent decades have seen the rapid development of ground-coupled heat pumps (GCHP) in the USA and Europe [START_REF] Lund | Geothermal (ground-source) heat pumps, a world overview[END_REF]). A GCHP is a conventional heat pump (HP) coupled with a ground heat exchanger (GHE): tubes buried in the soil in which circulates a calorific fluid. Using a stable temperature source at depth, a GCHP is more efficient than a conventional air-to-air HP and is more independent from climatic conditions. This efficiency must therefore be payed by expensive tube burying costs. Hence depth is often not a choice and a design study is necessary. As a consequence, many researchers deal with modeling and simulation of GHE for their dimensioning.

As reported by [START_REF] Lamarche | A new contribution to the finite-line model for geothermal boreholes[END_REF], the modeling of GHE involves two heat transfer problems. The first transfer is a partly convective heat flux between the fluid in the tubes and the wall of the borehole. It is mostly solved analytically and we find examples of this analytical resolution in [START_REF] Hollmuller | Utilisation des échangeurs air/sol pour le chauffage et le rafraîchissement des bâtiments Mesures in situ, modélisation analytique, simulation numérique et analyse systémique[END_REF] and [START_REF] Hellström | Ground Heat Storage Thermal Analysis of Duct Storage Systems[END_REF]. Hollmuller gives a nice demonstration for the temperature evolution inside the tubes, from which he deduces simple dimensioning rules. Hellström lays the basis for the line-source theory (LST) which is used in most GHE simulations.The second transfer is considered by most researchers to be a pure conductive heat flux between the borehole and the soil. This conduction phenomenon is mostly numerically treated. Here we study the influence of tubes' length for heating or cooling a particular building using a TRNSYS type developed by [START_REF] Giardina | Evaluation of ground coupled heat pumps for the state of wisconsin[END_REF].

Finally, we aim to focus on the role of soil's thermal proprieties (section 4). Indeed, these parameters have a great effect on the efficiency of GHE: we intend to show this with TRNSYS simulations. Here we discuss several methods for measuring these proprieties in situ or in laboratory. Then we investigate a relation between electrical and thermal resistivity empirically, shown by [START_REF] Singh | A generalised relationship between soil electrical and thermal resistivities[END_REF], which seems promising.

PHYSICAL PRINCIPLES OF HEAT EXTRACTION THROUGH GROUND-COUPLED HEAT PUMP

THE HEAT PUMP: A THERMODYNAMICAL MACHINE

A heat pump is a thermodynamical machine. Provided an electrical work W , it extracts a heat amount Q 1 from a cold source (whose temperature we note T 1 ) and transfers a heat amount Q 2 to a hot source, whose temperature is T 2 (see figure i). In the case of GCHP, the cold source is the fluid at the tubes output (i.e. T 1 = T f o ) and the hot source is the fluid to be warmed: distributed air, warming fluid, or domestic hot water. The coefficient of performance (COP) of such a machine is defined by:

COP = Q useful W = Q 2 W = - Q 2 Q 1 + Q 2 = T 2 T 2 -T f o (1)
If we consider the machine as ideal (Carnot's machine), we can express its COP with the 1 st and the 2 nd thermodynamical principles like in equation ( 1). This ideal COP is seen as the theoretical maximal COP of real machines. As a consequence, we must get a temperature T f o as large as possible at the tubes output to improve the COP. It means that the heat exchange between soil and tubes must be optimized. We can notice that cooling is free and doesn't require any electrical work in summer.

HEAT CONDUCTION AND THERMAL INERTIA IN THE GROUND

In a homogeneous soil, heat transfers obey to

∂ t T = α ∇ 2 T (2)
where α is the thermal diffusivity defined by α = λ/(ρC P ).

The heat law allows to calculate some orders of magnitude. For example, we can determine the time τ that is necessary to reach the permanent regime in a ground where is buried a GHE. We assume for this calculation that the soil's temperature is constant far from the GHE (i.e. at a characteristic distance L from the GHE). The heat law can thus be rewritten:

∆T τ = α ∆T L 2 (3)
and we deduce τ = L 2 /α. For example, the heat exchange between a vertical GHE and the surrounding soil with α = 10 -6 m 2 .s -1 reaches permanent regime at a distance of 5 m in 9 months and 19 days. After this time, the heat exchange is null and the system is inoperative. We see here the risk of using intensively such a heat exchanger: we must not pump heat continuously but intermittently. Inversely, the same type of calculation can determine the distance we must put between tubes when their operating duration is known. Finally, soil's thermal parameters have an effect on this calculation through the diffusivity α.

MODELING OF GROUND HEAT EXCHANGERS

To model GHE more precisely than with rough orders of magnitude, one must solve two distinct problems of heat transfer. [START_REF] Lamarche | A new contribution to the finite-line model for geothermal boreholes[END_REF] summarized these problems as follows:

(1) Heat transfer between the fluid in the tubes and the wall of the borehole. Due to the flow in the tubes, this transfer is partly convective. Most researchers solve this problem analytically making some approximations that we detail in subsection 2.1.1. Numerical approaches are used to deal with long-term evolution of the system (see section 2.2.1). (2) Heat transfer between the borehole and the soil by conduction. This problem needs a numerical resolution because the conduction is a transient phenomenon which takes place in a soil that is an heterogeneous and three-dimensional medium.

ANALYTICAL SOLUTIONS

On the basis of analytical solutions, we can understand the physics of heat transfers in GHE and it is possible to parameterize them. However, in order to get these solutions, we need to approximate the system and we cannot consider heterogeneous proprieties of the soil or complex dispositions of the tubes.

Temperature profile in undisturbed soil

Before we extract any heat, we must know the soil's temperature profile, as a function of depth and time. In order to determine this profile, we make the following assumptions: H1: the surface temperature can be expressed as

T S (z = 0,t) = T av + ∆T. cos [ω(t -t min )] ( 4
)
where ω is the frequency of interest (ω = 2 × 10 -7 rad.s -1 for a year and ω = 7.27 × 10 -5 rad.s -1 for a day), T av et ∆T are the average temperature and the amplitude of temperatures in a period τ = 2π/ω, and t min is the date (in s) of the minimal temperature in the period, H2: the soil is homogeneous. In particular, it has a constant thermal diffusivity α s (with depth, time and temperature), H3: there is no geothermal gradient in the ground (but only solar energy input), H4: there are no water flows, so that the heat exchanges are purely conductive. We can note that this is a very strong assumption. In real soils, water flows exist, vary depending on seasons and weather, induce convective heat exchanges and dramatically disturb the temperature profile we will calculate here.

To determine the temperature profile, we solve the heat equation (2). [START_REF] Kusuda | Earth temperature and thermal diffusivity at selected stations in the united states[END_REF] were first to solve this equation in the soil. Many authors applied their work, especially in TRNSYS type 556 developed by [START_REF] Giardina | Evaluation of ground coupled heat pumps for the state of wisconsin[END_REF]. The solution, drawn on figure ii can be written as follows:

T (z,t) = T av + ∆T. e -z/δ . cos ω(t -t min ) - z δ (5)
This expression shows that the ground's thermal inertia induces both a phase delay and a damping of the temperature wave in depth. Both effects can be characterized by the penetration coefficient

δ = 2 α s ω (in m)
This penetration coefficient can be interpreted as the depth where the surface signal is damped by a 1/e factor. It determines the particular depth z o where the temperature oscillation is in phase opposition with the surface temperature. Thus, in theory, horizontal tubes could extract heat stored in summer for use in winter if they are buried at a depth of z o = π δ.

Unfortunately, at this depth, the temperature damping is already e -π = 0.043, which means a 95.7% loss of the signal intensity. Although extracting heat stored in the summer during winter seems attractive, it is not beneficial. Nevertheless, the damping is still an interesting phenomenon: the temperature is increasingly stable with growing depth. The expression of undisturbed temperature (5) is thus a first relation for dimensioning a GHE.

2

oling effect produced at the evaporator as well as the heating effect produced enser. In these dual-mode GSHP systems, a reversing valve is used to switch eating and cooling modes by reversing the refrigerant flow direction. The third system is the ground loop in which water or an antifreeze solution exchanges he refrigerant and the earth. iciencies of GSHP systems are much greater than conventional air-source heat ems. A higher COP (coefficient of performance) can be achieved by a GSHP 2.1.2 Modeling of the temperature at the tubes output Subsection 1.1 demonstrated that greater the temperature T f o at the tubes output is, better is the COP of the heat pump. In his Ph.D. thesis, [START_REF] Hollmuller | Utilisation des échangeurs air/sol pour le chauffage et le rafraîchissement des bâtiments Mesures in situ, modélisation analytique, simulation numérique et analyse systémique[END_REF] determines an analytical solution for the temperature at the x-coordinate in the tube:

T f (x,t) = T f i + (T f o -T f i ). exp - 2πr o C P f D f h x . cos ω t - x v f - 2πr o C P f D f k x (6)
where h and k are functions of soil diffusive coefficients and of the convective coefficient h f of the fluid.

It is important to notice in expression ( 6) that the convection of the fluid shifts the temperature phase between the fluid and the ground. Therefore, the lower the fluid velocity is, the greater is this phase difference. Consequently, this phenomenon could solve the problem whereby temperature phase shift through the earth induces a damping (subsection 2.1.1): reducing the fluid flow in the tubes also shifts the output temperature phase (this technique is called low flow).

Moreover, according to [START_REF] Hollmuller | Utilisation des échangeurs air/sol pour le chauffage et le rafraîchissement des bâtiments Mesures in situ, modélisation analytique, simulation numérique et analyse systémique[END_REF], some values of h f can privilege phase difference against damping. A phase delay almost without damping has even been experimentally brought out at small scale and over a time period of a day with a laminar fluid flow between two concrete plates. Unfortunately, this configuration cannot be applied to heat an entire building. Nevertheless, it could offer possibilities for horizontal GHE systems.

Models of temperature in the ground around the tubes

The Line-Source Theory [START_REF] Hellström | Ground Heat Storage Thermal Analysis of Duct Storage Systems[END_REF]) proposes the following expression for the temperature in the soil around the borehole (located at the coordinate r = 0):

T (r, z,t) -T soil (z,t) = q 4πλ s +∞ r 2 4 αs t e -u u du (7)
where T soil is the undisturbed temperature given in subsection 2.1.1 and q is the heat flow per length unit along the tube (W/m). This model describes a global heat transfer along the entire tube that we consider to be a single section. In the earth, since the temperature difference between the soil and the fluid decreases along the tube (in the case of heat extraction), the heat flux is probably greater in the first sections of the tube than in the last ones. As a consequence, the heat flux which is calculated analytically with the LST often seems lower than this, which is numerically determined (see type 556, [START_REF] Giardina | Evaluation of ground coupled heat pumps for the state of wisconsin[END_REF]).

We can obtain the temperature profile in a ground where several boreholes are buried with the superposition principle of this type of solutions. However, we have to apply this superposition principle continuously i.e. at each time step of the calculation to attend the interference phenomena between each borehole. Thus, the modeling must be numerically processed.

NUMERICAL MODELING

Long-term interferences and thermal deviation in the soil

Lee and Lam (2008) used an implicit finite difference method to study the problem of thermal deviation in the soil that results from an overpumping of heat by several boreholes. To do this, the authors described the tubes with cylindrical slices. In each slice, they considered the mean temperatures T f of the fluid and T b of the borehole wall. They assumed that the heat transfer q between the fluid and the borehole obeys the law:

T f -T b = R b q (8)
where R b is the thermal resistance of the entire GHE (tubes, fill material and wall of the borehole). This approach draws an analogy between electric and thermal flux to describe the heat transfer, so that the resistance R b is considered to be one of the principal parameters of a GHE system. [START_REF] Hellström | Ground Heat Storage Thermal Analysis of Duct Storage Systems[END_REF] proposes an analytical expression for the resistance R b as a function of geometrical and thermal characteristics of the tubes. [START_REF] Lee | Computer simulation of borehole ground heat exchangers for geothermal heat pump systems[END_REF] numerically get the temperature profile by superposing the solutions that they obtain for several boreholes. Their results suggest that a finite difference scheme is insufficient to estimate the performance of an entire field of boreholes. They conclude that we need to discretize simultaneously each borehole of the field to improve the accuracy of dynamical studies.

Numerical tools for integrating GHE into building simulation

In energy simulation of buildings most of the physical processes are modeled with a set of algebraic equations, mostly non-linear, differential and integral. Several tools exist which solve these equations with different methods. The dynamic simulation systems are the more precise because they solve equations iteratively over time with a given time-step. Such a well-known system is TRNSYS from University of Wisconsin1 or Comfie/Pleiade3 from Ecole des Mines. For our simulations, we use the TRNSYS tool.

INFLUENCE OF THE SOIL'S THERMAL PARAMETERS

LITERATURE'S VALUES OF THERMAL PARAMETERS AND THEIR RELIABILITY

Some values of thermal conductivity λ, densities ρ and heat capacity C P are given in books by Turcotte and Schubert (2001) and by [START_REF] Sleep | Principles of Geophysics[END_REF], others are published on the Wikipedia website. We summarize them in the table i. About the reported data, we can make some remarks:

(1) According to these values, a diffusivity α = 10 -6 m 2 .s -1 , as we often find as default value in the simulations tools, would correspond to a granitic rock or to a soil composed by gravel, silt or watersaturated sand. But this value is not consistent with a dry soil or with clay. (2) Water increases the values of ρC P and λ. The increase of ρC P is due to the great heat storage capacity of water. The increase of λ whereas water is rather a better thermal insulator than the soil's mineral components, may be due to thermal bridges that water build between the soil's particles. At the contrary, in a dry soil, the interstices are filled with air, a very good insulator. Thus, water is a very important component of soils whose content we must precisely know. Moreover, as moisture content varies seasonally and daily, we can imagine that water-dependent thermal parameters vary with time as well.

(3) According to [START_REF] Fromentin | Recommandations pour la réalisation d'installations avec pieux échangeurs[END_REF], water has also a convective effect when in motion in soils. In a permeable layer, water can disperse the heat that is transferred from the tubes in the soil. On the one hand, this convective effect prevents a seasonal storage of heat in summer but on the other hand, it enables a quicker thermal reload in winter. The authors emphasized on the necessity to better understand heat transfers in porous media to model this effect. In this paper, we will neglect the convective phenomena, assuming that water does not move in soil.

Some thermal properties of complex soils were given by [START_REF] Giardina | Evaluation of ground coupled heat pumps for the state of wisconsin[END_REF] and used in TRNSYS type 556. They are summarized in table ii. Apart from the vague definition of the considered soils, without any information about composition, geologists would be astonished by the given densities: an heavy saturated soil is here as heavy as mantle rocks! Though, we saw in section 1.2 that soil's thermal parameters can largely impact calculations and dimensioning of the GHE. In the next section, we present TRNSYS simulation results to quantify the mistake due to erroneous thermal parameters values.

DIMENSIONING A GEOTHERMAL SYSTEM USING NUMERICAL SIMULATION EXPERIMENT

The previous remarks are supposed to have effects in dimensioning studies. In order to verify it, we simulate with the TRNSYS software the working of a commercialized product that contains:

• High-density polyethylene tubes, which must be horizontally buried.

• An hydraulic module with a 3-stages pump and input-output thermometers.

• An air/water heat exchanger.

• An 2-stages outdoor thermostat.

• Specific connections.

All components of the system already have fixed characteristics (thermal proprieties, pump flow, exchanger efficiency), so that the only parameters that are still to be determined is the length of the tubes and the soil's thermal parameters. As the constructor gives a value of around 1.2 m for the depth at which the tubes must be buried, we use this value to reduce the number of freedom degrees (although it may be varied to optimize the heat exchange).

Numerical simulation parameters are tested through the ground exchanger model implemented in TRNSYS type 556 in the climate of Saint Etienne, France, over a typical year. We tried to evaluate the maximum quantity of energy extractable for heating or cooling from the ground. This was done for several soil types corresponding to the parameters described in Table ii. Our results presented in Table iii show that in that range of usage, the main influent parameter is the soil's thermal conductivity while its density has a second role.

Thus, soil's parameters characterization and values choice for the numerical simulation is crucial. For instance, the water present in the soil will drastically change its conductivity and hence the thermal exchanges. One can note that the numerical tool's default values are somehow surprising. Indeed, while at low depth most soils are composed of sand and topsoil which correspond to Giardina's light soil, one can often found values corresponding to heavy saturated soils.

Figure iv shows the heat exchange as a function of the conductivity. For heat extraction, only low conductivity soils benefit from longer pipes. In the over hand, for cold extraction the pipe length is beneficial for all type of soils.

Figure iii obtained from simulations shows that for a given load to extract, soil's conductivity has significant influence on feasibility. Eventually, pipe length will allow for adjustment, especially with low conductivity. Thus the default values of the simulation tools can likely lead to system malfunction and mis-dimensioning for both heat and cold extraction.

For example, assuming that we want to extract 10 5 kJ/y from a light dry soil (λ = 0.35 W/m/K), figure iii shows that a calculation leads with a soil's conductivity λ = 1.17 W/m/K (so that α = 10 -6 .s -1 as often assumed) will suggest a pipe's length of 150 m (red dashed line). But in reality, 50 m would have been sufficient (red solid line): wecould have made economies. Inversely, simulating the GHE with an over-estimated soil's conductivity, we will not deduce enough length of pipe. Thus, the same heat exchange will induce a stronger temperature gradient, running the risk of a soil's thermal impoverishment on the long view.

EXPERIMENTAL DETERMINATION OF THE SOIL'S THERMAL PROPRIETIES

STATE AND OF THE ART AND ITS LIMITS

Most of studies on GHE modeling underline the importance of the soil's thermal proprieties to design efficient systems (e.g. [START_REF] Hollmuller | Utilisation des échangeurs air/sol pour le chauffage et le rafraîchissement des bâtiments Mesures in situ, modélisation analytique, simulation numérique et analyse systémique[END_REF], [START_REF] Chiasson | Advances in modeling of ground-source heat pump system[END_REF]). However, they rarely focus on this aspect. In the literature, there are few convenient methods to measure thermal parameters in soils at depth precisely.

Usually, their values are estimated from a rough composition evaluation and tabulated values as in table i. Three principal elements affect these thermal proprieties:

• The composition of the soils and rocks in which the GHE are buried. As we saw in Table i, page 6, clay, silt, sand and granite have different thermal conductivities, densities and specific heats. • The soil's moisture content. With its large C P , water has an effect on soil's thermal diffusivity. Moreover, water can induce seasonal variations of soil's proprieties. • Groundwater flows, which improve the soil's thermal reload through convective heat transfers. Today, we use three major methods to estimate these parameters. First, [START_REF] Mogensen | Fluid to duct wall heat transfer in duct system heat storages[END_REF] presents the Temperature Response Test (TRT). Recently, [START_REF] Mattsson | Advanced compact device for the in situ determination of geothermal characteristics of soils[END_REF] improved this method by reducing the size of the required equipment. For the TRT, we build a GHE into which we inject a known heat power. We then measure the temperature at the tube's input and output to deduce the soil's diffusivity and the thermal resistance R b of the borehole. With these parameters, we can calculate the heat transfer with the soil. With the output temperature, we can also calculate the COP of the heat pump. Although this method is efficient, it is obviously not convenient in the framework of a previous field study in order to determine the soil's proprieties for designing a GHE. Its application is only interesting to compare simulation results with field observations. A third method is to measure thermal proprieties of small samples in a laboratory and establish tables of values for different types of soils. We deduce the value of the thermal conductivity of a certain soil according to those of its components. But this approach does not consider the variability of soils.

A third method concerns surface measurements. [START_REF] Bristow | A small multi-needle probe for measuring soil thermal properties, water content and electrical conductivity[END_REF] give an example of a small probe which measures temperature, moisture content, electrical conductivity and thermal conductivity in the first ten centimeters of the ground. However, this kind of measurement is useful for agricultural studies but cannot access the soil at depth.

PERSPECTIVES WITH INDIRECT GEOPHYSICAL MEASUREMENTS

So far, thermal proprieties of soils seem difficult to measure, principally because of their dependence on composition and moisture content of the soils. As a consequence, [START_REF] Singh | A generalised relationship between soil electrical and thermal resistivities[END_REF] propose to measure these proprieties indirectly. They use a relation between thermal resistivity and electrical resistivity, which is easy to measure. Indeed, the authors show experimentally that both of these resistivities (we note them R E and R T respectively) depend on moisture content S r according to the following laws:

log R E = C 1 +C 2 log S r and log R T = C 3 +C 4 log S r (9)
where C 1 , C 2 , C 3 and C 4 are constants which depend on the type of soil. [START_REF] Bouchez | Relation entre les résistivités électrique et thermique des sols[END_REF] have verified these relations during an experimental lesson at ENS Lyon. Despite inherent experimental difficulties for measuring thermal conductivity, these relations seem us reliable and we can then deduce a relation between thermal and electrical resistivities:

log R E R o E = C 2 C 4 log R T R o T (10)
where, formally, log R o E = C 1 and log R o T = C 3 . In our study, equation ( 10) is under-determinated: Measuring the electrical resistivity of a soil is not sufficient to infer its thermal resisitivity. Even though we know the coefficients C 1 and C 4 from a measurement of electrical conductivity vs. water content, the R o T parameter is still unknown and only accessible through thermal measurements.

However, [START_REF] Singh | A generalised relationship between soil electrical and thermal resistivities[END_REF] deduced from their observations a direct relation between R E and R T :

log(R E ) = C R log(R T ) (11)
where C R only depend on the sand and gravel size fractions F of the soil: C R = 1.34 + 0.0085 × F.

Plotting the graph of log(R T ) vs. log(R E ) with consideration for the uncertainties due to the measurement, [START_REF] Bouchez | Relation entre les résistivités électrique et thermique des sols[END_REF] have some doubts about the reliability of such a linear regression. However, this kind of indirect measurement is promising to estimate the soil's thermal proprieties in situ. Assuming that for a given soil, thermal and electrical conductivities are both essentially governed by water content, we can apply the geophysical prospection methods developed in the field of hydrogeophysics (Rubin and Hubbard (2005) [START_REF] Robinson | Advancing process-based watershed hydrological research using near-surface geophysics: a vision for, and review of, electrical and magnetic geophysical methods[END_REF]) to profile, map or even monitor the electrical resistivity (or dielectric permittivity) at depth. Among these methods, we can cite:

• ground-penetrating radar (GPR, e.g. [START_REF] Huisman | Measuring soil water content with ground penetrating radar: A review[END_REF], [START_REF] Loeffler | Ground penetrating radar measurements in a controlled vadose zone: Influence of the water content[END_REF]), which initially measured the dielectric permittivity. Present research tend to access to the electrical conductivity through amplitude versus offset approaches (Deeds and Bradford ( 2002)) or full-waveform inversion (Klotzsche et al. ( 2010)) but it requires significant calculation power.

• electrical resistivity tomography (ERT) which can provide hundred-meter-long electrical resistivity profiles, up to a few dozen meters deep (e.g. [START_REF] Rings | Soil water content monitoring on a dike model using electrical resistivity tomography[END_REF], [START_REF] Koch | Joint interpretation of hydrological and geophysical data: Electrical resistivity tomography results from a process hydrological research site in the black forest mountains, germany[END_REF], see figure v). ERT requires a specific multi-electrodes instrumentation but the software RES2DINV1 enables to run the inversion on personal computers with an user-friendly interface. • electromagnetic induction (EMI, e.g. [START_REF] Abdu | Geophysical imaging of watershed subsurface patterns and prediction of soil texture and water holding capacity[END_REF]) which can quickly provide maps of apparent electrical conductivity, resultant value from the resistivities distribution at depth. Recently, [START_REF] Lavoué | Electromagnetic induction calibration using apparent electrical conductivity modelling based on electrical resistivity tomography[END_REF] improved the calibration of this method by coupling EMI values with ERT data. This advance open the possibility to get profiles of local values at depth with EMI as well. • electric self-potential (SP) measurements which can be linked to groundwater flows through electrokinetic relations (e.g. [START_REF] Trique | Fluid flow near reservoir lakes inferred from the spatial and temporal analysis of the electric potential[END_REF]).

Electrical conductivity and dielectric permittivity can then be converted into water content thanks to Archie's law [START_REF] Archie | The electrical resistivity log as an aid in determining some reservoir characteristics[END_REF]) and Topp's equation [START_REF] Topp | Electromagnetic determination of soil water content: measurements in coaxial transmission lines[END_REF]) respectively. However, these empirical laws are likely to depend on soil's charasteristics and thus only provide a water content scale rather than reliable water content values. Nevertheless, it allows to have an overview of the spatial variations of electrical proprieties and water content and to plane time-monitoring to detect seasonal variations and groundwater circulation. In an ideal world, it could be applied to adapt the GHE dimensioning and its working along the year knowing the seasonal soil's characteristics.

CONCLUSION AND PERSPECTIVES

We presented the modeling of GHE that has been made for many years and by many researchers, both with analytical and numerical methods. We then investigated a practical question which is rarely found in the literature: the influence of the soil's thermal parameters on GHE simulations. We used the TRNSYS software for the design of GHE test systems. Indeed, we showed that the required length of tubes depends strongly on the values of conductivity, and less on density and heat capacity of the soil. As a consequence, an effort must be made in further research to improve in situ measurements of soil's thermal parameters in the framework of previous studies for dimensioning GCHP systems. To do so, we really believe in the possibilities that geophysical prospection offers with GPR, EMI or ERT methods. With these techniques, we can indeed have profiles of dielectric permittivity, electrical conductivity, and thus water content at depth. An interesting interpretation of these data would then be to convert them into thermal conductivity data. This interpretation would involve a relation between electrical resistivity, thermal resistivity and water content in the heterogeneous, porous, and multi-phasic media that are soils. We must still investigate this relation precisely to make it applicable.

Figure 1 - 1 .

 11 Figure 1-1. Schematic of cycles in a GSHP system in (a) cooling mode and (b) heating mode.

  Figure ii : Temperature profile at March 21 st in a homogeneous soil, with α = 10 -6 m 2 .s -1 , T av = 12 • C, ∆T = 7.5 • C and t min = 0 (i.e. 1 st of January is the coldest day).

  Figure iii : Influence of pipe length on the heat exchanges obtained from the TRNSYS GHE simulation component for extreme conductivities, namely light dry soil against literature's default (erroneous) value.

  Figure v : ERT profile after Lavoué et al. (2010), converted into water content through a compact Archie's law (θ = 2.43σ 0.48 ).

Table i :

 i Values of thermal conductivity, heat capacity and diffusivity of some buildings, rocks and soils materials, dry or water-saturated (sat.).

	Air	0.062	1.25 × 10 -3	50
	Water	0.60	4.19	0.14
	Concrete	0.92	2.11	0.44
	Granite	2.2	2.13	1.03
	Clay	0.2 (dry) to 1.6 (sat.) 0.3 (dry) to 3.2 (sat.) 0.50 (sat.) to 0.67 (dry)
	Silt	0.2 (dry) to 2.5 (sat.) 0.60 (dry) to 2.4 (sat.) 0.33 (dry) to 1.0 (sat.)
	Sand	0.3 (dry) to 3.2 (sat.) 1.0 (dry) to 2.4 (sat.)	0.30 (dry) to 1.3 (sat.)
	Gravel	0.3 (dry) to 3.3 (sat.) 1.2 (dry) to 2.4 (sat.)	0.25 (dry) to 1.4 (sat.)

Materials λ (W.m -1 .K -1 ) ρC P (MJ.K -1 .m -3 ) α (.10 -6 m 2 .s -1 )

Table ii :

 ii Complex soil's parameters from[START_REF] Giardina | Evaluation of ground coupled heat pumps for the state of wisconsin[END_REF].

	Soil's type	Conductivity Diffusivity	Density	Specific heat
		(W /m/K)	(10 -7 m 2 /s) (10 3 kg/m 3 ) (10 2 J/kg/K)
	Heavy saturated soil	2.43	9.03	3.20	8.38
	Heavy damp soil	1.30	6.45	2.10	9.64
	Heavy dry soil	0.87	5.16	2.00	8.38
	Light damp soil	0.87	5.16	1.60	10.48
	Light dry soil	0.35	2.80	1.44	8.38

Table iii :

 iii Energy exchanged over a year, in typical Saint-Etienne climate external surface temperature, between the soil and an infinite load at 20 • C. J) (10 8 J) (10 8 J) (10 8 J) (10 8 J) (10 8 J)

		Heat	Cold	Heat	Cold	Heat	Cold
	(10 8 Pipe length	50 m	100 m	150 m
	Heavy saturated soil	0.80	-1.63	0.90	-2.08	0.90	-2.27
	Heavy damp soil	0.78	-1.06	0.95	-1.57	0.97	-1.84
	Heavy dry soil	0.79	-0.70	1.00	-1.15	1.07	-1.45
	Light damp soil	0.80	-0.70	1.00	-1.15	1.08	-1.51
	Light dry soil	0.96	-0.19	1.35	-0.38	1.69	-0.48
	m 2

TRNSYS homepage: http://sel.me.wisc.edu/trnsys/index.html

http://www.geoelectrical.com/index.php

Acknowledgements

We want to thank Jean-Jacques Roux and Kim Eui-Jong from the INSA of Lyon for sharing their experience in dealing with TRNSYS and more generally in modeling GHE.