Marius Cocou Dedicated To Professor

Nicolae Cristescu

A variational analysis of a class of dynamic problems with slip dependent friction

Keywords: Implicit inequalities, set-valued mapping, dynamic problems, pointwise conditions, Coulomb friction Mathematics Subject Classification (2010) : 35Q74, 49J40, 74A55, 74D05, 74H20

This work is concerned with the extension of some recent existence results proved for a class of nonsmooth dynamic frictional contact problem, to the case of a coefficient of friction depending on the slip velocity. Based on existence and approximation results for some general implicit variational inequalities, which are established by using Ky Fan's fixed point theorem, several estimates and compactness arguments, relaxed unilateral conditions with slip dependent friction between two viscoelastic bodies of Kelvin-Voigt type are analyzed.

Introduction

This paper deals with the analysis of a dynamic contact problem with some relaxed unilateral contact conditions, adhesion, and slip dependent pointwise friction, between two Kelvin-Voigt viscoelastic bodies.

The quasistatic unilateral contact problems with local Coulomb friction have been studied in [START_REF] Andersson | Existence results for quasistatic contact problems with Coulomb friction[END_REF][START_REF] Rocca | Existence and approximation of a solution to quasistatic Signorini problem with local friction[END_REF][START_REF] Rocca | Numerical analysis of quasistatic unilateral contact problems with local friction[END_REF] and adhesion laws based on the evolution of intensity of adhesion were investigated in [START_REF] Raous | A consistent model coupling adhesion, friction, and unilateral contact[END_REF][START_REF] Cocou | Existence results for unilateral quasistatic contact problems with friction and adhesion[END_REF]. Also, the normal compliance model, which can be seen as a particular regularization of the Signorini's conditions, has been considered by several authors, see e.g. [START_REF] Kikuchi | Contact Problems in Elasticity : A Study of Variational Inequalities and Finite Element Methods[END_REF][START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF][START_REF] Shillor | Models and Analysis of Quasistatic Contact[END_REF] and references therein.

A recent unified approach including unilateral and bilateral contact with nonlocal friction, or normal compliance conditions, in the quasistatic case and for a nonlinear elastic behavior, has been proposed in [START_REF] Badea | Internal and subspace correction approximations of implicit variational inequalities[END_REF].

In the dynamic case, viscoelastic contact problems with nonlocal friction laws were considered in [START_REF] Jarušek | Dynamic contact problems with given friction for viscoelastic bodies[END_REF][START_REF] Kuttler | Dynamic bilateral contact with discontinuous friction coefficient[END_REF][START_REF] Kuttler | Dynamic contact with Signorini's condition and slip rate depending friction[END_REF][START_REF] Eck | Unilateral Contact Problems -Variational Methods and Existence Theorems[END_REF][START_REF] Cocou | Existence of solutions of a dynamic Signorini's problem with nonlocal friction in viscoelasticity[END_REF][START_REF] Cocou | Analysis of a dynamic unilateral contact problem for a cracked viscoelastic body[END_REF] and the corresponding problems with normal compliance laws have been analyzed in [START_REF] Martins | Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws[END_REF][START_REF] Kikuchi | Contact Problems in Elasticity : A Study of Variational Inequalities and Finite Element Methods[END_REF][START_REF] Kuttler | Dynamic friction contact problems for general normal and friction laws[END_REF][START_REF] Chau | A dynamic frictional contact problem with normal damped response[END_REF][START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities[END_REF].

1 Dynamic frictionless problems with adhesion have been studied by several authors, see, e.g [START_REF] Sofonea | Analysis and Approximation of Contact Problems with Adhesion or Damage[END_REF] and references therein, and dynamic viscoelastic problems coupling unilateral contact, recoverable adhesion and nonlocal friction have been investigated in [START_REF] Cocou | A dynamic unilateral contact problem with adhesion and friction in viscoelasticity[END_REF][START_REF] Cocou | A class of implicit evolution inequalities and applications to dynamic contact problems[END_REF]. Using the Clarke subdifferential, various variational contact problems can be analyzed by using the theory of hemivariational inequalities, which represent a broad generalization of the variational inequalities to locally Lipschitz functions, see [START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities[END_REF][START_REF] Naniewicz | Mathematical Theory of Hemivariational Inequalities and Applications[END_REF][START_REF] Panagiotopoulos | Hemivariational Inequalities: Applications in Mechanics and Engineering[END_REF] and references therein.

Based on new abstract formulations and on Ky Fan's fixed point theorem, a static contact problem with relaxed unilateral conditions and pointwise Coulomb friction was studied in [START_REF] Rabier | Fixed points of multi-valued maps and static Coulomb friction problems[END_REF]. The extension of this interesting approach to an elastic quasistatic contact problem was considered in [START_REF] Cocou | Sur la modélisation mathématique de conditions unilatèrales en mécanique du contact[END_REF] and to a dynamic viscoelastic contact problem with slip independent coefficient of friction was investigated in [START_REF] Cocou | A class of dynamic contact problems with Coulomb friction in viscoelasticity[END_REF].

This paper extends the results presented in [START_REF] Cocou | A class of dynamic contact problems with Coulomb friction in viscoelasticity[END_REF] to the case of a coefficient of friction depending on the slip velocity, which enables to treat more realistic situations.

The paper is organized as follows. In Section 2 the classical formulation of the dynamic contact problem is presented. In Section 3 two variational formulations are given as a two-field problem. In Section 4 a more general evolution implicit variational inequality is considered and some auxiliary results are proved. Section 5 is devoted to the study of a fixed point problem, which is equivalent to the previous variational inequality. Using the Ky Fan's theorem, the existence of a fixed point is proved. In Section 6 this abstract result is used to prove the existence of a variational solution of the dynamic contact problem with slip dependent friction.

The applications presented in this paper concern the contact between two linear viscoelastic bodies but these results can be extended to more general constitutive laws, as, for example, the ones characterizing some elastoviscoplastic materials investigated in [START_REF] Cristescu | Viscoplasticity[END_REF].

Classical formulation

Let Ω α be the reference domains of R d , d = 2 or 3, occupied by two viscoelastic bodies, characterized by a Kelvin-Voigt constitutive law. Suppose that the bodies have Lipschitz boundaries

Γ α = ∂Ω α , α = 1, 2. Let Γ α U , Γ α F and Γ α C be three open disjoint sufficiently smooth parts of Γ α such that Γ α = Γ α U ∪ Γ α F ∪ Γ α C
and, to simplify the estimates, meas(Γ α U ) > 0, α = 1, 2. We shall assume the small deformation hypothesis and we shall use Cartesian coordinate representations.

Let y α (x α , t) denote the position at time t ∈ [0, T ], where 0 < T < +∞, of the material point represented by x α in the reference configuration, and u α (x α , t) := y α (x α , t)x α denote the displacement vector of x α at time t, with the Cartesian coordinates u α = (u α 1 , ..., u α d ) = (ū α , u α d ). Let ε α , with the Cartesian coordinates ε α = (ε ij (u α )), and σ α , with the Cartesian coordinates σ α = (σ α ij ), be the infinitesimal strain tensor and the stress tensor, respectively, corresponding to Ω α , α = 1, 2. The usual summation convention will be used for i, j, k, l = 1, . . . , d.

Assume that the displacement U α = 0 on Γ α U × (0, T ), α = 1, 2, and that the densities of both bodies are equal to 1. Let

f 1 = (f 1 1 , f 2 1 ) denote the given body forces in Ω 1 ∪ Ω 2 and f 2 = (f 1 2 , f 2 2 ) denote the tractions on Γ 1 F ∪ Γ 2 F .
The initial displacements and velocities of the bodies are denoted by

u 0 = (u 1 0 , u 2 0 ), u 1 = (u 1 1 , u 2 1
). Suppose that the solids can be in contact between the potential contact surfaces Γ 1 C and Γ 2 C which are parametrized by two C 1 functions, ϕ 1 , ϕ 2 , defined on an open and bounded subset Ξ of

R d-1 , such that ϕ 1 (ξ)-ϕ 2 (ξ) ≥ 0 ∀ ξ ∈ Ξ and each Γ α C is the graph of ϕ α on Ξ that is Γ α C = { (ξ, ϕ α (ξ)) ∈ R d ; ξ ∈ Ξ}, α = 1, 2.
Define the initial normalized gap between the two contact surfaces by

g 0 (ξ) = ϕ 1 (ξ) -ϕ 2 (ξ) 1 + |∇ϕ 1 (ξ)| 2 ∀ ξ ∈ Ξ
and suppose that this initial gap is sufficiently small. Let n α denote the unit outward normal vector to Γ α , α = 1, 2. We shall use the following notations for the normal and tangential components of a displacement field v α , α = 1, 2, of the relative displacement corresponding to v := (v 1 , v 2 ) and of the stress vector σ α n α on Γ α C :

v α (ξ, t) := v α (ξ, ϕ α (ξ), t), v α N (ξ, t) := v α (ξ, t) • n α (ξ), v N (ξ, t) := v 1 N (ξ, t) + v 2 N (ξ, t), [v N ](ξ, t) := v N (ξ, t) -g 0 (ξ), v α T (ξ, t) := v α (ξ, t) -v α N (ξ, t)n α (ξ), v T (ξ, t) := v 1 T (ξ, t) -v 2 T (ξ, t), σ α N (ξ, t) := (σ α (ξ, t)n α (ξ)) • n α (ξ), σ α T (ξ, t) = σ α (ξ, t)n α (ξ) -σ α N (ξ, t)n α (ξ),
for all ξ ∈ Ξ and for all t

∈ [0, T ]. Let g := -[u N ] = g 0 -u 1 N -u 2
N be the gap corresponding to the solution u := (u 1 , u 2 ). Using a similar method as the one presented in [START_REF] Boieri | Existence, uniqueness, and regularity results for the two-body contact problem[END_REF] (see also [START_REF] Cocou | Analysis of a dynamic unilateral contact problem for a cracked viscoelastic body[END_REF], [START_REF] Cocou | A class of implicit evolution inequalities and applications to dynamic contact problems[END_REF]) we obtain the following unilateral contact condition at time t in the set Ξ: [u N ] (ξ, t) = -g(ξ, t) ≤ 0 ∀ ξ ∈ Ξ. Let A α , B α denote two fourth-order tensors, the elasticity tensor and the viscosity tensor corresponding to Ω α , with the components A α = (A α ijkl ) and B α = (B α ijkl ), respectively. Assume that these components satisfy the following classical symmetry and ellipticity conditions:

C α ijkl = C α jikl = C α klij ∈ L ∞ (Ω α ), ∀ i, j, k, l = 1, . . . , d, ∃ α C α > 0 such that C α ijkl τ ij τ kl ≥ α C α τ ij τ ij ∀ τ = (τ ij ) verifying τ ij = τ ji , ∀ i, j = 1, . . . , d, where C α ijkl = A α ijkl , C α = A α or C α ijkl = B α ijkl , C α = B α ∀ i, j, k, l = 1, . . . , d, α = 1, 2.
Let µ = µ(ξ, uT ) be the slip rate dependent coefficient of friction and assume that µ : Ξ × R d → R + is a bounded function such that for a.e. ξ ∈ Ξ µ(ξ, •) is Lipschitz continuous with the Lipschitz constant, denoted by C µ , independent of ξ, and for every v ∈ R d µ(•, v) is measurable.

Let κ, κ : R → R be two mappings with κ lower semicontinuous and κ upper semicontinuous, satisfying the following conditions:

κ(s) ≤ κ(s) and 0 / ∈ (κ(s), κ(s)) ∀ s ∈ R, (2.1)
∃ r 0 ≥ 0 such that max(|κ(s)|, |κ(s)|) ≤ r 0 ∀ s ∈ R. (2.2)
We consider the following dynamic viscoelastic contact problem.

Problem P c : Find u = (u 1 , u 2 ) such that u(0) = u 0 , u(0) = u 1 and, for all t ∈ (0, T ), üα -div σ α (u α , uα ) = f α 1 in Ω α , (2.3)
σ α (u α , uα ) = A α ε(u α ) + B α ε( uα ) in Ω α , (2.4
)

u α = 0 on Γ α U , σ α n α = f α 2 on Γ α F , α = 1, 2, (2.5)
σ 1 n 1 + σ 2 n 2 = 0 in Ξ, (2.6)
κ([u N ]) ≤ σ N ≤ κ([u N ]) in Ξ, (2.7) |σ T | ≤ µ( uT ) |σ N | in Ξ and
(2.8)

|σ T | < µ( uT ) |σ N | ⇒ uT = 0, |σ T | = µ( uT ) |σ N | ⇒ ∃ ϑ ≥ 0, uT = -ϑσ T ,
where σ α = σ α (u α , uα ), α = 1, 2, σ N := σ 1 N and σ T := σ 1 T . Some contact and friction conditions, corresponding to particular κ and κ, with a general coefficient of friction, are presented in the following examples. Example 1. (Adhesion and friction conditions) Let s 0 ≥ 0, M ≥ 0 be constants, k : R → R be a continuous function such that k ≥ 0 with k(0) = 0 and define

κ(s) =    0 if s ≤ -s 0 , k(s) if -s 0 < s < 0, -M if s ≥ 0, κ(s) =    0 if s < -s 0 , k(s) if -s 0 ≤ s ≤ 0, -M if s > 0.

Example 2. (Friction condition)

In Example 1 we set k = s 0 = 0 and define

κ M (s) = 0 if s < 0, -M if s ≥ 0, κ M (s) = 0 if s ≤ 0, -M if s > 0.
The classical Signorini's conditions correspond formally to M = +∞.

Example 3. (General normal compliance conditions)

Various normal compliance conditions, friction and adhesion laws can be obtained from the previous general formulation if one considers κ = κ = κ, where κ : R → R is some bounded Lipschitz continuous function with κ(0) = 0, so that σ N is given by the relation σ N = κ([u N ]), see e.g. [START_REF] Cocou | A class of implicit evolution inequalities and applications to dynamic contact problems[END_REF], where the intensity of adhesion was also considered.

Variational formulations

We shall consider two different variational formulations of problem P c . We adopt the following notations:

H s (Ω α ) := H s (Ω α ; R d ), α = 1, 2, H s := H s (Ω 1 ) × H s (Ω 2 ), v, w -s,s = v 1 , w 1 H -s (Ω 1 )×H s (Ω 1 ) + v 2 , w 2
H -s (Ω 2 )×H s (Ω 2 ) ∀ v = (v 1 , v 2 ) ∈ H -s , ∀ w = (w 1 , w 2 ) ∈ H s , ∀ s ∈ R.
Define the Hilbert spaces (H, |.|) with the associated inner product denoted by (. , .), (V , . ) with the associated inner product (of H 1 ) denoted by . , . , and the closed convex cones L 2 + (Ξ), L 2 + (Ξ × (0, T )) as follows:

H := H 0 = L 2 (Ω 1 ; R d ) × L 2 (Ω 2 ; R d ), V := V 1 × V 2 ,
where

V α = {v α ∈ H 1 (Ω α ); v α = 0 a.e. on Γ α U }, α = 1, 2, L 2 + (Ξ) := {δ ∈ L 2 (Ξ); δ ≥ 0 a.e. in Ξ}, L 2 + (Ξ × (0, T )) := {η ∈ L 2 (0, T ; L 2 (Ξ)); η ≥ 0 a.e. in Ξ × (0, T )}.
Let a, b be two bilinear, continuous and symmetric mappings defined on

H 1 × H 1 → R by a(v, w) = a 1 (v 1 , w 1 ) + a 2 (v 2 , w 2 ), b(v, w) = b 1 (v 1 , w 1 ) + b 2 (v 2 , w 2 ) ∀ v = (v 1 , v 2 ), w = (w 1 , w 2 ) ∈ H 1 , where, for α = 1, 2, a α (v α , w α ) = Ω α A α ε(v α )•ε(w α ) dx, b α (v α , w α ) = Ω α B α ε(v α )•ε(w α ) dx. Assume f α 1 ∈ W 1,∞ (0, T ; L 2 (Ω α ; R d )), f α 2 ∈ W 1,∞ (0, T ; L 2 (Γ α F ; R d )), α = 1, 2, u 0 , u 1 ∈ V , g 0 ∈ L 2 + (Ξ)
, and define the following mappings:

J : L 2 (Ξ) × (H 1 ) 2 → R, J(δ, v, w) = Ξ µ(v T ) |δ| |w T | dξ ∀ δ ∈ L 2 (Ξ), ∀ v = (v 1 , v 2 ), w = (w 1 , w 2 ) ∈ H 1 , f ∈ W 1,∞ (0, T ; H 1 ), f , v = α=1,2 Ω α f α 1 • v α dx + α=1,2 Γ α F f α 2 • v α ds ∀ v = (v 1 , v 2 ) ∈ H 1 , ∀ t ∈ [0, T ].
Assume the following compatibility conditions:

[u 0N ] ≤ 0, κ([u 0N ]) = 0 a.e. in Ξ and ∃ p 0 ∈ H such that (p 0 , w) + a(u 0 , w) + b(u 1 , w) = f (0), w ∀ w ∈ V . (3.1)
For every ζ ∈ L 2 (0, T ; L 2 (Ξ)) = L 2 (Ξ × (0, T )), define the following sets:

Λ(ζ) = {η ∈ L 2 (0, T ; L 2 (Ξ)); κ • ζ ≤ η ≤ κ • ζ a.e. in Ξ × (0, T ) }, Λ + (ζ) = {η ∈ L 2 + (Ξ × (0, T )); κ + • ζ ≤ η ≤ κ + • ζ a.e. in Ξ × (0, T ) }, Λ -(ζ) = {η ∈ L 2 + (Ξ × (0, T )); κ -• ζ ≤ η ≤ κ -• ζ a.e. in Ξ × (0, T ) },
where, for each r ∈ R, r + := max(0, r) and r -:= max(0, -r) denote the positive and negative parts, respectively.

For each ζ ∈ L 2 (0, T ; L 2 (Ξ)) the sets Λ(ζ), Λ + (ζ) and Λ -(ζ) are clearly closed, convex and nonempty, because the bounding functions belong to the respective set. Since meas(Ξ) < ∞ and κ, κ satisfy (2.2), it follows that for all ζ ∈ L 2 (0, T ; L 2 (Ξ)) these three sets are also bounded in norm in L ∞ (Ξ × (0, T )) by r 0 , and in

L 2 (0, T ; L 2 (Ξ)) by r 1 = r 0 T 1/2 meas(Ξ) 1/2 .
A first variational formulation of the problem P c is the following.

Problem P 1 v : Find u ∈ C 1 ([0, T ]; H -ι )∩W 1,2 (0, T ; V ), λ ∈ L 2 (0, T ; L 2 (Ξ)) such that u(0) = u 0 , u(0) = u 1 , λ ∈ Λ([u N ]) and u(T ), v(T ) -u(T ) -ι, ι -(u 1 , v(0) -u 0 ) - T 0 ( u, v -u) dt + T 0 a(u, v -u) + b( u, v -u) -(λ, v N -u N ) L 2 (Ξ) dt (3.2) + T 0 {J(λ, u, v + k u -u) -J(λ, u, k u)} dt ≥ T 0 f , v -u dt ∀ v ∈ L ∞ (0, T ; V ) ∩ W 1,2 (0, T ; H), where 1 > ι > 1 2 , k > 0.
The formal equivalence between the variational problem P 1 v and the classical problem (2.3)-(2.8) can be easily proved by using Green's formula and an integration by parts, where the Lagrange multiplier λ satisfies the relation

λ = σ N . Let φ : (L 2 + (Ξ)) 2 × (V ) 2 → R be defined by φ(δ 1 , δ 2 , v, w) = -(δ 1 -δ 2 , w N ) L 2 (Ξ) + Ξ µ(v T ) (δ 1 + δ 2 ) |w T | dξ ∀ (δ 1 , δ 2 ) ∈ (L 2 + (Ξ)) 2 , ∀ v = (v 1 , v 2 ), w = (w 1 , w 2 ) ∈ V . (3.3) Since η ∈ Λ(ζ) if and only if (η + , η -) ∈ Λ + (ζ) × Λ -(ζ)
, it follows that the variational problem P 1 v is equivalent with the following problem.

Problem P 2 v : Find u ∈ C 1 ([0, T ]; H -ι )∩W 1,2 (0, T ; V ), λ ∈ L 2 (0, T ; L 2 (Ξ)) such that u(0) = u 0 , u(0) = u 1 , (λ + , λ -) ∈ Λ + ([u N ]) × Λ -([u N ]) and u(T ), v(T ) -u(T ) -ι, ι -(u 1 , v(0) -u 0 ) + T 0 {-( u, v -u) + a(u, v -u) + b( u, v -u)} dt + T 0 {φ(λ + , λ -, u, v + k u -u) -φ(λ + , λ -, u, k u)} dt ≥ T 0 f , v -u dt ∀ v ∈ L ∞ (0, T ; V ) ∩ W 1,2 (0, T ; H). (3.4)
The existence of variational solutions of the problem P c will follow from some general existence results that will be proved in the next sections.

Existence results for some variational inequalities

Let U 0 , (V 0 , . , . , . ), (U, . U ) and (H 0 , |.|, (. , .)) be four Hilbert spaces such that U 0 is a closed linear subspace of V 0 dense in H 0 , V 0 ⊂ U ⊆ H 0 with continuous embeddings and the embedding from V 0 into U is compact. Let B r (Ξ), B r (Ξ T ) denote the closed balls with center 0 and radius r in L ∞ (Ξ), L ∞ (Ξ T ), respectively, where Ξ T := Ξ × (0, T ) and r > 0.

Let a 0 , b 0 : V 0 × V 0 → R be two bilinear and symmetric forms such that

∃ M a , M b > 0 a 0 (u, v) ≤ M a u v , b 0 (u, v) ≤ M b u v , (4.1)
∃ m a , m b > 0 a 0 (v, v) ≥ m a v 2 , b 0 (v, v) ≥ m b v 2 ∀ u, v ∈ V 0 . (4.2) Let l : V 0 → L 2 (Ξ) and φ 0 : [0, T ] × (L 2 + (Ξ)) 2 × (V 0 ) 2 →
R be two mappings satisfying the following conditions:

∃ k 1 > 0 such that ∀ v 1 , v 2 ∈ V 0 , l(v 1 ) -l(v 2 ) L 2 (Ξ) ≤ k 1 v 1 -v 2 U , (4.3)
∀ t ∈ [0, T ], ∀ γ 1 , γ 2 ∈ L 2 + (Ξ), ∀ v, v 1 , v 2 ∈ V 0 , φ 0 (t, γ 1 , γ 2 , v, v 1 + v 2 ) ≤ φ 0 (t, γ 1 , γ 2 , v, v 1 ) + φ 0 (t, γ 1 , γ 2 , v, v 2 ), (4.4) φ 0 (t, γ 1 , γ 2 , v, θv 1 ) = θ φ 0 (t, γ 1 , γ 2 , v, v 1 ), ∀ θ ≥ 0, (4.5) φ 0 (t, γ 1 , γ 2 , v, w) = 0, ∀ w ∈ U 0 , (4.6)
φ 0 (0, 0, 0, 0, v) = 0, (4.7)

∀ r > 0, ∃ k 2 (r) > 0 such that ∀ t 1 , t 2 ∈ [0, T ], ∀ γ 1 , γ 2 , δ 1 , δ 2 ∈ L 2 + (Ξ) ∩ B r (Ξ), ∀ v 1 , v 2 , w 1 , w 2 ∈ V 0 , |φ 0 (t 1 , γ 1 , γ 2 , v 1 , w 1 ) -φ 0 (t 1 , γ 1 , γ 2 , v 1 , w 2 ) +φ 0 (t 2 , δ 1 , δ 2 , v 2 , w 2 ) -φ 0 (t 2 , δ 1 , δ 2 , v 2 , w 1 )| ≤ k 2 (r)(|t 1 -t 2 | + γ 1 -δ 1 L 2 (Ξ) + γ 2 -δ 2 L 2 (Ξ) + v 1 -v 2 U ) w 1 -w 2 U , (4.8) if (γ n 1 , γ n 2 ) ∈ (L 2 + (Ξ T )) 2 for all n ∈ N and (γ n 1 , γ n 2 ) ⇀ (γ 1 , γ 2 ) in (L 2 (0, T ; L 2 (Ξ))) 2 , then T 0 φ 0 (s, γ n 1 , γ n 2 , v, w) ds → T 0 φ 0 (s, γ 1 , γ 2 , v, w) ds ∀ v, w ∈ L 2 (0, T ; V 0 ).
(4.9)

Remark 4.1 i) Since by (4.5) φ 0 (•, •, •, •, 0) = 0, from (4.8), for w 1 = w, w 2 = 0, we have

∀ t 1 , t 2 ∈ [0, T ], ∀ γ 1 , γ 2 , δ 1 , δ 2 ∈ L 2 + (Ξ) ∩ B r (Ξ), ∀ v 1 , v 2 , w ∈ V 0 , |φ 0 (t 1 , γ 1 , γ 2 , v 1 , w) -φ 0 (t 2 , δ 1 , δ 2 , v 2 , w)| ≤ k 2 (r)(|t 1 -t 2 | + γ 1 -δ 1 L 2 (Ξ) + γ 2 -δ 2 L 2 (Ξ) + v 1 -v 2 U ) w U . (4.10)
ii) From (4.7) and (4.8), for

t 1 = t, t 2 = 0, δ 1 = δ 2 = 0 and v 1 = v, v 2 = 0 we derive ∀ t ∈ [0, T ], ∀ γ 1 , γ 2 ∈ L 2 + (Ξ) ∩ B r (Ξ), ∀ v, w 1 , w 2 ∈ V 0 , |φ 0 (t, γ 1 , γ 2 , v, w 1 ) -φ 0 (t, γ 1 , γ 2 , v, w 2 )| ≤ k 2 (r)(t + γ 1 L 2 (Ξ) + γ 2 L 2 (Ξ) + v U ) w 1 -w 2 U . (4.11) iii) If v n → v, w m → w in L 2 (0, T ; U ), (γ n 1 , γ n 2 ) ∈ (L 2 + (Ξ T ) ∩ B r (Ξ T )) 2 , for all n ∈ N, and (γ n 1 , γ n 2 ) ⇀ (γ 1 , γ 2 ) in (L 2 (0, T ; L 2 (Ξ))) 2 , then lim n,m→∞ T 0 φ 0 (s, γ n 1 , γ n 2 , v n , w m ) ds → T 0 φ 0 (s, γ 1 , γ 2 , v, w) ds, (4.12)
which can be proved by taking into account (4.11) in the following relations:

| T 0 {φ 0 (s, γ n 1 , γ n 2 , v n , w m ) -φ 0 (s, γ 1 , γ 2 , v, w)} ds| ≤ T 0 |φ 0 (s, γ n 1 , γ n 2 , v n , w m ) -φ 0 (s, γ n 1 , γ n 2 , v n , w)| ds + T 0 |{φ 0 (s, γ n 1 , γ n 2 , v n , w) -φ 0 (s, γ n 1 , γ n 2 , v, w)}| ds +| T 0 {φ 0 (s, γ n 1 , γ n 2 , v, w) -φ 0 (s, γ 1 , γ 2 , v, w)} ds| ≤ T 0 k 2 (r)( γ n 1 L 2 (Ξ) + γ n 2 L 2 (Ξ) + v n U ) w m -w U ds + T 0 k 2 (r)( γ n 1 L 2 (Ξ) + γ n 2 L 2 (Ξ) + v n -v U ) w U ds +| T 0 {φ 0 (s, γ n 1 , γ n 2 , v, w) -φ 0 (s, γ 1 , γ 2 , v, w)} ds|,
and passing to limits by using that (γ n 1,2 ) n are bounded and (4.9).

Assume that f 0 ∈ W 1,∞ (0, T ; V 0 ), u 0 , u 1 ∈ V 0 are given, and that the following compatibility condition holds: κ(l(u 0 )) = 0 and ∃ p 0 ∈ H 0 such that (p 0 , w) + a 0 (u 0 , w)

+ b 0 (u 1 , w) = f 0 (0), w ∀ w ∈ V 0 . (4.13)
We consider the following problem.

Problem Q 1 : Find u ∈ W 0 , λ ∈ L 2 (0, T ; L 2 (Ξ)) such that u(0) = u 0 , u(0) = u 1 , (λ + , λ -) ∈ Λ + (l(u)) × Λ -(l(u)) and u(T ), v(T ) -u(T ) U ′ × U -(u 1 , v(0) -u 0 ) + T 0 {-( u, v -u) + a 0 (u, v -u) + b 0 ( u, v -u)} dt + T 0 {φ 0 (t, λ + , λ -, u, v + k u -u) -φ 0 (t, λ + , λ -, u, k u)} dt ≥ T 0 f 0 , v -u dt ∀ v ∈ L ∞ (0, T ; V 0 ) ∩ W 1,2 (0, T ; H 0 ),
where

W 0 := C 1 ([0, T ]; U ′ ) ∩ W 1,2 (0, T ; V 0 ). The sets Λ + (ζ), Λ - (ζ
) and Λ(ζ) have the following useful properties, proved in [START_REF] Cocou | A class of implicit evolution inequalities and applications to dynamic contact problems[END_REF], see also [START_REF] Rabier | Fixed points of multi-valued maps and static Coulomb friction problems[END_REF].

Lemma 4.1 Let ζ ∈ L 2 (0, T ; L 2 (Ξ)) and (η 1 , η 2 ) ∈ Λ + (ζ) × Λ -(ζ). Then η 1 η 2 = 0 a.e. in Ξ T and there exists η ∈ Λ(ζ) such that η + = η 1 , η -= η 2 a.e. in Ξ T .
Based on the previous lemma, consider the following problem, which has the same solution u as the problem Q 1 , and the solutions λ 1 , λ 2 satisfy the relation

λ = λ 1 -λ 2 , where λ is a solution of Q 1 . Problem Q 2 : Find u ∈ W 0 , λ 1 , λ 2 ∈ L 2 (0, T ; L 2 (Ξ)) such that u(0) = u 0 , u(0) = u 1 , (λ 1 , λ 2 ) ∈ Λ + (l(u)) × Λ -(l(u)) and u(T ), v(T ) -u(T ) U ′ × U -(u 1 , v(0) -u 0 ) + T 0 {-( u, v -u) + a 0 (u, v -u) + b 0 ( u, v -u)} dt + T 0 {φ 0 (t, λ 1 , λ 2 , u, v + k u -u) -φ 0 (t, λ 1 , λ 2 , u, k u)} dt ≥ T 0 f 0 , v -u dt ∀ v ∈ L ∞ (0, T ; V 0 ) ∩ W 1,2 (0, T ; H 0 ). (4.14)
For the convenience of the reader, an existence and uniqueness result proved in [START_REF] Cocou | Analysis of a dynamic unilateral contact problem for a cracked viscoelastic body[END_REF] will be restated here under an adapted and more general form that will enable to study problem Q 2 .

Let β : V 0 → R and φ 1 : [0, T ] × V 3 0 → R be two sequentially weakly continuous mappings such that

β(0) = 0 and φ 1 (t, z, v, w 1 + w 2 ) ≤ φ 1 (t, z, v, w 1 ) + φ 1 (t, z, v, w 2 ), (4.15)
φ 1 (t, z, v, θw) = θ φ 1 (t, z, v, w), (4.16
)

φ 1 (0, 0, 0, w) = 0 ∀ t ∈ [0, T ], ∀ z, v, w, w 1,2 ∈ V 0 , ∀ θ ≥ 0, (4.17)
∃ k 3 > 0 such that ∀ t 1,2 ∈ [0, T ], ∀ u 1,2 , v 1,2 , w ∈ V 0 , |φ 1 (t 1 , u 1 , v 1 , w) -φ 1 (t 2 , u 2 , v 2 , w)| ≤ k 3 (|t 1 -t 2 | + |β(u 1 -u 2 )| + v 1 -v 2 U ) w , (4.18)
∃ k 4 > 0 such that ∀ t 1,2 ∈ [0, T ], ∀ u 1,2 , v 1,2 , w 1,2 ∈ V 0 , |φ 1 (t 1 , u 1 , v 1 , w 1 ) -φ 1 (t 1 , u 1 , v 1 , w 2 ) + φ 1 (t 2 , u 2 , v 2 , w 2 ) -φ 1 (t 2 , u 2 , v 2 , w 1 )| ≤ k 4 ( |t 1 -t 2 | + u 1 -u 2 + v 1 -v 2 U ) w 1 -w 2 . (4.19)
Let L ∈ W 1,∞ (0, T ; V 0 ) and assume the following compatibility condition on the initial data:

∃ p 1 ∈ H 0 such that (p 1 , w)+a 0 (u 0 , w)+b 0 (u 1 , w)+φ 1 (0, u 0 , u 1 , w) = L(0), w ∀ w ∈ V 0 . (4.20)
Consider the following problem. Problem Q 3 : Find u ∈ W 2,2 (0, T ; H 0 )∩W 1,2 (0, T ; V 0 ) such that u(0) = u 0 , u(0) = u 1 , and for almost all t ∈ (0, T )

(ü, v -u) + a 0 (u, v -u) + b 0 ( u, v -u) +φ 1 (t, u, u, v) -φ 1 (t, u, u, u) ≥ L, v -u ∀ v ∈ V 0 . (4.21)
We have the following existence and uniqueness result. The proof, which will be presented in a forthcoming paper, is based on a similar method to that used to prove the Theorem 3.2 established in [START_REF] Cocou | Analysis of a dynamic unilateral contact problem for a cracked viscoelastic body[END_REF] and on a useful estimate, see [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] or [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF], which, when applied to the spaces V 0 ⊂ U ⊆ H 0 , implies the following result: for every ǫ > 0 there exists C ǫ > 0 such that

u U ≤ ǫ u + C ǫ |u| ∀ u ∈ V 0 . ( 4
) ∈ (W 1,∞ (0, T ; L 2 (Ξ))) 2 ∩ (L 2 + (Ξ T )) 2 ∩ (B r (Ξ T )) 2 with γ 1 (0) = γ 2 (0) = 0, there exists a unique solution u = u (γ 1 ,γ 2 )
of the following evolution variational inequality: find u ∈ W 2,2 (0, T ; H 0 ) ∩ W 1,2 (0, T ; V 0 ) such that u(0) = u 0 , u(0) = u 1 , and for almost all t ∈ (0, T )

(ü, v -u) + a 0 (u, v -u) + b 0 ( u, v -u) +φ 0 (t, γ 1 , γ 2 , u, v) -φ 0 (t, γ 1 , γ 2 , u, u) ≥ f 0 , v -u ∀ v ∈ V 0 . (4.23)
Proof. We apply Theorem 4.1 to β = 0, L = f 0 and

φ 1 (t, z, v, w) = φ 0 (t, γ 1 (t), γ 2 (t), v, w) ∀ t ∈ [0, T ], ∀ z, v, w ∈ V 0 .
Since φ 0 satisfies (4.4), (4.5) and (4.7), one can easily verify the properties (4.15)-(4.17). Also, (4.13) and (4.17) imply the condition (4.20).

Using (4.8), for some k 2 = k 2 (r) > 0 we have

∀ t 1,2 ∈ [0, T ], ∀ u 1,2 , v 1,2 , w 1,2 ∈ V 0 , |φ 1 (t 1 , u 1 , v 1 , w 1 ) -φ 1 (t 1 , u 1 , v 1 , w 2 ) + φ 1 (t 2 , u 2 , v 2 , w 2 ) -φ 1 (t 2 , u 2 , v 2 , w 1 )| = |φ 0 (t 1 , γ 1 (t 1 ), γ 2 (t 1 ), v 1 , w 1 ) -φ 0 (t 1 , γ 1 (t 1 ), γ 2 (t 1 ), v 1 , w 2 )| +φ 0 (t 2 , γ 1 (t 2 ), γ 2 (t 2 ), v 2 , w 2 ) -φ 0 (t 2 , γ 1 (t 2 ), γ 2 (t 2 ), v 2 , w 1 )| ≤ k 2 (|t 1 -t 2 | + γ 1 (t 1 ) -γ 1 (t 2 ) L 2 (Ξ) + γ 2 (t 1 ) -γ 2 (t 2 ) L 2 (Ξ) + v 1 -v 2 U ) w 1 -w 2 U ≤ k 2 ((1 + C γ 1 + C γ 2 )|t 1 -t 2 | + v 1 -v 2 U ) w 1 -w 2 U , ≤ k 5 (|t 1 -t 2 | + v 1 -v 2 U ) w 1 -w 2 U ,
where C γ 1 , C γ 2 denote the Lipschitz constants of γ 1 , γ 2 , respectively, and

k 5 = k 2 (1 + C γ 1 + C γ 2 ). Thus, |φ 1 (t 1 , u 1 , v 1 , w 1 ) -φ 1 (t 1 , u 1 , v 1 , w 2 ) + φ 1 (t 2 , u 2 , v 2 , w 2 ) -φ 1 (t 2 , u 2 , v 2 , w 1 )| ≤ k 5 (|t 1 -t 2 | + v 1 -v 2 U ) w 1 -w 2 U ∀ t 1,2 ∈ [0, T ], ∀ u 1,2 , v 1,2 , w 1,2 ∈ V 0 , (4.24)
and, since by the continuous embedding V 0 ⊂ U there exists

C U > 0 such that w U ≤ C U w ∀ w ∈ V 0 , it follows that φ 1 satisfies (4.19) with k 4 = k 5 C U .
Taking in (4.24) w 1 = w, w 2 = 0, by (4.16) with θ = 0, we obtain

|φ 1 (t 1 , u 1 , v 1 , w) -φ 1 (t 2 , u 2 , v 2 , w)| ≤ k 5 (|t 1 -t 2 | + v 1 -v 2 U ) w U ∀ t 1,2 ∈ [0, T ], ∀ u 1,2 , v 1,2 , w ∈ V 0 , (4.25)
and using the continuous embedding V 0 ⊂ U , it follows that φ 1 satisfies (4.18) with k

3 = k 5 C U . inequality: u(0) = u 0 , u(0) = u 1 , u(T ), v(T ) -u(T ) U ′ × U -(u 1 , v(0) -u 0 ) + T 0 {-( u, v -u) + a 0 (u, v -u) + b 0 ( u, v -u)} dt (5.1) + T 0 {φ 0 (t, γ 1 , γ 2 , u, v -u + k u) -φ 0 (t, γ 1 , γ 2 , u, k u)} dt ≥ T 0 f 0 , v -u dt ∀ v ∈ L ∞ (0, T ; V 0 ) ∩ W 1,2 (0, T ; H 0 ).
The proof, based on Lemmas 4.2, 4.3 and on some compactness results established in [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF], is similar to the proof of Theorem 3.2 in [START_REF] Cocou | A class of dynamic contact problems with Coulomb friction in viscoelasticity[END_REF], so that will be not presented here. Now, let Φ :

(L 2 + (Ξ T ) ∩ B r 0 (Ξ T )) 2 → 2 (L 2 + (Ξ T )∩Br 0 (Ξ T )) 2 \ {∅} be the set-valued mapping defined by Φ(γ 1 , γ 2 ) = Λ + (l(u (γ 1 ,γ 2 ) )) × Λ -(l(u (γ 1 ,γ 2 ) )) for all (γ 1 , γ 2 ) ∈ (L 2 + (Ξ T ) ∩ B r 0 (Ξ T )) 2 , (5.2)
where u (γ 1 ,γ 2 ) is the solution of the variational inequality (5.1) which corresponds to (γ 1 , γ 2 ) by the procedure described in Theorem 5.1.

It is clear that if (λ 1 , λ 2 ) is a fixed point of Φ, i.e. (λ 1 , λ 2 ) ∈ Φ(λ 1 , λ 2 ), then (u (λ 1 ,λ 2 ) , λ 1 , λ 2 ) is a solution of the problem Q 2 .
Consider a new problem, which consists in finding a fixed point of the set-valued mapping Φ, called also multivalued function or multifunction, which will provide a solution of problem Q 1 .

The existence of a fixed point of the multifunction Φ will be obtained by using a corollary of the Ky Fan's fixed point theorem [START_REF] Fan | Fixed points and minimax theorems in locally convex topological linear spaces[END_REF], proved in [START_REF] Rabier | Fixed points of multi-valued maps and static Coulomb friction problems[END_REF] in the particular case of a reflexive Banach space. Note that since Y is a reflexive Banach space and D is convex, closed and bounded, we don't need to assume that Y is separable, see [START_REF] Browder | Nonlinear operators and nonlinear equations of eyolutionin Banach spaces[END_REF][START_REF] Rabier | Fixed points of multi-valued maps and static Coulomb friction problems[END_REF].

= | -(γ 1 -γ 2 , w 1N ) L 2 (Ξ) + Ξ µ(v 1T ) (γ 1 + γ 2 ) |w 1T | dξ +(γ 1 -γ 2 , w 2N ) L 2 (Ξ) - Ξ µ(v 1T ) (γ 1 + γ 2 ) |w 2T | dξ -(δ 1 -δ 2 , w 2N ) L 2 (Ξ) + Ξ µ(v 2T ) (δ 1 + δ 2 ) |w 2T | dξ +(δ 1 -δ 2 , w 1N ) L 2 (Ξ) - Ξ µ(v 2T ) (δ 1 + δ 2 ) |w 1T | dξ| ≤ |(γ 1 -γ 2 -δ 1 + δ 2 , w 1N -w 2N ) L 2 (Ξ) | +| Ξ (µ(v 1T ) -µ(v 2T )) (γ 1 + γ 2 ) (|w 1T | -|w 2T |) dξ| +| Ξ µ(v 2T )(γ 1 + γ 2 -δ 1 -δ 2 ) (|w 1T | -|w 2T |) dξ| ≤ ( γ 1 -δ 1 L 2 (Ξ) + γ 2 -δ 2 L 2 (Ξ) ) w 1N -w 2N L 2 (Ξ) + Ξ C µ |v 1T -v 2T | (|γ 1 | + |γ 2 |) |w 1T -w 2T | dξ + Ξ µ(v 2T ) (|γ 1 -δ 1 | + |γ 2 -δ 2 |) |w 1T -w 2T | dξ ≤ ( γ 1 -δ 1 L 2 (Ξ) + γ 2 -δ 2 L 2 (Ξ) ) w 1 -w 2 (L 2 (Ξ)) d +2rC µ Ξ |v 1 -v 2 | |w 1 -w 2 | dξ +M µ Ξ (|γ 1 -δ 1 | + |γ 2 -δ 2 |) |w 1 -w 2 | dξ ≤ ( γ 1 -δ 1 L 2 (Ξ) + γ 2 -δ 2 L 2 (Ξ) ) w 1 -w 2 (L 2 (Ξ)) d +2rC µ v 1 -v 2 (L 2 (Ξ)) d w 1 -w 2 (L 2 (Ξ)) d +M µ ( γ 1 -δ 1 L 2 (Ξ) + |γ 2 -δ 2 | L 2 (Ξ) ) w 1 -w 2 (L 2 (Ξ)) d ≤ (1 + M µ C tr )( γ 1 -δ 1 L 2 (Ξ) + γ 2 -δ 2 L 2 (Ξ) ) w 1 -w 2 H ι +2rC µ C 2 tr v 1 -v 2 H ι w 1 -w 2 H ι ,
and so (4.8) is satisfied with k 2 (r)=max(1 + M µ C tr , 2rC µ C 2 tr ), where M µ is an upper bound of µ.

Note that the same method can provide a unified approach to study more complex dynamic surface interactions, for which the evolution of the intensity of adhesion is governed by a variational inequality or a differential equation, see e.g. [START_REF] Sofonea | Analysis and Approximation of Contact Problems with Adhesion or Damage[END_REF], [START_REF] Shillor | Models and Analysis of Quasistatic Contact[END_REF], [START_REF] Cocou | A dynamic unilateral contact problem with adhesion and friction in viscoelasticity[END_REF], [START_REF] Cocou | A class of implicit evolution inequalities and applications to dynamic contact problems[END_REF].

Theorem 4 . 1

 41 Under the assumptions (4.1), (4.2), (4.15)-(4.20), there exists a unique solution to the problem Q 3 .

Definition 5 . 1

 51 Let Y be a reflexive Banach space, D a weakly closed set in Y , and F : D → 2 Y \ {∅} be a multivalued function. F is called sequentially weakly upper semicontinuous if z n ⇀ z, y n ∈ F (z n ) and y n ⇀ y imply y ∈ F (z).

Proposition 5 . 1 (

 51 [START_REF] Rabier | Fixed points of multi-valued maps and static Coulomb friction problems[END_REF]) Let Y be a reflexive Banach space, D a convex, closed and bounded set in Y , and F : D → 2 D \ {∅} a sequentially weakly upper semicontinuous multivalued function such that F (z) is convex for every z ∈ D. Then F has a fixed point.

Now, taking in (4.24) t 1 = t, t 2 = 0, u 1 = z, v 1 = v, u 2 = v 2 = 0, by (4.17) we have

(4.26)

As the embedding from V 0 into U is compact, from (4.25) and (4.26) it follows that φ 1 is weakly sequentially continuous. By Theorem 4.1 there exists a unique solution u = u (γ 1 ,γ 2 ) of the variational inequality (4.23).

Also, we have the following result, which is similar to Lemma 3.3 in [START_REF] Cocou | A class of dynamic contact problems with Coulomb friction in viscoelasticity[END_REF].

be the corresponding solutions of (4.23). Then there exists a constant

(4.27)

An equivalent fixed point problem

Since D(0, T ; L 2 (Ξ)) is dense in L 2 (0, T ; L 2 (Ξ)), which is classically proved by using the convolution product with suitable mollifiers, it follows that for every γ ∈ L 2 + (Ξ T ) ∩ B r 0 (Ξ T ), there exist r > 0 and a sequence (

Theorem 5.1 Assume that (4.1), (4.2), (4.4)-(4.9) and ( 4

be the solution of (4.23)

) n is strongly convergent in W 0 , its limit, denoted by u = u (γ 1 ,γ 2 ) , is independent of the chosen sequence converging to (γ

Proof. The proof is similar to the proof of Theorem 3.3 in [START_REF] Cocou | A class of dynamic contact problems with Coulomb friction in viscoelasticity[END_REF] but, for the convenience of the reader, we shall present it. By Lemma 4.

We apply Proposition 5.

The set D ⊂ (L 2 (0, T ; L 2 (Ξ))) 2 is clearly convex, closed and bounded. Since for each ζ ∈ L 2 (0, T ; L 2 (Ξ)) the sets Λ + (ζ) and Λ -(ζ) are nonempty, convex, closed, and bounded by r 0 , it follows that Φ(γ 1 , γ 2 ) is a nonempty, convex and closed subset of D for every (γ 1 , γ 2 ) ∈ D.

In order to prove that the multifunction Φ is sequentially weakly upper semicontinuous, let (

) and let us verify that (η 1 , η 2 ) ∈ Φ(γ 1 , γ 2 ). By Theorem 5.1, there exists a sequence (γ

where

is the solution of (4.23) corresponding to (γ n 1 , γn 2 ), u (γ n 1 ,γ n

2 ) is the solution of (5.1) corresponding to (γ n 1 , γ n 2 ) and to the procedure which enables to define Φ(γ n 1 , γ n 2 ). As (γ n 1 , γ n 2 ) ⇀ (γ 1 , γ 2 ) in (L 2 (0, T ; L 2 (Ξ))) 2 , Theorem 5.1 implies also that u (γ n 1 ,γ n 2 ) → u (γ 1 ,γ 2 ) in W 0 , and by (5.3) and the triangle inequality, we obtain

which may be rewritten as

for all n ∈ N, where

), and also under the following equivalent form:

for every measurable subset ω ⊂ Ξ T and for all n ∈ N.

Using (5.4), (4.3), the semi-continuity of κ and κ, the relation (2.2), the convergence property

, and passing to limits according to Fatou's lemma (see also [START_REF] Rabier | Fixed points of multi-valued maps and static Coulomb friction problems[END_REF]), we obtain

for every measurable subset ω ⊂ Ξ T , which implies (η 1 , η 2 ) ∈ Φ(γ 1 , γ 2 ).

6 Existence of a solution to the contact problem Theorem 6.1 Under the assumptions of Section 3 there exists a solution of the problem P 1 v .

Proof. We shall prove that there exists at least a solution (u, λ + , λ -) of the problem P 2 v which will provide a solution (u, λ) of the problem P 1 v with λ = λ +λ -.

We apply Theorem 5.2 to

, we obtain (4.1).

The condition meas(Γ α U ) > 0, the ellipticity properties of the coefficients A α ijkl , B α ijkl and the Korn's inequality imply that there exist m α a , m α b > 0 such that

and we obtain

where m a = min(m 1 a , m 2 a ), m b = min(m 1 b , m 2 b ). Also, the properties (4.3)-(4.7), (4.9) and (4.13) can be easily verified. Now, let C tr be a positive constant such that v (L 2 (Ξ)) d ≤ C tr v H ι for all v ∈ H ι . Using (3.3), the following estimates hold:


Jérémy Leconte
email: jeremy.leconte@lmd.jussieu.fr

Francois Forget

Benjamin Charnay

Robin Wordsworth

Alizée Pottier

Increased insolation threshold for runaway greenhouse processes on Earth like planets

Because the solar luminosity increases over geological timescales, Earth climate is expected to warm, increasing water evaporation which, in turn, enhances the atmospheric greenhouse effect. Above a certain critical insolation, this destabilizing greenhouse feedback can "runaway" until all the oceans are evaporated [START_REF] Simpson | Some studies in terrestrial radiation[END_REF][START_REF] Komabayashi | Discrete equilibrium temperatures of a hypothetical planet with the atmosphere and the hydrosphere of one component-two phase system under constant solar radiation[END_REF][START_REF] Ingersoll | The Runaway Greenhouse: A History of Water on Venus[END_REF][START_REF] Kasting | Response of earth's atmosphere to increases in solar flux and implications for loss of water from Venus[END_REF] . Through increases in stratospheric humidity, warming may also cause oceans to escape to space before the runaway greenhouse occurs [START_REF] Kasting | Runaway and moist greenhouse atmospheres and the evolution of earth and Venus[END_REF][START_REF] Kopparapu | Habitable Zones around Main-sequence Stars: New Estimates[END_REF] . The critical insolation thresholds for these processes, however, remain uncertain because they have so far been evaluated with unidimensional models that cannot account for the dynamical and cloud feedback effects that are key stabilizing features of Earth's climate. Here we use a 3D global climate model to show that the threshold for the runaway greenhouse is about 375 W/m 2 , significantly higher than previously thought [START_REF] Kopparapu | Habitable Zones around Main-sequence Stars: New Estimates[END_REF][START_REF] Goldblatt | Low simulated radiation limit for runaway greenhouse climates[END_REF] . Our model is specifically developed to quantify the climate response of Earth-like planets to increased insolation in hot and extremely moist atmospheres. In contrast with previous studies, we find that clouds have a destabilizing feedback on the long term warming. However, subsident, unsaturated regions created by the Hadley circulation have a stabilizing effect that is strong enough to defer the runaway greenhouse limit to higher insolation than inferred from 1D models. Furthermore, because of wavelength-dependent radiative effects, the stratosphere remains cold and dry enough to hamper atmospheric water escape, even at large fluxes. This has strong implications for Venus early water history and extends the size of the habitable zone around other stars.

Planetary atmospheres naturally settle into a thermal equilibrium state where their outgoing thermal emission balances the heating due to sunlight absorption. The resulting climate is stabilized by the fact that a temperature increase results in an enhanced thermal-emission cooling. When a condensable greenhouse gas is present at the surface, such as water on Earth, this stabilizing feedback is somewhat hampered by the destabilizing greenhouse feedback: evaporation, and thus water vapor greenhouse effect, increases with temperature, reducing the cooling. Fortunately, under present Earth conditions, this greenhouse feedback is both strong enough to maintain clement surface temperatures and weak enough for the climate to remain stable.

When solar heating becomes stronger, however, water vapor can become abundant enough to make the atmosphere optically thick at all thermal wavelengths [START_REF] Goldblatt | Low simulated radiation limit for runaway greenhouse climates[END_REF][START_REF] Goldblatt | The runaway greenhouse: implications for future climate change, geoengineering and planetary atmospheres[END_REF] . Thermal flux then originates from the upper troposphere only and reaches a maximum, ∼ 282 W/m 2 , independently from the surface temperature [START_REF] Nakajima | A study on the ?runaway greenhouse effect? with a one-dimensional radiative-convective equilibrium model[END_REF] . If the planet absorbs more than this critical flux, thermal equilibrium can be restored only by vaporizing all the water available and reaching high surface temperatures at which the surface starts to radiate at visible wavelengths [START_REF] Kasting | Response of earth's atmosphere to increases in solar flux and implications for loss of water from Venus[END_REF][START_REF] Goldblatt | Low simulated radiation limit for runaway greenhouse climates[END_REF] . This is the runaway greenhouse state.

Because it has mostly been studied through unidimensional atmosphere models, the aforementioned mechanism strongly relies on the assumption that the troposphere is saturated in water vapor. Furthermore, by construction, these studies could not account for spatial inhomogeneities in insolation and in the resulting water vapor and cloud distributions [START_REF] Kasting | Response of earth's atmosphere to increases in solar flux and implications for loss of water from Venus[END_REF][START_REF] Kopparapu | Habitable Zones around Main-sequence Stars: New Estimates[END_REF][START_REF] Goldblatt | Low simulated radiation limit for runaway greenhouse climates[END_REF] . To overcome these limitations, we have developed a three-dimensional (3D) global climate model fit to describe hot atmospheres in which water vapor can become a dominant species [START_REF] Wordsworth | Gliese 581d is the First Discovered Terrestrial-mass Exoplanet in the Habitable Zone[END_REF][START_REF] Forget | 3D modelling of the early martian climate under a denser CO2 atmosphere: Temperatures and CO2 ice clouds[END_REF][START_REF] Leconte | 3D climate modeling of close-in land planets: Circulation patterns, climate moist bistability, and habitability[END_REF][START_REF] Charnay | Exploring the faint young Sun problem and the possible climates of the Archean Earth with a 3D GCM[END_REF] (see Methods). The main challenges of such a model are threefold. First, the radiative transfer must be fast, yet able to describe accurately the spectroscopic properties of various gases in a wide temperature-pressure domain. Second, the modeling of the physical processes that are not specific to a given planet (convection, turbulence, etc) must rely on the fewest free parameters possible to ensure its validity in stringent conditions. Finally, for hot, moist atmospheres, the description of the water cycle (moist convection, cloud formation, etc) must take into account the fact that water can become a major constituent of the atmosphere. For these reasons, usual global climate models used to predict Earth climate are generally not suited for such studies.

Here we perform simulations of future Earth climate by running our baseline model for various (increasing) values of the solar constant until radiative balance is achieved. For the current solar flux (F ≈ 341 W/m 2 ), our generic model reproduces the energetic budget and the characteristics of our climate 13 (see Fig. 1). When the flux is increased, the planet undergoes a decrease in surface albedo which is due to the melting of the permanent polar ice caps and the reduced seasonal snow cover. Above ∼350 W/m 2 , only seasonal ice caps appear during the polar night. The amount of water vapor also increases. This results in a more efficient absorption of the incoming stellar light, but also in an enhanced greenhouse effect which tends to homogenize the surface temperatures. While continental surfaces can reach temperatures around 100 • C because of the intense solar and greenhouse heating, sea surface temperatures remain moderate with a small diurnal variation because they are thermodynamically controlled ). As visible in the hot case, unsaturated subtropical regions allow the atmosphere to emit more than this limit.

by latent heat cooling [START_REF] Pierrehumbert | Principles of Planetary Climate[END_REF] . Finally, above ∼375 W/m 2 , no thermal equilibrium can be found. Although surface temperature increases with time, thermal emission reaches a limiting value. This is the onset of the runaway greenhouse instability. This runaway greenhouse limit arises at higher fluxes than recently found by previous unidimensional studies [START_REF] Kopparapu | Habitable Zones around Main-sequence Stars: New Estimates[END_REF][START_REF] Goldblatt | Low simulated radiation limit for runaway greenhouse climates[END_REF] and confirmed by our 1D model (see Extended Data Fig. 1). To understand the mechanisms deferring runaway greenhouse on actual planets, we first analyze the radiative effect of clouds. While 1D simulations cannot capture spatial variations in cloud distribution properly, it has been suggested that clouds should have a stabilizing feedback against runaway greenhouse effect [START_REF] Kasting | Response of earth's atmosphere to increases in solar flux and implications for loss of water from Venus[END_REF][START_REF] Kasting | Runaway and moist greenhouse atmospheres and the evolution of earth and Venus[END_REF][START_REF] Kopparapu | Habitable Zones around Main-sequence Stars: New Estimates[END_REF] . This tentative conclusion was based on the fact that on present Earth, the net cloud radiative forcing is negative, meaning that the albedo increase due to low level clouds exceeds the greenhouse effect of high level clouds [START_REF] Ramanathan | Cloud-Radiative Forcing and Climate: Results from the Earth Radiation Budget Experiment[END_REF] . Several authors thus proposed that, because of the increased evaporation resulting from the warming, cloud thickness would increase and enhance cloud stabilizing effect [START_REF] Kasting | Response of earth's atmosphere to increases in solar flux and implications for loss of water from Venus[END_REF][START_REF] Kasting | Runaway and moist greenhouse atmospheres and the evolution of earth and Venus[END_REF][START_REF] Kopparapu | Habitable Zones around Main-sequence Stars: New Estimates[END_REF] .

As shown in Fig. 2c where we show the evolution of the radiative cloud forcing with insolation, our simulations suggest the opposite. This is due to a displacement of the cloud formation region toward higher altitudes (see Fig. 3). As a result, the temperature at the mean cloud emission level increases much less with insolation than the surface temperature. Even though the cloud optical depth increases, the greenhouse feedback of the clouds exceeds their albedo effect. At higher fluxes, clouds become thinner, which reduces both longwave and shortwave radiative cloud forcing. For the reason described above, however, the greenhouse effect of clouds prevails and the net radiative forcing tends to cancel out.

There are several reasons for the cloud vertical displacement. First, moist convection, and Hadley circulation become more intense as the insolation increases, extending the troposphere. Second, to form and persist, clouds need to be able to get rid of the latent heat released during condensation. Because of the infrared opacity increase of the atmosphere, the level at which clouds can efficiently cool radiatively rises. This may explain both the progressive disappearance of low level clouds and the small cloud-deck temperature change seen in Fig. 3.

Although the tendencies described above should be robust, the actual precise value of the cloud radiative forcing does depend on the assumptions made on the cloud microphysics. In our baseline model, for instance, we assume that the number density of cloud condensation nuclei, i.e. the number of cloud particles, per unit mass of air remains fixed to a value that is representative of modern Earth. As is expected, a larger mass mixing ratio of condensed water thus entails bigger cloud particles which precipitate more easily and have a smaller radiative effect. To explore the possible behavior of clouds and bracket the reality, we have conducted another set of simulations where it is the radii of cloud particles which is kept constant. This assumption results in smaller radii, and thus overestimates cloud optical depths, and shortwave and longwave cloud forcing. As can be seen in Extended Data Fig. 2, longwave and shortwave forcing indeed increase continuously with insolation. However, for the reasons mentioned above, the greenhouse effect of clouds eventually overcomes the albedo effect. Therefore, the fact that the cloud feedback is destabilizing under an extreme insolation seems a robust conclusion which comes to support results obtained in the context of anthropogenic global warming [START_REF] Soden | An Assessment of Climate Feedbacks in Coupled Ocean Atmosphere Models[END_REF] .

If cloud do not help stabilize climate against a solar flux increase, atmospheric dynamics does. This is due to the fact that when a parcel of moist air is heated without any source of moisture, its water vapor pressure decreases relative to the saturation pressure. As a result, Earth troposphere itself is not saturated everywhere, unlike what is often assumed in 1D models. Sub tropical regions receiving hot air from the Hadley cell that has been dried during its ascent in the tropics and compressed adiabatically during its descent, are a perfect example [START_REF] Sherwood | Relative humidity changes in a warmer climate[END_REF] (see Fig. 4). As they stabilize Earth tropics today, such dynamically unsaturated regions where water vapor greenhouse effect is reduced stabilize climate against runaway greenhouse by playing the role of "radiative fins" where the emission can exceed the maximum emission for a saturated atmosphere [START_REF] Pierrehumbert | Thermostats, Radiator Fins, and the Local Runaway Greenhouse[END_REF][START_REF] Ishiwatari | A Numerical Study on Appearance of the Runaway Greenhouse State of a Three-Dimensional Gray Atmosphere[END_REF] (see Fig. 1).

Furthermore, as seen in Fig. 4, these unsaturated regions have the interesting property of extending upward and poleward with the Hadley cell (see Fig. 4), and to get dryer when the insolation is increased [START_REF] Sherwood | Relative humidity changes in a warmer climate[END_REF] (see Extended Data Figs. 3 and4). Such unsaturated regions, that cannot be predicted by a 1D model, ensure that the infrared photosphere remains always lower than in the saturated case, yielding a more efficient cooling to space. This stabilizing feedback is the main reason why the runaway greenhouse limit predicted by 3D simulations is closer to the Sun than previously found. And it would be even closer if it were not for the positive cloud feedback (see Methods).

With this new estimation, the inner edge of the habitable zone for Earth like planet in the Solar System is pushed inward to ∼0.95 AUs which means that the Earth should not enter a runaway greenhouse state before at least 1 billion years [START_REF] Gough | Solar interior structure and luminosity variations[END_REF] . However, the question of whether or not the warming-induced increase of the stratospheric temperature and humidity could cause the loss of the oceans by atmospheric escape before the runaway greenhouse is triggered, as suggested by some 1D studies [START_REF] Kasting | Response of earth's atmosphere to increases in solar flux and implications for loss of water from Venus[END_REF][START_REF] Kopparapu | Habitable Zones around Main-sequence Stars: New Estimates[END_REF] , is still subject to debate. Our simulations tend to invalidate this scenario. This is due to both non-gray radiative effects [START_REF] Pierrehumbert | Principles of Planetary Climate[END_REF][START_REF] Arking | The influence of line shape and band structure on temperatures in planetary atmospheres[END_REF][START_REF] Wordsworth | Water loss from terrestrial planets with CO2-rich atmospheres[END_REF] (see Methods) and unsaturated regions that flatten the thermal profile in the troposphere. As a result, stratospheric temperature are much colder than anticipated and stratospheric humidity cannot reach the threshold needed for efficient water photodissociation and hydrogen escape to space (see Fig. 3).

While extending the habitable zone [START_REF] Kasting | Habitable Zones around Main Sequence Stars[END_REF] towards the star, more importantly, our results highlight the fact that global climate models are needed to understand subtle climate feedbacks resulting from the inhomogeneous insolation of planetary surfaces that cannot be modeled in 1D. Although adding complexity and uncertainties, these subtle effects must be accounted for in modeling real planets, and especially tidally-synchronized exoplanets where they are even more pronounced [START_REF] Leconte | 3D climate modeling of close-in land planets: Circulation patterns, climate moist bistability, and habitability[END_REF][START_REF] Pierrehumbert | A Palette of Climates for Gliese 581g[END_REF][START_REF] Yang | Stabilizing Cloud Feedback Dramatically Expands the Habitable Zone of Tidally Locked Planets[END_REF] . In particular, while our simulation suggest that Venus, if it had Earth rotation rate, would have been in a runaway greenhouse state since its formation, its slow retrograde motion (because it results in a long Solar day), or a small water inventory [START_REF] Leconte | 3D climate modeling of close-in land planets: Circulation patterns, climate moist bistability, and habitability[END_REF][START_REF] Abe | Habitable Zone Limits for Dry Planets[END_REF] , could completely change the picture.

Methods Summary

Our simulations were performed with an upgraded version of the LMD generic global climate model (GCM) specifically developed for the study of extrasolar planets [START_REF] Wordsworth | Gliese 581d is the First Discovered Terrestrial-mass Exoplanet in the Habitable Zone[END_REF][START_REF] Leconte | 3D climate modeling of close-in land planets: Circulation patterns, climate moist bistability, and habitability[END_REF][START_REF] Selsis | Thermal phase curves of nontransiting terrestrial exoplanets. I. Characterizing atmospheres[END_REF] and paleoclimates [START_REF] Forget | 3D modelling of the early martian climate under a denser CO2 atmosphere: Temperatures and CO2 ice clouds[END_REF][START_REF] Wordsworth | Global modelling of the early martian climate under a denser CO2 atmosphere: Water cycle and ice evolution[END_REF] . The model uses the 3D dynamical core of the LMDZ Earth GCM used in IPCC studies [START_REF] Hourdin | The LMDZ4 general circulation model: climate performance and sensitivity to parametrized physics with emphasis on tropical convection[END_REF] , based on a finite-difference formulation of the primitive equations of geophysical fluid dynamics.

General physical processes relevant for present Earth -including ground thermal inertia and albedo (see Extended Data Fig. 5), turbulent transport, dry convection, evaporation/condensation and precipitations -are parametrized in the most physically based way to ensure the robustness of the model under extreme conditions [START_REF] Charnay | Exploring the faint young Sun problem and the possible climates of the Archean Earth with a 3D GCM[END_REF] . For the present study, special care has been taken to treat properly situation where water vapor can become a major constituent of the atmosphere. As detailed in the online-only Methods, the radiative transfer and moist convection scheme have been validated against recent 1D results in the hot, pure water atmosphere regime (see Extended Data Fig. 1). A numerical scheme accounting for change in atmospheric mass and surface pressure with water vapor evaporation/condensation has also been implemented.

Methods 1 Numerical climate model

Our simulations were performed with an upgraded version of the LMD generic global climate model (GCM) specifically developed for the study of extrasolar planets [START_REF] Wordsworth | Gliese 581d is the First Discovered Terrestrial-mass Exoplanet in the Habitable Zone[END_REF][START_REF] Leconte | 3D climate modeling of close-in land planets: Circulation patterns, climate moist bistability, and habitability[END_REF][START_REF] Selsis | Thermal phase curves of nontransiting terrestrial exoplanets. I. Characterizing atmospheres[END_REF] and paleoclimates [START_REF] Forget | 3D modelling of the early martian climate under a denser CO2 atmosphere: Temperatures and CO2 ice clouds[END_REF][START_REF] Wordsworth | Global modelling of the early martian climate under a denser CO2 atmosphere: Water cycle and ice evolution[END_REF] . The model uses the 3D dynamical core of the LMDZ Earth GCM used in IPCC studies [START_REF] Hourdin | The LMDZ4 general circulation model: climate performance and sensitivity to parametrized physics with emphasis on tropical convection[END_REF] , based on a finite-difference formulation of the primitive equations of geophysical fluid dynamics.

Generalities

A horizontal resolution of 64 × 48, corresponding to resolutions of 3.75 • latitude by 5.625 • longitude, is used for the simulations. The vertical grid uses hybrid coordinates: a terrainfollowing σ coordinate system in the lower atmosphere (σ being equal to the pressure divided by the surface pressure), and pressure levels in the upper atmosphere. In this work, we used 30 layers, with the lowest mid-layer level at 4 m and the top mid layer level at 1 mb. The dynamical core is called every 90 s, the tendencies given by the physical parametrizations are updated every 15 min and the radiative transfer is called every 90 min.

The boundary layer dynamics are accounted for by the unstationary 2.5-level closure scheme of ref 30 , plus a convective adjustment which rapidly mixes the atmosphere in the case of dry unstable temperature profiles. Turbulence and convection mix energy (potential temperature), momentum (wind), and water (condensed and gaseous). A standard roughness coefficient of z 0 = 10 -2 m is used for both rocky and ocean surfaces for simplicity.

The evolution of surface temperature is governed by the balance between radiative, latent, and sensible heat fluxes (direct solar insolation, thermal radiation from the atmosphere and the surface, and turbulent and latent heat fluxes; see Sect. 1.3) and thermal conduction in the soil. The parameterization of this last process is based on an 13-layer soil model solving the heat diffusion equation using finite differences. The depth of the layers were chosen to capture diurnal thermal waves as well as the deeper annual thermal wave. A vertically homogeneous soil is assumed. The thermal inertia of continental surface is set to 2 000 Js -1/2 m -2 K -1 . To model the high thermal inertia due to mixing in the upper layers of the oceans, the thermal inertia of the oceans is set to 18 000 Js -1/2 m -2 K -1 . The albedo surface map used in this work is shown in Extended Data Fig. 5 1

.2 Radiative transfer

The method used to produce our radiative transfer model is similar to [START_REF] Wordsworth | Gliese 581d is the First Discovered Terrestrial-mass Exoplanet in the Habitable Zone[END_REF][START_REF] Leconte | 3D climate modeling of close-in land planets: Circulation patterns, climate moist bistability, and habitability[END_REF] . For a gaseous composition similar to the Earth (1 bar of N 2 with 376 ppmv of CO 2 and a variable amount of water vapor; CH 4 , O 2 , and O 3 have been discarded for more generality), we computed high-resolution spectra over a range of temperatures and pressures using the HITRAN 2008 database [START_REF] Rothman | The HITRAN 2008 molecular spectroscopic database[END_REF] . For this study we used temperature and pressure grids with values T = {110, 170, ..., 710} K, p = {10 -3 , 10 -2 , ..., 10 5 } mbar. The H 2 O volume mixing ratio could vary in the range {10 -7 , 10 -6 , ..., 1}. The H 2 O lines were truncated at 25 cm -1 , while the water vapor continuum was included using the CKD model [START_REF] Clough | Line shape and the water vapor continuum[END_REF] . We also account for opacity due to N 2 -N 2 collision-induced absorption [START_REF] Borysow | Collision-induced rototranslational absorption spectra of N2-N2 pairs for temperatures from 50 to 300 K[END_REF][START_REF] Richard | New section of the hitran database: Collision-induced absorption (cia)[END_REF] .

The correlated-k method was then used to produce a smaller database of coefficients suitable for fast calculation in a GCM. Thanks to the linearity of the Schwarzschild equation of radiative transfer, the contribution of the thermal emission and downwelling stellar radiation can be treated separately, even in the same spectral channel. We therefore do not assume any spectral separation between the stellar and planetary emission wavelengths. For thermal emission, the model uses 19 spectral bands, and the two-stream equations are solved using the hemispheric mean approximation [START_REF] Toon | Rapid calculation of radiative heating rates and photodissociation rates in inhomogeneous multiple scattering atmospheres[END_REF] . Absorption and scattering of the downwelling stellar radiation is treated with the δ -Eddington approximation within 18 bands. Sixteen points are used for the ḡ-space integration, where ḡ is the cumulated distribution function of the absorption data for each band. Rayleigh scattering by N 2 and H 2 O molecules is included using the method described in ref 36 with updated cross section for water [START_REF] Bucholtz | Rayleigh-scattering calculations for the terrestrial atmosphere[END_REF] .

Water cycle

In the atmosphere, we follow the evolution of water in its vapor and condensed phases. These tracers are advected by the dynamical core, mixed by turbulence and dry and moist convection. Much care has been devoted to develop a robust and numerically stable water cycle scheme that is accurate both in the trace gas (water vapor mass mixing ratio q v 1) and dominant gas (q v ∼ 1) limit. In particular, the atmospheric mass and surface pressure variation (and the related vertical transport of tracers through pressure levels) due to any evaporation/sublimation/condensation of water vapor is taken into account (see next section).

Cloud formation

Cloud formation is treated using the prognostic equations of ref 38. For each column and level, this scheme provides the cloud fraction f c and the mass mixing ratio of condensed water q c , which are both functions of q v and the saturation vapor mass mixing ratio q s . In addition, when part of a column reaches both 100% saturation and a super-moist-adiabatic lapse rate, moist convective adjustment is performed following ref 39, and the cloud fraction is set to unity. This moist convection scheme has been chosen instead of more refined ones because it is more robust for a wide range of pressure, at the cost of giving enhanced precipitations at the equator [START_REF] Frierson | The Dynamics of Idealized Convection Schemes and Their Effect on the Zonally Averaged Tropical Circulation[END_REF] . Furthermore, the moist-adiabat lapse rate has been generalized to account for the fact that water can be a dominant species. This yields

∂ ln T ∂ ln p moist = p p -p v (1 -q v ) R a + q v L v T q v c p,v + q a c p,a + q c c p,c + q v L v T p p-p v d ln p s d ln T , (1)
where p i , M i , c p,i , and R i ≡ R /M i are respectively the pressures, molar masses, specific heat capacity at constant pressure, and specific gas constants of the various phases (non condensable gas, or air, denoted with a subscript a; condensible gas, or vapor, denoted by v; condensed material, denoted by c); p and T being the total pressure and temperature, R the molar ideal gas constant, L v the specific latent heat of vaporization and p s the water vapor saturation pressure.

In our baseline model, the liquid/icy water (depending on the temperature) is assumed to condense on a number N c of activated cloud condensation nuclei (CCN) per unit mass of moist air. This number density of CCNs is assumed to be spatially uniform but to have a different value for liquid or ice cloud particles (the condensed water phase is determined by the local temperature). Indeed, on Earth, nucleation is found to be much less efficient in cold, high level clouds, resulting in larger particle radii in these ice clouds compared to liquid water clouds. This dichotomy is essential to recover the observed balance between longwave and shortwave cloud radiative forcings. To recover current Earth climate, the values used are 4 × 10 6 kg -1 for liquid water clouds and 2 × 10 4 kg -1 for ice cloud.

Then, the effective radius of the cloud particles is given by

r eff = 3 q c/ 4 π ρ c N c 1/3
, where ρ c is the density of condensed water (10 3 kg/m 3 for liquid and 920 kg/m 3 for ice). In the special case where radii have been kept fixed, the values used where 12 and 38 µm for water droplets and ice particles respectively. Precipitations are computed with the scheme given in ref 41. Since this scheme explicitly considers the dependence on gravity, cloud particle radii, and background air density, it should remain valid for a wide range of situations (see ref 13 for details). Finally, the total cloud fraction of the column is assumed to be the cloud fraction of the level with the thickest cloud, and radiative transfer is computed both in the cloudy and clear sky regions. Fluxes are then linearly weighted between the two regions.

Hydrology and evaporation

The ground is modeled using a simple bucket model with a maximal water capacity of 150 kg/m 2 . When the water amount exceeds this limit, the surplus is regarded as runoff and added to the oceans. The effect of latent heat release during solidification/melting of snow/ice at the surface is taken into account. On the surface, ice can also have a radiative effect by linearly increasing the albedo of the ground to A ice = 0.5 until the ice surface density exceeds a certain threshold (here 30 kg m -2 ).

Evaporation E (in kg/m 2 /s) is computed within the boundary layer scheme, using a bulk aerodynamic formula multiplied by a dryness coefficient (β which is zero for a dry surface and linearly increases towards 1 when the water amount reaches 75 kg/m 2 at the surface; β =1 over oceans). This yields

E = ρ CV β q s (T surf ) -q 1 v ( 2
)
where ρ is the mass density of air, V the wind speed above the surface, q s (T surf ) the water vapor mass mixing ratio at saturation at the surface, and q 1 v the mixing ratio in the first layer. The aerodynamic coefficient is given by C = κ/ ln(1 + z 1 /z 0 ), where κ = 0.4 is the Von Karman constant, z 0 the roughness and z 1 the altitude of the first level.

1.4 A numerical scheme for atmospheric mass redistribution due to the condensation of a non-trace gas.

In most GCM's, when a trace gas condenses in a grid box, the variation of the total gas mass is neglected and the mass mixing ratio variation is given by δ q v = δ m v /m g . Because water vapor is not a trace gas in our simulations, we must take into account this effect as well as the change of surface pressure (total mass of the atmosphere) entailed by the massive evaporation/condensation that can take place both at the surface and aloft. To that purpose, we developed a numerical scheme similar to the one used for CO 2 clouds on Mars [START_REF] Forget | CO2 Snowfall on Mars: Simulation with a General Circulation Model[END_REF] . This scheme can be used for any kind of condensing species without any assumptions on the mixing ratio of the vapor phase. The scheme is thus valid both in the trace and dominant gas limit as well as in intermediate regimes. This is done in two steps. First, the various routine describing the water cycle (evaporation at the surface, moist convection, cloud condensation, precipitations) compute the mass of vapor that has been added at a level k, δ m k v (>0 when vapor is created [START_REF] Simpson | Some studies in terrestrial radiation[END_REF] ), assuming that the gas mass in the layer is constant (usual approximation). Second, the change in the total pressure and mass mixing ratios due to the change in the total mass of gas (δ m k g ) is computed as follows. For an atmospheric column of N layers and area A, the change in the surface pressure is given by

δ p s = g A N ∑ k=0 δ m k v , ( 3
)
where g is the gravity. For the computation in each cell, the difficulty comes from the fact that we do not use lagrangian coordinates. As a result, evaporation/condensation occurring both at the surface and aloft induces artificial mass fluxes through horizontal levels that must be accounted for.

For hybrid coordinates, as used in many GCMs, layer k encompasses the matter between p k-1 2 and p k+ 1 2 , where p k ≡ σ k p s + γ k . The so-called σ coordinates are retrieved by setting σ k ≡ p k /p s and γ k = 0. Because both σ k and γ k are time independent, the variation in the gas mass in the layer during one timestep is given by

δ m k g = A g δ p k-1 2 -δ p k+ 1 2 = A g σ k-1 2 -σ k+ 1 2 δ p s . (4)
Besides, the variation in the gas mass in layer k is also linked to the mass flux through its interfaces (W k-1 2 being the mass flux in kg between level k -1 and k counted positive when 1 δ m 0 v is the mass of water vapor evaporated directly at the ground going upward) by

δ m k g = δ m k v +W k-1 2 -W k+ 1 2 (5)
Using Eqs (3), ( 4) and ( 5), we can get a recursive formula for the fluxes

W k+ 1 2 = W k-1 2 + δ m k v -σ k-1 2 -σ k+ 1 2 N ∑ k=0 δ m k v ( 6
)
Because W N+ 1 2 ≡ 0, we can compute the other fluxes, and especially W

1 2 = δ m 0 v , since σ 1 2 ≡ 1.
Once the total mass fluxes are known, one can get the variation of the mixing ratios of the various tracers (generically called q here) by considering the tracer mass budget

δ m k g q k = q k-1 2 W k-1 2 -q k+ 1 2 W k+ 1 2 + ε δ m k v , ( 7
)
where ε is equal to 1 if the tracer considered is the vapor phase of the condensing gas, -1 if we consider the condensed phase, and 0 for every other tracers. The q k-1 2 are the tracer mixing ratios transported through the interface σ k-1 2 . They are calculated using a "Van-Leer I" finite volume transport scheme [START_REF] Van Leer | Towards the ultimate conservative difference scheme: IV. A new approach to numerical convection[END_REF] .

Alternatively, the tracer mass conservation can be written exactly

δ m k g q k = m k g + δ m k g δ q k + q k δ m k g , (8)
where δ q k is the correction to be applied to the tracer mixing ratios. Combining Eqs. ( 5), (7), and (8), one can show that

δ q k = 1 m k g +δ m k g × q k-1 2 -q k W k-1 2 -q k+ 1 2 -q k W k+ 1 2 + ε -q k δ m k v , (9)
where the first terms represent transport through σ levels and the last term corresponds to enrichment/depletion due to the variation of the gas mass. Finally, apart from tracers, the mass of gas that is advected through levels also transports some energy and momentum. Once the mass fluxes have been computed, these exchanges are easily computed using

δ u k = 1 m k g +δ m k g u k-1 2 -u k W k-1 2 -u k+ 1 2 -u k W k+ 1 2 , (10)
where we have shown an example for the zonal speed u, but equations are the same for other quantities.

Model validation at high temperatures

To validate both the radiative transfer and the implementation of the moist adiabat in very hot atmospheres, we have developed a 1D "inverse climate modeling" version of our code, that has already been used in ref 22. In its spirit, this model is very similar to the codes developed in refs 4 and 6.

For a given surface temperature and background gas surface pressure, the vertical thermal and water vapor profiles are integrated upward following a moist adiabatic lapse rate (Eq. ( 1)) until a fixed stratospheric temperature is reached (here 200K). Once the profile is built, the two stream radiative transfer routine used in our full 3D GCM is used to compute the outgoing thermal radiation and the bond albedo of the planet. In these calculations, the most important assumptions are that 1) the planet is spherically symmetric with a fixed surface albedo, 2) the atmosphere is always saturated in water vapor, and 3) the atmosphere is cloud free.

The results of this model are shown in Extended Data Fig. 1 for a surface albedo of 0.25 to be comparable with the recent results of ref 7. We can see that the asymptotic behavior of the thermal flux as well as its quantitative asymptotic value (282 W/m 2 ) are in good agreement with recently published similar models with up-to-date spectroscopic data [START_REF] Kopparapu | Habitable Zones around Main-sequence Stars: New Estimates[END_REF][START_REF] Goldblatt | Low simulated radiation limit for runaway greenhouse climates[END_REF][START_REF] Pierrehumbert | Principles of Planetary Climate[END_REF] . Our albedo calculations show a little discrepancy of less than 0.02 with respect to the relevant cases of ref 7 (i.e. the pure water case and the transitional atmosphere case with 1b N 2 and preindustrial CO 2 ). These small differences are of little importance given the uncertainties on the primary mechanism controlling the albedo, i.e. the clouds.

Intrinsic three-dimensional effects: the role of unsaturated regions

To have more insight into the intrinsic differences between 1D and 3D simulations, we have performed a set of idealized simulations. In this numerical experiment, we run our 3D model in a configuration without topography and with a uniform surface albedo (The effect of ice albedo is not taken into account). Furthermore, we do not take into account the radiative effect of clouds. Because we want to understand the role of the dynamically driven distribution of humidity, we treat the whole surface as an infinite reservoir of water. In these simulations, we are left with only one major free parameter (the background atmosphere being kept similar to the Earth present atmosphere, see above) which is the surface albedo that has been fixed to 0.22 to recover a mean surface temperature similar to the Earth one under the present insolation. The value of the mean surface temperature as a function of insolation is presented in Extended Data Fig. 3a. For comparison, we have computed the surface temperature given by our 1D model for the same surface albedo. Because we do not take into account both cloud forcing and ice albedo in this set of 3D simulations, the only differences between the two models are due to dynamical effects and horizontal inhomogeneities in vapor and temperature distri-butions. One can readily see from Extended Data Fig. 3a that the 3D model always predict lower temperatures than the 1D one. As a result, runaway greenhouse occurs at much higher insolation.

The explanation to this fact can be found in Extended Data Figs. 3 and4. Indeed, as observed on Earth, subtropical regions of the troposphere are unsaturated [START_REF] Sherwood | Relative humidity changes in a warmer climate[END_REF][START_REF] Pierrehumbert | Thermostats, Radiator Fins, and the Local Runaway Greenhouse[END_REF] . As a result, these regions can emit more thermal flux than the asymptotic limiting flux (see Extended Data Fig. 4b). At a given mean surface temperature, the global thermal flux can thus be larger than in the fully saturated 1D case, as can be seen in Extended Data Fig. 3b. To retrieve comparable results with the 1D model, we have to decrease the relative humidity in the radiative transfer to 0.6 for present earth and ∼0.45 near the runaway greenhouse limit. Thus both the fact that subtropics are unsaturated and that they seem to get dryer under an increased insolation strongly stabilize the climate. Thus, as on Earth today, subtropics will continue to play the role of radiative fins [START_REF] Pierrehumbert | Thermostats, Radiator Fins, and the Local Runaway Greenhouse[END_REF] , deferring catastrophic consequences of a runaway greenhouse further in the future.

In the context of the runaway greenhouse, the effect of these unsaturations have been recognized in ref 19 where the authors performed a highly-idealized numerical experiment. They found that their 3D results could be mimicked in 1D by forcing a relative humidity of 0.6 or a little lower throughout the atmosphere. However, because of their use of a cloud-free, gray radiative transfer and of a very crude description of the water cycle, they were unable to make a quantitative assessment of the runaway greenhouse insolation limit for a more realistic cloudy, non-gray atmosphere. Indeed, this limit is directly determined by the gray opacity used, which is a free parameter. In terms of relative humidity, however, these authors find values and trends that are roughly similar to ours. Considering the very different setup involved (no ground thermal inertia, no seasonal or diurnal variations, ...), this suggests that, at least when looking at rough mean values, the mechanisms controlling the evolution of the relative humidity should be controlled by relatively simple processes [START_REF] Sherwood | Relative humidity changes in a warmer climate[END_REF] . On the contrary, even if their very different temperature profiles, due to the gray approximation (especially their warm stratospheres), should not affect the runaway greenhouse threshold, they certainly prevent any robust conclusion considering water escape.

Finally, another result of our simulations is to confirm that both cloud and ice albedo feedbacks are destabilizing in the runaway greenhouse context. Indeed, because we keep the surface albedo fixed to a higher value than Earth surface today to mimic both the present effect of ice and clouds, this set of simulations can be seen as a case where clouds and ice are present but have no feedback. In this case, as visible on Extended Data Fig. 3a, runaway occurs at a flux greater than ∼ 400 W/m 2 which is much higher than the ∼ 375 W/m 2 found when the cloud feedback is accounted for.

Stratospheric temperatures in non-gray atmospheres

One might be surprised by the fact that the stratospheric temperature found in our baseline model (see Fig. 3) can be much colder than the skin temperature, T skin = T eq /2 1/4 , where T eq ≡ [(1 -A)F /σ ] 1/4 is the equilibrium temperature, A the bond albedo, F the mean insolation and σ the Stefan-Boltzman constant.

The reason for the existence of a skin temperature in gray atmospheres is that there is, by construction, no radiative window by which the lower atmosphere and surface can cool to space without interacting with the upper atmosphere. The radiative equilibrium thus implies that the upward and downward emission from the optically thin upper atmosphere must both be equal to half the absorption of upward infrared flux, explaining the usual 2 -1/4 factor.

As extensively discussed in ref 14, in the real case, the radiation illuminating the upper atmosphere is depleted in the portion of the spectrum which is efficiently absorbed by the gas in the troposphere. The stratosphere, which, by definition, can efficiently absorb the upwelling infrared radiation only in this part of the spectrum, is thus poorly heated from below. However, it is still able to radiate efficiently, and is thus forced to cool to maintain radiative equilibrium balance. This does not threaten the global radiative balance of the whole atmosphere because even if the upwelling flux is small in the opaque regions of the troposphere, the surface and lower troposphere can cool through the transparent radiative windows.

A first quantitative, analytical estimate of the magnitude of this effect can be done using an idealized 1 band gas. As described in ref 14, this simple model already shows that the ratio T skin /T eq can be lower than inferred in the gray case. A much more comprehensive quantification of this phenomenon using an idealized (in the sense that they can arbitrarily choose the line shape) semi-analytical non-gray radiative model has been provided in ref 21. The main conclusions of their study was that in a non-gray atmosphere, the temperature increases with increasing mean opacity below a certain height and decreases with increasing mean opacity above that height. They also found that there is no lower limit to the temperature at the top of the atmosphere: it can arbitrarily approach zero as the width of the lines is decreased.

While the aforementioned studies demonstrate that there is no theoretical paradox in having (arbitrarily) cold stratospheres, let us now turn to the validation of our more realistic radiative transfer model in the case of our ozone-free Earth-like planet. Recently, several different climate models (mainly addressing early Earth climate), have computed temperature profiles for an ozone-free atmosphere under present insolation. Among these, ref 44 find a minimum stratospheric temperature of 140-145K, and ref 45 find it to be around 150 K. In a more idealized, 1D setup, ref 22 have published consistently computed 1D temperature profiles for N 2 -CO 2 -H 2 O atmospheres and they also find stratospheric temperatures around 140-150 K. Considering the fact that our estimation for the ozone-free stratospheric temperature (∼140 K) agrees well within the uncertainties with these numerous published results using different numerical models, we are confident that our non-gray radiative transfer model is suitable for the present application. Emitted thermal flux (a), bond albedo (b), and effective solar constant (with respect to the current solar constant; c) dependence on surface temperature with a 1D version of our GCM in the "reverse climate modeling" mode [START_REF] Kasting | Response of earth's atmosphere to increases in solar flux and implications for loss of water from Venus[END_REF][START_REF] Kopparapu | Habitable Zones around Main-sequence Stars: New Estimates[END_REF] . The dashed curve is the pure water case, and the solid curve is the case with a 1 bar N 2 background atmosphere including 376 ppmv of CO 2 . The surface albedo is 0.25. In all panels, filled dots stand for the idealized 3D set of aquaplanet simulations, and gray curves stand for the 1D model. In both cases, a uniform surface albedo of 0.22 is used. In the 1D case, the relative humidity in the radiative transfer is forced to be 1, 0.6 and 0.45 (solid, dashed and dotted curves respectively). Extended Data Figure 5: Surface albedo map used for the Earth baseline case. This map does not include the effect of the ice albedo which is computed directly by the GCM. The albedo of Greenland and Antartica, in particular, was set to 0.35. The altitude of these regions was, however, left unchanged, explaining in part the temperature contrasts around these areas in Fig. 1.

Extended Data Figure

Figure 1 :

 1 Figure1: Temperature and radiative budget for the Earth under two insolations. Map of the annual mean surface temperature for the model corresponding to present Earth (F = 341 W/m 2 ; a) and to a mean solar flux of 375 W/m 2 (b), just before runaway greenhouse is triggered. c and d: Zonally and annually averaged surface temperature (solid black), and absorbed (gray dashed) and emitted (gray dotted) fluxes. The red horizontal line shows the radiation limit on the emitted flux for a saturated water atmosphere (282 W/m 2 ). As visible in the hot case, unsaturated subtropical regions allow the atmosphere to emit more than this limit.

Figure 2 :

 2 Figure 2: Evolution of the mean surface temperature, planetary albedo, and cloud radiative forcing with the mean solar incoming flux. Curves start form present Earth conditions (∼ 341 W/m 2 ). a) The solid curve stands for the average surface temperature in our 3D baseline model and the dashed curve for the temperature in the 1D cloud-free model. b) The solid line is the average surface albedo and the dashed line is the effective planetary albedo (or Bond albedo). c) Solid, dotted and dashed curves are the shortwave, long wave and net radiative cloud forcing in the baseline model.

Figure 3 :

 3 Figure 3: Evolution of globally averaged vertical profiles. Mean vertical profiles of temperature (a), water vapor (b), and condensed water mixing ratios (c; in kg per kg of moist air) for 4 different insolations (solid: 341 W/m 2 ; long dashed: 353 W/m 2 ; dashed: 365 W/m 2 ; dotted: 375 W/m 2 ).

Figure 4 :

 4 Figure 4: Meridional distribution of relative humidity. Zonally and annually averaged relative humidity (see color bar) for models receiving a stellar flux equal to 341, 365 and 375 W/m 2 (from panel a to c). White contours are streamlines showing the Hadley circulation (solid and dashed curves depict clockwise and counter-clockwise rotation, respectively).

1 :

 1 Validation of the radiative transfer model at high temperature.

2 :

 2 Evolution of the cloud radiative forcing with the mean solar incoming flux for the scenario with fixed cloud particle radii. Solid, dotted and dashed curves are the shortwave, longwave and net radiative cloud forcing. Although less simulations have been run, the changes in the value of the slopes around 353 and 365 W/m 2 seem to have the same origin as the behavior change seen in Fig.2c(although they occur at different fluxes). These changes in cloud behavior might be due to the disappearance of both permanent ice caps (at lower fluxes) and seasonal snow covers (at higher fluxes).

m 2 c

 2 Extended DataFigure 3: Comparison between 1D and 3D cloud-free aquaplanet simulations. a, Mean surface temperature as a function of incoming stellar flux. b, Emitted thermal flux dependence on surface temperature. c, Water vapor column as a function os surface temperature.

F 375 W m 2 b

 2 Extended Data Figure4: Relative humidity and radiative budget for an idealized, cloudfree aquaplanet. a, Annually and zonally averaged relative humidity in a latitude altitude plane. b, Absorbed (gray dashed) and emitted (gray dotted) flux and surface temperature (solid black) distribution with latitude (annually and zonally averaged). The red line is the asymptotic limit infrared flux for a saturated atmosphere. Results are shown for the case receiving 375 W/m 2 .

Acknowledgments

We thank our referees, J. Kasting and Y. Abe, for their thorough review, and A. Spiga and F. Selsis for discussions. This work was supported by grants from Région Ile-de-France

Author contributions

J.L. developed the "high temperature/humidity" version of the generic global climate model, performed the calculations, and led the analysis and writing of the results. F.F. initiated the development of the generic global climate model and provided critical advice during analysis and writing. B.C worked on the development of the model and helped performing the comparison with current Earth climatology. RW developed the original version of the generic model, including the correlated-k radiative transfer scheme. A.P. performed comparison runs and sensitivity studies. All the authors commented on the manuscript.

Author Information

Author Information Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests. Readers are welcome to comment on the online version of the paper. Correspondence and requests for materials should be addressed to J.L. (jeremy.leconte@lmd.jussieu.fr).