Towards robust viscoelastic-plastic-damage material model with different hardenings / softenings capable of representing salient phenomena in seismic loading applications

Keywords: Thermodynamics with internal variables, Phenomenological approach, Continuum-discrete plasticity and damage models, Generalized standard material, Earthquake Engineering. 1

This paper presents the physical formulation of a 1D material model suitable for seismic applications. It is written within the framework of thermodynamics with internal variables that is, especially, very efficient for the phenomenological representation of material behaviors at macroscale: those of the representative elementary volume. The model can reproduce the main characteristics observed for concrete, that is nonsymetric loading rate-dependent (viscoelasticity) behavior with appearance of permanent deformations and local hysteresis (continuum plasticity), stiffness degradation (continuum damage), cracking due to displacement localization (discrete plasticity or damage). The parameters have a clear physical meaning and can thus be easily identified. Although this point is not detailed in the paper, this material model is developed to be implemented in a finite element computer program. Therefore, for the benefit of the robustness of the numerical implementation, (i) linear state equations (no local iteration required) are defined whenever possible and (ii) the conditions in which the presented model can enter the generalized standard materials class -whose elements benefit from good global and local stability properties -are clearly established. To illustrate the capabilities of this model -among them for Earthquake Engineering applications -results of some numerical applications are presented.

Introduction

The development of numerical computing facilities within the last decades has made it possible to carry out numerical large scale dynamic nonlinear analyses that take into account the nonlinear behavior of the materials, such as steel and concrete in reinforced concrete structures [START_REF] Bhattacharjee | Seismic cracking and energy dissipation in concrete gravity dams[END_REF][START_REF] Davenne | Efficient approaches to finite element analysis in earthquake engineering[END_REF][START_REF] Martinelli | Simulation of the shaking table test of a sevenstory shear wall building[END_REF][START_REF] Mousseau | Seismic performance of a full-scale, reinforced high-performance concrete building -Part II: analytical study[END_REF]. Multi-scale-like stategies for defining frame section behavior laws from information collected at the material level also emerged, e.g. for static ultimate load design [START_REF] Pham | Stress resultant model for ultimate load design of reinforced concrete frames -Combined axial force and bending moment[END_REF]. The implementation of robust refined material models in computer programs has thus become a key issue in structural engineering. Indeed, the prediction of the timehistory of some material properties provides useful information to predict the global structural seismic response: strength and stiffness degradations as well as residual deformation in the materials help for drawing conclusions about the residual capacity of the structure in post-earthquake conditions. In Earthquake Engineering, the physical modeling of damping is not an easy task because many dissipative phenomena occur at several levels. Refined material models are required for a reliable computation of the amount of material energy dissiped in seismic loading. Other dissipation sources at the level of the bounding between steel and concrete in RC structural elements have already been introduced in numerical static analyses [START_REF] Domínguez | Enhanced solid element for modelling of reinforced concrete structures with bond-slip[END_REF][START_REF] Ibrahimbegovic | On modeling of fire resistance tests on concrete and reinforced concrete structures[END_REF] and one can hope that they will be adapted to seismic analyses within the next few years. Since many computations are required for complete seismic analyses (need to consider several input motions) or in automatic parameters identification procedures [START_REF] Kucerova | Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures Part II: identification from tests under heterogeneous stress field[END_REF], the robustness of the material models should be considered with care.

This paper deals with modeling the behavior of concrete until collapse in dynamic loading applications so that it is robust once implemented in a computer program. The authors already presented a robust elasto-plastic-damage model and its implementation in finite element procedures [START_REF] Ibrahimbegovic | Coupled plasticity-damage model and direct stress interpolation[END_REF]. In addition here, in the physical formulation: (i) viscosity is added by appealing to the so-called Kelvin-Voigt viscoelastic model; (ii) the modeling of the softening part, which is introduced in the formulation by a discrete plastic or damage model that requires the definition of an enhanced strain field, is presented in detail and (iii) a clear distinction between tensile and compressive constitutive models is made. On the numerical side, the conditions under which the model can enter the generalized standard materials class [START_REF] Halphen | On generalized standard materials (in French)[END_REF], whose elements benefit from good robustness properties, are identified. This latter point, along with the model formulation that allows solving the local problem (the integration point level) without any iteration, except in tension, in a numerical nonlinear structural analysis [START_REF] Ibrahimbegovic | Coupled plasticity-damage model and direct stress interpolation[END_REF] justifies the use of the words "towards robust" in this paper title.

This material model is formulated in the framework provided by thermodynamics with internal variables (TIV ), which is very suitable for developing phenomenological material constitutive models with good numerical properties [START_REF] Maugin | The thermodynamics of nonlinear irreversible behaviors an introduction[END_REF]. A set of internal variables is chosen to describe nonlinear phenomena that occur in the material. Our approach is phenomenological because the nonlinear phenomena we aim at modeling for concrete are those observed on the strain-stress experimental macroscopic response of a concrete representative elementary volume (REV). However, internal variables often are connected to microscopic phenomena that have macroscopic effects such as in steel materials where microscopic dislocations motion lead to a macroscopic strain hardening behavior. More accurate representations of concrete could be provided for numerical simulations by, for instance, taking into account the heterogeneous nature of concrete [START_REF] Hautefeuille | Failure model for heterogeneous structures using structured meshes and accounting for probability aspects[END_REF]. The macroscopic phenomenological approach is however an efficient way to carry out large-scale structural numerical simulations with a rather refined material description.

The introduction of viscosity in a material model can be motivated by several purposes. (i) Viscosity in concrete can be associated to its loading ratedependent macroscopic behavior that results from several nano-and microscopic phenomena such as complex interaction between moisture and micro-structural solid skeleton, micro-cracking process, Stefan effect, micro-inertia of the material surronding the crack tip [START_REF] Pedersen | An analysis of dynamic fracture in concrete with a continuum visco-elastic visco-plastic damage model[END_REF]. (ii) Viscosity can also be associated to creep phenomena in highly dissipative materials such as asphalt concrete [START_REF] Panoskaltsis | The generalized Kuhn model of linear viscoelasticity[END_REF]. (iii) The introduction of viscosity leads to the explicit appearance of a time-step in the governing equations of the discretized problem. The possibility to determine a critical time-step then provides a tool to reguralize numerical problems that can be ill-posed in the presence of materials that exhibit strain softening behaviors [START_REF] Oliver | Theoretical and computational issues in modelling material failure in strong discontinuity scenarios[END_REF]. (iv) Viscosity can also be used to introduce a source of energy dissipation at the material level: it can lump the energy dissipation that physically comes from unknown or not (accurately) modeled physical nonlinear phenomena in the materials; this is what the authors have in mind when developing this concrete constitutive model.

The aim of this paper is to present the theoretical formulation of a refined 1D cyclic constitutive concrete model that has been developed in a framework leading to a robust implementation into finite element numerical procedures. We first present in section 2 the theoretical formulation of the constitutive concrete model; in section 3, the conditions in which the model can enter the generalized standard material class are examined for the benefit of robustness; in section 4, the mechanical problem of a concrete 1D structure in dynamic loading is formulated in a complementary variational form and the Euler-Lagrange equations are derived; before to conclude, section 5 is dedicated to the presentation of some numerical examples that illustrate some capabilities of the constitutive concrete model implemented in a finite element program.

2 Constitutive concrete model

Thermodynamics with internal variables

To represent the salient phenomena observed in the experimental cyclic response of a concrete REV (Fig. 1), a set α of internal variables is chosen. The phenomena we aim at modeling are: (i) loading-rate dependent behavior, (ii) strain hardening, (iii) strain softening, (iv) appearance of residual deformation, (v) stiffness and strength degradations, (vi) hysteresis loops. We consider that the phenomena (ii), (iv) and (vi) only occur in compression, that strain softening in compression is due to the localization of permanent deformation whereas strain softening in tension is due to the localization of deformation that completely disappear after unloading. Internal variables are associated to each of these phenomenologically identified mechanisms.

The set of internal variables we choose is defined in Tab. 1. • and • refer to continuum and discrete quantities. Concerning the continuum internal Fig. 1: [left] Experimental (adapted from [START_REF] Ramtani | Contribution to the modeling of the multi-axial behavior of damaged concrete with description of the unilateral characteristics[END_REF]) and [right] numerical (by using the proposed constitutive concrete model) cyclic behavior of concrete in quasi-static loading.

variables: ǭv and ǭp are the viscous and plastic deformations, ξp and λp are the strain-like variables that represent isotropic and kinematic plastic strain hardenings, D represents stiffness degradation and ξd isotropic damage strain hardening. Then, concerning the discrete internal variables, that is in the section where displacement localizes, ūp represents plastic deformation, ξp plastic strain softening in compression, D the stiffness degradation and ξd damage strain softening in tension.

t/c compression tension

α ǭv ǭp ξp λp D ξd ūp ξp D ξd A σ v σ qp κp σ 2 2 qd t qp t 2 2 qd PC E, η Ē Kp Hp Ē-1 Kd Ē Kp Ē-1 σ ∞ , a
Tab. 1: Set of internal variables α, corresponding affinities A and phenomenological coefficients PCs (parameters of the model).

A noteworthy advantage of this material model is that all its parameters have a clear meaning and can therefore be easily identified according to experimental curves. The complete procedure to identify a concrete law is presented in Section 5.1: the viscous parameter η can be seen as a material property and identified from free vibration tests; the tensile softening curve coefficient a is related to the fracture energy G F , a fracture mechanics concept [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF] that quantifies the total amount of energy that has to be furnished in tension to a concrete section between the time when displacement begins to localize and the time when the section is completely broken; all the other parameters can be identified from stress-strain curves of concrete specimen in quasi-static cyclic loading.

Governing equations of the constitutive model 2.2.1 Basic ingredients

We first define the three basic ingredients for developing a constitutive model in the framework of the TIV :

1. Split of the total strain into a viscoelastic, a plastic, a damage and another discrete part due to the likely localization of the displacement. A-space is devided into two parts: the terms involved in tension σ ≥ 0 are mentioned by the sign • + and those involved in compression σ < 0 by the sign • - (note that for a 1D case, tension and compression have a clear meaning):

u(x, t) = ū(x, t) + ū(t)H x(x) (1)
ǫ(x, t) = ∂u ∂x = ǭv + ǭp + ( ū⋄ + ūp )δ x where ǭ = ∂ ū ∂x = ǭv + ǭp is the continuum part of the strain and ūδ x = ( ū⋄ + ūp )δ x -with ū⋄ = ūd in tension and ū⋄ = ūe in compression -its discrete part introduced by the Dirac's function δ x that is the derivative of the Heaviside's function H x that introduces the displacement jump ū in the description of the displacement. ǭv corresponds to the deformation in a Kelvin-Voigt rheological model: a spring -elastic in tension (ǭ v = ǭe ) and damaging in compression (ǭ v = ǭd ) -and a dashpot in parallel. For simplicity, we assume here that the displacement jumps only take place in one section; extension to several localization sections is straightforward.

2. Helmholtz free energy functional:

ψ + (u, α) = ψe (ū + ) + ψd ( ū+ , D) + Ξd ( ξd ) δ x (2)
ψ -(u, α) = ψd (ū -, ǭp , D) + Ξp ( ξp ) + Λp ( λp ) + Ξd ( ξd ) + ψe ( ū-, ūp ) + Ξp ( ξp ) δ x 3. The nonlinear plasticity and damage mechanisms are activated for a positive value of the corresponding plasticity and damage criteria:

• For σ ≥ 0 (tension), the dissipative mechanisms that can be involved in the evolution process are activated according to the following criteria:

φv,+ (σ) = σ v = σ -σ e ≤ 0 (3) φd,+ (t, qd ) = t -(σ t u -qd ) ≤ 0
where t is the stress at the discontinuity and σ t u is the ultimate stress in tension.

• For σ < 0 (compression), the criteria functions are:

φv,-(σ) = -σ v = -σ -σ d ≤ 0 (4) φp,-(σ, qp , κp ) = |σ + κp | -(σ y -qp ) ≤ 0 φd,-(σ, qd ) = -σ -(σ f -qd ) ≤ 0 φp,-(t, qp ) = -t -(σ c u -qp ) ≤ 0
where σ y , σ f and σ c u are yield, fracture and ultimate stresses in compression.

State equations of the system

We can now give the expression of the energy dissipated by this viscoelasticplastic-damage model with different hardenings / softenings within a time unit. This can be computed as the difference between the total amount of energy imparted to the system and the amount of energy stored by the system (internal potential energy) during this time unit:

In tension:

Ḋ+ = σ ǫ+ + t u+ δ x -ψ+ (ū + , ū+ , D, ξd ) = (σ e + σ v ) ǫ+ + t u+ δ x - ∂ ψe ∂ǭ + ǫ+ - ∂ ψd ∂ ū+ u+ + ∂ ψd ∂ D Ḋ + ∂ Ξd ∂ ξd ξd δ x
We also have ǫ+ = ǫe = ǫv and u+ = ud because ūp cannot evolve in tension according to Eqs. ( 3) and ( 4). Thus:

Ḋ+ = σ e - ∂ ψe ∂ǭ e ǫe + σ v,+ ǫv,+ + t - ∂ ψd ∂ ūd ud δ x - ∂ ψd ∂ D Ḋδ x - ∂ Ξd ∂ ξd ξd δ x (5)
In compression:

Ḋ-= σ ǫ-+ t u-δ x -ψ-ū-, ū-, ǭp , ξp , λp , D, ξd , ūp , ξp = σ d + σ v ǫ-+ t u-δ x - ∂ ψd ∂(ǭ --ǭp ) ( ǫ--ǫp ) + ∂ ψd ∂ D Ḋ + ∂ Ξp ∂ ξp ξp + ∂ Λp ∂ λp λp + ∂ Ξd ∂ ξd ξd - ∂ ψe ∂( ū--ūp ) ( u--up ) + ∂ Ξp ∂ ξp ξp δ x
We also have ǫ-= ǫv + ǫp with ǫv = ǫd and u-= ue + up . Thus:

Ḋ-= σ d - ∂ ψd ∂(ǭ --ǭp ) ǫd + σ v ǫv,-+ σ ǫp - ∂ ψd ∂ D Ḋ - ∂ Ξp ∂ ξp ξp - ∂ Λp ∂ λp λp - ∂ Ξd ∂ ξd ξd + t - ∂ ψe ∂( ū--ūp ) ue δ x + t up δ x - ∂ Ξp ∂ ξp ξp δ x (6)
Note that ǭe , ūd , ǭd , ūe are not internal variables. Moreover, when there is no evolution of any internal variable, the dissipation is null. Therefore, according to Eqs. ( 5) and ( 6), we have the following state equations:

σ e = ∂ ψe ∂ǭ e , σ d = ∂ ψd ∂ǭ d , t + = ∂ ψd ∂ ūd , t -= ∂ ψe ∂ ūe (7)
To recover, for the thermodynamic forces, expressions that will then lead to an accurate reproduction of the experimentally observed behavior of a concrete REV, we define the different terms in the expression of the Helmholtz free energy as follows:

ψe = 1 2 ǭe Ēǭ e ψd = 1 2 ǭ--ǭp D-1 ǭ--ǭp (8) ψd = 1 2 ūd D-1 ūd Ξp = 1 2 ξp Kp ξp Ξd = -σ t u -σ ∞ ξd + 1 a e -a ξd Λp = λp Hp λp Ξd = 1 2 ξd Kd ξd ψe = ū--ūp Ē ū--ūp Ξp = 1 2 ξp Kp ξp
We define the state equations of the internal variables from the Helmholtz free energy functional that plays the role of a thermodynamic potential:

A i = - ∂ψ ∂α i (9)
Finally, the state equations of the system, that is the equations that characterize the state of a concrete material point at a given instant in A-space, are given in Tab. 2.

Tension

Compression

σ e = Ēǭ e σ d = Dǭ d t = D-1 ūd t = Ē ūe qd = -∂ Ξd ∂ ξd = (σ t u -σ ∞ ) 1 -e -a ξd qp = -∂ Ξp ∂ ξp = -Kp ξp 1 2 t 2 = -∂ ψd ∂ D κp = -∂ Λp ∂ λp = -Hpλp qd = -∂ Ξd ∂ ξd = -Kd ξd 1 2 σ 2 = -∂ ψd ∂ D qp = -∂ Ξp ∂ ξp = -Kp ξp
Tab. 2: State equations of the system.

Equations of evolution of the internal variables

Now, we determine the equations that govern the evolution of the internal variables of this concrete phenomenological constitutive law. We assume that the evolution is driven by the principle of maximum dissipation. With the state equations in Tab. 2 at hand we first rewrite the dissipation -Eqs. ( 5) and ( 6

To characterize the evolution of the system, except the viscous one, we extend the principle of maximum dissipation -which was first formulated for plasticity problems [START_REF] Hill | The Mathematical Theory of Plasticity[END_REF][START_REF] Lubliner | Plasticity Theory[END_REF] -to our model and thus compute the flow of the internal variables by assuming that they maximize the dissipation. Remembering that the criteria functions in Eqs. ( 3) and ( 4) must be satisfied, we thus have to solve a problem of constrained maximization that involves inequalities. Denoting f = (σ, t, q p , κ p , q d , qp , qd ), we have the following problem P:

max f Ḋ with constraints φ p (f ) ≤ 0 , φ d (f ) ≤ 0 , φp (f ) ≤ 0 , φd (f ) ≤ 0 (11)
Such a problem can be solved with the Lagrange's multipliers method [START_REF] Ibrahimbegovic | Nonlinear solids mechanics -Theoretical formulations and finite element solution methods[END_REF][START_REF] Strang | Introduction to Applied Mathematics[END_REF]. We define the Lagrangian as:

L(f, γ) = Ḋ -( γp φ p + γd φ p + γp φp + γd φd ) (12)
where γ = ( γp , γd , γp , γd ) is the set of Lagrange's multipliers associated to each constraint. According to the Kuhn-Tucker conditions, if the set f * = {f * i } is a solution of P, then there exists a unique set of Lagrange multipliers γ * = { γ * j } such that the following relations are verified for all i and for all j:

∂L(f * , γ * ) ∂f i = 0 , ∂L(f * , γ * ) ∂ γj ≥ 0 , γ * j ≥ 0 , γ * j φ j (f * ) = 0 ( 13
)
where the last three conditions are referred to as the loading / unloading conditions. Finally, the equations of evolution of the internal variables are given in Tab. 3. For the viscous internal variable ǭv , the evolution is expressed, in an associative form too, à la Perzyna [START_REF] Perzyna | Fundamental problems in viscoplasticity[END_REF]:

ǫv = φv η ∂ φv ∂σ ( 14
)
where η is the viscous parameter assumed here to be rate-insensitive although this might not be the case [START_REF] Pedersen | An analysis of dynamic fracture in concrete with a continuum visco-elastic visco-plastic damage model[END_REF]. Note that Eq. ( 14) is in accordance with the rheological relation in a dashpot σ v = η ǫv . Note that concrete is in general not an associative material. However, in the 1D context of this paper, we nevertheless assume that it is the case. All the equations of evolution of the internal variables thus are expressed in an associative form. The normality of the flow rule to the loading surface is in particular implied by the appeal to the maximum dissipation principle.

Generalized standard materials

Materials that belong to the generalized standard materials class benefit from nice local and global stability properties required for robust implementation in, for instance, finite element procedures. For other materials, some attention must be paid on the local or global stability properties. There are five conditions, in an isothermal context, to be fulfilled for a material to enter the generalized standard materials class [START_REF] Maugin | The thermodynamics of nonlinear irreversible behaviors an introduction[END_REF]:

1. Its Helmholtz free energy ψ(α) is convex in α = {α i }.

2. The laws of state for the internal variables are derived from the Helmholtz free energy that plays the role of a thermodynamic potential:

A i = -∂ψ ∂αi .
3. The intrinsic dissipation verifies: Ḋ = A i αi ≥ 0.

4. There exists a closed convex set K in the space of the associated variables A and A ∈ K.

5. The flux of the internal variables α i is normal to the boundary of K (normality law).

As shown in section 2.2, the material model presented in this paper satisfies the conditions 2) and 5). The sets K mentioned in the condition 4) are characterized according to the criteria functions φ(A) introduced in Eq. ( 3) and ( 4) and their boundary is defined by the equations φ(A) = 0. In our 1D context, these criteria functions clearly define closed and convex sets in A-space. About condition 1), we recall that a function of the form f (x) = cx 2 is convex only if c > 0. The convexity of the Helmholtz free energy is therefore not verified for a negative phenomenological coefficient K, which is required to represent softening.

Let prove that our model satisfies the condition 3). First, since η > 0, Ḋv = σ ǫv = η ǫv2 ≥ 0. For the other dissipative mechanisms:

Ḋ⋄ = σ ǫ⋄ + q ⋄ ξ⋄ = γ⋄ (σ ∂φ ⋄ ∂σ + q ⋄ ) (see Tab. ( 3))
= γ⋄ (φ ⋄ + σ y,f,u ) (see Eqs. ( 3) and ( 4)) [START_REF] Ibrahimbegovic | Nonlinear solids mechanics -Theoretical formulations and finite element solution methods[END_REF] According to the loading / unloading conditions Eq. ( 13), γ⋄ ≥ 0 and γ⋄ φ ⋄ = 0. Moreover, the material parameters verifies σ y,f,u > 0 and thus, finally, Ḋ⋄ ≥ 0. Fulfillment of the condition 3) ensures local stability properties (hysteretic loops builded clockwise).

To conclude this section, (i) because condition 1) is not satisfied with K < 0, the global stability of the numerical model has to be investigated: the tangent stiffness matrix could become non-invertible during nonlinear resolution procedure; moreover, (ii) the constitutive model presented here is rate-sensitive and the existence of a critical time-step must thus also be investigated; finally, (iii) A-space is divided into two parts, one for tension and the other for compression, which shall be considered in the numerical resolution phase. In this paper, we only focus on the physical formulation of the model, but it can be shown that these three points can be treated in an efficient way that makes the numerical implementation in finite element procedures robust at local and global levels.

4 Euler-Lagrange equations of a concrete 1D structure in dynamic loading

In this section, we move up from the formulation of the local constitutive model to the formulation of the structural mechanical problem. This latter, illustrated here in the 1D case, has to be adapted to the local constitutive model and requires the definition of an enhanced displacement field.

Enhanced displacement field kinematics

x x ūū H x ũ ũ u Fig. 2: Enhanced displacement field in a 1D structure. The solid, dashed and dotted lines represent the enriched displacement field u(x, t), the smooth displacement field ũ(x, t) and the continuum displacement field ū(x, t).

As shown in Fig. 2, the displacement field is written as the sum of a smooth linear part ũ(x, t) enhanced by an additional part ũ(x, t) to take into account the possible appearance of a displacement jump ū(t) in the 1D structure [START_REF] Garikipati | A study of strain localization in a multiple scale framework -the one-dimensional problem[END_REF][START_REF] Ibrahimbegovic | Combined hardening and softening constitutive model of plasticity: precursor to shear slip line failure[END_REF]:

u(x, t) = ũ(x, t) + (H x(x) -ϕ(x)) ū(t) ũ(x,t) (16)
which can also be written, in the form adapted for identification with Eq. ( 1), as:

u(x, t) = ũ(x, t) -ϕ(x) ū(t) ū(x,t) + ū(t)H x(x) (17)
Then ǫ(x, t) = ǭ(x, t) + ū(t)δ x with:

ǭ(x, t) = ∂ ū(x, t) ∂x = ∂ ũ(x, t) ∂x - dϕ(x) dx ū(t) (18)
In this 1D context, only the failure mode I is considered; we refer to [START_REF] Dujc | An embedded crack model for failure analysis of concrete solids[END_REF] for a presentation of the kinematics with embedded crack of quadrilateral 2D finite elements where several failure modes are considered.

Complementary Lagrangian variational formulation

The kinetic energy, assumed to depend only on the smooth part of the displacement, is:

T (ũ, t) = 1 2 Ω ρ ∂ ũ(x, t) ∂t 2 dΩ (19)
where the volumic mass ρ in the volume Ω is considered as constant. The total potential energy can be written:

U (ũ, ū, α, t) = Ω ψ(ũ, ū, α, t)dΩ -U ext (ũ, t) (20)
To introduce the stress fields σ and t as other independent variables of the problem and thus give a more general setting that can be helpful to derive the governing equations of the problem [START_REF] Markovic | Complementary energy based FE modelling of coupled elasto-plastic and damage behavior for continuum microstructure computations[END_REF][START_REF] Simo | Complementary mixed finite element formulations for elastoplasticity[END_REF], we appeal to the partial Legendre transformation of the Helmholtz free energy:

ψ + = σ e ǭe - 1 2 σ e Ē-1 σ e + t ūd - 1 2 t Dt + Ξd δ x (21)
ψ -= σ d ǭd - 1 2 σ d Dσ d + Ξp + Λp + Ξd + t ūe - 1 2 t Ē-1 t + Ξp δ x
Recalling that ǭe = ǭv in tension, ǭd = ǭv in compression, σ = σ v + σ e/d and according to Eqs. ( 1) and ( 18), we rewrite these equations with respect to ũ(x, t), ū(t), σ(x, t), t(t) and α:

ψ + = σ ∂ ũ ∂x - dϕ dx ū -ǭp -σ v ǭv - 1 2 (σ -σ v ) Ē-1 (σ -σ v ) (22) + t( ū -ūp ) - 1 2 t Dt + Ξd δ x ψ -= σ ∂ ũ ∂x - dϕ dx ū -ǭp -σ v ǭv - 1 2 (σ -σ v ) D(σ -σ v ) + Ξp + Λp + Ξd + t( ū -ūp ) - 1 2 t Ē-1 t + Ξp δ x
We denote L the Lagrangian of the mechanical system. It is defined as L = T -U and Lagrange's variational principle can be written as:

t2 t1 δL(ũ, ū, σ, t, t)dt = 0 ∀t 1 , ∀t 2 (23)
with δũ(x, t), δ ū(t), δσ(x, t), δt(t) null when t = t 1 and t = t 2 .

After some calculus (integration by parts,. . . ), we obtain the following set of governing equations of the system available for tension and compression:

Ω ρ ∂ 2 ũ ∂t 2 δũdΩ + Ω ∂ψ ∂ ũ δũdΩ -δU ext = 0 ( 24
) Ω ∂ψ ∂ ū dΩδ ū = 0 Ω ∂ψ ∂σ δσdΩ = 0 Ω ∂ψ ∂t dΩδt = 0
which leads, with Γ denoting the section in which the displacement is likely to occur and by remembering that neither ū nor t depends on the position x (they are only defined at the position x Γ of the section Γ), to the set of Euler-Lagrange equations of a concrete 1D structure written in Tab. 4.

Tension Compression

Ω ρ ∂ 2 ũ ∂t 2 δũ + ∂δ ũ ∂x σ dΩ -δU ext = 0 Ω ρ ∂ 2 ũ ∂t 2 δũ + ∂δ ũ ∂x σ dΩ -δU ext = 0 Ω dϕ dx σ -tδ x dΩ δ ū = 0 Ω dϕ dx σ -tδ x dΩ δ ū = 0 Ω δσ ∂ ũ ∂x -dϕ dx ū -ǭv -ǭp dΩ = 0 Ω δσ ∂ ũ ∂x -dϕ dx ū -ǭv -ǭp dΩ = 0 Γ ū -ūp -Dt dΓ δt = 0 Γ ū -ūp -Ē-1 t dΓ δt = 0
Tab. 4: Euler-Lagrange equations of a concrete 1D structure.

Either for tension or compresssion, the first equation in Tab. 4 enforces satisfying global force equilibrium, the second equation gives the condition of compatibility between the continuum σ and discrete t stresses, the third and fourth equations correspond to the weak form of the local continuum and discrete constitutive models.

To give another -also simplified in our viewpoint -expression of the set of the governing equations of the problem, one can take advantage of the dependancy between the continuum and the discrete stresses that appears in Tab. 4:

Ω dϕ dx σ -tδ x dΩ = 0 ⇒ t = t(σ).
Then, the expression of the Lagrangian is modified according to:

Ω ψ + dΩ = Ω σ ∂ ũ ∂x -ǭv -ǭp + -t ūp - 1 2 t Dt + Ξd δ x dΩ Ω ψ -dΩ = Ω σ ∂ ũ ∂x -ǭv -ǭp + Ξp + Λp + Ξd + -t ūp - 1 2 t Ē-1 t + Ξp δ x dΩ (25)
Moreover, the response of the localization section is supposed to be infinitely rigid before the stress reaches an ultimate value and therefore Ē-1 → 0. Finally, denoting A Γ the area of the section Γ where displacement localization can take place, the simplified form of these equations is written in Tab. 5. Only two unknown fields remain: the smooth displacement ũ(x, t) and the continuum stress σ(x, t).

Tension Compression

Ω ρ ∂ 2 ũ ∂t 2 δũ + ∂δ ũ ∂x σ dΩ -δU ext = 0 Ω ρ ∂ 2 ũ ∂t 2 δũ + ∂δ ũ ∂x σ dΩ -δU ext = 0 Ω δσ ∂ ũ ∂x -ǭv -ǭp Ω δσ ∂ ũ ∂x -ǭv -ǭp -ūp dt dσ δ x dΩ = 0 -ūp + Dt dt dσ δ x dΩ = 0 with t = t(σ) = 1 AΓ Ω dϕ dx σdΩ with t = t(σ) = 1 AΓ Ω dϕ dx σdΩ
Tab. 5: Final Euler-Lagrange equations of a concrete 1D structure.

Numerical applications

The numerical implementation has been done in the general purpose finite element computer program FEAP [START_REF] Taylor | FEAP: A Finite Element Analysis Program, User manual and Programmer manual[END_REF]. Because any external distributed loading is involved (self-weight is in particular neglected) in these applications, the concrete structure is always modeled with constant-stress bar elements. The 1D structure tested is presented in Fig. 3. Neglecting the concrete mass m = ρAL as compared to the added mass M , we compute the stiffness k, the fundamental pulsation ω and period T of this structure as:

M L x u ext (t) or F ext (t)
k = ES L = 1.4 GPa ; ω = k M = 44.7 s -1 ; T = 2π ω = 0.14 s (26)

Identification of a concrete law

All the parameters of the 1D cyclic concrete model presented in this paper have a clear physical meaning. It is shown in this section that (i) the viscous parameter η is related to the material critical damping ratio ξ, (ii) the tensile softening law parameter a is related to the fracture energy G F and (iii) all the other parameters can be directly identified from the experimental stress-strain response of a concrete specimen in quasi-static cyclic loading.

Identification of the viscous parameter η

We show here that the viscous parameter η can be interpreted as a material property that can be identified according to, among others, the experimental results of free low-amplitude -so that the structure remains elastic -vibration tests (Fig. 4):

1. First, denoting ξ the critical damping ratio and c the structural viscous parameter, if one assumes that ξ << 1, one can write [START_REF] Chopra | Dynamics of structures: Theory and applications to earthquake engineering[END_REF]:

c = 2ξmω (27)
ξ = 1% ξ = 2% ξ = 5%
Fig. 4: Free-vibrations viscoelastic response of the tested structure for various critical damping ratios.

2. Then we link the material viscous parameter η to the structural one c starting from the 1D local form of the equilibrium and the Kelvin-Voigt model constitutive equations:

ρü + ∂σ ∂x = 0 ; σ = E ∂u ∂x + η d dt ∂u ∂x (28)
For a tested structure with constant cross-section area S and strain field ( ∂u ∂x = u L ), the introduction of Eq. ( 28) 2 into Eq. ( 28) 1 and then the integration of the resulting local equation over the whole structure gives:

ρSL m ü + S L η c u + ES L k u = 0 ( 29
)
3. The Kelvin-Voigt viscoelastic model exhibits a hysteretic -energy dissipative -behavior in cyclic loading as illustrated in Fig. 5 for an excitation of the form ū = A • sin (Ωt) with, except for other indicated values, A = 0.3 mm, Ω = 157.08 s -1 and τ = 0.001 s. The amount of energy dissipated per cycle D cyc -the area of the elliptical loop -can be analytically written as [START_REF] Wang | Intrinsic damping: modeling techniques for engineering systems[END_REF]:

D cyc = πcΩA 2 (30)
The main drawback of the Kelvin-Voigt model is that D cyc is dependent on the loading frequency Ω. It is indeed not realistic in particular for seismic excitations [START_REF] Wang | Intrinsic damping: modeling techniques for engineering systems[END_REF]. This dependency is only negligeable in the cases where the response is primarily represented by the resonant pulsation, that is Ω ≃ ω. 4. Now, suppose that a concrete structure is excited by an impulsion that is weak enough not to active the nonlinear phenomena that are represented by our model, and that, assuming that the response is primarily represented by the fundamental pulsation ω, a critical damping ratio ξ = 1 % is measured (for instance thanks to the logarithmic decrement), we can then see η as a material property and identify it according to the following relation:

η = 2ξmωL S (31)
In the case of a tested structure identical to those in Fig. 3, we compute:.

η = 15.7 MPa.s ; τ = η E ≃ 4.5e -4 s (32)

Identification of the tensile softening parameter a

The tensile softening parameter a can be computed according to the fracture energy G F . G F is one of the parameters that characterize the softening stressseparation (crack opening) curve of the cohesive crack model [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF][START_REF] Bazant | Choice of standard fracture test for concrete and its statistical evaluation[END_REF] and represents the total amount of energy that has to be furnished in tension to a section between the time t loc when displacement begins to localize and the time t cri that we define as the time when the local softening problem becomes (numerically) ill-posed, because Kd < 0; physically, t cri has to coincide with the time when the section is completely broken: when all the energy G F is consumed. In the case of our concrete model, one can demonstrate that t cri is always defined and that:

G F = tcri t loc t udt = σ t u 2a ⇒ a = σ t u 2G F (33)
In addition, there exists an empirical relation [START_REF] Bazant | Statistical prediction of fracture parameters of concrete and implications for choice of testing standard[END_REF] that allows computing G F according to material characteristics:

G F = 2.5 α 0 σ c u 0.051 0.46
1 + d a 11.27 0.22 w c -0.30 [START_REF] Strang | Introduction to Applied Mathematics[END_REF] with σ c u in MPa and where α 0 = 1 for rounded aggregates, α 0 = 1.44 for crushed or angular aggregates, d a is the maximum aggregate size in mm and w c is the water-cement ratio. Finally, for a concrete with α 0 = 1.44, d a = 25 mm and w c = 0.5, we compute G F = 136 N.m -1 that leads, from [START_REF] Simo | An analysis of strong discontinuity induced by strain softening solutions in rate-independent solids[END_REF]

with σ t u = 3 MPa, to: a = 22060 m -1 (35)

Identification of the other parameters

Here is the procedure to identify the remaining parameters of the model from the stress-strain experimental response in Fig. 1 [left]:

1. Identify the Young's modulus for tensile and compressive parts. We consider here that they have the same value: E = Ē = 35 GPa.

2. Identify the set of stresses σ y , σ f , σ c u , σ t u so as to characterize the changes of the slope of the backbone curve. We choose: σ y = 3 MPa, σ f = 37 MPa, σ c u = 49 MPa and σ t u = 0.6 • σ c u = 3 MPa.

3. Identify simultaneously Kp and Hp so as they describe both the strain hardening phase of the backbone curve in the range σ y < σ < σ f and the shape of the local hysteresis loops.

4. Identify Kd so as to describe the remaining strain hardening phase of the backbone curve in the range σ f < σ < σ c u .

5. Identify Kp so as to describe the softening part of the backbone curve.

Fig. 1 [right] is plotted with the parameters indicated above. Note that although all the parameters have a clear physical meaning, they also are all connected and the results of a first identification process thus often need to be refined. Indeed, the slope of the hardening part of the backbone curve in the range

σ y < σ < σ f is C 1 = Ē Kp Ē+ Kp , in the range σ f < σ < σ c u is C 2 = C1 Kd C1+ Kd
and the slope of the softening part of the backbone curve is C 3 = C2 Kp C2+ Kp . However, in the absence of an automatic identification process [START_REF] Kucerova | Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures Part II: identification from tests under heterogeneous stress field[END_REF], following the procedure indicated above gives satisfying results within a few iterations. Note also that, although we observed that activating plasticity before damage leads to a better identification, it is possible to invert the roles of σ y and σ f .

Mesh objectivity

In finite element procedures, when displacement localizes, the solution of the mechanical problem can depend on the size of the mesh. In the formulation of our model, we introduced a strong discontinuity -a displacement jump - [START_REF] Ibrahimbegovic | Combined hardening and softening constitutive model of plasticity: precursor to shear slip line failure[END_REF][START_REF] Simo | An analysis of strong discontinuity induced by strain softening solutions in rate-independent solids[END_REF] and thus concentrated the localization in a zero-length zone. This method leads to a formulation that does not need any characteristic length and thus satisfies the mesh objectivity requirement, that is unicity of the solution, in a bar as illustrated in Fig. 6 where the quasi-static tensile response of the section where displacement localizes in the bar in Fig. 3 submitted to an imposed displacement is plotted.

Concerning the modeling of the tensile softening phase, the fracture energy is so small with resect to the elastic energy that it is numerically difficult to grasp the post peak tensile response of concrete bar. The curve in Fig. 6 has been plotted with 4,000 points, which cannot be achieved for instance in seismic applications. Note that the presence of viscosity helps the numerical representation of this tensile softening part.

Loading rate-sensitive response

The four main effects of the loading strain-rate on the response of concrete are: the increases of (i) the compressive strength [START_REF] Bischoff | Impact behavior of plain concrete loaded in uniaxial compression[END_REF][START_REF] Lu | Modelling of dynamic behaviour of concrete materials under blast loading[END_REF], (ii) the tensile strength [START_REF] Weerheijm | Tensile failure of concrete at high loading rates: new test data on strength and fracture energy from instrumented spalling tests[END_REF], (iii) the Young's modulus [START_REF] Weerheijm | Tensile failure of concrete at high loading rates: new test data on strength and fracture energy from instrumented spalling tests[END_REF] and (iv) the brittle behavior [START_REF] Dilger | Ductility of plain and confined concrete under different strain rates[END_REF]. For seismic excitations, the strain rates are comprised between 10 -5 s -1 and 10 -2 s -1 , which respectively corresponds to an increase of (i) the compression strength between 0 and 30% [START_REF] Bischoff | Impact behavior of plain concrete loaded in uniaxial compression[END_REF][START_REF] Lu | Modelling of dynamic behaviour of concrete materials under blast loading[END_REF], (ii) the tensile strength between 0 and 60% [START_REF] Weerheijm | Tensile failure of concrete at high loading rates: new test data on strength and fracture energy from instrumented spalling tests[END_REF] and (iii) the Young's modulus between 0 and 10% [START_REF] Weerheijm | Tensile failure of concrete at high loading rates: new test data on strength and fracture energy from instrumented spalling tests[END_REF].

The numerical applications presented in Fig. 7 [left] show that the Young's modulus computed for a loading strain rate of 10 -2 s -1 is about 5.7% larger than the one computed for a loading strain rate of 10 -5 s -1 . The increase of Young's modulus is thus represented by the model but a little underestimated here. However, this increase is related to the viscous parameter τ computed in these numerical applications from the damping ratio ξ = 1% whose value has been arbitrarily chosen: one can thus hope that, in reality, the right value for the damping ratio would lead to a correct modeling of the increase of elastic stiffness. The compressive response is presented in Fig. 7 [right]. The proposed concrete model has not been developed yet to represent by itself the rate-dependent tensile and compressive strengths, what can be done a priori by hand, as it is common practice in engineering where the parameters that define the strengths for concrete in quasi-static loading are majorated for seismis analyses. We focus on the compressive part only of the concrete law in Fig. 1. Thus a static compressive normal force in imposed to the bar before the seimic excitation.

In the absence of steel, once the strength degradation process has begun, the structure is no more capable of resisting the seismic force demand and global equilibrium can no more be satisfied; we therefore couple to the concrete bar an elastic steel bar with cross section area S s = 0.02 • S c = 8e -4 m 2 and Young's modulus E s = 210 GPa. The behavior law for concrete is those presented in Fig. 1 [right]. The loading pattern is shown in Fig. 8 and the global structural response in Fig. 9.

Local response and intrinsic dissipated energies

TIV provides a useful framework to quantify the material intrinsic dissipation. Indeed, the amount of continuum viscous, plastic and damage energies dissipated in all the structure wide and the discrete plastic and damage energies dissipated in the section where displacement localizes can be computed according to Eq. ( 10): Note that the total intrinsic dissipation is composed by a volumic and a surfacic -in the section Γ where displacement localizes -part. The local response of material points located outside and inside Γ is ploted in Fig 10 : the localization of the displacement in Γ leads to a larger amount of dissipated energy in this section. The amount and sources of intrinsic energy dissipation within the seimic excitation are detailed in Fig. 11. Around t = 11 s, a big amount of seimic energy is imparted to the structure (see Fig. 8); so as this latter do not collapse, it has to dissipate this input energy. The input energy is converted both into stored energy in the material and dissipated energy (irreversible mechanisms).

Ēv = T 0 Ω Ḋv dΩ dt = T 0 Ω σ v ǫv dΩ dt (36)

Conclusion and future work

In this paper, the theoretical formulation of a 1D cyclic constitutive concrete model is presented. The model is capable of representing most of the salient phenomena that can be experimentally observed in the stress-strain response of a concrete representative elementary volume (REV) -at macroscale -in cyclic 1D loading: brittle loading-rate dependent behavior in tension, quasibrittle loading-rate dependent behavior in compression with strain hardening, appearance of residual deformation, stiffness and strength degradations, local hysteretis. Although the details of the numerical implementation are not given in this paper, the model has been developed to be implemented in finite element numerical procedures and thus special attention is paid on its robustness: (i) the conditions under which the model enters the generalized standard materials class are determined and (ii) linear hardening / softening laws that make the resolution processes noniterative at the local -or numerical integration point -level are prefered to nonlinear ones. The Euler-Lagrange equations of a 1D structure, made of this concrete material, in dynamic loading are derived and finally some numerical applications are presented. They illustrate that: the viscous parameter can be considered as a material property, the objectivity of the mesh of the tested bar is satisfied when displacement localizes and the material REV thus exhibits a softening behavior, the model can represent the salient phenomena that occur in concrete in cyclic loading. The model is developed in the framework of thermodynamics with internal variables, that is that a set of internal variables is chosen to represent the experimentally observed response of a concrete REV in 1D cyclic loading. All the material parameters introduced have a clear physical meaning, what makes their identification simple. The model is developed in a purely phenomenological way at macroscale: for instance "strain hardening" means here "observed increase of the yield stress while strain increases" and has no relation with physical effects that occur at lower scales such as dislocations motion well known for steel materials. One could give a more refined description of the micro-or even nanoscopic behavior of concrete and one could hope that relating physical microor nanoscopic effects to the macroscopic response of a concrete REV would still improve the accuracy of the material models. However, our approach, albeit phenomenological and thus somehow arbitrary, can be efficiently used in common computer programs for large-scale structural numerical dynamic analyses. To that purpose, the implementation of the model in a fiber beam element is currently being achieved.

  ∂φ p ∂σ Ḋσ = γd ∂φ d ∂σ ξp = γp ∂φ p ∂q p ξd = γd ∂φ d ∂q d λp = γp ∂φ p ∂κ p γp ≥ 0 , φ p ≤ 0 , γp φ p = 0 γd ≥ 0 , φ d ≤ 0 , γd φ d = 0 Discrete plasticity Discrete damage up = γp ∂ φp ∂t Ḋt = γd ∂ φd ∂t ξp = γp ∂ φp ∂ qp ξd = γd ∂ φd ∂ qd γp ≥ 0 , φp ≤ 0 , γp φp = 0 γd ≥ 0 , φd ≤ 0 , γd φd = 0 Tab. 3: Equations of evolution of the internal variables.

Fig. 3 :

 3 Fig. 3: 1D concrete structure model used for the numerical applications. The cross-section area S = 0.04 m 2 , the length L = 1 m, the added mass M = 7.0e 5 kg, the volumic mass ρ = 2400 kg.m -3 and the concrete elastic modulus E = 35 GPa.

1 Ω 1 Ω 1 Fig. 5 :

 1115 Fig. 5: Hysteretic response of the viscoelastic model in quasi-static sine loading for several [top-left] viscous parameters, [top-right] loading amplitudes and [bottom] forcing pulsations.

Fig. 6 :

 6 Fig.6: Tensile response of the section where the discontinuity appears. n el is the number of elements of the regular mesh of the bar.

  -5 s -1 ǫ = 5 • 10 -1 s -1 ǫ = 1s -1

Fig. 7 :

 7 Fig. 7: Tensile [left] and [right] compressive concrete response for several loading rates.

Fig. 8 :

 8 Fig. 8: Loading pattern: static loading + seismic loading time histories.

Fig. 9 :

 9 Fig. 9: [left] Right-node displacement and [right] reaction-displacement timehistories.

2 t

 2 Ḋt + qd ξd dΓ dt

Fig. 10 :Fig. 11 :

 1011 Fig. 10: Local response of a material point [left] located and [right] not located in section Γ where displacement localizes.

Acknowledgements

The financial support provided by the French Ministry of Research, the Quebec Fund for Research on Nature and Technology (FQRNT), the Natural Science and Engineering Research Council of Canada (NSERC), and ENS-Cachan Invited Professor Funds is acknowledged.


Armel Oumbe
email: armel.oumbe@mines-paristech.fr

Lucien Wald
email: lucien.wald@mines-paristech.fr

A PARAMETERISATION OF VERTICAL PROFILE OF SOLAR IRRADIANCE FOR CORRECTING SOLAR FLUXES FOR CHANGES IN TERRAIN ELEVATION

This paper deals with the modelling of vertical profile of solar irradiance for correcting solar radiation data. An approximation of the influence of the terrain elevation on all-skies radiation and a parameterisation of vertical profile on clear-skies radiation are proposed. This parameterisation, called double-z fitting function, is validated using two radiative transfer models. The main applications of this work are: computing clear-sky irradiances from satellite images, extrapolation of irradiances from one site to another of different altitude and correction of databases derived from satellite images.

NOMENCLATURE

z, z 0 , z H , z 1 , z 2 altitude (km). Altitude is defined above mean sea level I 0 irradiance at the top of the atmosphere (W m -2 ) I(z), B(z) global, direct irradiance at altitude z on horizontal plane (W m -2 ) I c (z), B c (z) global, direct irradiance at altitude z on horizontal plane under clear sky (W m -2 ) f I , f B vertical profile irradiance fitting functions for global, direct irradiance A(z), A B (z) attenuation coefficients of the atmosphere for global, direct irradiance (unitless) ,  B parameters of the double-z fitting function (km -1 ) K c clear-sky index. It is equal to the ratio of I to I c (unitless) vis visibility at ground level (km)

INTRODUCTION

The surface solar irradiance (SSI) -also called downwelling shortwave irradiance -is an important element of the surface radiation budget. It is seldom measured and several works have developed methods for the assessment of the surface solar irradiance from images taken by geostationary or polar-orbiting satellites [START_REF] Ba | Satellite-Derived Surface Radiation Budget over the African Continent. Part I: Estimation of Downward Solar Irradiance and Albedo[END_REF][START_REF] Pinker | Effects of spatial sampling of satellite data on derived surface solar irradiance[END_REF]. The various Heliosat methods [3 -10], the methods of Mueller et al. [START_REF] Mueller | Rethinking satellite based solar irradiance modelling -The SOLIS clear sky module[END_REF], Perez et al. [START_REF] Perez | A new operational model for satellite-derived irradiances description and validation[END_REF], Raschke et al. [START_REF] Raschke | Solar Radiation Atlas of Africa[END_REF] or Stuhlmann et al. [START_REF] Stuhlmann | An improvement of the IGMK model to derive total and diffuse solar radiation at the surface from satellite data[END_REF] are based on the principle that the irradiance at a given pixel can be computed by reducing the irradiance that should be observed under clear-sky at that pixel by a certain amount that is a function of the properties of the clouds. Such satellite-derived irradiance data is currently used in various domains, ranging from climate to materials weathering [START_REF] Gschwind | Converting a successful research project into a sustainable service: the case of the SoDa Web service[END_REF].

The computation of the SSI from satellite images calls upon a digital terrain model (DTM) whose cell size fits that of the pixel. For example, the ESRA atlas [START_REF] Esra ; Greif | European Solar Radiation Atlas[END_REF] or the HelioClim-1 database [START_REF] Lef ‡vre | Using reduced data sets ISCCP-B2 from the Meteosat satellites to assess surface solar irradiance[END_REF] exploits the DTM TerrainBase [START_REF] Terrainbase | Worldwide digital terrain data, Documentation Manual, CD-ROM Release 1.0[END_REF] whose cell size is 5' of arc angle, i.e. approximately 10 km at mid-latitude. The size of the cell is even larger for the NASA-SSE database: 1• of arc angle [START_REF] Chandler | NASA climatological data for renewable energy assessment[END_REF] or for the ISIS database: 280 km [START_REF] Lohmann | Long-term variability of solar direct and global irradiance derived from ISCCP data and comparison with re-analysis data[END_REF]. Generally, these sizes are too large to describe changes in altitude with a sufficient accuracy in areas of steep relief; large discrepancies can be found between the mean altitude of a cell and the altitude of a particular site within this cell. For example, in Switzerland, the altitude of the measuring station at Saentis Mountain is 2490 m while that provided by the DTM TerrainBase is 1126 m, i.e. an underestimation of 1364 m. If uncorrected, this difference leads to an error in the irradiance provided by the ESRA and HelioClim databases. Note that in very steep relief, irradiance depends upon shadows cast on the sites by surrounding obstacles and not only on change in altitude; this is not the subject of this paper.

Wahab et al. [START_REF] Wahab | Assessing surface solar irradiance in northern africa desert climate and its long-term variations from meteosat images[END_REF] established, for monthly irradiation, that when the difference in altitude reached 200 m, the error on estimated irradiation is significant and conclude that the correction for altitude should be applied in a systematic way. The main objective of this paper is to propose a parameterisation of the vertical variation of SSI with altitude z. Most models for SSI under clear-sky take into account change in SSI with z, whether they are in analytical form [START_REF] Esra ; Greif | European Solar Radiation Atlas[END_REF][START_REF] Iqbal | An Introduction to Solar Radiation[END_REF][START_REF] Gueymard | SMARTS2, A simple model of the atmospheric radiative transfer of sunshine: algorithms and performance assessment[END_REF] or in numerical form such as libRadtran, LOWTRAN or MODTRAN. Laue [START_REF] Laue | The measurement of solar spectral irradiance at different terrestrial altitudes[END_REF] performed measurements of solar spectral irradiance at different altitudes. Blumthaler et al. [START_REF] Blumthaler | Increase in solar UV radiation with altitude[END_REF], Dvorkin and Steinberger [START_REF] Dvorkin | Modelling the altitude effect on solar UV radiation[END_REF] focused their work on the UV part of the radiation. Hofierka and Suri [START_REF] Hofierka | The solar radiation model for Open source GIS: implementation and applications[END_REF] propose an empirical relationship between the Linke turbidity factor and the altitude; the corrected factor may then be used in a clear-sky model such as ESRA. No model describing in an explicit manner how the total irradiance changes with z in clear-or all-skies was found in the literature. For practical reasons, it would be convenient to have a simple parameterisation, preferably in analytical form. This would allow e.g., easy implementation in spreadsheet or fast calculation for demonstration and education purposes. It would help in benefiting of the availability of DTM at high spatial resolution such as SRTM: 30 or 90 m [START_REF] Van Zyl | The Shuttle Radar Topography Mission (SRTM): a breakthrough in remote sensing of topography[END_REF] to increase the accuracy of the assessment of SSI at a given site. Other applications are possible and three of them are discussed hereafter.

This paper proposes a parameterisation of the vertical profile of the global and direct irradiances. We firstly establish that, in first order, the influence of elevation on SSI in all-skies is equivalent to that in clear-skies. This parameterisation requires solar irradiance at two different altitudes as inputs that can be provided by codes simulating the radiative transfer in atmosphere (RTM: radiative transfer model) or other means. The parameterisation is validated by the results produced by the two RTM: libRadtran and 6S. These RTM are well known and serve as references. The RTM libRadtran was originally developed for modelling UV irradiance [START_REF] Bernhard | Comparison of measured and modeled spectral ultraviolet irradiance at Antarctic stations used to determine biases in total ozone data from various sources[END_REF][START_REF] Balis | Study of the effect of different type of aerosols on UV-B radiation from measurements during EARLINET[END_REF]. Its accuracy has also been demonstrated for actinic fluxes [START_REF] Kylling | Spectral actinic flux in the lower troposphere: measurement and 1-D simulations for cloudless, broken cloud and overcast situations[END_REF][START_REF] Thiel | Influence of clouds on the spectral actinic flux density in the lower troposphere (INSPECTRO): overview of the field campaigns[END_REF] and total irradiance [START_REF] Mueller | Rethinking satellite based solar irradiance modelling -The SOLIS clear sky module[END_REF][START_REF] Lohmann | Long-term variability of solar direct and global irradiance derived from ISCCP data and comparison with re-analysis data[END_REF][START_REF] Ineichen | Comparison of eight clear sky broadband models against 16 independent data banks[END_REF]. The RTM 6S is largely used in processing and simulation of satellites images [START_REF] Vermote | Second simulation of the satellite signal in the solar spectrum (6S), 6S: An overview[END_REF][START_REF] Gascon | Radiative transfer model for simulating high-resolution satellite images[END_REF]. Contrary to libRadtran, 6S deals only with clear-sky case.

APPROXIMATION OF INFLUENCE OF TERRAIN ELEVATION ON ALL-SKIES

The main factors influencing the SSI under clear-sky I c (z) received at an altitude z are the ground albedo and optical path. The optical path is mainly influenced by clouds, aerosols and water vapour. Note that this paper deals with total irradiance, i.e., the irradiance received for all wavelengths less than 4 µm. The adjective total is omitted in the following. The ground albedo has an impact on the diffuse component of the irradiance; the greater the albedo, the greater the diffuse component. This impact decreases as the altitude z increases. The optical path is a function of the altitude. As the altitude increases, the thickness of the part of the atmosphere under concern decreases, the optical path decreases leading to a decrease in the attenuation of the radiation, and the downwelling irradiance increases. Besides the purely geometrical effect induced by the reduction of the geometrical path length, the increase in irradiance with z is due to the fact that atmospheric components influencing the radiation under clear skies -water vapour, aerosols -are mostly located in the lower parts of the atmosphere [START_REF] Liou | An Introduction to Atmospheric Radiation[END_REF].

On many occasions, one wants to know the SSI I(Q) at a site Q of altitude z(Q) but has only observations of SSI I(P) at a site P of altitude z(P). Practically, this is the case of a measuring site on a mountain side (P) and a site of interest located in the valley (Q) or reciprocally. This can be extended to the case where I(P) results from spatial interpolation of neighbouring measuring sites of similar altitude. In this application, the possible exploitation of the fitting functions to all-skies case is discussed. The optical state of the atmosphere, including cloud properties, is assumed the same at P and Q. This means that if P and Q were having the same altitude, the SSI I at P and Q would be equal. It is believed that this assumption is better verified for daily or monthly means of irradiance than for shorter periods of integration. If the clear-sky value I c (P) at P is known, e.g., by a clearsky model, the clear-sky index K c is defined as:

I(P) = K c (P) I c (P) (1) Thus, the change in I due to change in altitude z is given by I / z = K c I c / z + I c K c / z
(2) LibRadtran is exploited with many different atmospheric conditions to compute each term in Eq. 2. Figure 1 exhibits examples of the changes in SSI and K c as a function of the cloud optical depth (COD) for various altitudes. These curves are quasi-independent on the solar zenith angle, visibility, precipitable water, and atmospheric profile; the COD is the major influencing factor. As expected, I and K c decrease as the COD increases. For a given COD, I increases with the altitude. This is not the case for K c ; it remains quasi-constant for altitude from 0 km to 3 km. This means that K c / z is small. The relative contribution of the rightmost term in Eq. 2: I c K c / z to I / z is now further investigated. To that purpose, derivatives are written in a discrete form: I / z  (I(z 2 ) -I(z 1 )) / (z 2 -z 1 )

(3) and I c K c / z  [(I c (z 2 ) + I c (z 1 ))/2] [(K c (z 2 ) -K c (z 1 )) / (z 2 -z 1 )] (4) where z 2 = z 1 + 1 km. I / z reaches its maximum for COD close to 3 and then decreases as the COD increases. I c K c / z reaches its maximum for COD between 4 and 6, depending on the altitude. It then decreases as the COD increases. As a whole, the values of I c K c / z are approximately 10 times lower than I / z. This figure shows that the relative contribution of I c K c / z to I / z decreases as the altitude increases and as the COD increases beyond values 4-6. In the studied cases, this relative contribution may reach approximately 20 % but most often is less than 10 %. This means that in most cases, the term I c K c / z may be neglected at first order.

Eq. 2 becomes: (12) Thus, the function used to model the vertical profile in clear-sky can be used in all-skies conditions at first order.

I / z  K c I c / z (5) or in a discrete form [I(Q) -I(P)] / [z(Q) -z(P)]  [(K c (Q)+ K c (P))/2] [(I c (Q) -I c (P)) / (z(Q) -z(P))] (6) or [I(Q) -I(P)] / [z(Q) -z(P)]  (K c (P) [(I c (Q) - I c (P)) / (z(Q) -z(P))] (7) It follows: I(Q) -I(P)  K c (P) I c (Q) -K c (P) I c (P) (8) Using Eq. 1 for I(P) I(Q)  K c (P) I c (Q) (9)

PARAMETERISATION OF VERTICAL PROFILE IN CLEAR-SKY

As a consequence of the previous conclusion, at follows that the search for a parameterisation of vertical profile of SSI can be restricted to clear-sky case. Figure 3 shows an example of vertical profiles of global I c and direct B c irradiances. Both profiles are superimposed at high altitude where the scattering effects are small and consequently, the diffuse irradiance is small. Close to sea level (z = 0), both profiles differ because of the increasing role of the scattering effects. Most of the photons removed from the direct beam end up as diffuse irradiance. Global irradiance is therefore always larger than direct irradiance. This figure clearly shows that the direct and global components decrease with z. For example, there is a decrease of 30 W m -2 in direct irradiance between 1 km and 0.5 km. We are looking for analytical functions f I and f B that model the vertical profile of global I c (z) and direct B c (z) irradiances under clear-sky at any altitude z starting from known values of I c (z 0 ) and B c (z 0 ) at an altitude z 0 :

I c (z) = I c (z 0 ) f I (z 0 , z, I c (z 0 )) (13) B c (z) = B c (z 0 ) f B (z 0 , z, B c (z 0 ))
where f I is the quantity to be determined for the global irradiance I c , f B that for the direct component B c and I 0 the irradiance at the top of the atmosphere. Our work is based on the assumption that within the range of terrain altitude on the earth: 0 to 8 km, the vertical profile of irradiance follows an exponential form. This assumption is supported by the law of Beer-Bouguer-Lambert and by the ESRA model [START_REF] Esra ; Greif | European Solar Radiation Atlas[END_REF]. Fig. 3 illustrates the soundness of this assumption. Note that this assumption is more appropriate to the direct irradiance than to the global irradiance.

Knowing two values of clear-sky irradiances I c (z 0 ) and I c (z H ) at two different altitudes z 0 and z H allows to define the parameters of a fitting function in a deterministic way. This fitting function is called "double-z fitting function" because it needs as inputs the clear-sky irradiances at two altitudes. Two attenuation coefficients, A(z), for global and, A B (z), for direct irradiances are defined:

I c (z) = I 0 (1-A(z)) (14) B c (z) = I 0 (1-A B (z))
It is assumed that A(z) and A B (z) have the following forms:

A(z) = A(z 0 ) exp[-(z -z 0 )] (15) A B (z) = A B (z 0 ) exp[- B (z -z 0 )]
Then A(z 0 ) and  are determined from Eqs 14 and 15: [START_REF] Terrainbase | Worldwide digital terrain data, Documentation Manual, CD-ROM Release 1.0[END_REF] where z H has been empirically defined, in km, as: z H = max(3, z 0 +2) ( 18) after several trials. The same equations hold for A B (z 0 ) and  B , where B c (z H ) and B c (z 0 ) are the direct clear-sky irradiances at z H and z 0 : [START_REF] Wahab | Assessing surface solar irradiance in northern africa desert climate and its long-term variations from meteosat images[END_REF] This parameterisation would allow e.g., easy implementation in spreadsheet or fast calculation for demonstration and education purposes. An alternative to the analytical parameterisation is a method based on look-up tables or abacus. The irradiances I c (z) and B c (z) are computed by invoking libRadtran for a large number of cases that are representative of all expected cases. Great attention should be brought to the description of aerosol properties as they are presently not well known and have a crucial importance on the clear sky irradiance [START_REF] Lohmann | Long-term variability of solar direct and global irradiance derived from ISCCP data and comparison with re-analysis data[END_REF][START_REF] Breitkreuz | A case study to prepare for the utilization of aerosol forecasts in solar energy industries[END_REF]. An example of these representative cases can be the following:

A(z 0 ) = 1 -(I c (z 0 ) / I 0 ) (16)  = -ln[(I 0 -I c (z H ))/(I 0 -I c (z 0 ))] / (z H -z 0 )
A B (z 0 ) = 1 -(B c (z 0 ) / I 0 ) (19)  B = -ln[(I 0 -B c (z H ))/(I 0 -B c (z 0 ))] / (z H -z 0 )
 every 1 km in altitude, from 0 km to 9 km,  every 5 • in solar zenith angle, from 0 • to 75 •,  every 10 km in visibility, from 10 km to 110 km,  four types of aerosol in the lower 2 km of atmosphere: rural, maritime, urban and tropospheric,  every 10 kg m -2 in precipitable water, from 10 kg m -2 to 80 kg m -2 ,  six atmosphere profiles: tropical, midlatitude summer, midlatitude winter, subarctic summer, subarctic winter and U.S standard,  and every 0.2 in ground albedo, from 0 to 0.8. The results are stored in look-up tables. For any set of inputs, the irradiances are computed by interpolation of the known values stored in the look-up tables. Such a method would be similarly accurate and fast.

Both solutions have benefits and drawbacks. The analytical parameterisation is easy to handle and to implement. But, observations at two different altitudes of clear-sky irradiances at the same instant and same place are very rare. The look-up tables are cumbersome to compute and then to implement in the routine operations. For efficiency, the look-up tables -whose size would be approximately 200 Mbytes-should be stored in computer memory. The accuracy of the analytical parameterisation cannot be increased, except if the parameterisation is changed. The accuracy in the look-up tables can be improved by decreasing the quantification step in all inputs at the expenses of the increase of the initial computation and of the size of the look-up tables. If new descriptions of aerosol properties are available, the analytical parameterisation adapt itself to changing conditions while the look-up tables should be recomputed for the new conditions and reimplemented in the processing software.

COMPARISON WITH TWO RTM

The outcomes of this double-z fitting function are compared for several altitudes to the outcomes of libRadtran which serves here as a reference. LibRadtran was run for several input conditions, e.g., different solar zenith angle (0, 10, 20, 40, 50, 60, 65, 70, 75, 80) deg, four different aerosol types, different visibilities [START_REF] Gschwind | Converting a successful research project into a sustainable service: the case of the SoDa Web service[END_REF][START_REF] Wahab | Assessing surface solar irradiance in northern africa desert climate and its long-term variations from meteosat images[END_REF][START_REF] Dvorkin | Modelling the altitude effect on solar UV radiation[END_REF][START_REF] Kylling | Spectral actinic flux in the lower troposphere: measurement and 1-D simulations for cloudless, broken cloud and overcast situations[END_REF]40,50,60, 80, 100) km, different amounts of precipitable water [START_REF] Rigollier | The method Heliosat-2 for deriving shortwave solar radiation from satellite images[END_REF][START_REF] Gschwind | Converting a successful research project into a sustainable service: the case of the SoDa Web service[END_REF][START_REF] Kylling | Spectral actinic flux in the lower troposphere: measurement and 1-D simulations for cloudless, broken cloud and overcast situations[END_REF]60,80) kg m -2 , different standard atmospheres and different ground albedos (0.0, 0.2, 0.3, 0.6, 0.9). Similar runs were made with 6S. Similar conclusions and values were obtained, ensuring that the comparison is not biased by the use of a single RTM.

Fig. 4 illustrates the performance in the form of the relative deviation in percentage, between the vertical profiles (from 0 km to 9 km) obtained by the fitting function and those obtained by libRadtran, in the case where z 0 is set to mean sea level. As a whole, it is observed that the relative deviation is larger for the direct irradiance than for the global irradiance. There is a dependency of the deviation with the altitude with no clear trend. Quite often, the deviation decreases when the visibility increases but not always. The visibility and solar zenith angle are the most influencing factors on the vertical profile and therefore on the performance of a fitting function to reproduce the profile. The amount of precipitable water has little effect on the deviation. Generally, the relative deviation is less than 1 % (global) or 2 % (direct) for current visibilities, i.e. greater than 30 km (WMO 1981). It may amount to respectively 3 % and 5 % for very low visibilities of 15 km. For larger visibilities, the relative deviation is almost zero for altitudes from sea level up to 6-7 km (Fig. 4). Given current visibilities and solar zenith angles less than 75 •, the relative deviations compare favourably to typical deviations on ground measurements, which are 5 % to 7 % [START_REF]Meteorological aspects of the utilization of solar radiation as an energy source. Annex: World maps of relative global radiation[END_REF].

Since in practice, discrepancies on altitude are generally less than 1 km, it is of interest to evaluate the ability to correct the irradiance assessed incorrectly at a given altitude z 0 while the actual altitude is z = z 0 +1 km. Given the irradiance I c (z 0 ) and I c (z H ), the double-z fitting function is applied to compute an estimate I* c (z). This estimate is then compared to the correct value I c (z) given by libRadtran and the difference is computed. This comparison is performed for various atmospheric conditions and various z 0 . Table 1 shows example of the deviations reached by the double-z fitting function for different z 0 . One can see that the double-z fitting function exhibit small relative deviations, often less than 1 %, for global and direct irradiances. It is concluded that this function may be used to correct satellitederived SSI.

6

-0.2 (-0.0 %) -0.2 (-0.0 %) -0.1 (-0.0 %) -0.1 (-0.0 %)

7

-0.5 (-0.1 %) -0.5 (-0.1 %) -0.2 (-0.0 %) -0.2 (-0.0 %) Table 1. Error (in W m -2 ) between actual irradiance I c (z 0 +1) provided by libRadtran and the irradiance predicted by fitting functions. The relative deviation is indicated in brackets. Two visibilities are reported: 30 km (vis 30) and 60 km (vis 60). Solar zenith angle is 60•. Other conditions are those in Fig. 3.

CONCLUSION

This study demonstrates that it is possible to find fitting functions that reproduce the vertical profile of the global and direct irradiances under clear-sky with sufficient accuracy. The possibility of using look-up-tables for the same purpose is also discussed. It is shown that in first order, the influence of altitude on all-sky irradiance is equivalent to the influence of altitude on the corresponding clear-sky irradiance. It follows that any fitting function developed for clear-sky can be used in all-sky conditions. It is also found that the double-z fitting function might be used to correct irradiance with a sufficient accuracy in the case where irradiance is incorrectly assessed because of a deviation in altitude. This result can be exploited in the case of extrapolating irradiance measuring at a station to a close site of different altitude and also in the case of design a fast method for processing satellite data based on RTM.

Figure 2

 2 Figure 2 displays the comparison between these two terms as a function of the COD.

Figure 2 .

 2 Figure 2. Comparison between I / z (full lines and left axis) and I c K c / z (dashed lines, right axis) for different cloud optical depth and altitudes (z). Other inputs are similar to those in figure 1.

  If an analytical function f I , describing the vertical profile of the SSI in clear sky I c I c (z) = I c (z 0 ) f I (z 0 , z, I c (z 0 )) (10) is available, then the irradiance I(Q) at site Q is easily computed as: I(Q) = K c (P) I c (P) f I (z(P), z(Q), I c (P)) (11) or I(Q) = I(P) f I (z(P), z(Q), I c (P))

Figure 3 .

 3 Figure 3. Example of vertical profiles of direct and global irradiances. These profiles have been computedwith the radiative transfer model libRadtran. Main input parameters were solar zenith angle (30°), precipitable water (15 kg m -2 ), visibility (60 km), number of the day in the year (170), and ground albedo (0.2).

Figure 4 .Figure 5 .

 45 Figure 4. Relative deviation obtained with the double-z fitting function. global (direct) means global (direct) irradiance, vis means "visibility", and other conditions are those in Fig.3. z 0 is equal to 0 km.