Leonid Galtchouk

Sergey Pergamenshchikov

Adaptive sequential estimation for ergodic diffusion processes in quadratic metric. Part 1: Sharp non-asymptotic oracle inequalities

Keywords: primary 62G08; secondary 62G05, 62G20 2 adaptive estimation, asymptotic bounds, efficient estimation, nonparametric estimation, non-asymptotic estimation, oracle inequality, Pinsker's constant. 1

Adaptive sequential estimation for ergodic diffusion

processes in quadratic metric. Part 1: Sharp non-asymptotic oracle inequalities.

Introduction

Let (Ω, F, (F) t≥0 , P) be a filtered probability space on which the following stochastic differential equation is defined

dy t = S(y t ) dt + dw t , 0 ≤ t ≤ T , (1.1)
where (w t ) t≥0 is a scalar standard Wiener process, the initial value y 0 is a given constant and S(•) is an unknown function.

The problem is to estimate the function S from observations of (y t ) 0≤t≤T

in non-asymptotic setting, and to obtain sharp non-asymptotic bounds for a minimax quadratic risk. This problem is important for applications.

It seems that, for the first time, the problem of non-asymptotic parameter estimation for diffusion processes has been studied in [START_REF] Arato | On parameter estimation of a complex stationary gaussian process[END_REF] for wobling analysis of the axis of the equator. There, for a special diffusion process, the exact distribution of the ML-estimators of unknown parameters has been obtained

for any finite sample time T . Unfortunately, in the majority cases, when the sample time is finite, it is difficult to study classical estimators such as LS-estimators or ML-estimators since they are non-linear functionals of observations.

The paper [START_REF] Liptser | Statistics of a random process[END_REF] has shown that many difficulties in non-asymptotic parameter estimation for one-dimensional diffusion processes can be overcome by the sequential approach. It turns out that the sequential ML-estimator is more advantageous than the usual one. In particular, it is possible to calculate non-asymptotic bounds for quadratic risk in the sequential procedure.

By making use of the sequential approach non-asymptotic parameter estimation problems have been studied in [START_REF] Konev | Sequential Estimation of the Parameters of Random Processes with Continuous Time[END_REF], [START_REF] Galtchouk | On Sequential Estimation of Parameters in Continuous-Time Stochastic Regression[END_REF] for multidimensional diffusion processes and recently in [START_REF] Galtchouk | On Sequential Estimation of Parameters in Semimartingale Regression Models with Continuous Time Parameter[END_REF] for multidimensional continuous and discrete time semimartingales. In the paper [START_REF] Konev | On Truncated Sequential Estimation of the Parameters of Diffusion Processes[END_REF] a truncated sequential method has been developed for parameter estimation in diffusion processes.

The sequential approach to nonparametric minimax estimation problem of the drift coefficient in ergodic diffusion processes has been developed in [START_REF] Galtchouk | Estimateurs séquentiels à noyaux et aux polynômes locaux pour le problème d'estimation nonparamétrique de la dérive d'un processus de diffusion[END_REF]- [START_REF] Galtchouk | Nonparametric sequential minimax estimation of the drift in diffusion processes[END_REF], [START_REF] Galtchouk | Asymptotic efficient sequential kernel estimates of the drift coefficient in ergodic diffusion processes[END_REF]. The papers [START_REF] Galtchouk | Estimateurs séquentiels à noyaux et aux polynômes locaux pour le problème d'estimation nonparamétrique de la dérive d'un processus de diffusion[END_REF], [START_REF] Galtchouk | Sequential nonparametric adaptive estimation of the drift coefficient in diffusion processes[END_REF], [START_REF] Galtchouk | Nonparametric sequential minimax estimation of the drift in diffusion processes[END_REF] deal with sequential pointwise kernel estimators of the drift coefficient. A non-asymptotic upper bound was obtained for absolute error risks. The estimators yield also the optimal convergence rate as the sample time T → ∞. In the paper [START_REF] Galtchouk | Sequential nonparametric adaptive estimation of the drift coefficient in diffusion processes[END_REF] it is shown that this procedure is minimax and adaptive in the both cases when either the smoothness is known or unknown. The same type of the sequential kernel estimators is used in the paper [START_REF] Galtchouk | Nonparametric sequential estimation of the drift in diffusion via model selection[END_REF] for the nonparametric estimation in the L 2 -metric of the drift coefficient via model selection. A non-asymptotic upper bound for the quadratic risk is proved. The procedure is minimax and adaptive in the asymptotic setting as well. A sequential asymptotically efficient kernel estimator is constructed for pointwise drift estimation in [START_REF] Galtchouk | Asymptotic efficient sequential kernel estimates of the drift coefficient in ergodic diffusion processes[END_REF].

A new method using the concentration inequality for the nonlinear LSestimators analysis has been proposed in [START_REF] Birgé | Gaussian model selection[END_REF], [START_REF] Massart | A Non-asymptotic theory for model selection[END_REF] for homoscedastic nonparametric gaussian regression. A non-asymptotic upper bound for a quadratic risk (the Oracle inequality) was obtained which is best over the LS-estimators.

This method was developed in [START_REF] Fourdrinier | Improved selection model method for the regression with dependent noise[END_REF] for any estimators to regression models with depending noises for which a non-asymptotic Oracle inequality has been proved. Further development of this method was done recently in [START_REF] Galtchouk | Efficient adaptive nonparametric estimation in heteroscedastic regression models[END_REF] for heteroscedastic nonparametric regression models with a varying variance depending on unknown regression function. In the last paper weighted LSestimators were used to obtain the non-asymptotic Oracle inequality.

This paper deals at the first time with the non-asymptotic estimation of the drift coefficient S in adaptive setting for the quadratic risk

R( ŜT , S) = E S ŜT -S 2 , (1.2)
where ŜT is an estimator of S, S 2 = b a S 2 (x)dx , a, b are some real numbers, b -a ≥ 1. Here E S is the expectation with respect to the distribution law P S of the process (y t ) t≥0 given the drift S.

Our main goal is to construct by sequential methods an adaptive estimator Ŝ * of the function S in (1.1) and to show that the quadratic risk of this estimator satisfies a sharp non-asymptotic oracle inequality with respect to a family of weighted LS-estimators.

The sharp non-asymptotic oracle inequality with respect to a estimators family ( Ŝα , α ∈ A) means that the principal term in the upper bound for the quadratic risk has a coefficient which closes to one, i.e.

R( ŜT , S)

≤ (1 + D(ρ)) min α∈A R( Ŝα , S) + B T (ρ) T , (1.3)
where the function D(ρ) → 0 as ρ → 0 and B T (ρ) is some slowly varying function, i.e. for any δ > 0 and ρ > 0, lim

T →∞ B T (ρ) T δ = 0 . (1.4)
The goal of this paper is to obtain this inequality for model (1.1). To this end we pass to a discrete time regression model by making use of the truncated sequential procedure introduced in [START_REF] Galtchouk | Sequential nonparametric adaptive estimation of the drift coefficient in diffusion processes[END_REF], [START_REF] Galtchouk | Nonparametric sequential estimation of the drift in diffusion via model selection[END_REF] and [START_REF] Galtchouk | Asymptotic efficient sequential kernel estimates of the drift coefficient in ergodic diffusion processes[END_REF] and we come to the Γ-regression scheme, i.e. the sequential kernel estimator at the point x k on the set Γ has the following representation

Y k = S(x k ) + ζ k . (1.5)
sequence (ζ k ) k≥1 has a complex structure, i.e.

ζ k = σ k ξ k + δ k , (1.6)
where (σ k ) k≥1 is a sequence of observed random variables, (δ k ) k≥1 is a sequence of bounded random variables and (ξ k ) k≥1 is a sequence of i.i.d. random variables N (0, 1) which are independent of (σ k ) k≥1 .

In this paper we propose a new method to construct adaptive estimators that is based on the approach proposed in [START_REF] Galtchouk | Efficient adaptive nonparametric estimation in heteroscedastic regression models[END_REF] and on the penalisation method proposed in [START_REF] Baron | Risk bounds for model selection via penalisation[END_REF]. In our method we propose to add an penalty term with a small coefficient into the quadratic loss functional used to find an estimator.

This new functional criterion with the penalty term allows us to obtain a sharp oracle inequality and the above small coefficient provides that the principal term of the oracle inequality has a coefficient closes to one. The new criterion allows us to overcome difficulties bound up with depending noises appearing after the discrete Fourier transformation to model (1.5).

To estimate the function S in model (1.5) we make use of the estimator family ( Ŝα , α ∈ A), where Ŝα is a weighted least squares estimator with the Pinsker weights. For this family similarly to [START_REF] Galtchouk | Efficient adaptive nonparametric estimation in heteroscedastic regression models[END_REF], we construct a special selection rule, i.e. a random variable α * with values in A, for which we define the selection estimator as Ŝ * = Ŝα * . We show (see Theorem 3.2 ) that this selection rule for any 0 < ρ < 1/6 and any T ≥ 32 satisfies the oracle inequality (1.3) with

D(ρ) = 11ρ + 7ρ 2 + 3ρ 3 1 -6ρ (1.7)
and with some function B T (ρ) satisfying property (1.4).

We make use of the non-asymptotic analysis method from [START_REF] Galtchouk | Efficient adaptive nonparametric estimation in heteroscedastic regression models[END_REF] based on the non-asymptotic investigation proposed in [START_REF] Baron | Risk bounds for model selection via penalisation[END_REF] for the family of leastsquares estimators and developed in [START_REF] Fourdrinier | Improved selection model method for the regression with dependent noise[END_REF] for some other families of estimators.

In the second part of this investigation we prove the asymptotic efficiency of the proposed procedure that is we show that the oracle inequality (1.3) provides the minimal quadratic risk which coinsides with the Pinsker constant. Juste in this sense the inequality (1.3) is sharp.

The paper is organized as follows. In the next Section we give the definitions of functional classes and the passage from diffusion model (1.1) to a discrete time regression model. In Section 3 the sharp non-asymptotic oracle inequality is given for the regression model (Theorem 3.1). In Section 4 we prove Theorem 3.1 and Theorem 3.2. Appendix contains the proofs of some technique results.

Passage to a discrete time regression model

To obtain a good estimate of the function S, it is necessary to impose some conditions on the function S which are similar to the periodicity of the deterministic signal in the white noise model. One of conditions which is sufficient for this purpose is the assumption that the process (y t ) in (1.1) returns to any vicinity of each point x ∈ [a, b] infinite times. The ergodicity of (y t ) provides this property.

Let L > 1. We define the following functional class :

Σ L = {S ∈ Lip L (R) : sup |z|≤L |S(z)| ≤ L ; ∀|x| ≥ L , ∃ Ṡ(x) ∈ C such that -L ≤ Ṡ(x) ≤ -1/L} , (2.1)
where

Lip L (R) = f ∈ C(R) : sup x,y∈R |f (x) -f (y)| |x -y| ≤ L .
First of all, note that if S ∈ Σ L , then there exists the ergodic density

q(x) = q S (x) = exp{2 x 0 S(z)dz} +∞ -∞ exp{2 y 0 S(z)dz}dy (2.2)
(see,e.g., Gihman and Skorohod (1972), Ch.4, 18, Th2). It is easy to see that this density is uniformly bounded in the class (2.1), i.e.

q * = sup x∈R sup S∈Σ L q S (x) < +∞
and bounded away from zero on the interval [-b * , b * ], i.e.

q * = inf |x|≤b * inf S∈Σ L q S (x) > 0 , (2.3)
where b * = 1 + |a| + |b|. Moreover, we note that the functions from Σ L are uniformly bounded, i.e.

s * = sup a≤x≤b sup S∈Σ L S 2 (x) < ∞ . (2.4)
To obtain the oracle inequality we pass to discrete time regression model (1.5) by the same way as in [START_REF] Galtchouk | Nonparametric sequential estimation of the drift in diffusion via model selection[END_REF]. We start with the partition of the interval [a, b] by points (x k ) 1≤k≤n defined as

x k = a + k n (b -a) (2.5) with n = 2[(T -1)/2] + 1 ([x]
is the integer part of x). At each point x k we will estimate the function S by the sequential kernel estimator from [START_REF] Galtchouk | Nonparametric sequential estimation of the drift in diffusion via model selection[END_REF]- [START_REF] Galtchouk | Nonparametric sequential minimax estimation of the drift in diffusion processes[END_REF].

We fix some 0 < t 0 < T and we set

τ k = inf{t ≥ t 0 : t t 0 Q y s -x k h ds ≥ H k } (2.6)
and

S * k = 1 H k τ k t 0 Q y s -x k h dy s , (2.7)
where

Q(z) = 1 {|z|≤1} , h = (b -a)/(2n
) and H k is a positive threshold. From

(1.1) it is easy to obtain that

S * k = S(x k ) + ζ k .
The error term ζ k is represented as the following sum of the approximated term B k and the stochastic term:

ζ k = B k + 1 H k ξ k ,
where

B k = 1 H k τ k t 0 Q y s -x k h (S(y s ) -S(x k ))ds , ξ k = 1 √ H k τ k t 0 Q y s -x k h dw s .
Taking into account that the function S is lipschitzian, we obtain the upper bound for the approximated term as

|B k | ≤ L h . (2.8)
It is easy to see that the random variables (ξ k ) 1≤k≤n are i.i.d. normal N (0, 1).

Moreover, in [START_REF] Galtchouk | Asymptotic efficient sequential kernel estimates of the drift coefficient in ergodic diffusion processes[END_REF] it is established that the efficient kernel estimator of type (2.7) has the stochastic term distributed as N (0, 2T hq S (x k )). Therefore for the efficient estimation at each point x k by kernel estimator (2.7) we need to estimate the ergodic density (2.2) from the observations (y t ) 0≤t≤t 0 . We put

qT (x k ) = max{q(x k ) , T } , where T is positive, 0 < T < 1, q(x k ) = 1 2t 0 h t 0 0 Q y s -x k h ds .
(2.9)

Now we choose the threshold H k in (2.6)-(2.7) as

H k = (T -t 0 )(2q T (x k ) -2 T )h .
We suppose that for any T ≥ 32 the parameters t 0 = t 0 (T ) and T satisfy the following conditions

16 ≤ t 0 ≤ T 2 and √ 2 t -1/8 0 ≤ T ≤ 1 . (2.10) Moreover, lim T →∞ T = 0 , lim T →∞ T T t 0 (T ) = ∞ . (2.11)
Finally, we assume also for any δ > 0 and µ > 0, lim

T →∞
T µ e -δ √ t 0 = 0 .

(2.12)

For example, one can take, for T ≥ 32,

t 0 = max{min{ln 4 T , T /2} , 16} and T = √ 2 t -1/8 0 .
We set now

Γ = { max 1≤k≤n τ k ≤ T } and Y k = S * k 1 Γ . (2.13)
Therefore on the set Γ one has the discrete time regression model (1.5) with

δ l = B l and σ 2 l = n (T -t 0 )(q T (x l ) -2 T /2)(b -a) . (2.14)
To study properties of the set Γ we introduce the following functions

(c, L) = 2c e L 3 + 2q * q * + q * ∞ 0 e -z 2 /L+2Lz dz , K * (L) = L max y 2 0 , L 2 (1 + 2L 2 + 2L) 2 + y 2 0 + 1 and γ(c, L) = 1 54 2 (c, L)K * (L) . ( 2

.15)

In Appendix A.1 we show the following result.

Proposition 2.1. Suppose that the parameters t 0 and T satisfy the conditions (2.10)-(2.12). Then sup

S∈Σ L P(Γ c ) ≤ T 1 {q * <2 T } + 16T e -γ 1 √ t 0 1 {q * ≥2 T } := Π T (2.16)
with

γ 1 = γ(2, L) min(q 2 * , 1/8) .
Note that the conditions (2.10)-(2.12) imply directly that for any m > 0 lim

T →∞ T m Π T = 0 .
(2.17)

Oracle inequalities

In this section we consider the estimation problem for the following discrete time regression model. We suppose that on some set Γ ⊆ Ω we have the following nonparametric regression

Y l = S(x l ) + ζ l , 1 ≤ l ≤ n , (3.1)
where ζ l = σ l ξ l + δ l and n is odd. The points (x l ) 1≤l≤n are defined in

Ŝ -S 2 n = ( Ŝ -S, Ŝ -S) n = b -a n n l=1 ( Ŝ(x l ) -S(x l )) 2 .
We suppose that (ξ i ) 1≤i≤n is a sequence of i.i.d. random variables N (0, 1), (σ l ) 1≤l≤n is a sequence of observed positive random variables independent of (ξ i ) 1≤i≤n and bounded by some nonrandom constant

ς n ≥ 1, i.e. max 1≤l≤n σ 2 l ≤ ς n . (3.2)
Concerning the random sequence δ = (δ l ) 1≤l≤n we suppose that

E S δ 2 n < ∞ . (3.3)
Note that the model (1.5) is a particular case of (3.1) with Γ defined in (2.13) and δ l = B l . Indeed, in view of the inequality (2.6), we obtain that, for any function

S ∈ Lip L , max 1≤l≤n |B l | ≤ Lh = L(b -a) 2n .
Therefore, in this case

δ 2 n ≤ L 2 (b -a) 3 4n 2 . ( 3.4)
Moreover from (2.14), in view of condition (2.10), we obtain that

max 1≤l≤n σ 2 l ≤ 4 (b -a) T ≤ max 4 (b -a) T , 1 = ς n . (3.5)
We make use of the trigonometric basis (φ j ) j≥1 in L 2 [a, b] with

φ 1 = 1 √ b -a , φ j (x) = 2 b -a T r j (2π[j/2]v(x)) , j ≥ 2 , (3.6)
where the function T r j (x) = cos(x) for even j and T r j (x) = sin(x) for odd j; v(x) = (x -a)/(b -a). This basis is orthonormal for the empiric inner product generated by the sieve (2.5) , i.e. for any 1 ≤ i, j ≤ n,

(φ i , φ j ) n = b -a n n l=1 φ i (x l ) φ j (x l ) = Kr ij , (3.7)
where Kr ij is Kronecker's symbol.

By making use of this basis we apply the discrete Fourier transformation to (3.1) and we obtain the Fourier coefficients

θj,n = b -a n n l=1 Y l φ j (x l ) , θ j,n = b -a n n l=1 S(x l ) φ j (x l ) .
From (3.1) it follows directly that these Fourier coefficients satisfy the fol-

lowing equation θj,n = θ j,n + ζ j,n (3.8)
with

ζ j,n = b -a n ξ j,n + δ j,n ,
where

ξ j,n = b -a n n l=1 σ l ξ l φ j (x l ) and δ j,n = b -a n n l=1 δ l φ j (x l ) .
Notice that inequality (3.3) and the Bouniakovskii-Cauchy-Schwarz inequality imply that

|δ j,n | ≤ δ n φ j n = δ n . (3.9)
We estimate the function S on the sieve (2.5) by the weighted least squares estimator

Ŝλ (x l ) = n j=1 λ(j) θj,n φ j (x l ) 1 Γ , 1 ≤ l ≤ n , (3.10)
where the weight vector λ = (λ(1), . . . , λ(n)) belongs to some finite set Λ into R n , the prime denotes the transposition. We set for any x ∈ [a, b]

Ŝλ (x) = Ŝλ (x 1 )1 {a≤x≤x 1 } + n l=2 Ŝλ (x l )1 {x l-1 <x≤x l } . (3.11)
Let ν be the cardinal number of the set Λ. We suppose that the components of the weight vector λ = (λ(1), . . . , λ(n)) are such that 0 ≤ λ(j) ≤ 1. By setting λ 2 = (λ 2 (1), . . . , λ 2 (n)) we define the following sets

Λ 1 = {λ 2 , λ ∈ Λ} and Λ 2 = Λ ∪ Λ 1 . (3.12)
We set

χ j (λ) = 1 {λ j >0} , χ * = max λ∈Λ n j=1 χ j (λ) and µ n = χ * √ n . (3.13)
Moreover, for the basis functions (3.6) we shall use the following upper bound

λ * = max λ∈Λ 1 sup a≤x≤b | n j=1 λ(j)φ j (x)| , (3.14)
where φ j = (b -a)φ 2 j -1. Now we have to write a rule to choose a weight vector λ ∈ Λ. It is obviously, that the best way is to minimize with respect to λ the empiric squared error

Err n (λ) = Ŝλ -S 2 n ,
which in our case is equal to

Err n (λ) = n j=1 λ 2 (j) θ2 j,n -2 n j=1 λ(j) θj,n θ j,n + n j=1 θ 2 j,n . (3.15)
Since θ j,n are unknown, we need to replace the term θj,n θ j,n by some its estimator which we choose as

θj,n = θ2 j,n - b -a n s j,n with s j,n = b -a n n l=1 σ 2 l φ 2 j (x l ) . (3.16)
For this substitution in the empiric squared error one has to pay a penalty.

Finally, we define the cost function by the following way

J n (λ) = n j=1 λ 2 (j) θ2 j,n -2 n j=1 λ(j) θj,n + ρP n (λ) , (3.17)
where the penalty term we define as

P n (λ) = (b -a)|λ| 2 s n n , s n = 1 n n l=1 σ 2 l , |λ| 2 = n j=1 λ 2 (j)
and 0 < ρ < 1 is some positive constant which will be chosen later. We set λ = agrmin λ∈Λ J n (λ) (3.18) and define an estimator of S at the model (3.1) as

Ŝ * (x) = Ŝλ (x) for a ≤ x ≤ b . (3.19)
In the Oracle inequality we make use of the following function

Υ * n (ρ) = υ * 1 (ρ)ς n (1 + µ n )ν + υ * 2 (ρ)ς n λ * + υ * 3 (ρ)χ * δ * n (3.20) with δ * n = nE S 1 Γ δ 2 n , υ * 1 (ρ) = 16(b -a)(6 + 3ρ) ρ , υ * 2 (ρ) = 4ρ(b -a) 3 and υ * 3 (ρ) = 4(6 + ρ) ρ .
Now with the help of Υ * n (ρ) we introduce the function which describes the second term in the upper bound of the Oracle inequality :

Ψ n (ρ) = ψ 1 (ρ)Υ * n (ρ) + ψ 2 (ρ)ς n (2ν + ρ 2 λ * ) + ψ 3 (ρ)χ * δ * n , (3.21)
where

ψ 1 (ρ) = 1 -2ρ 2ρ(1 -6ρ) , ψ 2 (ρ) = b -a 2ρ(1 -6ρ)
and ψ 3 (ρ) = 4 + ρ 2ρ(1 -6ρ) .

Theorem 3.1. For any n ≥ 3 and 0 < ρ < 1/6, the estimator Ŝ * satisfies the following Oracle inequality

E S Ŝ * -S 2 n ≤ 1 + 3ρ 1 -6ρ min λ∈Λ E S Ŝλ -S 2 n + Ψ n (ρ) n + κ n (ρ) + ρ 1 -6ρ S 2 n P S (Γ c ) , (3.22)
where

κ n (ρ) = ( √ 3 + 1) ρ 1 -6ρ µ * n and µ * n = (b -a)ς n χ * n .
Remark 3.1. It is clear that to obtain the asymptotically optimal variance in the case of unknown regularity one can take, for example, the parameter ρ going to zero like

ρ = ρ n = O 1 ln γ n (3.23)
for some γ > 0.

Now we consider the adaptive setting, i.e. we suppose the function S belongs to the Sobolev space

W k r = {S ∈ Σ L ∩ C k-1 ([a, b]) : k j=0 S (j) 2 ≤ r} ,
where the regularity exposant k ≥ 1 and the Sobolev radius r > 0 are unknown.

We consider the weight family Λ introduced in [START_REF] Galtchouk | Efficient adaptive nonparametric estimation in heteroscedastic regression models[END_REF]. For any 0 < ε < 1 we define the set

A ε = {1, . . . , k * } × {t 1 , . . . , t m } ,
where

k * = [1/ √ ε], t i = iε, m = [1/ε 2 ]
and we take ε = 1/ ln n for any n ≥ 3.

For any α = (β, t) ∈ A ε we define the weight vector λ α = (λ α (1), . . . , λ α (n))

with

λ α (j) =      1 , for 1 ≤ j ≤ j 0 , 1 -(j/ω α ) β + , for j 0 < j ≤ n , (3.24)
where

j 0 = j 0 (α) = [ω α / ln n] + 1, ω α = (A β t n) 1/(2β+1) and A β = (b -a) 2β+1 (β + 1)(2β + 1) π 2β β .
Therefore,

Λ = {λ α , α ∈ A ε } (3.25) and ν = k * m ≤ ln 5/2 n. Note that in this case n j=1 χ j (λ) = n j=1 1 {λ α (j)>0} ≤ ω α .
Therefore,

χ * ≤ max α ω α ≤ (b -a)π -2/3 (6n ln n) 1/3 .

Therefore, for any

α ∈ A ε sup a≤x≤b | n j=1 λ 2 α (j)φ j (x)| ≤ 2 + 2 • 2 k * +1 + 2 2k * +1 . lim n→∞ λ * n δ = 0 .
Moreover, if we choose the parameter ρ as in (3.23) and if we assume as in For this procedure we show the following Oracle inequality.

Theorem 3.2. Assume that S ∈ W k r with unknown parameters r > 0 and k ≥ 1. Then the procedure Ŝ * with n = 2[(T -1)/2] + 1, for any T ≥ 32 and 0 < ρ < 1/6, satisfies the following inequality

R( Ŝ * , S) ≤ (1 + ρ) 2 (1 + 3ρ) 1 -6ρ min α∈A ε R( Ŝα , S) + B T (ρ) n , (3.28)
where

B T (ρ) = (1 + ρ)(Ψ n (ρ) + g T ) + (1 + ρ) 2 L 2 (b -a) 3 ρ(1 -6ρ)n with g T = n(κ n (ρ) + s * (b -a)) Π T .
Moreover, for any 0 < ρ < 1/6 and any δ > 0, lim

T →∞ B T (ρ) T δ = 0 . (3.29)
Remark 3.2. Note that to obtain the inequality (3.28) we make use of the method proposed in [START_REF] Galtchouk | Efficient adaptive nonparametric estimation in heteroscedastic regression models[END_REF] and based on the nonasymptotic approach from [START_REF] Baron | Risk bounds for model selection via penalisation[END_REF] and [START_REF] Fourdrinier | Improved selection model method for the regression with dependent noise[END_REF].

Remark 3.3. In Part 2 of this research we will show the obtained oracle inequality is sharp in the sensee that it provides the asymptotic efficiency of the estimator Ŝα * . Moreover, we will find the Pinsker constant and will prove that the asymptotic quadratic risk of the above estimator coincides with this constant.

Proofs

Proof of Theorem 3.1

First of all, note that on the set Γ we can represent the empiric squared error

Err n (λ) by the following way

Err n (λ) = J n (λ) + 2 n j=1 λ(j)θ j,n + S 2 n -ρ P n (λ) (4.1)
with θ j,n = θj,n -θ j,n θj,n . From (3.8) we find that

θ j,n = θ j,n ζ j,n + b -a n ξj,n + 2 b -a n ξ j,n δ j,n + δ 2 j,n ,
where ξj,n = ξ 2 j,n -s j,n . We set

∆(λ) = b -a n n j=1 λ(j) ξj,n , ∆(λ) = 2 b -a n n j=1 λ(j) ξ j,n δ j,n . (4.2)
We remind that for any > 0 and for any x, y

2xy ≤ εx 2 + ε -1 y 2 . (4.3)
Therefore, for some 0 < < 1 with help inequality (3.9) we can estimate the last term as

|∆(λ)| ≤ b -a n n j=1 λ 2 (j) ξ 2 j,n + 1 n j=1 χ j δ 2 j,n ≤ b -a n n j=1 λ 2 (j) s j,n + ∆(λ 2 ) + 1 χ * δ 2 n ,
where the vector λ 2 is defined in (3.12). Taking into account that

n j=1 λ 2 (j) s j,n -|λ| 2 s n = 1 n n l=1 σ 2 l n j=1 λ 2 (j) φ j (x l ) ≤ ς n λ * (4.4)
we get the following upper bound for ∆(λ), i.e.,

|∆(λ)| ≤ P n (λ) + ∆(λ 2 ) + (b -a)ς n λ * n + 1 χ * δ 2 n . (4.5)
Now putting

M (λ) = n j=1 λ(j) θ j,n ζ j,n and ∆ 1 (λ) = n j=1 λ(j) δ 2 j,n (4.6)
we rewrite (4.1) by the following way

Err n (λ) = J n (λ) + 2M (λ) + 2 ∆(λ) + 2∆(λ) + 2∆ 1 (λ) + S 2 n -ρ P n (λ) (4.7)
Note now that, we can represent ξj,n as

ξj,n = 2(b -a) n l=2 u j,l ξ l + b -a n n l=1 σ 2 l φ 2 j (x l ) ξl = 2U j,n + V j,n , (4.8)
where ξl = ξ 2 l -1 and

u j,l = 1 n σ l φ j (x l ) l-1 r=1 σ r φ j (x r )ξ r .
To study the function ∆(λ) we set

U (λ) = 1 √ s n n j=1 λ(j) U j,n 1 {s n >0} and V (λ) = n j=1 λ(j) V j,n ,
where λ(j) = λ(j)/|λ|. Therefore, we can represent ∆(λ) as

∆(λ) = 2 b -a n U (λ) P n (λ) + b -a n V (λ) .
In Appendix A.2 we show that sup

λ∈R n E S U 2 (λ) ≤ 2ς n . (4.9)
Therefore, in view of definition (3.12) we find that

E S sup λ∈Λ 2 U 2 (λ) ≤ λ∈Λ 2 E S U 2 (λ) ≤ 4ς n ν . (4.10)
Moreover, taking into account that

E S V 2 j,n ≤ 4σ 2 * n ,
we get that sup

λ∈Λ 2 E S |V (λ)| ≤ 2 ς n µ n
and similarly to (4.10) we obtain that

E S sup λ∈Λ 2 |V (λ)| ≤ 4 ς n µ n ν .
Therefore, if we set

∆ * = sup λ∈Λ 2 U 2 (λ) + sup λ∈Λ 2 |V (λ)| , 20
then we get that

E S ∆ * ≤ 4 ς n (1 + µ n )ν . (4.11)
Moreover, taking into account that P n (λ 2 ) ≤ P n (λ) we obtain that for any λ ∈ Λ 2 and 0 < < 1

| ∆(λ)| ≤ P n (υ λ ) + (b -a)(1 + ) n ∆ * (4.12)
with υ λ = λ for λ ∈ Λ and υ λ = ( λ(1), . . . , λ(1)) for λ ∈ Λ 1 . We remind that the set Λ 1 is defined in (3.12). Therefore, from (4.5) replacing by 2 we can get for any λ ∈ Λ the following upper bound for ∆(λ)

|∆(λ)| ≤ 2 P n (λ) + (b -a)(1 + ) n ∆ * + (b -a)ς n λ * n + b -a χ * δ 2 n .
Moreover, directly from (3.9) we find that

sup λ∈Λ |∆ 1 (λ)| ≤ χ * δ 2 n .
Therefore, setting

∆ 2 (λ) = 2 ∆(λ) + 2∆(λ) + 2∆ 1 (λ) ,
we can estimate this term for any λ ∈ Λ as

|∆ 2 (λ)| ≤ 6 P n (λ) + 2(b -a)(1 + 3 ) n ∆ * + 2 (b -a) n ς n λ * + 2(1 + ) χ * δ 2 n . (4.13)
Now from (4.1) we obtain that for some fixed λ 0 from Λ

Err n ( λ) -Err n (λ 0 ) = J( λ) -J(λ 0 ) + 2 M ( θ) + ∆ 2 ( λ) -ρP n ( λ) -∆ 2 (λ 0 ) + ρP n (λ 0 ) ,
where θ = λ -λ 0 . Therefore by definition of λ in (3.18) and putting = ρ/6

we obtain that on the set Γ

Err n ( λ) ≤ Err n (λ 0 ) + 2 M ( θ) + Υ n (ρ) + 2ρP n (λ 0 ) , (4.14)
where

Υ n (ρ) = 4(b -a)(6 + 3ρ) nρ ∆ * + 4ρ(b -a) 3n ς n λ * + 4(6 + ρ) ρ χ * δ 2 n .
From (4.11) it follows that

E S Υ n (ρ) ≤ Υ * n (ρ) n , (4.15)
where the function Υ * n (ρ) is defined in (3.20). Now we study the second term in the right-hand part of the inequality (4.14). Firstly for any weght vector λ ∈ Λ we set ϑ = λ -λ 0 . Then we decompose this term as

M (ϑ) = √ b -a Z(ϑ) + N (ϑ) (4.16)
with

Z(ϑ) = n -1/2 n j=1
ϑ(j) θ j,n ξ j,n and N (ϑ) = n j=1 ϑ(j) θ j,n δ j,n .

We define now the weighted discrete Fourier transformation of S, i.e. we set

S λ = n j=1 λ(j) θ j,n φ j . (4.17)
It is easy to see that in this case

E S Z 2 (ϑ) = b -a n 2 n l=1 E S σ 2 l   n j=1 ϑ(j)θ j,n φ j (x l )   2 ≤ ς n d(ϑ) , (4.18)
where

d(ϑ) = S ϑ 2 n n .
Moreover, note that we can represent ϑ(j) = λ(j) -λ 0 (j) as On the set Γ we have

ϑ(j) = ϑ(j)χ * j with χ * j = max(1 {λ(j)>0} , 1 {λ 0 (j)>0} ) .
S ϑ 2 n -Ŝϑ 2 n = n j=1 ϑ 2 (j)(θ 2 j,n -θ2 j,n ) ≤ -2 n j=1 ϑ 2 (j) θ j,n ζ j,n = -2 √ b -aZ 1 (ϑ) -2N 1 (ϑ) , (4.22)
where

Z 1 (ϑ) = 1 √ n n j=1 ϑ 2 (j)θ j,n ξ j,n and N 1 (ϑ) = n j=1 ϑ 2 (j) θ j,n δ j,n .
By the same way, taking into account that |ϑ(j)| ≤ 1, we find that on the set Γ

E S Z 2 1 (ϑ) ≤ ς n d(ϑ) .
Similarly, we set

Z * 1 = sup ϑ∈Θ Z 2 1 (ϑ) d(ϑ)
and we obtain that on the set Γ

E S Z * 1 ≤ νς n . (4.23)
Similarly to (4.17) we estimate the second term in (4.22) as

2N 1 (ϑ) ≤ ρ S ϑ 2 n + 2χ * δ 2 n ρ .
Therefore, on the set Γ

S ϑ 2 n ≤ Ŝϑ 2 n + 2ρ S ϑ 2 n + (b -a)Z * 1 nρ + 2χ * δ 2 n ρ , i.e. S ϑ 2 n ≤ 1 1 -2ρ Ŝϑ 2 n + 1 (1 -2ρ)ρ (b -a)Z * 1 n + 2χ * δ 2 n . ( 4
(ϑ) ≤ 2ρ 1 -2ρ Ŝϑ 2 n + 1 ρ(1 -2ρ) (b -a)Z * 2 n + 2χ * δ 2 n ≤ 4ρ(Err n (λ) + Err n (λ 0 )) 1 -2ρ + 1 ρ(1 -2ρ) (b -a)Z * 2 n + 2χ * δ 2 n . Therefore (4.14) implies that Err n ( λ)1 Γ ≤ 1 + 2ρ 1 -6ρ Err n (λ 0 )1 Γ + 1 -2ρ 1 -6ρ Υ n (ρ)1 Γ + 1 ρ(1 -6ρ) (b -a)Z * 2 n + 2χ * δ 2 n 1 Γ + 2ρ(1 -2ρ) 1 -6ρ P n (λ 0 )1 Γ , E S Err n ( λ)1 Γ ≤ 1 + 2ρ 1 -6ρ E S Err n (λ 0 )1 Γ + 1 -2ρ n(1 -6ρ) Υ * n (ρ) + 2νς n (b -a) + 2χ * δ * n nρ(1 -6ρ) + ρ(1 -2ρ) 1 -6ρ E S 1 Γ P n (λ 0 ) .
By applying here Lemma A.1 with = 2ρ we get that

E S Err n ( λ)1 Γ ≤ 1 + 3ρ 1 -6ρ E S Err n (λ 0 )1 Γ + Ψ n (ρ) n + κ n (ρ) P S (Γ c ) ,
where the functions Ψ n (

-S 2 ≤ (1 + ρ) Ŝ -S 2 n + (1 + ρ -1 )L 2 (b -a) 3 1 n 2 and Ŝ -S 2 n ≤ (1 + ρ) Ŝ -S 2 + (1 + ρ -1 )L 2 (b -a) 3 1 n 2 .
Moreover, inequality (2.4) implies

S 2 n ≤ s * (b -a)
for any S from Σ L . Therefore inequality (3.22) implies directly (3.28). Hence Theorem 3.2.

A Appendix

A.1 Proof of Proposition 2.1

Fisrt note that for this proposition we make use the consentration inequality from [START_REF] Galtchouk | Uniform concentration inequality for ergodic diffusion processes[END_REF]. We remind it. (ψ(y t ) -M S (ψ))dt .

We will denote by ∆ t (ψ) = ∆ 0,t (ψ). Note that now by Theorem 3.2 from where

ψ h,k (y) = 1 h Q y -x k h .
It is easy to see that this function satisfies condition (A.1) with c = 2.

Moreover, notice now that we can represent the integral

( 2 . 5 )

 25 with b -a ≥ 1. The function S is unknown and to be estimated from observations Y 1 , . . . , Y n . The quality of an estimator Ŝ will be measured by the empiric squared error

  the sequence (δ l ) 1≤l≤n satisfies inequality (3.4) we deduce that lim n→∞ χ * δ * n = 0 . (3.26) By Lemma 6.2 from [11] we get for any integer m ≥ 0 and a ≤ x ≤ b sup x) ≤ 2 m+1 .

( 3 . 5 )

 35 that the parameter ς n is slowly increasing in n, i.e. for any δ > 0 lim n→∞ ς n n δ = 0 we obtain that for any δ > 0 lim n→∞ Ψ n (ρ) n δ = 0 . (3.27) Now we consider the estimation problem (1.1) via model (1.5). By denoting Ŝα = Ŝλ α We make use of estimating procedure (3.18)-(3.19) with the weight coefficients set (3.25), i.e. denoting Ŝα = Ŝλ α we set Ŝ * = Ŝα with α = agrmin α∈A ε J n (λ α ) .

Therefore, by inequality ( 4 . 3 )

 43 with = ρ, we can estimate N (ϑ) as follows 2N (ϑ) = 2 n j=1 ϑ(j) θ j,n χ * j δ j,n from (4.18) it follows that E S Z * ≤ νς n . (4.21) Now we estimate the first term in the right-hand part of the inequality (4.20).

  For any c > 0 and for any function ψ : R → R

[ 13 ]E S e γ∆ 2 T

 132 for any function ψ(•) satisfying (A.1) we havesup (ψ) ≤ 4 , (A.2)where the constant γ = γ(c, L) is defined by(2.15). Therefore taking into account the inequality|∆ t 0 ,T (ψ)| ≤ |∆ T (ψ)| + |∆ t 0 (ψ)| ,we obtain that for any µ > 0 and for any ψ satisfying (A.1)sup 1≤t 0 ≤T sup S∈Σ L P S (|∆ t 0 ,T (ψ)| ≥ 2µ) ≤ 8 e -γµ 2 . (A.3)Now we estimate the probability of Γ c asP S (Γ c ) ≤ n k=1 P S (τ k > T ) .By the definition of τ k in (2.6) one hasP S (τ k > T ) = P S T t 0 ψ h,k (y t ) dt < (T -t 0 )(2q T (x k ) -2 T ) ,

1 - 1 qλ 2 8 )-a n n l=1 σ 2 l φ 2 jλ 2 λ 2 λ 2 2 P

 112822222 k (y t ) dt = (T -t 0 ) M S (ψ h,k ) + √ T ∆ t 0 ,T (ψ h,k ) ,whereM S (ψ h,k ) = +∞ -∞ ψ h,k (y) q S (y) dy = S (x k + hz) dz .Taking into account that2q T (x k ) -M S (ψ h,k ) = 1 √ t 0 ∆ t 0 (ψ h,k ) (A.4)we obtainP S (τ k > T ) ≤ P S (q(x k ) < T ) + P S ∆ t 0 ,T (ψ h,k ) < T -t 0 √ T 1 √ t 0 ∆ t 0 (ψ h,k ) -)(λ(j) -1)θ j,n ξ j,n and I 2 = n j=1 (j)δ j,n ξ j,nwe get on the set Γ the following lower bound for the impirical riskErr n (λ) ≥ b -Notice that the random variable ξ j,n ∼ N (0, s j,n ) (conditionally with respect to σ 1 , . . . , σ n ) withs j,n = b (x l ) ≤ ς n .ThusE S ξ 4 j,n ≤ 3σ 2 * . (A.9)Let us consider the first term in (A.8). We haveE S 1 Γ n j=1 (j) ξ 2 j,n = E S n j=1 (j) s j,n -E S 1 Γ c n j=1 (j) ξ 2 j,n .Therefore, by inequalities (4.4) and (A.9) we get thatE S 1 Γ n j=1 λ 2 (j) ξ 2 j,n ≥ |λ| 2 E S s n 1 Γ -ς n λ * -√ 3ς n χ * P S (Γ c ) . (A.10)Moreover, taking into account that E S I 1 = 0 we estimate E S I 1 1 Γ as|E S I 1 1 Γ | = |E S I 1 1 Γ c | S (Γ c ) ≤ S n √ ς n χ * P S (Γ c ) . (A.11)

The second author is partially supported by the RFFI-Grant 04-01-00855.

Moreover, thanks to the conditions (2.10), the last probability can be estimated as follows :

Finally we obtain that

Let us estimate now the first probability in the left hand of this inequality.

Taking into account that M S (ψ h,k ) ≥ 2q * we get

Therefore if T ≤ q * /2 then the inequality (A.2) implies

Now, by applying in (A.5) the inequalities (A.3) and (A.6) we obtain that for T ≤ q * /2

where γ 1 is defined in (2.16). Therefore taking into account that n ≤ T we obtain the following upper bound

Hence the condition (2.11) implies Proposition 2.1.

A.2 Proof of Inequality (4.9)

By putting α l = n j=1 λ(j)u j,l and taking into account that the random variables (ξ k ) 1≤k≤n are independent of (σ k ) 1≤k≤n we obtain that

where

Therefore the orthonormality of the functions (φ j ) implies that

Now by making use of this inequality in (A.7) we get (4.9).

A.3 Technical lemma

Lemma A.1. For any n ≥ 1 and 0 < < 1

where µ * n is defined in (3.22).

Proof. Indeed, by the definition of Err n (λ) on the set Γ we have

.

(4.3), i.e. for 0

Therefore from (A.8) we get that

diffusion processes.


Tiejun Wang
email: twang3@unl.edu

Xunhong Chen

Anh Minh Tang

Yu-Jun Cui

On the Use of the Similar Media Concept for Scaling Soil Air Permeability

Keywords: Similar Media Concept, Scaling Factor, Soil Air Permeability, Air Permeability Model

)-KA, were selected to test this method using two independent datasets.

For the first dataset that included k a measured for river sediments, all of the four models were able to delineate the spatial distribution of k a with a reference curve of k a and a set of scaling factors. Specifically, the MQ model gave the least satisfactory results due to the less flexibility of its form; whereas, there were no significant differences in the performances for the HT, BC, and KA models. For the second dataset that contained k a measured for agricultural soils, the overall performance of the four models for scaling k a deteriorated, largely due to the alterations in the microscopic structures of soil samples caused by repacking and compression of soil samples.

Nonetheless, as the first attempt, this study shows the viability of using the similar media concept for scaling k a . The merit of this method resides in the fact that the spatial variations of moisture conditions and soil properties can be simultaneously included for analyzing the spatial distribution of k a . With a reference curve of k a and the distribution of scaling factors, this method would be particularly suitable for modeling subsurface gas transport.

Introduction

As one of the key factors that control subsurface gas transport and exchange of gas across the soil-atmosphere interface, soil air permeability (k a ) describes the ability of a soil to transmit gas.

The rising interest in k a is manifested by its broad range of applications in various fields, including greenhouse gas emission [START_REF] Ball | The influence of gas transport and porosity on methane oxidation in soils[END_REF][START_REF] Conen | A re-examination of closed flux chamber methods for the measurement of trace gas emissions from soils to the atmosphere[END_REF], landfills [START_REF] Jain | Air permeability of waste in a municipal solid waste landfill[END_REF][START_REF] Wu | Field air permeability and hydraulic conductivity of landfilled municipal solid waste in China[END_REF], soil vapor extraction systems (Poulsen et la., 1998;[START_REF] Farhan | Interaction of soil air permeability and soil vapor extraction[END_REF], crop growth [START_REF] Lipiec | Quantification of compaction effects on soil physical properties and crop growth[END_REF][START_REF] Barrios | Fallow management for soil fertility recovery in tropical Andean agroecosystems in Colombia[END_REF]. Moreover, due to easy operation and cost effectiveness, k a measured near field capacity has been used to predict soil saturated hydraulic conductivity [START_REF] Loll | Predicting saturated hydraulic conductivity from air permeability: Application in stochastic water infiltration modeling[END_REF][START_REF] Iversen | Air and water permeability in differently textured soils at two measurement scales[END_REF].

Under natural conditions, k a is affected by a number of soil factors (e.g., air-filled porosity and pore size distribution) and moisture conditions, all of which show various degrees of spatial variations. As such, field measured k a exhibits significant spatial variability [START_REF] Iversen | Field application of a portable air permeameter to characterize spatial variability in air and water permeability[END_REF].

For application and modeling purposes, it is thus crucial to evaluate the spatial distribution of k a and its controlling factors. However, compared to relevant researches on soil hydraulic properties, only few studies are available on the spatial distribution of k a , which nonetheless provided valuable insights into the understanding of the spatial pattern of k a [START_REF] Poulsen | Spatial and temporal dynamics of air permeability in a constructed field[END_REF][START_REF] Iversen | Field application of a portable air permeameter to characterize spatial variability in air and water permeability[END_REF][START_REF] Iversen | Runoff modelling at two field slopes: Use of in situ measurements of air permeability to characterize spatial variability of saturated hydraulic conductivity[END_REF]. Geostatistical techniques were mainly used in previous studies to analyze the spatial pattern of k a . Based on the results of variograms, [START_REF] Iversen | Field application of a portable air permeameter to characterize spatial variability in air and water permeability[END_REF] found that the range of k a for sandy soils was larger than the one for a loamy soil, probably due to the difference in the depositional processes of those two types of soils. Although geostatistical techniques have been proven to be powerful tools for investigating naturally occurred phenomena, there are certain shortcomings for those techniques as pointed out by [START_REF] Henley | Geostatisticscracks in the foundations?[END_REF]. Most notably, the underlying processes associated with studied targets cannot be explicitly considered in geostatistical techniques, which rather rely on statistical models for examining the spatial correlations of the targets.

With respect to k a , the main issue of using geostatistical techniques stems from the impact of soil moisture on k a . It has been well known that k a is highly dependent on soil moisture levels, but the spatial distribution of soil moisture cannot be directly included in the geostatistical analyses of k a . To some degree, this is analogous to study the spatial distribution of soil unsaturated hydraulic conductivity without specifying moisture conditions. Therefore, precautions were usually taken in previous field studies on the spatial distribution of k a . As conjectured by [START_REF] Iversen | Field application of a portable air permeameter to characterize spatial variability in air and water permeability[END_REF], when soil moisture contents reach field capacity, the air flow takes place in the majority of soil pores; thus, the impact of moisture on k a can be neglected with moisture contents near and below field capacity. Although the assumption made by [START_REF] Iversen | Field application of a portable air permeameter to characterize spatial variability in air and water permeability[END_REF] is useful for studying the spatial distribution of k a under dry conditions, it may fail at regions with wet climates or shallow groundwater tables. For application purposes (e.g., landfills and soil vapor extraction systems), the inclusion of moisture in analyzing the spatial distribution of k a is also inevitable. Therefore, it is desirable to seek alternative methods to assess the spatial distribution of k a under the influence of soil moisture.

Along the line of above thinking, the scaling scheme based on the similar media concept may provide a promising approach to investigating the spatial distribution of k a with the consideration of moisture conditions. First introduced by [START_REF] Miller | Physical theory for capillary flow phenomena[END_REF], the similar media concept is based on the assumption that the internal geometry of similar media only differs by microscopic length scales that can be characterized by scaling factors. The purpose of this scaling approach is to coalesce a range of functional relationships into a single curve through scaling factors that depict the spatial distribution of those functional relationships. More specifically, scaling factors are used to relate soil properties at a given location to the mean properties at a reference location, which are invariant of moisture conditions. This scaling method has been widely used in soil hydrology to characterize spatial variability of soil hydraulic properties and soil hydraulic functions with associated model parameters [START_REF] Warrick | Scaling field measured soil hydraulic properties using a similar media concept[END_REF][START_REF] Hopmans | A comparison of various methods to scale soil hydraulic properties[END_REF][START_REF] Shouse | Scaling of near-saturated hydraulic conductivity measured using disc infiltrometers[END_REF][START_REF] Zavattaro | Use of similar media scaling to characterize spatial dependence of near-saturated hydraulic conductivity[END_REF][START_REF] Hendrayanto | Scaling hydraulic properties of forest soils[END_REF][START_REF] Tuli | Simultaneous scaling of soil water retention and unsaturated hydraulic conductivity functions assuming lognormal pore-size distribution[END_REF]. With a reference functional relationship and the distribution of scaling factors, this method is particularly suitable for modeling purposes [START_REF] Peck | Effect of spatial variability of soil hydraulic properties in water budget modeling[END_REF][START_REF] Kabat | The scaling characteristics of soil parameters: from plot scale heterogeneity to subgrid parameterization[END_REF][START_REF] Salvucci | Limiting relations between soil moisture and soil texture with implications for measured, modeled and remotely sensed estimates[END_REF][START_REF] Oliveira | Simulation of solute transport in a heterogeneous vadose zone describing the hydraulic properties using a multistep stochastic approach[END_REF]. Given the similarities between water flow and gas transport in soils, one can expect the viability of applying this scaling approach for investigating the spatial distribution of k a under the influence of soil moisture.

To our knowledge, this research was the first attempt to extend the use of the similar media concept for scaling k a . The main objective of this study was to examine the feasibility of this approach using two datasets collected from USA and France. Four air permeability models were selected to delineate the functional relationship between k a and saturation degree of air. The results of this study demonstrated the feasibility of using the similar media concept for scaling k a , which also opened the door for utilizing this method for simulating gas transport in soils.

Materials and Methods

Similar Media Concept in Soil Hydrology

The use of the similar media concept for scaling soil water retention and hydraulic conductivity curves is well documented in the literature (see the review by [START_REF] Vereecken | Upscaling hydraulic properties and soil water flow processes in heterogeneous soils: a review[END_REF]. So, only a brief overview is given here. Based on the similar media concept, it is assumed that the microscopic structures (e.g., tortuosity, and relative particle size and pore size distributions) of similar soils are identical and only differ by microscopic length scales that can be characterized by scaling factors [START_REF] Peck | Effect of spatial variability of soil hydraulic properties in water budget modeling[END_REF][START_REF] Warrick | Scaling field measured soil hydraulic properties using a similar media concept[END_REF]. The scaling factor (α) is thus defined by the ratio of the microscopic characteristic length of a soil (λ) to the characteristic length of a reference soil (λ m ):

(1) where i=1, 2, …, L is the location of the soil and the subscript m denotes the reference soil. By its definition, α is invariant of soil moisture conditions and only dependent upon the location of the soil. As such, the soil water retention curve at any location can be scaled to the reference water retention curve through α:

(2) where h is the soil matric potential, θ is the volumetric moisture content, and the subscript w denotes water. The scaling relationship of the soil hydraulic conductivity curve can be written as:

(3) where K w is the hydraulic conductivity. Due to the fact that soil porosity may vary across locations, instead of θ, the saturation degree of moisture (S w ) is usually used [START_REF] Warrick | Scaling field measured soil hydraulic properties using a similar media concept[END_REF]:

(4) where ϕ is the soil porosity.

To describe the reference curves of h m and K w,m , soil water retention and hydraulic conductivity models (e.g., Brooks-Corey model and van Genuchten model) have been used, although polynomial functions of h m (S w ) and K w,m (S w ) were also adopted [START_REF] Warrick | Scaling field measured soil hydraulic properties using a similar media concept[END_REF].

By optimizing the differences between measured and calculated h(S w ) and K w (S w ) across all the locations, the model parameters for the reference curves of h m and K w,m , and the scaling factor α a,i at each location can be obtained [START_REF] Hopmans | A comparison of various methods to scale soil hydraulic properties[END_REF][START_REF] Tuli | Simultaneous scaling of soil water retention and unsaturated hydraulic conductivity functions assuming lognormal pore-size distribution[END_REF].

Extension of the Similar Media Concept for Scaling Soil Air Permeability

Although the scaling theory based on the similar media concept has been widely used in soil hydrology, there is a surprising lack of studies on its application for assessing the spatial distribution of k a . Strictly speaking, soil water permeability (k w ) should be used in the scaling procedure (Eq. ( 3)), as regardless of fluid properties, soil permeability, whether it is k w or k a , represents the intrinsic properties of soils to transmit fluids and thus reflects the microscopic structures of soils. However, given that the fluid of interest remains under the same conditions (e.g., temperature), the fluid properties do not change across locations and conductivities (e.g., K w ) can be used in Eq. ( 3). Therefore, one can write a similar scaling relationship for k a based on Eq. ( 3):

(5)

where S a =(1-S w ) is the saturation degree of air and the subscript a denotes air. For the same reason for using S w in scaling h(S w ) and K w (S w ), S a was used in this study instead of the volumetric soil-air content. Note that α a is also invariant of soil moisture conditions and only dependent upon the location of the soil.

To delineate the reference relationship of k a,m , four air permeability models were selected [START_REF] Ghanbarian-Alavijeh | Comparison of the predictions of universal scaling of the saturation dependence of the air permeability with experiment[END_REF]. The first model was proposed by [START_REF] Millington | Transport in porous media[END_REF] (denoted as MQ hereafter):

(6)
where k o is the air permeability at the porosity ϕ.

The second model was developed by Hunt (2005) (HT), which is based on the continuum percolation theory: (7) where S a,t is the percolation threshold used by [START_REF] Ghanbarian-Alavijeh | Comparison of the predictions of universal scaling of the saturation dependence of the air permeability with experiment[END_REF]. Note that the HT model can be reduced to Eq. ( 6) by assuming S a,t =0.

The third model was based on the hydraulic conductivity function of [START_REF] Brooks | Hydraulic properties of porous media[END_REF]

) (BC): (8)
where γ is the pore size distribution index, and l is the tortuosity-pore connectivity factor and assumed to be 2 based on the percolation theory [START_REF] Hunt | Percolation Theory for Flow in Porous Media[END_REF].

The last model was taken from [START_REF] Kawamoto | Gas transport parameters in the vadose zone: Development and tests of power-law models for air permeability[END_REF]

(KA): (9)
where (1+3/η) represents the combined effects of tortuosity and connectivity of air-filled pores.

In [START_REF] Kawamoto | Gas transport parameters in the vadose zone: Development and tests of power-law models for air permeability[END_REF], the measured k a at -100 cm H 2 O was used instead of k o . Essentially, the MQ model is a simplified version of Eq. ( 9) with η=3.

In parallel to previous studies on scaling h(S w ) and K w (S w ), the scaling factors α a,i and the model parameters (given in Tables 2 and3) for describing the reference curves of k a,m were obtained by optimizing the differences between measured and calculated k a . The following objective function was used: (10) where M(i) is the total number of observed k a at location i, is the jth observed k a corresponding to S a,j at location i, and k a,m (S a,j ) is the calculated k a at S a,j from Eq. ( 6) to Eq.( 9).

The generalized reduced gradient method was used to optimize Eq. ( 10) with the geometric mean of the scaling factors assumed to be unity [START_REF] Hopmans | A comparison of various methods to scale soil hydraulic properties[END_REF][START_REF] Zavattaro | Use of similar media scaling to characterize spatial dependence of near-saturated hydraulic conductivity[END_REF][START_REF] Hendrayanto | Scaling hydraulic properties of forest soils[END_REF]:

(11)
To evaluate the performance of each air permeability model, the root mean square error (RMSE) was used:

(12

)
where n is the total number of observed k a across all the locations, o k is either observed k a,i (S a,j ) or scaled k a,i (S a,j ) (i.e., observed k a,i (S a,j )×α a,i 2 ), and e k is the corresponding k a,m (S a,j ) calculated from the reference curve of k a,m .

Datasets

Two independent datasets collected from USA and France were used in this study to test the feasibility of applying the similar media concept for scaling k a . The first dataset (denoted as Dataset I) included measured air permeability data using sediment samples taken from the channel of the Platte River near Kearney, NE (Figure 1). The streambed sediments at the site are moderately-sorted, and mainly consist of sands and gravels with very low organic matter contents (Table 1). A total of 25 samples with a length of 20 cm for each sample were collected using transparent polycarbonate tubes with an inside diameter of 5.1 cm. The tubes were vertically pressed to the depth of 20 cm into the sediments. The distance between neighboring sampling locations was 1 m. After taking out the sample from the river channel, the bottom of the tube was wrapped by a figberglass screening to prevent sediments from falling off. A falling head permeability test was carried out on site to determine the saturated hydraulic conductivity of the sediments inside the tube [START_REF] Chen | Statistical and geostatistical features of streambed hydraulic conductivities in the Platte River, Nebraska[END_REF]. Water temperatures were also recorded during the test to calculate water dynamic viscosity, which was used to calculate k w .

After the measurements of saturated hydraulic conductivities, the sediment samples were brought back to the laboratory and air dried at room temperatures. On Day 2 and every three days afterwards until Day 17, air permeability tests were performed using a portable air permeameter (Figure 2), which resulted in 6 pairs of data points of k a and S a for each tested sample. To measure k a , the polycarbonate tubes were connected to the air permeameter, which mainly consisted of an air tank, an air pressure regulator, a rotameter that measured the flow rate of air through the sediments, and a manometer that recorded the pressure difference across the length of the sediments. Darcy's Law was then used to calculate k a [START_REF] Kirkham | Field method for determination of air permeability of soil in its undisturbed state[END_REF]. Due to the slippage of gas molecules near the pore wall, the flow of gas in porous media may exhibit non-Darcian flow that leads to higher k a [START_REF] Klinkenberg | The permeability of porous media to liquids and gases[END_REF]. Therefore, additional tests were performed to examine the applicability of Darcy's Law by measuring the pressure difference across the sample length under varying flow rates. The linearity of the results proved the applicability of Darcy's Law and confirmed the air permeameter was operational for measuring k a . The merit of the air permeameter used in this study is that it allows to measure both saturated hydraulic conductivity and k a of the same sample without the need to remove the sediments from the permeameter for repacking. Thus, it minimizes the possible disturbance of the sample structure. After each air permeability test, samples were immediately weighed to determine moisture contents and corresponding volumetric soil-air contents. On Day 17, all of the samples were dried in the oven for 24 h at the temperature of 105 °C to determine bulk density and porosity. The results of the grain size analyses showed that the texture of the sediments was sand according to the USDA classification (Table 1).

The the impact of S a on k a , the results from four sediment samples in Dataset I (numbered as #1, #2, #3, and #4) are plotted in Figure 3. As expected, k a increased with increasing S a or with decreasing moisture levels; however, the k a (S a ) relationship also differed among different samples. More specifically, k a showed similar ranges for #1 (13.7 to 41.5 µm 2 ) and #2 (13.2 to 44.5 µm 2 ) with S a ranging from 0.40 to 0.81 and 0.56 to 0.86, respectively (Figures 3a and3b).

In comparison, k a for #3 (11.3 to 142.3 µm 2 ) and #4 (4.5 to 111.5 µm 2 ) varied over an order of magnitude with S a increasing from 0.46 to 0.90 and 0.42 to 0.89, respectively (Figures 3c and3d). Clearly, a single curve of k a (S a ) is not sufficient to describe spatially distributed k a , even at plot scales (e.g., a reach scale in this study). This presents a great challenge for simulating subsurface gas transport with spatially varied k a and associated model parameters. The same challenge is also faced by soil hydrologists to simulate water movement in vadose zone, where soil hydraulic properties may vary substantially across landscapes. One of the solutions to tackle this problem in soil hydrology relies on the use of scaling factors within the framework of the similar media concept [START_REF] Kabat | The scaling characteristics of soil parameters: from plot scale heterogeneity to subgrid parameterization[END_REF][START_REF] Salvucci | Limiting relations between soil moisture and soil texture with implications for measured, modeled and remotely sensed estimates[END_REF][START_REF] Oliveira | Simulation of solute transport in a heterogeneous vadose zone describing the hydraulic properties using a multistep stochastic approach[END_REF].

To have a thorough view of the effectiveness of using the similar media concept for scaling k a , the results of unscaled and scaled k a from Dataset I are plotted in Figure 4 with the reference curves of k a,m for all of the air permeability models. For the purpose of clarity, the results of scaled k a for each model are also plotted separately in Figure 4. The scaled k a was calculated using the obtained scaling factors (e.g., Eq. ( 5)). The obtained fitting parameters and RMSE values for each model are given in Table 2. Note that the results of two sediment samples were removed from the analyses, as they exhibited irregular patterns between k a and S a (e.g., k a may decrease with increasing S a ). Those irregular patterns were suspected to be the consequence of the collapse of root holes in the samples during the drying process. As expected, a positive trend existed between unscaled k a and S a . When S a was lower than 0.7, unscaled k a largely remained below 60 µm 2 , while unscaled k a varied about 8 times from approximately 20 to 160 µm 2 when S a reached above 0.8. Figure 5a further shows that a single curve of k a (S a ) is not sufficient to delineate the distribution of the unscaled k a as attested by the k a,m curves. By comparison, with the aid of scaling factors, the scaled k a coalesced along the k a,m curves with greatly reduced scatter for all of the four air permeability models. Note that scaled k a was systematically more scattered at the lower range of S a for all of the models. It was due to the fact that more weights were given to higher values of k a in Eq. ( 10) during the optimization process. If k a was replaced by log k a in Eq. ( 10), the scatter of scaled k a at the lower range of S a could be reduced; however, the overall performance of the models deteriorated. Compared to the other three k a,m curves, the one for the HT model slightly deviated; whereas, the resulting scaled k a was comparable for all of the four models (Figure 4b). Among all of the results, the RMSE value of the scaled k a for the MQ model was highest with k a,m overestimating scaled k a at the lower end of S a and underestimating scaled k a at the higher end of S a (Figure 4c). This was because without additional fitting parameters other than k o , the form of the MQ model was the least flexible. By contrast, the performances of the other three k a,m curves improved at both ends of S a with the BC model showing the most satisfactory results (Figures 4d-4f); however, there appeared no significant differences in the performances of those three air permeability models for scaling k a (Table 2), indicating the applicability of those models for delineating the distribution of k a in the studied sediments.

The values of obtained k o varied among the four models (Table 2). The fitted k o in the KA model was highest with the value of 106.0 µm 2 , while the one in the MQ model was only 92.5 µm 2 . As a result, the scaled k a for the MQ model was lowest under the same level of S a ; however, the differences in the scaled k a were minimum among different models. The value of the obtained percolation threshold S a,t in the HT model was 0.168, which was comparable to the 0.1 threshold proposed by [START_REF] Ewing | Dependence of the electrical conductivity on saturation in real porous media[END_REF]. In addition, the pore size distribution index γ (γ =3.067 in this study) in the BC model fell within the range for sandy soils (Ghanbarian-Alavijeh and [START_REF] Ghanbarian-Alavijeh | Comparison of the predictions of universal scaling of the saturation dependence of the air permeability with experiment[END_REF]. In summary, Figure 4 shows that for the studied sediments, the spatially distributed k a could be described by a reference curve of k a,m and a set of scaling factors, which attests the feasibility of using the similar media concept for scaling k a under the influence of moisture. With a known reference curve of k a,m and the distribution of scaling factors, this method would be particularly suitable for modeling subsurface gas transport.

As previously explained, the spatially distributed k a is mostly controlled by moisture conditions and soil properties, both of which may vary at different spatiotemporal scales.

Essentially, the impacts of moisture levels and S a on k a can be determined through the reference curve of k a,m ; whereas, the spatial variations of soil properties are embedded in the distribution of scaling factors, as scaling factors only depend on soil intrinsic properties (i.e., the microscopic characteristic length). The probability plots of the scaling factors (both α a,i and lnα a,i ) are given in Figure 5 for the HT model. The resulting probability plots were similar for all of the four models, so only the ones for the HT model are presented here. The respective mean and standard deviation was 1.03 and 0.23 for α a,i , and 0 and 0.23 for lnα a,i . The values of the standard deviations for α a,i and lnα a,i were smaller than previously reported values derived from soil water retention and hydraulic conductivity data [START_REF] Warrick | Scaling field measured soil hydraulic properties using a similar media concept[END_REF][START_REF] Hopmans | A comparison of various methods to scale soil hydraulic properties[END_REF][START_REF] Hendrayanto | Scaling hydraulic properties of forest soils[END_REF], probably due to the smaller number of sediment samples used in this study. Although those previous studies showed a log normal distribution of scaling factors, the Kolmogorov-Smirnov Test did not reject the null hypothesis that α a,i and lnα a,i were normally distributed at the significance level of 0.01. Thus, future studies are needed to include a larger number of samples, and to compare the distributions of scaling factors drawn separately from soil water retention and hydraulic conductivity data and from air permeability data.

To test the applicability of the derived reference curves of k a,m , the values of k a on Day 17

were calculated using the obtained scaling factors and measured S a based on Eq. ( 5). As hypothesized by [START_REF] Iversen | Field application of a portable air permeameter to characterize spatial variability in air and water permeability[END_REF], when soil moisture contents reach field capacity, the air flow takes place in the majority of soil pores; thus, k a measured at or near field capacity can be used as a proxy of soil water permeability k w to predict saturated hydraulic conductivity. On Day 17, the average value of soil moisture contents for all of the sediment samples was 0.04, which was quite comparable to the observed field capacity of sandy soils (soil moisture data can be accessed from the High Plains Regional Climate Center at http://www.hprcc.unl.edu/awdn/soilm/). Based on the above reasoning, the calculated k a on Day 17 is plotted against k w derived from measured saturated hydraulic conductivities in Figure 6.

Again, only the calculated k a for the HT model is shown here as a demonstration. It can be seen from Figure 6 that the calculated k a was very close to k w . On average, the calculated k a only underestimated k w by 11%. Figure 6 further corroborates the feasibility of using scaling factors within the similar media framework for describing spatially distributed k a .

Dataset II

The previous section has demonstrated the feasibility of using the similar media concept for scaling k a in river sediments. To assess the viability of this method in agricultural soils, analyses were carried out for Dataset II. A mean value of 0.02 for S a,t in the HT model was obtained by Ghanbarian-Alavijeh and Hunt (2012) by fitting individual curves of k a to soil samples in Dataset II. Thus, to avoid the negative values of (S a -S a,t ) in Eq. ( 7), S a,t was fixed to be 0.02 for Dataset II. In addition, an unusually high value of γ in the BC model was obtained from the optimization, which was caused by fitting the k a,m curve to comparatively low values of k a from one soil sample (the data points with S a between 0 and 0.2 shown in Figure 7a). Moreover, no physical constraints for the fitting parameters were considered in the optimization process and k a (S a ) in Eq. ( 8) became insensitive to γ for large γ values. Thus, instead of using this high value of γ, the class-averaged value of γ for sandy loam was used, which was derived from soil water retention data [START_REF] Rawls | Estimation of soil water properties[END_REF].

The unscaled and scaled k a from Dataset II are presented in Figure 7, and the obtained fitting parameters and RMSE values are given in Table 3. Overall, the unscaled k a from Dataset II exhibited a much wider distribution from 0.5 to 1980.8 µm 2 . Some of the remolded samples (e.g., the one with the maximum k a of 1980.8 µm 2 ) also showed different ranges of k a (Figure 7a), which was probably caused by the repacking processes. As expected, the improvement of the scaled k a from the MQ model was least satisfactory (Figure 7c). By comparison, the performances of the HT and BC models were considerably improved; however, the obtained k a,m tended to overestimate most of the scaled k a (Figures 7d and 7e Using the same remolded samples from Dataset II, Ghanbarian-Alavijeh and [START_REF] Ghanbarian-Alavijeh | Comparison of the predictions of universal scaling of the saturation dependence of the air permeability with experiment[END_REF] showed that the HT model was better than the other three models (e.g., MQ, BC, and KA) to describe the k a (S a ) relationship for individual samples; however, no significant improvement of the HT model was found in this study for scaling k a . One of the possible reasons for the different conclusions about the model performances was due to the differences in the optimization procedures. In Ghanbarian-Alavijeh and [START_REF] Ghanbarian-Alavijeh | Comparison of the predictions of universal scaling of the saturation dependence of the air permeability with experiment[END_REF], the k a (S a ) curves were optimized for individual samples and some of the fitting parameters (e.g., γ in the BC model) were obtained through optimizing soil retention data rather than soil permeability data.

In comparison to the results from Dataset I, the overall performance of the four selected air permeability models for scaling k a was less satisfactory for Dataset II, which can be attributed to two reasons. First, for the soil samples in Dataset II, soil structures might have been altered

during the repacking processes. Secondly, in order to achieve different moisture levels, soil samples were compressed using an oedometer (Tang et al., 2011), which unavoidably changed the microscopic structures of soils. The repacking processes and the compression of soil samples might have led to variations in the microscopic structures among different soil samples and thus the less satisfactory performance of using the similar media concept for scaling k a . Nonetheless, Figure 7 shows that the use of scaling factors was able to narrow down the scatter of k a around the k a,m curves for the HT, BC, and KA models. In terms of the distribution of scaling factors from Dataset II, the respective mean and standard deviation for the HT model was 1.01 and 0.14 for α a,i , and 0 and 0.13 for lnα a,i . Again, the values of the standard deviations for α a,i and lnα a,i were smaller than the previously reported values derived from soil water retention and hydraulic conductivity data [START_REF] Warrick | Scaling field measured soil hydraulic properties using a similar media concept[END_REF][START_REF] Hopmans | A comparison of various methods to scale soil hydraulic properties[END_REF][START_REF] Hendrayanto | Scaling hydraulic properties of forest soils[END_REF].

Conclusions

The use of the scaling factors based on the similar media concept was tested for scaling air permeability k a using two independent datasets in this study. With a reference curve of k a,m and a set of scaling factors, this method was shown to be able to delineate the spatial distribution of k a for the first dataset, which included k a measured under different moisture conditions by drying river sediment samples. For the second dataset that contained k a measured for agricultural soils, however, the use of the similar media concept for scaling k a was less successful. It was most likely due to the alterations in the microscopic structures of soil samples caused by repacking and compression of soil samples. The merit of this method resides in the fact that the spatial variations of moisture conditions and soil properties can be simultaneously included for analyzing the spatial distribution of k a . At any given location, the impact of S a and thus moisture levels on k a can be quantified by the reference curve of k a,m ; whereas, soil properties are embedded in the scaling factors. Since this study was the first attempt to apply the similar media concept for scaling k a , more studies are needed to test this method on soils with different textures and to examine the impacting factors that control the distribution of scaling factors. It would be also useful to compare scaling factors derived independently from air permeability data, and from soil water retention and hydraulic conductivity data, which may further elucidate the physical meanings of those scaling factors. Future studies are also needed to incorporate this method in modeling subsurface gas transport.

List of Figures

  second dataset (denoted as Dataset II) was taken fromTang et al. (2011) and contained nine remolded soil samples from an experimental forest site at Le Breuil located in central France. The reason for choosing Dataset II was that several air permeability models, including the HT model, were successfully tested on those soil samples(Tang et al., 2011;[START_REF] Ghanbarian-Alavijeh | Comparison of the predictions of universal scaling of the saturation dependence of the air permeability with experiment[END_REF]. However, unlike the air permeameter used in Dataset I, a specially designed oedometer was adopted byTang et al. (2011) to measure k a (S a ) by vertically compressing soil samples, which led to varying moisture contents under different compressional stresses. Except for one sample with 6 data points, there were 8 pairs of k a and S a for each sample. The texture of the soils was sandy loam and their physical properties are reported It has been well known that besides soil properties, k a is also affected by moisture conditions and thus S a . With the drying process, the mean value of S a for Dataset I increased from 0.49 on Day 2 to 0.87 on Day 17. Consequently, the mean value of k a was almost quadrupled from 21.0 µm 2 on Day 2 to 83.6 µm 2 on Day 17, as more smaller pores opened up for transmitting air when the sediments became drier. The results from Dataset I again illustrate the importance of including moisture conditions in analyzing the spatial distribution of k a . To further demonstrate

  ), owing to the fact that the fitting parameter (or shape factor) S a,t in the HT model and γ in the BC model were fixed during the optimizations. Without the need to fix the fitting parameter η in the KA model, the k a,m curve seemed to more reasonably delineate the distribution of the scaled k a over the entire range of S a (Figures 7f), although it should be noted that a very high value of k o was obtained for the KA model in order to match the data points from the sample with the maximum k a of 1980.8 µm 2 .

Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6

 123456 Figure 1 The sampling location for Dataset I near Kearney, NE

Figure 7 Figure 2 Figure 3 Figure 4 Figure 6

 72346 Figure7Unscaled and scaled air permeability (k a ) with the reference curves of k a,m for the selected air permeability models from Dataset II: (a) unscaled k a , (b) scaled k a for all of the four air permeability models, (c) scaled k a for theMillington and Quirk (1960) model, (d) scaled k a for theHunt (2005) model, (e) scaled k a for the[START_REF] Brooks | Hydraulic properties of porous media[END_REF] model, and (f) scaled k a for theKawamoto et al. (2006) model

Figure 7

 7 Figure 7 Unscaled and scaled air permeability (k a ) with the reference curves of k a,m for the selected air permeability models from Dataset II: (a) unscaled k a , (b) scaled k a for all of the four air permeability models, (c) scaled k a for the Millington and Quirk (1960) model, (d) scaled k a for the Hunt (2005) model, (e) scaled k a for the Brooks and Corey (1964) model, and (f) scaled k a for the Kawamoto et al. (2006) model

Acknowledgements

Iain Nicolson Audubon Center at the Rowe Sanctuary in Nebraska allowed the authors to access the Platte River channel for collecting streambed cores and conducting permeameter tests. The authors would like to thank P. Hanson of University of Nebraska-Lincoln for providing assistance in the lab analyses. The research was partially supported by National Basic Research Program of China (2014CB744700).

List of Tables

 Table 1 Summaries [START_REF] Rawls | Estimation of soil water properties[END_REF]


Kenneth L Kuttler

Yves Renard

Meir Shillor

K L Kuttlera

Y Renardb

M Shillorc

Models and simulations of dynamic frictional contact of a thermoelastic beam

We investigate a mathematical model for the dynamic thermomechanical behavior of a viscoelastic beam that is in frictional contact with a rigid moving surface. Friction is modeled by a version of Coulomb's law with slip dependent coefficient of friction, taking into account the frictional heat generation. We prove the existence and uniqueness of the weak solution, describe an algorithm for the numerical solutions and present results of numerical simulations, including the frequency distribution of the noise generated by the stick/slip motion. We also show that when the surface moves too fast there are no steady solutions and therefore the system is thermally unstable.

Introduction

We investigate the dynamic behavior of a one-dimensional thermomechanical system which consists of a beam that is in frictional contact with a rigid obstacle. Our interest lies in the well-posedness of the model and in the stick-slip behavior of the contacting end when the coefficient of friction depends on the relative slip velocity. There is considerable interest in dynamic contact with friction, both mathematically and in applications, especially, since the noise and vibration characteristics of engineering systems are recognized as possible obstacles in product marketing.

We consider a 'simple' setting, so that we can concentrate on the behavior that the solutions exhibit, without having to deal with technical mathematical difficulties. Indeed, the purpose is to gain insight into the possible types of behavior of systems in frictional contact, and use this insight in two-and three-dimensional problems which are much more complicated, both mathematically and numerically.

Often, in engineering literature, the coefficient of friction is assumed to jump from a static value to a smaller dynamic value (see e.g. [START_REF] Rabinowicz | The intrinsic variables affecting the stick-slip process[END_REF] and the survey [ 11] and references therein). Models with a variable coefficient of friction are recent in the mathematical literature (see e.g. [1,5,[START_REF] Kuttler | Set-valued pseudomonotone maps and degenerate evolution inclusions[END_REF][START_REF] Renard | Singular perturbation approach to an elastic dry friction problem with a non-monotone coefficient[END_REF] and references therein). Here, we assume that the coefficient of friction is slip dependent, but the dependence is continuous.

The isothermal problem of frictional contact between a beam and a rigid foundation can be found in [2] where a number of different contact conditions were investigated, and the existence of weak solutions established. In addition, we consider here the thermal aspects of the process. We model the contact with a version of Coulomb's law which was derived recently in [ 1 0] with a slip dependent friction coefficient. The considerable heat generated by frictional contact is taken into account. Indeed, it is found that if too much heat is generated, the problem does not have a steady solution, which indicates thermal instability (see e.g. [3,[START_REF] Yeo | Finite element analysis of thermoelastic contact stability[END_REF] and references therein).

The classical model is presented in Section 2. It consists of the beam equation for a viscoelastic material and the heat equation. The problem is nonlinear and coupled because of the contact and frictional heat generation boundary conditions. Since, in general, there are no classical solutions to problems with contact conditions, we present a weak formulation in Section 3. We state our existence and uniqueness result in Theorem 3.1. The proof is based on a sequence of solutions to a family of linearized and regularized problems. The existence of unique weak solutions to these problems follows from an abstract result of Kuttler [6], which has been generalized and improved in [START_REF] Kuttler | Set-valued pseudomonotone maps and degenerate evolution inclusions[END_REF]. Then, we obtain a priori estimates on these solutions and pass to the limit when the regularization parameter tends to zero. In Section 4 we describe a Crank-Nicholson type of algorithm for the thermoelastic problem. Then, we present four different results of numerical simulations. These are: (i) constant coefficient of friction I.L and constant foundation velocity v*; (ii) constant I.L and oscillatory v*; (iii) slip dependent I.L and oscillatory v*; and (iv) temperature dependent I.L and oscillatory v*. The last case is not covered by our theory. For each of the simulations, we present the Fourier transform of the contacting end's velocity which shows the possible noise characteristics of the solution. In Section 6 we obtain the steady solution of the problem when the frictional heat generation is not too large, since then there is no steady state and the system is thermally unstable. We summarize our findings in Section 7. The results are very interesting and indicate that the problem warrants a more thorough investigation.

The model

In this section, we present the model for the process. The physical setting, depicted in Fig. I, is as follows. A beam is rigidly attached at its left end; its right end is in frictional contact with a moving surface. The area-center of gravity of the beam in its (stress free and isothermal) reference configuration coincides with the interval 0 ~X~ I. For T > 0, we set fly= (0, I) X (0, T); u = u(x, t) and (} = (J(x, t), defined on nT' represent the vertical displacements and the temperature of the beam, respectively.

We employ the Kelvin-Voigt law of viscoelasticity for the shear stress

-er= -er(x, t) = c 2 u",(x, t) + dc 2 u,,,,(x, t),
where c is proportional to the elastic modulus, d (d ;:. 0), is the viscosity coefficient, and subscripts denote partial derivatives. When d = 0 the material is simply elastic. Let f denote the density (per unit length) of vertical applied forces. The equation of motion, in dimensionless form, takes the form u 11 -er, = f, and therefore,

(2.1)

The energy equation can be written, in dimensionless form, as

(} 1 - k(},, = 0 in flT , (2.2)
where k is the (scaled) coefficient of thermal diffusivity. The initial conditions take the form for 0 ~ x ~ I. Here u 0 , v 0 , and 0 0 represent the beam's initial displacement, velocity and temperature, respectively.

u(x, 0) = u 0 (x) , u,(x, 0) = V 0 (x), O(x, 0) = 0 0 (x), (2.3
The beam is rigidly attached to the left wall, and has the wall's temperature, denoted by

0 1 = 0 1 (t), thus u(O, t) = 0, u_(O, t) = 0, 0(0, t) = (}/t), (2.4)
where 0 ~ t ~ T. We allow for time dependence of 0 1 in case it is periodic, although below we simply set it to be the constant ambient temperature.

Our main interest lies in the frictional contact process at the right end. The surface is assumed to move with vertical velocity v* = v*(t), which may be time dependent. A simple argument shows that the contact pressure that results from the horizontal thermal expansion of the beam is given by

p(t) = (a f O(x, t) dx) + , (2.5)
where f+ = max{f, 0}, and a is proportional to the coefficient of thermal expansion. We use ( • ) + to make sure that when the beam shrinks and loses contact because of cooling the pressure vanishes.

Let f..L denote the coefficient of friction, then f..Lh is the friction bound, where h is given by

h = h(t) = p(t)(l -f3p(t)) + .
(2.6)

Here, f3 is a very small positive constant related to the wear and hardness of the surface. This version of the friction bound has been derived in [10]. When f3 = 0 we obtain the usual Columb's law.

We are interested, in particular, in the case when the coefficient of friction depends on the relative slip between the end and the moving surface, thus This is often the case in applications, and for recent results and further references see e.g. [I ,5,7 ,9].

The modified Coulomb's friction condition is

la( I, t)l ~ f..Lh(t), la(!, t)l = f..Lh(t) then u/1, t)-v*(t) = -Aa(l, t), la( I, t)l < f..Lh(t) then u,( I, t) = v *(t) .
(2.7)

When the shear stress a equals ± f..Lh the end is in a slip motion, and the shear direction is opposite to the slip (here A~ 0 is an unknown multiplier); when lal < f..Lh the end sticks to the obstacle and moves with it, which is called stick motion.

The sum of the moments acting on the free end is assumed zero. The thermal boundary condition is chosen as

kO,(l, t) = f..Lh(t)lu,(l, t)-v*(t)l-h,,(O(l, t)-O"(t)).
(2.8)

The first term on the right-hand side represents the frictional heat generated during slip motion and is proportional to the pressure and the slip; the second term represents the thermal interaction between the end and the surface which has temperature 0" = o,,(t).

We summarize the classical problem as: Find a pair {u, 8} such that

u, + c\u,,,, + du'"'') = f in f2r, 0, -kO,, = 0 in {}T , u(x, 0) = u 0 (x), u,(x, 0) = v 0 (x), 0 ~ x ~ I , O(x, 0) = 0 0 (x), 0 ~x ~I, u(O, t) = 0, 0(0, t) = o, ' u,(O, t) = 0 u"( I, t) + du,,, (I, t) = 0, kO,( I, t) = f..Lh(t)lu,( I, t)-v *(t)l-h"(O( I, t) -O"(t)) , (2.9)
(2.10)

(2.11)

(2.12)

(2.13)

(2.14)

(2.15)

(2.16)

10"(1, t)l,;; t-th(t),

I£T(I, t)l = t-th(t) then U 1 (i, t)-v*(t) = -A£T(I, t),
10"( I, t)l < t-th(t) then U 1 ( I, t) = v *(t) for 0,;; t,;; T, where h is given by (2.6) and p by (2.5).

(2.17)

Note that the coupling between the mechanical and thermal parts is through the boundary conditions (2.16) and (2.17) which introduce nonlinearities into the problem.

We assume that t-t = t-t(r) is a smooth bounded function that depends on the magnitude of the slip

r = lu 1 -v*l• It satisfies t-t(O) = f-to > 0
, is decreasing and tends asymptotically to t-t"' > 0. For the sake of convenience, we assume that rt-t '(r) is bounded, (2.18) so that v(r) == rt-t(r) is Lipschitz continuous. The case when t-t is a graph with a vertical segment at the origin is more complicated and will be considered in the sequel.

The friction condition precludes, in general, the existence of smooth solutions to the problem, therefore, we turn to a weak formulation for it.

Existence of the weak solution

We consider a weak abstract formulation of problem (2.9)-(2.17) and prove that it has a unique weak solution. Without loss of generality we assume 0 1 = 0. We define ) and ;JC= L \0, T; H), and let'}': V~~ denote the trace on the right boundary of [0, I] given by yu = u(i ). Next, we define the operator A :V~ V' by

V= {u E H\0, I) : u(O) = u,(O) = 0}, H = L \0, I), 1/= L \0, T; V
(Au, v) = { 1 u,,v" d.x.

Jo

Moreover, let E={O EH 1 (0, I): 0(0) =0}, 'if:=L 2 (0, T;E), and define the operators B, N(t,•): E~E' by

(BO, TJ) = {I 0, TJ, d.x, Jo (N(t, 0), TJ) = -h"(yO-O"(t))'YTJ,
where, with a slight abuse of notation '}' is also the trace of E onto the right endpoint, i.e. yO = 0( I). For 0 E ;Jt and v E /1 we define g(O, yv) E 'if:'= L \0, T; £') by

(g(O, '}'V)(t), TJ) = t-t(lyv(t)-v*(t)l)lyu(t)-v*(t)lh(O)(t)'}'TJ
where for 0 E H. It is straightforward to show that when f3 > 0.

LEMMA 3. 1. There exists a positive constant K such that lh(O)I,;; K, and lh(O)-h(TJ)I,;; KIO-TJif/, lt-t(rl )rl -t-t(r2)r21,;; Kh -r2l• Next, we introduce ~ = ~(t) = £T( I, t), the contact shear, as an additional dependent variable.

The abstract weak form of problem (2.9)-(2.17) is: Find a triplet {u, 0, ~} such that

OE'$, O'E'$', u,u'E'V, u"E'V, ~Ee(O,T), u" + c\Au + dAu') + y*~ = f, 0' + kBO +NO= g(O, yu'), 0(0) = 0 0 EH, u(O) = u 0 E V, u'(O)=v 0 EH,
and for all w E 'V,

We have, under the above assumptions on the data, THEOREM 3.2. There exists a unique solution of problem (3.1)-(3.5).

(3.1)

(3.2) (3.3) (3.4) (3.5)
We prove the existence of a solution by considering a sequence of approximate problems obtained by retarding the arguments of 0 and yu' in g. We define the retardation operator r,, for E > 0, by

() ={w(t-E) if t 3 E, r,w t -0 .f I t< €.
Then, the approximate problem is: Find a triplet {u,, 0,, ~.} such that

0, E 'I;, 0~ E "$', u,, u~ E 'V, u'~ E 'V', ~.EL \0, T), u'~ + c\Au, + dAu:) + y*~. =f, o: + kBO, +NO,= g(r.(}., r, yu:), 0,(0) = 0 0 EH , u,(O) = u 0 E V, u~(O) = V 0 EH,
and for all w E 'V,

(3.6) (3.7) (3.8)
(3.9)

(3.10)

It follows from the results in [2] and Lemma 3.1 that there exists a solution to (3.6)-(3.10) for each small E > 0.

Indeed, first the existence follows for 0,;;;; t,;;;; E, and then the solution is constructed step by step, for step length E, up tot= T.

We proceed to obtain a priori estimates on the approximate solutions. For the sake of simplicity we omit the subscript E, and use C as a positive constant independent of E. We multiply (3.8) by 0, integrate over (0, t) and perform routine estimates to obtain

I L (N(s, O(s)), O(s)) ds I , ; ; ; ; L (C + C([yO(s)i + [Oh(s)[))[yO(s)i ds ,;;;; C + CL [yO(s)[ 2 ds ,;;;; C + CL [[O(s)li~ ds ,
where F is an intermediate space such that y maps F continuously into lffi and E embeds compactly into F, e.g. F = H 31 \0, 1). It follows from the compactness of the embedding that the term is majorized by [[~. ds , (3.11) where, by the definition of E, IIOII~ = [0,1~. Now, using the definition of g and v(r) = J.L(r)r, we find I L (g(r/}, r,(yu')) e) ds I ~ f v(i(yu')(s-E)v*(s)i)h(e(s-E))iye(s) ds which is no larger than (3.28) Now, (g(rA, r,(yu~)), 77) =iT ,u(lr.(yu:)-v*(t)l)lr,(yu~)-v*(t)lh<rA)y7J dt = J.T ,u(lr,(yu')-v*(t)l)lr,(yu')-v*(t)lh(r,O,)y7]dt + g(E), where,

C + C 8 i' [O(s)[~ ds + 8 i' i[O(s)
= L-• v(l( yu')(s)-v*(s + t:)i)h(e(s))!ye(s + t:)l ds ~ C L-• <ICru')(s)l + C)lye(s + t:)l ds ~ C + C J' lru'(s)l 2 ds + C J' lre(s)l 2 ds, 0 0 C + C 8 f
g(E) ~ r iv(ir,(yu~)-v*(t)i)-v(lr,(yu')-v*(t)i)ih(r,O,)y7J dt ~ K Lipv(LT lr,(yu~)r,(yu')l)lr7JI dt,

where Lip v is the Lipschitz constant of v, which converges to zero by (3.23 ).

The first term on the right-hand of (3.29) converges to f T ,u(lyu'-v*(t)l)iyu'-v*(t)lh(O)y7J dt = (g(O, yu'), 7]), ()

(3.29) by Lemma 3.1, the continuity of translations in L 2 and the dominated convergence theorem. Therefore, the limits (3.17)-(3.28) are sufficient to pass to the limit in (3.6)-(3.10) and obtain (3.1)-(3.4). It remains to verify that~ satisfies (3.5). But, this follows from (3.10), the continuity of translations in L 2 , the Lipschitz continuity of h and the pointwise limits together with the dominated convergence theorem. This proves the existence part of the theorem. PROOF OF UNIQUENESS. Suppose {ui, Oi, (), where here and below i = l, 2, are two solutions to the problem. Then, for the mechanical part of the problem, we obtain

l c2 fr 2lu;(t)-u;(t)l~ + 2llu 1 (t)-U 2 (t)il~ + d 0 llu;(s)-u;(s)ll~ ds = -fr ( ~1 (s) -~2(s))( yu; (s) -yu;(s)) ds . () (3.30)
We consider the last term in (3.30) and note that (3.10) implies f r ~2(s)(yu;(s)yu;(s)) ds ~ fr ,u(lyu;(s)-u*(s)i)h( 0

2 )(s) () ()
X (lyu;(s) -v *(s) + yu; (s) -yu ;(s)l -lyu;(s) -v*(s)l} ds, and

-f 1 ~1 (s)( yu; (s)-yu;(s)) ds () = fr ~1 (s)(yu;(s)-yu;(s)) ds ~ fr ,u(lyu;(s)-v*(s)i)h(0 1 )(s) () ()
Therefore, the last term in (3.30) is bounded by

X (lyu;(s) -v *(s)l-lru; (s)-v *(s)l) ds
~ fr (,u(lyu;(s)-v*(s)l)-,u(lyu;(s)-v*(s)l))ih(0 1 )(s)llru;(s)-yu;(s)i ds

~ C f 1 lru;(s)-yu; (s)l 2 ds + C f 1 10 1 (s)-0 2 (s)l~ ds, 0 0
which in turn is bounded, by the compactness of the embedding of V into a space on which the trace map is continuous, by ~ 8 L llu;-u;ll~ ds +CoL iu;-u; I~ ds +CL 10 1 -0 2 1~ ds.

~ C 0 f 1 iu;(s)-u;(s)i~ ds + C f 1 I0 1 (s)-0 2 (s)l~ ds. 0 0 Now, Gronwall's inequality implies iu;(t)-u;Ct)l~ + llu 1 (t)-U 2 (t)il~ + Cd L llu;(s)-u;(s)ll~ ds
~CL 10 1 (s) -0 2 (s)l~ ds .

(3.31)

Letting F = H 31 \0, 1), it follows that the trace map y maps F continuously into IR. Therefore, we have

Since h, is Lipschitz, i(N(t, 0 1 ) -N(t, 0 2 ), 0 1 -0 2 }1 ~ Ciy0 1 -y0 2 i 2 ~ EIIOI -Ozil~: + c.lol -Ozi~ .
Now, we consider g, and find

Lemma 3.1 and the assumption that v* is bounded imply l(g(Op yu;)-g(Oz, yu;), OI-Oz>l ~ C{lru;-yu;llrOI-yOzi

+ iv<iru;-v*I)IIOI -OziHiyOI -yOzi} ~ C{lru;-yu;l 2 + EIIOI-Ozll! + c.lol-Ozi~ +(I+ llu;II~Ml -Ozi~ + EIIOI -Ozil~: + c.lol -Ozl~}.
Therefore, when E is small enough,

IOI(t)-0 2 (t)i~ + f 1 110 1 -0 2 11!-ds ~ C f 1 lru;-yu;l 2 ds + f' (C, + llu;II~M 1 -0 2 1~ ds. 0 0 0
Gronwall's inequality now yields (3.32)

We apply the same method to the thermal boundary condition at the contact point, and this leads to a nonlinear equation.

The friction condition can be divided into two: the stick and the slip conditions, which are of different type. This leads to two nonlinear systems: one for the solution with the stick condition and the other one for the slip solution. When the solution of the continuous model is unique, we expect that at each time step only one of these two holds true. The stick solution is valid if lu(l, t)l < ;..th{t), and the slip solution holds when lu(l, t)l = ;..th(t).

To solve the system, which contains two coupled nonlinear equations corresponding to the friction and thermal conditions, we apply a fixed point method. At each time step we make a first guess for the contact pressure, the friction force and the value of(}~••"+ 1 by taking the values from the previous time step. Then, we solve two linear systems and compute the corresponding values of the contact pressure, friction force and

(} N .u+ I d h . b . h , '
, an t en we Iterate y repeatmg t e process.

Numerical simulations

We now describe a number of numerical simulations of the problem. The purpose of these simulations is to show different types of behavior that the model solutions exhibit. We present only four simulations, a more thorough investigation will be performed in the sequel. Our main interest lies in the stick-slip phenomenon, and the values of the coefficients were chosen accordingly, and not necessarily for a specific material or setting.

For each simulation we show the sliding velocity, the temperature and the Fourier transform of the velocity at

~ Or-~~~~--~~--~--~~--~~ -~ -0.2 ~ -0.4 § -0.6 (J ~ -0.8 §-1.0 -~ -1.2 g a) -1.4 > ~-1.6
] -1.8 the contact point. The last is meant to indicate the possible noise that the slip/stick motion of the contacting end produces. The constants are chosen as c = 10, k = I, h" = I, a= 0.01, {3 = 0 and () 1 = ()" = 400.

The first simulation, presented in Fig. 2, is a simple case where the friction coefficient is constant, J.L = 0.9, and the velocity of the foundation is also constant v* = I. The initial condition is u 0 = 0. In the first period, the beam remains stuck to the surface and the system is under loading. The beam begins to slide when the shear stress threshold, given in the friction law, is reached. Then, the solution is oscillatory about the stationary solution, with a slow convergence toward the stationary solution.

The end oscillations are large, and so are the temperature oscillations, both due to the relatively large coefficient of friction. There is a sharp peak in the frequency domain near 6 Hz, the natural first frequency of this beam.

In the second simulation, presented in Fig. 3, the rigid foundation oscillates with velocity v*(t) = 1 + 2 sin( 2 ; r).

In this particular case, the frequency is equal to 3Hz (T = 113). This leads to a stick-slip phenomenon. Indeed, the end undergoes rapid oscillations with periods of stick, where it moves with the obstacle. The frequency distribution is richer, and more interesting. Higher modes are excited, but the three main peaks are below 10Hz, the smaller ones above it.

The third simulation, shown in Fig. 4, deals with the case of slip dependent coefficient of friction and constant velocity of the rigid foundation v* = l. The slip dependence of J.L is shown in the fourth graph of Fig. decreasing coefficient of this form is often used to model the stick-slip behavior (see e.g. [START_REF] Rabinowicz | The intrinsic variables affecting the stick-slip process[END_REF] or [START_REF] Renard | Singular perturbation approach to an elastic dry friction problem with a non-monotone coefficient[END_REF]), and we observe that this leads to this kind of behavior in the present model, too. Indeed, even a gently decreasing J.L leads to a fairly 'noisy' behavior, which usually is assumed to be caused by a discontinuous coefficient of friction, with different static and dynamic values. The oscillations of the end and of the temperature are noticeable.

In the final example, depicted in Fig. 5, the friction coefficient is chosen to depend on the temperature at the contact point (the dependence is shown in the fourth graph of Fig. 5). The temperature leads to two opposite effects: reduction of the friction force due to the decrease in the friction coefficient and the beam's thermal expansion which increases the contact pressure, and thus the friction force. In the example presented in Fig. 5, there is a stick-slip behavior similar to the one with a slip velocity dependent friction coefficient, but with a lower frequency of the oscillations.

Clearly, even the 'simple' one-dimensional system exhibits varied and interesting types of behavior.

Steady states

We investigate the steady solutions to problem (2.9)-(2.17), with variable J.L. But we use the usual statement of Coulomb's law where h = p and the friction bound is J.LP (i.e. we set f3 = 0). We find that when the heat generation is too large, the system does not have a steady state. We note that in the above setting, only slip can take place in the steady state, which leads to the simplification of conditions (2.17). It follows from the model (2.9)-(2.17), when all the functions are time independent, that the steady problem is: Find a pair { 0 = ii(x), ii = ii(x)} such that This solution is valid only when the frictional heat generation is not too large. Indeed, assume first that v(r) = rJL(r) is monotone increasing, and let v** be the solution of I h" + k = 2 aJL(V**)v**, then, for all v*<v** the solution (5.6) and (5.7) holds. But we note that A~x as v*~v**• If v(r) is monotone decreasing then the steady state exists for v * > u **• Finally, if v(r) is neither, there may be intervals of u* for which the steady solution exists. While intuitively obvious, these results are of considerable applied interest since they determine the long time behavior of the solutions.

The stability properties of the steady states, when existing, are an open question, and so is the rate of convergence of the solution to the steady state.

Conclusions

We presented a 'simple' one-dimensional model for the thermomechanical behavior of a beam that is in frictional contact with a rigid moving obstacle. The problem is very interesting, and the numerical simulations indicate that the system exhibits different types of behavior. Our main interest was in the stick-slip phenomenon, which can be seen for slip dependent coefficient of friction, but the dependence was rather weak, while in many engineering investigations the coefficient of friction is assumed to jump from a static value to a smaller dynamic value. We plan to investigate such problems, but mathematically they are harder and pose more challenge. Moreover, the temperature oscillations indicate that the thermal effects must be included in any consideration of such processes; when the heat generation is too large, there are no steady states and the system is thermally unstable.

This research is a step in our investigation into the vibration and noise characteristics of systems in frictional contact. Clearly, we have only touched the tip of the iceberg.

  Fig. I. The deflected beam in frictional contact.

  Fig. 2. The case with f.1. = const. and v * = const.

Fig. 3 .Fig. 4 .

 34 Fig. 3. The case with periodically oscillating v*.

  ' lu'(s)i~ + 8 f' llu'(s)ll~ ds + 8 f'11e(s)ll~ ds + C 8 f' le(s)l~ ds. Ieo I~+~ J: lle(s)ll~ ds ~ C + C f' le(s)l~ ds + Co f' lu'(s)l~ + 8 f' llu'(s)ll~ ds.

				(3.12)
	()	0	0	0
	(3.11) and (3.12) imply that for 8 small enough,
	~ le(t)l~-~ 0	0	0	(3.13)
	An application of Gronwall's inequality yields
				(3.14)
	0	()	0
				(3.15)
	0		0
	Also, it is straightforward to show that
				(3.16)
	u'~~u"	weakly in /I' ,		(3.18)
	u. ---7 u	strongly in C(O, T; W) ,		(3.19)
	u: ---7u'	strongly in L 2 (0, T; W),		(3.20)
	u:(t) ---7 u'(t) in W a.e. t,		(3.21)
	u. ---7 u	weak* in L x(O, T; V),		(3.22)
	u: ---7u'	weakly in /I,		(3.23)
	(. ---7 (	weakly in L 2 (0, T) ,		(3.24)
	e: ---7 e I	weakly in "6' ,		(3.25)
	e. ---7 e	weakly in "6 ,		(3.26)

le(t)l~ + f' lle(s)ll~ ds ~ C + C f' lu'(s)l~ + C f' llu'(s)ll~ ds.

Now, it follows from similar estimates on u and (3.14) that le(t)l~ + f' l!e(s)ll~ ds + lu'(t)l~ + llu(t)ll~ + f' llu'(s)ll~ ds ~C. Therefore, there exists a subsequence {u,, e., ~J of solutions to (3.6)-(3. 10), for a sequence E ---7 0, such that u:---7u' weak* in L x(O, T; H), (3.17) e. ---7 e strongly in L \0, T; F), (3.27) e. c t) ---7 ec t) in F for a.e. t E [0, T] .

J; ( ~1 -~2 )yw dt ~ 0, for all w E 'V. This proves the theorem. D

Numerical algorithms and simulations

We describe an algorithm for the problem and present numerical simulations of the solutions. The numerical algorithm is based on a Crank-Nicholson finite differences scheme. For the sake of simplicity we consider the elastic case, when d = 0, and when there are no applied forces, f = 0.

We denote by V the column vector (u/x, t), u,Jx, t), O(x, t)), then the evolution equations can be written as a a 2

at V= ax2 (DV), with Let vJ,n = (u;•", u~: 1 , oi•") approximate V at time t" = n !::..t and position xi= j Ax (j = 0, ... , N,), the Crank- Nicholson scheme reads as follows:

An implicit scheme was chosen since such a scheme is stable, and permits to avoid numerical instabilities (see e.g. [4]) due to the friction condition, which may be considered as a stiffness condition (in fact nondifferentiable and multivalued).

The boundary conditions are of two kinds. The conditions u(O, t) = 0, 8(0, t) = 8 1 and uxx(l, t) = 0, can be approximated as o . ., -o u, -'

We use the virtual point method to discretize the thermal boundary condition and the condition u,(O, t) = 0. Indeed, for the condition u,(O, t) = 0, we introduce the virtual point u -I•", and we have

if u is sufficiently smooth. Thus, for an order two approximation we use u -I•" = u 1

•". Applying now the Crank-Nicholson scheme yields u. c, > . . ,"+l -~ul,n+l =uo,., +~ul,., ---------------------- (5.5)

Here, v* is the constant sliding velocity, for convenience chosen in the upward (positive) direction, p = const., and for the sake of simplicity we set c = I, 8 1 = 0 and f = 0. Also, we recall that p = a f 1 : {}(x) dx.

A straightforward calculation shows that the unique solution of the steady state problem is