Azam Marie-Hélène
email: marie-helene.azam@cerema.fr

Morille Benjamin

Bernard Jérémy

Musy Marjorie

Rodriguez Fabrice

A new urban soil model for SOLENE-microclimat: review, sensitivity analysis and validation on a car park

Keywords: Urban soil model, Heat transfer, Soil surface temperature, SOLENE-Microclimat, Urban Heat Island. Highlights

The main purpose of this study is to evaluate an urban soil model that will accurately reproduce the heat flux into urban soil, which has an influence on the urban heat island effect, for typical urban land use such as a car park. After a complete literature review, a sensitivity study is carried out on a large number of parameters: material properties, layer size, deep boundary condition, and convective heat transfer coefficient. The model's ability to reproduce heat conduction transfer is validated via a measurement campaign performed on an asphalt car park during hot days. The mean daily RMSE between estimated and observed surface temperature is 0.86 0 C, and 0.72 0 C, 0.58 0 C, 0.26 0 C and 0.13 0 C respectively at 5cm-, 10cm-, 34cm-and 50cm-depths. Performances obtained using different node distributions are

Introduction

In conditions of global warming, the development of cities must be carried out considering the urban heat island (UHI) phenomenon [START_REF] Oke | Boundary layer climates[END_REF] as a serious environmental issue. This phenomenon has several consequences on outside comfort and on building energy needs. In order to mitigate the UHI, it is necessary to identify its causes and to quantify the impact of mitigation solutions. Measurements campaigns are useful to evaluate the UHI, but linking it to the influence of modifications in urban form or urban planning choices is quite tricky. For this purpose, numerical simulation is a powerful tool. Several models under development simulate the UHI phenomenon and its consequences.

Different scales are considered, depending on the application intended: for example TEB [START_REF] Masson | A physically-based scheme for the urban energy budget in atmospheric models[END_REF] or ARPS-VUC [START_REF] Tavares | Modelling the impact of green infrastructures on local microclimate within an idealized homogeneous urban canopy[END_REF] are more suitable for city-scale applications while models like SOLENE-microclimat [START_REF] Bouyer | Modelisation et simulation des microclimats urbains-Etude de l'impact de l'amenagement urbain sur les consommations energetiques des batiments[END_REF], Envi-met [START_REF] Yang | Evaluation of a microclimate model for predicting the thermal behavior of different ground surfaces[END_REF] and EnviBatE [START_REF] Gros | Modélisation de la demande énergétique des bâtiments à l'échelle d'un quartier[END_REF] are more appropriate for the district scale.

For a given scale, each tool may have one specific feature among many others: EnviBatE [START_REF] Gros | Modélisation de la demande énergétique des bâtiments à l'échelle d'un quartier[END_REF] is designed to study the energy demand of a group of buildings, SOLENE-microclimat [START_REF] Bouyer | Modelisation et simulation des microclimats urbains-Etude de l'impact de l'amenagement urbain sur les consommations energetiques des batiments[END_REF] focuses on outdoor comfort and on the impact of urban climate on indoor comfort, and ENVImet [START_REF] Yang | Evaluation of a microclimate model for predicting the thermal behavior of different ground surfaces[END_REF] is dedicated to outdoor comfort.

All those models have in common the fact that they represent several physical mechanisms: radiative fluxes, thermal fluxes and fluid dynamics.

Furthermore, the representation of these phenomena is essential to accurately calculate the soil surface temperature, which is the key to the interaction between the soil and the urban environment (radiative and sensible fluxes).

The heat flux stored and released by urban material respectively during dayand night-time is one of the main causes of UHI development. This heat flux is greater in urban than in rural areas due to the high inertia of the materials used. The simulation of heat transfer in facades and in the soil are therefore of the highest importance.

Consequently, there is a need to enhance the ability of urban soil models to simulate heat fluxes in complex heterogeneous urban contexts. The SOLENE-microclimat model (whose efficiency has already been proven) is being developed and validated to answer this need. It is a complex model made up of several modules [START_REF] Musy | The use of solene-microclimat model to assess adaptation strategies at the district scale[END_REF], implying that validation should be performed on each module individually. This article presents one of these validation steps.

The main purpose of this article is to modify and validate the existing The open space is chosen to get away from the constraints of an urban environment: diffuse and reflected solar radiation depending on masks, the albedo of the surrounding surfaces, long-wave radiation exchanges with surrounding surfaces, etc.

As the accuracy of the model is often linked to its numerical cost, the additional aim of this study is to develop optimized discretization in order to reach a compromise between accuracy loss and computational efficiency.

The results of the study are divided into three parts:

• sensitivity study on model parameters,

• model validation using a centimetric grid in order to evaluate the model's ability to reproduce conductive heat transfer into the soil,

• calculation of the accuracy loss caused by different node distributions. The model's performance is then compared to that of models identified in the literature. Finally, performance is analysed for a whole year of simulation.

This exhaustive study provides accurate information about the reliability of the SOLENE-microclimat soil model.

State of the art

Existing models

In the literature on microclimate models, soil representation is rarely fully described. However, soil models are also used in other fields such as:

• Geothermal energy applications,

• Road applications: pavement sustainability or frost forecasts,

• Hydrology and interaction between soil, vegetation, and atmosphere.

Those other domains have the advantage of proposing a different point of view on the way to model heat transfer in the ground. In Table 1 the articles used for the following literature review are shown together with the characteristics of the soil models. Even though the field of application is different, all the articles presented have the common objective of predicting surface temperature or ground heat flux. Depending on the application, the physical mechanisms modelled are not the same. In addition to conductive heat flux, moisture flow is often modelled. This is the case for several applications that use bare-soil and vegetation covers [START_REF] Qin | Numerical solution of a complete surface energy balance model for simulation of heat fluxes and surface temperature under bare soil environment[END_REF][START_REF] Saito | Effects of meteorological models on the solution of the surface energy balance and soil temperature variations in bare soils[END_REF][START_REF] Herb | Ground surface temperature simulation for different land covers[END_REF][START_REF] Asaeda | The subsurface transport of heat and moisture and its effect on the environment: a numerical model[END_REF]. This is useful to estimate the water availability for vegetation or to adjust the thermal properties of the soil depending on humidity content. Nevertheless, for impervious surfaces, moisture flux is mostly neglected [START_REF] Herb | Ground surface temperature simulation for different land covers[END_REF].

The soil model presented in this paper is designed for an impervious surface in an urban environment. Urban grounds are heterogeneous and made up of different layers characterized by large differences in their physical properties. The model must therefore be able to take several layers into account. In the literature, soil is modelled by either a homogeneous or a heterogeneous column. But for a given area, the size of each layer and its physical properties are not accurately known. By simplification, half of the soil models presented here consider a homogeneous column of soil (Table 1).

Nevertheless, in order to accurately simulate the conduction flux all along the vertical axis, the soil profile should be consistent with reality, which implies considering a heterogeneous soil.

Thermal properties are either set to experimental data (i.e.: to better represent measurements the albedo is defined as the ratio between incident and reflected solar radiation), or calibrated.

However, the impact of thermal properties has rarely been investigated.

Best (1998, [START_REF] Best | A model to predict surface temperatures[END_REF]) and Herb et al. (2008, [13]) studied the influence of material characteristics (diffusivity, specific heat, etc.) on surface temperature.

According to Best (1998 [9]) and Herb et al. (2008 [13]), the emissivity and the thermal conductivity of the pavement have the most influence on surface temperature, while the characteristics of the soil underneath have little influence.

In the case of a homogeneous soil column, there are several possibilities for obtaining ground temperature variation. Among others, we may cite analytical solutions (Fourier analysis [START_REF] Best | Determining the optimal soil temperature scheme for atmospheric modelling applications[END_REF][START_REF] Chow | Estimation of soil temperature profile in hong kong from climatic variables[END_REF][START_REF] Ho | A soil thermal model for remote sensing[END_REF][START_REF] Jacovides | On the ground temperature profile for passive cooling applications in buildings[END_REF][START_REF] Mihalakakou | On the application of the energy balance equation to predict ground temperature profiles[END_REF][START_REF] Mihalakakou | On estimating soil surface temperature profiles[END_REF][START_REF] Ozgener | A practical approach to predict soil temperature variations for geothermal (ground) heat exchangers applications[END_REF][START_REF] Swaid | The prediction of impervious ground surface temperature by the surface thermal time constant (sttc) model[END_REF]), the empirical method [START_REF] Chow | Estimation of soil temperature profile in hong kong from climatic variables[END_REF][START_REF] Diefenderfer | Model to predict pavement temperature profile: development and validation[END_REF][START_REF] Mihalakakou | On estimating soil surface temperature profiles[END_REF], or even the force-restore method [START_REF] Lin | On the force-restore method for prediction of ground surface temperature[END_REF].

For a heterogeneous soil column, only the numerical method can be used to accurately estimate surface temperature and ground temperature for several depths. Twelve of the models described use the finite difference method with an implicit scheme, [START_REF] Asaeda | The subsurface transport of heat and moisture and its effect on the environment: a numerical model[END_REF][START_REF] Best | A model to predict surface temperatures[END_REF][START_REF] Best | Determining the optimal soil temperature scheme for atmospheric modelling applications[END_REF][START_REF] Bouyer | Modelisation et simulation des microclimats urbains-Etude de l'impact de l'amenagement urbain sur les consommations energetiques des batiments[END_REF][START_REF] Gros | Modélisation de la demande énergétique des bâtiments à l'échelle d'un quartier[END_REF][START_REF] Herb | Ground surface temperature simulation for different land covers[END_REF][START_REF] Hermansson | Mathematical model for calculation of pavement temperatures: comparison of calculated and measured temperatures[END_REF][START_REF] Ho | A soil thermal model for remote sensing[END_REF][START_REF] Masson | A physically-based scheme for the urban energy budget in atmospheric models[END_REF][START_REF] Qin | Numerical solution of a complete surface energy balance model for simulation of heat fluxes and surface temperature under bare soil environment[END_REF][START_REF] Saito | Effects of meteorological models on the solution of the surface energy balance and soil temperature variations in bare soils[END_REF][START_REF] Yang | Evaluation of a microclimate model for predicting the thermal behavior of different ground surfaces[END_REF] except Nowamooz et al. (2015, [20]) who use an explicit scheme. The soil model presented in this paper is based on an implicit finite difference method. The following section presents in detail all possible parametrizations in the case of the finite difference method.

Parametrization of finite difference models

For a given problem-solving method, different choices can be made regarding the node distribution (discretization), boundaries and initial conditions (Table 2).

Vertical discretization

Depending on the author, different vertical discretizations are used. Ho (1987, [START_REF] Ho | A soil thermal model for remote sensing[END_REF]) and Qin et al. (2002, [22]) propose uniform layout nodes, whereas most of the articles present a denser node distribution near the surface than deeper in the ground. Saito and Simunek (2009, [START_REF] Saito | Effects of meteorological models on the solution of the surface energy balance and soil temperature variations in bare soils[END_REF]) use a denser distribution near the surface and also near the interface between two layers. However, the choice of node distribution is rarely clearly justified. Only Best et al. (2005 [10]) studied the behaviour and accuracy of the model for different vertical discretizations, comparing the numerical solution with the analytical one. The present study also aims to perform a sensitivity analysis on the discretization in order to know how much accuracy is lost when reducing the number of nodes. Few sensitivity studies have been performed on this point.

Asaeda and Ca (1993, [START_REF] Asaeda | The subsurface transport of heat and moisture and its effect on the environment: a numerical model[END_REF]) studied the influence of grid size. As expected, the thinner the resolution of the grid, the more precise the surface temperature. But this accuracy gain is obtained to the detriment of calculation duration. So as Best et al. (2005 [10]) suggested, a compromise should be found between accuracy and execution time.

Boundary conditions

To solve one dimension heat conduction equation with the finite difference method, two boundary conditions are required. The bottom condition taken at a defined depth can either be a zero flux or a constant temperature (Table 2). If the simulation concerns a short period and the last node of the grid is below the diurnal amortization depth, the type of condition has little influence on the surface temperature [START_REF] Best | A model to predict surface temperatures[END_REF].

As the purpose of most articles is to predict surface temperature, the upper boundary condition is defined by the surface energy balance. This balance is composed of radiative fluxes, the sensible flux, the latent flux, and the conductive heat flux into the soil. The long and short-wave radiation fluxes are always calculated, except for the validation process for which measurements are used whenever they are available.

As urban soils have a low albedo, the temperature gradient between the surface and the air can be very great, especially during clear days. This phenomenon leads to natural convection (driven by buoyancy forces). Nevertheless, in the surface energy balance, most authors only take into account the forced convective mode (driven by wind forces), whereas Herb et al.

(2008, [START_REF] Herb | Ground surface temperature simulation for different land covers[END_REF]) proposes considering both forced and natural modes.

Both Hermansson (2004, [14]) and Herb et al. (2008, [13]) assumed that the method used to calculate the convective heat transfer coefficient is a critical point. To test this assumption, they both carried out a sensitivity study on the coefficients of an empirical formula. Herb et al. (2008, [13]) show that the modification of convection parameters can sometimes have a significant impact on surface temperature: for a 10% increase in the convection parameters, the average surface temperature variation is about 0.24 0 C. This might be explained by low wind speed and by a high-temperature gradient between the surface and the air. Hermansson (2004, [14]) dissociates winter and summer periods since the temperature gradient between the surface and the air is smaller during winter. He proposes to use two sets of parameters regarding solar radiation conditions. Depending on the season, convection losses and wind velocity are balanced by coefficients.

Initial conditions

In general, little information is given on the initial conditions (temperature profile) except for the validation process. Most authors use experimental values (Table 2). When this information is not available, several alternatives are proposed: a constant temperature profile [START_REF] Bouyer | Modelisation et simulation des microclimats urbains-Etude de l'impact de l'amenagement urbain sur les consommations energetiques des batiments[END_REF][START_REF] Ho | A soil thermal model for remote sensing[END_REF] or an exponential profile between the deep soil and the surface is set [START_REF] Best | A model to predict surface temperatures[END_REF][START_REF] Ho | A soil thermal model for remote sensing[END_REF]. However, according to (1987, [START_REF] Ho | A soil thermal model for remote sensing[END_REF]), this parameter has little impact on surface temperature

Methodology of the study

Description of the proposed soil model

The soil model presented in this paper is designed for an impervious surface such as a pavement coating. Only heat transfer is therefore taken into account (moisture transfer is neglected). The soil model is 1-D, defined as a one-dimensional soil column where each layer has its own characteristics.

In an unsteady state, temperature fluctuation is calculated from Equation 1, which is an application of the heat equation for a one-dimensional problem.

BT Bt " α soil . B 2 T Bx 2 (1)
α soil : thermal diffusivity of the soil [m 2 s ´1]
The problem is solved by a finite difference method using an electrical analogy. Heat resistances represent the resistance to heat transfer through a ground layer and heat capacities, the heat storage capacity of a ground layer.

They are defined in Figure 1.

The soil model is composed of n nodes. The energy balance equation is calculated at the surface (node i = 0, Equation 2), and then for each following node i Ps0 : n ´1s (Equation 3) until the last one which includes the bottom boundary condition (i " n, Equation 4).

T surf ace ´Tair R c `Tsurface ´T1 R 1 `Ce dT surf ace dt " R net ´LE (2)
T i ´Ti´1 R i `Ti ´Ti`1 R i`1 `Ci dT i dt " 0 (3)
T n ´Tn´1 R n´1 `Tn ´T8 R n `Ci dT n dt " 0 (4) R net : net radiation [W.m ´2] LE: latent heat flux [W.m ´2] H: sensible heat flux [W.m ´2]
C e : heat capacity of the surface layer [J.m ´2.K ´1]

C i : capacity of the layer at the node i

[J.m ´2.K ´1] R c : convection resistance [m 2 .K.W ´1] R i : heat resistance of the layer between the node i ´1 and i [m 2 .K.W ´1]
T surf ace : surface temperature [K]

T air : air temperature [K]

T i : temperature of the node i [K] According to this method, any node distribution and any boundary condition depth may be used. Most of the authors [START_REF] Asaeda | The subsurface transport of heat and moisture and its effect on the environment: a numerical model[END_REF][START_REF] Best | A model to predict surface temperatures[END_REF][START_REF] Best | Determining the optimal soil temperature scheme for atmospheric modelling applications[END_REF][START_REF] Yang | Evaluation of a microclimate model for predicting the thermal behavior of different ground surfaces[END_REF] in the literature work with a centimetric grid when accuracy is required. A model with one node per centimetre is used (see Section 3.1.3).

The top node is located at the ground surface and the deepest one at a depth of 1m. At this depth, the temperature is supposed to be constant over a day. A new value will be set for each day. More details are given in Section 3.1.1.

Figure 1 illustrates the way thermal properties are defined. Each layer of material is assumed to be homogeneous and isotropic and has its own characteristics which are considered to be constant over time.

As the node distribution is defined aside from the soil layers and characteristics, the nodes are not automatically located at the interface between two layers.

Deep boundary condition

In deep soil, the temperature is taken to be constant over a day. In the case of homogeneous soil, an analytic solution can be used to calculate the temperature for any depth z and any time step t. If the surface temperature is considered to be sinusoidal, the analytic solution follows Equation 5. The parameters T ma , A a and t 0 are respectively the mean, the amplitude, and the phase of a day surface temperature signal.

T pz, tq " T ma `Aa . expp´z zd a q sinpw a pt ´t0 q ´z zd a q (5)

T ma : mean annual temperature [ 0 C]
w a : annual beat w a " 2 ˚π{31536000 [rad.s ´1]
t 0 : day of the year where the surface temperature was the coldest From a certain depth, the daily signal is completely damped. The damping depth depends on the soil characteristics through the parameter z da , defined by Equation 6. The depth from where the signal is damped at 95 or 99 % can be estimated from Equation 7.

z da " c 2.α sol w a (6)
Apzq " A j . expp ´z z da q (7)

A j : half daily amplitude of the climatic thermal wave at the surface [ 0 C] For a range of materials (asphalt, concrete, bare-soil), the most diffusive is marble. For this material, the depth corresponding to a daily damping of 99% is 0.89cm. Beyond a meter in depth, the temperature is assumed to be constant throughout the day whatever the type of ground. For this reason, the bottom node is located below this depth.

Upper boundary: Heat flux across a ground surface

The upper boundary condition is defined by the energy balance at the ground surface (Equation 8).

R net " Q cond `H `LE (8)
with R net , Q cond , H, LE previously defined.

Radiative flux

The net radiative flux is the balance of all radiative fluxes at the soil surface. It is the sum of short-wave radiation and long-wave radiation. The historical SOLENE radiative model computes radiative transfers, including long-wave radiation, inter-reflexion and shading effects [START_REF] Musy | The use of solene-microclimat model to assess adaptation strategies at the district scale[END_REF]. So this input data does not need to be calculated in future simulations.

Convective heat flux

The heat flux exchanged between the surface and a moving fluid can be expressed with the Equation 9:

H " h c pT air ´Tsurface q (9)

h c : convective heat transfer coefficient [W.m ´2.K ´1]
T air : air temperature [K]

T surf ace : surface temperature [K]

In order to calculate this flux, the convective heat transfer coefficient is required. In the literature, this coefficient is always a function of wind speed.

Linear or power law functions are used, or correlations using dimensionless numbers ( [START_REF] Palyvos | A survey of wind convection coefficient correlations for building envelope energy systems' modeling[END_REF], [START_REF] Defraeye | Convective heat transfer coefficients for exterior building surfaces: Existing correlations and cfd modelling[END_REF]). For urban applications, the first solution is often used.

The simplest is a linear relation of the wind speed (Equation 9). Depending on the situation (i.e. surface texture, wind velocity, windward/leeward surface, etc.), Palyvos (2008, [START_REF] Palyvos | A survey of wind convection coefficient correlations for building envelope energy systems' modeling[END_REF]) suggests around forty combinations for a and b coefficients. For a horizontal surface, with low winds (V air ă 5m{s), several coefficients are proposed (Table 3).

h c " a `b ˚Vair (10)
N u " epaRe b `cGr d q f (11) N u " h c .L c λ f luid (12)
Re: Reynolds number

Gr: Grashof number N u: Nusselt number L c : characteristic length [m] λ f luid : thermal conductivity of the fluid [W.m ´1.K ]
This type of correlation can be used to define a convective heat transfer coefficient as a function of the flow modes:

• Natural or free convection: air flow driven by buoyancy forces,

• Forced convection: air flow driven by wind forces,

• Mixed convection: when the air flow is created by both wind and buoyancy forces.

In the case of urban applications, wind speed is often low and the temperature gradient great. Free convection may then become predominant.

A comparison between the MacAdams formula and the correlation method is made in order to choose the most suitable method for an urban climate application (Section 4.1.1).

Node distribution

Authors in the literature work with centimetric grid when accuracy is required ( [START_REF] Asaeda | The subsurface transport of heat and moisture and its effect on the environment: a numerical model[END_REF][START_REF] Best | A model to predict surface temperatures[END_REF][START_REF] Best | Determining the optimal soil temperature scheme for atmospheric modelling applications[END_REF][START_REF] Yang | Evaluation of a microclimate model for predicting the thermal behavior of different ground surfaces[END_REF]). In order to gauge the accuracy of the proposed grid different grid sizes are compared to the analytic solution (Equation 5) for a homogeneous soil of asphalt-concrete (α " 1.04 ˚10 ´6m 2 .s ´1).

The analytic solution is compared to a finite difference model with a sinusoidal temperature signal imposed at the surface (Equation 13). The bottom boundary condition is taken as a fixed temperature equal to the mean annual temperature. The model is initialized with a temperature profile calculated by the analytic solution. The simulation is run for 6 days with a time step of 900 seconds.

T s ptq " T ma `Aa . sinpw a pt ´t0 qq

Different grid sizes are tested from 0.0025m to 0.025m with steps of 0.0025m. For each grid size, the RMSE is calculated at the last time step, between the analytic solution and the finite difference model.

Presentation of the measurement campaign

The measurement campaign provided data that are used either as model inputs or as validation data to evaluate the performance of the model and its accuracy. This campaign especially focused on surface and air temperatures and on heat flux measurements during a warm summer period. Artificial rain events were created during the campaign but the present study focuses only on dry weather periods.

The experiment site is located near Nantes (France) within the IFSTTAR center of Bouguenais, and consisted of a 2, 500m 2 asphalt car park (Figure 3(b)). The soil structure is composed of a 5cm-asphalt layer, a 45cm-ballast layer and an altered mica-schist natural soil underneath. Throughout the observations available for this campaign, this study focuses on the following variables, all observed in the middle of the car park described in the Table 5 and in Figure 3(a):

• surface and ground temperatures: vertical profile at depths of 0, 1, 10, 3, 4, 5, 6, 10, 15, 104, 34, 50 and 75 cm;

• wind speed and direction;

• humidity and air temperature;

• convective heat fluxes;

• radiation components.

Physical quantity

Instrument type Height The meteorological conditions are presented in Figure 4 and Figure 5 .

Soil model setup

In order to model the car park described above, the model parameters need to be set for the current situation.

Deep boundary condition

The deep boundary condition is calculated from Equation 5. The value for each of the parameters (T ma , A a and t 0 ) is set according to air temperature signals. The values used, derived from measurements recorded over 4 years in three locations of the city of Nantes (France). The mean annual temperature is T ma " 2.50 0 C , the yearly half amplitude of the daily mean temperature is A a " 1.86 0 C and the phase shift t 0 " 0 days.

Upper boundary condition

The measurement campaign provided data that are used as model inputs for the calculation of the surface energy balance on the upper boundary condition (Equation 8).

The sensible heat flux is calculated from the measured air temperature.

For the radiative budget measurements are used. Latent heat flux is not considered for now, as in this step of the model only dry days are modeled.

In Equation 8, LE " 0.

Calibration of the material' properties

As seen before, the thermal characteristics of the soil are often unknown From the measured temperature gradient, changes in the soil thermal properties within the first layer were identified (0-1cm, 1-5cm). The total capacity of the first layer was then estimated with the surface energy budget measured, the ground heat flux measured at 3.5cm and the temperature gradient. The thermal characteristics were then adjusted by an iterative procedure reducing the difference between the measured and the simulated temperature at the surface and various depths. The thermal characteristics values were considered acceptable when the RMSE on each ground temperature had the same level of magnitude as the uncertainty of the sensor (0.2 0 C).

The calibrated material characteristics are summarized in Table 6.

As the material properties were calibrated with the calculated convective heat flux, the proposed thermal properties compensate for the error between the calculation and measurement of the convective heat flux. After calibration, the model correctly reproduces the conducted heat flux into the different layers of the soil as shown in Figure 6. In this article, the indicator chosen is the RMSE (Root Mean Square Error). As the authors used this, it will be easier to compare model accuracy.

In order to determine whether the dynamics of heat storage are properly reproduced, the RMSE will be calculated at the surface and at several depths. For grid optimization, the choice was made to work with annual simulation. First, an analytic solution is used to better understand ground temperature dynamics. In this way, two sinusoidal temperature signals are applied to the surface: a three-hour period (corresponding to weather change or shadows created by buildings during a day) and a day period (corresponding to the day-night cycle). The profiles obtained from the analytical solution are presented on the right-hand side of Figure 7. The closer we get to the deep condition, the more linear the profile. With the finite difference method, the more linear the profile, the less dense the distribution of nodes needs to be.

According to this result, node distribution is different in the three zones: the criteria used to identify each of them are given in Table 7.

Three node distributions are then proposed to better represent these ground temperature profiles (Figure 7). The final distribution profiles are consistent with those proposed in the literature: the node density is higher in the upper layer of the ground than in the lower layers. 4. Results

Sensitivity study

Some assumptions made for the model parametrization may affect the surface temperature calculation. Some authors carried out a sensitivity study on the surface and material parameters, on the convective heat transfer coefficient or on the grid size. However, none of these studies has compared the relative influence of all the affecting parameters.

In the following section, a model sensitivity study is performed on the surface temperature on June 5 th and 6 th . Two kinds of parameters are used for this study: those used for the calibration step (soil thermal characteristics and the convective heat transfer coefficient) and two additional parameters which have been identified to be relevant: the deep boundary conditions and the size of the layers. Parameter sensitivity is studied regarding the order of magnitude of its uncertainty.

The sensitivity study is presented from the most influential parameter to the least influential: the convective heat transfer coefficient, followed by the material characteristics, the layer size, and the deep boundary condition.

Convective heat transfer coefficient

Several methods are presented to estimate the convective heat transfer coefficient h in Section 3.1.2. Two methods are first compared to select the most suitable:

• MacAdams (1954, [START_REF] Mcadams | Heat transmission, 3rd[END_REF]) formula: linear function of the wind speed (coefficients a and b equal to 3.8 and 5.7 respectively).

• Correlation equation with dimensionless numbers and different characteristics length (1, 10, 50 m).

Figure 8 illustrates the comparison. For the correlation method, after calculating the dimensionless coefficients, the convection mode varies over time:

it is mostly mixed during the night and forced during the day. Free convection is not represented by the first formula, whereas it is not integrated into the MacAdams formula. The correlation method is closer to the measured heat flux throughout the comparison period. The correlation method used to calculate the convective heat transfer coefficient (Section 3.1.2), is based on results obtained by Tain and Petit (1989, [START_REF] Taine | Transferts thermiques: mécanique des fluides anisothermes: cours et données de base[END_REF]) applied for a horizontal flat plate. In order to apply the similitude theory, a characteristic length should be set which is most of the time defined as the distance from the leading edge. Applied to our case, this characteristic length is difficult to define. The influence of several lengths (1, 10 and 50 m) is tested (Figure 8). Among the tree values tested, it appears that the 1m-characteristic length correlation method produces the lowest RMSE (24.27W.m ´2 ).

For each value, the correlation method remains better than the MacAdams method. As a result, for the following sub-sections, the correlation method is chosen with a characteristic length of 1m. The resulting convective heat transfer coefficient varies from 2.6 to 12.6 W.m ´2.K ´1.

In order to see the influence of the characteristic length value on the sensible heat flux, the RMSE calculated with the 1m-characteristic length is compared to the daily mean sensible heat flux ( 61W.m ´2). The RMSE represents 40% of the mean experimental value.

According to these levels of magnitude, we performed the sensitivity study varying the convective heat transfer coefficient up to 40% of its initial value.

The aim is here to quantify the influence of this coefficient on the surface temperature. The results are presented in Table 8.

Indicator h `40% h ´40%

Maximum error ( 0 C) 3.70 5.46

Mean error ( 0 C) 1.89 2.78 RMSE ( 0 C) 2.12 3.14

Sensitivity of the layer definition

Modification of soil thermal conductivity, soil density and the thickness of the layers is performed one by one for each layer. The magnitude of the modifications and the layers concerned by the modification are described in Table 9. The surface temperature modification caused by each material property change is also presented in this table.

Parameter

Modified layer Maximum difference at the surface ( 0 C)

Thermal conductivity 0 ă 0, 1 λ `´10% 1 ă 0, 7 2 ă 0, 7 3 ă 0, 1 Density 0 ă 0, 07 ρ `´5% 1 ă 0, 3 2 ă 0, 3 3 ă 0, 07
Size of the layer 1 ă 0, 04 e `´0, 01m 2 ă 0, 15 The height of each soil layer is not constant over the depth. Each height is therefore roughly estimated. To evaluate this lack of accuracy, the influence of the biggest layers (1 and 2) is investigated (Table 9).

The temperature change associated with layer size modification has the same magnitude as that associated with soil property modification of layers 0 and 3. However, it is negligible compared to the temperature change associated with soil property modification of layers 1 and 2.

Deep boundary condition

The temperature imposed as a deep boundary condition may have an influence on the surface temperature. If the deep temperature is overestimated by one degree, the surface temperature increases by only 0.05 0 C.

Model ability to reproduce heat conduction transfer: validation

The ability of the model to properly reproduce the physical phenomenon is evaluated in this section. In order to quantify the uncertainty due to the model itself (i.e. the accuracy of the physical phenomenon representation), the model is first evaluated comparing temperature estimation to experimental data. For this purpose, the model with a centimetric grid is used. As was noticed in Section 3. To quantify the accuracy of the model in reproducing the temperature variation at several depths, Qin (2002, [START_REF] Qin | Numerical solution of a complete surface energy balance model for simulation of heat fluxes and surface temperature under bare soil environment[END_REF]) divide the RMSE by the amplitude of the signal. In fact, as the amplitude decreased, the relative error increased. Going deeper in the ground the amplitude decreased with the RMSE but in proportion this error increased. In fact, at the surface, the RMSE represents only between 10% to 4.7% of the amplitude, while at a depth of 34 cm and 50cm it represents between 6% to 11.7% and between 8.9% to 104.7% respectively. The error remains under the uncertainty of the temperature measurement.

Influence of node distribution

The ability of the model to reproduce the physical phenomenon (centimetric grid model compared to the experimental data) is compared to the accuracy loss due to the reduction of the number of nodes. For the 8/10/4 node model, the reduction in the number of nodes has a very negligible influence. Even though the total error of this model with experimental data (0.62 0 C or 1.30 0 C) increases, this demonstrates that a higher number of nodes might not be so relevant, since the uncertainty does not make it possible to say whether model accuracy would be higher or not.

In all cases, the number of nodes is not the main cause of the error. and minimum peak are slightly underestimated.

The underestimated daily maximum peak and minimum trough are due to the reduction in the number of nodes that worsens the representation of the heat transfer into the ground. In fact, this induces a time shift of heat conduction, the main influence of which appears when its sign changes.

For the 4/6/2 and for the 4/3/2 node model, the accuracy loss due to the reduction in the number of nodes remains lower than the model accuracy (centimetric grid RMSE). The absolute error of the model is sometimes lower using these node distributions than the centimetric distribution. The reason is that the model itself overestimates the surface temperature whereas as a lower number of nodes underestimates it (Figures 11 and12 ).

Finally, the reduction of the number of nodes has consequences on the reproduction of the daily peaks. This confirms that the correct representation Simulations are run during the summer period during sunny days with scattered clouds. The surface temperature rises from 41 0 C to 70 0 C, except for Best (1998, [START_REF] Best | A model to predict surface temperatures[END_REF]), who worked with a lower temperature.

According to the results presented in In general, one needs to exercise care in considering the results of this comparison for several other reasons:

• the surface temperature used as the reference temperature for each study comes from surface temperature measurement. Uncertainty on such measurement might be of the same magnitude as the RMSE observed, which makes comparisons between models difficult.

• the comparison of the model's performance is given for studies in similar conditions. Nevertheless, they cannot be exactly the same. The location is not the same, the weather conditions differ, and the surface type or at least the soil composition is not exactly the same.

Even though we conclude that the models proposed in this paper seem to present better performance than those of the literature, the difference must be put into perspective with application conditions. All those models should be applied to a single case study.

Performance according to meteorological data

The accuracy of the different grids according to meteorological data is evaluated comparing numerical profile of an annual simulation (hourly time step) for each node distribution (8/10/4, 4/6/2, 4/3/2). For those simulations, meteorological data recorded during one year (2010) at the Pin Sec station in the city of Nantes are used. These data were collected by the ONEVU (Observatoire Nantais des EnVironnements Urbains -urban environment observatory of the IRSTV [START_REF] Mestayer | The experimental campaign fluxsap 2010: climatological measurements over a heterogeneous urban area[END_REF]). Among the available observations, the following data are used as input in the model: air temperature, pressure and humidity, global and IR radiative flux, and wind velocity and direction. Table 13 presents the results of the comparison between the surface temperature model for each of the node distribution compared to the centimetric grid model. In the case of the 8/10/4 nodes model, the maximum absolute

Conclusion

The main purpose of this study is to propose a model that accurately reproduces the heat storage flux into urban ground as well as changes in surface temperature.

Some failings were identified in the literature review, pointing up the need:

• to perform an overall sensitivity analysis,

• to investigate the advantage of using different convection flow modes (forced, mixed, natural)

• to justify the chosen mesh,

• to assess models not only for surface temperature calculation.

To test the robustness of our model regarding its parameters, a complete sensitivity study was carried out on the model parameters: all have a negligible impact except the convective heat transfer coefficient.

Special attention has been paid to how the convective heat transfer coefficient is calculated. The chosen method makes it possible to take into account the different kinds of convection flow modes which is necessary for an urban application. Nevertheless, the sensitive heat flux can vary by 40%, leading to an uncertainty up to 3.14 0 C on the surface temperature RMSE.

After calibration of the different parameters, soil model (using a centimetric grid) accuracy is evaluated according to a measurement campaign carried out in a large asphalt car park during three sunny days with scattered clouds.

The RMSE between estimated and observed temperatures is calculated for surface temperature and ground temperature at several depths. Surface temperature RMSE is 0.86 0 C; RMSE for temperature at 34cm of depth is 0.26 0 C.

The results validate the ability of the model to reproduce heat storage into the ground.

Three node distributions are proposed on the basis of the analysis of the different temperature profiles throughout the depth. They all are dedicated to any kind of impervious urban surfaces. However, for less diffusive soils, the user is free to define other distributions in order to adapt the model to his application.The accuracy of the models varies from 0.62 0 C to 1.30 0 C. These performances are better than those of models from the literature applied under quite similar conditions.

Finally, application of the models over a whole year shows that only a few days are represented with an accuracy worse than 0.

  SOLENE-microclimat soil model. The modifications are based on a review of the literature and on a sensitivity analysis. Validation is performed using temperature and flux measurement carried out on an open car park.

Figure 1 :

 1 Figure 1: Schematic representation of the soil model: representation of node distribution, heat resistances and capacities, description of a common cross-section of an urban soil column with diffusive material layers in shade of grey and underneath natural soil in brown (1 column, print in color)

A

  a : annual half amplitude of the climatic thermal wave at the surface [ 0 C] zd a : damping depth with an annual beat [m]

V 9

 9 air : wind speed [m.s ´1] Cristofari et al. (2006) 5.67 3.86

Figure 2 :Figure 2

 22 Figure 2: Changes in RMSE for different grid sizes(1 column, print in color)

( a )Figure 3 :

 a3 Figure 3: Illustration of the ROSURE measurement campaign

Table 5 :

 5 Captor used and the corresponding measuresThe data were collected with a 1 min time step except for the sonic anemometer (0.1 s). The final data were averaged to 15 min time steps. From the whole measurement period, only days without artificial precipitation were used. Four days were selected with different meteorological conditions: three sunny days with scattered clouds (June 5 th , 12 th , 13 th ) and one clear sky day (6 th of June). The clear sky day was used to calibrate the model (6 th of June) while the other days were used to evaluate the model (June 5 th , 12 th , 13 th ).

Figure 4 :Figure 5 :

 45 Figure 4: Net Radiation measured at a height of 1m and Air Temperature measured at a height of 2m from June 5 th to 14 th (2 columns, print in color)

  and must be adjusted to properly represent fluxes and temperatures variations. The characteristics of the soil layers are calibrated according to the observed soil profile, reducing the difference between the measured and simulated surface temperature, with the centimetric grid. Data acquired on June 6 th are used for calibration.The temperature profile is initialized from ground temperatures measured on June 6 th at midnight. The deep boundary condition, corresponding to the ground temperature at a depth of 75 cm is set according to experimental data. At the surface, the boundary condition is calculated from the energy balance with the measured radiative heat flux and calculated convective heat flux (with correlation method formula and a characteristic length of 1m, see: Section 4.1.1 and Figure8).Albedo and emissivity of the surface are calculated from short and long wave radiations measured during the period of interest. For the albedo, the mean diurnal value is 0.173 (for reflected short-wave radiation flux K up ą 20W.m ´2), and the mean emissivity value over the period is 0.965. Soil composition and thermal properties are not known with any accuracy. Their values were not measured during the campaign. As a result, only approximate soil composition was available.Cohard et al. (2017, [29]) tried to estimate the first layers characteristics from the temperature and the ground heat flux measurement. They obtained the following thermal diffusivity α " 1.04˚10 ´6m 2 .s ´1 for a fixed capacity (with ρ " 920kg.m ´3 and Cp " 2050J.kg ´1.K ´1). The resulting thermal conductivity is 1.96W.m ´1.K ´1.

  Yang et al. (2013,[START_REF] Yang | Evaluation of a microclimate model for predicting the thermal behavior of different ground surfaces[END_REF]) noticed that thermal characteristics of the asphalt layer may vary with depth. Due to asphalt compaction, the layer density and the asphalt proportion are not constant along the depth, altering the thermal properties.Yang et al. (2013, [5]) chose to divide the asphalt layer into several layers to which different properties were attributed. The same phenomenon is observed through analysis of the temperature signal within the asphalt layer. For this reason, the asphalt layer is divided into two layers (of 1cm and 4cm respectively).

Figure 6 :

 6 Figure 6: Comparison of simulated and measured temperatures at the surface and at several depths on June 6 th (1 column, print in color)

3. 5 .

 5 Node distribution definition methodology Optimization of the node distribution represents a major challenge. The principle consists in reducing the number of nodes until a situation is reached where calculation time and lack of accuracy are minimal. of the daily signal for the most common material for this part ( soil, concrete, stone) 3 1.0 Amortization at 95% of the daily signal for the material with the highest diffusivity

Figure 7 :

 7 Figure 7: Three proposed node distributions (On the left-hand side the three node distributions, and on the right-hand side the temperature profiles obtained from the analytic solution.) (2 columns, print in color)

Figure 8 :

 8 Figure 8: Comparison between the measured sensible heat fluxes at a height of 2m and two calculation methods: the correlation method and the Mac Adams formula. (2 columns, print in color)

4 ,

 4 the evaluation of soil model performance based only on surface temperature comparison is one of the lacks identified in the literature. In this study, results are compared on the basis of the temperature at the surface and at several depths in the soil.The model was evaluated during June 5 th , 12 th , and 13 th . The model correctly reproduces the heat flux conduced into the different layers of the soil as shown in Figures9 and 10. The maximum difference between the measured and simulated temperature at the surface is 1.67 0 C, with an RMSE of 0.71 0 C for June 5 th ; 1.16 0 C with an RMSE of 0.60 0 C for June 12 th and 2.45 0 C with an RMSE of 1.28 0 C for June 13 th . The mean daily RMSE between the estimated and the observed surface temperature is 0.86 0 C, and 0.72 0 C, 0.58 0 C, 0.26 0 C and 0.13 0 C at 5cm-, 10cm-, 34cm-and 50cm-depths respectively (June 5 th ,12 th and 13 th ).

Figure 9 :Figure 10 :

 910 Figure 9: Comparison of simulated and measured temperatures at the surface and at several depths, June 5 th (1 column, print in color)

Figure 11 and

 11 Figure 11 and 12 compare the surface temperature changes of the different models with that measured. Firstly, all the models properly represent the change of surface temperature with time. Nevertheless, the daily maximum

Figure 11 :Figure 12 :

 1112 Figure 11: Surface temperature calculated with the different grids compared to that measured (June 5 th ).(1 column, print in color)

  Two of the authors highlight that the high amplitude of surface temperature recorded during the day is harder to represent than the surface temperature during the night.Yang et al.(2013, [5]) described the fact that the model fits well with the measurement at night time and in the morning, but when the temperature rises during the afternoon the model underestimates the temperature with an average o3.5 0 C.Best (1998, [9]) also has better results during the night (RMSE 0.83 0 C) than during the day (RMSE 1.26 0 C). This analysis is consistent with the assumption made previously: clear, hot days, characterized by high surface temperature amplitude, are the most difficult days to simulate.

  5 0 C. Most days are even reproduced with an accuracy better than 0.2 0 C. The investigation demonstrates that surface temperatures during clear and sunny days are the most difficult to reproduce. This paper provides a complete overview of soil model performance, comparing it with experimental data, centimetric grid model, and literature results. Because of measurement uncertainty, better performances would be difficult to obtain and especially to assess. Nevertheless, comparison of the model performances with results from the literature would require a benchmark. The model presented, now validated, is ready to include moisture transfer and in particular the evaporation of water at the soil surface, to help to properly assess the effect on local climate and outdoor comfort of moistening techniques.

Table 1 :

 1 Model application and type of soil

	Article	Type of application B, GE Roads SVA	UM	Surface Previous	Imprevious	Type of coating Pavement Bare-soil	Vegetation	Soil column Homogeneous	Heterogeneous
	Asaeda and Ca (1993,[8])		x	x		x	x			x
	Best (1998, [9])	x		x		x	x	x	x
	Best and al. (2005, [10])	x		x			x		x
	Bouyer (2009, [4])		x	x		x	x	x		x
	Chow and al. (2011, [11])	x		x				x	irrelevant
	Diefenderfer and al. (2006, [12])	x			x	x			irrelevant
	Gros et al. (2015, [6])		x		x	x				x
	Herb et al. (2008, [13])	x		x		x	x	x		x
	Hermansson (2004, [14])	x			x	x				x
	Ho (1987, [15])	x			x	x	x		x
	Jacovides (1996, [16])	x		x			x	x	x
	Lin (1980, [17])	x		x			x		x
	Masson (2000, [2])		x		x	x				x
	Milhalakakou and al. (1997, [18])	x		x			x	x	x
	Milhalakakou and al. (2002, [19])	x		x			x	x	x
	Nowamooz and al. (2015, [20])	x		x				x		x
	Ozgener and al. (2013, [21])	x		x			x		x
	Qin and al. (2002, [22])	x		x			x			x
	Saito and Simunek (2009, [23])	x		x			x			x
	Swaid and Hoffman (1989, [24])	x			x	x				x
	Yang and al. (2013, [5])		x	x		x	x	x		x

B, GE: Buildings, geothermal energy SVA: Interaction soil, vegetation, atmosphere UM: Urban Microclimate

Table 2 :

 2 Parametrization of finite difference models

	EV: Experimental values NI: No Information CFD: Computational Fluid Dynamics

Table 3 :

 3

	Convection mode Flow regime a	b	c	d	e	f
	Free	laminar	0	0	0.49 1/4 1	1
		turbulent	0	0	0.13 1/3 1	1
	Mixed	laminar	1	3/2 0,57 3/5 0,68 1/3
		turbulent	1	12/5 12,1 1	0,03 1/3
	Forced	laminar	0,56 1/2 0	0	1	1
		turbulent	0,03 4/5 0	0	1	1

Coefficients a and b for a flat surface low wind speed

Methods based on correlations that use dimensionless numbers (Reynolds, Grashof, and Nusselt) also exist. For flat surfaces Morille (2012,

[START_REF] Morille | Élaboration d'un modèle du climat distribué à l'échelle de l'abri et de la plante en cultures ornementales sous serres: analyse des transferts de masse et de chaleur, bilans énergétiques[END_REF]

) presents coefficients depending on the flow regime (Table

4

).

Table 4 :

 4 Coefficients as a function of the convection mode for a flat surface

Table 6 :

 6 Calibrated characteristics of the soil

	For this reason

2.88W.m ´1.K ´1, diffusivity ranges from 4.4 ˚10 ´7 to 14.4 ˚10 ´7m 2 .s ´1, and specific heat from 879 to 1364J.kg ´1.K ´1.

Table 7 :

 7 Criterion for the choice of size for each zone

Table 8 :

 8 Influence of the convective heat transfer coefficient on the surface temperature (June 5 th and 6 th )

Table 9 :

 9 Influence of soil characteristics (June 5 th and 6 th )

Table 11

 11

	presents

Table 10 :

 10 Evaluation of the centimetric grid model according to the experimental data the total error due to the model and node distribution (first and third line) and the part of the error which is due to the reduction in the number of nodes (second and fourth line).

	Date	Model	1pt/cm 8/10/4 4/6/2 4/3/2
	June 5 th	RMSE with experimental	0.71	0.71	0.65	0.65
		data ( 0 C)
		RMSE with centimetric grid	-	0.03	0.13	0.36
		( 0 C)
	June 12 th Root mean square error	0.60	0.62	0.65	0.85
		with experimental data ( 0 C)
		RMSE with centimetric grid	-	0.03	0.15	0.42
		( 0 C)
	June 13 th RMSE with experimental	1.28	1.30	1.10	1.26
		data ( 0 C)
		RMSE with centimetric grid	-	0.03	0.44	0.60
		( 0 C)

Table 11 :

 11 Evaluation of the model with a reduced numbers of nodes

Table 12

 12 , the centimetric grid model is more accurate than the other models from the literature. This is a logical result as the centimetric model uses a high number of nodes.

	Using a grid with 15 nodes, Herb et al. (2008 [13]) obtained lower perfor-
	mances (1.58 0 C) than our model. The application conditions being similar,
	it can be affirmed that our model is at least as accurate as that of Herb
	et al. (2008 [13]) but with more optimized discretization. The difference
	in performance can mainly be explained by thinner discretization near the
	surface.

The model presented by

Malys (2012, [32]

) has only 4 nodes whereas the other models have more than 14 nodes. Consequently, it is the least accurate one, even though the model is applied with fair solicitations: a low amplitude of surface temperature variation (19 0 C) and low maximum temperature (41 0 C). Using a higher number of nodes with good distribution seems to be essential to obtain good performance. However, this result should be treated with caution since Malys (2012,

[START_REF] Malys | Évaluation des impacts directs et indirects des façades et des toitures végétales sur le comportement thermique des bâtiments[END_REF]

) does not use the same ground surface type as the other authors.

Table 13 :

 13 Distribution of the nodes and accuracy

Acknowledgements

This research work was carried out within the scope of the EVA Project, funded by the ADEME (French Environment and Energy Management Agency) under Contract No. 1216C0037. The authors are grateful to the ADEME and Véolia for their financial support of this study, as well as to IFSTTAR, LHEEA and ONEVU for providing us with the experimental data.

of the surface temperature requires that the heat flux transfer into the soil be correctly represented.

Discussion

This part is devoted to discussing the overall accuracy of the model.

The results are firstly compared to those presented in the literature. The accuracy of the model according to meteorological data is then evaluated with an annual simulation.

Comparison with other models accuracy

The performance of the models presented above may be compared to the performance of other models. From the literature the results of the following authors are chosen: Yang et al. (2013 [5]), Herb et al. (2008 [13]), Best (1998 [9]) , Malys, (2012 [32, 7]). Note that Malys, (2015, [START_REF] Malys | Évaluation des impacts directs et indirects des façades et des toitures végétales sur le comportement thermique des bâtiments[END_REF][START_REF] Musy | The use of solene-microclimat model to assess adaptation strategies at the district scale[END_REF]) applied the model proposed by Bouyer (2009, [4]). Those models are chosen because their simulation conditions (weather conditions, node distribution, type of ground surface, indicator calculated for performance evaluation) are the closest to those of the present study. Conditions of simulation and results are summarized in Table 12.

Three of the four articles ( [START_REF] Yang | Evaluation of a microclimate model for predicting the thermal behavior of different ground surfaces[END_REF][START_REF] Best | A model to predict surface temperatures[END_REF][START_REF] Herb | Ground surface temperature simulation for different land covers[END_REF]) have in common the type of surface studied and a gradual distribution of nodes: thinner at the surface than bellow. The last one, Malys (2012, [START_REF] Malys | Évaluation des impacts directs et indirects des façades et des toitures végétales sur le comportement thermique des bâtiments[END_REF][START_REF] Musy | The use of solene-microclimat model to assess adaptation strategies at the district scale[END_REF]), evaluated the previous SOLENEmicroclimat soil model with a grass surface and a mesh with only 4 nodes. difference is 0.58 0 C and for the 4/3/2 node model, 2.18 0 C. The mean absolute difference for the first node distribution is only 0.1 0 C whereas that of the 4/3/2 distribution is 0.33 0 C, which remains in the order of magnitude of the measurement uncertainty.

For each case, the frequency of appearance of the mean daily error between the centimetric grid and the reduced grid model is calculated. Each day is ranked in a different performance class to illustrate the differences between the node distribution (Table 14).

Mean daily error Low : E ă 0.


Jérémy Bernard

:,; ID Auline Rodler

Benjamin Morille

Xueyao Zheng

How to design a park and its surrounding urban morphology to optimize the spreading of cool air?

Keywords: park cool island, urban cooling, urban morphology, micro-climate simulations

Green areas induce smaller increases in the air temperature than built-up areas. They can offer a solution to mitigating the Urban Heat Island impacts during heat waves since the cool air generated by a park is diffused into its immediate surroundings through forced or natural convection.

The purpose of this study is to characterize the effect of several variables (park size, morphology of surrounding urban area and wind speed) on the spreading of cool air. A parametric study is performed to run computational fluid dynamic simulations. The air temperature entering the computational domain was set at 35°C, and the 2-m high surface included within the 34°C isotherm was defined as an indicator of cool air spreading. The effects of park shape and orientation are negligible in comparison with its size effects. The number of buildings is better correlated with the cooled surface area than the typical urban parameters identified in the literature (i.e. building density, aspect ratio or mean building height). Since the number of buildings is obviously related to the number of streets, this result suggests that the greater the number of streets around a park, the wider the area that cool air spreads.

Introduction

According to the Intergovernmental Panel on Climate Change (IPCC) reports, heat wave intensity and frequency should increase in the coming decades [START_REF] Revel | Climate change, impacts and vulnerability in Europe[END_REF], which may lead to situations of outdoor and indoor discomfort as well as major health impacts. During the summer of 2003, European countries recorded 70,000 excess deaths, attributing them to an unusually hot summer [START_REF] Robine | Death toll exceeded 70,000 in Europe during the summer of 2003[END_REF].

Urban areas are especially vulnerable when facing such problems, and this for two reasons. First, heat waves are exacerbated in cities when compared to their surroundings due to the Urban Heat Island (UHI) phenomenon [START_REF] Oke | Boundary layer climates[END_REF]. Second, urban areas concentrate most of the human population (i.e. 70% in European Union countries - [START_REF]Eurostat regional yearbook 2014: Focus on european cities[END_REF]).

In order to cool cities, several solutions are investigated: evaporative techniques (fountains, water ponds, street watering), green techniques (grass, trees, green facades or green roofs, etc.), and material techniques (reflective, water retentive, etc.). All of these techniques have shown cooling potential, but further research is still needed to improve their performance [START_REF] Santamouris | Passive and active cooling for the outdoor built environment-Analysis and assessment of the cooling potential of mitigation technologies using performance data from 220 large scale projects[END_REF]. This article focuses on the use of parks as cooling solutions.

Several studies have shown that parks may create cool air during the day [START_REF] Barradas | Air temperature and humidity and human comfort index of some city parks of Mexico City[END_REF] but also at night [START_REF] Jauregui | Influence of a large urban park on temperature and convective precipitation in a tropical city[END_REF]. Moreover, this cool air may be transported to the neighborhoods surrounding parks [START_REF] Chang | Effects of urban parks on the local urban thermal environment[END_REF][START_REF] Spronken-Smith | Scale modelling of nocturnal cooling in urban parks[END_REF]. The ability of a park to cool the air and spread this cooler air to its vicinity depends mainly on three factors: park characteristics, wind speed conditions, and urban morphology of the peripheral areas. The characteristics of a park (tree species and density, presence of water ponds, soil types, etc.) are key parameters in explaining its cooling potential [START_REF] Spronken-Smith | Scale modelling of nocturnal cooling in urban parks[END_REF]. The cool air produced within the park is transported into the surrounding areas by diffusion, advection and convection. For the sake of simplification, the term spreading will be used in the following sections. The spreading of cool air varies in intensity depending on both wind conditions and the urban layout of the park surroundings. This article focuses solely on the aspect of cool air spreading (i.e. on the effect of wind speed and urban layout on the level of spreading).

In 1991, Jauregui [START_REF] Jauregui | Influence of a large urban park on temperature and convective precipitation in a tropical city[END_REF] established that park cooling may be measured up to one park's width away from its boundaries, which is consistent with several other studies [START_REF] Spronken-Smith | Scale modelling of nocturnal cooling in urban parks[END_REF][START_REF] Cao | Quantifying the cool island intensity of urban parks using ASTER and IKONOS data[END_REF][START_REF] Ca | Reductions in air conditioning energy caused by a nearby park[END_REF][START_REF] Shashua-Bar | Vegetation as a climatic component in the design of an urban street: An empirical model for predicting the cooling effect of urban green areas with trees[END_REF]. However, the methodology used to assess this distance was not clearly stated, which makes the results difficult to replicate. To explain the relationship between air temperature drop and distance from the park, Shashua-Bar and Hoffman [START_REF] Shashua-Bar | Vegetation as a climatic component in the design of an urban street: An empirical model for predicting the cooling effect of urban green areas with trees[END_REF], Doick et al. [START_REF] Doick | The role of one large greenspace in mitigating London's nocturnal urban heat island[END_REF] proposed an asymptotic nonlinear model, which was then used to define the Park Cool island Distance (PCD) as the 'distance [from the park] where 10% of the UHI is still present' [START_REF] Doick | The role of one large greenspace in mitigating London's nocturnal urban heat island[END_REF]. Despite its interest, this method has been applied based on air temperature observations for a single park, hence the model is only valid for this specific park and cannot be used to draw general conclusions. We have previously seen that the production of cool air depends on park characteristics, whereas its spreading depends on both wind conditions and the layout of the urban surroundings.

The spreading of cool air is affected by wind speed value [START_REF] Chandler | The Climates of London[END_REF]. Under light winds, the temperature difference observed between the park and the streets is responsible for the advection phenomenon [START_REF] Upmanis | The influence of green areas on nocturnal temperatures in a high latitude city (Göteborg, Sweden)[END_REF][START_REF] Jansson | Near surface climate in an urban vegetated park and its surroundings[END_REF]. In the presence of high wind speeds, this phenomenon becomes negligible. Skoulika et al. [START_REF] Skoulika | On the thermal characteristics and the mitigation potential of a medium size urban park in Athens, Greece[END_REF] showed that for a wind speed above 5.5 m/s, the PCD decreases linearly as wind speed value increases, whereas Doick et al. [START_REF] Doick | The role of one large greenspace in mitigating London's nocturnal urban heat island[END_REF], Upmanis et al. [START_REF] Upmanis | The influence of green areas on nocturnal temperatures in a high latitude city (Göteborg, Sweden)[END_REF], Oke et al. [START_REF] Oke | The micrometeorology of the urban forest [and discussion[END_REF] considered the PCD to be negligible whenever wind speed exceeds 2.3 m/s, 5 m/s and 6 m/s, respectively. According to Doick et al. [START_REF] Doick | The role of one large greenspace in mitigating London's nocturnal urban heat island[END_REF],

this wind speed threshold is likely a function of street geometry, but they did not study this point in further detail.

Chandler [START_REF] Chandler | The Climates of London[END_REF] demonstrated that some urban structural parameters characterizing the surroundings of a park play a key role in the spreading of cool air. Chang and Li [START_REF] Chang | Effects of urban parks on the local urban thermal environment[END_REF] analyzed the air temperature gradient around 60 parks relative to the urban canyon dimensions. Their results agreed with the previous literature regarding the average cooling distance from the park, i.e. approx. one park width. They were unable however to establish any relationship between cool air spreading and the aspect ratio or building density. This shortcoming might be explained by the use of experimental data, which make it difficult to differentiate the influence of street dimensions on both radiation trapping and the reduction of cool air flow originating from the park. Moreover, the performance of the shelter protecting the air temperature sensor from the sun was very sensitive to solar radiation conditions [START_REF] Lacombe | WMO Field Intercomparison of Thermometer Screens/Shields and Humidity Measuring Instruments[END_REF][START_REF] Lacombe | Results of the WMO intercomparison of thermometer screens/shields and hygrometers in hot desert conditions[END_REF]. A temperature measurement conducted with a shelter located in a street exposed to direct solar radiation may be overestimated, thus complicating comparisons with a measurement performed in a shaded street.

The objective of this article is to better understand the phenomenon of cool air drainage from a park to its surroundings, particularly the effects of wind conditions and street dimensions on the spreading of air. Simulation rather than observation is used in order to better control the key variables, such as wind speed, urban layout, and both park and urban thermal characteristics. Simulation is also appropriate to facilitate obtaining a temperature field, thus avoiding having to deploy a large number of sensors, which could lead to numerous observation issues [START_REF] Oke | Initial guidance to obtain representative meteorological observations at urban sites[END_REF].

Methodology

Definition of urban form parameters

A numerical approach is used herein to establish relationships between cool air spreading and both wind speed and urban form parameters. The urban form parameters are selected from a literature review, namely: building height, building density, and aspect ratio [START_REF] Chang | Effects of urban parks on the local urban thermal environment[END_REF][START_REF] Doick | The role of one large greenspace in mitigating London's nocturnal urban heat island[END_REF][START_REF] Perini | Effects of vegetation, urban density, building height, and atmospheric conditions on local temperatures and thermal comfort[END_REF]. For purposes of this analysis, a square park has been designed, and the street width (W), building footprint area (S B ) and reference area (S re f ) are all shown in Figure 1. The aspect ratio (H{W) and building density (D B ) are calculated by means of Equations 1 and 2:

H{W " H W ( 1
)
where H is the building height

D B " S B S re f ´Spark (2)
Since this investigation concerns the effect of the urban form on air spreading, the scenario used as a reference is represented by a park on its own (the reference area contains a park but no building and consists of the entire square surface area).

Model and implementation

The park is assumed to produce cool air that will then be spread over the urban area. Depending on the morphology of the urban area, radiation trapping may vary in intensity and thus interfere with the cool air originating from the park. To focus solely on air spreading and simplify the simulation analysis, we assume herein that the urban surfaces do not exchange any energy with the air, whereas the park surface cools the air at a constant energy rate of 300W{m 2 . This initial assumption is a major one: during a summer, Rodriguez [START_REF] Rodriguez | FluxSAP -A collaborative experimental campaign on water and energy fluxes in urban areas and the relation with the vegetation : the case of a Nantes district[END_REF] measured in Nantes (city in western France) both the latent and sensible heat fluxes according to an eddy-covariance method. At noon, they observed an average sensible heat flux value of 250W{m 2 for a highly urban neighborhood (i.e. vegetation density less than 10% in a 200-m buffer circle around the station) versus just 100W{m 2 for a more highly vegetated area (76% vegetation). The second assumption is much closer to reality, at least for cities located in the Cfb 1

or Dfb 2 climate zone [START_REF] Peel | Updated world map of the Köppen-Geiger climate classification[END_REF]. The cooling flux of 300W{m 2 corresponds to the maximum latent heat flux measured around noon in the urban forested park of Chicago [START_REF] Mcpherson | Chicago's urban forest ecosystem: results of the Chicago Urban Forest Climate Project[END_REF] and is consistent with observations recorded in the City of Nantes in a fairly green urban neighborhood (measurement performed using the eddy-covariance method with a vegetation density in the 200-m buffer circle around the station of 76%) [START_REF] Rodriguez | FluxSAP -A collaborative experimental campaign on water and energy fluxes in urban areas and the relation with the vegetation : the case of a Nantes district[END_REF]. The results of our study may not be replicable for any park type (differences in the proportion of grass, tree, concrete, etc.) or climate zone (differences in the amount of solar radiation, air humidity, air temperature, etc.) since the cooling intensity may be affected [START_REF] Santamouris | Passive and active cooling for the outdoor built environment-Analysis and assessment of the cooling potential of mitigation technologies using performance data from 220 large scale projects[END_REF]. Several other latent heat flux values (100 and 200W{m 2 ) have been tested, but the natural convection effect was harder to observe. We thus decided to apply the method for the most sensible case (300W{m 2 ). The air temperature was also set very high (35°C) to reflect the peak temperature that Europe has had to face during its most recent heat waves.

To assess the spreading phenomenon of the cooled air generated by a park, the Computational Fluid Dynamics (CFD) tool Code_Saturne is used in the environment of SOLENE-microclimat (Morille et al. [START_REF] Morille | SOLENE-microclimate: a tool to evaluate envelopes efficiency on energy consumption at district scale[END_REF], Musy et al. [START_REF] Musy | The use of SOLENE-microclimat model to assess adaptation strategies at the district scale[END_REF]), a numerical tool dedicated to urban climate modeling. Thanks to Code_Saturne, the airflow is computed by resolving Navier-Stokes equations (momentum, mass continuity, energy, species transport) using a k ´ turbulence model. The tool configuration is similar to that described in Malys [START_REF] Malys | Évaluation des impacts directs et indirects des façades et des toitures végétales sur le comportement thermique des bâtiments[END_REF], except that the buoyancy forces have been modeled in order to obtain realistic airflow when natural convection is predominant (i.e. under low wind speed and high temperature differences). These forces are modeled using the Boussinesq approximation, which considers air density to be a function of air temperature.

ρ " ρ r ¨p1 ´β ¨pT ´Tr qq

where:

T r is the reference air temperature (" 300K)

T the air temperature ρ r the air density at temperature T r (" 1.18kg ¨m´3 ) ρ the air density

β the coefficient of thermal expansion (" 1 T r )
The computational domain represents a 100-ha square area discretized using 270,000 tetrahedral meshes (Figure 2). The mesh size varies within the domain. The minimum size of 1.5 m for the tetrahedron side is set in the region of interest (the park and its close surroundings under 2 m high). A sensitivity analysis has been performed in order to evaluate the effect of mesh size on air temperature and wind speed values calculated at a height of 2 m. Several meshing scenarios with different mesh sizes have been tested. For each case, air temperature and wind speed have been calculated according to the following procedure:

1. An interpolation is performed to calculate the air temperature of the vertex from the air temperature in the tetrahedron 2. A 2-m high horizontal plane on the ground is generated and intersects the tetrahedron vertex 3. The air temperature at these points is the result of a linear interpolation between the segment ends.

Results of this analysis revealed that between the chosen scenario (270,000 tetrahedrons) and the most accurate one (507,000 tetrahedrons), the temperature error never exceeded 0.5°C inside the park and 0.02°C outside the park.

For each time step, the calculation stopped under a convergence criterion that needed to be set. This threshold was chosen based on several simulations, where both park size and wind speed varied but not wind direction. For each simulation, 2,000 iterations were performed and the temperature convergence parameter always stabilized above 10 ´4. This threshold was thus chosen as the convergence criterion.

Indicators of air spreading efficiency

To evaluate the spreading efficiency according to wind speed and urban form, performance indicators prove to be necessary. Several studies employed the concept of PCD, which is defined as the distance where the cooling induced by the park is still noticeable. In this study, the cooling effect of a park is considered to end where the 2-m high temperature drops below 1°C from the initial temperature. The PCD is then defined as the mean distance between the park boundary and the 34°C isoline at a height of 2 m above ground. The distance between each point of the isoline and the park is calculated differently depending on the specific location (Figure 3). Concerning zones 1, 2, 3 and 4, it is defined as the distance between the point and the closest corner of the park. As for zones 5, 6, 7 and 8, it is defined as the minimum distance between the point and a side of the park. The PCD however is a 1-dimensional indicator: in the case of a long but narrow isotherm shape, the average distance of the 34°C isotherm can be high, whereas the surface benefiting from the cooled air may be small. The surface of the area included in the 34°C isoline is then also calculated in order to overcome the shortcomings of the PCD. Yet these results show that both indicators exhibit similar behavior for all of the following analyses. For the sake of simplification, only the results obtained with the cooled surface area will be presented and analyzed.

Parametric study

The cooling surface is impacted by several kinds of parameters. Three of them will be investigated herein: park size (with the shape remaining square), the urban form of its surrounding areas, and wind conditions.

The park width is set at 50, 100 and 200 m. One simulation is run for each park width without any adjacent buildings.

Several Urban Form (UF) scenarios are adopted to investigate the effect of urban parameters on the spreading of cool air. Each scenario is established to define a building height, building size and street width in order to obtain a building density and aspect ratio values (Table 1). Values are chosen so as to highlight the separate effect of H{W, D B and H on the cool air spreading, but also to optimize the number of simulations to be carried out. Seven Urban Forms are proposed (Table 2). First, the building density influence is investigated by changing either the building size (UF 1 to 3 -method 1) or the street width (UF 2 to 4 -method 2). Three density levels are obtained and defined as low, medium and high. Next, the aspect ratio influence is investigated by adopting either a low density (UF 1 and 4)

or high density (UF 3 and 5). Lastly, the height is varied while the building density remains constant at a low value (UF 1 and 6) or high value (UF 3 and 7). For each scenario, the resulting number of buildings is also given in the Table. The influence of wind speed is investigated for each of these scenarios: values are selected every 0.5 m/s from 0.5 m/s to 5 m/s. The scenario combinations derived from urban form and wind speed lead to a total of 70 simulations. For all simulations, wind crosses the park along its diagonal (Figure 3). The effect of wind direction and park shape have been investigated, but only the conclusions will be briefly presented in the following section.

Results

Influence of park size

For each park size (50, 125 and 200 m wide), the variation in the cooled surface area (i.e. where temperature is lower than 34°C) is evaluated for various wind speeds. It is obvious that a larger park will cool a larger area. The cooled surface is then divided by the park area in order to obtain a cooling efficiency (Ce) indicator. This value may be interpreted as an equivalent surface cooled by each square meter of park.

The simulations show that Ce increases with park size or as wind speed decreases (Figure 4). The first observation implies that one square meter of a large park is more efficient than one square meter of a small park. The second observation may be attributed to the assumption made for the simulation: the cooling rate of the air located in the park is constant regardless of its air temperature or wind speed. One result of this assumption is that the longer the air stays in the park, the cooler it leaves the park. Thus, when wind speed is very low (0.5 m/s), the air temperature difference between the park and its surroundings is high (8°C). In this case, the air temperature of 34°C is easily reached and natural convection dominates: the cool air flow is emitted in all directions, including the opposite direction of the incoming wind, and the 34°C isotherm is quite large (Figure 5(a)). For higher wind speeds (2 m/s), the air temperature difference between the park and its surroundings is lower (2°C), in which case the air temperature barely reaches 34°C and spreading is mainly driven by the wind. As a result, the flow is channeled in a single direction and the 34°C isotherm is small (Figure 5(b)). Overall, upon analysis of the temperature distribution around a park, we can conclude that the lower the wind speed, the further the spreading of cooled air in the direction opposite the incoming air flow. The shape and orientation of the park (with respect to wind speed) have been investigated. For a given surface area, they both exert a very limited influence compared to park size (which is why the results are not shown in this study).

Overall, two main observations can be drawn:

• The larger the park, the greater the cooling intensity

•

The higher the wind speed, the smaller the cooled surface area.

In conclusion, the longer the air stays in the park, the cooler it is and the larger the size of the 34°isotherm.

Influence of the neighboring district

Adding buildings around the park has a direct impact on the cooled surface shape. To investigate the influence of the neighboring district on cool air spreading, the larger park is used (200 m wide).

The surface generated by the 34°C isotherm has nearly the same size for both the reference urban form (without any buildings) and the UF 1 (with buildings). However, the former is located further downstream than the latter (Figure 6). The buildings can be seen as contributing to slowing the wind speed. In agreement with the previous analysis, this finding leads to increasing the natural convection influence, diffusing the air further toward the incoming airflow direction. The influence of the urban form of the district surrounding the park is now investigated. Three parameters are considered: building density, aspect ratio, and building height. As described in the Methodology section, the building density is increased using two distinct methods (either the building size is modified -method 1 -or the street width is modified -method 2). These two methods are considered separately and results are shown on separate plots. As previously observed, the cooled surface decreases as wind speed increases under all scenarios (Figure 7). For a given wind speed, when building density is increased by increasing the building size, the cooled surface area decreases (Figure 7(a)). In contrast, when building density is increased by increasing street width, the cooled surface area expands (Figure 7(b)). Building density cannot therefore be used as the lone parameter explaining the influence of the built environment on the cooled surface of a park.

The effect of the aspect ratio on the cooled surface is tested by decreasing the street width. In this manner, the aspect ratio is increased, thus keeping building density and building height constant. For a given wind speed, the size of the cooled surface increases as the aspect ratio increases (as well for a given building density regardless of its value 8. Thirdly, the building height is increased for a given building density with an increasing aspect ratio. For high building density (Figure 9(a)) and at a given wind speed, the cooled surface area is nearly similar when the building height increases from 10 m to 17 m. For the low building density case (Figure 9(b)), it increases as the building height rises from 4 m to 10 m. According to an initial analysis, we could assume that building height, like building density, is not a key parameter in explaining the cooled surface, while the aspect ratio is a key parameter. It is also possible that the effect of building height only appears at a certain threshold value. This assumption will be further analyzed in the following section. The three urban form parameters (building density, aspect ratio and height) investigated herein do not impact the air spreading pattern in a way that could have been expected:

• For a given building height and aspect ratio, the building density increase enhances or reduces the spreading of cooled air, depending on how the density is modified. Density does not seem to be a key parameter affecting air spreading.

• For a given building height and building density, the aspect ratio appears to enhance the cooling process, whereas it is not often identified in the literature as a key parameter.

• Higher buildings may increase air spreading under a certain threshold value, which has not yet been identified (note that this threshold may be dependent on wind speed conditions).

Discussion

The results presented above lead to investigating in greater detail the simulation results and then discussing them. This discussion section will focus on two points:

• An understanding of the influence of building height on air spreading based on the configuration

•

The identification of other parameters that could be relevant for the air spreading characterization.

The air temperature field for UF 3 and UF 7 (buildings 10 m and 4 m high, respectively, in a high-density district due to large buildings -few streets) is presented in Figure 10. In the low building height case (UF 7), the cool air diffuses over the buildings (Figure 10(a)). In the tall building case (UF 3), the buildings seem to be sufficiently high to avoid air from spreading above them (Figure 10(b)). As a result, the cool air is channeled into the streets, reaching a further distance from the park and thus covering a wider area. Above a certain building height threshold, the cool air can no longer spread over the buildings.

Further increasing the building height would then have no effect on the cooled surface since the air is only spreading under the urban canopy.

The issue of building height relevance calls into question the relevance of two other indicators, namely: urban fragmentation (e.g. number of buildings), and street width. We previously showed that the aspect ratio is positively correlated with the cooled surface area. The aspect ratio is defined as the ratio of building height to street width. Since past results have demonstrated that building height only affects the cool air spreading in very specific configurations, we might wonder whether street width is an influential parameter. Urban fragmentation was not initially identified as an important parameter as regards the cool air spreading. However, the number of buildings varies along the chosen urban form and thus its cooling potential warrants investigation. The cooled surface area is plotted versus both the number of buildings (Figure 12(a)) and the street width (Figure 12(b)). The number of buildings appears to be much more closely correlated than street width with the cooled surface area. A large variability in the cooled surface area exists for the UF with 8-m width streets, thus making this parameter insignificant (Figure 12(b)). In Figure 12(a), the cooled surface area clearly increases as the number of buildings or streets increases. Moreover, the distance to the regression line remains low (except when the number of buildings is very small).

Ultimately, urban fragmentation (here the number of buildings) seems to be a more relevant variable than building density, aspect ratio or building height for characterizing the cooled air spreading from a park through a district.

Conclusion

The parameters impacting the spreading of cool air produced by a park have been investigated

herein. Simulations were performed to observe the cooling distance (surface included inside the 1°C temperature decrease isoline) induced by the park and by the form of its surrounding urban district.

Some of the studied parameters are intrinsic to the park (shape and size), while others describe the urban morphology of the surrounding neighborhoods (building density, building height and aspect ratio). The simulations were run for various wind speeds and under a constant cooling flux of 300 W/m² generated by the green area.

The analysis of results leads to the following conclusions:

• A larger park will cool more efficiently: it will generate a greater surface area of cooled air per square meter than a smaller park.

•

The building density variation can enhance or reduce the spreading of cooled air. This parameter therefore is irrelevant as regards cool air spreading.

• For a given building height and building density, the aspect ratio seems to enhance the cooling process, even though it has not often been identified in the literature as a key parameter.

• Building height does not exert any influence on air spreading above a certain threshold value.

•

The number of buildings seems to be a very relevant parameter for characterizing the cooled air spreading from a park through a district. This parameter and all other fragmentation indicators (the number of buildings is an indicator of fragmentation intensity, but the direction and shape of fragmentation may also be of interest) should be further analyzed by the research community.

In this paper, a constant cooling flux was considered in order to focus solely on the spreading of the cool air being produced. Some bias was thus introduced when using this method for analysis (if air remains for an infinite time in the park, its temperature should be infinitely negative). In reality, the flux should have been modified by the temperature difference between the ground and the air and by the wind speed. This hypothesis could be revised in future work by employing a coupled heat transfer simulation (thermo radiative and CFD) in order to generalize the methodology proposed herein and consider a cooling flux that varies in both space and time.

This paper has highlighted new findings but it is based on strong hypotheses. Further investigation is thus required (empirical data should be used) to verify that these assumptions do not affect the veracity of the main findings. However, only a few studies regarding the topic of cool park air spreading have so far been performed using a numerical approach. The academic approach proposed in this paper has led to both a better understanding of the physical phenomena involved and an identification of the key parameters affecting cool air spreading. For these reasons, it is complementary to the experimental studies available in the literature.

Figure 1 .

 1 Figure 1. Study area and parameter definition -note that the reference surface consists of the entire square surface area

1 Temperate

 1 Without dry season Warm summer 2 Cold Without dry season Warm summer The mesh size then increases linearly up to the limits of the domain (10 m at surface level, and 50 m at the top of the domain).

Figure 2 .

 2 Figure 2. Meshed computational domain.

Figure 3 .

 3 Figure 3. Zoning (dashed lines and numbers) and example of calculation (red lines) used to calculate the distance between any point on the 34°C isoline and the park boundary.

Figure 4 .

 4 Figure 4. Cooled surface evolution with wind speed

Figure 5 .

 5 Figure 5. 2 m high air temperature distribution inside and around a park with a width of 125 m: (a) Exposed to a wind speed of 2 m/s. (b) Exposed to a wind speed of 0.5 m/s -the white line is the 34°C isoline

Figure 6 .

 6 Figure 6. 34°C isotherm at a 0.5 m/s wind speed for the reference urban form (no buildings -orange line) and for the UF 1 (red line) -the small white squares represent building footprints

Figure 7 .

 7 Figure 7. Influence of wind speed on the cooled surface for various building density values: (a) When building density is increased modifying the building size (b) When building density is increased modifying the street width

Figure 8 .

 8 Figure 8. Influence of wind speed on the cooled surface for different aspect ratio values. (a) For high building density (b) For low building density

Figure 9 .

 9 Figure 9. Influence of wind speed on the cooled surface for various building heights. (a) For a low building density district (b) For a high building density district

Figure 10 .

 10 Figure 10. Vertical slice of the air temperature field for a high building density district. (a) UF 7 : buildings are 4 m high. (b) UF 3 : buildings are 10 m high -the white line is the 34°C isotherm

Figure 11 .

 11 Figure 11. Vertical slice of the air temperature field for a low building density district. (a) UF 1 : buildings are 10 m high. (b) UF 6 : buildings are 17 m high) -the white line is the 34°C isotherm

Figure 12 .

 12 Figure 12. Cooled surface area evolution versus: (a) the number of buildings; (b) street width -wind speed set at 2.5m{s

Table 1 .

 1 Parameters set for each of the seven scenarios investigated

	UF height (m) building width (m) street width (m)	D B (%)	H/W number of buildings
	1	10	8	8	25 (low)	1.25	672
	2	10	18	8	48 (medium)	1.25	260
	3	10	96	8	85 (high)	1.25	12
	4	10	18	18	25 (low)	0.56	160
	5	10	18	5	95 (high)	2.00	360
	6	17	8	8	25 (low)	2.12	672
	7	4	96	8	85 (high)	0.50	12

Table 2 .

 2 Geometry of each urban form

	UF1 / UF6	UF2	UF3 / UF7
	UF4	UF5

Author Contributions: This work was been conducted within the scope of Xueyao Zhang's internship. The work program was managed by Jérémy Bernard, Benjamin Morille and Auline Rodler. Xueyao Zhang studied the bibliography and carried out all of the simulation runs. Jérémy Bernard, Benjamin Morille and Auline Rodler wrote the paper and contributed to this effort in equal measure.

Conflicts of Interest:

The authors declare no conflict of interest."