Sébastien Blein
email: sébastienblein@meteo.fr

Romain Roehrig

Aurore Voldoire

Parametrizing the mesoscale enhancement of oceanic surface turbulent fluxes: A physical-statistical approach

The meso-scale enhancement of surface turbulent fluxes at the air-sea interface is driven by the meso-scale surface wind speed variability, especially the gustiness velocity and the meso-scale wind speed standard variation. This study proposes a parameterization of these two variables. A large dataset based on the operational 2.5 km-AROME convection permitting model is used in a coarse-graining framework, to quantify various quantities that are subgrid at the scale of a 100-km resolution global circulation model grid cell. This provides a learning dataset to help build the parameterization. The analysis of two case studies of intense wind speed meso-scale variability, combined with the literature review, provides a physically-based set of twelve potential predictors, accounting for the convection activity and the large-scale dynamics. The least absolute shrinkage and selection operator then frames a penalized multivariate linear regression approach to objectively identify the most relevant predictors. Five predictors are selected for predicting the gustiness velocity: the updraft mass flux at the lifting condensation level, the density current spreading velocity, the large-scale horizontal shear and divergence and the large-scale wind speed. The parameterization of the meso-scale wind speed standard deviation requires an additional predictor, namely the cold-pool object aggregation index. The proposed parameterization performs significantly better than the previously published parameterizations and is able to capture 80%, 99% and 93% of the meso-scale enhancement of the momentum, sensible heat and latent heat fluxes, respectively. In the perspective of a global circulation model implementation, in which some predictors may be unavailable, simpler versions of the parameterization, i.e. involving fewer predictors, are also discussed.

Introduction

The exchanges of energy, water and chemical compounds at the ocean-atmosphere interface are critical for the Earth system and its variability (e.g., [START_REF] Trenberth | Atmospheric circulation climate changes[END_REF][START_REF] Timmermann | El niño-southern oscillation complexity[END_REF]1 Zhang 2005;[START_REF] Gulev | North atlantic ocean control on surface heat flux on multidecadal timescales[END_REF]. In a large part, they are driven by the turbulence intensity within the atmosphere surface layer (e.g., [START_REF] Geernaert | Air-Sea Exchange: Physics, Chemestry and Dynamics[END_REF]. The associated surface turbulent fluxes (e.g., the surface momentum flux τ , the sensible heat flux H and the latent heat flux LE) provide boundary conditions to either the atmosphere or the ocean, especially within the context of Earth system modelling. In atmospheric models, these fluxes have to be parameterized as they are governed by motions of smaller scales than the model resolution. Most parameterizations (e.g., [START_REF] Zeng | Intercomparison of bulk aerodynamic algorithms for the computation of sea surface fluxes using TOGA COARE and TAO data[END_REF], for a review) are based on the Monin-Obukhov similarity theory (MOST - [START_REF] Monin | Basic laws of turbulent mixing in the surface layer of the atmosphere[END_REF], which is applied throughout the so-called bulk formulas.

The MOST describes the turbulence within the atmospheric surface layer, namely the layer where the turbulent flux vertical variations are negligible. Its depth is typically of the order of a few meters to a few tens of meters. The surface layer turbulence is driven locally by the wind shear and the stratification. The MOST also assumes horizontal homogeneity of the surface and atmospheric state over the targeted area. It therefore quantifies the locally-generated turbulence at a scale smaller than those possibly associated with meso-scale motions such as boundary-layer structures generated by convective processes (typically from a few hundreds of meters to a few tens of kilometres). However, turbulence generated by processes outside from the surface layer often significantly contributes to surface fluxes at large scale (typically between 5 and 100 km, [START_REF] Sun | Scale dependence of air-sea fluxes over the western equatorial pacific[END_REF] and systematically acts as a flux enhancement (e.g. [START_REF] Redelsperger | A parameterization of mesoscale enhancement of surface fluxes for large-scale models[END_REF]. This contribution, which will be referred here to as the meso-scale flux enhancement, remains sub-grid in current Global Circulation Models (GCM -resolution of the order of 100 km), and thus needs to be parameterized [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF].

Above ocean, the research community has mainly focused on convective activity as a firstorder contributor to the meso-scale flux enhancement. [START_REF] Redelsperger | A parameterization of mesoscale enhancement of surface fluxes for large-scale models[END_REF] propose a scale separation, between the scale of boundary-layer free convection (typically a few hundreds of meters) and that of deep convection (typically a few tens of kilometres). [START_REF] Beljaars | The parametrization of surface fluxes in large-scale models under free convection[END_REF] suggests that other types of meso-scale motion not resolved by GCMs contribute to surface fluxes. [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF] confirm that purely dynamical meso-scale processes, such as convergence lines or orographic perturbations, significantly contribute to GCM-scale flux.

Since the 1990's, several parameterizations of the meso-scale flux enhancement have been proposed, based on a gustiness approach, which add a wind correction in the bulk formulas [START_REF] Godfrey | On the turbulent fluxes of buoyancy, heat and moisture at the air-sea interface at low wind speeds[END_REF][START_REF] Miller | The sensitivity of the ECMWF model to the parameterization of evaporation from the tropical oceans[END_REF][START_REF] Beljaars | The parametrization of surface fluxes in large-scale models under free convection[END_REF][START_REF] Fairall | Bulk parameterization of air-sea fluxes for tropical ocean-global atmosphere coupled-ocean atmosphere response experiment[END_REF][START_REF] Jabouille | Modification of surface fluxes by atmospheric convection in the TOGA COARE region[END_REF][START_REF] Mondon | A study of a fair weather boundary layer in TOGA-COARE: Parameterization of surface fluxes in large scale and regional models for light wind conditions[END_REF][START_REF] Redelsperger | A parameterization of mesoscale enhancement of surface fluxes for large-scale models[END_REF][START_REF] Williams | A physically based parametrization for surface flux enhancement by gustiness effects in dry and precipitating convection[END_REF][START_REF] Zeng | Parameterization of wind gustiness for the computation of ocean surface fluxes at different spatial scales[END_REF][START_REF] Hourdin | LMDZ6a: The atmospheric component of the IPSL climate model with improved and better tuned physics[END_REF]. For example, [START_REF] Jabouille | Modification of surface fluxes by atmospheric convection in the TOGA COARE region[END_REF] suggested to relate the gustiness velocity U g to the precipitation rate, as a proxy of the convective downdraft and cold pool impacts on the surface layer variability. [START_REF] Emanuel | Development and evaluation of a convection scheme for use in climate models[END_REF], [START_REF] Williams | A physically based parametrization for surface flux enhancement by gustiness effects in dry and precipitating convection[END_REF], [START_REF] Hourdin | LMDZ6a: The atmospheric component of the IPSL climate model with improved and better tuned physics[END_REF] followed a similar but more direct approach, without any proxy. [START_REF] Redelsperger | A parameterization of mesoscale enhancement of surface fluxes for large-scale models[END_REF] studied three different parameterizations of U g , based on the rainfall rate, the updraft mass flux and the downdraft mass flux, respectively. [START_REF] Zeng | Parameterization of wind gustiness for the computation of ocean surface fluxes at different spatial scales[END_REF] argued for the relevance of the cloud fraction. More recently, [START_REF] Bessac | Stochastic parameterization of subgrid-scale velocity enhancement of sea surface fluxes[END_REF] and [START_REF] Bessac | Scale-aware spacetime stochastic parameterization of subgrid-scale velocity enhancement of sea surface fluxes[END_REF] chose both the GCM-resolved wind speed and the rainfall rate as two complementary U g predictors. The authors also addressed the stochastic component of the meso-scale flux enhancement. In the context of parameterizing sea salt or dust emissions, a few other studies propose a parametric statistical approach to model the nearsurface wind speed meso-scale variability. The wind speed distribution is parameterized according to a given family of distributions (e.g., Weibull) and model variables are used to derive its parameters. Then, based on a stochastic sampling approach, it is used to compute the desired emission flux (e.g. [START_REF] Ridley | Toward resolution-independent dust emissions in global models: Impacts on the seasonal and spatial distribution of dust[END_REF][START_REF] Zhang | Quantifying the impact of sub-grid surface wind variability on sea salt and dust emissions in cam5[END_REF] For the deterministic component of the meso-scale flux enhancement, which is the main focus of the present study, the literature review unveils at least three main potential caveats in current parameterizations:

• the gustiness approach is assumed valid, even though [START_REF] Beljaars | The parametrization of surface fluxes in large-scale models under free convection[END_REF] warned that this intuitive and pragmatic approach might not be optimal;

• convective activity is considered as the only driver of meso-scale flux enhancement;

• predictors used to parameterize U g are chosen in a rather heuristic and univariate manner, though based on intuitive physically-based conceptual models. This may not be optimal.

The recent availability of kilometer-scale long atmospheric simulations, which explicitly represent meso-scale circulations occurring at the scale of a few kilometers, provided a unique and large dataset to [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF] to address the first two limits, at least for the scales resolved by such simulations. On the one hand, the authors showed that the wind correction between the true wind speed (scalar average) and the GCM-resolved wind speed (vector average) in the bulk formula is able to account for most of the meso-scale enhancement of the sensible and latent heat fluxes. In that sense, the gustiness velocity U g approach, which makes use of the quadratic wind difference is valid. For the momentum flux, another contribution associated to the wind speed meso-scale variance must be accounted for (10% contribution on the momentum flux enhancement on average above all the studied region during one month). On the other hand, [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF] illustrated that other processes are responsible for meso-scale heterogeneity which significantly contributes to the GCM-scale flux, such as orographic or coastal features or large-scale circulations.

The Blein et al.'s dataset now provides the opportunity to more systematically quantify the contribution of many GCM-scale parameters to the GCM-scale fluxes, in a multivariate framework. In particular, while physical intuition yields a wide number of potential predictors for the meso-scale flux enhancement, statistical methods help to more objectively select and combine the most relevant and less redundant ones. This trade-off between purely physical approaches (focusing in a limited amount of physical processes) and purely statistical approaches (such as machine-learning based on raw model parameters) seeks to benefit from the physical processes knowledge while accounting for a potential high complexity degree, which is enhanced, for instance, by the processes co-occurrence and interaction. Such an approach hopefully paves the way to parameterizations of a wider use, while offering further insight in the underlying processes.

This article is organised as follow. Section 2 formalizes the meso-scale flux enhancement background based on [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF], especially introducing useful notations for the whole present work. Section 3 presents the coarse-grained numerical framework at the basis of this study and the statistical method used for the predictor selection and to design a skillful parameterization. The reference dataset main statistical properties are documented in Section 4. Section 5 details the proposed parameterization set up and its testing. Its performance for representing the meso-scale flux enhancement is also presented in this Section. Section 6 intercompares the present parameterization with previous ones. With the objective of preparing GCM implementation, simpler versions of the present parameterization, involving fewer predictors, are also discussed in Section 7. Finally, Section 8 concludes the present work.

Surface fluxes and meso-scale enhancement a. Surface fluxes

The MOST bulk formulas reads:

F φ (∆U, ∆φ) = -A φ C φ (∆U, ∆φ) ∆φ∆U (1)
where: F φ is either the momentum flux norm ∥τ ∥, the sensible heat flux H or the latent heat flux LE. A φ stands for ρ a , ρ a c pa or ρ a L v , where ρ a is the near-surface air density, c pa the moist air specific heat and L v the latent heat of vaporization. C φ is either C D , C H and C E , namely the transfer coefficients for momentum (drag coefficient), heat and humidity, respectively. ∆U is the wind speed difference between the near-surface atmospheric value (the first model level) and the surface value. ∆φ is the difference between the near-surface atmospheric value and the surface value of the potential temperature θ, the specific humidity q or the horizontal wind speed U . As sea surface currents are not considered in this study, ∆U will be replaced by U in the remainder of this paper. The transfer coefficients C φ are generally defined from empirical stability functions (e.g., [START_REF] Businger | Flux-profile relationships in the atmospheric surface layer[END_REF]) and roughness lengths. The latter depends on the sea state, itself depending on the turbulence within the surface layer [START_REF] Fairall | Bulk parameterization of air-sea fluxes: Updates and verification for the COARE algorithm[END_REF]. Iterative methods are therefore often used to estimate the transfer coefficients [START_REF] Liu | Bulk parameterization of air-sea exchanges of heat and water vapor including the molecular constraints at the interface[END_REF][START_REF] Zeng | Intercomparison of bulk aerodynamic algorithms for the computation of sea surface fluxes using TOGA COARE and TAO data[END_REF][START_REF] Fairall | Bulk parameterization of air-sea fluxes: Updates and verification for the COARE algorithm[END_REF][START_REF] Edson | On the exchange of momentum over the open ocean[END_REF]). If the sea state is not available, known or considered, transfer coefficients only depend on U and ∆φ.

b. Meso-scale surface flux enhancement Equation 1 can be used to compute surface fluxes over a wide area such as a GCM grid cell, as long as the horizontal homogeneity hypothesis is fulfilled. If meso-scale motions exist at subgrid scale, the MOST should not be applied directly:

F φ ̸ = F φ = F φ U , ∆φ (2)
where F φ is the true average of the surface flux over the grid cell and F φ is the surface flux computed from GCM-resolved variables. For the latter, U = U = ∥(u, v)∥ is the GCM-resolved (vector average) wind speed and ∆φ represents either U , ∆θ = ∆θ or ∆q = ∆q. The meso-scale flux enhancement is thus defined as: [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF] show that the true GCM-scale fluxes F φ can be well approximated by:

F φMS = F φ -F φ (3)
∥τ ∥ ≃ F U U + δU + ρ a C D U + δU σ 2 U (4) H ≃ F θ U + δU, ∆θ (5)
LE ≃ F q U + δU, ∆q (6)
with δU = U -U = ∥U∥ -U . For the momentum flux, a second correction term is needed to account for the contribution of the meso-scale wind speed variance σ 2 U = U ′2 [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF],prime indicates local deviation to the grid-cell average). When agregated over one month and over a large Indian Ocean domain, approximations 4, 5 and 6 lead to average flux bias reductions from -5.5% to -0.7%, -3.8% to -0.2% and -2.2% to 0% for momentum, sensible heat and latent heat, respectively. When focusing only on cases of intense meso-scale flux enhancement (greater than 10% of the GCM-computed value), [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF] found that the approximations (4-6) reduce the biases from -25.8% to -2.9%, -23.4% to -1.6% and -22.3% to 0.9% for momentum, sensible heat and latent heat, respectively. At this stage, these approximations imply that δU and σ 2 U are the two parameters to be parameterized.

c. Gustiness approach and relevant variables to be parameterized

In contrast to δU , the studies mentioned in the introduction proposed to parameterize the gustiness velocity U g defined as

U 2 = U 2 + U 2 g (7)
and to use U = U 2 + U 2 g in place of U = U +δU in Equations 4, 5 and 6. The two approaches are equivalent in the end as they both lead to mathematically describe U . The only difference lies in the chosen parameter to be parameterized. Both were tested and better scores (as those presented in the analysis of this paper) are found when using the gustiness velocity (not shown). We therefore only present hereafter the results for the gustiness approach. However, as, in contrast to other studies, [START_REF] Bessac | Scale-aware spacetime stochastic parameterization of subgrid-scale velocity enhancement of sea surface fluxes[END_REF] directly parameterize δU , we also briefly discuss in Section c how the change in the parameterized field may impact our results.

Similarly σ U leads to a more skillful parameterization than σ 2 U . Results for σ U are thus synthesized in section e. Note that none of previous studies tackled this component of the meso-scale turbulent surface stress enhancement.

Numerical framework and methods a. Protocol

A coarse-graining approach is used in order to quantify and predict the surface meso-scale wind speed variability at the scale of a GCM grid cell (100 km). The high-resolution simulation (2.5 km) is provided by the operational forecasts from the convection-permitting model (CPM) AROME, in operations at Météo-France, the French weather service (see next subsection for its description). For each GCM-like grid cell (a 100 km×100 km square), the subgrid information (40 × 40 CPM grid points) is thereby explicitly available from the CPM simulation, considered here as a reference. Especially, it allows to diagnose the surface wind speed variability (e.g., U g and σ U ) as well as many other parameters characterizing e.g., convection or large-scale circulations that will be considered to predict the meso-scale wind speed variability (see Section a for the predictors used in the present work). As a reminder (see [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF]), U g is computed following:

U g = U 2 -U 2 (8) with U = √ u 2 + v 2 and U = √ u 2 + v 2 .
u and v are the CPM zonal and meridional wind components, respectively. Bars indicate the horizontal average over the CPM grid-points within the considered GCM-like grid cell. Similarly,

σ U = (U -U ) 2 , with U = √ u 2 + v 2 .

b. Numerical model

The CPM AROME is operating for numerical weather prediction (NWP) activities over several domains covering France and some of the French overseas regions in the tropical belt [START_REF] Faure | Operational implementation of the AROME model in the tropics: multiscale validation of rainfall forecasts[END_REF]. The AROME configuration used in this study has a horizontal resolution of 2.5 km and a 90-level vertical grid with 33 levels located below the 2000-m altitude and a first level at 5 m above the surface [START_REF] Brousseau | Improvement of the forecast of convective activity from the AROME-france system[END_REF]. The non-hydrostatic fully compressible Euler equation system is solved by the AROME dynamical core described in [START_REF] Bubnová | Integration of the fully elastic equations cast in the hydrostatic pressure terrain-following coordinate in the framework of the ARPEGE/aladin NWP system[END_REF]. The AROME physics is detailed in [START_REF] Seity | The AROME-france convective-scale operational model[END_REF] and [START_REF] Brousseau | Improvement of the forecast of convective activity from the AROME-france system[END_REF]

c. Datasets

Two domains which cover mainly open ocean areas are used in this study: (i) the Indien domain which includes La Réunion Island in the tropical Indian Ocean (see Figure 1) and (ii) the Antilles domain which covers the West Indies in the Caribbean Sea and a large fraction of the tropical Atlantic (see Figure S1 in the Supplementary Material). Only GCM-like grid cells with a 100% sea fraction are included in the analysis.

One month of simulation is used on each domain. The month of January 2017 is selected for the Indien domain, as it corresponds to a convectively-active period (expected to be a significant source of meso-scale variability, e.g., [START_REF] Redelsperger | A parameterization of mesoscale enhancement of surface fluxes for large-scale models[END_REF][START_REF] Williams | A physically based parametrization for surface flux enhancement by gustiness effects in dry and precipitating convection[END_REF]). The month of August 2017 is selected for the Antilles domain, also because of its significant convective activity. Hourly model outputs between 12 h and 36 h lead time are used so that meso-scale motions are properly established.

The month of data available for the Indien domain provides 356 280 spatio-temporal GCM-like samples (503 sea points every hour). In this dataset, hereafter referred to as D All Indien , 75 % of the GCM-like samples are randomly sampled and dedicated to the statistical model training (D

D All

Antilles . No difference is observed on the parameterization skills and behaviour. However, for the sake of the clarity and depending on the step in our following analysis, only one of these three testing datasets is used when computing parameterization skill scores. It is clearly indicated whenever relevant. Only the spatial distribution (monthly averaged maps) of the parameterization output biases is based on D All

Indien in order to get a homogeneous temporal sample number for calculating the local statistics for each GCM-like grid cells (e.g., Figure 6). But these diagnostics are only used to illustrate the spatial behaviour of the parameterization.

d. Statistical modelling: penalized linear regression

A penalized linear multi-variate regression is used in order to predict the surface wind speed variability parameters (U g and σ U , referred to as y in this section) from a predefined set of p potential predictors (x j ) j=1,...,p .

For a given set of observed states of size N , the statistical model reads:

y = Xβ + ε. (9)
y is the N -vector of centred and normalized observed responses (y i = 0 and σ y i = 1 for i = 1, . . . , N ), X is the N ×p observed predictor value matrix (also centred and normalized per column), β contains the unknown coefficients of the linear model and ε is the residual vector. β coefficients are usually computed from the ordinary least square (OLS) method.

While efficient for minimizing the prediction bias, a direct applications of the linear regression is likely to lead to over-fitting, especially when the number of a priori predictors is large. The challenge is therefore to achieve a parsimonious model by selecting the most relevant predictors. The least absolute shrinkage and selection operator (LASSO, Tibshirani 1996) is used here to penalize the most complex models (i.e. those with a large number of predictors). In practice, this leads to the cancellation of certain coefficients in β (sparse β). The LASSO problem can be written as:

argmin β 1 2N ∥y -Xβ∥ 2 2 + λ ∥β∥ 1 (10)
where ∥•∥ 1 and ∥•∥ 2 are the l 1 and l 2 norms, respectively. λ is the regularization parameter, which defines the degree of penalization. For each value of λ, the resolution of Equation 10 provides a solution β(λ). Note the l 1 norm allows to get some coefficient of β being exactly equal to zero. The degree of penalization increases with λ, λ = 0 meaning no penalization at all, and λ = 1 corresponding to the highest degree of penalization (i.e. β = 0). The dependence of β j to λ draws the so-called LASSO path (see e.g, Figure 4a for an application in the present work). The LASSO procedure is performed based on the Scikit-learn Python package [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF].

The choice of λ is crucial. A K-fold cross-validation technique [START_REF] Hastie | Statistical Learning with Sparsity: The Lasso and Generalizations[END_REF]) is used to choose λ and thus the most parsimonious linear model. The D Training Indien dataset is randomly divided into K = 100 independent groups. Each group is used once for model testing and the LASSO procedure is run K times on the K -1 remaining groups (model training). The average and standard deviation of the K mean-squared prediction errors (MSE) can be computed for each λ. These are then used to select the λ that minimises the average MSE within one standard deviation (the "onestandard-error rule", see [START_REF] Hastie | Statistical Learning with Sparsity: The Lasso and Generalizations[END_REF]: we keep the λ, which has the smallest value of -log(λ), and for which the average MSE reaches the minimum average value of MSE over all λ plus one MSE standard deviation (see Figure 4 for an application).

The LASSO-based penalization introduces a bias in the resulting linear regression, and therefore is used only to select the relevant predictors. The sparse β is re-calculated through the non-biased OLS procedure using only the LASSO-retained predictors. The intercept is kept in order to compensate, on average, for potential missing predictors that could explain, for instance, a non zero U g when all predictors are zeros.

Note that higher-order predictors (e.g., polynomial functions of the predictors introduced hereafter) or the use of transfer functions (e.g., logarithm) have been considered but not retained as they did not significantly improved the model skills without including a much larger number of predictors (not shown). A simple sparse multivariate linear model based directly on the predictors is therefore preferred.

Wind speed variability characteristics a. Preliminary analysis

In January, southeasterly trade winds prevail on the eastern half of the Indien domain (Fig. 1a). They can reach a wind speed of 8 m s -1 or higher. The near-surface convergence (or minimal wind speed) line positions the InterTropical Convergence Zone (ITCZ, e.g., [START_REF] Waliser | A satellite-derived climatology of the ITCZ[END_REF] in the northern part of the domain and on the eastern part of the Mozambique Channel (see also [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF], for further details). The GCM-resolved wind speed in the D All Indien dataset has a unimodal spatio-temporal distribution with an average wind speed of 5.7 m s -1 and a maximum near 12 m s -1 (Fig. 1b). The distributions of the meso-scale wind speed properties U g and σ U are similar, with a significant positive skewness (Figures 1d,f). The northeastern quarter of the domain contains a large pattern of meso-scale variability (Figures 1c,e), which is generated by two tropical depressions crossing the region during January 2017. In the Mozambique Channel, the meso-scale wind speed variability is associated with regular deep convection activity, specifically along northern coast of Madagascar, where organized deep convection is triggered every day through coastal interactions and spreads offshore [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF]. The region west and southwest of La Réunion island is subject to a regular orographic perturbation (Von Karman vortex street) on the lee-side of the island which generates wind speed perturbations near the surface. In the rest of the domain, sporadic mesoscale features occur, such as isolated deep convective systems, meteorological fronts, squall lines or orographic perturbations. This spatial distribution is similar to the one of the meso-scale flux enhancement (see Figures 3a, 3c and 3e in [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF], in agreement with the fact that the wind speed variability is the main contributor to the meso-scale flux enhancement (Equations 4 to 6).

b. Meso-scale processes and first-guess predictors

Several types of meteorological meso-scale processes generate meso-scale surface wind speed variability (see previous paragraph and [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF]. The idea of an objective classification of these processes in order, for instance, to develop a fully process-based meso-scale wind speed parameterization might be attractive. It has however been discarded due to the difficulties to easily isolate a given process from the others, as superimposition of several mechanisms is often observed. Such a classification is thus left for future work and an intermediate approach is followed hereafter. Yet, we will still use two almost canonical examples to illustrate the involved mechanisms and the parameterization behaviour: (i) a case of isolated convective cells (Figures 2a, c ande) and (ii) a case of a synoptic front (Figures 2b, d andf). Remind that these two purely illustrative examples were hard to find and thus do not necessarily picture the complexity of most of the events of enhanced meso-scale wind speed variability in our dataset.

Case (i) occurred on 21 January 2017 in the middle of the Mozambique Channel. The surface wind speed horizontal distribution and streamlines exhibit the occurrence of several isolated convective cells that develop from around 0200 LT (local time) and last about 8 hours under a roughly homogeneous large-scale northerly wind field (Figure 2a). The meso-scale wind speed variability components U g and σ U increase together at dawn as the convection cells deepen. They in turn modify the surface wind inducing a maximum of meso-scale variability at 0700 LT (Figure 2a, see also the vertical dashed line in the time series in panels c and e). The time evolution of the rainfall rate R, the density current spreading velocity c and the updraft mass flux at the condensation level M up , which are first guess predictors inspired from previously published studies (see later for their detailed definition), are also shown in Figure 2e. Their signal roughly follow the wind speed variability signal. The updraft mass flux increases first -as expected from the life cycle of a convective cell (Houze 2018) -and seems to be the only potential predictor able to explain the meso-scale wind speed variability during the first 2 hours.

Case (ii) occurred on 4 January 2017 in the northeastern quarter of the Indien domain. It consists of a tropical depression centred a few hundreds of kilometres south-west of the selected grid-cell, and which generates a synoptic front whose signature is visible in the wind speed field and streamlines of Figure 2b. For this case, the updraft mass flux signal is weak, with almost no rainfall during the first 12 hours of the enhanced meso-scale wind speed event. Case (ii) thus evidences a potential role of the large-scale dynamics in generating meso-scale wind speed variability, as suggested by the large-scale horizontal wind convergence and wind shear emphasized by the streamline (see also LSS and LSD on Figure 2f, defined in the next section).

Note finally that these two examples, in particular Case (ii), highlight that the gustiness velocity U g and the wind standard deviation σ U do not necessarily peak at the same time.

Parameterization of the meso-scale flux enhancement

In this section, a list of potential predictors is introduced to be used for building a penalized linear regression model of the wind gustiness U g and then of the meso-scale wind standard deviation σ U , following Section d. The obtained statistical models thus objectively select the most relevant predictors and provide the basis of a meso-scale flux enhancement parameterization.

a. Potential predictors

The prior identification of the potential predictors for the meso-scale wind variability is subjective though critical. It is guided here following both earlier studies and the authors' physical intuition about the processes that are likely to contribute to this variability. Note that the analysis of several events such as those presented in Section b further validated the value of some predictors. Mainly two main families of processes are considered in the following, namely processes associated with convection and processes induced by the large-scale dynamics. For each family, a set of scalar metrics is defined to quantify the magnitude of the processes. Note that a given metric may ultimately involve several processes at the same time. Besides, as the derived parameterization is aimed to be ultimately implemented in a GCM, the choice of the metrics is partly driven by what can be available in such models (see also Section 7). Similarly, as most GCM parameterization are uni-dimensional, each metric is supposed to be local in the sense that it only depends on the properties of the atmospheric column above the grid cell in which the meso-scale surface wind variability is diagnosed. Non-local processes such as those induced in the lee side of small islands with orography are thus unlikely to be accounted for (see also Section d).

The 12 potential predictors retained hereafter are listed in Table 1 and further discussed in the two next sub-sections. Their detailed definition is given in Appendix 8. Figure S2 of the supplementary material also provides the probability density function of each predictor over the domain and period considered in the present study, as it may be useful in future studies for e.g., parameterization comparison and implementation.

M up (kg m -2 ) Updraft mass flux X X X X SCAI up (‰) Updraft objects aggregation index X X X X c (m s -1 ) Cold pool spreading velocity X X X X SCAI cp (‰) Cold pool objects aggregation index X X X X R (mm h -1 ) Rainfall rate X X X X clt (-) Total cloud cover X X X X X CAPE (J kg -1 ) Convection available potential energy X X X CIN (J kg -1 ) Convection inhibition X X X Dynamics U (m s -1 ) GCM-resolved wind speed X X X X X dU 850 950 (m s -1 )
Wind shear between 850 and 950 hPa X X X

LSS (s -1 ) GCM-resolved horizontal shear X X X X X LSD (s -1 )
GCM-resolved horizontal wind divergence X X X X X associated updraft mass flux M up is calculated at the lifting condensation level. The choice of a convective updraft detection without any vertical velocity threshold will be discussed in Section c. With this definition, the downdraft mass flux at the same level is exactly the opposite of M up , thus it is not used in the following. The spatial organisation of the convective updrafts may be of importance (e.g., one large updraft vs. many small updraft) and quantified by the simple convective aggregation index (SCAI up [START_REF] Tobin | Observational evidence for relationships between the degree of aggregation of deep convection, water vapor, surface fluxes, and radiation[END_REF]).

Cold pools have been identified as significant contributors to the meso-scale wind variability (e.g., [START_REF] Jabouille | Modification of surface fluxes by atmospheric convection in the TOGA COARE region[END_REF][START_REF] Williams | A physically based parametrization for surface flux enhancement by gustiness effects in dry and precipitating convection[END_REF][START_REF] Hourdin | LMDZ6a: The atmospheric component of the IPSL climate model with improved and better tuned physics[END_REF]. Their intensity is measured by the velocity scale associated to their bulk negative buoyancy (see e.g., [START_REF] Jabouille | Modification of surface fluxes by atmospheric convection in the TOGA COARE region[END_REF][START_REF] Williams | A physically based parametrization for surface flux enhancement by gustiness effects in dry and precipitating convection[END_REF][START_REF] Grandpeix | A density current parameterization coupled with emanuel's convection scheme. part i: The models[END_REF]. The spatial organisation of the cold pools may also be of importance. It is quantified by the simple convective aggregation index SCAI cp , similarly to the convective updrafts.

Finally, a few more standard convection parameters are introduced, namely (i) the surface rainfall rate R averaged over the last hour similarly to [START_REF] Redelsperger | A parameterization of mesoscale enhancement of surface fluxes for large-scale models[END_REF] and [START_REF] Bessac | Stochastic parameterization of subgrid-scale velocity enhancement of sea surface fluxes[END_REF], (ii) the total cloud cover clt as suggested by [START_REF] Zeng | Parameterization of wind gustiness for the computation of ocean surface fluxes at different spatial scales[END_REF], and (iii) the Convective Available Potential Energy (CAPE) and the Convection INhibition (CIN), which characterized the convective instability of the atmospheric column.

ii. Dynamical predictors The large-scale horizontal dynamics is prone to generate meso-scale wind variability, possibly before or without convective activity (see Section b and [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF]. In this study, it is quantified by the large-scale (or GCM-scale) horizontal wind shear modulus (LSS) and the large-scale horizontal wind divergence (LSD).

Vertical wind shear may generate meso-scale eddies and is often used to forecast gust intensity in operational weather prediction models. Therefore, the large-scale, low-level wind shear between 950 hPa and 850 hPa (dU 850 950 , e.g., [START_REF] Bechtold | Parametrization of convective gusts[END_REF] is also considered as a potential predictor.

Finally, following [START_REF] Bessac | Stochastic parameterization of subgrid-scale velocity enhancement of sea surface fluxes[END_REF], the horizontal GCM-resolved wind speed U is introduced to further characterize the large-scale dynamics.

b. Correlations between potential predictors and the gustiness velocity

Before performing the penalized linear regression, the relationship between the parameters listed in the previous section and the gustiness velocity U g is illustrated, based on the D All Indien dataset (Figure 3). U g exhibits the highest Pearson correlations with the updraft mass flux M up (0.87) and the cold pool spreading velocity c (0.81). The associated scatter plots confirm such linear relationships (bottom row on Figure 3), only up to values of U g near 3.5 m s -1 . Note that this range contains 99 % of the data. Aggregation indices of updrafts (SCAI up ) and cold pools (SCAI cp ) also show moderate correlations with U g (0.64 and 0.7, respectively), but the relationships appear weakly linear or display quite a large scatter. Regarding the rainfall rate R, its correlation with U g is also non negligible (0.66), but this dependence to U g is weak up to 2 m s -1 (92 % of the data). For the remaining predictors, only the two dynamical parameters LSS and LSD have correlations larger than 0.5 (0.61 and 0.55, respectively).

If this analysis provides a first insight on the relationship between the potential predictors and U g , their noisiness exhibits that none of these predictors can explain alone and appropriately the spatial and temporal variability of U g in the dataset D All Indien . It thus clearly advocate for a multivariate approach. Besides, the several high correlations possibly hide redundant physical information between predictors, that may lead to overfitting. And even a weakly-correlated potential predictor may convey useful additional information, not accounted for by other predictors. For instance, M up and c are highly correlated (0.83, Figure 3), and thus it remains unclear whether both should be kept in a multivariate regression (the statistical predictor selection in the next section will prove that they should). In contrast, the dynamical parameters LSS and LSD, which are moderately correlated with U g appear rather independent of the other predictors, and thus likely provide useful independent information.

c. Statistical modelling of the gustiness velocity

The statistical modelling presented in Section d is now applied to build a multi-variate linear regression of the gustiness velocity U g , based on the prior predictors detailed in Section a. The LASSO path (Figure 4a) unveils three parameters of prominent weight for all the penalization range: M up , c and LSS. All the remaining parameters have weights lower than ±0.1 whatever the degree of penalization. They thus provide a weakly-relevant contribution to U g or include information that is redundant to M up , c or LSS. The relative importance of the different predictors weakly varies as a function of the regularization parameter, which supports the robustness of the statistical model. The spread of the ( β j ) j=1,...,p coefficients as computed from the K-fold cross-validation is also negligible (hardly visible error bars in Figure 4a, which indicate the K-fold standard deviation).

The optimum degree of penalization given by the cross-validation (vertical dotted line in Figure 4) defines the sparse β vector of the LASSO linear regression, which thus contains five non-zero coefficients: the updraft mass flux M up , the cold pool spreading velocity c, the large-scale horizontal wind shear LSS, the large-scale horizontal wind divergence LSD and the GCM-resolved wind speed U . The five non-zero coefficients β j of the optimal U g statistical model are finally computed using OLS on the D Training Indien dataset and given in Table 2. Relevant predictors include both convective and dynamical parameters. The predictors with the two largest weights are M up and c, and also those with the highest Pearson correlation with U g (Figure 3). Even if the correlation between M up and c is high (0.83), both of them are necessary for a skillful prediction of U g : This probably emphasizes the complementarity of these predictors, which can be active at different stages of the convective system life cycle (e.g., [START_REF] Houze | 100 years of research on mesoscale convective systems[END_REF]). On the one hand, weak updrafts are likely to occur before any cold pool has formed, while, on the other hand, cold pools are likely to remain after the updrafts that generated them have vanish. Besides, cold pool may travel a few hundreds kilometers from the place they were triggered. Given the correlation analysis of Section b, the three other selected parameters, namely LSS, LSD and U would probably not have been chosen, as their correlation with U g was moderate or weak. The statistical method thus reveals that these dynamical parameters bring relevant physical information to predict U g , thereby that the large-scale dynamics is able to generate meso-scale variability of the surface wind. Finally, the updraft mass flux M up has the largest weight in the multi-variate regression. As mentioned above, updrafts are detected as CPM grid cell with a positive vertical velocity, while many previous studies used a threshold of about 1 m s -1 on vertical velocity (e.g., [START_REF] Zipser | Cumulonimbus vertical velocity events in GATE. part II: Synthesis and model core structure[END_REF][START_REF] Igau | Updraft and downdraft cores in TOGA COARE: Why so many buoyant downdraft cores[END_REF]. Several values of this threshold were tested, but our least constrained updraft definition remains clearly superior in terms of derived multi-variate regression (not shown). The upward mass flux generated by vertical velocities between 0 and about 1 m s -1 thus significantly contributes to the meso-scale surface wind variability. So far, the physical understanding of this issue remains elusive to the authors and requires further investigation beyond the scope of the present work. Nevertheless, the visual inspection of a few cases reveals organized features of these weak vertical velocities, thus Antilles dataset. The coefficients are associated to the raw predictors (uncentred and unnormalized) and are presented by decreasing LASSO-weight from the left to the right. A column common factor (power of 10) is only detailed on the column title for the sake of clarity. The first row presents the complete parameterization. Rows 2 to 6 present the coefficients of the simpler versions of the U g multivariate linear regression and the two last rows present two simpler versions, with a pragmatic selection of a predictor subset (see Section 7, "w.o. c" meaning without the parameter related to the cold pools). For each version, the grey shading indicates the predictors that are removed from the initial version. they are not just noise. They may be artefact of the convection-permitting model at the origin of our dataset, but also be associated to the model grey zone of convection processes (e.g., congestus), which are partly resolved and partly unresolved. Similar datasets based on higher-resolution models may help better understand this issue.

β c (×10 -1 ) β LSS (×10 4 ) β LSD (×10 3 ) β U (×10 -2 ) β 0 (×10 -
As mentioned in Section c, a parameterization of U g is preferred to δU , as the resulting parameterization provides better scores (not shown). Nevertheless, if δU is considered as the field to be parameterized, the predictors selected by the LASSO procedure slightly differ, consistently with the results of [START_REF] Bessac | Stochastic parameterization of subgrid-scale velocity enhancement of sea surface fluxes[END_REF] (case n = 1 in their terminology): the rainfall rate R is now selected while the LSD is removed. The resolved wind speed U also becomes a more important parameter.

Finally, note that the present parameterization applied to our dataset does not generate any negative value of U g (U g is by definition positive), even though no such constraint is added in the parameterization. Nevertheless, when used in a GCM, a lower bound of zero is likely to be needed.

d. Testing of the gustiness velocity parameterization

The skills of the proposed U g parameterization is now assessed, using the independent dataset D Test

Indien . The parameterized gustiness velocity is noted U g . Figure 5a shows an overall good behaviour in comparison to the reference, i.e. the parameterization is unbiased and the regression coefficient of determination R 2 = 0.846 is close to one. The distribution shape of the reference gustiness (Figure 1d) is also well accounted for up to about U g = 3 m s -1 (see inset in Figure 2, but for the σ u parameterization 5a) which is the value of quantile 98%. The residual analysis also indicates a good behaviour of the parameterization with only a slight negative bias tendency for large gustiness values, above 2-3 m s -1 (Figure 5b). On the inter-quantile range 2 nd to 98 th (evidenced by the orange lines in 5b), the residual distribution is close to a Gaussian distribution (see inset in Figure 5b). Only the lowest and highest quantile (approximately below the 2 nd and above the 98 th ) evidence an associated tails of the distribution which exhibits more weight than a Gaussian distribution. The origin of the significant negative biases for the largest gustiness values mentioned earlier is discussed in the next paragraph. This shape of the residual distribution highlights limited predictive skill of the parameterization for high gustiness velocities, possibly because of non-linear relationships not accounted for in the present framework, additional relevant metrics not considered here, or poor sampling of these regimes by the current training dataset. The use of the entire dataset D All Indien shows that the derived parameterization is weakly-biased over the Indian Ocean (∼ 0.1 m s -1 ), except over the West and Southwest of La Réunion Island where the parameterization bias reaches -0.8 m s -1 (Figure 6a). The latter region also corresponds to that of highest RMSE and highest frequency of occurrence of residuals being below the 2 nd percentile (Figure 6b). When removing these few GCM grid cells from the available dataset, the negative tendency of the residuals almost vanishes (not shown). The illustration in [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF] suggests that most of these deficiencies are likely attributed to the streets of Von Karman vortices continuously generated during the 21 first days of the time period through an interaction between the prevailing northeasterlies and the high orography of La Réunion Island. As a result, large surface wind speed perturbations occur in the lee side of the island. Such a process is clearly non-local (i.e. induced by processes occurring in remote GCM grid cells), and thus cannot be properly captured by our prior predictors. The definition of the corresponding predictors is left for future work. The proposed parameterization thus captures well the gustiness velocity generated by most of the open ocean processes. This is confirmed when using the D All Antilles dataset as another testing dataset (Figure S1 in the supplementary material). It is reminded that a small intercept is found in the parameterization ( β 0 = 1.92 × 10 -1 , Table 2) and is interpreted as a compensation effect, on average, of non-local missing predictors when all other predictors are zeros. Figure 5: Joint probability distribution (a) between the parameterized gustiness velocity U g and its targeted reference U g and (b) between the parameterization residual ϵ = U g -U g and the reference U g . On panel (b), contours are equally spaced by 0.29 in the logarithmic space starting from 1.3 × 10 -6 , and the two horizontal orange lines indicate the 2 nd and 98 th percentiles of the residual distribution. The inset in both panels shows the quantile-quantile plot (black bullets: quantiles 1% to 99%, grey bullets: quantiles 0.1%, 0.5%, 99.5% and 99.9% and dashed grey line: "1 : 1" line) of (a) the U g distribution against the reference U g distribution (see Figure 1d for the U g distribution) (b) the residual ϵ distribution against the corresponding normal distribution (ϵ N , the normal distribution of ϵ, scaled by the standard deviation of the samples and centered on the samples mean). Only the testing dataset D Test Indien is used here.

β Mup (×10 1 ) β c (×10 -1 ) β LSS (×10 3 ) β SCAIcp (×10 -2 ) β LSD (×10 3 ) β U (×10 -2 ) β 0 (×10 -
(a) (b)

Figure 6: Spatial distribution based on the entire D All Indien dataset of (a) the median of the proposed U g parameterization bias U g -U g (in m s -1 ) and (b) the root mean square error (RMSE colours) and the frequency of occurrence of residuals below the 2 nd percentile of the residual distribution (contours -see also the lowest horizontal orange line in Figure 5b). x and y axis are expressed in ×10 3 km. The parameterization behaviour is further analyzed using the two cases introduced in Section b. The predicted U g matches well the temporal evolution of the reference U g both for the case of isolated convective cells and that of a large-scale synoptic front (Figure 7). In the isolated convective cells case, U g is mainly governed by the updraft mass flux M up and the cold pool spreading velocity c. The cold pool contribution starts about 2 hours after the strong increase of the updraft mass flux contribution, in agreement with the life cycle of convective systems (e.g., Houze 2018). The other input parameters of the parameterization have a negligible contribution. In the front case, convective updrafts significantly contributes all along the gustiness event. As for the first case, cold pools generate wind variability during the second part of the event. But here, the large scale dynamics (LSS and LSD) contribution is significant, the large-scale shear contributing in the same way as convective updrafts, especially during the first part of the event. In both cases, the contribution of the GCM-resolved wind speed U is fairly weak, but always acts as a gustiness damping.

e. Model for σ U

A similar work is performed for parameterizing the meso-scale wind standard deviation σ U , which is required to fully capture the meso-scale momentum flux enhancement. While there is no mathematical relationship between σ U and U g , both parameters are highly correlated (0.84). A simple parameterization of σ U based only on U g would thus be possible. Nevertheless, as discussed later in this section, a multivariate regression proves to be more skillful. The processes generating gustiness and wind variance are assumed similar, so that the 12 predictors of Section a are used to statistically model σ U . The corresponding LASSO path is shown on Figure S3 and the retained predictors and associated weights are provided in Table 2. The five predictors retained for U g are also retained for σ U and a sixth one appears relevant, namely the spatial organisation of cold pools (SCAI cp ). The correlation between the parameterized wind standard deviation σ U and its targeted reference (0.88) is higher than that between σ U and U g , meaning that the multivariate parameterization performs better than a parameterization based on U g only. This is even clearer when comparing the RMSE of the multivariate and U g -only parameterizations (0.228 vs 0.465).

f. Back to the meso-scale flux enhancement

To conclude this section on the added value of the derived parameterization, we come back to the key variables for a GCM, namely the turbulent surface fluxes, and assess the ability of the parameterization to capture the meso-scale flux enhancement. The parameterized meso-scale flux enhancement is identified by the symbol. The reference meso-scale flux enhancement (F φMS , Figures 8a, f andi) shows a similar regional distribution than the gustiness velocity U g (Figure 1c) and the wind speed variance σ 2 U (not shown but similar to that of σ U on Figure 1e). This is expected as these two parameters drive the meso-scale flux enhancement [START_REF] Blein | Meso-scale contribution to air-sea turbulent fluxes at GCM scale[END_REF]. In a "perfect parameterization" framework, in which the true U g and σ U (i.e. as derived from the convection-permitting simulation) is used to compute the meso-scale flux enhancement (noted F PP φMS , PP standing for perfect parameterization), the obtained bias is mostly negligible everywhere on the domain (Figures 8e, h andk). Note that, in the case of the momentum flux, a U g -only perfect parameterization leads to rather large biases and only 57% of the meso-scale momentum flux enhancement is captured on average (Figure 8c). The implementation of the proposed parameterization allows to catch most of the meso-scale flux enhancement for the three surface fluxes (Figures 8d, g andj): on average over the month of data, 80%, 99% and 93% of the meso-scale momentum, sensible heat and latent heat flux enhancements are caught, respectively. As expected from the rather large bias of the parameterized U g (Figure 6b) and σ U , especially near La Réunion Island, surface fluxes also exhibit significant biases there.

Comparison with published U g parameterizations a. U g inter-comparison

The introduction of this work mentions several published parameterizations of the gustiness velocity U g . The associated formulations are synthesized in Table 4. They are applied to the entire test datasets (D Test

Indien

D All

Antilles ) to assess their performance against that of the proposed parameterization (hereafter referred to as the BRV parameterization, according the authors' names). For a fair comparison, the empirical coefficient used in two of them (E1999 and WAPE) is tuned using the training dataset D Training Indien . The published parameterizations show a wide variety of behaviours (Figure 9). In four of them (R2000-R, R2000-M up ,R2000-M down and Z2002), the proposed types of function lead to an unrealistic saturation value at high U g . Several parameterizations (R2000-R and WAPE) also strongly underestimate the frequency of occurrence of weak U g . The B2019 parameterization is unable to generate weak U g values, below 0.5-0.6 m s -1 . In contrast the E1999 and WAPE parameterizations perform rather well over most of the range of U g . In any case, the BRV parameterization performs better than any of the previously-published parameterization (weakest heteroscedasticity and best shape of the residuals). We suppose that the use of a multivariate framework helps the parameterization render a larger range of behaviours. For instance, Figure 7 shows that meso-scale wind variability can arise prior to a precipitating event (and sometimes even without any), so that parameterizations based on R or c only cannot catch these regimes. The skillfull behaviour of the E1999 and WAPE parameterizations is consistent with our use of a linear approach, which includes similar predictors.

The behaviour of the available parameterizations is further synthesized in a Taylor diagram [START_REF] Taylor | Summarizing multiple aspects of model performance in a single diagram[END_REF]) and a normalized bias versus normalized RMSE plot (Figure 10). If some parameterizations show a reasonable U g variability (standard deviation of R2000-M up , R2000-M down , WAPE and E1999), their correlation with the reference is always weaker and their RMSE larger than the BRV parameterization. These diagrams thus highlight the step improvement provided by the present BRV U g parameterization. The spatial distribution of each parameterization bias is also provided in Figure S5 and mainly reveals a horizontally homogeneous bias sign for the most biased parameterization (R2000-M up , R2000-M down , R2000-R, WAPE). Contrarily, the B2019 parameterization have a homogeneous positive bias above regions of weak U g and a negative bias above regions of large U g . Both E1999 and BRV have a less organised residual bias. Besides, as expected, no parameterization is yet able to capture the regional pattern on the lee side of La Réunion Island.

b. Meso-scale flux enhancement

Published parameterizations of U g , together with the BRV parameterization developed here (thus incorporating both U g and σ U ), are now compared in terms of their capability in representing the meso-scale flux enhancement. Taylor diagrams and scatter plots between the normalized bias and the normalized RMSE are provided for the meso-scale enhancement of the momentum, sensible heat and latent heat fluxes on Figures 11 to 13. The monthly-mean spatial patterns of each predicted meso-scale flux enhancement bias are available in the Supplementary Material (Figures S6,S7 and S8). Consistently with the results presented in the previous subsection for U g , the BRV parameterization outperforms in representing the meso-scale enhancement of the three fluxes (weakest bias and RMSE and highest correlation with the reference). The E1999 parameterization exhibits scores often close to those of the BRV-(U g ) parameterization (including only the U g contribution for the momentum flux). The Z2002 parameterization has also weak biases, but poorly captures the variability of the meso-scale flux enhancement. This is also the case for the B2019 parameterization which has similar biases to the BRV-(U g ) parameterization. The B2019 parameterization would probably be improved when considering the subgrid standard deviation of the wind speed.

Simpler versions of the proposed parameterization

The present study developed a parameterization of the meso-scale flux enhancement throughout a parameterization of U g and σ U in a framework where the subgrid-scale variability can be assessed and quantified. A subjective step provided a set of prior physically-based predictors and a statistical method helped identify those that were key to build a multivariate linear regression for each of the two parameters to be predicted (five for U g and six for σ u ). However, in the perspective of a GCM implementation of the present parameterization, all retained predictors will not necessarily be available. Therefore, the performance of simpler versions of the parameterization, which would involve fewer predictors is assessed here. These versions are built by sequentially removing the predictor with the lowest weight in the regression model and re-conducting an OLS optimization using the fewer predictors. For instance, it is expected that the LSS and LSD, though probably easily computable in atmospheric dynamical cores, might not be easily accessible to the GCM physical component. Besides, parameterizations without the cold pool spreading velocity predictor are also proposed, as this parameter is generally not available in GCMs (To our knowledge, only the LMDZ GCM incorporate a cold pool parameterization from which c can be derived, e.g., [START_REF] Hourdin | LMDZ6a: The atmospheric component of the IPSL climate model with improved and better tuned physics[END_REF]. Parameterizations of this predictor may however be considered for implementation (e.g., [START_REF] Grandpeix | A density current parameterization coupled with emanuel's convection scheme. part i: The models[END_REF][START_REF] Pantillon | A parameterization of convective dust storms for models with mass-flux convection schemes[END_REF]. The RMSE and R 2 scores of these intermediate versions are provided in Table 2 and 3, while their corresponding joint probability distributions are available in the Supplementary Material (Figures S9 andS10). The bias of each version is always 4) together with the BRV parameterization derived from Section a-d, using the test datasets D Test Indien D All Antilles . On each panel, the color shading indicates the joint probability distribution between the gustiness velocity simulated by a given parameterization and its targeted reference value, the gray line being the x = y diagonal. The top-right sub-plot exhibits the scatter plot of the residuals, i.e the bias of the parameterized U g against its reference value as a function of the reference U g . The mean bias (in m s -1 ), the determination coefficient R close to zero, as a result of the OLS optimization (not shown). As expected, the full version outperforms the simpler ones. As already emphasized with the E1999 parameterization, the use of the updraft mass flux only (M up ) provides a version that seems already skillful. And as expected from the LASSO path (Figure 4), the successive addition of the cold pool spreading velocity (c) and the large-scale horizontal shear (LSS) provides two major step improvements of the parameterization.

(a) (b) (c) (d) (e) (f) (g) (h)
In particular, the large-scale shear enables the parameterization to capture the weakest but frequent values of U g , below 0.3 m s -1 (Figure S9b andc).

Conclusions

Blein et al. ( 2020) recently argued that the meso-scale enhancement of turbulent surface fluxes at the air-sea interface could be parameterized using the concept of gustiness velocity (e.g., [START_REF] Jabouille | Modification of surface fluxes by atmospheric convection in the TOGA COARE region[END_REF], together with a term including the meso-scale wind speed variance. The present study follows this work and develops a parameterization for these two parameters, through the gustiness velocity U g and the meso-scale wind speed standard deviation σ u .

To develop such a parameterization, a coarse-graining framework based on the operational 2.5km AROME convection-permitting model and gathering groups of 40×40 grid-cells (100×100 km, mimicking a GCM resolution) is used. This numerical set up allows to quantify explicitly the mesoscale (or subgrid in the context of a GCM) variability, thus both the targeted parameters U g and σ u and any relevant predictor, which accounts for the processes which generate such a variability. The available dataset provides two months of hourly forecasts over two oceanic domains covering the tropical Indian Ocean and the western tropical Atlantic Ocean respectively. This large dataset is used in a statistical framework to determine multivariate regressions of U g and σ u , which then serve as the roots of a parameterization of the meso-scale enhancement of surface fluxes. The parameterization sensitivity to the coarse resolution scale may be addressed in a future work.

The prior analysis of two case studies of intense wind speed meso-scale variability indicates that the prior predictors should at least account for both convection processes (e.g., updrafts and cold pools) and large-scale dynamical features that can imprint meso-scale wind variability. The identified processes are often combined in complex ways, so that a statistical approach is preferred to a purely mechanistic one. This prior analysis, combined with the input of previous studies addressing the parameterization of surface gustiness, provides a comprehensive set of twelve potential predictors. The least absolute shrinkage and selection operator (LASSO) then frames a penalized multivariate linear regression approach to objectively keep the most relevant metrics for predicting U g and σ U .

Five predictors are finally selected for U g : the convective updraft mass flux at the condensation level, the cold pool spreading velocity, the large-scale horizontal shear and divergence and the largescale wind speed. The σ U parameterization requires a sixth predictor, namely the cold-pool objects aggregation index. The parameterizations of U g and σ U are described by the regression coefficients listed in Table 2. The proposed parameterizations exhibit high skills and their error distributions marginally deviate from a Gaussian distribution (below approximativelly the 2 nd percentile). The associated negative residual (i.e. an underestimate of the wind-speed variability) is attributed to non-local processes which result from the interaction between the large-scale flow and a substantial topography (e.g. La Réunion Island). Such processes probably can be accounted for, but this is left for future work. Note also that some atmospheric processes may also act on the near-surface wind mesoscale variability with some delay, so that time-lagged predictors may also be of interest for such a parameterization.

The proposed U g parameterization performs significantly better than the previously published parameterizations, as it exhibits no bias, the weakest RMSE and the highest correlation with the reference data. Besides, it also better captures the weak values of the gustiness velocity and does not exhibit any unrealistic saturation behaviour. Note that the σ U parameterization, which does not have any equivalent in the literature, clearly improves the prediction of the momentum flux enhancement. We suggest that the use of a multivariate framework helps our parameterization to render a larger range of behaviour, and thus partly explains its better skills. Finally, the resulting parameterization is able to capture 80%, 99% and 93% of the meso-scale enhancement of the momentum, sensible heat and latent heat fluxes, respectively. In the perspective of the implementation in a GCM, in which some predictors may be unavailable, simpler versions of the propsoed parameterization, i.e. with fewer predictors, are also proposed. Though less skillful than the complete parameterization, the simpler versions perform significantly better than the previously published ones.

Because of the variety of situations explored in our training dataset, the parameterization derived here does not hopefully depend too much on it. Nevertheless, the weights between the selected predictors may depend on the way the convection-permitting model AROME represents or parameterizes some key processes playing in e.g., convection (turbulence, microphysics, resolution). Therefore, similar studies based on different models and domains are encouraged to further consolidate the results presented here, and possibly introduce other relevant predictors accounting for processes not seen in our dataset. Besides, because of the horizontal resolution of AROME, a wide range of processes at play in the boundary layer (e.g., thermals) is still excluded from the present work. This will require a significant endeavour in the future.

Finally, the next step is evidently the implementation and testing of the proposed parameterization in a GCM. Besides a possible practical selection of the predictors, it will require validation of the parameterization input parameters as simulated by the GCM (some of them will clearly rely on other parameterizations such as convection). The implementation may also require specific tuning if one wants to achieve similar meso-scale enhancement effects as quantified here and thus compensate for some of the model errors. In this perspective, the exact value of the regression coefficients provided here should mainly be taken as indicative of the relative weight between processes and, to some extent, may serve as tuning parameters. detection as seen from visual inspection (not shown). This may be related to the AROME model representation of cold pools.

To characterize the cold pool organisation within a GCM-like grid cell, the simple convective aggregation index of [START_REF] Tobin | Observational evidence for relationships between the degree of aggregation of deep convection, water vapor, surface fluxes, and radiation[END_REF] is adapted, following the same definition as for the convective updrafts (see Equation 12) and yielding the metric SCAI cp .

c. Other convection parameters

A few other more standard metrics characterizing the convection intensity are also considered:

• The rainfall rate R (in mm h -1 ) being the average over the last hour and over the GCM-like grid cell of surface liquid precipitation.

• The total cloud cover clt being the instantaneous total cloud cover averaged over the GCMlike grid cell, based on the total cloud cover computing by the CPM AROME. The inclusion of this parameter as a potential predictor follows the study of [START_REF] Zeng | Parameterization of wind gustiness for the computation of ocean surface fluxes at different spatial scales[END_REF].

• The Convective Available Potential Energy CAPE and the Convective INhibition CIN computed from the GCM-resolved thermodynamic state assuming no dilution and a pseudoadiabatic transformation of the lifted air parcel.

d. Large-scale horizontal wind shear and divergence

The large-scale horizontal shear (LSS, in s -1 ), is defined as the modulus of the horizontal wind shear calculated at the AROME first vertical level (5 m) from the four adjacent GCM-like grid cells, based on a centered difference scheme:

Figure 1 :

 1 Figure 1: Monthly average (denoted by ⟨•⟩, left column, x and y axis are expressed in ×10 3 km) on the Indien domain (D All Indien dataset) and the corresponding spatio-temporal probability density function (right column) of: (a) and (b) the GCM-resolved horizontal wind speed at z=5 m ( U , in m s -1 ); (c) and (d) the gustiness velocity (U g , in m s -1 ); (e) and (f) the meso-scale wind speed standard deviation (σ U , in m s -1 ).

  i. Convection-related predictors Convective drafts are commonly characterized by their mass flux. In this study, convective updraft are detected based on a positive vertical velocity and the

Figure 2 :

 2 Figure 2: Two events of enhanced meso-scale wind speed variability over the Indien domain: (left) isolated convective cells in the middle of the Mozambique Channel on 21 January 2017 and (right) meteorological front in the northeastern quarter of the domain on 4 January 2017. For the two cases, (a) and (b) display the horizontal wind speed (color, m s -1) and streamlines at the AROME first vertical level (5 m) at 0700 LT and 1400 LT, respectively, (c) and (d) the time series of the GCMresolved wind speed U , the gustiness velocity U g and the wind speed standard deviation σ U (in m s -1 ) for the corresponding GCM-like grid cell highlighted by the blue square, and (e) and (f) the time series for a few GCM-scale predictors chosen as a first guess: the rainfall rate R (mm h -1 ), the density current spreading velocity c (m s -1 ), the updraft mass flux at the condensation level M up (kg m -2 s -1 ), the large-scale horizontal wind shear LSS (s -1 ) and the large-scale horizontal wind divergence LSD (s -1 ). Vertical dashed line on (c-f) indicates the time of the snapshot of (a-b).11

Figure 3 :

 3 Figure3: Matrix of the Pearson correlation coefficients between all potential predictors and the gustiness velocity U g (filled circle). The bottom row display the joint probability density function (shading) between each predictor and U g , over which the associated U g -bin-averaged scatter plot and standard deviation are superimposed (red dots and bars, U g bins being of 0.5 m s -1 ). The analysis is based on the D All Indien dataset. For the sake of clarity, the CAPE, CIN, LSS and LSD parameters are multiplied by 10 -3 , 10 -2 , 10 5 and 10 5 , respectively.

Figure 4 :

 4 Figure 4: (a) LASSO path of each predictor weight as a function of the degree of penalization λ (logarithmic scale). Error bars on each path, which represent the standard deviation of the weight as derived from the K-fold cross-validation, are so small that they are hardly visible. (b) MSE of each fold of the K-fold cross-validation as a function of the degree of penalization λ (logarithmic scale -gray lines). The black line and error bars exhibit the mean and ±1 standard deviation of all folds. The orange horizontal shading shows ±1 standard deviation around the fold-mean MSE minimum (black arrow). On both panels, the vertical dotted line indicates the degree of penalisation from which the minimum of fold-mean MSE is achieved within ±1 standard deviation, and thus its optimum value.

Figure 7 :

 7 Figure 7: Time series of the parameterized U g (solid black thick line), the reference U g (dotted black thick line) and the different contributing terms of the parameterization (thin lines) for (a) the case of isolated convective cells and (b) that of a large-scale synoptic front discussed in Section b.

Figure 8 :

 8 Figure 8: Monthly average of the meso-scale flux enhancements and associated parameterization biases for momentum (a-e), sensible heat (f-h) and latent heat (i-k): (a, f and i) reference meso-scale flux enhancement ⟨F φMS ⟩, (b, d, g and j) bias of the parameterized meso-scale flux enhancement based on the parameterization derived from Sections a-e, i.e. ⟨ F φMS -F φMS ⟩ and (c, e, h and k) bias of the meso-scale flux enhancement based on a perfect paramterization using the true U g and σ U (⟨F PP φMS -F φMS ⟩). For the momentum flux, two parameterization biases are considered: (b) and (c) only includes the U g contribution ( τ MS (U g ) and τ PP MS (U g )) and (d) and (e) includes the contributions of both U g and σ U ( τ MS (U g , σ U ) and τ PP MS (U g , σ U )).The ⟨•⟩ operator indicates a time average. x and y axis are expressed in ×10 3 km.

Figure 9 :

 9 Figure 9: Intercomparison of the previously-published U g parameterizations (see Table4) together with the BRV parameterization derived from Section a-d, using the test datasets D Test

  2 and the RMSE (in m s -1 ) are also indicated in the upper-left corner of each panel. The considered parameterizations are (a) E1999, (b) R2000-M up , (c) R2000-M down , (d) R2000-R, (e) Z2002, (f) WAPE, (g) B2019 and (h) the BRV parameterization.

Figure 10 :

 10 Figure 10: Statistical scores of previously-published U g parameterizations together with the present BRV parameterization, using the test datasets D Test Indien D All Antilles : (a) Taylor diagram (Taylor 2001, radial distance: standard deviation with the reference represented by the filled black circle along the x-axis; angular position: Pearson correlation coefficient; grey iso-contours: centered RMSE) and (b) scatter plot between the bias and the RMSE of each parameterization, both normalized by the mean reference U g . A common way to read the Taylor diagram on panel (a) is the following: the closer the marker of a parameterization to the reference filled black circle on the x-axis, the better.

Figure 11 :

 11 Figure 11: Same as Figure 10 but for the meso-scale enhancement of surface momentum flux as predicted by the intercompared parameterizations. For the BRV parameterization developed in the present work, two versions are indicated, a first one considering only the use of the gustiness velocity (BRV-(U g )) and a second one including further the impact of the subgrid wind speed standard deviation (BRV-(U g + σ u )).

Figure 12 :Figure 13 :

 1213 Figure 12: Same as Figure 11 but for the meso-scale enhancement of sensible heat flux.

  large-scale horizontal wind divergence (LSD, in s -1 ) at the AROME first vertical level is computed as: Large-scale, low-level vertical wind shearThe large-scale, low-level vertical wind shear dU 850 950 is estimated from the GCM-scale horizontal wind speed difference between 850 hPa and 950 hPa: dU 850 950 = U (p = 850 hPa) -U (p = 950 hPa) (16)

Table 1 :

 1 List of scalar metrics considered to predict the meso-scale wind speed variability. Metrics that have not yet been used in a published parameterization are indicated in bold font. Examples of meteorological features in which the associated processes are likely involved is also given in the right side of the table.

Table 2 :

 2 Coefficients and scores of the multivariate linear regressions of U g . Coefficients are estimated through an ordinary least square optimization using the D Training Indien datasets, while the RMSE and R 2 scores are calculated using the D Test

	1 )	RMSE	R 2

Table 3 :

 3 1 ) Same as Table

	RMSE	R 2

acknowledgements

The authors acknowledge the thoughtful comments of Adam Monahan and an anonymous reviewer, which clearly helped in improving an earlier version of this manuscript. The work of S. Blein was supported by the French national research agency (Agence National de la Recherche) through the ANR-COCOA project "COmprehensive Coupling approach for the Ocean and the Atmosphere" (grant: ANR-16-CE01-0007) and the European Union Horizon 2020 H2020-CRESCENDO project "Coordinated Research in Earth Systems and Climate: Experiments, Knowl-edge, Dissemination and Outreach" funded under the programm SC5-01-2014: Advanced Earth-system models; grant agreement No 641816. This work was also supported by a national funding by the Agence Nationale de la Recherche within the framework of the Investissement d'Avenir program under the ANR-17-MPGA-0003 reference. The INSU-LEFE DEPHY project "Développement et Évaluation des PHYsiques des modèles atmosphériques" allowed for helpful interactions with participating researchers. The authors also thank Aurélien Ribes, Thomas Rieutor and Alix Rigal for their interesting feedbacks.

E1999

Ug = βM down ρσ d

M down in kg m -2 s -1 σ d : downdraft fractional area, fixed at 0.05 β: constant parameter, originally fixed at 10.0, and here tuned at 2.46.

Redelsperger et al. (2000)

R2000-R Ug = log(1.0 + 6.69R -0.476R 2 ) R in cm day -1 R2000-Mup Ug = log (1.0 + 386.6Mup with:

Ψ N,2 is a stochastic residual (centered Gaussian distribution with a variance σ 2 Ψ that scales as a power law of the coarsening scale N (in degree, here N = 1.15 • ). µ = 0.42 and λ = -0.31 Table 4: Synthesis of state-of-the-art U g parameterizations.

Appendix 1: Predictor definitions

This appendix details the computation of the predictors introduced in Table 1. As a reminder (see Section 2), the operator • stands for the average over all the CPM grid cells contained in the considered GCM-like grid cell.

a. Convective updrafts

Updrafts are simply detected as CPM grid cells with a positive vertical velocity. This provides a fraction α u of the GCM grid cell covered by the updrafts and a mean updraft vertical velocity w u . The associated updraft mass flux M up is defined as:

with w the GCM-scale vertical velocity. M up is calculated at the lifting condensation level obtained from the GCM-scale thermodynamics variables. The level of free convection was also tested but did not improve the results (not shown). Previous studies have often used a strictly positive threshold on the vertical velocity to detect convective updraft (e.g., [START_REF] Redelsperger | A parameterization of mesoscale enhancement of surface fluxes for large-scale models[END_REF], but the testing of a few thresholds between 0 to 1 m s -1 led to statistical models of significantly weaker quality (see discussion in Section c). Note also that the downdraft mass flux has also been used in previous studies [START_REF] Emanuel | Development and evaluation of a convection scheme for use in climate models[END_REF][START_REF] Redelsperger | A parameterization of mesoscale enhancement of surface fluxes for large-scale models[END_REF]. In the present study, this is simply the opposite of M up and therefore it is not considered. The meso-scale organisation of convective updrafts is quantified by the simple convective aggregation index (SCAI) proposed by [START_REF] Tobin | Observational evidence for relationships between the degree of aggregation of deep convection, water vapor, surface fluxes, and radiation[END_REF]. It accounts for both the number of convective clusters N and the clusters clumping throughout the mean distance D between each cluster mass center. The SCAI up , expressed in per thousand, is defined as:

with N max the maximum cluster number (here fixed at 800, half the number of GCM-subgrid points) and L is the characteristic length of the GCM-like grid cell (here 100 km). Clusters are defined by an object identification based on watershed algorithm (Scikit-image Python package, van der Walt et al. 2014).

b. Cold pools

Following Grandpeix and Lafore (2010), a velocity scale c is introduced for the horizontal spreading of cold pools (or wakes) due solely to their bulk negative buoyancy. It is related to the Wake Available Potential Energy (WAPE -see also [START_REF] Jabouille | Modification of surface fluxes by atmospheric convection in the TOGA COARE region[END_REF][START_REF] Williams | A physically based parametrization for surface flux enhancement by gustiness effects in dry and precipitating convection[END_REF]; Grandpeix and Lafore 2010) as

where h is the cold pool top height, θ v is the virtual potential temperature and δθ v = θ cp v -θ v with θ cp v the virtual potential temperature averaged over the cold pools. Cold pools are detected between the surface and the layer where the local buoyancy (θ v -θ v ) becomes weaker than -0.008 m s -2 . This threshold is consistent with those used by Tompkins (2001, -0.005 m s -2 ) or Feng et al. (2015, -0.003 m s -2 ). It is slightly more restrictive as weaker values led to erroneous cold pool