N Le Touz

J Dumoulin

J-M Piau

Solar hybrid road: from concept to modeling and lab scale mockup experiments

Keywords: Solar hybrid road, multi-physics, finite element method, thermal heat transfer, energetic system

We present in this paper the concept of solar hybrid road and the numerical model studied and develop to compute its thermal state at any time step. The aim of these new structures is to collect part of the solar radiation during the summer period, to prevent from icing at the surface during the winter period. As it is reversible, it will help also to prevent from a too high structure temperature (in particular closed to the surface) in summer season. Two types of modifications, compared to traditional structures are presented. The first one concerns the insertion of a porous layer allowing a heat fluid to exchange thermal energy. The second one concerns the use of a semitransparent material at pavement surface allowing the solar radiation to penetrate deeper in the structure. The numerical model of this multi-physics problem develop with the finite element method is presented and discussed. A validation on two test cases is proposed. First results on energy harvesting evaluation for few locations (i.e. climatic conditions) in France are presented and discussed. Finally, conclusion and perspectives are proposed.

Nomenclature

Introduction

Car traffic on roads represents at the same time economic, safety and ecologic stakes, which are particularly important during extreme environmental conditions, as heat waves or cold spells. In particular, winter conditions can entail black ice occurrence and/or snow accumulation on roads, increasing the risk of car crash and reducing the economic activities. The solution of pouring salt on roads submitted to such environmental conditions allows to maintain traffic, but have the drawback to be expensive. In fact, in North America, about 14 million tons of salt are used annually (Transportation Research Board (1991)). In France, this quantity reaches 1 million tons (Parienté (2011)). Because of the cost of salt, which is about 100 € per ton, from [START_REF] Schepman | Neige : combien coûte le régime de déneigement au sel ?, 20minutes.fr[END_REF], using salt for deicing has a non-negligible economic impact for the society, in addition to the ecological impact on the fauna and the flora, according to [START_REF] Sanzo | Effects of road de-icing salt (NaCl) on larval wood frogs (Rana sylvatica)[END_REF], due to flow of brine to rivers and ground water. This ecological impact can be evaluated to 450 M€ in Germany according to [START_REF] Huckestein | Externe Effekte des Verkehrs in Deutschland[END_REF]. Salt increases also corrosion on bridges, what entails a decrease in the lifetime of the structure [START_REF] Stewart | Time-dependant reliability of deteriorating reinforced concrete bridge decks[END_REF]).

In order to avoid these disagreements, there is also an interest to investigate alternative solutions to deice roads.

One of these solutions is to heat the surface of the road (The Hybrid Solar Road ( 2015)). Different methods exist to heat the pavement: some electrical resistance can be buried in the pavement. Oil-or gas based road heating system can also be brought into play [START_REF] Eugster | Road and Bridge Heating Using Geothermal Energy[END_REF]). However, heating roads need a large amount of energy. This high consumption can represent, for 20 square meters of pavement, the annual consumption in energy for hot water of a family [START_REF] Chen | Study of ice and snow melting process on conductive asphalt solar collector[END_REF]).

Due to amount of energy required for de-icing and cost consideration, alternative solutions have been studied and developed, in particular ones based on geothermal approaches. Such solution allows to use the energy stored in the ground during the summer when freezing conditions occur. It has moreover the advantage of needing less maintenance [START_REF] Chen | Study of ice and snow melting process on conductive asphalt solar collector[END_REF]). For the main part of these systems, the heat fluid flows in tubes buried below the pavement. We propose here a solution with a fluid flowing in a porous layer below the surface layer. To increase the ability of this hybrid road to store energy, a semi-transparent layer is used in order that the solar radiation can penetrate deeply inside the structure.

In this paper, we present a numerical model developed and studied to calculate evolution of the temperature field in the structure under known environmental conditions. It is used to establish annual energy recovery performances of these hybrid roads for various climatic conditions. In our numerical model three physical phenomena have been integrated: thermal diffusion in the whole structure, hydraulic convection in the porous layer and radiative heat transfer in the semi-transparent layer. Numerical resolution of this multi-physic problem is performed with the finite element method. We first introduce the concept of hybrid solar road, then the numerical model developed for each physic, then we detail the numerical coupling and present a comparison with measured data obtained on a laboratory mock-up. Finally, we present and discuss first results obtained on energetic performances.

Hybrid solar road concept

We focus here on the study of a multilayer system for which the two layers closest to the surface of the road are modified, compared to traditional road structures Figure 2 represents a sectional view of the studied pavement. Layer 1 is opaque or semi-transparent and waterproof. Layer 2 is porous: a fluid flows inside along the width of the road, with the effect of the imposed slope. Layer 3 and 4 are opaque and waterproof. This fluid flow allows to collect a part of solar radiation during summer period. On the contrary, under winter conditions, it allows to heat the surface of the pavement and also to prevent from icing [START_REF] Asfour | Hydrothermal modeling of porous pavement for its surface de-freezing[END_REF]).

The structure is also submitted to three coupled thermal phenomenon: heat diffusion and radiative transfer in the semi-transparent layer, heat diffusion and hydraulic convection in the porous layer and heat diffusion anywhere else. Numerical model of these three physical phenomena and their coupling are detailed in the next section.

For numerical computation, the studied structure has a length equal to 4 m, a height equal to 1.5 m. To reduce computation time, the structure is considered in 2 dimensions when conducting numerical simulations over a year. Thermal properties considered for the pavement layers are reported in table 1. Optical properties for the semi-transparent material vary with considered wavelength. For calculus, we use absorption coefficient given by [START_REF] Ping | Transient radiative-conductive heat transfer in flat glasses submitted to temperature, flux and mixed boundary conditions[END_REF] along each spectral band. These values are reported in table 2. Semi-transparent pavement is a diffusive medium; we take a scattering coefficient equal to 10 m -1 .

Multi-physic numerical model studied and developed

Our model combined three Physics using Finite Element Modelling approach and is presented hereafter.

Thermal diffusion

Thermal diffusion occurs in the whole structure. We consider here that the surface of the structure is submitted to convective exchanges with the ambient air, radiative exchanges with the atmosphere and solar radiation. Moreover, we suppose that there is no precipitation and no water nor snow at the surface. Other boundaries are supposed to be adiabatic. In particular, the depth of the considered domain is supposed to be large enough so that thermal exchange at its opposite boundary can be neglected. Heat diffusion equation to solve is reported in system (1)

� 𝜌𝜌𝐶𝐶 𝑝𝑝 𝜕𝜕𝑇𝑇 𝜕𝜕𝜕𝜕 = ∇ ⋅ [𝑘𝑘∇𝑇𝑇] + 𝑞𝑞 on Ω 𝑘𝑘∇𝑇𝑇 ⋅ 𝑛𝑛 �⃗ = Φ 𝑠𝑠 + ℎ(𝑇𝑇 𝑎𝑎 -𝑇𝑇) + 𝜀𝜀𝜎𝜎(𝑇𝑇 𝑐𝑐 4 -𝑇𝑇 4
) on the surface 𝑘𝑘∇𝑇𝑇 ⋅ 𝑛𝑛 �⃗ = 0 on the other boundaries (1) Modelization of thermal diffusion with the finite element method is well-documented, for instance in Gartling and Reddy (2010) or [START_REF] Taine | Transferts thermiques, introduction aux transferts d'énergie[END_REF].

Hydraulic convection

The porous layer, with a porosity written 𝜑𝜑, is supposed to be fully saturated in fluid. Fluid flow can be described by the Darcy law: the fluid velocity 𝑢𝑢 �⃗ is written as equation ( 2):

𝑢𝑢 �⃗ = -𝐾𝐾𝜌𝜌𝐾𝐾 𝜇𝜇 ∇ � 𝑝𝑝 𝜌𝜌𝐾𝐾 + 𝑧𝑧� (2)
Where 𝐾𝐾 is the permeability (in m 2 ), 𝐾𝐾 is the acceleration of gravity and 𝜇𝜇 the dynamic viscosity (in Pa.s). Two temperature fields are studied: a temperature field for the solid phase, noted 𝑇𝑇, defined in the whole domain Ω, and a temperature field for the fluid phase 𝑇𝑇 𝑓𝑓 , only defined on the porous domain Ω 𝑓𝑓 .

Heat equation is formulated at a microscopic scale, then the use of a representative elementary volume allows to get a coupled system of equation involving mean values, on each representative elementary volume, involving both temperature fields. The system of equations is written, from [START_REF] Bejan | Porous and Complex Flow Structures in Modern Technologies[END_REF], or [START_REF] Nield | Convection in Porous Media[END_REF]:

(1 -𝜑𝜑)𝜌𝜌𝐶𝐶 𝑝𝑝 𝜕𝜕𝑇𝑇 𝜕𝜕𝜕𝜕 = (1 -𝜑𝜑)∇ ⋅ [𝑘𝑘∇𝑇𝑇] + ℎ 𝑓𝑓𝑠𝑠 𝑎𝑎 𝑓𝑓𝑠𝑠 �𝑇𝑇 𝑓𝑓 -𝑇𝑇� (3) 𝜑𝜑𝜌𝜌 𝑓𝑓 𝐶𝐶 𝑝𝑝𝑓𝑓 𝜕𝜕𝑇𝑇 𝑓𝑓 𝜕𝜕𝜕𝜕 + 𝜌𝜌 𝑓𝑓 𝐶𝐶 𝑝𝑝𝑓𝑓 𝑢𝑢 �⃗ ⋅ ∇𝑇𝑇 𝑓𝑓 = 𝜑𝜑∇ ⋅ �𝑘𝑘 𝑓𝑓 ∇𝑇𝑇 𝑓𝑓 � + ℎ 𝑓𝑓𝑠𝑠 𝑎𝑎 𝑓𝑓𝑠𝑠 �𝑇𝑇 -𝑇𝑇 𝑓𝑓 � (4)
Where the index 𝑓𝑓 refers to the fluid phase, no index refers to the solid phase and ℎ 𝑓𝑓𝑠𝑠 is an exchange coefficient given by ( 5), from [START_REF] Nield | Convection in Porous Media[END_REF]:

1 ℎ 𝑓𝑓𝑠𝑠 = 𝑑𝑑 𝑝𝑝 𝑁𝑁𝑢𝑢 𝑓𝑓𝑠𝑠 𝑘𝑘 𝑓𝑓 + 𝑑𝑑 𝑝𝑝 𝛽𝛽𝑘𝑘 ( 5
)
𝑎𝑎 𝑓𝑓𝑠𝑠 is the specific surface (surface by volume unit), given by

𝑎𝑎 𝑓𝑓𝑠𝑠 = 6(1-𝜑𝜑) 𝑑𝑑 𝑝𝑝
, 𝑑𝑑 𝑝𝑝 is the aggregate size, 𝑁𝑁𝑢𝑢 𝑓𝑓𝑠𝑠 the Nüsselt number and 𝛽𝛽 a constant.

For the solid phase, boundary conditions are those described in (1). Moreover, for internal boundaries, the temperature field is supposed to be known. Two external boundaries are considered: with the fluid input and with the fluid output. The solid phase is supposed to be submitted to a Neumann boundary condition at the input, denoted by Γ 𝑓𝑓0 : 𝑘𝑘∇𝑇𝑇 ⋅ 𝑛𝑛 �⃗ = ℎ 𝑓𝑓 �𝑇𝑇 𝑓𝑓0 -𝑇𝑇�, with 𝑇𝑇 𝑓𝑓0 the input fluid temperature. At the output, we suppose an adiabatic boundary condition, as also used by [START_REF] Asfour | Hydrothermal modeling of porous pavement for its surface de-freezing[END_REF]. By written (𝜙𝜙 𝑖𝑖 ) a basis of 𝐿𝐿 2 �Ω 𝑓𝑓 �, temperature fields can be written 𝑇𝑇 = ∑ 𝜙𝜙 𝑖𝑖 𝑇𝑇 𝑖𝑖 𝑖𝑖 and 𝑇𝑇 𝑓𝑓 = ∑ 𝜙𝜙 𝑖𝑖 𝑇𝑇 𝑓𝑓,𝑖𝑖

𝑖𝑖

. Because of the advection term 𝜌𝜌 𝑓𝑓 𝐶𝐶 𝑝𝑝𝑓𝑓 𝑢𝑢 �⃗ ⋅ ∇𝑇𝑇, the equation ( 4) is hyperbolic, so solving numerically this equation with finite element method can entail spurious oscillations on the numeric solution (Johnson (1987). Meinköhn ( 2009)). A Petrov-Galerkin formulation, which was detailed by [START_REF] Brooks | Streamline Upwing/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes Equations[END_REF] and Johnson (1987) allows to prevent from these non-physical oscillations by adding an upwind term 𝑢𝑢 �⃗ ⋅ ∇𝜙𝜙 𝑖𝑖 to the test function 𝜙𝜙 𝑖𝑖 . The choice of the damping parameter 𝛿𝛿 is mainly due to the mesh size. Discussion about this parameter is proposed in [START_REF] Brooks | Streamline Upwing/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes Equations[END_REF], [START_REF] Codina | Comparison of some Finite Element Methods for Solving the Diffusion-Convection-Reaction Equation[END_REF], [START_REF] Donea | Finite Element Methods for Flow Problems[END_REF], and [START_REF] Shakib | A new Finite Element Formulation for Computational Fluid Dynamics: X. The Incompressible Euler and Navier-Stokes Equations[END_REF]. For any 𝜙𝜙 𝑗𝑗 , the Petrov-Galerkin formulation give us the equations ( 6) and ( 7):

� �(1 -𝜑𝜑) �� 𝜌𝜌𝐶𝐶 𝑝𝑝 𝜙𝜙 𝑖𝑖 𝜙𝜙 𝑗𝑗 dΩ Ω 𝑓𝑓 � 𝜕𝜕𝑇𝑇 𝑖𝑖 𝜕𝜕𝜕𝜕 + (1 -𝜑𝜑) �� ∇𝜙𝜙 𝑖𝑖 𝑘𝑘∇𝜙𝜙 𝑗𝑗 dΩ Ω 𝑓𝑓 + � ℎ 𝑓𝑓 𝜙𝜙 𝑖𝑖 𝜙𝜙 𝑗𝑗 dΓ Γ 𝑓𝑓0 � 𝑇𝑇 𝑖𝑖 � 𝑖𝑖 = �� ℎ 𝑓𝑓 𝜙𝜙 𝑗𝑗 𝑇𝑇 𝑓𝑓0 dΓ Γ 𝑓𝑓0 � + � ��� ℎ 𝑓𝑓𝑠𝑠 𝑎𝑎 𝑓𝑓𝑠𝑠 𝜙𝜙 𝑖𝑖 𝜙𝜙 𝑗𝑗 dΓ Γ 𝑓𝑓 � �𝑇𝑇 𝑓𝑓,𝑖𝑖 -𝑇𝑇 𝑖𝑖 �� 𝑖𝑖 , ∀𝑗𝑗 (6)
� �𝜑𝜑 �� 𝜌𝜌 𝑓𝑓 𝐶𝐶 𝑝𝑝𝑓𝑓 𝜙𝜙 𝑖𝑖 �𝜙𝜙 𝑗𝑗 + 𝛿𝛿𝑢𝑢 �⃗ ⋅ ∇𝜙𝜙 𝑗𝑗 �dΩ Ω 𝑓𝑓 � 𝜕𝜕𝑇𝑇 𝑓𝑓,𝑖𝑖 𝜕𝜕𝜕𝜕 𝑖𝑖 + �𝜑𝜑 � ∇𝜙𝜙 𝑖𝑖 𝑘𝑘∇𝜙𝜙 𝑗𝑗 dΩ + � �𝜙𝜙 𝑗𝑗 + 𝛿𝛿𝑢𝑢 �⃗ ⋅ ∇𝜙𝜙 𝑗𝑗 �𝜌𝜌 𝑓𝑓 𝐶𝐶 𝑝𝑝𝑓𝑓 𝑢𝑢 �⃗ ⋅ ∇𝜙𝜙 𝑖𝑖 dΩ Ω f Ω 𝑓𝑓 � 𝑇𝑇 𝑓𝑓,𝑖𝑖 � = � ��� ℎ 𝑓𝑓𝑠𝑠 𝑎𝑎 𝑓𝑓𝑠𝑠 𝜙𝜙 𝑖𝑖 �𝜙𝜙 𝑗𝑗 + 𝛿𝛿𝑢𝑢 �⃗ ⋅ ∇𝜙𝜙 𝑗𝑗 �dΩ Ω 𝑓𝑓 � �𝑇𝑇 𝑖𝑖 -𝑇𝑇 𝑓𝑓,𝑖𝑖 �� 𝑖𝑖 , ∀𝑗𝑗 (7)

Radiative transfer in the semi-transparent layer

We present in this part the modeling of thermal diffusion and radiative transfer coupling effects for the semitransparent layer. This one is transparent for solar radiation, and opaque for high wavelength radiation. The aim of this section is to compute the intensity of radiation in the semi-transparent layer in order to get the thermal source. Radiative transfer in the layer can be quantified by computing the radiative intensity 𝐼𝐼 𝜈𝜈 , which is governed by the radiative transfer equation introduced by [START_REF] Chandrasekhar | Radiative Transfer[END_REF]:

𝜗𝜗 ⋅ ∇𝐼𝐼 𝜈𝜈 (𝑥𝑥, 𝜗𝜗) = -(𝜅𝜅 𝜈𝜈 + 𝜎𝜎 𝜈𝜈 )𝐼𝐼 𝜈𝜈 (𝑥𝑥, 𝜗𝜗) + 𝜅𝜅 𝜈𝜈 𝐵𝐵 𝜈𝜈 �𝑇𝑇(𝑥𝑥)� + 𝜎𝜎 𝜈𝜈 4𝜋𝜋 � 𝑝𝑝 𝜈𝜈 �𝜗𝜗 ̃, 𝜗𝜗�𝐼𝐼 𝜈𝜈 �𝑥𝑥, 𝜗𝜗 ̃�d𝜗𝜗 4𝜋𝜋 (8)
Where 𝑥𝑥 denotes the position (𝑥𝑥 ∈ Ω), 𝜗𝜗 refers to the direction (𝜗𝜗 ∈ 𝑆𝑆 2 , the unit sphere). The transitional term is neglected, the evolution of phenomenon being slow [START_REF] Howell | Thermal Radiation Heat Transfer[END_REF]).

Radiative transfer equation is defined both in position along 𝑥𝑥 and in direction along 𝜗𝜗. The domain, here Ω × 𝑆𝑆 2 , is also in five dimensions. The boundaries need to be divided into two parts: input boundary Γ -and output boundary Γ + , which are defined as:

Γ + = {(𝑥𝑥, 𝜗𝜗) ∈ 𝜕𝜕Ω × 𝑆𝑆 2 such as 𝑛𝑛 �⃗ ⋅ 𝜗𝜗 ≥ 0} (9) Γ -= {(𝑥𝑥, 𝜗𝜗) ∈ 𝜕𝜕Ω × 𝑆𝑆 2 such as 𝑛𝑛 �⃗ ⋅ 𝜗𝜗 < 0} (10)
These boundaries are shown, for a given direction 𝜗𝜗, in Figure 3.

Some methods in literature allows to solve the radiative transfer equation, as finite volume, Monte Carlo, spherical harmonics [START_REF] Howell | Thermal Radiation Heat Transfer[END_REF], Modest ( 2003)). We propose here to focus on the finite element method presented by [START_REF] Kanschat | Solution of Radiative Transfer Problems with Finite Elements[END_REF] and [START_REF] Richling | Radiative Transfer with Finite Element. i. Basic Method and Tests[END_REF] for radiative transfer applied in astrophysics. The finite element method is here transposed and adapted for the case of a semi-transparent solid media with small dimensions. Two discretization are required for directional and spatial dependences. Directional discretization is performed with the discrete ordinate method. For a directional mesh, based on the unit sphere, with 𝑁𝑁 𝑆𝑆 nodes 𝜗𝜗 𝑖𝑖 located at the center of elements of weight 𝑤𝑤 𝑖𝑖 , integral terms on directions can be approximated, for any variable 𝑢𝑢 continuous along the unit sphere, by the expression in equation ( 12):

� 𝑢𝑢(𝜗𝜗) 4𝜋𝜋 𝑑𝑑𝜗𝜗 = � 𝑤𝑤 𝑖𝑖 𝑢𝑢(𝜗𝜗 𝑖𝑖 ) 𝑁𝑁 𝑆𝑆 𝑖𝑖=1 ( 12
)
This discretization is equivalent to a discontinuous Galerkin method with constant interpolation functions on the elements. In order to get a mesh with the most uniform allocation of direction nodes, we propose to use a mesh based on a regular polyhedron with triangular faces, as the icosahedron suggested by [START_REF] Richling | Radiative Transfer with Finite Element. i. Basic Method and Tests[END_REF]. Each direction of the mesh linked the center of the unit polyhedral to the centroid of a face. Weights being equal to the projected area of the face on the unit sphere, all weights are equal to

4𝜋𝜋 𝑁𝑁 𝑆𝑆
. The number of directions can be increased recursively by divided each triangular face into four similar triangles, this number of directions being also multiplied by four at each iteration.

Applying the discrete ordinate method to the equations ( 11) and ( 12) allows to obtain the system (13) and the boundary conditions ( 14). ( 13) and ( 14) are systems of 𝑁𝑁 𝑆𝑆 equations, defined only on space:

𝜗𝜗 𝑘𝑘 ⋅ ∇𝐼𝐼 𝜈𝜈 (𝑥𝑥, 𝜗𝜗 𝑘𝑘 ) = -(𝜅𝜅 𝜈𝜈 + 𝜎𝜎 𝜈𝜈 )𝐼𝐼 𝜈𝜈 (𝑥𝑥, 𝜗𝜗 𝑘𝑘 ) + 𝜅𝜅 𝜈𝜈 𝐵𝐵 𝜈𝜈 �𝑇𝑇(𝑥𝑥)� + 𝜎𝜎 𝜈𝜈 4𝜋𝜋 � 𝑤𝑤 𝑖𝑖 𝑝𝑝 𝜈𝜈 (𝜗𝜗 𝑙𝑙 , 𝜗𝜗 𝑘𝑘 )𝐼𝐼 𝜈𝜈 (𝑥𝑥, 𝜗𝜗 𝑙𝑙 ) 𝑁𝑁 𝑆𝑆 𝑙𝑙=1 , ∀(𝑥𝑥, 𝑘𝑘) ∈ Ω × [|1, 𝑁𝑁 𝑆𝑆 |] (13) 𝐼𝐼 𝜈𝜈 (𝑥𝑥, 𝜗𝜗 𝑘𝑘 ) = 𝐾𝐾 𝑜𝑜𝑜𝑜𝑜𝑜 (𝑥𝑥, 𝜗𝜗 𝑘𝑘 ) + 𝐾𝐾 𝑖𝑖𝑖𝑖 (𝑥𝑥, 𝜗𝜗 𝑘𝑘 , 𝐼𝐼 𝜈𝜈 ) , ∀(𝑥𝑥, 𝜗𝜗 𝑘𝑘 ) ∈ Γ - (14)
To apply the finite element method on this coupled system of equations, we introduce the following notations:

• 𝑊𝑊 Ω = {𝑢𝑢 ∈ 𝐿𝐿 2 (Ω) such as ∀𝜗𝜗 ∈ 𝑆𝑆 2 , 𝜗𝜗 ⋅ ∇𝑢𝑢 ∈ 𝐿𝐿 2 (Ω × 𝑆𝑆 2 )} • 𝜕𝜕Ω 𝑘𝑘 -= {𝑥𝑥 ∈ 𝜕𝜕Ω such as 𝑛𝑛 �⃗(𝑥𝑥) ⋅ 𝜗𝜗 𝑘𝑘 < 0} • 𝑊𝑊 0,𝑘𝑘 = {𝑢𝑢 ∈ 𝐿𝐿 2 (𝜕𝜕Ω 𝑘𝑘 -)}
We define inner products 〈⋅,⋅〉 Ω on Ω and 〈⋅,⋅〉 𝑘𝑘 -on 𝜕𝜕Ω 𝑘𝑘 -such as:

〈𝑢𝑢, 𝑣𝑣〉 Ω = � 𝑢𝑢(𝑥𝑥)𝑣𝑣(𝑥𝑥)𝑑𝑑Ω Ω , 𝑢𝑢, 𝑣𝑣 ∈ 𝑊𝑊 Ω (15) 〈𝑢𝑢, 𝑣𝑣〉 𝑘𝑘 -= � 𝑢𝑢(𝑥𝑥)𝑣𝑣(𝑥𝑥)𝑑𝑑𝑆𝑆 𝜕𝜕Ω 𝑘𝑘 - , 𝑢𝑢, 𝑣𝑣 ∈ 𝑊𝑊 0,𝑘𝑘 (16)
The systems of equations ( 13) and ( 14) can also be rewritten under weak form with these notations, under a basis (𝜑𝜑 𝑖𝑖 ) of 𝐿𝐿 2 (Ω): 〈𝒜𝒜 𝜈𝜈 𝐼𝐼 𝜈𝜈 (⋅, 𝜗𝜗 𝑘𝑘 ), 𝜑𝜑 𝑖𝑖 〉 Ω = 〈𝑓𝑓 𝑘𝑘 , 𝜑𝜑 𝑖𝑖 〉 Ω (17) 〈ℬ 𝜈𝜈 𝐼𝐼 𝜈𝜈 (⋅, 𝜗𝜗 𝑘𝑘 ), 𝜑𝜑 𝑖𝑖 〉 𝑘𝑘 -= 0 (18) With 𝒜𝒜 𝜈𝜈 and ℬ 𝜈𝜈 two linear operators. As for the hydraulic convection equation ( 4), radiative transfer equations discretized in directions (13) are hyperbolic if there is scattering in the media. To avoid spurious oscillations, a regularization needs to be applied. We overcome this problem with a Petrov-Galerkin formulation. Moreover, we inject boundary conditions (18) as weak constraints in the system (17). The systems ( 17) and ( 18) are then rewritten as:

〈𝒜𝒜 𝜈𝜈 𝐼𝐼 𝜈𝜈 (⋅, 𝜗𝜗 𝑘𝑘 ), 𝜑𝜑 𝑖𝑖 + 𝛿𝛿𝜗𝜗 𝑘𝑘 ⋅ ∇𝜑𝜑 𝑖𝑖 〉 Ω + 〈ℬ 𝜈𝜈 𝐼𝐼 𝜈𝜈 (⋅, 𝜗𝜗 𝑘𝑘 ), 𝜑𝜑 𝑖𝑖 〉 𝑘𝑘 -= 〈𝑓𝑓 𝑘𝑘 , 𝜑𝜑 𝑖𝑖 〉 Ω , ∀𝑘𝑘 ∈ [|1, 𝑁𝑁 𝑆𝑆 |] (19)
For a spatial mesh with 𝑁𝑁 Ω nodes, we get a linear system with 𝑁𝑁 Ω × 𝑁𝑁 𝑆𝑆 unknowns for each frequency interval. Absorption of radiation and heat production in the semi-transparent media occur by two ways. In the material itself, absorption of radiation entail a volumic heat source denoted here by 𝑆𝑆 𝑣𝑣𝑜𝑜𝑙𝑙 , and at boundaries between semitransparent and opaque media, there is also absorption of radiation which can be written under the form of a heat flux 𝑆𝑆 𝑓𝑓𝑙𝑙𝑜𝑜𝑓𝑓 . Both two components can be written as equations ( 20) and ( 21) as suggest by from Modest (2003).

𝑆𝑆 𝑣𝑣𝑜𝑜𝑙𝑙 = -∇ ⋅ �� � 𝐼𝐼 𝜈𝜈 (𝑥𝑥, 𝜗𝜗)𝜗𝜗𝑑𝑑𝜗𝜗𝑑𝑑𝜗𝜗 4𝜋𝜋 +∞ 𝜈𝜈=0 � (20)
𝑆𝑆 𝑓𝑓𝑙𝑙𝑜𝑜𝑓𝑓 = � �(1 -𝜌𝜌 𝜈𝜈 ) � 𝐼𝐼 𝜈𝜈 �𝑥𝑥, 𝜗𝜗 ̃� 𝜗𝜗⋅𝑖𝑖 �⃗>0 𝜗𝜗 ̃⋅ 𝑛𝑛 �⃗𝑑𝑑𝜗𝜗 -𝜖𝜖 𝜈𝜈 𝜋𝜋𝐵𝐵 𝜈𝜈 �𝑇𝑇(𝑥𝑥)�� 𝑑𝑑𝜗𝜗 +∞ 𝜈𝜈=0 (21)

Solving the coupled problems : Thermal diffusion, Hydraulic convection and Radiative transfer

Temperature field in the structure is computed dynamically from environmental data (sky temperature, air temperature and solar radiation). At each temporal iteration, radiative transfer equation is first solved (system (19) discretized in space), allowing to get the intensity field, then source terms are computed with ( 20) and ( 21).

The heat sources terms are then injected in the heat diffusion equation ( 1), which is coupled with the hydraulic convection equations ( 6) and ( 7). A Crank-Nicholson scheme is used for the temporal resolution.

Validation of the numerical model

Two test cases are presented hereafter.

Solar road with opaque surface in established steady state

We focus in this section on the case of an opaque surface layer, for an established steady state. The aim is here to get an estimation of the minimal input temperature of the fluid allowing to maintain the surface temperature higher than 0°C. Two methods are used: the numerical model presented in previous sections and an analytical formulation. This formulation is first described. The studied road structure is presented in figure 5 (left). We consider in this paragraph a one dimension problem: temperature in layers depends only from the distance to the injection point. Three layers are supposed to constitute the road structure: the porous layer is located between two impermeable layers. We suppose moreover that temperature in the middle of layers is equal to the mean temperature between the extremities of the layer and that there is a local thermal equilibrium. Under these assumptions, it can be shown that the temperature of the fluid can be written as an explicit function of the distance 𝑠𝑠 to the injection point:

𝑇𝑇 𝑓𝑓 (𝑠𝑠) = 𝑇𝑇 𝑓𝑓 (0)𝑒𝑒 -𝑠𝑠 𝜅𝜅 + 𝑇𝑇 ∞ �1 -𝑒𝑒 -𝑠𝑠 𝜅𝜅 � (22)
With : The temperature of the surface layer 𝑇𝑇 𝑠𝑠 can then be written as:

• 𝑇𝑇 ∞ = ℎ 𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐
𝑇𝑇 𝑠𝑠 (𝑠𝑠) = 𝑟𝑟 𝑑𝑑𝑟𝑟𝑠𝑠 𝑇𝑇 𝑓𝑓 (𝑠𝑠) + ℎ 𝑐𝑐𝑜𝑜𝑖𝑖𝑣𝑣 𝑇𝑇 𝑎𝑎𝑖𝑖𝑟𝑟 + ℎ 𝑟𝑟𝑎𝑎𝑑𝑑 𝑇𝑇 𝑠𝑠𝑘𝑘𝑠𝑠 + Φ 𝑟𝑟 𝑝𝑝𝑠𝑠 + ℎ 𝑐𝑐𝑜𝑜𝑖𝑖𝑣𝑣 + ℎ 𝑟𝑟𝑎𝑎𝑑𝑑 (23)
By denoting 𝐿𝐿 the width of the road, the input temperature of the fluid allowing to maintain a surface temperature higher than 0°C can be calculated:

𝑇𝑇 𝑓𝑓,𝑚𝑚𝑖𝑖𝑖𝑖 (0) = -𝑒𝑒 𝐿𝐿 𝜅𝜅 � 1 𝑟𝑟 𝑝𝑝𝑠𝑠 �ℎ 𝑐𝑐𝑜𝑜𝑖𝑖𝑣𝑣 𝑇𝑇 𝑎𝑎𝑖𝑖𝑟𝑟 + ℎ 𝑟𝑟𝑎𝑎𝑑𝑑 𝑇𝑇 𝑠𝑠𝑘𝑘𝑠𝑠 + Φ� + 𝑇𝑇 ∞ �1 -𝑒𝑒 -𝐿𝐿 𝜅𝜅 �� (24)
In a first approach, we consider here that solar road structure is studied under constant environmental conditions, what can be performed with a climatic chamber. We suppose although that the sky temperature is equal to the air temperature and that there is no solar radiation. Under these assumptions, equation ( 24) can be rewritten:

𝑇𝑇 𝑓𝑓,𝑚𝑚𝑖𝑖𝑖𝑖 (0) = -𝑒𝑒 𝐿𝐿 𝜅𝜅 � ℎ 𝑐𝑐𝑜𝑜𝑖𝑖𝑣𝑣 + ℎ 𝑟𝑟𝑎𝑎𝑑𝑑 𝑟𝑟 𝑝𝑝𝑠𝑠 + 1 -𝑒𝑒 -𝐿𝐿 𝜅𝜅 � 𝑇𝑇 𝑎𝑎𝑖𝑖𝑟𝑟 (25)
A Darcy velocity value equal to 2.18 × 10 -4 m.s -1 , corresponding to a slope of 1% and a permeability of 2.18 × 10 -2 m.s -1 (Asfour ( 2016)), is chosen. We consider the structure presented in figure 5 and the thermal parameters of table 1. Two wind speed are studied: 𝑉𝑉 𝑤𝑤 = 0 m.s -1 and 𝑉𝑉 𝑤𝑤 = 3.1 m.s -1 , with a convective exchange coefficient given by ℎ 𝑐𝑐𝑜𝑜𝑖𝑖𝑣𝑣 = 5.6 + 4𝑉𝑉 𝑤𝑤 , also used by [START_REF] Asfour | Hydrothermal modeling of porous pavement for its surface de-freezing[END_REF]). Results for input fluid temperature both with (25) and numerical simulation are given in table 3.

Table 3. Minimal temperature (in °C) for input fluid allowing to maintain a positive temperature along the surface The comparison between analytic and numerical formulation is shown in figure 5 (right). The input temperature with the analytic formulation is close to the computed solution, with a difference of about 10%. It can be noticed that the analytical formulation always under evaluates the input temperature compared to the computed temperature. Such difference may be explain by the thermal diffusion in the second dimension that is not considered. There is no temperature variation in surface layer along the depth, so the surface temperature corresponds more to the temperature at the middle of the surface layer rather than the temperature at its the surface, which is lower.

Wind speed Method 𝑇𝑇 𝑎𝑎𝑖𝑖𝑟𝑟 = -10°C 𝑇𝑇 𝑎𝑎𝑖𝑖𝑟𝑟 = -8°C 𝑇𝑇 𝑎𝑎𝑖𝑖𝑟𝑟 = -6°C 𝑇𝑇 𝑎𝑎𝑖𝑖𝑟𝑟 = -4°C 𝑇𝑇 𝑎𝑎𝑖𝑖𝑟𝑟 = -2°C 0

Solar road with semi-transparent surface: comparison with data acquired on a laboratory mock-up

In order to study the validity of our approach for the semi-transparent surface layer, a mock-up has been built (figure 6). A 1500 W halogen lamp at 3000 K allows to simulate the solar heat flux and the input fluid temperature is measured with a thermocouple. A FLIR A65 infrared camera allows to measure the surface temperature of the structure. A meteorological station and a pyranometer measure the relative humidity, ambient temperature and the solar radiation flux density. The halogen lamp is switched on during some hours, allowing to get a permanent established regime. Knowledge of environmental conditions enables us to apply the numerical model presented in section 2. A comparison between numerical and experimental data can then be led. The difference of surface temperature between the numerical model and measurements from the infrared camera is presented in figure 7 (right). On a central area, the difference of temperature between experiments and numerical simulations is about 1 °C. The numerical model drives to results close to the measures. Photovoltaic panels being present in the mock-up can explain a part of the observed temperature difference.

Energy harvesting performance evaluation: preliminary results and analysis

In this section, we determine the energetic performance of solar hybrid roads for energy harvesting over a year and for few climatic conditions in France. Annual weather data are used for eight locations in France, corresponding to the main climatic conditions. The schematic sectional view of the structure was shown in figure 2. The fluid is injected in the porous layer at a constant temperature, equal to 15°C in the presented numerical simulations. Environmental conditions entail temporal variations of the temperature field in the structure and also for the temperature of the fluid at the output of the structure. The velocity of the fluid flow is supposed to be maintained constant, equal to 1 × 10 -4 m.s -1 . Thermal and optical properties are supposed to be independent from temperature and from mechanical evolution of the structure over time.

We consider here only the ability of the structure to recover energy, that is to say that we consider only the moments when the output temperature of the fluid is higher than the input temperature. This variation of temperature between the input and the output is supposed to allow energy recovering when positive. Simulations are performed on more than one year of data, so that initial conditions are erased. Since output temperature during winter period is most of the time lower than the input temperature, only summer period is considered for energy recovering, from May to October. Efficiency of this modified structure can be evaluated by calculating the ratio between the energy harvested by the system and the cumulated solar radiation during the studied period. Results obtained are shown in table 3. This table also shows that the structure with a semi-transparent pavement allows to collect more energy than the one with the opaque pavement. For opaque surface wearing course, the gain is between 21 and 29% of the recoverable energy. Efficiency values presented in table 3 show that the gain of efficiency by using a semitransparent pavement surface layer is between 6.5 and 7% of the incoming solar radiation for all these locations. These data of recoverable energy also depend highly from the width of the road. In fact, the efficiency of this heat exchanger decrease as the length of the flow increases. This deterioration of the performances of the structure with its width is shown in figure 8 (right). This figure represents the ratio of evolution of the energy recovery term 𝐸𝐸 computed from equation ( 22) for some lengths from the input area on the energy recovery at the output of the structure, after 4 meters. This figure points out the fact that half of the contribution for energy recovering occurs less than 1.4 meter after the injection of the heat fluid. On the same way, 90% of the contribution occurs after about 3 meters, at 75% of the width of the structure. For dimensioning the structure, it is also important to take into account this decrease in performances, the same length yielding more energy if located near the injection area. There are however no significant differences between both pavements.

Conclusion

In this study, we have presented the concept of solar hybrid road. In order to get an estimation of energetic performance of such structures, multi-physics numerical model has been studied and developed, implying thermal diffusion, hydraulic convection and radiative transfer. Each physic has been described and finite element method has been applied on each of them, with variations of the Galerkin formulation for hydraulic convection and radiative transfer equations. A numerical coupling has been then explained in order to realize numerical simulation relying on environmental data. Our numerical model has been compared with data coming from two test cases: first, a comparison with an analytical expression in steady state has been performed with an opaque surface layer. Finally, numerical data from the model are compared with data acquired on a mock-up with a semi-transparent surface layer submitted to a radiative heat flux. These two test cases show that the numerical model is coherent with particular test case allowing analytical solution or experimental results obtained in laboratory conditions. First results obtained for energy harvesting performances have been presented and discussed.

Future works will refine energy harvesting performance evaluation and also address heating efficiency (during winter season).

Figure 1 :

 1 Figure 1: Schematic view of a modified structure from a 3D CAD Figure 2: Sectional schematic view of modified road structure

Figure 3 :

 3 Figure 3: Sub-boundaries for the studied domain Boundary conditions are required only in the input boundaries, due to the structure of the radiative transfer equation. Radiative intensity at the input boundaries comes from external sources such as solar radiation and internal contributions due to reflections: 𝐼𝐼 𝜈𝜈 (𝑥𝑥, 𝜗𝜗) = 𝐾𝐾 𝑜𝑜𝑜𝑜𝑜𝑜 (𝑥𝑥, 𝜗𝜗) + 𝐾𝐾 𝑖𝑖𝑖𝑖 (𝑥𝑥, 𝜗𝜗, 𝐼𝐼 𝜈𝜈 ), for (𝑥𝑥, 𝜗𝜗) ∈ Γ -(11)

Figure 5 :

 5 Figure 5: Studied structure for established steady state (left image) -Minimal fluid temperature to prevent surface from icing as a function of air temperature (right graph)

  𝑇𝑇 𝑎𝑎𝑎𝑎𝑎𝑎 +𝑟𝑟 𝑎𝑎𝑎𝑎𝑟𝑟 𝑇𝑇 𝑠𝑠𝑘𝑘𝑠𝑠 +Φ ℎ 𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐 +ℎ 𝑎𝑎𝑎𝑎𝑟𝑟 , • 𝜅𝜅 = 𝜌𝜌 𝑓𝑓 𝑐𝑐 𝑓𝑓 ℎ 𝑟𝑟𝑎𝑎 𝑜𝑜�𝑟𝑟 𝑝𝑝𝑠𝑠 +ℎ 𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐 +ℎ 𝑎𝑎𝑎𝑎𝑟𝑟 � 𝑟𝑟 𝑝𝑝𝑠𝑠 (ℎ 𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐 +ℎ 𝑎𝑎𝑎𝑎𝑟𝑟 ) , • 𝑟𝑟 𝑝𝑝𝑠𝑠 = 2𝑘𝑘 𝑠𝑠 𝑘𝑘 𝑝𝑝 ℎ 𝑠𝑠 𝑘𝑘 𝑝𝑝 +ℎ 𝑝𝑝 𝑘𝑘 𝑠𝑠 . • Φ refers to the solar radiation, • 𝑇𝑇 𝑎𝑎𝑖𝑖𝑟𝑟 to the air temperature, • 𝑇𝑇 𝑠𝑠𝑘𝑘𝑠𝑠 to the sky temperature, • ℎ 𝑠𝑠 to the thickness of the surface layer, • ℎ 𝑝𝑝 to the thickness of the porous layer, • 𝑘𝑘 𝑠𝑠 to the thermal conductivity of the surface layer, • 𝑘𝑘 𝑝𝑝 to the thermal conductivity of the porous layer, • 𝑢𝑢 to the Darcy velocity of the fluid, • ℎ 𝑐𝑐𝑜𝑜𝑖𝑖𝑣𝑣 to the convective exchange coefficient • ℎ 𝑟𝑟𝑎𝑎𝑑𝑑 to the radiative exchange coefficient.

Figure 6 :

 6 Figure 6: Experimental set-up An example of thermal image obtained with the infrared camera is presented in figure 7.

Figure 7 :

 7 Figure 7: Example of infrared image captured by the infrared camera (left image) -Comparison of measured temperature with the numerical simulation (middle image with zoom on right)

  Energy recovering capacity is defined as: 𝐸𝐸 = �𝜌𝜌 𝑓𝑓 𝐶𝐶 𝑝𝑝,𝑓𝑓 ℎ 𝑃𝑃𝑜𝑜 𝑢𝑢(𝑇𝑇 𝑠𝑠𝑜𝑜𝑟𝑟𝑜𝑜𝑖𝑖𝑠𝑠 -𝑇𝑇 𝑠𝑠𝑖𝑖𝑜𝑜𝑟𝑟𝑠𝑠𝑠𝑠 ) + 𝑑𝑑𝜕𝜕 𝑜𝑜 (22) With (⋅) + referring to the positive part and ℎ 𝑃𝑃𝑜𝑜 to the thickness of the porous layer. The eight locations studied are submitted to various solar radiation conditions during this period, with cumulated solar contributions between 4047 kWh.m -2 in Nancy and 5524 kWh.m -2 in Carpentras, as shown in figure 8 (left).

Figure 8 :

 8 Figure 8: Cumulated solar radiation (May-October period) on left -Evolution of the part of final energy recovery along the width of the structure for opaque and semi-transparent pavements (right)

Table 1

 1

	. Thermal properties of the layers of the structure and fluid, from Asfour et al. (2016), Incropera
	and DeWitt (1990) and Taine et al. (2014)
		Layer 1	Layer 2	Layer 3	Layer 4	Fluid
	𝑘𝑘 [W.m -1 .K -1 ]	0.85	1.03	1.40	1.10	0.60
	𝜌𝜌 [kg.m -3 ]	2700	2360	2620	2300	1000
	𝐶𝐶 𝑝𝑝 [J.kg -1 .K -1 ]	840	780	860	1000	4180

Table 2

 2

	. Optical properties of the semi-transparent media
	𝜆𝜆 [µm]	𝜆𝜆 < 0.5	0.5 ≤ 𝜆𝜆 < 2.7 2.7 ≤ 𝜆𝜆 < 4.5	4.5 ≤ 𝜆𝜆 < 50	𝜆𝜆 ≥ 50
	𝜅𝜅 𝜆𝜆 [m -1 ]	0	10	1000	5000	+∞
	𝜎𝜎 𝜆𝜆 [m -1 ]	10	10	10	10	10

Table 3 .

 3 Results for 4m width opaque and semi-transparent pavement

	Location	Solar cumulated	Recoverable	Efficiency	Recoverable	Efficiency	Energy gain
		input (Op	energy (Op	(Op	energy (ST	(ST	[kWh/m 2 ]
		pavement)[kWh/m 2 ]	pavement)	pavement )	pavement)	pavement)
			[kWh/m 2 ]	[%]	[kWh/m 2 ]	[%]
	Agen	4615	1440	31,2	1744	37,7	+304 (+21%)
	Carpentras	5524	1812	32,8	2188	39,6	+376 (+21%)
	La Rochelle	4973	1288	25,9	1657	33,2	+369 (+29%)
	Macon	4164	1299	31,2	1581	37,8	+282 (+22%)
	Nancy	4047	1117	27,6	1394	34,4	+277 (+25%)
	Nice	5437	1805	33,2	2201	40,4	+396 (+22%)
	Rennes	4336	1084	25,0	1386	31,9	+302 (+28%)
	Trappes	4312	1190	27,6	1470	34,1	+280 (+24%)

Acknowledgements

French Ministry of "Transition écologique et solidaire" for supporting part of this work under grant agreement DGITM N° 17/389. Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org)


Loïc Louison
email: loic.louison@guyane.univ-ag.fr

Abdennebi Omrane
email: aomrane@guyane.univ-ag.fr

Harry Ozier-Lafontaine
email: harry.ozier-lafontaine@antilles.inra.fr

Delphine Picart
email: dpicart@bordeaux.inra.fr

Modeling plant nutrient uptake : mathematical analysis and optimal control

Keywords: 2010 Mathematics Subject Classification. Primary: 35M32, 35Q92, 49J20, 49K20; Secondary: 49M15, 65M06 NTB-System, transport-diffusion system, nutrient uptake, optimal control, numerical analysis

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1. Introduction. Plant nutrition has always been a main challenge in crop production but chemical fertilizers can be no longer the only sources of nutrients for crop growth and development. Switching from conventional intensive cropping systems to ecologically ones is necessary (see Griffon [START_REF] Griffon | Nourrir la Planète[END_REF]). The awareness of the dependency to fossil fuels for the production of chemical inputs and the enactment of new rules for a healthier environment, have led to other agroecological alternatives in fertilizing practices. Plant growth is strongly linked to the amount of soil nutrients absorbed from its roots. These nutrients are produced naturally, they are present in the groundsoil at various levels of concentration. They may also be provided by humans in the form of chemical fertilizers, or by a secondary companion crop plants. It is the case for the nitrogen fixation and transfer from nitrogen fixing crops to cash crops in mixed cropping systems (see Jalonen et al. [START_REF] Jalonen | Root exudates of a legume tree as a nitrogen source for a tropical fodder grass[END_REF][9], Daudin and Sierra [START_REF] Daudin | Spatial and temporal variation of below-ground N transfer from a leguminous tree to an associated grass in an agroforestery system[END_REF]).

Modeling is required to understand the dynamics and the nutrient transfer mechanism and to define and test scenarii for the optimization of crop nutrition and sustainable systems. The root nutrient uptake as well as in solute movement in the soil are well explained by Tinker and Nye in [START_REF] Tinker | Solute Movement in the Rhizosphere[END_REF], where they describe the uptake and nutrient motion processes from a biological and chemical point of view, and where they give the model of nutrient uptake and transfer using partial differential equations (PDE) known as the Nye-Tinker-Barber (NTB) system. Plant nutrition models that have emerged in the 60s with, among others, the works by Nye [START_REF] Nye | The effect of the nutrient intensity and buffering power of a soil, and the absorbing power, size and root hairs of a root, on nutrient absorption by diffusion[END_REF] and by Nye and Marriott [START_REF] Nye | A theoretical study of the distribution of substances around roots resulting from simultaneous diffusion and mass flow[END_REF]. See also the work by Itoh and Barber [START_REF] Itoh | A numerical solution of whole plant nutrient uptake for soil-root systems with root hairs[END_REF], or by Cushman [START_REF] Cushman | Nutrient transport inside and outside the root rhizosphere: Generalized model[END_REF] who suggests a general framework of the model of nutrient uptake by roots in which he adds a term called source or sink that models either the increase or the decrease in solute concentration w.r.t time and space.

More recently, Roose [START_REF] Roose | Mathematical Model of Plant Nutrient Uptake[END_REF][20] used the NTB model in order to reflect in a more accurate way the morphology of the root system (modeling of root growth, root hair, mycorrhizae, ..) and the spatio-temporal dynamics of the solute in the soil. We also mention a recent work by Ptashnik [START_REF] Ptashnyk | Derivation of a macroscopic model for nutrient uptake by hairy-roots[END_REF], where the author studied a process of nutrient uptake by a single root branch using the asymptotic expansion method or a similar work by Schnepf et al. [START_REF] Schnepf | Impact of growth and foraging strategies of arbuscular mycorrhizal fungi on plant phosphorus uptake[END_REF].

In the above mentioned works, no attempt was made concerning the question of optimal control for crop models. In this article we focus on the modeling of plant nutrient uptake by roots from multiple nutrient sources, and more precisely, on the determination of the optimal amount of nutrients that must be provided in the soil so that the plant grows in the best conditions to satisfy its demand.

This article is built as follows: in Section 2 we present the model of nutrient uptake and we provide a proof of existence and uniqueness of its solution. In Section 3, we study the optimal control question and we give a characterization of the optimal solution. Section 4 is devoted to some simulations where we compute the optimal nutrient input for a given plant's growth in time.

Statement of the problem.

Let Ω ∈ R 3 be the part of the soil close to the root called the rhizosphere, of regular boundary Γ = Γ 1 ∪ Γ 2 , such that Γ 1 ∩ Γ 2 = ∅ as shown in Figure 1. Here, Γ 1 represents the root surface and Γ 2 plays the role of the rhizosphere frontier to the rest of the soil. During a time t ∈ (0, T ), the nutrient transport-diffusion and uptake by the roots is here described by the following Nye-Tinker-Barber (NTB) system:

                 α ∂c ∂t + q • ∇c -D∆c = 0 in Q :=]0, T [×Ω, div q = 0 in Q, (D∇c -1 2 qc) • n = Ic K on Σ 1 :=]0, T [×Γ 1 , (D∇c -1 2 qc) • n = 0 on Σ 2 :=]0, T [×Γ 2 , c(0, x) = c 0 (x) in Ω, (1)
where c = c(t, x) is the concentration of nutrient density at time t with position x. The coefficient α = b + θ is a constant, where b is the buffer power and θ the liquide saturation, the vector q = q(t, x) is the Darcy flux, and D the diffusion coefficient, which is a positive constant. The function h(c) = Ic/K is the Michaelis-Menten function, and represents the inflow density to the root surface, where I is the maximum uptake constant, and K is the Michaelis-Menten constant. 2.1. Existence of a solution to the NTB problem. This part is devoted to the weak formulation to problem (1) and to the proof of the existence of a positive unique solution c(t, x). We will consider the function t → c(t; .) ∈ V such that V ⊂ H 1 (Ω). More precisely, consider the Hilbert space:

V = ψ ∈ H 1 (Ω); ψ | Γ 2 = 0 ,
endowed with the H 1 (Ω) norm and equivalently with its semi-norm ∇ψ L 2 (Ω) . We also consider L 2 (Ω) as the pivot space i.e., V ⊂ L 2 (Ω) ⊂ V .

Lemma 2.1. The problem (1) has the equivalent weak formulation given by the following. Given

c 0 ∈ V , find c : t ∈ [0, T ] → c(t) ∈ V such that :    α d dt Ω c(t) ψ dx + a(t; c, ψ) = 0 a.e. t ∈]0, T [, ∀ ψ ∈ V, c(0, x) = c 0 (x), (2)
where

a(t; c, ψ) = 1 2 Ω q • (ψ∇c -c∇ψ) dx + D Ω ∇c ∇ψ dx - Γ1 Ic K ψ dx. (3)
Remark 1. We suppose that the concentration c(t; .) belongs the set of continuous

functions in time c ∈ L 2 (]0, T [; V ) ∩ C(]0, T [; L 2 (Ω)). The function c(t) is considered as an element of V which represents the function x → c(t, x) for all t ∈ [0, T ].
We point out that the regularity hypothesis on time is natural. The weak formulation (2)-( 3) is easy to get, and we refer to Allaire [START_REF] Allaire | Numerical Analysis and Optimization[END_REF] or to Louison [START_REF] Louison | Analysis and Optimal Control of Transport-Diffusion Problems of Incomplete Data: Agroecology Application to Nutrient Uptake in Polluted Soils[END_REF] for details.

Proposition 1 (Existence). We suppose that the flux |q| is uniformly bounded q ∈ L ∞ (Q). Then, there is a unique solution c ∈ V to the problem (2)-(3).

Proof. We show first that the bilinear form a is continuous. We have by the Cauchy-Schwartz inequality:

|a(t; c, ψ)| ≤ (q ∞ + D) c V ψ V + c L 2 (Γ1) ψ L 2 (Γ1)
where

q ∞ = q L ∞ (Q) . Now since ϕ | Γ 2 = 0, for ϕ = c, ψ, we have ϕ 2 L 2 (Γ1) = 2 Ω ϕ ∇ϕ dx ≤ 2 ϕ L 2 (Ω) ∇ϕ L 2 (Ω) ≤ 2 ϕ 2 V .
Hence, there is a positive constant

C > 0 such that |a(t; c, ψ)| ≤ C c V ψ V with C = q ∞ + D + 2,
which implies that a is continuous. Now we show that a is coercive in the sense of evolution problems. We have:

a(t; c, c) = D Ω |∇c| 2 dx - Γ1 I K |c| 2 dx ≥ D ∇c 2 L 2 (Ω) -c 2 L 2 (Γ1) . (4)
From the Young inequality we obtain

c 2 L 2 (Γ1) ≤ γ c 2 L 2 (Ω) + 1 γ ∇c 2 L 2 (Ω)
. Hence, the bilinear form a satisfies the conditions of J.-L. Lions's Theorem (Cf. Lions-Magenes [START_REF] Lions | Problèmes Aux Limites non Homogènes et Applications[END_REF]), and we have Gårding's inequality:

a(t; c, c) ≥ D - 1 γ c 2 V -γ c 2 L 2 (Ω) ∀ c ∈ V, (5)
where γ is a positive constant chosen such that we have γD > 1. The hypothesis of the Lions theorem are then satisfied, and (1) has a unique weak solution c ∈ V .

Positivity of the solution.

As usual, we denote by c = c + -c -, for the solution c ∈ V , where c + and c -are the classical nonnegative parts of c. We have the following result:

Lemma 2.2 (Positivity). Let be c = c(x, t) the solution to the NTB system. Suppose that c 0 ≥ 0 and that c | Σ 1 ≥ 0, then c is non negative and we have:

c(T, .) ≥ 0, ∀ T > 0.
Proof. We will show that c -= 0. We multiply the NTB system by c -and we integrate by parts over Q. Since (∂ t c + )c -= (∇c + )c -= (∆c + )c -= 0, we have for each integral :

-α T 0 Ω ( ∂c - ∂t )c -dxdt = - α 2 Ω |c -(T )| 2 -|c -(0)| 2 dx ,
and with the Gauss fromula div q = 0 in Ω,

- Q (q.∇c -)c -dxdt = - T 0 1 2 Γ (q.c -)c -.n dx dt,
and,

D Q ∆c -c -dxdt = D Σ ∇c -c -.n dxdt -D Q ∇c -2 dxdt.
thanks to the Green formula. Adding the three integrals we find:

- α 2 Ω |c -(T )| 2 -|c -(0)| 2 dx + Σ1 I K |c -| 2 .n dxdt = D Q |∇c -| 2 dxdt ≥ 0. Now since c - | Σ 1 = 0 we obtain c -(T ) L 2 (Ω) ≤ c -(0) L 2 (Ω) = 0, ∀ T ≥ 0.
3. Optimal control for the NTB model. We here study the optimal control for the NTB system (1). Nutrients come from the second plant by exudates through Γ 2 , and there, we put a control function. First, we rewrite the NTB system using simple notations:

         A c = 0 in Q :=]0, T [×Ω, B c . n = I K c on Σ 1 =]0, T [×Γ 1 , B c . n = w on Σ 2 =]0, T [×Γ 2 , c(0, x) = 0 in Ω, (6)
where we have:

A = α ∂ ∂t + q.∇ -D∆, and B = D∇ - 1 2 q.
Here, w = -v where v is a positive control function. It corresponds to the addition of nutrient in the soil. It is prescribed in the rhizosphere frontier Γ 2 as an inflow nutrient.

The existence of a solution to (6) can be deduced from the one of Proposition 1, since we have : Corollary 1. Under the hypothesis of Proposition 1, we suppose that w ∈ L 2 (Σ 2 ). Then, there is a unique solution c ∈ V to the weak problem (2) with c 0 = 0, and with a(t; ., .) in ( 5) replaced by :

b(t; c, ψ) = a(t; c, ψ) - Γ2 w ψ dx, ψ ∈ V. (7)
The goal in this paper is to characterize the control which minimizes the cost function:

J(v) = c(v) -c 2 L 2 (Σ1) + N v 2 L 2 (Σ2) , (8)
where c is the observation given for L 2 (Σ 1 ), and where N > 0 is a positive constant.

The control function v is taken in the region L 2 (Σ 2 ) of Ω. The control problem is similar to a regional control problem as introduced by El Jai et al. [START_REF] Jai | Regional controllability of distributed systems[END_REF], and El Jai [START_REF] Jai | Analyse régionale des systèmes distribués[END_REF] for example.

Existence of the optimal control.

Proposition 2. There exists a unique u ∈ U = L 2 (Σ 2 ) such that

J(u) = inf v∈U J(v).
Proof. The existence of the control is easy to prove since the L 2 norm in J is continuous and coercive. Moreover, we have

J(v) ≥ J(0) = 0, ∀ v ∈ L 2 (Σ 2 ),
hence there exists a positive constant m > 0 such that m = inf

v∈L 2 (Σ2)) J(v).
Let v n be a minimizing sequence satisfying m = lim

n→∞ J(v n ). For every n ∈ N, n > n 0 , J(0) ≤ J(v n ) ≤ m + 1.
Hence, there exists a constant d ≥ √ m + 1 independent of n such that

v n L 2 (Σ2) ≤ d, and c(v n ) L 2 (Σ1) ≤ d.
Therefore v n u ∈ L 2 (Σ 2 ) weakly, and c(v n ) c(u) weakly in L 2 (Σ 1 ) because c is a continuous function.

Finally, we deduce from the strict convexity of the cost function J that u is unique.

3.2.

Characterization of the optimal control. We give a first characterization of the optimal control u ∈ L 2 (Σ 2 ) in the following lemma: Lemma 3.1. The optimal control u of the problem ( 6)-( 8) satisfies to :

Σ1 (c(u) -c) c(w).n dxdt + N Σ2 uw.n dxdt = 0, ∀w ∈ L 2 (Σ 2 ). ( 9
)
Proof. We have the Euler-Lagrange equality satisfied by the optimal control u:

lim λ→0 J(u + λw) -J(u) λ = 0, ∀w ∈ L 2 (Σ 2 ). Then J(u + λw) -J(u) λ = λ||c(w)|| 2 L 2 (Σ1) + 2λ c(u) -c, c(w) L 2 (Σ1) +2N u, w L 2 (Σ2) + λN ||w|| 2 L 2 (Σ2) . (10)
Hence,

lim λ→0 J(u + λw) -J(u) λ = 2 c(u) -c, c(w) L 2 (Σ1) + 2N u, w L 2 (Σ2) = 0,
which is the desired equality [START_REF] Jalonen | Transfer of nitrogen from a tropical legume tree to an associated fodder grass via root exudation and common mycelial networks[END_REF].

Γ h(c) u

Γ 2 1 Figure 2.
The Michaelis-Menten uptake function h(c) = Ic K prescribed on the root surface Γ 1 . The boundary control u(t, x), which represents an addition of nutrient into the rhizosphere region is prescribed on Γ 2 . Now, we give the optimality system (SO) for the optimal control u. We have the following theorem: Theorem 3.2 (Optimality system). The optimal control u for the problem ( 6)-( 8) is characterized by the triplet {u, c(u), p(u)} solution to the optimality system:

       A c = 0 A * p = 0 in Q, B c . n = I K c B * p . n = c(u) -c + I K p on Σ 1 , B c . n = -u B * p . n = 0 on Σ 2 , c(0, x) = 0 p(T, x) = 0 in Ω, (11)
and by the adjoint equation :

p + N u = 0 in L 2 (Σ 2 ). ( 12
)
Proof. We introduce the function p = p(t, x) solution to the adjoint problem :

       A * p = 0 in Q, B * p . n = c(u) -c + I K p on Σ 1 , B * p . n = 0 on Σ 2 , p(T, x) = 0 in Ω. ( 13
)
where

A * = -α ∂ ∂t -q.∇ -D∆, B * = -D∇ - 1 2 q.
We multiply the first equation of ( 13) by c(w), and we integrate by parts over Q. We first have :

-α Q ∂p ∂t c(w) dxdt = α Q ∂c ∂t (w) p dxdt. ( 14
)
With the Gauss property we have q.∇p = div (qp)-p div q = div (qp). We integrate over x, we have on one hand :

Q q.∇p c(w) dxdt = Q div (qp)c(w) dxdt = Q div (qp)c(w) -qp.∇c(w) dxdt = Σ qpc(w) ndxdt - Q qp.∇c(w) dxdt (15)
using the Gauss theorem. On the other hand, using the Green formula we have :

- Q D∆p c(w) dxdt = - Σ D∇p c(w).n dxdt + Σ D∇c(w) p.n dxdt - Q D∆c(w) p dxdt.
We resume by

0 = Q A * p c(w) dxdt = Q α ∂c(w) ∂t + q.∇c(w) -D∆c(w) p dxdt - Σ D∇p + 1 2 qp c(w).n dxdt + Σ D∇c(w) - 1 2 qc(w) p.n dxdt = - Σ1 c -c(u) + I K p c(w).n dxdt + Σ1 I K c(w)p.n dxdt - Σ2 wp.n dxdt = - Σ1 c -c(u) c(w).n dxdt - Σ2 wp.n dxdt. Hence, c(u) -c, c(w) L 2 (Σ1) = w, p L 2 (Σ2) . (16)
Replacing in (9), we finally have :

p + N u, w L 2 (Σ2) = 0, ∀ w ∈ L 2 (Σ 2 ).
4. Numerical results. To ease reader's reading and understanding in the interpretation of the below results we make the computations in one dimension of space.

In this situation the spatial motion of nutrients is produced along an horizontal axis ]x 0 , x Nx+1 [ on which the point at the right extrema corresponds to the uptake point at the root surface (Γ 1 = x Nx+1 ) and the point at the left extrema corresponds to the input point of nutrients (Γ 2 = x 0 In this case the optimal control problem is written by

[P ] =                              min v∈L 2 (Σ2) J(v)
where

J(v) = T 0 [c(t, x Nx+1 ) -c(t, x Nx+1 )] 2 dt + N v 2 L 2 (Σ2) |x Nx+1 -x 0 |
and where c satisfies the nutrient uptake model ( 6), with Σ 1 = (0, T ) × {x Nx+1 } and Σ 2 = (0, T ) × {x 0 }.

The function c is a given nutrient uptake dynamics at the root surface and v is the searching control that solve the problem [P ].

The optimal control problem [P] is numerically solved by a Quasi-Newton method whose the algorithm is available on Matlab tools or on Absoft software (see Picart and Ainseba [START_REF] Picart | Parameter identification in multistage population dynamics model[END_REF] who used the method for hyperbolic systems). Here we used the latter with the BCONF routine which is well adapted for our problem. The solution of NTB-system is numerically computed with a Difference Element scheme of Euler kind that is reported in Annex. The parameters of the NTB-system are chosen so that a spatial motion is observed on the time and space domain (0, 5) × [0, 1] that is q = 0.5, D = 0.1, α = 1, I = 0.001, K = 1, with a space-discretization step equal to 0.1, and a time-discretization step equal to 0.01. We numerically test our problem with the function c of Figure 3. The root nutrient uptake is like a sigmoid curve that increases with time until to reach a maximum nutrient uptake value of 2.8. The function c is obtained by solving NTBsystem with a given control function that is represented in dotted curve on Figure 4 and called exact control function. The input of nutrients is a step-function equal to one during the first 4 units of time and to 0 after that. The numerical result consists in finding for this given nutrient uptake dynamics c the optimal control solution of problem [P ] that allows to reach it. For this, we initialize the Quasi-Newton algorithm with the step-function represented by the double-dotted curve of Figure 4. The optimal control solution is then obtained when the value of the cost function J(v) at the optimal control solution is close to zero.

The optimal solution is given by the solid curve of Figure 4. The L 2 -error of the difference between the estimated solution and the exact one (dotted curve) is of order 1.5910 -3 meaning that the QN's algorithm has converged toward the exact solution. In this example the numerical optimal solution is unique.

We conclude this work by a final remark:

Remark 2. In this paper we proposed a model of nutrient uptake by roots whose origin is either the soil or root exudates from a secondary plant [START_REF] Griffon | Nourrir la Planète[END_REF]. This model derives from the NTB model which describes the transport and diffusion of nutrients through the soil to the root surface in general. These initial theoretical and numerical results are encouraging for the study and development of alternative methods to fertilizing plants through the soil or the contribution of service plants for example. 5. Annex. Let Ω =]x 1 , x Nx [ be the domain of the space variable, with N x the space-point number. Let x i = i∆x be the space discretization point, for i = 1, N x , with ∆x the space discretization step. The time discretization point is defined by t n = n∆t, with n ≥ 0 and ∆t the time discretization step. The numerical solution of the NTB-system of ( 6) is given by the following scheme c n+1 i -c n i ∆t = F (c n i , x i , t n ), for i from 1 to N x and n ≥ 0, where

F (c n i , x i , t n ) = - q α c n i -c n i-1 ∆x + D α c n i+1 -2c n i + c n i-1 ∆x 2 .
We get a system of equations that can be rewritten as the following matrice form

C n+1 = ∆t α AC n ,
where A is given by The boundary conditions at the point x 0 and x Nx+1 are respectively approximated by

         α ∆t -q ∆x -2D
c n x0 = c n 2 + ∆x v n - q 2 c n 1 , c n x Nx +1 = c n x Nx-1 + ∆x h n Nx + q 2 c n
Nx , for all time n ≥ 0.

Figure 1 .

 1 Figure 1. The rhizosphere Ω is delimited by the two boundaries Γ 1 and Γ 2 , where Γ 1 represents the root surface and where Γ 2 the boundary of the rhizosphere where the nutrients enter.

Figure 3 .

 3 Figure 3. Nutrient uptake at the root surface with respect to time. This curve is used to define the c function in problem [P ].

Figure 4 .

 4 Figure 4. Nutrient input from exudates with respect to time. The dotted curve is the control function used to construct the function c by solving the NTB-system. The double dotted curve is the initial control used in the Quasi-Newton algorithm. The solid curve is the numerical solution of [P ].


N Le Touz

J Dumoulin

J.-M Piau

MULTI-PHYSICS FEM MODEL OF SOLAR HYBRID ROADS FOR ENERGY HARVESTING PERFORMANCE EVALUATION IN PRESENCE OF SEMI-TRANSPARENT OR OPAQUE PAVEMENT SURFACE LAYER

Keywords: Computational methods, Energy efficiency, Multi-physics heat transfer, Solar hybrid road, Energy harvesting, Finite element model, Porous media

We present in this paper the concept of solar hybrid road and focus on the thermal performances of such system. Main differences between these multi-layer structures with traditional road structures come from the pavement surface, which can be opaque or semi-transparent and the porous sub-layer, submitted to fluid flow. These structures aim at collecting solar energy during the summer season, what helps also to prevent from a too high temperature and urban heat island effect. A finite element model is presented to couple thermal diffusion, hydraulic convection and radiative transfer. This numerical model allows to compute the temperature field for different weather conditions and also to evaluate the thermal performances of the system. Annual simulations are performed and a comparison between two surface layer solutions for different locations and climates is presented and discussed.

INTRODUCTION

As part of the energy transition, numerous studies have been led on renewable energy sources and energy management. Because of the increase of car traffic, road structures represent a large surface area that is still growing. These structures are in direct contact with the environment and subject to weather conditions. During summer conditions, road are submitted to heavy thermal load through solar and convective heat flux. The increase of temperature in the structure can provoke a degradation of mechanical performances entailing rutting and also contribute to the rise of urban heat island effect [START_REF] Hasebe | Thermoelectric generators using solar thermal energy in heated road pavement[END_REF][START_REF] Mallick | Harvesting energy from asphalt pavements and reducing the heat island effect[END_REF]. During winter conditions, the surface temperature can decrease below 0 • C what may entail black ice formation, and consequently safety and economic problems. Moreover, current practice consisting in salting the road entails environmental problems.

In order to avoid these problems, we propose to control the temperature of the road structure by a reversible system able to recover energy during the summer and bringing energy during cold periods. Traditional approaches to build such embedded energetic systems rely on a heat fluid flowing in tubes buried below the surface layer of the pavement [START_REF] Liu | Modelling snow melting on heated pavement surfaces, part i: model development[END_REF][START_REF] Mallick | Harvesting energy from asphalt pavements and reducing the heat island effect[END_REF]. In this study, we present a solar hybrid road in which the fluid flows into a porous layer below the surface layer. To increase the ability of this hybrid road to harvest energy, a semi-transparent layer is used in order to favour solar radiation penetration in the structure.

IHTC16-23896

In this paper, computation of the temperature field in the structure with a multi-physics finite element model is first presented. This numerical model is then applied to calculate the energy harvesting during the summer period.

HYBRID SOLAR ROAD CONCEPT

This concept is based on a multilayer system for which the two layers closest to the surface have been modified, compared to traditional road structures. The multilayer system is shown in figure 1. Thermal properties of materials considered, in this study, for numerical simulations are reported in table 1. For optical properties, we take in a first approach properties of glass given by [START_REF] Ping | Transient radiative-conductive heat transfer in flat glasses submitted to temperature, flux and mixed boundary conditions[END_REF] and estimate them thanks to experimental mock-up [START_REF] Letouz | Infrared thermography applied to the study of heated and solar pavements: from numerical modeling to small scale laboratory experiments[END_REF]. They are reported in table 2.

Table 2 Optical properties of the semi-transparent layer, from [START_REF] Ping | Transient radiative-conductive heat transfer in flat glasses submitted to temperature, flux and mixed boundary conditions[END_REF] for absorption coefficient κ λ

[µm] λ < 0.5 0.5 ≤ λ < 2.7 2.7 ≤ λ < 4.5 4.5 ≤ λ < 50 λ ≥ 50 κ λ [m -1 ] 0 10 1000 5000 +∞ σ λ [m -1 ]
10 10 10 10 10

MULTI-PHYSIC NUMERICAL MODEL

Our model combining thermal diffusion, hydraulic convection and radiative transfer is presented hereafter

Thermal diffusion

The structure is submitted to convective exchange with the air, radiative exchange with the sky and solar radiation. Heat equation can be written as system 1. Depth of the structure is supposed to be large enough so

IHTC16-23896

that thermal exchanges can be neglected except at the surface S, shown in figure 2, where Robin boundary conditions are considered.

     ρC p ∂T ∂t = ∇ • [k∇T ] + q k∇T • n = Φ s + h conv (T a -T ) + h rad T sky -T on the surface 0 for other boundaries (1)
With ρ the density, C p the heat capacity, k the thermal conductivity, h conv the convective heat exchange coefficient, h rad the radiative exchange coefficient, T a the air temperature and T sky the sky temperature. The radiative exchange coefficient is given from sky and surface temperatures by h rad = εσ(T 2 surf +T 2 sky )(T surf +T sky ), where ε is the surface emissivity and σ the Stefan constant.

Finite element formulation of this problem has been widely studied, more information can be found in [START_REF] Gartling | The Finite Element Method in Heat Transfer and Fluid Dynamics[END_REF][START_REF] Taine | Transferts thermiques, introduction aux transferts d'nergie[END_REF].

Hydraulic convection

Two temperature fields are considered in the porous layer, for solid and fluid phases, respectively written T and T f . Equations on these temperature fields can be obtained from Navier-Stokes equations integrated over an elementary representative volume. With the hypothesis of a stationary flow, a low and uniform pore speed and a saturated porous layer [START_REF] Bejan | Porous and Complex Flow Structures in Modern Technologies[END_REF][START_REF] Nield | Convection in porous media[END_REF], we have:

(1 -ϕ) ρ s C p,s ∂T ∂t = (1 -ϕ) ∇ • [k∇T ] + h (T f -T ) (2)
ϕρ f C p,f ∂T f ∂t + ρ f C p,f u • ∇T f = ϕ∇ • [k f ∇T f ] + h (T -T f ) (3)
Where h is an exchange coefficient dependent of the pore and flow characteristics [START_REF] Nield | Convection in porous media[END_REF], u is the fluid velocity and ϕ is the porosity.

Boundary conditions for fluid temperature are supposed to be Dirichlet for the input surface Γ f,in and adiabatic anywhere else.

Because of the advection term in equation ( 3), Galerkin method applied to the weak formulation of ( 2) and (3) can entail spurious oscillations [START_REF] Johnson | Numerical solutions of partial differential equations by the finite element method[END_REF]. To avoid that, a streamline upwind Petrov Galerkin (SUPG) formulation is used, allowing to add numerical diffusion [START_REF] Brooks | Streamline Upwind/Petrov-Galerkin Formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes Equations[END_REF][START_REF] Letouz | Etude numerique de la resolution du couplage convection/radiation/diffusion dans une structure de chaussée hybride[END_REF].

Radiative transfer

Radiative transfer equation [START_REF] Howell | Thermal radiation heat transfer[END_REF] is presented in equation ( 4) :

ϑ • ∇I ν = -(κ ν + σ ν ) I ν + κ ν I b,ν + σ ν 4π 4π p ν ( θ, ϑ)I ν d θ ( 4
)
Where x is the spatial position, ϑ the direction, I ν the intensity, I b,ν the blackbody intensity, κ ν the absorption coefficient, σ ν the scattering coefficient and p ν the phase function.

Two discretizations are applied to the radiative transfer equation: angular and spatial. For angular discretization, discrete ordinates method is used on an icosahedron mesh, which can be refined recursively. A system of semidiscrete equations, each with an advection term, is also obtained. As for the hydraulic convection, SUPG

IHTC16-23896

method is used for spatial discretization [START_REF] Kanschat | Solution of radiative transfer problems with finite elements[END_REF][START_REF] Letouz | Etude numerique de la resolution du couplage convection/radiation/diffusion dans une structure de chaussée hybride[END_REF]. Basis function ϕ of Galerkin weak formulation are replaced by ϕ + δϑ • ∇ϕ for Petrov-Galerkin weak formulation, where δ is a coefficient function of the mesh size.

A linear system is then solved. Inversion of this system allows to get the intensity field over the semi-transparent domain, and also the energy that is absorbed in.

Materials that are used here are semi-transparent for main wavelengths of solar radiations, but opaque for higher wavelengths, as those of the ambient radiation. In particular, blackbody emission at ambient temperature for main wavelengths of solar radiation are negligible, respect to solar radiation. Source terms can also be computed for a unit solar heat flux, independently of the temperature field.

Numerical coupling

These three finite element models are coupled in order to compute the evolution of the temperature field over the whole structure by knowing the climatic conditions. At each time step, knowledge of solar radiation and temperature field allow to compute, from the linear system of radiative transfer equation, the intensity and then the energy absorbed in the semi-transparent media. Heat and hydraulic convection are then solved by taking into account the thermal heat source induced by solar radiation. A Crank-Nicolson scheme allows finally to compute the temperature field at the next time step [START_REF] Letouz | Etude numerique de la resolution du couplage convection/radiation/diffusion dans une structure de chaussée hybride[END_REF][START_REF] Letouz | Infrared thermography applied to the study of heated and solar pavements: from numerical modeling to small scale laboratory experiments[END_REF].

ENERGY HARVESTING PERFORMANCES EVALUATION

In this section, we first present an overview of thermal performances of solar hybrid roads for some locations in France. From annual weather database, we calculate the temperature field distribution with time and also the output fluid temperature. The energy harvesting is calculated from the fluid temperature rise between input T f,in and output T f,out :

E = t ρ f C p,f uS (T f,out -T f,in ) + dt (5)
Where S is the section of the flow, u = 1 × 10 -4 m.s -1 is the Darcy velocity and (•) + refers to the positive part.

We consider a length of flow equal to 4 m. The studied period for energy harvesting comes from May to October. In fact, because of thermal inertia, this is the period when the temperature of the structure is the higher during the year in France. Results for energy harvesting are presented in table 3.

Table 3 Results for energy harvesting For all these locations, a non negligible part of solar radiation can be recovered, this part being more important when solar radiation and/or average temperature are higher. Using a semi-transparent surface layer allows to obtain an increase in energy harvest between 54 and 76 kWh.m -2 for studied locations, corresponding to an efficiency increase of about 7 % by reference to the incident radiation. However, it should be noticed that the energy which is recovered per square meter decreases with the length of the unidirectional porous exchanger.

Location total solar irradiation [kWh/m 2 ] E (opaque) [kWh/m 2 ] η (opaque) [%] E(ST) [kWh/m 2 ] η(ST) [%] δE [kWh/m 2 ]
To illustrate that, we present the efficiency distribution as a function of length for Agen in figure 5.

Fig. 5 Efficiency as a function of length in Agen

In a second step, our model was applied to weather data from Naples (Italy), Montreal (Canada) and Bergen (Norway). In Naples, where mean weather data are closed to Nice, the efficiency increases from 41.5 % with opaque surface layer to 48.3 % with semi-transparent surface layer. In Montreal, where the mean temperature

IHTC16-23896

is lower than in France while solar irradiance is closed to Agen, efficiency is lower, with respectively 29.0 and 35.9 %. In Bergen, both mean temperature and irradiance are lower than in France, entailing low efficiencies, equal to 17.6 and 23.2 %. It confirms results and trends observed at France level.

CONCLUSION

In this paper, we have presented a concept of solar hybrid road and its multi-physics finite element model to solve the coupling between heat diffusion, hydraulic convection and radiative transfer. This numerical model was used to calculate the thermal behavior of such structure in different weather conditions. Energy harvesting in different locations in France is presented and results shows that a non negligible amount of solar energy can be recovered. Using a modified semi-transparent surface layer allows to increase the average energy harvesting of round 7 % by reference to the average solar irradiance in France. In future works a command law to control the surface temperature by acting on the input fluid temperature will be studied. Finally, an improvement of our model would be to take into account surface precipitations in the boundaries conditions.

Fig. 1

 1 Fig. 1 Sectional schematic view of modified road structure Fig. 2 Interactions with the environment

  semi-transparent and opaque surface layer and δη is the efficiency increase given respect to the incoming total solar irradiation. Results obtained for the French geographical distribution in energy harvesting and efficiency are shown respectively in figures 3 and 4.

Fig. 3 Fig. 4

 34 Fig. 3 Energy harvesting for opaque and semi-transparent surface layer

Table 1

 1

	Thermal properties of layers and fluid (from [1, 7, 16])
		layer 1 layer 2 layer 3 layer 4 fluid
	k [W.m -1 .K -1 ]	0.85	1.03	1.40	1.10	0.60
	ρ [kg.m -3 ]	2700	2360	2620	2300 1000
	C p [J.kg -1 .K -1 ]	840	780	860	1000 4180

  Efficiencies η are given as the energy harvest over the cumulated solar irradiation. Subscripts opaque and ST refers respectively to the cases of opaque and semi-transparent surface layer. δE is the difference between

								δη
								[%]
	Agen	903	339	37.5	399	44.2	60	6.7
	Carpentras	1034	408	39.5	479	46.3	71	6.8
	La Rochelle	961	309	32.2	381	39.6	72	7.4
	Mâcon	822	312	38.0	367	44.7	55	6.7
	Nancy	786	267	34.0	321	40.8	54	6.8
	Nice	1035	424	41.0	500	48.3	76	7.3
	Rennes	838	261	31.1	320	38.2	59	7.1
	Trappes	833	275	33.0	330	39.6	55	6.6


N Le Touz

M Marchetti
email: mario.marchetti@cerema.fr

J Dumoulin

L Peiffer

A Escal

L Ibos

M Fois

P Bourson

Appreciation of the delay in the benefits of the thermal energy released by PCM in civil engineering structures

To optimize de-icers use in winter on roads and airports, some sustainable solutions were promoted in a global project entitled "Roads of the Future". PCM into pavements were then considered to rely on the released enthalpy to delay the occurrence of black ice or snow accumulation. Coupling a spectroscopic monitoring of the PCM behavior and thermography on pavement surface indicated a significant delay between the phase change and the thermal manifestation of the released energy. Considering an accurate description of PCM according to cooling weather phenomena, these elements might provide valuable information for numerical models.

Introduction

Maintaining mobility on roads and the possibility for aircrafts to land and to take-off is a major issue in poor weather conditions as encountered in winter. Many efforts have been developed over the years to evaluate the economical costs of winter maintenance [START_REF] Petty | Federal Highway Administration winter road Maintenance Decision Support System (MDSS): Recent enhancements and refinements[END_REF], its environmental impacts [START_REF] Suaire | Transfer, Exchanges and Effects of Road Deicing Salts in a Detention Pond Treating Road Water[END_REF] or to forecast these adverse situations [START_REF] Khalifa | Appreciation of the traffic effects on the RST by infrared thermography[END_REF][START_REF] Marchetti | Methodology to forecast road surface temperature with Principal Components Analysis and Partial Least-Square regression. Application to a urban configuration[END_REF][START_REF] Khalifa | Accounting for anthropic energy flux of traffc in winter urban road surface temperature simulations with TEB model[END_REF][START_REF] Khalifa | Individual contributions of anthropogenic physical processes associated to urban traffic in improving the road surface temperature forecast using TEB model[END_REF]. Still the amounts of de-icers are increasing, conducting to consider a modification of the infrastructure itself to make it compliant with the avoidance of ice occurrence and snow accumulation, as performed on the Serso bridge in Switzerland [START_REF] Eugster | Road and Bridge Heating Using Geothermal Energy. Overview and Examples[END_REF]. Previous research indicated the benefits of a multilayered asphalt concrete [START_REF] Marchetti | Infrared thermography applied to surface temperature survey of enhanced pavement sample addressing self-de-icing functionality[END_REF] and phase change materials (PCM) included into a pavement [START_REF] Marchetti | Comparative study in the identification of liquid to solid transition phase with DSC, Raman spectra analysis and chemiometrics methods applied to phase change materials used for icing-delay in civil engineering infrastructures[END_REF]. These inclusions are meant to maintain a constant surface temperature during the phase transition process. The phase transition did indeed take place but often not at the temperature provided by the PCM manufacturer. Furthermore, a delay is expected between the moment the phase transition occur and its thermal detection at the surface. Coupling a spectroscopic monitoring of the PCM behavior and a thermographic one of the pavement surface allowed to quantifying subcooling effect and this delay in the benefits of the energy released during the phase change. These elements could therefore be implemented in numerical models to properly appreciate to what extent these PCM could contribute to a sustainable winter maintenance.

Coupling infrared thermography and Raman spectroscopy to pavement monitoring

The main idea was to identify to what extent the energy released during the solid to liquid phase transition of a PCM embedded into a pavement was able to increase the pavement surface temperature. But because of the PCM volume and the poor and heterogeneous thermal characteristics, the determination of the time necessary for this available released energy to reach the surface was difficult to appreciate. This time was difficult to identify for two reasons. First, the PCM being embedded into a non transparent material, the phase change can not be observed. The second aspect is related to subcooling phenomena that can hardly be predicted since they depend on the experimental conditions. So as to identify the moment the phase transition is taking place, a Raman probe was inserted within the PCM container, the spectral signature being different according to the phase of the material.

Pavement, PCM samples and test facility

Two pavement samples were manufactured, and were parallelepipeds of 10 cm x 18 cm x 50 cm, of semi-granular asphalt concrete (granular materials with a bitumen matrix). This material is very common either on national roads, highways and taxiways. The first one was a pristine asphalt concrete one used as reference. Some wood inclusions of different heights were placed into the second sample during the manufacturing to generate voids where PCM inclusions in appropriate metallic containers were placed (figure 1a). One of the inclusions is then located 3.9 cm below the surface and the second one 2.8 cm. All pavement samples sides were surrounded by extruded polystyrene to provide a thermal insulation of the samples and to avoid lateral losses.

. Pristine pavement sample and PCM metallic containers before installation of pavement sample with voids wi t h T -type thermocouples (a) and pavements samples with FLIR A320G infrared camera (b)

Both pavement samples were then placed on the floor of a climatic chamber before being submitted to a dry air temperature variations (-10°C ≤ Tair ≤ +10°C) representative of a winter in France. The test lasted five days, excluding a preliminary thermal conditioning at 10°C for 24 hours. Several thermocouples were attached to both pavement samples and PCM containers to monitor temperature variations during the test.

The selection of PCM relied on some requirements. The first one was to have a liquid to solid phase transition temperature compliant with the winter maintenance context on civil engineering infrastructures. A maximum of 5°C was then selected. The second aspect to be considered are the chemical characteristics of the materials, i.e. the potential corrosivity with respect to metals, or their flash points. Therefore, a first PCM inclusion consisted in an aqueous saline solution and with a freezing point of -5°C. The second one was a Rubitherm ® SP5, an inorganic material with a freezing temperature in the 2-7°C range according to its datasheet.

Thermographic and spectroscopic monitoring

Infrared images were generated with a FLIR A325 (320 x 240 IRFPA, LWIR, pitch 25 µm), with a NEDT of 0.1°C. The camera was placed in the climatic chamber with an incidence of roughly 20° with respect to the pavement samples surfaces, at a distance of 1.5 m (figure 1b).

The spectroscopic monitoring of the aqueous saline solution as PCM inclusion was obtained with an i-Raman from BWTek ® operating with a 532 nm laser light, and with a spectral resolution of 2 cm -1 . An immersion probe was sealed into the container with the PCM, the signal being collected through an optical fibre, the passage being under the pavement samples. One of the great advantage of such experimental configuration is the total absence of external light. Within these conditions, the Raman spectral signature of the PCM was not containing any fluorescence artefacts.

Both instruments were synchronized, with an infrared and Raman spectrum collection every 3 minutes. The over 2000 infrared images were analysed with ThermaCam Researcher ® and its associated OLE Automation implemented to extract surface temperature of some transects over pavement samples. Both transects and over 2000 Raman spectra were then analysed by means of principal component analysis (PCA) with Unscrambler ® X 10.1 software, using a NIPALS algorithm, as detailed in the literature [START_REF] Marchetti | Comparative study in the identification of liquid to solid transition phase with DSC, Raman spectra analysis and chemiometrics methods applied to phase change materials used for icing-delay in civil engineering infrastructures[END_REF] to identify the PCM phase transition. The same tool was applied to thermal images considering the surface temperature increase will generate an anomaly that a PCA analysis will identify.

Finite elements and numerical simulations

We propose in this section to perform numerical simulations to compare results issued from measurements with numerical ones. A finite element model was developed to compute the temperature field in the structure from the knowledge of air temperature evolution at boundary. Two pavements samples were studied: a pavement with two inclusions, Rubitherm ® SP5 and brine solution, shown in figure 2, and a healthy pavement without inclusion and with the same dimensions than experimental pavement samples.

Thermal diffusion in the structure is governed by the heat equation ( 1):

𝜌𝜌𝜌𝜌 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 -∇ ⋅ (𝑘𝑘∇𝜕𝜕) = 0 (1)
Where 𝜌𝜌 is the density, 𝜌𝜌 the massic heat capacity and 𝑘𝑘 the thermal conductivity. We suppose that thermal exchanges of the structure with the outside occur at the surface. Other boundaries are supposed to be adiabatic. Boundary conditions are as follows:

𝑘𝑘∇𝜕𝜕 ⋅ 𝑛𝑛 �⃗ = ℎ(𝜕𝜕 ∞ -𝜕𝜕) (2)
With ℎ an exchange coefficient and 𝜕𝜕 ∞ the temperature far from the structure. Exchanges between the studied domain and the environment take into account convective exchanges with the air and radiative exchanges. Radiative exchanges are supposed to be linearized and environment temperature equal to the measured air temperature.

Applying the finite element method to equations ( 1) and ( 2) allows to get a linear system [START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF]. For the time dependence, we used a Crank-Nicolson scheme.

Fig. 2. Schematic view of the domain

We suppose that thermal properties of the asphalt concrete and insulating material are constant. From [START_REF] Asfour | Hydrothermal modeling of porous pavement for its defreezing[END_REF], we use 𝜌𝜌 = 2450 kg.m -3 , 𝜌𝜌 = 920 J.kg -1 .K -1 and 𝑘𝑘 = 1.40 W.m -1 .K -1 for asphalt concrete and from [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF] 𝜌𝜌 = 55 kg.m -3 , 𝜌𝜌 = 1210 J.kg -1 .K -1 and 𝑘𝑘 = 0.027 W.m -1 .K -1 for insulating material (extruded polystyrene).

Because of the phase change of both brine and SP5, thermal properties of the other materials cannot be considered as constant over the studied temperature range, and particularly along the temperature range of phase change.

Thermal data for SP5 are provided by the supplier [13] and are shown in figure 3. Thermal conductivity is moreover supposed to vary linearly with the temperature.

For the brine solution, we consider a two-phase mixture. The salt quantity is equal to 𝐶𝐶 𝑖𝑖 ,0 = 0.080 kg per kg of water, below the saturation value. For a decreasing temperature from a positive temperature, the mixture is first liquid, then, the eutectic temperature is reached and ice begin to form. Salt concentration in the remaining brine solution also increases. Ice formation continues until the eutectic point, at -20.6°C where the NaCl solution reaches the saturation. If temperature continues to decrease, mixture is composed of ice and solid sodium chloride dihydrate (NaCl-2H2O). Phase change occurs also between the eutectic point at -20.6°C and eutectic temperature. We consider here a concentration in NaCl equal to 80 g per kg of water. From values given by [START_REF] Aguirre-Puente | Contribution à l'étude du changement de phase de solutions aqueuses[END_REF], the eutectic temperature 𝜕𝜕 𝑙𝑙 (in °C) can be interpolated as a function of NaCl concentration 𝐶𝐶 (in kg per kg of water):

𝜕𝜕 𝑙𝑙 (𝐶𝐶) = -55.3𝐶𝐶 -42.05𝐶𝐶 2

Thermal capacity can be computed from infinitesimal variations of enthalpy 𝑑𝑑𝑑𝑑 that come from variations of temperature for ice and brine and latent heat enthalpy from phase change of brine:

𝑑𝑑𝑑𝑑 = 𝑓𝑓 𝑖𝑖 𝜌𝜌 𝑖𝑖 𝑑𝑑𝜕𝜕 + 𝑓𝑓 𝑏𝑏 𝜌𝜌 𝑏𝑏 𝑑𝑑𝜕𝜕 + 𝐿𝐿𝑑𝑑𝑓𝑓 𝑏𝑏 (4)
Where 𝑓𝑓 𝑖𝑖 is the massic fraction of ice in the mixture, 𝑓𝑓 𝑏𝑏 = 1 -𝑓𝑓 𝑖𝑖 the massic fraction of brine, 𝐿𝐿 the latent heat, 𝜌𝜌 𝑖𝑖 the massic capacity of ice and 𝜌𝜌 𝑏𝑏 the massic capacity of brine.

Values for 𝜌𝜌 𝑖𝑖 ,𝜌𝜌 𝑏𝑏 and 𝐿𝐿 are given by [START_REF] Aguirre-Puente | Contribution à l'étude du changement de phase de solutions aqueuses[END_REF] as a function of temperature for 𝜌𝜌 𝑖𝑖 and 𝐿𝐿 or concentration for 𝜌𝜌 𝑏𝑏 .

From the definition of brine fraction, we have:

𝑑𝑑𝑓𝑓 𝑏𝑏 = - 𝑓𝑓 𝑏𝑏 2 𝐶𝐶 𝑖𝑖,0 𝑑𝑑𝐶𝐶 (5)
And 𝑑𝑑𝐶𝐶 can be computed from (3) as 𝑑𝑑𝐶𝐶 =

𝑑𝑑𝑇𝑇 𝑙𝑙 𝑑𝑑𝑇𝑇 𝑑𝑑𝜕𝜕
We suppose here that thermal capacity is constant over temperature range of 1°C. The thermal capacity over the temperature range is obtained by integrating 𝑑𝑑𝑑𝑑 over the temperature range:

𝜌𝜌(𝜕𝜕 1 ≤ 𝜕𝜕 < 𝜕𝜕 2 ) = 1 𝜕𝜕 2 -𝜕𝜕 1 � 𝑑𝑑𝑑𝑑 𝑇𝑇 2 𝑇𝑇 1 (6)
As for data from SP5 supplier, we work with temperature ranges of 1°C for the brine solution. Outside of the phase change temperature range, that is to say to temperatures higher than -4.7°C.

We have the partial enthalpy distributions shown in figure 3. The enthalpy distribution for SP5 has an hysteresis: values for enthalpy distribution are different along the direction of temperature variation.

For thermal conductivity, it is supposed to be constant for SP5, equal to 0.2 W.m -1 .K -1 , and function of both temperature and NaCl concentration for brine and ice. From [START_REF] Levy | The thermal conductivity of commercial brines and seawater in the freezing range[END_REF]:

𝑘𝑘 𝑖𝑖𝑖𝑖𝑏𝑏 = 2.25 -0.006875𝜕𝜕(1 -0.0159𝜕𝜕) 𝑘𝑘 𝑏𝑏 𝑏𝑏𝑖𝑖 𝑏𝑏𝑏𝑏 = (1 -0.14𝐶𝐶)�0.562 + 0.001967𝜕𝜕(1 -0.00338𝜕𝜕)� ( 8
)
Air temperature is shown in figure 4. In a first approach, we consider that the environmental temperature is at air temperature for radiative exchanges.

Fig. 4. Air temperature as a function of time

Results obtained by numerical simulation were conducted the air temperature evolution measured during experiments and presented in figure 4. Results obtained are presented and discussed in the following paragraph.

Results and discussion

Identification of the PCM phase transition

Raman and thermographic data acquisition conducted each to the generation of over 2000 respectively spectra and thermal images, covering the whole duration of the test. The phase change in the PCM could either be detected through a change in the Raman spectra signature, and through a temperature change detected at the surface of the pavement sample containing the inclusions.

When considering the surface temperature evolution monitored by the infrared camera (figure 5a), it appeared a clear delay of over 10 hours before the surface temperature reached more or less the air temperature. This clearly indicated the thermal inertia of the material. The temperature variation with time also indicated two moments with some heterogeneities in its evolution, indicated in figure 5a, These heterogeneities do correspond, with a timeshift, to thermal anomalies detected on the surface (figure 5b), which could be attributed to the phase change, with respectively a reduction in temperature decrease and increase. This is consistent with the fact that in one case some heat is released during the liquid to solid phase transition, while some energy is absorbed to the solid to liquid phase transition.

PCA results on thermal images (a) and Raman spectra (b) for the whole test duration for the saline solution inclusion

The main interest was to quantify the time shift between the phase transition as detected by the Raman spectrometer and its thermal manifestation at the pavement surface. This was performed with a PCA conducted on data regarding a temperature decrease followed by a plateau including the phase transition temperature. Results are illustrated in figure 7. As indicated, the time shift between the phase transition detected spectroscopically and its thermal occurrence on pavement surface was of nearly 9 hours. This time shift was also identified for the fourth temperature decrease and seventh plateau. This time shift is quite large, and this can easily be explained by the low thermal conductivity of the material the pavement samples were made of. This low thermal conductivity is also caused by the porosity of the pavement sample.

Comparison between measurements and thermal numerical simulations

Temperature field computed from the model is compared with measurements realized with infrared thermography and with thermocouples inserted at surface level and inside samples. Numerical and experimental probes locations are reported in figure 8.

Fig. 8. Sectional view of pavement samples and probes locations with and without inclusion

Temperature at the surface of inclusions and for healthy pavement from measures and computation are shown in figure 9. Even if the influence of phase change is particularly visible close to the inclusions, the temperature at the surface is also impacted by it. As for figure 8, the computed temperature at the surface of the pavement is shown in figure 11. Finally, a first comparative analysis between reconstructed thermal image sequence issued from numerical simulations and experimental ones is proposed.

Simulated thermal image sequence are obtained by interpolating the surface temperature field evolution, issued from 3D numerical simulations, and using a 1 mm spatial step interpolation. Two sequences were generated from numerical calculations, one for the specimen with the two PCM inclusions and one for the specimen with no PCM inclusion. Using the same decomposition, the experimental thermal image sequence, has been spatially stich in two thermal sequences, one for with the two PCM inclusions and one for the specimen with no PCM inclusion.

A principal component analysis by using Singular Value Decomposition [START_REF] Rajic | Principal component thermography for flaw contrast enhancement and flaw depth characterization in composite structures[END_REF] was applied to the first air temperature decay. Results obtained for second empirical orthogonal function (EOF) maps are presented in figure 15. For the specimen with two inclusions, (b) and (d) difference of behaviour depending of the nature of PCM materials can be observed. Difference of behaviour is clearly visible in second EOF map calculated using numerical simulations and more difficult to observe on experimental one. Two perturbations are affecting the experimental results, one is due to convective effect in the climatic chamber used for experiments (see upper part of map) and the second one is induced by thermocouples glued at specimen surface. Furthermore, numerical simulations consider that the specimen is homogeneous, but experiments are affected by noise induced by the mix of aggregates and binder.

Conclusion

The objective of this research was to investigate to what extent the solid to liquid phase transition of Phase Change Material (PCM) inclusions into pavement were able to generate enough heat to be considered for the avoidanc e of ice occurrence and reducing snow accumulation. One asphalt concrete received two PCM inclusions and was submitted to controlled air temperature cycles representative of winter. The surface temperature variations were compared to another pristine asphalt concrete sample. The main concern was to identify the moment where the phase transition took place to clearly identify the timeshift between the energy release and its availability at the pavement surface. To do so, a spectroscopic Raman monitoring was specifically dedicated to the PCM, while an infrared camera recorded the variations of the surface temperatures of both pavement samples.

The analysis of both Raman spectra and infrared images by means of a multivariate data analysis tool clearly allowed to detecting both the phase transition and its surfacic thermal effect. The timeshift reached nine hours in the selected experimental conditions.

A new set of experiments will be planed in a near future to confirm trends observed by using SVD analysis. In particular, duration of constant cold and hot periods will be extended to favour a complete phase change in PCM material, along with an enhancement of the thermal signature effect.

Fig. 3 .

 3 Fig. 3. Partial enthalpy distributions for brine (left) and SP5 (right) Density of RT5HC is supposed to evaluate linearly, with 𝜌𝜌 = 880 kg.m -3 for 𝜕𝜕 = -15°C and 𝜌𝜌 = 760 kg.m -3 for 𝜕𝜕 = +20°C. For brine solution and ice, we use data from [Levy 1982]: 𝜌𝜌 𝑏𝑏𝑏𝑏𝑖𝑖𝑏𝑏𝑏𝑏 = �1 + 0.689𝐶𝐶 (1 -0.0056𝜕𝜕)� × 10 3 𝜌𝜌 𝑖𝑖𝑖𝑖𝑏𝑏 = 0.917(1 + 0.000173𝜕𝜕) × 10 3 (7)

Fig. 5 .

 5 Pavement surface temperature evolution above and on the saline solution-PCM inclusion (a) and thermal image of the pavements samples shortly before an air temperature increase (b)Such a phase change could also indeed be detected through a PCA on both spectral and Raman data. PCA results of this data is presented in figure6. For both thermographic and spectroscopic data, results indicated a structure depending on the air temperature decreases, plateaus or increases.

Fig. 7 .

 7 PCA results on Raman spectra (a) and thermal images (b) for second temperature decrease and the third plateau

  InfraRed Thermography Conference, 25 -29 June 2018, Berlin, Germany 6

Fig. 9 .

 9 Fig. 9. Thermocouple measurements and computed surface temperature above inclusions The effects of the inclusions on the temperature field are visible when the internal temperature is closed to values of state change. Temperature field computed from our model is shown on a vertical sectional view for 𝜕𝜕 = 25 hours in figure 10 on the left. The difference with the pavement without inclusion is also illustrated in this figure on the right.

Fig. 10 .

 10 Fig. 10. Effects of inclusions on temperature field on sectional views for 𝑦𝑦 = 9 cm at 𝜕𝜕 = 25 h. (left: temperature field, right: difference with a healthy pavement)

Fig. 11 .

 11 Fig. 11. Effects of inclusions on temperature field at the surface at 𝜕𝜕 = 25 h. (left: temperature field, right: difference with a healthy pavement)This same difference on a line at the surface for 𝑦𝑦 = 9 cm along the studied time range is plotted in figure12. The influence of phase change is visible when the difference of temperature reaches extrema. It can also be seen in this figure that only one phase change occurs for brine inclusion, when the air temperature is stabilized at about -10°C between 𝜕𝜕 = 61 h and 73 h.

Fig. 12 .

 12 Fig. 12. Difference of temperature between pavement with two inclusions and healthy pavement at the surface for 𝑦𝑦 = 9 cm as a function of 𝑥𝑥 position and time. MCP inclusion is located for 𝑥𝑥 between 0.123 and 0.186 m and brine inclusion between 0.311 and 0.375 m Computed temperatures are compared with measured at thermocouples. Differences above PCM and brine inclusions are shown in figure 13. Difference above healthy pavement is shown in figure 14. The absolute differenc e remains for the three cases below 1.5°C during most of the experience time.

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Difference between measured and computed temperatures above PCM inclusion (left) and brine inclusion (right)

Fig. 15 .

 15 Second empirical orthogonal function maps (EOF2) calculated for the first air temperature decay: Experimental results (left (a) no inclusion, left (b) 2 inclusions) Numerical simulations results (right (c) no inclusion, right (d) 2 inclusions)For the specimen with no inclusion (a) and (c) same trends are observed in the second EOF maps. Anyway, some undesirable convective effects are present in the experimental map (upper right).

Acknowledgments

The authors will like to thank the French Ministry of "Transition Écologique et Solidaire" for supporting part of this work under grant agreement DGITM N 17/389.

Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org)


PhD Jérémie Pourchez
email: pourchez@emse.fr

PhD Valérie Forest

E-cigarettes: from nicotine to cannabinoids, the French situation

published or not. The documents may come

The promotion of vaping to reduce tobacco-related harms is debated worldwide. E-cigarettes could play a role in public health if they are shown to be effective and safe cessation aids. In January 2018, two comments illustrated the substantial controversy regarding e-cigarettes in the UK [START_REF] Bareham | E-cigarettes: further flavours of controversy within the controversy[END_REF][START_REF] Flint | The irresponsible promotion of e-cigarettes and Swaptober[END_REF] . Despite this debate, vaping has become the most popular cessation aid in the UK; the British situation is therefore being carefully watched by public health bodies around the world.

By contrast, things are moving more slowly in France in terms of attitudes towards e-cigarettes.

In 2016, 34.5% of French people aged 15-75 smoked tobacco, while only 3.3% used ecigarettes [START_REF] Pasquereau | Tobacco and e-cigarette in France: levels of consumption according to the preliminary results from the 2016 Health Barometer[END_REF] . In February 2016, the French High Council for Public Health stated that vaping was "a tool for reducing the risks of smoking" and "a smoking cessation aid for smokers" [START_REF] Hcsp | The risks and benefits of electronic cigarettes for the general population[END_REF] .

However, it also stated that e-cigarettes "may act as a way to start smoking" and "run a risk of renormalising tobacco consumption" [START_REF] Hcsp | The risks and benefits of electronic cigarettes for the general population[END_REF] . In May 2017, the French National Cancer Institute concluded that "health risks are expected to be reduced for tobacco smokers who switch to ecigarettes". However, the health effects of prolonged use of e-cigarettes are not yet known, which is why their use is discouraged among non-smokers. In May 2017 the French National Public Health agency claimed that a high proportion of French vapers were either current smokers (48.0% daily, 41.2% occasionally) or former smokers (10.8%) [START_REF] Pasquereau | Tobacco and e-cigarette in France: levels of consumption according to the preliminary results from the 2016 Health Barometer[END_REF] , with a tendency to quit tobacco smoking or to reduce the number of cigarettes they smoked every day by at least half [START_REF] Pasquereau | Electronic cigarettes, quit attempts and smoking cessation: a 6-month follow-up[END_REF] .

While controversy about the role of vaping in smoking cessation continues to rage, a new debate is growing around e-cigarettes and cannabinoids. Recently, many French media outlets have published content denouncing the vaping of cannabis derivatives. In November 2017, "Le Parisien" newspaper headlined: "Cannabis: the worrying success of the electronic joint" [http://www.leparisien.fr/societe/cannabis-l-inquietant-succes-du-joint-electronique-28-11-2017-7419270.php]. This trend deserves further exploration. First, e-cigarettes are increasingly being used as a new device for inhaling  9 -tetrahydrocannabinol (THC), the primary psychoactive drug of the cannabis plant. In France the sale, transport and production of THCcontaining products for personal use remains prohibited [START_REF] Legifrance | de l'article R. 5132-86 du code de la santé publique pour le cannabis[END_REF] . Concomitantly, cannabidiol (CBD) vaping is becoming more common in France in the second half of 2017, with anecdotal reports that sales of CBD e-liquid have been increasing. As CBD is considered non-psychoactive, pure CBD-based e-liquids of synthetic origin, as well as CBD products with THC content lower than

Cannabinoids are increasingly used in France for recreational purposes (e.g. smoked cannabis or CBD vaping) and for medical applications (e.g. to relieve some symptoms of multiple sclerosis). Based on lessons learned from the controversy surrounding nicotine vaping, the risks and benefits of cannabinoid vaping should be discussed. Even if e-cigarettes can be considered as a promising harm reduction strategy of smoked cannabis, as for nicotine vaping the potential issues are: exposure to airborne toxicants, limited evidence to support the efficacy in smoking cessation, and a possible gateway to start smoking. There is an urgent need for a scientific debate to inform political and regulatory decisions, and the sale of cannabinoid-based e-liquids, health warnings and better consumer information.

Conflict of interest

The authors report no conflicts of interest.