Mohamad El Houssami

Aymeric Lamorlette

Dominique Morvan

Rory M Hadden

Albert Simeoni
email: asimeoni@wpi.edu

Framework for submodel improvement in wildfire modeling

Keywords:

An experimental and numerical study was carried out to assess the performance of the different submodels and parameters used to describe the burning dynamics of wildfires. A multiphase formulation was used and compared to static fires of dried pitch pine needles of different bulk densities. The samples were exposed to an external heat flux of 50 kW/m 2 in the FM Global Fire Propagation Apparatus and subjected to different airflows, providing a controlled environment and repeatable conditions. Submodels for convective heat transfer, drag forces, and char combustion were investigated to provide mass loss rate, flaming duration, and gas emissions. Good agreement of predicted mass loss rates and heat release rates was achieved, where all these submodels were selected to suit the tested conditions. Simulated flaming times for different flow conditions and different fuel bulk densities compared favorably against experimental measurements. The calculation of the drag forces and the heat transfer coefficient was demonstrated to influence greatly the heating/cooling rate, the degradation rate, and the flaming time. The simulated CO and CO 2 values compared well with experimental data, especially for reproducing the transition between flaming and smoldering. This study complements a previous study made with no flow to propose a systematic approach that can be used to assess the performance of the submodels and to better understand how specific physical phenomena contribute to the wildfire dynamics. Furthermore, this study underlined the importance of selecting relevant submodels and the necessity of introducing relevant subgrid-scale modelling for larger scale simulations.

Introduction

Understanding the burning behavior of wildland fuels is essential to develop a complete understanding of wildfire spread. For the last decades, numerous studies have been attempted to better describe the fire dynamics of solid fuels [START_REF] Fernandez-Pello | The solid phase[END_REF][START_REF] Babrauskas | Ignition Handbook[END_REF][START_REF] Drysdale | An Introduction to Fire Dynamics[END_REF] , including polymers [START_REF] Tewarson | Flammability of plastics-I. Burning intensity[END_REF][START_REF] Ohlemiller | Modeling of smoldering combustion propagation[END_REF][START_REF] Sohn | Transient modeling of thermal degradation in non-charring solids[END_REF][START_REF] Tewarson | Combustion efficiency and its radiative component[END_REF][START_REF] Stoliarov | Prediction of the burning rates of non-charring polymers[END_REF] , cellulosic materials [START_REF] Nakabe | Ignition and transition to flame spread over a thermally thin cellulosic sheet in a microgravity environment[END_REF][START_REF] Olson | Experimental observations of spot radiative ignition and subsequent three-dimensional flame spread over thin cellulose fuels[END_REF] , and different species of wood [START_REF] Parker | Prediction of the heat release rate of wood[END_REF][START_REF] Boonmee | Glowing ignition of wood: the onset of surface combustion[END_REF][START_REF] Blasi | Modeling chemical and physical processes of wood and biomass pyrolysis[END_REF][START_REF] Fateh | Multi-scale modeling of the thermal decomposition of fire retardant plywood[END_REF] . The transfer of this knowledge to describe the burning behavior of wildland fuels, which are natural polymers, is not trivial because they differ from classical solid fuels, as they often consist in highly porous media where the solid phase is very small compared to the gas phase. The main challenge arises from the porous nature of the fuels and its coupling with the transient surrounding environment, which strongly influences the heat transfer, degradation rates, and the burning dynamics [START_REF] Grishin | A mathematical modelling of forest fires and new methods of fighting them[END_REF] . Indeed, during combustion of a highly porous bed, the heat captured by the solid phase is a result of the energy balance from the flame feedback by radiation and from the fresh air entrained into the fuel layer by convection, but also reacting with the pyrolysis products and with the char residue. As a consequence, these homogeneous and heterogeneous combustions contribute to transfer energy to the solid phase. This behavior is different from that observed for solid fuels for which the radiative and convective contributions of the flame can be considered to be transferred into the solid by conduction [START_REF] Torero | Heat and mass transfer in fires: scaling laws, ignition of solid fuels and application to forest fires[END_REF] . This highly porous configuration is typically found in pine needle beds, which are often found on the forest floors. They can typically be found accumulated near structures in the Wildland Urban Interface (WUI), increasing the fire risk [START_REF] Mell | The wildland-urban interface fire problem -current approaches and research needs[END_REF] . As the need to understand the burning dynamics of wildland fuels is a matter of great urgency in order to improve wildfire assessment tools, many studies have been performed over the last decades. Among those, Dimitrakopoulos and Papaioannou [START_REF] Dimitrakopoulos | Flammability assessment of mediterranean forest fuels[END_REF] used the cone colorimeter to determine the ignition time and the flammability of typical Mediterranean forest fuels and proposed four flammability classes. Schemel et al . [START_REF] Schemel | A calorimetric study of wildland fuels[END_REF] reported a calorimetric study of pine needle beds using the FM Global Fire Propagation Apparatus (FPA) Volumetric rate of charring, pyrolysis, and vaporization which it was shown that calculating the heat release rate (HRR) by means of calorimetry for pine needles can be reinforced by the use of mass loss rate and by knowing the heat of combustion in well ventilated test conditions. Bartoli et al . [START_REF] Bartoli | Determination of the main parameters influencing forest fuel combustion dynamics[END_REF] demonstrated that permeability drives the burning dynamics in porous beds and that the energy released increases with permeability. Additionally it was shown that for a given permeability, the fuel species have an influence on times to ignition and duration of flames. Mindykowski et al . [START_REF] Mindykowski | Piloted ignition of wildland fuels[END_REF] conducted experiments in the same experimental setup to investigate the piloted ignition time of litters composed of oven Mediterranean pine species and kermes oak leaves.

The porous bed can be described through CFD modeling using a multiphase formulation that allows the fire rate of spread (ROS) [START_REF] Porterie | Firespread through fuel beds: modeling of wind-aided fires and induced hydrodynamics[END_REF][START_REF] Porterie | Predicting wildland fire behavior and emissions using a fine-scale physical model[END_REF][START_REF] Morvan | Numerical study of the behaviour of a surface fire propagating through a firebreak built in a mediterranean shrub layer[END_REF] , and temperature field to be estimated [START_REF] Morvan | Modeling the propagation of a wildfire through a mediterranean shrub using a multiphase formulation[END_REF] . For instance, Padhi et al. [START_REF] Padhi | Computational investigation of flame characteristics of a non-propagating shrub fire[END_REF] described the flow field generated due to the steady burning of a shrub using the multiphase formulation. Velocity and temperature profiles were analyzed in the continuous flame region and thermal plume region and it was highlighted that the porous nature of the shrub affects air entrainment into the shrub, thus affecting the nature of the fire plume. Hence, it is important to evaluate the model's capabilities to simulate air entrainment and air mixture in the fuel bed. The use of the multiphase formulation is not new, it has been applied in numerous studies to describe the interaction between the atmospheric boundary layer with a canopy problems using isothermal flows [START_REF] Shaw | Large-eddy simulation of turbulent flow above and within a forest[END_REF][START_REF] Patton | Large-eddy simulation of windbreak flow[END_REF][START_REF] Shaw | Canopy element influences on resolved-and sub-grid-scale energy within a large-eddy simulation[END_REF][START_REF] Watanabe | Large-eddy simulation of coherent turbulence structures associated with scalar ramps over plant canopies[END_REF][START_REF] Katul | One and two-equation models for canopy turbulence[END_REF][START_REF] Foudhil | A fine-scale K-Epsilon model for atmospheric flow over heterogeneous landscaped[END_REF][START_REF] Gavrilov | Numerical simulation of coherent turbulent structures and passive scalar dispersion in canopy sublayer[END_REF] . However, the multiphase formulation is usually applied to simulate large-scale wildfires under complex conditions without thoroughly pre-validating or verifying if the selected submodels (i.e. the closure models) are well adapted to the fuel conditions and to the surrounding environment in which they are used, and if they allow the relevant physics to be captured [START_REF] Grishin | A mathematical modelling of forest fires and new methods of fighting them[END_REF][START_REF] Larini | A multiphase formulation for fire propagation in heterogeneous combustible media[END_REF][START_REF] Morvan | FIRESTAR: a physically based model to study wildfire behaviour[END_REF][START_REF] Mell | Numerical simulation and experiments of burning douglas fir trees[END_REF][START_REF] Consalvi | Numerical study of piloted ignition of forest fuel layer[END_REF] . For instance, degradation rates are provided from TGA analysis at a much smaller scale (microscopic) than reality [START_REF] Grishin | A mathematical modelling of forest fires and new methods of fighting them[END_REF] , and drag coefficients are derived from wind tunnel tests for simple geometries at ambient temperatures [START_REF] Incropera | Fundamentals of heat and mass transfer[END_REF] . In a recent study, Hoffman et al . [START_REF] Hoffman | Evaluating crown fire rate of spread predictions from physics-based models[END_REF] evaluated crown fire spread rate using FIRETEC [START_REF] Linn | Studying wildfire behavior using FIRETEC[END_REF] and WFDS [START_REF] Mell | WFDS preliminary user's guide[END_REF] (both include a multiphase approach) and compared the results to a compilation of wildfire observations in North American forests [START_REF] Alexander | Evaluating a model for predicting active crown fire rate of spread using wildfire observations[END_REF] . It was noted that the ability to test the model was limited by a lack of environmental and fuel data such as errors associated with point-to-point comparison, which is a limitation of large scale fires. This aspect is essential because the rate of spread predictions from detailed physics-based models are sensitive to small variations in both the spatial pattern of the fuels and the environment. Therefore, this study recommended further assessment of detailed physics-based models, particularly by providing additional data regarding fuel and environmental characteristics (i.e. wind). This is why it is necessary to establish a framework tailored to the development of fire modeling with the multiphase approach. It allows to study elementary aspects of the problem and to gradually move towards complexity (by changing the fuel or environmental properties). This methodology follows a building block approach to model development and facilitates a better understanding of forest fuel flammability and of its corresponding fire dynamics.

This study includes well documented fire experiments that are conducted in a controlled environment, providing precise measurements for different fuel and ambient conditions, to quantify the influence of the parameters on the models numerical predictions. As these conditions change, some parameters and submodels in the numerical model are no longer valid and need to be adjusted. This parametric study is often overlooked at larger scale due to the excessive uncertainties caused by the large variations of fuel properties (i.e. special distribution, fuel moisture content), the unstable wind (i.e. gusts), and the strong coupling between submodels, which makes it impossible to pinpoint exactly which submodel is not behaving physically. In this study, experiments were undertaken in the FM Global Fire Propagation Apparatus (FPA) [START_REF]Reaction to fire tests -measurement of material properties using a fire propagation apparatus. International Organization for Standardization[END_REF] , similarly as in [START_REF] Schemel | A calorimetric study of wildland fuels[END_REF][START_REF] Bartoli | Determination of the main parameters influencing forest fuel combustion dynamics[END_REF] where it was demonstrated to proving repeatable conditions. Pine needle beds of varying porosity were subjected to a radiative heat and various values of airflows to observe the burning behavior under different conditions relevant to wildfire spread. Pine needle beds were used as a reference fuel because they allow repeatable fuel bed properties under laboratory settings. Mass loss rates, heat release rates, flaming times, gas emissions, and temperature fields were used to compare the experiments to the simulations. Particular attention was given to forced flow conditionsextending previous work on natural convection [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] -because in wildfires, the flow is usually highly unstable and not well defined. Consequently, local variations in the flow can have important impacts on the local fire regime. However, the flow is controlled and can be well characterized in the FPA. This framework, with the appropriate modifications supports the development of large scale CFD modeling by providing inputs and indications for the necessary subscale modeling.

Methods

Experimental details

Experiments were conducted by burning beds of pine needles in the FPA [START_REF]Reaction to fire tests -measurement of material properties using a fire propagation apparatus. International Organization for Standardization[END_REF] ( Fig. 1 ). The apparatus provides a controlled environment with respect to air flow and repeatable, uniform heating conditions. Samples were exposed to a radiative heat flux of 50 kW/m 2 applied to the top surface of the sample, throughout the entire experiment mimicking a strong flame feedback from a larger fire surrounding the sample. This value of received heat flux is representative of typical fire front propagation [START_REF] Silvani | Fire spread experiments in the field: temperature and heat fluxes measurements[END_REF] . Dead pitch pine ( Pinus rigida ) needles were packed in cylindrical porous baskets (67%% opening fraction) of 12.6 cm diameter and 3.0 cm depth following the protocol used in [START_REF] Thomas | An experimental study evaluating the burning dynamics of pitch pine needle beds using the FPA[END_REF] . Fuel moisture content (FMC) was determined by conditioning needles at 60 °C for 24 h. The average FMC was 7%% of the dry weight. The surface to volume ratio ( σ s ) is 7295 m -1 , and the density ( ρ s ) is 607 kg/m 3 .

Airflow was introduced at the inlet under the sample and passed through and around the porous sample. Three inlet fluxes of air were used No Flow (NF), Low Flow (LF) (50 L/min, corresponding to a free stream velocity of 6.67 cm/s), and High Flow (HF) (200 L/min, corresponding to 26.8 cm/s). A hot wire anemometer (Kimo ® AMI301 with a 0.01 m/s resolution) was used to estimate the averaged velocity of the flow penetrating, and the flow circling around the porous sample. Even if the flow conditions used in this study are relatively low compared to the mean velocities that can be found in wildfires [START_REF] Nelson | Flame characteristics of wind-driven surface fires[END_REF][START_REF] Pitts | Wind effects on fires[END_REF] , it is important to assess the model's performance under these specific conditions, which could also be found locally during a fire. Following this framework, other flow regimes and more moderate flows could be studied.

Mass loss rate (MLR) was derived from measurements of mass using a load cell (0.001 g resolution) and the exhaust gases were analyzed for composition (CO, CO 2 and O 2 ) using non-dispersive infrared and paramagnetic techniques at a sampling rate of 2 Hz. Flaming time was measured with a 0.5 s error margin and the end of flameout time was reached when the CO concentration dropped under 4 ppm.

Three bulk densities were tested: 23, 30, and 40 kg/m 3 , as shown in Fig. 2 . The corresponding masses and porosities are listed in Table 1 . 40 kg/m 3 is the maximum density that can fit in the sample basket without over compression and breaking of the needles. 23 kg/m 3 corresponds to the minimal load that needs to be used to avoid significant reflection of the basket when heated [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] . These values are in the same range as typical fuel loads measured in the forest [START_REF] Clark | Fire management and carbon sequestration in pine barren forests[END_REF] . More details of the experimental setup are provided in [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] . All experimental results were averaged over 3 repetitions or more, and the error bars represent the standard deviation for all experiments.

Computational modeling

In this study, simulations were carried out using ForestFire-FOAM (FFF) [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] . It has been previously demonstrated that FFF allows modeling the fire-induced behavior of a porous, reactive and radiative medium [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] . Conservation equations (mass, momentum, and energy) are solved for an averaged control volume at a scale sufficient to contain both coexisting gas and solid phases, and considering strong coupling between the phases [START_REF] Larini | A multiphase formulation for fire propagation in heterogeneous combustible media[END_REF] . Processes such as drying, pyrolysis, and char combustion are described through temperature-dependent interaction between the solid and gas phases. The gas phase combustion is represented using an extension of the Eddy Dissipation Concept (EDC), in which the characteristic time scale of fuel-air mixing is different under turbulent and laminar flow conditions [START_REF] Ren | Large eddy simulation of turbulent vertical wall fires supplied with gaseous fuel through porous burners[END_REF] . This combustion model was used in a study that focused on fully-developed turbulent region of a wall flame where the diffusion time used for the laminar flow was a placeholder for the near wall region and was not verified. It is recognized that the validity of the modified EDC model for laminar flame studies still needs to be established, but this one was used in the absence of a better and simple model for laminar flames in LES. A complete description of FFF is available in [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF][START_REF] Morvan | Physical modelling of fire spread in grasslands[END_REF] and the governing equations are presented in the Appendix. The computational domain representing the FPA in a two dimensional configuration is sketched in Fig. 3 . The overall domain is a rectangle of 5 m wide and 2 m high. The mesh was composed primarily of hexahedral cells. Based on a mesh sensitivity study, a mesh resolution of 0.001 × 0.001 m 2 was determined as providing converging results in the vicinity of the sample and near the lamps. The mesh is stretched beyond the zone of interest until the boundaries reducing computational time without affecting the results. Adaptive time step is calculated based on Courant-Fredrichs-Lewy (CFL) number [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF] , with a maximum value of 0.7 to achieve temporal accuracy and numerical stability.

Sub-models for heat transfer, degradation, combustion, and gas emission have been implemented to allow better representation of these processes. However, in order to appropriately simulate the conditions presented above, specific submodels had to be adjusted. These submodels are presented hereafter and all other submodels that are not presented are discussed in the previous study [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] .

Drag forces

In the multiphase model, the drag force acting on the solidphase is often [START_REF] Porterie | Firespread through fuel beds: modeling of wind-aided fires and induced hydrodynamics[END_REF][START_REF] Porterie | Predicting wildland fire behavior and emissions using a fine-scale physical model[END_REF][START_REF] Mell | Numerical simulation and experiments of burning douglas fir trees[END_REF][START_REF] Zhou | Modeling of marginal burning state of fire spread in live chaparral shrub fuel bed[END_REF] approximated using the correlation proposed by Clift et al . [START_REF] Clift | Bubbles, drops and particles[END_REF] for estimating the drag force coefficient ( C D ) of spheres.

C D = 24 1 + 0 . 15 R e 0 . 687 Re × 3 8 ; 1 < Re < 10 0 0 (1)
To more closely represent the cylindrical geometry of the pine needles, the drag coefficient was approximated using a pseudo fluid model for arrays of emerging circular cylinders [START_REF] Cheng | Calculation of drag coefficient for arrays of emergent circular cylinders with pseudofluid model[END_REF] :

C D = 11 R e -0 . 75 + 0 . 9 1 -exp - 10 0 0 Re +1 . 2 1 -exp - Re 4500 ; 0 . 02 < Re < 2 × 10 5 (2)
The use of a correlation for arrays of cylinders is justified by the "sheltering" effect that diminishes the drag downstream an element [START_REF] Raupach | Drag and drag partition on rough surfaces[END_REF][START_REF] Mueller | LES modeling of flow through vegetation with applications to wildland fires[END_REF] . As element spacing decreases, the bulk drag coefficient decreases [START_REF] Nepf | A model for diffusion within an emergent plant canopy[END_REF] . Hence, correlations established for isolated elements are no longer effective because of this strong sheltering effect.

Convective heat transfer

The term representing the contribution of convective heat transfer Q (S) CONV between gas and solid phases is described as follows in the gas phase energy balance equation (in Appendix):

Q ( S ) CONV = α s σ s χ ( T -T s ) (3)
With T and T s the gas and solid phase temperatures, respectively. The heat transfer coefficient χ is estimated from the Nusselt number [START_REF] Incropera | Fundamentals of heat and mass transfer[END_REF] :

Nu = χ D k = CR e m P r n (4)
where, k is the air thermal conductivity (W/m.K) and D the equivalent diameter of a fuel particle (cylindrical needle) approximated as 4/ σ s . D is also used as the characteristic length for calculating the Reynolds number. As demonstrated in [START_REF] Lamorlette | Analytical quantification of convective heat transfer inside vegetal structures[END_REF] , despite the complex vegetative structure, a simple correlation (i.e. Eq. 4 ) is appropriate to represent the convective term, as long as the coefficients C, m and n are adapted to represent the observed packing effect, or "sheltering" effect [START_REF] Raupach | Drag and drag partition on rough surfaces[END_REF] . Since samples are prepared by stacking All the cited correlations are widely applied in studies that are using the multiphase approach [START_REF] Mell | Numerical simulation and experiments of burning douglas fir trees[END_REF][START_REF] Consalvi | Numerical study of piloted ignition of forest fuel layer[END_REF][START_REF] Morvan | Physical modelling of fire spread in grasslands[END_REF][START_REF] Zhou | Modeling of marginal burning state of fire spread in live chaparral shrub fuel bed[END_REF] , but there are no indications on which ones are more suitable.

Char oxidation

Previous work [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] has shown that using a 1st order Arrhenius submodel underestimates char oxidation rate ( ˙ ω char ) due to the low induced mass flux of oxygen. Consequently, the char oxidation reaction could not be sustained after flameout, despite the constant lamp radiation. In the present study, the same model was tested again, but this time in forced flow conditions to verify if it is sufficient to initiate char combustion without using more detailed modeling or additional forcing term depending on the Reynolds number to account for blowing effects [START_REF] Porterie | Predicting wildland fire behavior and emissions using a fine-scale physical model[END_REF] .The char oxidation rate is estimated from the following Arrhenius equation, where E CHAR is the energy activation of char oxidation and K CHAR is the preexponential coefficient, all derived from TGA analysis [START_REF] Grishin | A mathematical modelling of forest fires and new methods of fighting them[END_REF] .

˙

ω char = α s σ s υ O 2 α g ρY O 2 K CHAR e -E CHAR R T s ( 5
)
Char oxidation represents an important source of heat released during smoldering. Hence, it is considered that any material that forms char during its thermal decomposition can potentially sustain a smoldering process [START_REF] Ohlemiller | Smoldering combustion[END_REF] . This heterogeneous oxidation is incomplete and emits a higher yield of CO than the gas-phase combustion [START_REF] Rein | Carbon emissions from smouldering peat in shallow and strong fronts[END_REF] . By implementing the model of Evans et al. [START_REF] Evans | Combustion of wood charcoal[END_REF] into the multiphase formulation, the contribution of char combustion to the energy equation of the solid phase (in Appendix) becomes: [START_REF] Sohn | Transient modeling of thermal degradation in non-charring solids[END_REF] with α sg == 0.5, assuming that 50%% of the heterogeneous combustion takes place at the surface of the solid phase and the other 50%% takes place in the gas phase [START_REF] Porterie | Firespread through fuel beds: modeling of wind-aided fires and induced hydrodynamics[END_REF][START_REF] Porterie | Predicting wildland fire behavior and emissions using a fine-scale physical model[END_REF][START_REF] Larini | A multiphase formulation for fire propagation in heterogeneous combustible media[END_REF][START_REF] Morvan | Physical modelling of fire spread in grasslands[END_REF] . The split function ( ϕ) allows prediction of the ratio of CO and CO 2 concentrations produced during smoldering [START_REF] Mendes | Numerical simulation of reactive transfers in spouted beds at high temperature: application to coal gasification[END_REF] . However, the gas phase combustion submodel (EDC) is designed to oxidize all CO mixed with oxygen, including the CO produced by smoldering combustion. But in reality not all CO is oxidized, particularly that originated from smoldering, due to the low temperature. To separate between the two types of CO, an "inert" CO phase was introduced as a smoldering product.

α sg Q ( s ) char = [ ( 2 ϕ -1 ) h C O 2 + 2 ( 1 -ϕ ) h CO ] ˙ ω char
ϕ = 2 + CO / C O 2 2 CO / C O 2 + 2 (7)

Results and discussion

The estimation of the drag coefficient ( Eq. ( 2) ) affects the mass and the velocity of the flow that entered the sample. A misrepresentation of this value can lead to a misrepresentation of the total drag force resulting from the interaction between the solid phase and the gas phase. This results in incorrect heat transfer estimation, and poorly represented burning dynamics. Regarding the drag force coefficient estimation, the outcome of the aforementioned submodels ( Eq. ( 1) and ( 2) ) are presented in Fig. 4 . The Reynolds ( Re ) number was calculated using the kinematic air viscosity of dry air at 300 K (1.568 × 10 -5 m 2 s -1 ), resulting in Re LF == 2.3 and Re HF == 9.3. For Re > 4, both submodels provide similar estimations, but for lower Reynolds numbers ( Re < 4), corresponding to flows lower than 13.14 cm/s (or HF/2), the value of C D is higher using Clift et al. model for a sphere [START_REF] Clift | Bubbles, drops and particles[END_REF] . However, flow velocity fields were successfully simulated using the pseudofluid model at both low flow (LF) and at high flow (HF), matching the measured velocities on top, and around the sample for tests without combustion ( Fig. 5 a andb). An overestimation of the drag force leads to an underestimation of the flow in the fuel bed, which affect the burning dynamics through cooling, air mixing, combustion rate, char oxidation rate, among others. When modeling large scale forest fires, these correlations are usually applied to estimate the mean drag force generated from both wind/litter and wind/tree interactions in a control volume larger than 1 m × 1 m [START_REF] Linn | Studying wildfire behavior using FIRETEC[END_REF][START_REF] Morvan | Numerical simulation of the interaction between two fire fronts in grassland and shrubland[END_REF] . By doing so, the flow and consequently, the fire behavior of the litter are misrepresented. Hence, vegetation elements producing different drag coefficients have to be separated.

Once ignition occurred, the observed flame was mostly laminar at the base, as shown in Fig. 6 a. Transient behavior was also observed ( Fig. 6 b), but only in the intermittent zone of the flame, which was mostly due to the entrained air.

The difference between the three convective heat transfer coefficients outcomes presented in Fig. 7 χ > 40 W/m 2 K. However, this value is high and competes with the radiative heat transfer source term; thereby preventing the solid temperature from rising, degradation to occur, and ultimately ignition to happen. Moreover, these correlations were initially reported for moderate flow temperatures [START_REF] Irvine | Advances in heat transfer[END_REF] and were not verified at high temperatures.

On the other hand, coefficient proposed by Colburn et al .

[38] corresponded to χ ≤ 20 W/m 2 K, providing a more moderate heating/cooling rate. We can justify that these coefficients are more adapted to account for the sheltering effect produced when elements are close together and act as a bulk. Even though, this value is relatively low, it was demonstrated that the packing of needles generally causes the heat transfer coefficient to decrease [START_REF] Lamorlette | Analytical quantification of convective heat transfer inside vegetal structures[END_REF] .

The reason why the other correlations worked in previously cited studies is probably because the energy balance was compensated by an overestimation of the radiative heat transfer. In contrast, the radiative heat transfer of the FPA configuration was carefully examined and quantified [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] .

Figure 8 a illustrates the predicted temperature field of the corresponding experiment during flaming (HF, 40 kg/m 3 ). Temperature ranges were in the same order of magnitude as other flames from wildland fuels [START_REF] Dupuy | Fires from a cylindrical forest fuel burner: combustion dynamics and flame properties[END_REF][START_REF] Tihay | Experimental study of laminar flames obtained by the homogenization of three forest fuels[END_REF] . The shape of the simulated flame is consistent with the experimental one ( Fig. 6 a). The evolution of the temperature in the solid phase is presented in Fig. 8 b at the top of the sample, the middle and the bottom using the correlation proposed by Colburn et al . [START_REF] Incropera | Fundamentals of heat and mass transfer[END_REF] . The temperature on top, in the middle, and in the bottom increased at a slower rate than observed in NF conditions [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] before ignition due to the convective cooling. After the initial ignition, the middle temperature promptly increased. At 12 s, the middle temperature reached 800 K, followed by an increase in the temperature rate, corresponding to a localized ignition. As the flame propagated downwards, the bottom temperature increased until it was fully involved in the flame. Temperature measurements inside the fuel bed were not obtained in this study due to experimental limitations. Two main issues occurred with temperature measurements with forced flow. First, the vibration of the thermocouples inside the fuel bed had an influence on the pine needle arrangement in the sample, which it is important great care in having similar ones for every experimental procedure, in order to have better repeatable conditions. Second, the wiring of the thermocouples in the FPA combustion chamber obstructs the inflow. However, the simulated temperature evolution is consistent with temperature measurements conducted for Pinus strobus ( σ s == 10,788 m -1 ) in the FPA under similar testing conditions (heat flux of 50 kW/m 2 ) [START_REF] Simeoni | Flammability studies for wildland and wildland-urban interface fires applied to pine needles and solid polymers[END_REF] .

Measured and simulated flaming times are plotted in Fig. 9 for different bulk densities and inlet flows. Figure 9 shows that flaming time increased with increasing bulk density. Flaming time was longer as there was more fuel to burn. It also increased with decreasing flow velocity, and the flames were maintained the longest (45 s) under natural convection at maximum bulk density. The forced flow enhanced the mixing pyrolysis gases with air and increased the burning rate. For low bulk densities, the influence of the forced flow was not as dominant as for higher bulk densities: for all flow conditions, flaming times varied between 13 and 17 s. Since the sample was very porous (96%%), the induced air easily penetrated and provided enough oxygen to obtain well-ventilated combustion conditions. Flaming time linearly increased with bulk density for various flow rates since the air contribution is directly related to an increase in the combustion rate, as the latter becomes limited by the fuel available. In contrast, under natural convection the relation between flaming time and bulk density is non linear, because the combustion rate is limited by both the available fuel and the oxygen, especially for high bulk densities (92%%). For clarity, simulations are only plotted for NF and HF, which are the lower and upper bounding conditions. Also, good agreement was found for LF. Simulations slightly over predicted the flaming time by a few seconds. It is a direct consequence from using a highly efficient gas phase combustion model, where all available fuel was consumed when mixed with Oxygen. In reality the efficiency rate is not maximum (due to air dilution and cooling).

One of the main advantages of numerical simulations is the ability of separating the processes involved in the evolution of the solid phase, such as in Fig. 10 , in which, the measured and simulated time evolution of the mass sample (normalized by the initial mass) were plotted for a bulk density of 40 kg/m 3 and HF. The simulated mass loss was split into a dry fraction, a moisture content fraction, and a char fraction. Piloted ignition was observed experimentally. However, the pilot (in Fig. 6 a) was not simulated because the combustion model implies that flaming combustion occurs as soon as fuel and oxidizer are mixed. As ignition occurred faster (5 s) in the simulation, all the curves were shifted to synchronize both ignition times. This is not surprising because piloted ignition is a marginal event that any small variation in the experiment or in the numerical condition can influence it. More importantly, we are interested in comparing the combustion phases following ignition. One can observe that the simulated (total) mass loss followed the mean experimental mass loss closely. The different dotted lines correspond to the three simulated fractions, giving more insight on the mass loss mechanisms. In the beginning, the mass fraction was only composed of moisture content and dry fraction. Then the moisture content evaporated during the first 13 s of the simulation. More dry material was lost (via pyrolysis) when ignition occurred, during which the dry fraction decreased steeply, and the char fraction formed (corresponding to flaming). The transition to smoldering combustion occurred after 30 s (near flameout time), during which the char oxidation reaction was dominant until the end.

Mass loss rates (MLR) were calculated from experimental and numerical mass loss data and are shown in Fig. 11 for two bulk densities: 23 kg/m 3 ( Fig. 11 a) and 40 kg/m 3 ( Fig. 11 b) at HF. A 5 s moving average was used to smooth experimental curves. In the first 10 s of both cases, the smaller slope corresponds to pyrolysis and evaporation that occurred before ignition. When ignition occurred, the slope increased steeply. The subsequent decrease corresponds to flameout (vertical line around 30 s), where after mainly char oxidation occurred. The numerical results over predicted the maximum value by around 20%%, but followed the same trend as in the experiments under both conditions.

The radiation attenuation coefficient can drastically change influence the radiation distribution in the fuel depth and affect the distribution of the mass loss rate. In these simulations, the extinction coefficient was kept as 4/ α s σ s , similarly as in a previous study [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] . However, it was demonstrated in Acem et al . [START_REF] Acem | Analytical determination and numerical computation of extinction coefficients for vegetation with given leaf distribution[END_REF] that the estimation of this theoretical coefficient could be could be affected by the geometrical shape by 10%%.

For better comparison between the different bulk densities and flow conditions, MLR were averaged for each configuration, and are presented separately during flaming and after flameout in Fig. 12 a and b, respectively. During flaming ( Fig. 12 a), one can observe that the average MLR increased with an increasing inlet flow for all tested conditions because, generally the flow enhances the mixing of pyrolysis gases with air, increasing the combustion rate. The following observations can be made:

-Under HF conditions, the average value increased almost linearly with the bulk density. This is a result of the direct influence of the airflow on the combustion rate. -For NF conditions, the average MLR increases at 30 kg/m 3 , because there is more fuel to burn. However, it decreased at the highest bulk density, mostly due to the low ventilation conditions induced by natural convection, which slowed down the reaction. -As for LF condition, the mean MLR remains non-monotonic as a result of the competition between radiation heating more of the solid phase with increasing bulk density (increasing α s ), and the decreasing oxygen available to sustain the reaction. -At the lowest bulk density (23 kg/m 3 ), which corresponds to the highest porosity ( Table 1 ), the mean MLR is similar regardless of the flow conditions, because air can easily penetrate and provide well-ventilated combustion conditions.

-For higher bulk densities, the sample was more compact, causing more contact between pine needles, and more shadowing effect. Therefore, there was less contact with the gas phase. As the extinction length ( ∼ 4/ α s σ s ) [START_REF] Albini | A model for fire spread in wildland fuels by radiation[END_REF] changes from 13.7 mm to 7.8 mm, between 23 and 40 kg/m 3 , respectively, radiation cannot penetrate deeper in the sample body. The optical thickness defined as the ratio between the depth of the sample (3 cm) and the extinction length shows that the sample is optically thick (2.19 and 3.84 for 23 and 40 kg/m 3 , respectively). -When the total amount of radiation is distributed on a shallow layer the limited factor is the amount of fuel available for pyrolysis, However, when the same amount of energy is distributed deeper in the fuel (i.e. low bulk density), more mass undergoes pyrolysis. Ultimately, pyrolysis of a denser fuel will be limited by the heat transferred downwards, which is correlated to the bulk density.

-Between experiments and simulations with and without flow, we can highlight the existence of two regimes. The oxygen limited oxygen regime (no wind) and the fuel limited regime (with wind) [START_REF] Morvan | Modeling of one dimensional fire spread in pine needles with opposing air flow[END_REF] . Because the mean MLR is not a monotonic function of the density in NF and LF conditions, we can determine that the combustion is limited by oxygen supply (via saturation).

Overall, simulations are consistent with measurements, regardless of the peak values that are slightly overestimated numerically, as shown in Fig. 11 a and b, and in [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] .

Concerning the average MLR after flameout presented in Fig 12 b, we can safely assume that mass loss is driven by smoldering. Since the latter occurs at very slow reaction rate, the average MLR is one order of magnitude smaller than the one of flaming, which makes it difficult to distinguish tendencies. Globally, one can observe that the average MLR increases with bulk density and with the flow, because there is more fuel and more air, respectively. This is consistent with Eq. ( 5), where the char oxidation rate is proportional to the oxygen concentration and to the bulk density. The numerical predictions in HF compared well with measurements, as a result of the proper flow and heat transfer estimation. However, results under NF conditions did not match the experimental results because the char oxidation was not sustained due to the low air supply in natural convection, as explained in [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] . Hence, the reaction only occurred in the few seconds after flameout then dropped, resulting in a high average MLR compared to the experiments, which lasted longer. The smoldering time in NF could be enhanced by either testing a more comprehensive model for estimating the char combustion rate, or by accounting for the increase of the surface to volume ratio in the solid phase (which is constant) as a result of the char pores that are formed in the fuel.

The measured and simulated Heat Release Rate (HRR) are presented in Fig. 13 for two bulk densities at HF. Computed curves were time-shifted so that they coincide with experimental curves when ignition occurred. The peak HRR was slightly over estimated for both cases, a result of the high MLR peaks. Nevertheless, the overall trend is matched and the total heat released is similar. The peak HRR at 40 kg/m 3 ( Fig. 13 b) is at approximately 11 kW, which is comparable to results in [START_REF] Bartoli | Determination of the main parameters influencing forest fuel combustion dynamics[END_REF][START_REF] Simeoni | Flammability studies for wildland and wildland-urban interface fires applied to pine needles and solid polymers[END_REF] for Aleppo pine ( Pinus halepensis ) tested under similar conditions. The transition between flaming and smoldering (at flameout) was better predicted numerically for the higher bulk density ( Fig. 13 b). However, for the lower bulk density ( Fig. 13 a), an abrupt transition can be noticed only in the simulation. This behavior can be justified by the larger amount of fuel burning (i.e. higher solid fraction) available to increase the heat transfer from solid to solid, and to better sustaining char oxidation after flameout.

For better illustration of the influence of the porosity in the HRR, averaged values for measured and simulated HRR were normalized over their corresponding initial fuel mass ( HRR / m 0 ) and are presented during flaming and after flameout, in Fig. 14 a and b, respectively. The normalization allows better comparison of the energy rate released between different bulk densities. In Fig 14 a, HRR / m 0 globally increases with the flow. This is consisted with the observed MLR ( Fig. 12 a) and with the shorter flaming times for higher flow conditions. Essentially, Higher HRR is reached for shorter flaming time due to the better mixing, and consequently the enhanced the combustion rate. The observed HRR / m 0 trends to decrease with an increasing bulk density, which is also consistent with the aforementioned observations about the competition between radiation penetration and oxygen availability. As for the numerical predictions, these tendencies are well matched at high flow but cannot be obtained under NF conditions due to the low combustion rate estimated [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] .

Regarding HRR / m 0 after flameout ( Fig. 14 b), overall the measured values are also an order of magnitude lower than during flaming. They decrease with an increasing bulk density, since more air can reach the reacting solid phase for lower bulk densities. HRR / m 0 is consistently slightly higher at HF, especially for the higher bulk density. Simulations of HF conditions matched the measurements. However, for NF condition, a high average MLR was calculated after flameout ( Fig. 12 b). Consequently, the calculated HRR / m 0 was overestimated compared to measurements only because a higher HRR was released over a shorter period of time. However, the total energy released is smaller, since not all the char was oxidized in NF.

Figure 15 a, b, and c describe the evolution of the fuel and oxygen mass fractions inside the fuel bed under HF conditions at ignition, 5 s and 15 s after ignition, respectively. It can be noted that the propagation occurs from top to bottom, which is in agreement with the experimental observations. Moreover, the evolution of the mass fractions (step) are representative of a typical diffusion flame, meaning that no partial mixing in involved. It is clear that the oxygen is completely consumed in the first cell of the reaction zone.

On the other hand, Fig. 15 d,e, and f represent the evolution of the mass fraction in NF conditions. In Fig. 15 d, the conditions at ignition are similar to those of encountered with HF ( Fig. 15 a). As the flame propagates slower in NF, Fig. 15 e and f show the mass fractions 20 and 40 s after ignition. It appears that the entire fuel sample is embedded in the fuel rich zone and that there is no sufficient oxygen inside it to ensure flaming combustion. Anecdotal observations of experiments using pine needles in the cone calorimeter under natural convection allowed to see the flame burning below the sample, which is characteristic of the behavior of a fuel rich diffusion flame and confirms the prediction from Fig. 15 e andf. This behavior was also observed in the FPA in natural convection but it was harder to see because of the light of the heaters. the transition from flaming to smoldering (second phase) occurs in a shorter time than in forced flow ( Fig 16 ). This behavior is mainly due to the additional oxygen supplied inside the fuel bed by the forced flow, which allows flaming and smoldering to occur simultaneously. In contrast, smoldering occurs mainly after flameout in natural convection. The differences in CO curves between HF and NF ( Figs 16 and17 ) demonstrate the changing in behavior in the combustion process. CO concentration is an indicator of the two stages of combustion when correlated to visual observations. The first step increase in Fig. 17 is due to the ignition of the sample. A steady production of CO follows, which corresponds to the flaming stage (first phase) with possibly an oxidation of CO in the flame. The consequent increase is due to the decrease of the flame and the initiation of smoldering combustion. The two processes overlap during this stage (second phase). After flameout, CO production peaks due to the sole occurrence of smoldering combustion (third phase). In contrast, the steady state of the CO disappears in HF. A faster combustion time was observed with embers starting to glow before flameout. These tendencies confirm the measurements and the observations described by Schemel et al . [START_REF] Schemel | A calorimetric study of wildland fuels[END_REF] .

Concerning the numerical predictions, all the main tendencies were captured qualitatively by the model, even though the quantitative predictions were not accurate. In both configurations, the simulations overestimated the CO 2 production during flaming because the entire pyrolysis products were consumed in the gas phase due to the nature of the combustion model. In Fig. 16 , smoldering was not activated during flaming, therefore, no CO was detected in the simulation. In reality, not all CO from pyrolysis is entirely oxidized, due to dilution with air and cooling. Indeed, the concentration of measured CO increased as the flame was quenching and more smoldering was occurring. An extinction model could be added to account for the unburned gases and better represent the incomplete combustion in the flame, especially at a larger scale [START_REF] Vilfayeau | Numerical simulation of underventilated liquid-fueled compartment fires with flame extinction and thermally-driven fuel evaporation[END_REF] . However, its absence did not severely influence the burning dynamics of the pine needle bed (which was the focus of this work), as it was mainly driven by the incident radiation coming from the lamps. In Fig. 17 , numerical CO is detected even before flame out time, which is due to an early activation of the smoldering process, as it was observed experimentally. Hence, both burning dynamics in NF and HF are captured by the model. The experimental peak value for CO was 25 ppm for HF, whereas it was 35 ppm for the simulation. This is also due to the combustion model that overestimated the combustion efficiency in the gas phase and had little influence on the burning dynamics. An adjustment will certainly be needed for large fires, where extinction is likely play a role, especially if the model is to be used to estimate fire emissions.

Conclusion

The goal of the presented framework was to assess the relevance and the performance of the submodels used to close CFD models particularly with the multiphase approach for wildfires. Some submodels were successfully adapted in a specific and realistic range of conditions, which allowed improving burning rate estimations and better describing the underlying physics. Moreover, the importance and the implications of using appropriate submodels were also demonstrated. By comparing simulations to experiments with various fuel bulk densities and various inlet flows, at a constant applied heat flux with a well defined fuel, it was shown which submodels need to be appropriately defined in order to provide acceptable predictions. Because of the strong coupling between the different submodels, the proposed framework of using controlled experiments and matching simulations at a small scale is necessary for choosing the adequate submodels. This will ensure that the physics involved in the burning of the fuel are correctly captured before predicting fire spread at a larger scale. The main results are summarized as follows:

• The estimation of the drag force coefficient affects the flow field inside the porous fuel, consequently it affects the associated burning dynamics during flaming and smoldering.

The proposed correlations can be used for subgrid modeling and included in large scale simulations to represent the drag force through litter beds separately from the one for trees. • For the convective heat transfer coefficient, we have verified that Colburn coefficients are appropriate for these experimental conditions, whereas other cited coefficients can prevent the fuel from igniting. • Under forced flow conditions, the char oxidation model based on a single step Arrhenius equation is sufficient for characterizing smoldering combustion, and especially mass loss after flameout. However, this model is not adapted under natural convection, and further investigation is necessary. • The additional split function in the char oxidation allows predicting an acceptable gas emission and adequately modeling the transition between flaming and smoldering emission. • Overall, the selected gas phase combustion model allowed a successful prediction of the burning dynamics of the solid phase in the chosen range of conditions. However, it needs further improvements for supporting emission estimations.

Appendix. Governing equations in FFF

The gas phase is governed by a set of transport equations representing the mass balance equations:

∂ α g ρ ∂t + ∂ α g ρ ˜ u j ∂ x j = ( 1 -ν ash ) ˙ ω char + ˙ ω v ap + ( 1 -ν char ) ˙ ω pyr ( 9
)
where ∼ is the Favre filter operator, ν char and ν ash are the char and ash fractions, respectively. For forest fuel ν char ∼ 0.3 and ν ash ∼ 0.05 [START_REF] Grishin | A mathematical modelling of forest fires and new methods of fighting them[END_REF][START_REF] Kollman | Principles of wood science and technology[END_REF] . The momentum equation is defined as: [START_REF] Olson | Experimental observations of spot radiative ignition and subsequent three-dimensional flame spread over thin cellulose fuels[END_REF] ν is the molecular viscosity, ν t is the sub grid scale viscosity, and p is the pressure. F D represents the drag force source term resulting from the interaction between the gas flow and the solid phase, and is defined as:

∂ ρ ˜ u i ∂t + ∂ ( α g ρ ˜ u i ˜ u j ) ∂ x j = - ∂ p ∂ x i + ∂ ∂ x j × α g ρ( ν+ ν t ) ∂ ˜ u i ∂ x j + ∂ ˜ u j ∂ x i - 2 3 ∂ ˜ u k ∂ x k δ i j + α g ρg i -F D
F D = α g ρC D α s σ s 2 | u | u i ( 11
)
Since thermal equilibrium is not assumed between the solid fuel particle and gaseous phase, the temperature in the solid phase is solved separately, in the following equation:

C ( S ) p α s ρ s d T s dt = Q ( S ) ± -h v ap ˙ ω v ap -h pyr ˙ ω pyr -α sg Q ( S )
char [START_REF] Boonmee | Glowing ignition of wood: the onset of surface combustion[END_REF] with Q (S) ± the energy balance on the solid phase exchanged with the gaseous phase by convection and radiation. C p (S) is the specific heat capacity of the solid phase, h char , h pyr , h vap are the heat of reaction for charring, pyrolysis, and evaporation, respectively. The time evolution of the fuel is characterized by the variation of its dry, water, and char mass fractions. They can be described by the following three ordinary differential equations: ω char [START_REF] Grishin | A mathematical modelling of forest fires and new methods of fighting them[END_REF] The global mass balance equation for the solid phase is: d ( α s ρ s ) dt = ( ν char -ν soot -1 ) ˙ ω pyr -˙ ω char -˙ ω v ap [START_REF] Torero | Heat and mass transfer in fires: scaling laws, ignition of solid fuels and application to forest fires[END_REF] Assuming that the solid consumption is only due to char combustion. The balance equation for the solid fraction is resolved following:

d α s ρ s ϕ ( S )
d α s dt = - ˙ ω char ρ s (17)
Finally, the energy balance of the gas phase is written as:

∂ α g ρ ˜ h ∂t + ∂ α g ρ ˜ u j ˜ h ∂ x j = D p Dt + ∂ ∂ x j α g ρ α D + ν t P r t ∂ ˜ h ∂ x j -Q ( S ) CONV + Q RAD + Q COMB -( 1 -α sg ) Q ( s ) char ( 18
)
where h is the enthalpy, α D is the thermal diffusivity (considering a unity Lewis number approximation). Q COMB and Q RAD are the source term for combustion and radiation in the gas phase, respectively, as described in [START_REF] Houssami | Experimental and numerical studies characterizing the burning dynamics of wildland fuels[END_REF] .

Fig. 1 .

 1 Fig.1. Overview of the Fire Propagation Apparatus (FPA).

Fig. 2 .

 2 Fig. 2. Pitch pine needle samples with (a) 23 kg/m 3 ; (b) 30 kg/m 3 ; (c) 40 kg/m 3.

Fig. 3 .

 3 Fig. 3. Computational domain of the FPA (not to scale).

  pine needles over each other, one can assume that the formulation of the convective heat transfer coefficient is similar to that of array of staggered cylinders in cross flow. Many values for C, m , and n are proposed to represent the Nusselt ( Nu ) number in these specific or similar conditions[START_REF] Irvine | Advances in heat transfer[END_REF][START_REF] Khan | Convection heat transfer from tube banks in crossflow: analytical approach[END_REF] . Values for C, m , and n were proposed by Žukauskas[START_REF] Irvine | Advances in heat transfer[END_REF] ( C == 1.04; m == 0.4; n == 0.36), Colburn[START_REF] Incropera | Fundamentals of heat and mass transfer[END_REF] ( C == 0.33; m == 0.6; n == 1/3), Incropera and De-Witt[START_REF] Incropera | Fundamentals of heat and mass transfer[END_REF] ( C == 0.683; m == 0.466; n == 1/3), and others[START_REF] Khan | Convection heat transfer from tube banks in crossflow: analytical approach[END_REF] for flow across banks of staggered tubes, for 10 or more rows of tubes and for different Reynolds regimes.

Fig. 4 .

 4 Fig. 4. Drag force coefficient estimation using different submodels.

  is non negligible. For a typical HF inlet, the measured flow inside the fuel bed is approximately 0.2 m/s, corresponding to Re HF == 9.3. The correlations proposed by Zukauskas et al. and DeWitt et al. give convective coefficients

Fig. 5 .Fig. 6 .

 56 Fig. 5. Measured and simulated air velocity on the peripheral free space (gap) and on top of the fuel (sample) in the FPA at HF for (a) 23 kg/m 3 and (b) 40 kg/m 3 .

Fig. 7 .

 7 Fig. 7. Convective heat transfer coefficient estimation.

Fig. 8 .

 8 Fig. 8. (a) Temperature field during flaming (b) Simulated solid temperature at the top, middle and bottom of the sample corresponding to the three blue sampling points in (a), for HF and 40 kg/m 3 . Vertical line marks ignition. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

Fig. 9 .

 9 Fig. 9. Measured and simulated flaming time for different bulk densities and inlet flow. Filled markers: experiments, empty markers: simulations, error bars: standard deviations for experiments.

Fig. 10 .

 10 Fig. 10. Experimental and simulated mass loss, including simulated mass loss of dry, water, and char phases. Vertical dotted line corresponds to experimental and simulated ignition times. Solid vertical lines correspond to simulated and experimental flameout times respectively. (HF -40 kg/m 3 ).

Fig. 11 .

 11 Fig. 11. Mass loss rates (MLR) at HF for different bulk densities (a) 23 kg/m 3 and (b) 40 kg/m 3 . Error bars: standard deviation. Vertical line: flameout time.

Fig. 12 .

 12 Fig. 12. Averaged values for measured and simulated mass loss rates (a) during flaming (b) after flameout.

Fig. 13 .

 13 Fig. 13. Measured and simulated Heat Release Rates at HF for different bulk densities (a) 23 kg/m 3 and (b) 40 kg/m 3 . Error bars: standard deviation. Vertical line: flameout time.

Fig. 14 .

 14 Fig. 14. Averaged values for measured and simulated heat release rates (HRR) per unit of initial fuel mass (a) during flaming (b) after flameout.

Fig. 15 .

 15 Fig. 15. Fuel and oxygen mass fractions along the vertical and at the center of the sample (hatched area) during flaming. Figures a, b, and c correspond to distributions at ignition, 5 s and 15 s, respectively for HF conditions. Figures d, e, and f correspond to distributions at ignition, 20 s and 40 s, respectively for NF conditions.

  Fig. 15. Fuel and oxygen mass fractions along the vertical and at the center of the sample (hatched area) during flaming. Figures a, b, and c correspond to distributions at ignition, 5 s and 15 s, respectively for HF conditions. Figures d, e, and f correspond to distributions at ignition, 20 s and 40 s, respectively for NF conditions.

  Fig. 15. Fuel and oxygen mass fractions along the vertical and at the center of the sample (hatched area) during flaming. Figures a, b, and c correspond to distributions at ignition, 5 s and 15 s, respectively for HF conditions. Figures d, e, and f correspond to distributions at ignition, 20 s and 40 s, respectively for NF conditions.

Figures 16 and 17

 17 present the experimental and simulated time evolution of CO and CO 2 generation for a fuel bulk density of 40 kg/m 3 in HF and NF, respectively. The experimental values were measured in the exhaust duct, similarly the simulated values were obtained at the outlet boundary. Useful information can be extracted from these figures and three phases can be defined. During the first phase, which occurs just after ignition, only flaming is involved as CO 2 production increases. During the second phase, flaming and smoldering occur simultaneously. Finally, during the third phase which is after flameout, only heterogeneous combustion occurs. By comparing Figs 16 and 17 , it can be noted that the time of each phase can be very different. In natural convection ( Fig 17 ),

Fig. 16 .

 16 Fig. 16. Measured and simulated CO and CO 2 production (HF -40 kg/m 3 ). The vertical lines represent the experimental (red dashed) and the numerical (black solid) flameout times. Plots are shifted to have ignition at 0 s. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

Fig. 17 .

 17 Fig. 17. Measured and simulated CO and CO 2 production (NF -40 kg/m 3 ). The vertical lines represent the experimental (red dashed) and the numerical (black solid) flameout times. Plots are shifted to have ignition at 0 s. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

d

  α s ρ s ϕ ( S ) char dt = ( ν char -ν soot ) ˙ ω pyr -ν ash ν char + 1 ˙

  , during

	Nomenclature
	C		Convective heat transfer constants
	C D		Drag force coefficient
	D		Equivalent diameter
	E CHAR	Energy activation for char oxidation
	HRR	Average Heat Release Rate
	k		Air thermal conductivity
	K CHAR	pre-exponential coefficient for char oxidation
	m 0	Initial mass
	m, n	Convective heat transfer coefficient constants
	Nu		Nusselt number
	P		Pressure
	Pr Q (s ) char Q (S) CONV	Prandtl number Char oxidation source term Convective heat transfer source term
	R		Ideal gas constant
	Re		Reynolds number
	T, T s	Gas and solid phase temperatures
	U		Velocity
	Y O 2	Oxygen mass fraction
	Greek symbols
	α g , α s		Volume fraction of gas and solid phases
	ρ	h CO , h C O 2	Heat of combustion of CO and CO 2 Gas density
	σ s			Surface to volume ratio
	ϕ ϕ (S)		Split function

char , ϕ

(S)

DRY , ϕ (S) H 2 O Mass fraction of char, dry pine and water in solid phase χ Convective heat transfer coefficient ˙ ω char , ˙ ω pyr , ˙ ω v ap

Table 1

 1 Summary of bulk densities and porosities.

	Mass [g]	Bulk density [kg/m 3 ]	Porosity αg [-]
	8.7	23	0.96
	11.4	30	0.95
	15.0	40	0.93

Acknowledgments

The authors gratefully acknowledge the Fire Dynamics Group at FM Global for the help provided in the development of Forest-FireFOAM. Insights offered by G. Accary are most appreciated. This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso ( ANR-10-EQPX-29-01 ) of the program « Investissements d'Avenir » supervised by the French Agence Nationale pour la Recherche.


Stanislas Rebaudet

Karolina Griffiths
email: karolinagriffiths@hotmail.com

Mazard Trazillio

Anne-Gaelle Lebeau
email: aglebeau@unicef.org

Aaron A Abedi

Gregory Bulit
email: gbulit@unicef.org

Renaud Piarroux
email: renaud.piarroux@ap-hm.fr

Jean Gaudart
email: jean.gaudart@univ-amu.fr

Cholera spatial-temporal patterns in Gonaives, Haiti: from contributing factors to targeted recommendations

Keywords: Cholera, Socio-economic risk factors, Environmental risk factors, Spatial analysis, Water availability, Sanitation

Introduction: Gonaives, in the Artibonite region, is one of the most affected areas by the cholera epidemic in Haiti. Five years on, the epidemic persists and further information is needed to guide water and sanitation strategies for cholera elimination.

Methods: The number of cholera cases, socio-economic and environmental characteristics of Gonaives city neighbourhoods were included in an unsupervised classification approach and a general additive model was used to estimate standardized incidence ratios, adjusted on locations.

Results: Five classes of neighbourhoods were identified in Gonaives. Three classes were defined as high-risk in comparison with the area with the lowest incidence ratio: coastal flood-prone poor neighbourhoods (Class 5) with frequent open defecation and wide access to a shallow water table 7.04]); a coastal trade area (Class 3) with good access to a brackish and shallow water table likely contaminated by numerous latrines (SIR of 2.5 IC95%[1.05; 5.95]); and the hillside poor residential neighbourhoods (Class 4) with poor water availability and frequent open defecation (SIR=3.2 IC95% [1.25; 8.18]).

Conclusion:

Water and sanitation prevention strategies should not be homogenous. Within an urban area, prevention strategies should be planned according to the specific socio-environmental and geographical context of each neighbourhood.

 Hillside areas are high-risk associated with poverty, poor water availability and open defecation.

 Development of water supply networks and safe sanitation practices must be neighbourhood-specific.

A C C E P T E D M

A N U S C R I P T

Introduction

Since October 2010, Haiti has been affected by the most important cholera epidemic in recent history. A dramatic number of cases were seen with the first phase of the epidemic in October-December 2010, with daily case counts exceeding 4 000 in mid-December [START_REF] Barzilay | Cholera Surveillance during the Haiti Epidemic -The First 2 Years[END_REF]. This was followed by further large peaks of cases in mid-2011, late-2011 and in 2012. Despite the overall trend of decreasing number of cases of cholera over the past five years, cholera persists with epidemic waves in the rainy season separated by residual cases in between peaks. As of the 14 th April 2016, a total of 776 661 cases and 9 287 deaths have been reported in Haiti (MSPP, 2016(MSPP, , p. 2016)).

It is now well established that the cholera epidemic originated from a camp of UN Nepalese peacekeepers near Mirebalais, in the Centre Department, and initially spread along the Artibonite river, both the longest river in Haiti and an essential water source for local communities [START_REF] Piarroux | Understanding the cholera epidemic, Haiti[END_REF][START_REF] Cravioto | Final Report of the Independent Panel of Experts on the Cholera Outbreak in Haiti[END_REF][START_REF] Hendriksen | Population genetics of Vibrio cholerae from Nepal in 2010: evidence on the origin of the Haitian outbreak[END_REF][START_REF] Frerichs | Nepalese origin of cholera epidemic in Haiti[END_REF][START_REF] Gaudart | Spatio-Temporal Dynamics of Cholera during the First Year of the Epidemic in Haiti[END_REF].

The cholera epidemic quickly propagated from the Centre and Artibonite departments, with rapid transmission country wide over several weeks in October-November 2010 [START_REF] Gaudart | Spatio-Temporal Dynamics of Cholera during the First Year of the Epidemic in Haiti[END_REF][START_REF] Blackburn | Household-Level Spatiotemporal Patterns of Incidence of Cholera, Haiti, 2011[END_REF][START_REF] Allan | High-resolution spatial analysis of cholera patients reported in Artibonite department, Haiti in 2010-2011[END_REF]. Not only was the Artibonite department among the first hit, it also continues to be one of the most affected. After the capital, Port-au-Prince, the administrative department of Artibonite saw the highest number of cases (119 519, as of 14 th April 2016) (MSPP, 2016(MSPP, , p. 2016)). Gonaives, with an estimated population of around 356 000 (IHSI (Institut Haitien de Statistique et d 'Informatique) et al., 2015, p. 20156), is both the principal town of the Artibonite department and an area severely hit by cholera. A retrospective [14.9-24.4], an estimated two-fold increase on previously reported mortality rates (Luquero et al., 2016).

Previous studies of cholera in urban settings have identified patterns in environmental and socio-economic characteristics associated with cholera cases, which include socio-economic status, such as income level [START_REF] Penrose | Informal urban settlements and cholera risk in Dar es Salaam, Tanzania[END_REF] and hygiene conditions, such as proximity to refuse dumps [START_REF] Osei | Spatial dependency of V. cholera prevalence on open space refuse dumps in Kumasi, Ghana: a spatial statistical modelling[END_REF]. Other studies have analysed the protective effects of improved sanitation, such as safe water supply sources, hand-washing with soap [START_REF] Kanungo | Clinical, epidemiological, and spatial characteristics of Vibrio parahaemolyticus diarrhea and cholera in the urban slums of Kolkata, India[END_REF] and water treatment [START_REF] Dunkle | Epidemic cholera in a crowded urban environment, Port-au-Prince, Haiti[END_REF] at the individual level.

Environmental characteristics, such as proximity to lakes, have been shown to be significant in a spatial analysis [START_REF] Bompangue | Cholera epidemics, war and disasters around Goma and Lake Kivu: an eight-year survey[END_REF], although few studies use geographical coordinates, especially within an urban area. Studies have highlighted this importance of research at the local level, where there can be large differences of cholera transmission within a city, between different communes or informal settlements [START_REF] Penrose | Informal urban settlements and cholera risk in Dar es Salaam, Tanzania[END_REF][START_REF] Azman | Urban cholera transmission hotspots and their implications for reactive vaccination: evidence from bissau city, Guinea bissau[END_REF].

Studies at the city level may hide intra-urban differences in social and environmental conditions [START_REF] Sur | Factors associated with reported diarrhoea episodes and treatment-seeking in an urban slum of Kolkata, India[END_REF][START_REF] Sasaki | Spatial analysis of risk factor of cholera outbreak for 2003-2004 in a peri-urban area of Lusaka, Zambia[END_REF]. Furthermore, it is essential to identify the areas of a city most affected, as they may act as drivers of the epidemic across the city and should be targeted first for control strategies and vaccination programmes [START_REF] Azman | Urban cholera transmission hotspots and their implications for reactive vaccination: evidence from bissau city, Guinea bissau[END_REF]. This therefore indicates the need to study the association of environmental and socio-economic At the request of UNICEF in 2013, the authors have been involved in epidemiological studies of the cholera outbreak, helping to understand the dynamics of cholera by bringing together the environmental context and population characteristics with a spatial-temporal analysis. Here, we present an epidemiological study on cholera in Gonaives, the objective of which is to provide information on the socio-economic and environmental characteristics associated with the number of cases of cholera and describe the profiles of the neighbourhoods most affected by cholera. This will help guide the middle and long-term prevention strategies to reduce the vulnerability of the population and reduce the risk of cholera transmission.

Gonaives: situation and characteristics

The commune of Gonaives covers an area of 620.5 km 2 . Over three quarters of its 356,000 estimated inhabitants (IHSI (Institut Haitien de Statistique et d'Informatique) et al., 2015) live in the coastal city, which represents the third urban concentration of Haiti, and is located on the main road connecting Port-au-Prince Metropolitan Area to Cap-Haitien. Local rainfall is scarce in the commune, but the city is situated at the bottom of a large watershed, where runoff water collects and drains to the sea. Due to these geographical characteristics, Gonaives presents an increased vulnerability to hydro-meteorological events, such as floods and cyclones in comparison to other surrounding areas. The city of Gonaives is a major transit route for people travelling to and from the south of Haiti. Alongside the cholera epidemic, other recent socio-economic concerns for the city include the displacement of people after the 2010 earthquake, riots, and increasing populations in shantytowns with poor living conditions. The commune of Gonaives is subdivided into Gonaives city plus five communalsections -1 st Pont Tamarin, 2 nd Bassin, 3 rd Petite Rivière de Bayonnais, 4 th Poteaux, 5 th Labranle (Figure 1 Panel B). Gonaives city has spread rapidly over the past decades, and is currently split into 28 neighbourhoods (Supplementary Database), of which 20 of the neighbourhoods have been partially mapped by Gonaives city council (Supplementary Figure S1).

Data collection

The remaining eight neighbourhoods that had not been previously mapped were delimited by field investigations, through interviews and confirmation with local habitants (Supplementary Database). Addresses of cholera suspected cases seen

Statistical analysis

Univariate analysis was performed to analyse the relationships between cholera incidence and the collected socio-economic and environmental characteristics, as listed in the supplementary Database, using the non-parametric Kruskal-Wallis test. Continuous variables were transformed into categorical using quartile discretization. In order to classify neighbourhoods according to their characteristics, hierarchical classification on principal components of a multiple correspondence analysis (MCA) was used, as previously described [START_REF] Lê | FactoMineR: An R Package for Multivariate Analysis[END_REF][START_REF] Husson | Principal component methods -hierarchical clustering -partitional clustering: why would we need to choose for visualizing data?[END_REF]. All the variables were included in this analysis. A multiple correspondence analysis (MCA) was the first pre-processing step to explore the socio-economic and environmental characteristics on the mixed data set, which takes into consideration relationships between variables. Then, the coordinates of each variable in the first 20 principal components, which summarise 99.22% of the information, were used to perform a Hierarchical Ascendant Classification (HAC). This method provides classes according to the socio-economic and environmental profile using an objective non-supervised classification technique, allowing the classification to be independent from the number of cases of cholera. Furthermore, using the first 20 principal components

A C C E P T E D M A N U S C R I P T
10 is a way of removing the noise from the data and therefore renders the classification more stable [START_REF] Husson | Principal component methods -hierarchical clustering -partitional clustering: why would we need to choose for visualizing data?[END_REF]. A sensibility analysis was performed in order to evaluate the stability of the classification through different component numbers (15 to 25). A generalized additive model (GAM) was then applied to compare these classes with the number of cholera cases:

Log(CholeraCase) = Cluster + s(Lon,Lat )+ offset(log(Pop)) ~ NegBin
A Negative binomial model was performed, taking into account over-dispersion.

The log population was used as an offset in order to estimate the standardised incidence ratio associated with each class. A bivariate thin plate spline, s(Lon,Lat), based on geographical coordinates of each neighbourhood centroid was used to take spatial autocorrelation into account, as previously described by Simon Wood [START_REF] Wood | Generalized Additive Models: An Introduction with R[END_REF]. The level of significance was fixed at α=0.05.

Statistical analysis was performed using the software programmes R v3.3.0 (The R Foundation for Statistical Computing, Vienna, Austria), and the mgcv and FactoMineR packages. GIS analyses and mapping was performed using QGIS v2.12.1 (QGIS Development Team, 2016. QGIS Geographic Information System.

Open Source Geospatial Foundation Project. http://qgis.osgeo.org). Background satellite images used in the maps were obtained from Google Maps. Figures were formatted with GIMP v2.8.16 (The GIMP Team, 2016, www.gimp.org).

Ethical considerations

Cholera counts and addresses were anonymously extracted from registries and aggregated at the neighbourhood scale. This study was approved by the Haitian Ministry of Public Health and Population. and March 2016 (Supplementary Figure S2). However, spatial distribution of cholera was not homogeneous, as 50% of neighbourhoods provided less than 20, or more than 104 suspected cases (Supplementary Database & Supplementary Figure S2). The HAC on MCA (Supplementary Results & Figure 2) provided five classes for socio-economic and environmental profile characteristics. The first class (Class 1) was characterized by residential neighbourhoods of intermediate altitude (mean altitude=23m), a higher economic status, the use of latrines and modern toilets and no marketplaces. It corresponded to seven neighbourhoods at the city periphery, with an incidence ratio of 56.1 cases /10,000 inhabitants (min-max, 18.5-144.0) (Figure 3).

Results

A C C E P T E

The second class (Class 2) was characterized by a lower altitude (mean=9m) but no flooding risk, an intermediate distance to the sea (mean=1,809m), a shallow and fresh water table, an intermediate economic status, a main use of latrines, and a moderate water availability with a functional DINEPA network. It

A C C E P T E D M

A N U S C R I P T 14 represented 10 neighbourhoods in the middle of the city, with a mean incidence ratio of 61.0 cases /10,000 inhabitants (min-max, 6.9-132.1) (Figure 3).

The third class (Class 3) was characterized by a low altitude (mean=6m), a short distance to sea (mean=614m), a poor sanitation with a major use of open defecation and unimproved latrines, a high water availability with access to a shallow and brackish water table, a residential and commercial use of mainly two floors buildings, and large markets. It represented 3 coastal neighbourhoods at the city centre, with a mean incidence ratio of 104.0 cases/10,000 inhabitants (min-max, 88.2-114.3) (Figure 3).

The fourth class (Class 4) was characterized by a higher altitude (mean=41m), a low water availability with a deep water table, a low socio-economic status, and a prevailing use of open defecation. It represented four hillside neighbourhoods at the city northern periphery, with a mean incidence ratio of 166.9 cases /10,000 inhabitants (min-max, 93.4-546.1) (Figure 3).

The fifth class (Class 5) was a flood-risk area characterized by a low altitude (mean=5m), a short distance to sea (mean=631m), a high water availability with a shallow water table, a poor economic status, a poor sanitation with a major use of open defecation and unimproved latrines. It represented 4 coastal neighbourhoods divided in two part (north and south), with 119.1 cases /10,000 inhabitants (min-max, 99.3-160.3) (Figure 3). Compared to Class 2 (the intermediate area with the lowest incidence ratio), the regression analysis (GAM) demonstrated that, adjusting for the geographical location, the peripheral higher economic status class (Class 1) was not significantly different (SIR=0.76 IC95%[0.38; 1.49], p-value=0.43). The hillside poor area with no water availability and open defecation (class 4) was the area most at risk, with a SIR of 3.20 (IC95%[1.25;8.18], p-value=0.015). The coastal flood-prone and poor area (Class 5), with high access to a shallow water table, was significantly at risk with a SIR of 3.15 (IC95%[1.41; 7.04], p-value=0.005).

A C C E P T E D M A N U S C R I P T

A C C E P T E D M A N U S C R I P T 16

The central coastal trading area with good access to a brackish water table and the main use of latrines (class 3), was also at risk with a SIR of 2.49 (IC95%[1.05; 5.95], p-value=0.039).

Discussion

We highlight the heterogeneity of the mechanisms causing cholera persistence in the city of Gonaives, a major transmission hotspot in the Artibonite department and more generally in Haiti. Amongst the neighbourhoods of Gonaives we identified five different classes based on their socio-economic and environmental profiles. The classes exhibited significantly different cholera incidence rates. Three high-risk classes were identified -poor hillside neighbourhoods, poor neighbourhoods of the coastal floodplain and coastal trade neighbourhoods. Compared to the lowest risk Class 2, the highest risk class (Class 4), incorporated the northern hillside poor neighbourhoods, which have limited access to water and frequent open defecation. There is no functional water supply, and water has to be purchased downhill from private pumps or connections to DINEPA network. There are few latrines, and open-air defecation is the most common practice. The other at-risk classes (Class 3 and 5) are the low-lying coastal neighbourhoods, with increased risk of flooding and a high mainly brackish water table. Due to poor maintenance of limited existing water networks, there is no reliable water supply resulting in the use of uncontrolled water supplies, such as local shallow wells. Their contamination is at higher risk due to their close proximity to numerous porous latrines with poor maintenance and sparse dislodging. The low-lying terrain increases flooding risk, and essential water drains are clogged by rubbish. Whereas Class 5 includes poor residential neighbourhoods, important trading activities are present in the Class 3 neighbourhoods of the coastal centre, with large markets where poor hygiene practices, such as dirty hands and contaminated fruit and vegetables, may contribute to disease transmission. The relationship with water is complex: our main results showed that cholera risk is significantly concentrated in an area with high water availability and in another area with low water availability. The problem remains the same, poor access to clean drinking water.

To our knowledge, this is the first study that incorporates a statistical spatial analysis of cholera risk factors with various socio-economic and environmental characteristics, obtained from available geographical databases and field studies, in an urban context. This wide-angle approach enables us to have an overall vision of the patterns of cholera socio-environmental risk factors in the city of Gonaives. However, we note several limitations of our study. The number of suspected cases of cholera, rather than confirmed cases, was used in this study, as only 10% of suspected cases were sampled for confirmation culture, which is in accordance with WHO recommendations (Word Health Organization (WHO) and Global Task Force on Cholera Control, 2010). Secondly, some of the variables included in the analysis are continuous, such as economic status, functionality of the water network, or drinking water availability. However, due to local constraints, the field investigation was only able to observe qualitative information with no systematic sampling of households. Furthermore, classification methods identify variables that exhibit differences between the neighbourhoods. However, there may be variables that play a relevant role in cholera transmission, but we were unable to identify in the analysis due to homogeneity in each neighbourhood. We performed an ecological study, not an individual population study. This introduces classical ecological fallacy, which therefore means we cannot infer a direct individual risk from these neighbourhood risk profiles [START_REF] Richardson | Ecological correlation studies[END_REF].

A C C E P T E D M

A N U S C R I P T 19 However, our results do not aim to describe individual risk factors of cholera transmission, which have already been assessed by several case control studies [START_REF] Dubois | Epidemic cholera in urban Zambia: hand soap and dried fish as protective factors[END_REF][START_REF] Dunkle | Epidemic cholera in a crowded urban environment, Port-au-Prince, Haiti[END_REF][START_REF] Kanungo | Clinical, epidemiological, and spatial characteristics of Vibrio parahaemolyticus diarrhea and cholera in the urban slums of Kolkata, India[END_REF][START_REF] Morof | Community Mortality from Cholera: Urban and Rural Districts in Zimbabwe[END_REF] including a study conducted in Gonaives in 2011 [START_REF] Grandesso | Risk factors for cholera transmission in Haiti during inter-peak periods: insights to improve current control strategies from two case-control studies[END_REF]. Additionally, our study could not capture important factors associated with the dynamics of cholera transmission, such as population movements [START_REF] Bengtsson | Using mobile phone data to predict the spatial spread of cholera[END_REF], population immunity [START_REF] Jackson | Seroepidemiologic Survey of Epidemic Cholera in Haiti to Assess Spectrum of Illness and Risk Factors for Severe Disease[END_REF], rainfall [START_REF] Eisenberg | Examining rainfall and cholera dynamics in Haiti using statistical and dynamic modeling approaches[END_REF] or mitigation activities implemented by cholera response teams. We aimed to focus on neighbourhood characteristics associated with cholera incidence. In our univariate analysis, few risk factors seem to be significantly associated with cholera incidence, probably due to the heterogeneity of cholera-associated socio-environmental characteristics in the city and the limited number of neighbourhoods (28). Furthermore, our methods were chosen to take into account possible highly correlated variables, an important risk of biased estimates during the multivariate analysis. We implemented a dimensional reduction through the use of other statistical methods for large numbers of variables, such as multiple correspondence analysis and hierarchical classification. The identified socio-environmental classes were significantly associated with cholera incidence, which highlights the importance of this approach. Rather than describing the individual risk patterns, our results propose a classification of cholera distribution patterns in one of the most affected communes in the country, in order to implement adapted elimination strategies for the specific socio-economic contexts.

These results have already been incorporated into a contingency plan promoted by several governmental and non-governmental stakeholders; Pan-American Health Organisation (PAHO). This contingency plan is essential to raise funds from international funding agencies, in order to prioritise the extension and/or rehabilitation of water supply networks, as identified in the National Cholera Elimination Plan. Our results helped to identify the neighbourhoods in which to prioritise these middle and long-term prevention strategies (Figure 4). Within the hillside neighbourhoods, the priority is to construct a new water network, which does not currently exist, with enough pressure to ensure clean water supply uphill. Strategies to improve hygiene practices should focus on the washing of hands and the reduction of open defecation practices by promoting the construction of improved latrines. Within the coastal neighbourhoods, improving water quality should focus on the extension and rehabilitation of existing water supply networks. An efficient chlorination system should be implemented. This water supply has to be affordable, in order to reduce the use of domestic wells that are currently easily contaminated due to the high water table and high flooding risk. Furthermore, new latrines that are not in contact with the water table should be constructed in the flood-risk neighbourhoods, alongside the regular dislodging of excrement. To further limit the risk of flooding, it is important to regularly clear water drains.

Education should focus on improving the hygiene practices and conditions in the markets.



  Patterns of cholera socio-environmental risk factors were identified in Gonaives City, Haiti.  Coastal flood plains are high-risk associated with easy access to a contaminated brackish water table.

  survey undertaken by Médecins Sans Frontières, identified Gonaives as one of the main affected towns in 2010-2011, with an estimated crude mortality rate of 19.1 deaths per 1000 person years 95% CI

  with the cholera incidence in Gonaives City, taking into account the geographical coordinates of each commune. In October 2012, the Ministry of Public Health and Population (MSPP) in Haiti launched a National Plan for the elimination of cholera (2013-2022) (Republic of Haiti et al., 2013). It integrates prevention strategies at the short, middle and long term along several themes: water and sanitation, epidemiological surveillance, health promotion and behavioural change, and clinical care of infected patients in healthcare institutions. One of the main aims is to decrease the risk of transmission by rehabilitating existing drinking water supply networks and constructing new waste water collection and treatment systems in the main towns and peri-urban areas. In order to prioritise the areas requiring interventions in Gonaives, it is vital to consider the cholera epidemic in the local infrastructural context. This helps to identify the most vulnerable neighbourhoods at risk of cholera and adapt control and prevention strategies.

  2013 to 22 March 2016 were collected from registry books at the three cholera treatment institutions in Gonaives (Hôpital Secours-Providence in 2013-2015, CDI Raboteau in 2015, then the new Hôpital Secours-Providence in 2015-2016) and surrounding communes of Ennery and L'Estère. The addresses were attributed to the correct neighbourhood with local and experienced staff from ACF non-governmental organization (Action Contre la Faim) and IOM inter-governmental organization (International Organization for Migrations). Neighbourhoods were characterized using the 2008 Local Urbanism Plan of Gonaives provided by the city council; a geographic information system (GIS) analysis; and field studies performed with staff from ACF NGO, who have worked in the city on WASH (water, sanitation and hygiene promotion) and cholera response programs for many years. Based on known risk factors associated with cholera transmission in urban settings, available data and particularities of Gonaives City, the following socio-economic and environmental characteristics were collected for each neighbourhood and are detailed in the supplementary Database: mean and minimum altitude, distance from the sea, flooding risk, water table depth and salinity; population, socio-economic status; drinking water availability (from domestic wells, water pumps or DINEPA water network); prevailing type of sanitation including open defecation practices, unimproved latrines or modern improved toilets; size and cleanliness of markets; presence and traffic activity of motorised taxi stations; main building use and number of floors. Geographical analysis was performed with the use of shapefiles of communes and communal administrative boundaries (obtained from the National geospatial information center, http://www.cnigs.ht, accessed 1 August 2016), a shapefile of Gonaives neighbourhoods (drawn for this study), A C C E P T E D M A N U S C R I P T 9 an altitude raster (magic.csr.utexas.edu/Event/10Haiti/Elevation.cfm, accessed 1 August 2016) and a raster of the population estimated in 2015 (LandScan, http://web.ornl.gov/sci/landscan/, accessed 1 August 2016). Finally, satellite estimates of daily-accumulated rainfall in Gonaives (TRMM_3B42_daily v7, on the commune centroid) were extracted from NASA website (http://gdata1.sci.gsfc.nasa.gov/daac-bin/G3/gui.cgi?instance_id=daily_p_comp, accessed 1 August 2016).

  2013 to 22 March 2016, 2,567 patients with suspected cholera living in the commune of Gonaives were seen in cholera treatment institutions in Gonaives City or in the surrounding communes. During this period, 268 of these patients were sampled for cholera culture confirmation at the National Laboratory of Public Health, and 51% of these stool cultures were positive for Vibrio cholera O1. Most cases were seen in 2013 and 2014, with two main peaks in November 2013 and September 2014, which did not appear to be related to increase rainfall (Figure 1 Panel A). Conversely, 2015 exhibited a relative lull period from April to August, although staff strikes at the Hôpital Secours-Providence and local riots with the robbery of a register book in CDI Raboteau resulted in missing data during this period. Overall, 71% of these 2,567 suspected cases originated from the city of Gonaives, where the cumulated incidence ratio between November 2013 and March 2016 was 96.8 cases /10,000 inhabitants. The monthly number of total suspected cases per communal section is shown in Figure 1 Panel A and B.

Figure

  Figure 1: (A) Epidemic curve of cholera in the commune of Gonaives with

  At the neighbourhood level, variables significantly associated with cholera incidence in univariate non-parametric analysis included mean altitude (p-value=0.02), drinking water availability (p-value=0.01), prevailing open defecation (p-value=0.008), prevailing use of unimproved latrines (p-value=0.023), socio-economic status (p-value=0.04), water table depth (p-value=0.01) and the presence of a functioning DINEPA water network (p-value=0.005).

Figure 2 :

 2 Figure 2: (A) Hierarchical Ascendant Classification (HAC), and (B) Multiple

Figure 3 :

 3 Figure 3: (A) Cholera epidemic curve for each of the five classes of

  authorities, the Ministry of Public Health and Population, the National directorate for potable water and sanitation (DINEPA), UNICEF and the

Figure 4 :

 4 Figure 4: Profile characteristics (blue boxes) of the three high-risk

Acknowledgements

Johane Occilien, Paul Adrien and Jacques Boncy, Ministry of Public Health and Population, Gonaives, Haiti Jude Pierre, National Directorate for Potable Water and Sanitation, Gonaives, Haiti Ytiende Eliezaire, International Organization for Migrations, Gonaives, Haiti

Funding

This work was supported by the United Nations Children's Fund (UNICEF)

In conclusion, factors associated with the cholera spatial distribution are heterogeneous across neighbourhoods of the city of Gonaives, Haiti. Therefore, prevention strategies should not be homogenous. We provide an important example of the need to adapt prevention strategies according to the specific socio-environmental context of each neighbourhood within the city. Not only will these strategies help mitigate and eliminate cholera transmission in Gonaives and the surrounding communes, they will also have an impact on other diseases sharing the same mode of transmission, such as shigellosis, typhoid fever, hepatitis A etc. More generally, this will help improve the livelihood of the population in the third largest city in Haiti. A similar approach is currently being implemented by UNICEF and its partners in other cities affected by cholera in Haiti, such as Cap-Haitien and Saint-Marc. Harnessing the information from spatial analysis, observational field studies and geographical databases provides an adaptive contextualised approach that may be generalised to other urban areas affected by cholera worldwide.

Declaration of interest