Marie Mazoyer
email: marie.mazoyer@meteo.fr

Romain Roehrig

Fanny Duffourg

Olivier Nuissier

| Olivier Nuissier

Simulating a Mediterranean heavy-precipitating event with parametrized convection: role of subgrid-scale topography

HAL is

Getting the right precipitation probability density function in the Mediterranean area is a challenge for most global and regional climate models with parametrized convection. In particular over land, the intensity of heavy precipitating events is often underestimated. In the present study, we provide a process-based analysis of the representation by the CNRM-ALADIN63 regional climate model of one of these events, which occurred in the southeast of France on November 1-2, 2008.

The CNRM-ALADIN63 model, when run in a configuration where the large-scale dynamics is nudged towards that of the ERA-Interim reanalysis, is first shown to capture appropriately the location and intensity of the heavy-1 precipitating event. Then, using a reference convection-permitting simulation of the same event and a conditional sampling approach to identify and characterize convective updraughts, the ability of the model convection parametrization to capture further convective details is assessed. The model misses the occurrence of the updraught mass flux largest values, despite a significant and systematic overestimation of the updraught vertical velocity. The area covered by convective updraught is in fact found severely underestimated, suggesting an inappropriate approach for the convective closure: the event occurred along the foothills of the Massif Central, where the subgrid-scale features of the topography highly interact with the impinging large-scale flow and thereby drive the position and the large area fractions of convective updraughts, at least during the mature phase of the event. A preliminary topographybased convective closure is proposed and implemented in the model to further assess this hypothesis. The results confirm that the parametrization deficiencies can be significantly reduced with a proper inclusion of the subgrid-scale topographic features. [START_REF] Bresson | Idealized numerical simulations of quasi-stationary convective systems over the northwestern mediterranean complex terrain[END_REF]Hagen et al., 2011a) to the less-than-1-km small scales (e.g., [START_REF] Fuhrer | Dynamics of orographically triggered banded convection in sheared moist orographic flows[END_REF].

K E Y W O R D S
The ability of the impinging flow to be lifted by topographic features is often quantified using the Froude number (in particular its moist version). The latter has been shown relevant to predict the location and potential stationary property of the triggered convective systems (i.e. down-or up-stream, e.g., [START_REF] Miglietta | Numerical simulations of sheared conditionally unstable flows over a mountain ridge[END_REF][START_REF] Chen | Orographic effects on a conditionally unstable flow over an idealized threedimensional mesoscale mountain[END_REF].

Until recently, the climate-scale analysis of HPEs in the Mediterranean region, as in many other regions of the world, has mostly relied on global and regional climate models (GCM and RCM, respectively), in which convection is parametrized as a subgrid-scale process (horizontal resolution ranging from 10 to 150 km -e.g., [START_REF] Randall | Beyond deadlock[END_REF][START_REF] Ruti | Med-cordex initiative for mediterranean climate studies[END_REF][START_REF] Mascaro | Evaluation of precipitation from euro-cordex regional climate simulations in a small-scale mediterranean site[END_REF]. Convection parametrization is often pointed as key in many severe weaknesses of these "coarse" resolution models, responsible of large-scale biases in the atmospheric circulation (e.g., [START_REF] Bellucci | The double-itcz syndrome in coupled general circulation models: The role of large-scale vertical circulation regimes[END_REF][START_REF] Dirmeyer | Simulating the diurnal cycle of rainfall in global climate models: resolution versus parameterization[END_REF]Van Weverberg et al., 2018) or in the distribution of precipitation at the daily or hourly timescales [START_REF] Fumière | Extreme rainfall in mediterranean france during the fall: added value of the cnrm-arome convection-permitting regional climate model[END_REF][START_REF] Demott | Convective precipitation variability as a tool for general circulation model analysis[END_REF]. During the last decade, convection-permitting models have emerged as potential tools to go beyond the limits of convection-parametrized models, first at the regional scale (e.g., [START_REF] Kendon | Realism of rainfall in a very high-resolution regional climate model[END_REF][START_REF] Prein | A review on regional convection-permitting climate modeling: Demonstrations, prospects, and challenges[END_REF][START_REF] Coppola | A first-of-its-kind multi-model convection permitting ensemble for investigating convective phenomena over europe and the mediterranean[END_REF], and now at the global scale (e.g., [START_REF] Stevens | The added value of large-eddy and storm-resolving models for simulating clouds and precipitation[END_REF]. Although these convection-permitting models offer new and relevant pathways to understand convective extremes, their variability and their response to various forcings, they remain numerical tools that are particularly computationally intensive. As such, they will not be able to tackle several scientific or society-relevant questions about present and future climate in a foreseeable future, which require particularly long-lasting simulations or large ensembles of simulations (e.g., [START_REF] Deser | Uncertainty in climate change projections: the role of internal variability[END_REF][START_REF] Kay | The community earth system model (cesm) 21 large ensemble project: A community resource for studying climate change in the presence of internal climate variability[END_REF]. Therefore there is still a critical need to understand current convection parametrization weaknesses and find novel ways to overcome them. Recent works exemplify that this is achievable, with notable breakthroughs in the representation of the diurnal cycle of convection over land (e.g., [START_REF] Rio | Shifting the diurnal cycle of parameterized deep convection over land[END_REF][START_REF] Rio | Ongoing Breakthroughs in Convective Parameterization[END_REF][START_REF] Bechtold | Representing Equilibrium and Nonequilibrium Convection in Large-Scale Models[END_REF]. The present work positions along this main objective, bridging together convection-permitting and convection-parametrized simulations to assess the ability of a convection parametrization to capture the main features of a Mediterranean HPE and understand which components of the scheme are likely deficient.

This study focuses on a representative Mediterranean HPE to compare a convection-parametrized simulation based on the regional climate model CNRM-ALADIN63 [START_REF] Nabat | Modulation of radiative aerosols effects by atmospheric circulation over the euro-mediterranean region[END_REF] with a convectionpermitting simulation. The focus is mostly on the representation of convection processes and variables internal to the convection parametrization. The role of subgrid-scale topographic features is investigated, and its relevance for convection parametrization is argued.

The article is organized as follows. Section 2 presents the HPE selected for the present study. Section 3 then details the experimental setup. Section 4 evaluates the representation of the selected HPE in the two simulations and section 5 follows with a more detailed analysis of the representation of convection by the convection-parametrized simulation, using the convection-permitting simulation as a reference. Section 6 further investigates the role of subgrid-scale topographic features in the representation of convection, thereby suggesting a preliminary pathway to improve convection parametrization. Section 7 develops, implements and tests a topography-based convective closure in CNRM-ALADIN63 to further assess how sugrid-scale topographic features may indeed be responsible for some of the parametrization deficiencies. Section 8 finally concludes the present study.

| CASE OVERVIEW

Among the events analysed in [START_REF] Duffourg | Origin of the moisture feeding the heavy precipitating systems over southeastern france[END_REF], we consider the HPE that occurred in the south of France over the Cévennes mountainous range on November 1-2, 2008. [START_REF] Duffourg | Origin of the moisture feeding the heavy precipitating systems over southeastern france[END_REF] provide a global overview of the large-and meso-scale conditions in which the event occurred, and therefore only a brief description is given here. Atmosphériques à la Neige - [START_REF] Durand | A meteorological estimation of relevant parameters for snow models[END_REF][START_REF] Quintana-Segui | Analysis of near-surface atmospheric variables: Validation of the safran analysis over france[END_REF] rainfall analysis (Fig. 1b), the HPE led to daily accumulations of more 200 mm day -1 on November 1, 2008, over several 8x8-km 2 pixels along the Cévennes foothills. Several rain-gauge stations of the Météo-France network recorded more than 400 mm day -1 on that day [START_REF] Duffourg | Origin of the moisture feeding the heavy precipitating systems over southeastern france[END_REF]. This scenario is typical of HPEs occurring over southeastern France [START_REF] Nuissier | A numerical study of three catastrophic precipitating events over southern france. i: Numerical framework and synoptic ingredients[END_REF]. The strong low-level wind impinging the Cévennes likely explains the absence of a large cold pool below the MCS which would have help it to migrate upstream, over the plain area [START_REF] Miglietta | Numerical simulations of sheared conditionally unstable flows over a mountain ridge[END_REF][START_REF] Bresson | Idealized numerical simulations of quasi-stationary convective systems over the northwestern mediterranean complex terrain[END_REF]. This HPE is thus fully relevant for our study as (i) it is representative of a large portion of HPE occurring over the region, and (ii) it does not mix orographic triggering and cold pool mechanisms. Besides, because we are mostly interested in the early stage of the convective system, we mostly focus on its first day, namely on November 1, 2008.

| EXPERIMENTAL DESIGN

3.1 | The regional climate model CNRM-ALADIN63

The convection-parametrized simulation is based on the regional atmospheric model CNRM-ALADIN63 [START_REF] Nabat | Modulation of radiative aerosols effects by atmospheric circulation over the euro-mediterranean region[END_REF]. CNRM-ALADIN63 is used over the Med-CORDEX domain [START_REF] Ruti | Med-cordex initiative for mediterranean climate studies[END_REF], at a 50-km horizontal resolution. The vertical grid encompasses 91 levels from the surface to about 80 km, with a higher resolution in the atmospheric boundary layer. The model atmospheric physics is fully described in [START_REF] Roehrig | The cnrm global atmosphere model arpege-climat 6.3: description and evaluation[END_REF], and we only briefly detail here the convection parametrization.

Convection is parametrized according to the work of [START_REF] Piriou | An approach for convective parameterization with memory: Separating microphysics and transport in grid-scale equations[END_REF] and [START_REF] Guérémy | A continuous buoyancy based convection scheme: one-and three-dimensional validation[END_REF]. This massflux-based scheme considers a single updraught, which aims at representing dry, shallow and deep convection in a continuous way. The convective profile treats dilution processes through organized and turbulent entrainment and detrainment (e.g., [START_REF] Tiedtke | A Comprehensive Mass Flux Scheme for Cumulus Parameterization in Large-Scale Models[END_REF]. Organized entrainment and detrainment processes are based on the buoyancy-sorting approach of [START_REF] Bretherton | Relationships between water vapor path and precipitation over the tropical oceans[END_REF] while turbulent entrainment and detrainment solely depend on the updraught vertical velocity (stronger for weak updraughts, weaker for strong updraughts). A prognostic equation of the updraught vertical velocity (w u ) following the work of [START_REF] Simpson | Models of precipitating cumulus towers[END_REF] is used. It classically includes an aerodynamic drag and a virtual mass parameter in the buoyancy term (see [START_REF] Guérémy | A continuous buoyancy based convection scheme: one-and three-dimensional validation[END_REF]. The scheme closure uses a relaxation of the Convective Available Potential Energy (CAPE) to provide the final value of the updraught mass flux. The latter is then used to scale the fractional area profile α u , given that the updraught vertical velocity is already computed from its own equation. The CAPE relaxation timescale is proportional to the time needed by a buoyant parcel to travel from the base to the top of the convective cloud. Finally, the convective microphysics fully mirrors the large-scale microphysics [START_REF] Lopez | Implementation and validation of a new prognostic large-scale cloud and precipitation scheme for climate and data-assimilation purposes[END_REF]. Namely, it describes the same microphysics variables (cloud liquid water, cloud ice, rain and snow) as the large-scale microphysics scheme, with the same representation of microphysical processes. This consistency allows to consider entrainment of microphysical condensates in the updraught or their detrainment in the environment (for further detail see [START_REF] Piriou | An approach for convective parameterization with memory: Separating microphysics and transport in grid-scale equations[END_REF][START_REF] Roehrig | The cnrm global atmosphere model arpege-climat 6.3: description and evaluation[END_REF].

With respects to surface processes, CNRM-ALADIN63 embarks the SURFEX version 8 platform, as described in [START_REF] Nabat | Direct and semi-direct aerosol radiative effect on the mediterranean climate variability using a coupled regional climate system model[END_REF] (for further detail, see also [START_REF] Roehrig | The cnrm global atmosphere model arpege-climat 6.3: description and evaluation[END_REF][START_REF] Decharme | Recent Changes in the ISBA-CTRIP Land Surface System for Use in the CNRM-CM6 Climate Model and in Global Off-Line Hydrological Applications[END_REF].

| The convection-permitting model Meso-NH

The convection-permitting simulation is based on the non-hydrostatic numerical model Meso-NH [START_REF] Lafore | The meso-nh atmospheric simulation system. part i: Adiabatic formulation and control simulations[END_REF][START_REF] Lac | Overview of the meso-nh model version 5.4 and its applications[END_REF]. It is run at a 2.5-km resolution. Its ability to simulate Mediterranean HPEs, either in realistic or idealized frameworks, has been shown for instance in [START_REF] Nuissier | A numerical study of three catastrophic precipitating events over southern france. i: Numerical framework and synoptic ingredients[END_REF], [START_REF] Duffourg | Origin of the moisture feeding the heavy precipitating systems over southeastern france[END_REF] and [START_REF] Martinet | Fine-scale numerical analysis of the sensitivity of the hymex iop16a heavy precipitating event to the turbulent mixing-length parametrization[END_REF]. The Meso-NH configuration used in the present study is fully described in [START_REF] Duffourg | Origin of the moisture feeding the heavy precipitating systems over southeastern france[END_REF]. The three two-way nested domains can be seen on Fig. 1. The ERA-Interim reanalysis [START_REF] Dee | The ERA-Interim reanalysis: configuration and performance of the data assimilation system[END_REF] is first downscaled at a 40-km resolution over a large domain encompassing the whole Mediterranean Sea, the Eastern Atlantic Ocean, Europe and North Africa (full domain displayed on Fig. 1a). The 10-km-resolution nested domain covers Western Europe and the Western Mediterranean Sea (larger array on Fig. 1a). Finally, the third domain, at a 2.5-km resolution, covers the Northwestern Mediterranean Sea and the southeast of France (smaller array on Fig. 1a, which is also approximately the domain shown on Fig. 2). The two-way nesting approach has been shown to perform often better than the one-nesting approach (e.g., Warner et al., 1997;[START_REF] Harris | An idealized comparison of one-way and two-way grid nesting[END_REF][START_REF] Gao | Improving agcm hurricane structure with two-way nesting[END_REF][START_REF] Madhulatha | Impact of different nesting methods on the simulation of a severe convective event over south korea using the weather research and forecasting model[END_REF]. Though, it may introduce inconsistencies due to different grids, thereby amplifying sometimes noisy adjustments on the coarse domains, especially in the regions along the internal nested boundaries [START_REF] Liu | Simultaneous nested modeling from the synoptic scale to the les scale for wind energy applications[END_REF]. Visual inspection of several instantaneous fields of the Meso-NH simulations did not reveal significant noise near the nested domain boundaries. Besides, the large distance between domain boundaries limits the propagation of such noise up to the Cévennes area which is our main region of interest.

| Simulation set-up

For both models, the initial and lateral conditions are provided by the ERA-Interim reanalysis at its native horizontal resolution (about 75 km, see [START_REF] Dee | The ERA-Interim reanalysis: configuration and performance of the data assimilation system[END_REF].

The Meso-NH simulation starts 4 days before the HPE (October 28, 2008, 0600 UTC), with only the 40-km and 10-km domains activated. The 2.5-km domain is activated the day before the HPE triggered, namely on October 31, 2008 at 1800 UTC (see [START_REF] Duffourg | Origin of the moisture feeding the heavy precipitating systems over southeastern france[END_REF], for further details). The Meso-NH simulation provides instantaneous hourly outputs of various 3D fields.

As for the Meso-NH simulation, the CNRM-ALADIN63 simulation is initialized from the ERA-Interim reanalysis on October 31, 2008 at 1800 UTC for the atmosphere. The surface initial state is provided by a previous long-term simulation with the same model configuration, including spectral nudging of wind vorticity and divergence, surface pressure, temperature and humidity towards the ERA-Interim reanalysis, as described in [START_REF] Nabat | Modulation of radiative aerosols effects by atmospheric circulation over the euro-mediterranean region[END_REF]. This ensures an initial state consistent with the surface model physics, and consistent with the simulated period. For the present simulation, in contrast to [START_REF] Nabat | Modulation of radiative aerosols effects by atmospheric circulation over the euro-mediterranean region[END_REF], spectral nudging [START_REF] Radu | Spectral nudging in a spectral regional climate model[END_REF] towards the ERA-Interim reanalysis is activated only for the wind vorticity and divergence and the surface pressure. This helps to better constrain the synoptic chronology of the simulation, while weakly perturbing the convection scheme. The nudging is maximal above 700 hPa up to the model top, null below 850 hPa, and has a linear transition these two altitudes. Spectral nudging has the advantage to select the scales to be nudged: it is maximal for scales above about 400 km, null below 300 km, and varies linearly between these two spatial scales. Finally, the nudging timescale depends on the variable: 6 hours for wind vorticity, 24 hours for surface pressure, and 48 hours for wind divergence.

| THE HPE AS SIMULATED BY MESO-NH AND CNRM-ALADIN63

| Precipitation

Figs. 2a-c exhibit the accumulation of precipitation from November 1, 2008, 0000 UTC to November 2, 2008 0000 UTC for the reference SAFRAN product, and for the CNRM-ALADIN63 and Meso-NH simulations. SAFRAN and Meso-NH data are conservatively interpolated onto the CNRM-ALADIN63 horizontal grid. SAFRAN precipitation maxima (∼ 108 mm day -1 ) are located on the southern slopes of the Massif Central mountainous range, and surrounded by weaker but significant precipitation amount upwind along the mountain foothills and downwind over the main mountainous area. Both simulations broadly capture the precipitation pattern, especially the two pixels of maximum precipitation. These two maxima have similar amplitude in the two simulations. Note that Meso-NH also captures the precipitation pattern of the HPE at 2.5 km, though with a probable slight underestimate of the highest rain rates (242 at 2.5 km compared to 245 mm day -1 in SAFRAN at 8 km -not shown). CNRM-ALADIN63 precipitation spreads upwind and downwind over a larger area (i.e. more precipitation over the Massif Central foothills and the main mountainous area), while the Meso-NH simulation slightly underestimates precipitation along the Massif Central foothills. Overall the two simulations exhibit high performance in capturing the HPE mean spatial footprint.

| Properties of the regional atmospheric flow

Several metrics are now introduced to quantify the properties of the regional atmospheric flow on November 1, 2008, at 1200 UTC, just before the HPE initiation, with an emphasis on the convective and orographic instabilities of the flow (Fig. 2d-l). In particular, the properties of the flow impinging the mountainous range are key for determining the location and intensity of the HPE (e.g., [START_REF] Bresson | Idealized numerical simulations of quasi-stationary convective systems over the northwestern mediterranean complex terrain[END_REF]. For instance, according to [START_REF] Chu | Effects of orography on the generation and propagation of mesoscale convective systems in a two-dimensional conditionally unstable flow[END_REF] and [START_REF] Chen | Orographic effects on a conditionally unstable flow over an idealized threedimensional mesoscale mountain[END_REF], the moist Froude number (Fr w ) partially controls the location and propagation of a convective system in the presence of mountainous area. A high Fr w favours the mechanical triggering of convection over the mountain range as the flow can overpass it, while a low Fr w favours blocking conditions, dynamical convergence along the mountain foothills and thus convection there. The Froude number is defined in appendix A, following [START_REF] Chen | Orographic effects on a conditionally unstable flow over an idealized threedimensional mesoscale mountain[END_REF]. In particular, a height scale of 1200 m, relevant for the Massif Central is used (e.g., [START_REF] Bresson | Idealized numerical simulations of quasi-stationary convective systems over the northwestern mediterranean complex terrain[END_REF]. Of course, the Froude number is not sufficient to determine the location of any HPE (e.g., [START_REF] Chu | Effects of orography on the generation and propagation of mesoscale convective systems in a two-dimensional conditionally unstable flow[END_REF]: a large CAPE may enhance upstream propagation if the orographic forcing is sufficiently large to lift the incoming air parcels above the level of free convection (LFC). The Fr w , CAPE and LFC are shown on Fig. 2 for ERA-Interim, CNRM-ALADIN63 and Meso-NH. The three datasets indicate that the impinging flow has a rather high Fr w , greater than one (Fig. 2d-f), indicating that mechanical lifting can occur over the Massif Central. Together with rather the low altitude of the LFC (Figs. 2j-l), at least in CNRM-ALADIN63 and Meso-NH, this is consistent with precipitation maxima located along the first slopes of the mountainous range. Convective instability is present, though weak, with a CAPE of only a few hundreds of J kg -1 (Figs. 2g-i). The CAPE is slightly higher in CNRM-ALADIN63 than in Meso-NH, which may indicate that the CNRM-ALADIN63 convective scheme has not been active enough in the morning. Overall, CNRM-ALADIN63 and Meso-NH are rather consistent, while ERA-Interim exhibits particularly low CAPE and elevated LFC. Such a difference remains so far unclear and the radiosoundings in the region on that particular day do not have enough data to compute and validate these parameters. Nevertheless, the qualitative consistency between CNRM-ALADIN63 and Meso-NH allows us to compare these two simulations in more details.

| CNRM-ALADIN63 REPRESENTATION OF CONVECTION

This section qualitatively evaluates the internal variables of the CNRM-ALADIN63 convection scheme, using the convection-permitting simulation as a reference, and thereby identify whether CNRM-ALADIN63 captures the main HPE features for the right reasons or due to some hidden compensating errors. Convective updraught are identified and characterized in the convection-permitting simulation, based on a threshold on the vertical velocity (e.g., [START_REF] Lemone | Cumulonimbus vertical velocity events in gate. part i: Diameter, intensity and mass flux[END_REF][START_REF] Khairoutdinov | Large-eddy simulation of maritime deep tropical convection[END_REF]. There is no consensus on the threshold value to be used, but 1 m s -1 seems appropriate to characterize the core updraughts of convective systems that explain most of the non-local vertical transport of mass, energy and moisture. The sensitivity of our results to this threshold value was explored and even though they quantitatively differ, the comparison between CNRM-ALADIN63 and Meso-NH remains qualitatively similar, and this does not impact our main conclusions. Following [START_REF] Derbyshire | Sensitivity of moist convection to environmental humidity[END_REF], we consider only cloudy updraught (positive liquid or ice water content), as it helps us remove most of orographic gravity waves. Finally, to compare the CNRM-ALADIN63 parametrized updraught with those explicitly simulated by Meso-NH, the detection criteria is applied on the Meso-NH simulation at its native resolution (2.5 km) and the properties (e.g., mass flux, vertical velocity, area fraction) are aggregated over each cell of the CNRM-ALADIN63 horizontal grid. Note that the CNRM-ALADIN63 updraughts are assumed to be fully sub-grid, i.e. the model resolved vertical velocity does not include any direct contribution from them. The size of the convective systems compared to the model horizontal resolution may however question this convection parametrization paradigm, but this remains beyond the scope of this work.

The model behaviour on the grid cell having the maximum accumulated precipitation on November 1, 2008 (44.25 • N, 3.66 • E) is illustrated on Fig. 3. The CNRM-ALADIN63 convection parametrization significantly underestimates the convective mass flux (Fig. 3a-b), especially from 1200 UTC. The convective parametrization also simulates only shallow convection in the morning, while in Meso-NH, deep yet weak convection occurs. Nevertheless, this shallow convection provides significant precipitation, more than in Meso-NH (not shown). The updraught mass flux decomposition in CNRM-ALADIN63 also exhibits a significantly different balance between updraught vertical velocity and area fraction (Figs. 3c-f). On the one hand, the CNRM-ALADIN63 updraught area fraction remains mostly below 6%, while reaching much larger values in Meso-NH, up to 38% near 4 km above sea level. On the other hand, the CNRM-ALADIN63 updraught vertical velocity is strongly overestimated, reaching up to 10 m s -1 in the mid-troposphere, while, in contrast, Meso-NH does not simulates mean updraught vertical velocities higher than 3 m s -1 (note that the 2.5-km Meso-NH resolution remains rather coarse to fully capture the detail of convective updraughts, thereby inducing some uncertainty in our estimates of their vertical velocity). Thus the CNRM-ALADIN63 convection scheme compensates the area fraction small values with stronger updraught vertical velocity, albeit still resulting in an underestimated convective mass flux. The overestimated updraught vertical velocity is consistent with an overestimated buoyancy within the updraught (Figs. 3g-h). The latter behavior may be related to the moister updraught environment in the CNRM-ALADIN63 simulation than in the Meso-NH simulation (not shown), but also to too weak entrainment rates at the edges of the updraught.

Beyond this single grid cell, the model behavior is further documented on Fig. 4 for all the grid cells within a larger domain covering the Cévennes area (indicated by the thick black line on Fig. 2b) and for the whole day of November 1, 2008. It is synthesized using the bi-variate probability density function between the altitude and a given updraught properties (mass flux, vertical velocity, area fraction, buoyancy).

Consistently with Fig. 3a-b, CNRM-ALADIN63 does not capture the strongest values of the mass flux.

CNRM-ALADIN63 also simulates a bi-modal distribution with a main mode below 0.02 kg m -2 s -1 and another weaker one near 0.15 kg m 2 s -1 . In contrast, Meso-NH simulates a more continuous distribution from the smallest to the highest values. Stronger mass flux values are also found in the lower troposphere, under 2 km, indicating enhanced entrainment of environmental air near the base of convective cells. As in Fig. 3e, the updraught vertical velocity in CNRM-ALADIN63 is most of time strong, with a maximum frequency of occurrence of 8 m s -1 between 2 and 4 km (Fig. 4i), while it mostly remains below 3-4 m s -1 in Meso-NH. Once again, this compensates in CNRM-ALADIN63 the small updraught area fractions (Fig. 4d-f).

The bi-modal shape of the mass flux distribution is explained by that of the area fraction distribution, which is not seen in Meso-NH.

In Summary, CNRM-ALADIN63 exhibits two main errors that compensate each other: on the one hand an overestimation of the updraught vertical velocity and, on the other hand, an underestimation of the updraught area fraction. In the following, we focus on the area fraction (thus the closure function of the parametrization) and assess to what extent the sub-grid properties of the topography may contribute to its strong underestimate by the CNRM-ALADIN63 convection parametrization. The updraught vertical velocity errors may be related to errors in the convective profile, especially with respect to too weak entrainment rates.

This will be addressed in a future work.

| THE ROLE OF SUBGRID-SCALE TOPOGRAPHY

This section now aims at assessing the relationship between the small-scales (i.e. subgrid-scale) features of the regional topography and the area covered by updraughts. Our reference topograpy is that of the reference simulation, namely the 2.5-km-resolution topography of Meso-NH.

Fig. 2 provides a broad view of the topography at the scale of the south-east of France, while Figs. 5b-c document the slope and orientation of the topography (see appendix B for their computation). Slopes are mostly of the order of 4 to 8% over the Cévennes and Massif Central, with a few places above 12%, close to the mountain foothills. There, their orientation mostly falls between 50 • and 150 • , meaning that the slopes point to the north and the northwest. The mechanical forcing for southerly and southeasterly air masses is thus high, as in the present case. Fig. 5a quantifies the location where convection is particularly active on November 1, 2008, based on the accumulated Meso-NH mass flux on that day. Convectively-active regions concentrate over small areas, qualitatively anchored near the crests at the top of steep slopes, approximately aligned with the large-scale low-level wind. Thus, together with the conditional instability of the flow, the dynamical lifting by the small-scale topographic features significantly drives the location and coverage of the most convectively-active areas (see also e.g., [START_REF] Smith | The role of orography in the regeneration of convection: A case study from the convective and orographically-induced precipitation study[END_REF]. Thermally-driven circulation on steep slopes may also be of importance (e.g., Hagen et al., 2011b). However, the simulated buoyancy flux over the Massif Central is mostly negative on that day (not shown), thus rather inhibiting convection than driving it. This is likely expected for such a case as warm and moist air from the Mediterranean Sea flows over a mountainous surface that cooled during the previous night.

CNRM-ALADIN63 resolves part of the topography and thus part of its impact on the flow. We therefore identify the subgrid-scale features of the topography (i.e. those which scale is below 50 km) to assess to what extent they can drive the subgrid-scale convective processes in the Meso-NH simulation, in particular the updraught area, and can be used in the CNRM-ALADIN63 convection parametrization. Fig. 6 quantifies the subgrid-scale heterogeneities of the topography: the Meso-NH topography is first smoothed at 50 km using a simple 21-point running average operator. Then, for any property of the topography, we compute its anomaly against the same property computed from the smoothed topography. The slope and slope orientation anomalies are illustrated on Figs. 6a,b, respectively. The comparison between Fig. 5a and Fig. 6 indicates that locations favourable to convection initiation and maintenance are qualitatively related to subgrid-scale features of the topography, namely small-scale slopes that are significantly steeper than the resolved topography and well aligned with the large-scale flow. This is not captured by the CNRM-ALADIN63 resolution. Fig. 6 further shows that the main CNRM-ALADIN63 grid cells receiving the highest precipitation amount simulated in either CNRM-ALADIN63 or Meso-NH (see Fig. 2b-c) are also those having large subgrid-scale variability of the topography altitude and slopes. Such large subgrid-scale variability probably increases the number of convective updraughts in the CNRM-ALADIN63 grid cells, which, in turn, may explain the larger updraught area fractions indicated by Meso-NH but missed by CNRM-ALADIN63.

To further assess this relationship and provide possible insight in the parametrization of the effect of the topography subgrid-scale variability onto convection (also subgrid), we adopt the parametrization paradigm.

The latter assumes a clear separation between the resolved and parametrized processes. In our case, the interactions between the resolved flow and the resolved topography are thus supposed to be handled by the CNRM-ALADIN63 dynamical core, while the interaction between the resolved flow (its dynamics and thermodynamics) and the parametrized convection only involves the topography anomalies with respect to the resolved topography. Fig. 7 exhibits the properties of the topography anomalies computed with respect to the topography smoothed at the CNRM-ALADIN63 50-km resolution. Namely, Fig. 7 shows the distributions of the anomalous topography slope and slope orientation with respect to the large-scale wind direction at 2400 m above sea level, the wind components being horizontally averaged over the relevant CNRM-ALADIN63 grid cell, considering only the Meso-NH grid points falling in the eight CNRM-ALADIN63 grid cells indicated with the thick black lines on Figs. 5 and 6. The same distributions but conditionally-sampled to the occurrence of an updraught at 2400 m during at least one hour of the analysed day are superimposed with a red line. The conditionally-sampled distributions are accumulated over both space and time and consist of 13 725 updraught occurrences (a pixel can thus be counted several times during the day). This allows to account for any change in wind direction along the analysed day. The updraught occurrence is significantly favoured over steep slopes (Fig. 6a.), oriented within about ±30 • of the large-scale wind directionat an angle near -30-30 • , thus consistently aligned with the south-southeasterly large-scale wind. Of course, small-scale topography does not explain all the details of the updraught occurrence, as heterogeneities in the flow instability, the convective organisation, etc, are also likely key players.

Finally, in the perspective of accounting for subgrid-scale topographic features in a convective parametrization, we quantify the fraction α oro of a CNRM-ALADIN63 grid cell, which is likely contributing the updraught triggering or maintenance. Any Meso-NH grid cell meeting the following criteria is supposed to contribute to this area fraction: slope of the anomalous topography s ≥ 4.0% and orientation of the anomalous topography within ±30 • of the 50-km-scale wind direction at 2400 m above sea level. This fraction is computed for each of the eight 50-km grid cells indicated on Fig. 5 and 6 with thick black squares, and for each hour of November 1, 2008. This represents 192 estimates of α oro that can be compared to the updraught area fraction α u computed from the Meso-NH simulation (Fig. 8). The comparison is thus made only on the main grid cells where convection occurs and where subgrid-scale topographic features are likely to play a role, as these are obviously not the unique ingredient controlling convection properties. Fig. 8 confirms that highest values of α u also correspond to the highest values of α oro . The opposite is however not true, especially during the early stage of the convective event. This likely indicates that other processes are acting during its triggering phase (e.g., small scale heterogeneities in the flow instability or in the surface buoyancy flux). To summarize, the convectively-active topographic fraction α oro contributes to the large updraught area fractions observed in Meso-NH, and missed by CNRM-ALADIN63, especially during the mature phase of the event. This clearly indicates that subgrid-scale topographic features must be accounted for in the convection parametrization closure.

| TOWARDS A TOPOGRAPHY-BASED CONVECTIVE CLOSURE IN CNRM-ALADIN63

To further assess how the subgrid-scale topography may be at the origin of the CNRM-ALADIN63 convection parametrization deficiencies in this HPE context, a preliminary implementation of this effect is performed.

The Meso-NH 2.5-km resolution topography is used to compute several descriptors of the subgrid-scale topography within the CNRM-ALADIN63 grid cells covering the southeast of France. Namely, considering the topography anomalies with respect to the 50-km smoothed topography, we compute the fraction of the CNRM-ALADIN63 grid cells which have slopes greater than 4%. This fraction is then distributed over 12 bins of direction, each bin being 30 • wide. These 12 topography parameters are read at the beginning of the simulation. Then, at each model timestep, the fraction of any CNRM-ALADIN63 grid cell, orographically active for convection, is computed as the cell fraction αoro having a slope greater than 4% and a slope direction within ±30 • of the model low-level wind direction (computed from the wind components averaged over the first 200 m above the surface). This is consistent with the fraction α oro defined in the previous section. Given the updraught area fraction α CAPE provided by the CAPE closure, the updraught area fraction α u is updated as follows:

α u = [1 -f (U BL ) ] α CAPE + f (U BL )k oro αoro . ( 1
)
U BL is the mean low-level wind magnitude computed from the wind components averaged over the first 200 m above the surface. f is a smooth function of the wind speed which is equal to 0 below 1 m s -1 and to 1 above 5 m s -1 (a squared sinus in between). k oro is a tuning parameter. The depth over which the wind is averaged probably needs to be further tuned. This tuning of the parametrization is not addressed in the present work.

With this parametrization of the subgrid topography effect on the convection closure, two new simulations are run with CNRM-ALADIN63, testing two values of k oro : 2 and 4. A value of 1 was also tested, but, as expected from Fig. 8, it does not allow to reach area fractions comparable to those estimated in the Meso-NH simulation, and it is therefore not discussed hereafter. The precipitation averaged over the event is weakly impacted compared to the standard CNRM-ALADIN63 configuration (not shown). The maximum precipitation rate is also slightly increased by about 10 mm day -1 from k oro = 0 to k oro = 2 (similarly from k oro = 2 to k oro = 4 -not shown). Fig. 9 illustrates the impact of the new closure formulation on the convection parametrization internal variables on the same grid cell as Fig. 3. The updraught area fraction is significantly increased to 8-9% for k oro = 2 and to 12-13% for k oro = 4. In both cases, the updraught mass flux is increased, mainly due to the updraught area fraction increase. For k oro = 2, the increased updraught vertical velocity also contributes. In contrast, the updraught vertical velocity decreases in the case of k oro = 4, while convection becomes shallower and more intermittent. The enhanced mass flux likely makes convection more efficient to stabilise the troposphere, as seen with the updraught buoyancy. Except that the parametrized updraughts are shallower during the core of the HPE, the configuration with k oro = 4

simulates updraughts more comparable with those estimated with Meso-NH (Fig. 3). This is statistically confirmed for the CNRM-ALADIN63 grid cells covering the Cévennes area (Fig. 10).

In summary, this preliminary implementation of a topography-based convective closure emphasises the role of the subgrid-scale features of the topography in driving some of the CNRM-ALADIN63 deficiencies in the representation of convection during the HPE: the parameterised updraught area fractions, mass fluxes, and vertical velocities are larger, stronger and weaker, respectively, when subgrid-scale topography effect are accounted for, making them more consistent with our reference convection-permitting simulation.

However, the updraughts become shallower during the core of the event, highlighting that the subgrid-scale orography does not explain all the parametrization deficiencies. On the one hand the present topographybased convective closure remains preliminary and requires further development and consolidation. On the other hand, while a new tuning of the parametrization is likely required, the entrainment/detrainment parametrization component probably needs further evaluation in this HPE context.

| CONCLUSIONS

The present study evaluates the representation of convection by the regional climate model CNRM-ALADIN63 during a heavy precipitating event, which occurred on November 1-2, 2008 over the mountainous area of the Cévennes in the southeast of France. For that purpose, a reference simulation of the same event by the convection-permitting model Meso-NH is used to diagnose several properties of convection. As documented by [START_REF] Nuissier | A statistical downscaling to identify the large-scale circulation patterns associated with heavy precipitation events over southern france[END_REF] and [START_REF] Duffourg | Origin of the moisture feeding the heavy precipitating systems over southeastern france[END_REF], the studied event resulted from a conditionally-unstable air flow of Mediterranean origin impinging the foothills of the Massif Central, in the Cévennes area. For November 1, 2008, the SAFRAN analysis indicates a daily precipitation maximum of about 246 mm, located on the southern slopes of the Massif Central mountainous range, and surrounded by weaker but still significant precipitation along the upwind side of the mountain range. The two simulations appropriately capture the HPE spatial footprint and intensity at the scale of the CNRM-ALADIN63 horizontal resolution. They are also consistent in terms of impinging flow properties: low-CAPE instability, level of free convection between 500 and 1000 m above sea level, moist Froude number greater than unity indicating that the flow is able to go over most of the Massif Central rather than being constrained to go around. The thermodynamical and dynamical (orographic) instabilities of the flow (CAPE, LFC and Froude number) are slightly higher in the CNRM-ALADIN63 simulation than in the Meso-NH simulation. The similarities between the two simulations emphasize the ability of 50-km-resolution regional climate models with parametrized convection to capture the coarse properties of a specific heavy precipitating event in the southeast of France, at least when the large-scale circulation is constrained through model nudging.

The evaluation of the CNRM-ALADIN63 convection scheme is pursued using a conditional sampling approach based on vertical velocity and cloud liquid or ice water occurrence to identify convective and cloudy updraughts in the Meso-NH simulation and characterize their bulk properties at the scale of the CNRM-ALADIN63 grid cells (e.g., [START_REF] Lemone | Cumulonimbus vertical velocity events in gate. part i: Diameter, intensity and mass flux[END_REF][START_REF] Khairoutdinov | Large-eddy simulation of maritime deep tropical convection[END_REF][START_REF] Derbyshire | Sensitivity of moist convection to environmental humidity[END_REF]. The range of updraught mass flux values is similar in the two simulations, but the CNRM-ALADIN63 simulation misses its largest values, even though it overestimate the updraught vertical velocity. In fact, the area fraction covered by convective updraught is severely underestimated, and this cannot be fully compensated by the vertical velocity overestimation.

As a result, the role of subgrid-scale topographic features in driving part of the large updraught area fraction simulated by Meso-NH is further investigated. Unsurprisingly, it is first found that, in the context of this particular event, the strong slopes of the topography, well aligned with the large-scale impinging flow concentrates most of the convective vertical transport over November 1, 2008. For a 50-km resolution, most of theses slopes are subgrid-scale. They are then shown to explain a large proportion of the highest updraught area fractions simulated by Meso-NH.

To further assess the ability of these subgrid-scale features of the topography to dynamically impact and potentially improved the CNRM-ALADIN63 parameterised convection, a preliminary topography-based convective closure is developed, implemented and tested in the model. Partly inspired by [START_REF] Leung | A subgrid parameterization of orographic precipitation[END_REF], a bin approach combining slope and slope direction is taken, as to reduce the computational burden of keeping detailed information about the high-resolution topography within a climate model. When this new closure activates in the case of strong horizontal wind, it yield updraughts covering larger fractions of the model grid cell than those predicted by the CAPE closure. Though their vertical velocity is weaker due to more active convection during the first stages of the event which stabilises the atmospheric column, the updraught mass flux is overall increased. As a result, the three convective variables become more consistent with those estimated from the Meso-NH reference simulation. The main (new) caveat is shallower convective towers, either calling for a new tuning of the convection parametrization or for further revisit of its other components (e.g., entrainment/detrainment).

The new convective closure is preliminary but encouraging. Further development is needed to assess it in a wider set of conditions. In particular the consideration of the flow instability with respect to orographic lifting is probably key to quantify the subgrid-scale topography slopes that are truly active for triggering and maintaining convective cells. The wider availability of regional and global convection-permitting simulations should help refine this topography-based convective closure and make it more robust. In particular, other interactions between the topography and convection, such as the thermally-driven circulations along slopes, are known to be of paramount importance.

In the present work, our reference convection-permitting simulation has only a 2.5-km resolution, which likely remains a bit coarse to fully capture the detailed structure of the convective updraughts and the transition from boundary-layer thermals to deep convective cells. As a consequence, the initiation of deep convection might be slightly oversensitive to orographic features in our reference simulation. Simulations at higher resolution of a similar heavy-precipitating event would help better disentangle the role of the topography from other boundary-layer processes in triggering deep convection.

AC K N O W L E D G E M E N T S
We thanks the three anonymous reviewers for their helpful and constructive comments.

A | MOIST FROUDE NUMBER

Following [START_REF] Chen | Orographic effects on a conditionally unstable flow over an idealized threedimensional mesoscale mountain[END_REF] and [START_REF] Bresson | Idealized numerical simulations of quasi-stationary convective systems over the northwestern mediterranean complex terrain[END_REF], the non-dimensional moist Froude number Fr w reads:

Fr w = ∫ hm 0 U (z )d z h m ∫ hm 0 N v (z )d z (2)
where U is the wind speed, N v the Brunt-Väisälä frequency for stable moist air (taken to 0 in case of unstable air) computed using the virtual air temperature and h m is a reference height taken to 1 200 m to be representative of the Massif Central mountain range. When the topography is higher than this reference height, the Froude number is simply not calculated.

B | SLOPE AND SLOPE ORIENTATION

Considering a topography described by the height h (x, y ) above the Earth geoide, its slope s (x, y ) is computed as:

s (x, y ) = 100 × ∂h (x, y ) ∂x 2 + ∂h (x, y ) ∂ y 2 (3)
The multiplication by 100 allows us to express the slope in percent (%). The slope orientation γ (x, y ) reads: (d-f), (g-i) and (j-l) same as (a-c) but for the updraught area fraction, vertical velocity (in m s -1 ) and buoyancy (in m s -2 ), respectively. See text for details on the computation of these convective variables in the Meso-NH simulations.

γ (x, y ) = tan -1 ∂h (x, y ) ∂ y - ∂h (x, y ) ∂x (

Fig. 1a exhibits

  Fig. 1a exhibits the large-scale situation over Western Europe on November 1, 2008, at 0600 UTC. The 500-hPa geopotential height indicates a large upper-level trough, which extends southwards over the Iberian Peninsula. Close to the surface, the pressure low centred over Spain drives a strong south-westerly flow over the Gulf of Lion, which becomes almost southerly when entering over the land of the southeast of France. Such a synoptic configuration is referred to as a Cyclonic South Westerly pattern in the climatological study of Nuissier et al. (2011). The southerly moist flow that impinges the Cévennes is weakly conditionally unstable (see Figs. 2g-i) and leads to the triggering of a quasi-stationary MCS on the Cévennes foothills at about 1000 UTC on November 1, 2008. The MCS lasts until the day after and stops at approximately 0600 UTC on November 2, 2008. As shown by the SAFRAN (Système d'Analyse Fournissant des Renseignements

  This research has received funding from the French National Research Agency (ANR) project REMEMBER (grant ANR-12-SENV-001). It was also supported by the DEPHY2 project, funded by the French national program LEFE/INSU and the GdR DEPHY. The authors would like to thank S. Somot (CNRM) for useful discussions and comments along this work.

  4)Numerically, the atan2 operator is used to get first an angle between -180 • and 180 • . This angle is then converted in the range 0 • -360 • to follow a conventional wind orientation. The direction increases from 0 • to 360 • clockwise, 0 • and 90 • pointing to the northward and eastward directions, respectively. F I G U R E 1 (a) ERA-Interim 500-hPa geopotential height (grey contours, in m), mean sea-level pressure (colors, in hPa) and 925-hPa wind (vectors, in m s -1 ) on November 1, 2008, at 0600 UTC. (b) SAFRAN precipitation accumulated from November 1, 2008, 0000 UTC to November 2, 2008, 0000 UTC. Black contours indicate the topography (one contour every 400 m starting at 400 m. The black arrays indicate the intermediate 10-km and innermost 2.5-km Meso-NH domains. U R E 2 Top row: Precipitation (in mm) regridded on the CNRM-ALADIN63 grid accumulated from November 1, 2008, 0000 UTC to November 2, 2008, 0000 UTC, for the SAFRAN analysis (left), CNRM-ALADIN63 (center) and Meso-NH (right). Second row: Moist Froude number and 925-hPa wind on November 1, 2008, 0000 UTC for ERA-Interim (left), CNRM-ALADIN63 (center) and Meso-NH (right). Third and fourth row: same as second row but for the CAPE (in J kg -1 ) and the LFC (in m), respectively. Grey shading indicates missing values, i.e. grid cell where the computed diagnostics are not meaningful. Black contours indicate the topography (one contour every 400 m starting at 400 m). The black thick rectangle on panel b indicates the Cévennes domain used for Fig. 4.

  U R E 3 Temporal evolution on November 1, 2008 of the updraught mass flux (a-b, in kg m -2 s -1 ), area fraction α u (c-d), vertical velocity (e-f, in m s -1 ) and buoyancy (g-h, in m s -2 ) for the CNRM-ALADIN63 grid cell with the maximum daily accumulated precipitation on November 1, 2008 (44.25 • N, 3.66 • E). The left and right columns correspond to the CNRM-ALADIN63 and Meso-NH simulations, respectively. The altitude is provided in km above sea level. See text for details on the computation of these convective variables in the Meso-NH simulations. U R E 4 Bivariate probability density function (PDF) of the altitude (in km above sea level) and the updraught mass flux (in kg m -2 s -1 ) for (a) CNRM-ALADIN63 and (b) Meso-NH based on all CNRM-ALADIN63 grid cells within the domain indicated with thick black lines on Fig. 2b and every hour on November 1, 2008. (c) emphasizes the PDF of the updraught mass flux at 2400 m for the two simulations.

  U R E 5 (a) Updraught mass flux at 2400 m above sea level simulated by Meso-NH and accumulated from November 1, 000 UTC to November 2, 2008, 000UTC (in kg m -2 ). See text for the updraught identification in the Meso-NH simulation. (b) Slope in % and (c) slope orientation in degrees north of the Meso-NH 2.5-km-resolution topography. A slope orientation of 0 • indicates that the height of topography increases from south to north (thus northward slope). 90 • corresponds to an eastward slope. The 925-hPa horizontal wind simulated by Meso-NH at 1200 UTC on November 1, 2008 is superimposed with arrows. The eight squares indicated by the thick black lines indicate the CNRM-ALADIN63 grid cells used for further diagnostics (Figs. 7 and 8).

  U R E 6 (a) Slope (%) and (b) slope orientation (same convention as as on Fig.5c) of the topography anomalies computed against the 50-km smoothed Meso-NH topography, as computed using a 21-point running average operator. The eight squares indicated by the thick black lines indicate the CNRM-ALADIN63 grid cells used for further diagnostics (Figs.7 and 8).

  density function (PDF) of (a) the slope and (b) slope orientation of the topography anomalies with respect to the 50-km smoothed Meso-NH topography. The slope orientation is adjusted with respect to the 50-km-scale wind direction at 2400 m above sea level. 0 • indicates that the slope and the mean low-level wind are aligned. Only the Meso-NH grid cells within the eight squares indicated by thick black lines on Fig.6are considered to compute the raw PDF (black line). The red line corresponds to the PDF when those Meso-NH grid points are conditionally sampled to the occurrence of an updraught at 2400 m above sea level (see text for details). The grey shading indicates the 95% confidence interval that can be used to compare the raw and conditionally-sampled distributions. It is estimated using a bootstrap approach. U R E 8 Scatter plot between α oro and α u based on the eight CNRM-ALADIN63 grid cells indicated by thick black lines on Fig.6and for each hour of November 1, 2008 (192 points). α u is the Meso-NH updraught area fraction and α oro is the convectively-active fraction of subgrid-scale topographic features as defined by the criteria of a slope s ≥ 4% and a slope orientation within ±30 • of the mean low-level wind direction. Colors emphasise the time UTC on November 1, 2008 (in hour) when the area fractions are computed. U R E Same as Fig.3but for the CNRM-ALADIN63 simulations in which the topography-based convective closure is activated (see Section 7 for details). Left column: k oro = 2; right column k oro = 4 U R E 1 0 Same as Fig.4but for the CNRM-ALADIN63 simulations in which the topography-based convective closure is actived (see Section 7 for details). (a), (c), (e) and (g): k oro = 2: (b), (d), (f) and (h):k oro = 4

W. I., Qian, Y., Berg, L. K., Liu, Y., Huang, M., Ahlgrimm, M., Forbes, R., Bazile, E., Roehrig, R., Cole, J., Merryfield, W., Lee, W. S., Cheruy, F., Mellul, L., Wang, Y. C., Johnson, K. and Thieman, M. M. (2018) CAUSES: Attribution of Surface Radiation Biases in NWP and Climate Models near the U.S. Southern Great Plains. Journal of Geophysical Research: Atmospheres, 123, 3612-3644. Warner, T. T., Peterson, R. A. and Treadon, R. E. (1997) A tutorial on lateral boundary conditions as a basic and potentially serious limitation to regional numerical weather prediction. Bulletin of the American Meteorological Society, 78, 2599-2618. URL: https://journals.ametsoc.org/view/journals/bams/78/11/1520-0477_ 1997_078_2599_atolbc_2_0_co_2.xml.

Funding information

French National Research Agency (ANR),

French national program LEFE/INSU and


Arun Dhaka

Olivier Jeannin

Emmanuel Aubert

Enrique Espinosa

Marc Fourmigué

Ie-Rang Jeon

Ie- Rang Jeon

N-Chlorobenzimidazoles as efficient and structurally diverse amphoteric halogen bond donors in crystal engineering

published or not. The documents may come   L'archive ouverte pluridisciplinaire

The ability of amphoteric N-chlorobenzimidazoles to self-associate into 1D chains through strong and linear N-Cl•••N halogen bond interactions is demonstrated. Less polarisable Cl atom is strongly activated thanks to the intramolecular amphoteric character and the intermolecular cooperativity effect. The obtained family of compounds featuring different substitution patterns opens perspectives toward the elaboration of macroscopic polar structures.

During the past decades, the development of molecular solids with new applications has been achieved through supramolecular approaches that rely on our understanding of intermolecular interactions 1,2 i.e. hydrogen bond. 3 In this context, the halogen bond (XB) 4 has been established as an efficient structuring tool in crystal engineering strategies owing to its advanced fundamental understanding, inherent directionality and most significantly predictability. 5,6 Strong activation of iodine in classical XB donors like diiodotetrafluorobenzne 7 and diiodoperfluoroaalkanes 8 allows them to competently participate in co-crystal formation with an electron rich partner. However, difficult activation of a sigma-hole in XB donors involving rather less polarizable halogen atoms find limited use of lighter elements i.e. F < Cl << Br in supramolecular synthesis. On the other hand, molecular assemblies recognized through selfassociation ability, require a self-complementary system bearing both an activated donor and an acceptor site. 9,10 Such amphoteric molecules for XB remain limited, which might be attributed to the difficult integration of donor and acceptor sites in a molecule at a desired position. For instance, linear X-cyanoacetylene were demonstrated to organize into infinite 1D chains through short X•••N contacts with distances 2.970 , 2.956, and 2.930 Å respectively for X = Cl, Br, and I (Scheme 1a). [11][12][13] Crystal structure of a neutral diaryl-5-iodo-1,2,3-triazole was shown to exhibit an antiparallel arrangement of molecular tapes formed through I•••N contacts. 14 Moreover, examples of self-associating XB donors resulting into a dimeric [15][16][17] and trimeric [18][19][20] structures in the solidstate are also known. A self-complementary system that features structural diversity is of importance as it allows us to induce the structural polarity in the solid-state. For example, moving from X-CC-CN to perfluorinated X-C6F4-CN based system promoted a reduced interaction between the XB bonded 1D chains that favoured a co-parallel/non-centrosymmetric arrangement over an antiparallel/centrosymmetric arrangement. 21 Interestingly, X-C6F4-CN (X= I, Br) resulted into 1D chain motif through X•••N contacts whereas X-C6F4-CN (X= F, Cl) crystallized differently favouring CN••• interactions (Scheme 1b). This example indeed elegantly demonstrates the tendency of lighter halogen atoms to involve into interactions that are less structure directing. Therefore, use of XB based self-associating systems in crystal engineering strategies, especially those involving lighter halogen atoms (X = F, Cl) remains very challenging. This can be attributed to the more demanding activation of a sigma-hole on highly electronegative atoms as well as the synthetic complexity or chemical instability of such systems. Benzimidazole derivatives bearing both hydrogen bond donor and acceptor sites were recently reported to systematically selfassociate into hydrogen bonded polar chains (Scheme 1c). 22 Inspired from this work, we considered replacing the N-H with N-Cl functionality anticipating that activation of a strong sigma-hole on chlorine atom would similarly favour the molecular assembly into halogen bonded polar chains (Scheme 1d). In this manuscript, we report synthesis of structurally diverse N-chlorobenzimidazole derivatives that upon crystallization self-assemble into the sought after XB bonded polar chains though linear and short N-Cl•••N contacts. Calculation of electrostatic potential (ESP) in a series of Nchlorobenzimidazole derivatives (Fig. S3) revealed presence of a strong electron depletion on the Cl atom and an electron concentration on the nucleophilic N atom (see Table 1). This indeed confirm a favourable self-association ability of our system. Scheme 2 Synthesis strategy with atom numbering (a) and Nchlorobenzimidazole derivatives synthesized in this study (b). Following a reported procedure, 23 functionalization by chlorine was achieved in a single step by treating commercially available benzimidazoles with sodium hypochlorite (NaOCl) solution (Scheme 2) to afford a white crystalline powder of

N-chloro-2- methylbenzimidazole (2•CH3), N-chloro-2-bromobenzimidazole (2•Br), N-chloro-2-mchlorobenzimidazole (2•Cl), and N-chloro-5,6- dichloro-2-methylbenzimidazole ((5,6-Cl)2•CH 3
) in moderate yields (30-60%). Slow evaporation of ethyl acetate solution containing these newly synthesized derivatives yielded colourless crystals. To our surprise, the only known crystal structure for Nchlorobenzimidazole is 1-chlorobenzimidazole itself (R = R' = H), albeit incorrect. 19 The structure was reported indeed with an extra H atom on the sp2 N3 atom of the heterocycle.

Compound 2•CH3 crystallizes in the monoclinic system, a centrosymmetric space group P21/c with each molecule situated in general position (Fig. S1a). Moreover, we observed an antiparallel arrangement of XB chains that are primarily driven by strong electrophilic•••nucleophilic interactions with a Cl•••N distance of 2.648 Å, corresponding to a reduction ratio (RR) of 0.80 relative to the sum of their van der Waals radii (Fig. 1). The XB directionality in this motif is manifested by a highly linear Cl•••N contact with ND-Cl•••NA = 176.3° (Table 1). These structural details indeed highlight the excellent ability of systematic self-association driven by Cl•••N XB interactions. Apart from that, we also observe that these XB chains running in the opposite direction are interconnected through secondary intermolecular interactions such as C-H•••N and ••• interactions (Fig. S1b). These benzimidazole-based XB systems offer a large structural diversity that allowed us to think of having different substituents at the imidazole unit (C2 position) as well as on the phenyl ring (C5 and C6 positions). To this end, we synthesized 2•Cl and 2•Br derivatives, which crystallize in the monoclinic and the orthorhombic systems, respectively, with space groups P2 1 /c and Pcab. Their crystal structures revealed a similar centrosymmetric arrangement of XB bonded chains driven through equally short and linear Cl•••N contacts (Fig. 2, Table 1).

The Cl atom at C2 position in 2•Cl does not participate in any interaction the molecular chains communicate with each other mostly through ••• interactions (Fig. 3a). Interestingly, the more polarizable Br atom in 2•Br participates in directional Br••• Please do not adjust margins Please do not adjust margins interactions between the neighbouring chains allowing them to organize orthogonal to each other (Fig. 3b). Such recurrent X••• interactions (X= Br, I) have been well established in literature that possibly occur when an electrophilic halogen is looking for a partner in the absence of a nucleophilic atom (N, O, S etc). 24,25 Overall, these results further prove the robust nature of Cl•••N contacts observed in the different derivatives of N-chlorobenzimidazole that strongly promote their systematic selfassociation.

This motivated us to now investigate the effect of possible substituents present on the phenyl ring. Therefore, we synthesized (5,6-Cl)2•CH 3 that has two chlorine atoms positioned in a symmetrical manner at C5 and C6 of the phenyl ring. As shown in Figure 4, the structure reveals a similar antiparallel arrangement of XB molecular chains formed through Cl•••N contacts of comparable strength (RR = 0.84). In this case, the molecular chains communicate with each other through various modes of secondary interactions i.e. weak type-I Cl•••Cl contacts, Cl•••H and ••• interactions (Fig. S2). It is evident that structural diversity in N-chlorobenzimidazolebased system allows us to introduce different interactions between these 1D chains to tune their relative orientation with each other without rupturing the self-association ability. Moreover, the persistent occurrence of Cl•••N halogen bonds in these derivatives suggests strong activation of a sigma-hole on Cl atom that effectively participate in molecular recognition. With these four structures in hand, we took the opportunity to look at the correlation between activation in donor/acceptor and shortening of the Cl•••N distance (RR). In this series, ESP calculations for isolated molecules in gas phase reveal the weakest activation of Cl (VS,max, Table 1) in 2•CH3 and the strongest in (5,6-Cl)2•CH 3 whereas that of nucleophilic nitrogen (V S,min , Table 1) follow exactly the opposite trend. This is a consequence of electron withdrawing power of different substituents present at C2 or C5/C6 position, i.e. a decrease in the nucleophilic power of NA/ND might enhance the electrophilic character of Cl. Moreover, comparable values of RR relate well to similar values of the electrostatic potential difference V = VS,max -VS,min, which integrates the permanent electrostatics at both donor and acceptor sites in monomers (Table 1).

To have a further insight on strong activation of ND-chlorine atom in this family, we compared the electrostatic potential values of halogen atoms on other positions as well. As can be seen from Table S2 and Figure S3, the sigma-hole of chlorine at ND position is strongly activated compared to the C2-chlorine (+87% increase), whereas the chlorine atoms at positions C5 and C6 are strongly deactivated. Moreover, C2-bromine atom is more polarized than C2-chlorine atom (+31% on VS,max), which explains the involvement of the C2bromine in Br••• interactions in contrast to the C2-chlorine revealing no particular interaction. It is worth to note that the C2-Cl/Br atoms still have a significantly lower electrostatic potential than the ND-Cl atom which is particularly well activated. The difference in activation of chlorine depending on its position is also clearly visible on its integrated topological atomic charge: whereas ND-Cl is positively charged (+0. Such strong chlorine activation at ND position is indeed due to the amphoteric character of the imidazole ring. Calculations performed on two model compounds, where NA in 2•Cl is replaced by either the CH-or CH2-group, show a significant decrease of the VS,max value at Cl at ND position along the series (+29,42, +25,51 and +16,99 kcal/mol for N, CH and CH2 groups, respectively). Furthermore, in addition to the amphoteric intramolecular effect, the intermolecular cooperativity effect in single chains assembled by ND-Cl•••NA interactions appears to be also significant. Increasing the n number of molecules in the chain affords a monotonic increase in the VS,max value, ranging from +29,22 (n = 1) to +38.15 kcal/mol (n = 6) (Table S3). In contrast, such ND-X•••NA (X=Cl, Br, I) interactions were shown to exhibit an additive effect on a self-assembly of XB dimers and trimers. 26 Lastly, to check the effect of interchain ••• interactions, we compared two anti-parallel chains assembled by N D -Cl•••N A and ••• interactions (formed around inversion centres in the crystal structure). Consequently, this interaction induces the decrease of VS,max values (for instance, from +29.22 to +27.84 kcal/mol for n = 1) while maintaining the cooperative effect with the increase of chain length (from +27.84 (n = 1) to +32.94 kcal/mol (n = 5)) (Table S4).

In conclusion, we have presented here synthesis and singlecrystal structures of a new class of neutral self-complementary XB systems based on N-chlorobenzimidazole. In this series, strong activation of a sigma-hole on rather less polarizable Cl-atom primarily drives the molecular assemblies into halogen bonded 1D chains through short and directional Cl•••N contacts. With this respect, it should be pointed the synergetic intramolecular effects of electron withdrawing groups and the amphoteric character of the molecule, to which is added the intermolecular effect of cooperativity. Moreover, with the excellent self-association ability, this system also offers a vast structural diversity that allows introducing different substituents to tune the communication between the chains. These results open perspectives towards supramolecular synthesis of halogen bonded solid-state structures with a macroscopic polar arrangement.

This research was supported by the French National Research Agency Grant ANR-17-CE07-0025. A. D, O. J., M. F, and I.-J. thank the CDIFX in Rennes for the use of X-ray diffractometer. The EXPLOR mesocenter is thanked for providing access to their computing facility (project 2021CPMXX2483).

Scheme 1

 1 Scheme 1 Literature examples of self-associating XB 1D systems (1a & 1b). Hydrogen bonding in benzimidazoles leading to 1D polar chains (1c), and halogen bonding in N-chlorobenzimidazoles illustrated in this work (1d).

Fig. 1 Table 1

 11 Fig. 1 Solid-state arrangement of 2•CH3 with antiparallel arrangement of XB chains. The XB motif (RR = 0.80) is shown in green dotted lines. Hydrogens on the phenyl ring have been omitted for clarity.Table 1 ESP values and solid-state structural details of Nchlorobenzimidazoles. Units: VS,max, VS,min and V = VS,max-VS,min in kcal/mol, angle (ND-Cl•••NA) in degrees, and RR(NA•••Cl) is unit less. Compound VS,max (Cl) VS,min (NA) V RR angle 2•CH3 2•Br 2•Cl (5,6-Cl)2•CH3

Fig. 2

 2 Fig. 2 Molecular arrangement with antiparallel arrangement of XB chains in 2•Cl (a) 2•Br (b). The XB motifs (RR = 0.81 and 0.80) are shown in green dotted lines. Hydrogens on the phenyl rings have been omitted for clarity.

Fig. 3

 3 Fig. 3 Secondary intermolecular interactions in crystal structure of 2•Cl exhibiting ••• interactions (a) and 2•Br exhibiting Br••• interactions (b). Hydrogens on the phenyl rings have been omitted for clarity.

Fig. 4

 4 Fig. 4 Molecular arrangement with antiparallel arrangement of dipolar XB chains in (5,6-Cl)2•CH3. The XB motif (RR = 0.84) is shown in green dotted lines. Hydrogens on the phenyl rings have been omitted for clarity.

  14 in (5,6-Cl)2•CH3; +0.16 in 2•Cl), chlorine at positions C2 and C5,6 are negatively charged (-0.14 in 2•Cl; -0.21 in (5,6-Cl)2•CH3).

Conflicts of interest

There are no conflicts to declare

Notes and references


J E Bergez
email: jbergez@toulouse.inra.fr

Using DEXi functions

This document explains how to use the DEXi_fct R code in order to work with DEXi-tree models under R (DEXi_fct) The code is in constant evolution. For more information,

A Node is a R object made of -id="numeric", #unique sequential id of the node -name="character", #name of the node -isLeaf="logical", #is it a leaf? -isLeafAndAggregated="logical", #is this leaf also an aggregated node ! -children="character", #list of the names of the node's children -sisters="character", #list of the names of the node's sisters -mother="character", #name of the node's mother -aggregation="matrix", #if aggregated node, table of aggregation -Proba="numeric", #Estimated weight of aggregation. If Leaf set basically to uniform -Depth="numeric", #Depth of the node -Twin="numeric", #In case of multiple leaves, give the id of the other leaves -CondiProbaList="list", -rangeScale="numeric", #range scale -scaleLabel="character", #Labels of the different scales -nodePath="character" #Node path from root to leaf

To access an attribute Y of the object X under R use "@" : X@Y e.g. aNode@name give you the name of the node aNode

Basic functions

To create a DEXi tree under R modelName<-"MASC2.0.dxi" #The name of the DEXi file version <-"DEXi_fct_1.9.r" #The name of the current file of function #Load the functions to read and manipulate the DEXi model source(paste("functions",version,sep="/")) #Depends of the structure of your directory ## Creating a generic function for 'print' from package 'base' in the global environment #read the structure of the model MT<-xmlTreeParse(paste("data",modelName,sep="/"),useInternalNodes=T)

listDEXi<-createTree(MT) DEXi <-listDEXi[[1]]
The function createTree returns a list in case of multiple trees (for example when integrating the satellite trees in the model). Then one choose the tree he wants to work with To view the tree: three possibilities: print, show, describe -print(aTree): gives a short summary of the tree

## < 1 > Z: Contribution au developpement durable ## < 2 > -Y : Dimension economique ## < 3 > --Y : Resultats economiques ## < 4 > ---X : Rentabilite ## < 5 > ---Y : Autonomie economique ## < 6 > ----X : Independance economique ## < 7 > ----X : Efficience economique ## < 8 > ---X : Surcout en materiel ## < 9 > --Y :
Capacite productive a long terme ## < 10 > ---Y : Maitrise de la fertilite physico-chimique ## < 11 > ----X : Maitrise du statut acido-basique du sol ## < 12 > ----X : Maitrise de l etat structural du sol ## < 13 > ----X : Maitrise de la fertilite phosphopotassique ## < 14 > ---Y : Maitrise des bioagresseurs ## < 15 > ----X : Maitrise des maladies et ravageurs

To manipulate the tree

You can create subtree: createSubTree <-function(aTree,nodeName) -aTree: the name of the main tree -nodeName: the name of the node to cut the tree You can change the aggregation table: createAggregationMatrix <function(aNode,expectedWeight,nbTables=1,popSize=50,iters=50) -aNode: the node that will be modified -expectedWeight: an array of weights -nbTables: the number of matrix to be returned -popSize and iters: parameters for the algogen

myWeights <-c(0.2,0.2,0.6) newAggregation <-createAggregationMatrix(DEXi@Nodes[[1]],myWeights,5) newTree <-DEXi newTree@Nodes[[1]]@aggregation <-newAggregation[[1]] getEstimatedWeights(newTree@Nodes[[1]]) ## xDimension economique xDimension sociale ## 0.2 0.2 ## xDimension environnementale ## 0.6
You can change the uniformity of the leaf modalities: changeLeavesUniformity<function(aTree,aDistributionList)

DEXi@Nodes[[4]]@Proba ## [1] 0.25 0.25 0.25 0.25 DEXi@Nodes[[6]]@Proba ## [1] 0.25 0.25 0.25 0.25 myDistributionList <-list() myDistributionList[[1]] <-c(0.5,0.2,0.2,0.1) ; names(myDistributionList)[1] <-"Rentabilite" myDistributionList[[2]] <-c(0.1,0.2,0.3,0.4) ; names(myDistributionList)[2] <-"Independance economique" newTree<-changeLeavesUniformity(DEXi,myDistributionList)

To simulate

You can create options: createOptions <-function(aTree,nbOptions=1,aSeed=-1) -aTree: the name of the Tree -nbOptions: the number of option to be created -aSeed: in case one wants to create the same series of option unless random values are choses option <-createOptions(DEXi,nbOptions=3) head(option)

## [,1] [,2] [,3] ## Rentabilite 2 3 4 ## Independance economique 3 4 2 ## Efficience economique 2 1 3 ## Surcout en materiel 2 2 1 ## Maitrise
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] ## Contribution au developpement durable 2 5 3 4 1 1 2 ## Dimension economique 4 4 4 3 1 1 3 ## Resultats economiques 2 3 4 3 1 1 3 ## Rentabilite 3 3 4 4 1 1 2 ## Autonomie economique 1 3 3 2 2 2 3 ## Independance economique 2 2 2 3 1 1 4 ## [,8] [,9] [,10] ## Contribution au developpement durable 2 3 1 ## Dimension economique 2 5 1 ## Resultats economiques 1 3 1 ## Rentabilite 2 3 1 ## Autonomie economique 1 2 2 ## Independance economique 1 3 1
You can save the scenarios: saveScenarios <-function(scenarios,file="scenarios.csv")

The structure of the created file is:

"" "V1" "V2" "V3" "V4" "V5" "V6" "V7" "V8" "V9" "V10" "Contribution au developpement durable" 1 1 2 1 1 2 3 4 2 2 "Dimension economique" 2 2 3 2 2 2 3 5 2 3 "Resultats economiques" 3 4 2 3 2 1 3 4 2 2 "Rentabilite" 3 4 2 3 3 1 3 4 1 1 "Autonomie economique" 2 3 3 3 1 2 3 3 4 4 "Independance economique" 1 4 2 4 2 1 2 3 4 4

You can view a scenario (barplot): showScenario<function(aScenario,aTree,isLabelY=TRUE,isPar=T) -aScenario: scenario to be shown -aTree: the tree -isLabelY: a boolean to draw the Y-axis -isPAR: a boolean to format the graph

showScenario(as.matrix(scenarios[,1]),DEXi)

You can compare some scenario (radar plot): compareScenario <function(aTree,theScenarios,listNodes) -aTree: the tree -theScenarios: the scenarios to be compared -listNodes: the variables to be represented compareScenario(DEXi,scenarios,c("Contribution au developpement durable","Dimension economique","Dimension sociale","Dimension environnementale"))

To perform sensitivity analysis of the tree

You can run an ANOVA -First check the requested time: infoAOV <-function(aTree,iTest=50) -Then if "adequate" run the anova AOV_DEXi<-function(aTree): two models are tested: additif and power 2 -And then show results: showAOV <-function(aAOV_DEXi, main.show=T, nb.plot=8, beside=T,las=1,...)

ShowMC(DEXi@Nodes[[2]],MC)

## 1 2 3 4 5 ## 0.17 0.35 0.22 0.17 0.09

ShowMC(DEXi@Nodes[[19]],MC)

## 1 2 3 ## 0.39 0.28 0.33 #Write a file that contain all the results obtained analysing the random options selected for the MC analysis write.table(MC,file=paste("output/","MC results_.csv",sep=""),sep=",",row.names=T,col.names=NA)

You can run a Sensitivity index computation -Perform the analysis and save the file (to be used with Excel macro):SI_DEXi <function(aTree,fileName="SI_out.csv",isFile=T) -And then show the results for a given aggregated node showSI<-function(aTree,aSI) SI <-SI_DEXi(DEXi,paste("output","SI_out.csv",sep="/"),TRUE)

  name: Contribution au developpement durable ## Number of attributes: 65 ## Number of aggregated attributes: 26 ## Number of true leaves (no multiple, no aggregated): 39 ## Maximum depth: 6 ## List of repeated aggregated nodes: Non ## No multiple leaves ## No Leaf-Aggregated Leaf -show(aTree): represents the hierarchy of the tree and give information on repeated leaves. In case of multiple leaves or leaf-aggregated nodes ID of similar nodes are provided show(DEXi)

  root to node: Contribution au developpement durable -> ## Is it a leaf: FALSE ## Is is a leaf-aggregated: FALSE ## Mother: Root ## Sisters: None ## Children: Dimension economique Dimension sociale Dimension environnementale ## Estimated weights: 33.33 33.33 33.33

  subTree <-createSubTree(DEXi,"Dimension economique") print(subTree) ## Root name: Dimension economique ## Number of attributes: 20 ## Number of aggregated attributes: 8 ## Number of true leaves (no multiple, no aggregated): 12 ## Maximum depth: 4 ## List of repeated aggregated nodes: Non ## No multiple leaves ## No Leaf-Aggregated Leaf You can get the estimated weights of an aggregation table: getEstimatedWeights <function(aNode)

  You can perform a One-At-a-Time analysis -Perform the analysis: OAT<-function(aTree,option) -And then show the results: showOAT<-function(nodeName,aResults,aTree) optionOAT<-createOptions(DEXi,1) results<-OAT(DEXi, optionOAT) showOAT("Contribution au developpement durable",results,DEXi)