Meso‐scale contribution to air–sea turbulent fluxes at GCM scale

Parametrizations of sea surface turbulent fluxes used in general circulation models (GCMs) assume horizontal homogeneity of atmospheric properties at the grid‐cell scale. The present study assesses the contribution of the meso‐scale (i.e., subgrid) to the grid‐scale surface fluxes, for GCM resolution ranging from 20 to 200 km, and thus quantifies the associated GCM surface flux error. A coarse‐graining method allows for an a priori analysis of the subgrid information. It is based on an atmospheric reference dataset produced by the convection‐permitting operational model AROME. The meso‐scale relative contribution to GCM‐scale fluxes exceeding 10% is shown to have large regional patterns, with large values (up to 90%) and high frequency of occurrence (up to 76% of the time). These meso‐scale motions are not necessarily due to convective activity but also occur frequently under dynamical perturbation conditions. Contributions to surface fluxes, at both the GCM scale and the meso‐scale, are disentangled through a Reynolds decomposition. It is found that temperature and humidity meso‐scale heterogeneities do not contribute much to the GCM‐scale fluxes. The subgrid dynamical processes are the main meso‐scale contribution and consist of two parts. The first one represents the wind magnitude and wind direction heterogeneities and corresponds to the so‐called gustiness approach. It is clarified that the contribution of the gustiness wind in the transfer coefficients cannot be neglected. The second part is the contribution of the wind speed subgrid variance, which contributes up to 10% of the GCM‐scale momentum flux.

circulation is sensitive to all types of turbulent exchange at the surface (Trenberth, 1995). The latent heat flux, associated with evaporation, is a key component of the global hydrological cycle (Cao et al., 2015). Heat fluxes also strongly impact the sea surface temperature (SST), whose prediction capability by climate modeling systems is critical for the reliability of climate variation predictions (Seager et al., 1995). In essence, most of the climate variability modes -such as the El Nino Southern Oscillation (ENSO; Timmermann et al., 2018), the Madden-Julian Oscillation (Zhang, 2005) or the West African monsoon (Roehrig et al., 2013;Rodríguez-Fonseca et al., 2015) -are interactive ocean-atmosphere processes, in which surface turbulent fluxes are often key. For instance, the momentum flux (or wind stress) enhances the upper ocean mixing, contributes to the generation of wind-driven currents (key component in the ENSO), and thus has a direct effect on the SSTs (Chen et al., 1994). The surface latent heat flux is also key at smaller scales, such as in tropical cyclone genesis (Zhou et al., 2014;Gao et al., 2019).
The surface turbulent fluxes modelling is commonly based on the Monin-Obukhov similarity theory (MOST; Monin and Obukhov, 1954), in addition to the concept of bulk formula and transfer coefficient. These bulk formulae estimate the surface fluxes from resolved (i.e., timeand space-averaged) quantities and the theoretical transfer coefficients. The MOST is widely used in surface flux parametrizations (e.g., Zeng et al., 1998), which have been improved by focussing on the estimation of transfer coefficients, to which climate models have been shown to be very sensitive (Torres et al., 2018). Challenges are especially found in regimes of weak and strong winds (Kara et al., 2000).
However, the MOST has fundamental hypotheses, which a priori restrict its application: the surface should be flat and homogeneous, the range of stability should be rather moderate, the surface layer should be thick enough to include the first modelling level and the mean atmospheric state should be stationary and horizontally homogeneous. The two latter points emerge because the MOST is meant to describe the turbulent scales of the surface layer, but does not account for wider scale (e.g., meso-scale) motions, whose contribution might be significant for the quantification of surface fluxes at large scale (Sun et al., 1996). Most state-of-the-art numerical weather prediction (NWP) and climate models consider these assumptions as valid at the scale of their time step and of their grid cell size. At best, additional parametrizations are introduced to partly eliminate these assumptions (e.g., the orographic form drag induced by subgrid orography, Beljaars et al., 2004). Among these restrictive conditions, the present study seeks to address the horizontal homogeneity hypothesis as it is often violated at the typical scale of a model grid cell. Atmospheric models running at coarse resolution are of course the most sensitive to this issue, and therefore, this study focuses on the scale of a few tens to hundreds of kilometres.
Based on the Taylor's frozen turbulence approximation (Taylor, 1938), the question of subgrid horizontal heterogeneity is supposed equivalent to that of time series temporal variability. In the 1980s, some researchers had already questioned the assumption of stationarity, with the aim of estimating monthly mean fluxes from long time series (Esbensen and Reynolds, 1981;Hanawa and Toba, 1985;Ledvina et al., 1993). Then, facing the need to represent the effects of meso-scale variability in numerical models, the so-called gustiness approach was introduced in the 1990s. It was first proposed by Godfrey and Beljaars (1991) and Miller et al. (1992), to represent the velocity perturbation at the surface due to the presence of free convection. It is a general simplification of the problem and relies on a quadratic correction of the mean wind speed used in surface flux bulk formulae by a gustiness velocity U g . This approach was then further explored in order, for instance, to better identify subgrid velocity scales (Mahrt and Sun, 1995;Vickers and Esbensen, 1998), to parametrize the wind variability induced by the boundary-layer free convection (Beljaars, 1995;Fairall et al., 1996;Mondon and Redelsperger, 1998;Redelsperger et al., 2000), and to represent the effects of the deep convection on surface fluxes (Jabouille et al., 1996;Redelsperger et al., 2000;Williams, 2001). However, Beljaars (1995) warned that this approach was intuitively and pragmatically introduced and might not be optimal. Today, to our knowledge, the relevance of this approach remains an open question: should we correct only the wind? Is the quadratic formulation the good one? Finally, the parametrization of the gustiness velocity deserves much care. All published works link this parameter to convection processes only (e.g., Redelsperger et al., 2000;Williams, 2001). However, as mentioned in Beljaars (1995), many other meso-scale motions, which remain subgrid in large-scale models, might influence surface fluxes (e.g., meso-scale features of synoptic fronts or tropical cyclones, circulations induced by islands or coastlines). At least their possible impacts at larger scales need to be quantified. Recently, stochastic approaches have been proposed to parametrize the gustiness velocity (Bessac et al., 2019).
Following Redelsperger et al. (2000), it is relevant to distinguish two categories of subgrid (or meso-) scale processes, which fall outside MOST. The smaller scale, typically below 1 km, can be associated with the boundary-layer free convection. The larger scale corresponds to deep convection, but may cover many other meso-scale processes (e.g., Mahrt and Sun, 1995). In the present work, this latter scale is addressed, which will be generically referred to as meso scale in the following. More specifically, we tackle the following two main research questions: (a) what is the impact of meso-scale features on larger scale fluxes?, and (b) how relevant is the pragmatic gustiness approach? Thanks to recent model and computational capacity developments, we are now able to address these questions using models that explicitly represent the relevant meso-scale motions. Here, the operational and convection-permitting NWP model AROME (Seity et al., 2011;Brousseau et al., 2016) is used to quantify the impact of meso-scale atmospheric heterogeneities on the surface fluxes at the scale of a general circulation model (GCM) grid cell. This meso-scale contribution, if neglected in GCMs, corresponds to the error made when applying directly common surface flux parametrizations in GCMs.
The paper is organised as follow. Section 2 provides some background about surface flux parametrizations. It also defines the total contribution of meso-scale heterogeneities to larger-scale surface fluxes, typically at the scale of a GCM grid cell. Section 3 then introduces the numerical framework used to build a large dataset (one month of hourly data at 2.5 km resolution covering two large domains over the tropical Indian Ocean and the Caribbean Sea) which serves as a basis to quantify the meso-scale contribution in Section 4. Both the ocean basin-scale mean quantities and more specific events of strong subgrid-scale heterogeneities are addressed. Section 5 details the Reynolds decomposition used to disentangle the various meso-scale contributions to the GCM-scale surface fluxes and to identify the main governing terms. This provides the basis for an a priori flux reconstruction which is evaluated. Section 6 positions our results in the context of previous studies, then Section 7 concludes this work and elaborates on a few perspectives.

Flux definition
Locally, turbulent eddies generate fluctuations of wind components, potential temperature and specific humidity around their respective mean value obtained when averaging them on a characteristic time-scale T (which must comply with the Reynolds averaging rules; Reynolds, 1895). T is characteristic of the time that is needed for the largest turbulent eddy to travel across an observation point. It separates turbulent and non-turbulent scales (e.g. Sun et al., 1996) and it is related to the length-scale D of the largest turbulent eddy. At an altitude z above the surface, D is of the order of z. The covariances of these fluctuations contribute to the turbulent fluxes of momentum ‖ ‖ = , sensible heat H and latent heat LE. As they act at smaller scales than the vertical grid resolution of any atmospheric model, they have to be parametrized. Near the surface, the turbulent fluxes are commonly described by the MOST. Larger eddies, which originate further from the surface layer, can interact with the local turbulence and contribute to the surface layer fluxes (McNaughton and Brunet, 2002). These larger scales are not considered in the MOST, as it assumes a homogeneous surface and horizontally homogeneous atmospheric averaged fields (at the length-scale D). This theory defines the surface layer as a layer across which the fluxes are considered as constant and thus equal to their surface value which are driven by the friction velocity u * , the temperature turbulent scale * and the humidity turbulent scale q * . Surface fluxes then read (first equality): where a is the near-surface air density, and c pa and L v are the specific heat of moist air and the latent heat of vaporization, respectively. The MOST introduces universal functions for describing the turbulence statistics and mean states. These functions depend on the non-dimensional parameter = z∕L, where z is a reference height and L the Monin-Obukhov length (Monin and Obukhov, 1954). The MOST flux-gradient relationships can be vertically integrated over the surface layer and hence provide the similarity profiles of wind velocity, temperature and humidity. This integration is the basis of the well-known bulk formulae used in surface flux parametrizations (Equation (1), second equality, where C D , C H and C E are transfer coefficients for momentum, sensible heat and latent heat, respectively; Δ = a − s , with representing either U, or q and the subscripts a and s represent the atmospheric value at the reference height and the surface value, respectively).
Even though the MOST (and thus surface flux parametrizations) deals with the wind stress module, assuming that its direction remains constant in the surface layer, atmospheric models need the stress zonal and meridional components x and y , respectively. The usual approach is to assume that the stress direction is given by the surface wind difference ΔU: For the sake of clarity, most of the paper focusses on the behaviour of ‖ ‖ = instead of ( x , y ). The approach by stress components is discussed whenever relevant, and more specifically in Section 5.2.4. More details are provided in the Supporting Information. All variables in Equation (1) (second equality) are implicitly averaged over the characteristic time-scale T. The magnitude of the wind relative to the surface reads: The following generic form of the bulk formulae (1) will be used below: where F is either F U = , F = H or F q = LE. A U , A and A q stand for a , a c pa and a L v , respectively. Flux formulation (4) is crucial in the context of numerical models as it allows for the computation of surface fluxes from mean (or resolved) quantities. The transfer coefficients C mainly depend on roughness lengths and on empirical stability functions (e.g., Businger et al., 1971). Over ocean, the roughness lengths depend on the sea state and thus on the turbulence within the surface layer (e.g. from u * ; Fairall et al., 2003). As a result, iterative numerical methods are often used.
To estimate transfer coefficients, various parametrizations have been developed. Based on the algorithm of Liu et al. (1979) and updated following the TOGA-COARE (Tropical Ocean-Global Atmosphere -Coupled Ocean-Atmosphere Response Experiment, Webster and Lukas, 1992) data analysis, the series of COARE algorithms has been continuously developed (Fairall et al., 1996;Fairall et al., 2003;Edson et al., 2013). Other algorithms have been published and analysed in the intercomparison study of Zeng et al. (1998). All of these parametrizations are based on the bulk approach introduced above. When no wave information is provided, transfer coefficients C thus depend only on the resolved quantities (ΔU, Δ and Δq), so that the fluxes can be written as:

Contribution of the meso-scale to the surface fluxes
The bulk formulation (Equation (5)) can be used to compute the surface fluxes over a wide area, as long as horizontal homogeneity is guaranteed. Hence, if meso-scale motions are resolved at the scale of the numerical model grid cell, the MOST is a priori valid and can be used to compute surface fluxes in the model. However, if meso-scale motions are subgrid, the MOST cannot be applied directly over the whole area, as these meso-scale motions are generally generated out of the surface layer, and consequently are not represented by the theory. Nevertheless, the MOST can still be used at more local scale, before the fluxes are aggregated at the larger scale. For the sake of simplicity, henceforth surface fluxes computed over an area that respects the horizontal homogeneity hypothesis will be called local fluxes.
At the scale of a GCM grid cell (surface S), the total flux is the average F of the local fluxes over the grid cell. It will be referred to as the sampling method (SM) flux and represents the true GCM-scale flux. For any variable X, the average operator reads: In the following, the variables computed from scalar average quantities X will be labelled with thêsymbol. The variables computed from vector average quantities (which are the GCM-computed quantities) will be labelled with thẽsymbol. For the sake of clarity, thêsymbol and thẽsymbol will be replaced by the symbol when they are equivalent.
Generally, GCMs do not compute F as defined above, but rather approximate it by the GCM-computed fluxF , using the resolved parametersΔ , namelyΔU, Δ (=Δ ) and Δq (=Δq): where the transfer coefficientC reads: The GCM wind magnitudeΔU is defined, following the vector average method (VAM), as: The contribution of the meso-scale (MS) to the SM surface fluxes thus reads: and represents the opposite of the error made when applying directly the MOST to a surface flux parametrization in GCMs.

Wind magnitude
It is useful to highlight that the GCM wind speedΔU differs from the wind speed averaged on the GCM grid cell, which is defined by the scalar average method (SAM) as: By definition,ΔU is always lower than or equal to ΔU. The gustiness velocity U g , mentioned in the introduction and intuitively proposed by Miller et al. (1992), is defined from the SAM and VAM results as:

Statistics
The horizontal variability of any given variable X is quantified by its spatial standard deviation X and the associated normalized standard deviation: If ⟨ ⟩ is the monthly average operator, the normalized root mean square error (NRMSE) of the GCM-computed flux is defined as:

Occurrence of meso-scale contribution to flux
The local occurrence of meso-scale contribution to flux is defined as the percent of time with meso-scale relative contribution greater than or equal to 10%.

NUMERICAL FRAMEWORK
In order to assess GCM subgrid heterogeneities and to quantify the meso-scale contribution (Equation (10)) to the GCM-scale surface fluxes (or SM fluxes), a coarse-graining numerical framework has been set up. It is based on a pair of simulations of different resolution but run with the same model and forcing. The coarsest resolution considered is equivalent to a GCM resolution and the finest one is supposed to provide the local (or GCM subgrid) information explicitly. This section describes this numerical framework.

Numerical model
The surface model SURFEX (SURFace EXternalisée; Masson et al., 2013) of the French research and operational community is used. SURFEX is implemented in all Météo-France atmospheric models, from the research model Meso-NH (Lac et al., 2018) to the global climate model CNRM-CM (Voldoire et al., 2019), including operational Météo-France NWP models. SURFEX provides the lower boundary conditions to the atmospheric systems, especially the surface turbulent fluxes. In the present work, SURFEX is run offline, forced at 5 m above the surface by a given atmospheric forcing (Section 3.2). We only consider ocean grid cells. Surface fluxes over ocean can be computed using two different iterative parametrizations: the COARE 3.0 algorithm (Fairall et al., 2003) and the exchange coefficients from unified multi-campaigns estimates (ECUME) algorithm (Belamari and Pirani, 2007;Seity et al., 2011). The iterative COARE 3.0 algorithm uses a formulation of the roughness length that depends on the friction velocity u * . This dependency is guaranteed by including both the contribution of a smooth flow limit and of a wind-wave induced roughness, through the Charnock expression (Charnock, 1955). In contrast, the ECUME parametrization does not include any formulation for the roughness lengths. It directly links the neutral transfer coefficients to the mean velocity, temperature and humidity through empirical functions. Both surface flux parametrizations have options to include the effect of boundary-layer free convection through a gustiness approach (e.g. Redelsperger et al., 2000). These options are not activated in the present work, in order to keep physical consistency with the model that generates the atmospheric forcing (Section 3.2), and which is tuned for operational forecast applications. The activation of such gustiness effects is unlikely to modify our conclusions, since we focus here on the scales resolved by our reference model.

Atmospheric forcing
The atmospheric forcing dataset is provided by the convection-permitting model (CPM) AROME (Seity et al., 2011;Brousseau et al., 2016), which has been operating for NWP activities for the last decade over several domains covering France and some of the French overseas territories (Faure et al., 2020). In the present work, we focus on two oceanic domains covering La Réunion island in the tropical Indian Ocean ( Figure 1) and the West Indies in the Caribbean Sea and a fraction of the tropical Atlantic (Figure 13 below). These AROME configurations have a horizontal resolution of 2.5 km and a 90-level vertical grid (Brousseau et al., 2016), with 33 levels below 2000 m height. Only the atmospheric fields of the first model level (5 m above the surface) are used. The AROME dynamical core solves the non-hydrostatic fully compressible Euler equation system (Bubnová et al., 1995) and thus explicitly resolves most of deep convection features. The subgrid shallow convection is parametrized with the Pergaud et al. (2009) scheme. The subgrid turbulence is tackled with a prognostic turbulent kinetic energy (TKE) equation (Cuxart et al., 2000) and uses the mixing-length formulation of Bougeault and Lacarrere (1989), which allows for a rather realistic subgrid dissipation. AROME is thus also able to capture most other meso-scale circulations, such as those induced by orography and islands and those occurring within fronts or tropical depressions. For the two domains of interest, the initial state and the lateral boundary forcing are provided by the high-resolution model (IFS) of the European Centre for Medium Range Weather Forecasts (ECMWF). The SST field is derived from the Operational Sea Surface Temperature and Ice Analysis (OSTIA) product computed from satellite observations (Donlon et al., 2012). It is constant through each forecast run. The AROME physics are described in Seity et al. (2011) and Brousseau et al. (2016). The AROME model accuracy, in particular in terms of convection process representation (e.g., size and life cycle of convective cells), was improved a few years ago, thanks to the currently used vertical grid (Brousseau et al., 2016). In a model intercomparison study, Field et al. (2017) showed that AROME and other CPMs, run at equivalent resolution, do agree reasonably with observations in term of radiative (broadband) fluxes for a convective case. Faure et al. (2020) also assessed the realism of the AROME configurations used in the present study, in terms of rainfall rate and organisation and showed good performances at fine scales against radar and rain-gauge observations. Hourly output between the 12 to 36 hr lead times are used so that meso-scale circulations are properly established.

Protocol
The numerical framework used below thus makes use of the SURFEX platform. The reference experiment (SFX-CPM) is performed using the same resolution as the original AROME forcing fields (2.5 km). SFX-CPM simulations provide fine-scale fluxes that allow for the true GCM-scale (or reference SM) flux calculation. Then, several experiments are run at coarser resolution, using coarse-grained forcing fields from the AROME model output (fields are averaged over GCM-type grid cells). They still cover the same domain as the one of the SFX-CPM reference simulation. Two different methods are used to deal with the wind velocity. The first one corresponds to what is naturally used in a GCM. It is based on the VAM, which averages wind components (wind magnitudẽ ΔU; Section 2.2). The corresponding simulations are noted SFX-GCM and therefore provide the GCM-computed flux.
The second approach is based on the SAM, which directly averages the wind magnitude (ΔU; Section 2.3). The corresponding simulations are noted SFX-SAM. Each simulation is further described with a dx-X label, where X stands for the targeted coarse resolution, which covers 20, 50, 100, 150 and 200 km. In this study, sea surface currents are neglected: U s = 0 and ΔU = U a , and the notation U will be used for the sake of simplicity. Table 1 synthesises input and output of each simulation. The COARE 3.0 parametrization is used to compute all the surface fluxes used in the following analysis. A sensitivity test with the ECUME parametrization is also done and the associated results are discussed in Section 5.3.1.
Note that the numerical framework accounts for the horizontal variability of SSTs, thus conditioning Δ and Δq. However, the variability of s and q s is an order of magnitude lower than their counterpart in the atmosphere. The horizontal variability of the sea surface properties is beyond the scope of the present study, but the meso-scale variability of atmospheric fields is expected to prevail in that of U, Δ and Δq.
The proposed framework thus consists of an a priori approach: for each coarse grid cell of SFX-GCM and SFX-SAM simulations, the SFX-CPM fine-resolution simulation provides explicit subgrid information. It will be used in Section 5.1 to assess the various meso-scale contributions (from the horizontal heterogeneities) to the GCM-scale surface fluxes.

Case-study
The following results are robust across the two domains and periods tested here (January 2017 for the Indian Ocean domain and August 2017 for the Atlantic Ocean domain). Therefore, we mainly emphasize the results for the AROME Indian Ocean domain ( Figure 1) in the next sections. Results over the AROME Atlantic Ocean domain are briefly discussed in Section 5.3.3. During austral summer, the intertropical convergence zone (ITCZ) tends to shift southward of the Equator, especially in the southwestern Indian Ocean (e.g., Waliser and Gautier, 1993). The southeasterly trade winds over most of the southeastern part of the domain (Figure 1b) encounter the northwesterlies north of Madagascar, enhancing large-scale convection, as underlined by the cumulative rainfall ( Figure 1a). The month of January 2017 is selected for its important convective activity, which is expected to be a significant source of meso-scale heterogeneities (e.g., Redelsperger et al., 2000;Williams, 2001). In the northeastern quarter of the domain (Figure 1a), the rainfall pattern is mainly the signature of two tropical depressions crossing the region during this month. There are also several westerly-propagating squall lines crossing the domain east of Madagascar coasts. Over the Mozambique Channel, Madagascar influences the dominant winds, leading to the ITCZ southward shift there. Over the southern half of the channel, deep convection activity occurs regularly, while the northern half is continuously subject to organised deep convection (influenced by the northwest coast of Madagascar). This latter process will be further analysed in the next section. Figure 2 (light blue histogram) further documents the wind speed distribution that is encountered over the domain and period considered here. It is characterised by a unimodal distribution with an average of 5.8 m⋅s −1 . Configurations with wind speed under 1 m⋅s −1 or above 10 m⋅s −1 are not much sampled.
As the focus is on air-sea interactions, any coarse simulation grid cell with a non-zero land fraction is masked and removed from the analysis (grey cells on Figure 1). Finally, a total of 359,280 spatio-temporal GCM-scale grid cells is available (503 sea points every hour over 1 month) for the dx-100 coarse resolution simulations. Except in Section 5.3.1, the following focuses on this dx-100 resolution.

Mean contribution
The monthly mean meso-scale contribution ⟨ F MS ⟩ to the GCM-scale (dx-100) surface fluxes (Figures 3a, c e), as defined by Equation (10), are approximately correlated to mean rainfall patterns ( Figure 1a). They are also consistent among the three fluxes, although they cover wider regions for momentum. Regions of larger mean contribution are also those where the contribution is most frequently 10% higher than the GCM-scale flux (contours in Figure 3a, c e).
The northeastern quarter of the domain is characterized by a large and frequent meso-scale contribution to the momentum flux: it can reach -0.011 N⋅m 2 and occurs up to 60% of time. The meso-scale relative contribution reaches 20% of the average local SM-flux ( Figure S1a in the F I G U R E 2 Histogram (left axis) of the mean wind speed U for the entire dataset (light blue), for the MS ∕ ≥ 10% subdataset (turquoise), and for the H MS ∕H ≥ 10% subdataset (blue). Box plot (right axis; whiskers show 5 to 95%iles; box shows interquartile range; line shows median; dot shows mean) of wind difference U −Ũ (black), gustiness velocity U g (red) and wind speed standard deviation U (yellow) for the different wind ranges. All diagnostics are based on the dx-100 coarse resolution simulation Supporting Information). Magnitude and occurrence of the meso-scale contribution to the heat fluxes are of lesser importance in this region (Figures 3c,e). A significant meso-scale contribution also occurs fairly frequently in a wide region that covers the southern part of the Mozambique Channel and in the region southeast of Madagascar, even though it is less important than in the northeast quarter of the domain.
The meso-scale contribution is also large in regions of smaller extension. It is particularly high and frequent west of La Réunion for the momentum flux in the lee side of the island. It reaches -0.0127 N⋅m 2 , which is about 20% of the average local GCM-scale flux ( Figure S1a), 76% of time. Regarding sensible and latent heat fluxes, local maxima of the meso-scale contribution are mostly northwest of Madagascar. There, they reach 2.2 W⋅m 2 (∼ 15% of the local average SM-flux; Figure S1b), with an occurrence of 46% for the sensible heat flux, and 16 W⋅m 2 (∼ 10% of the local average SM-flux; Figure S1c), with an occurrence of 44% for the latent heat flux. The processes associated with these large contributions are analysed in detail in Section 4.2.
A comparison between GCM-computed fluxes (from SFX-GCM simulation) and reference SM-fluxes (or GCM-scale fluxes) is further presented for the entire dataset in Figures 3b, d, f. GCM-computed fluxes are systematically lower in amplitude than SM-fluxes, which indicates that the meso-scale contribution systematically enhances the fluxes. It is the most frequent in the low-to-moderate range of the fluxes, but note that large flux events are less sampled in the present dataset. The error can be large even for weak fluxes, thus the meso-scale contribution can be important even under weak wind regimes. For the three fluxes, the regional mean meso-scale relative contribution to fluxes remains weak: 5.5, 3.8 and 2.2% for momentum, sensible heat and latent heat fluxes, respectively.

Example of meso-scale processes enhancing GCM-scale surface fluxes
A comprehensive analysis of the processes at play in the meso-scale organisation of the flow is clearly beyond the scope of this paper. Nevertheless, it is helpful to illustrate a few specific situations during which the meso-scale contribution to the GCM-scale momentum and heat fluxes (and thus to the GCM-computed flux errors) is large, and thus identify at least a few possible underlying mechanisms. To some extent, the choice of the situations detailed below is subjective, although a systematic analysis of many GCM-like grid cells reveals some representative mechanisms which can be encountered over the present region.

Impact of La Réunion orography
La Réunion orography generates a wake process whose impact spreads over a few GCM cells west and southwest of the island, and generates a regional pattern of meso-scale contribution (Figures 3a,c,e) and thus of GCM-computed flux error.
In the following, we focus on the grid cell just west of the island (Figure 4) where the meso-scale contribution to the momentum flux is the largest (Figure 3a). The interaction between the easterly trade winds and the high orography of the island frequently generates a von Kármán vortex street on the lee side of the island (three snapshots in Figures 4a, b, c). Such a process of vortex street generation has been widely studied for other case-studies (e.g. Grubišić et al., 2015). Starting from a purely dynamical meso-scale perturbation of the large-scale flow early in the night, the atmospheric situation becomes convectively active, with the triggering of deep convection near the low-level convergence lines induced by the vortices (cold pool and precipitation patterns; Figure 4, centre and right columns). Such vortices have already been shown able to induce updraughts that may trigger cloud formation (e.g.   Ito and Niino, 2016) and enhance precipitation (Chung and Kim, 2007). Table 2 (Figure 4f), which slightly reduce the wind meso-scale variability (NSTD(U) = 0.33) while increasing significantly the temperature variability (NSTD(Δ ) = 0.30). However the concomitance of the two processes leads to a relative meso-scale contribution to the momentum and sensible heat fluxes similar to those occurring during the early convectively inactive stage ( MS ∼ 45% and H MS ∼ 15%). The latent heat flux behaves similarly to the sensible heat flux (Table 2), with a weak contribution of moisture meso-scale fluctuations.
A large meso-scale contribution to the momentum flux occurs almost every day of January 2017 over this GCM-like grid cell just west of La Réunion (Figure 5a, red shading). By the end of January, easterlies became northeasterlies, so that the von Kármán vortex street moved out from the considered GCM-like grid cell, which yields a more homogenous wind field. Over the whole month, the meso-scale contribution to the momentum flux is highly correlated with the meso-scale variability of the wind, as quantified by NSTD(U) (correlation of 0.8; Figure 6), and with the normalized wind difference (U −Ũ)∕U (correlation of 0.95).
The meso-scale contribution to the sensible heat flux is more intermittent across the month (Figure 5c), although highly correlated to that to the momentum flux (0.85; Figure 6). The meso-scale variability of the temperature difference (NSTD(Δ )) generally peaks when the local convergence lines induced by the von Kármán vortices trigger deep convection, precipitation and ultimately cold pools (moderate correlation of 0.51 between precipitation and NSTD(Δ ). However, it remains weakly correlated to the meso-scale contribution to the sensible heat flux (0.31; Figure 6), which is mostly linked to the wind variability (correlation of 0.50 with NSTD(U) and of 0.92 with (U −Ũ)∕U; Figure 6). The meso-scale contribution to the latent heat flux behaves similarly (correlation of 0.46 with NSTD(Δq) and of 0.99 with (U −Ũ)∕U; Figure 6)

Convergence lines and squall lines offshore of northwest Madagascar
The region to the northwest of Madagascar is characterised by a large meso-scale contribution to the momentum flux (25% of the SM or GCM-scale flux; Figure S1a), which occurs frequently, about 73% of the time. It also corresponds to the highest meso-scale contribution to heat fluxes (15 and 10% of the SM or GCM-scale sensible and latent heat fluxes, respectively; Figures S1b and S1c). In order to describe the mechanism that is responsible for the regional maximum, we focus on the GCM-like grid cell which has the highest meso-scale contribution to the heat fluxes over this region (Figure 3). Figure 7a shows the time series of the relative meso-scale contribution to the three surface fluxes for this grid cell, while Figure 8 highlights some of the relevant fields (U, Δ and H) for a specific event occurring on 22-23 January 2017. Note that the correlation figure for this case, equivalent to Figure 6, is not shown, as the results are similar to those of the previous case.   In the region, convection triggers almost every day of January and organizes as a squall line, parallel to the coast, which propagates from the land to the ocean. The December-January-February climatology of this region shows convective rainfall that is likely driven from this coastal feature (Bergemann et al., 2015; their figure 3a,c) and during night-time, the precipitation pattern is shifted offshore (their Figure 5e). In the early evening (2100 local time, UTC+0300; Figure 8a), the wind blows from the ocean to the land (probably a superimposition between a land-sea breeze and the synoptic wind). It is rather homogeneous, both in terms of direction and magnitude (also Figure 7b). As Δ has also a weak meso-scale variability (Figure 8d and 7c), the sensible heat flux weakly varies within the GCM-like grid cell ( Figure 8g) and the meso-scale contribution to the GCM-scale fluxes is negligible (≤ 3% for the three fluxes; Table 3 and Figure 7a). Six hours later (0300 local time), the convergence line has moved from land to ocean (Figure 8b). Colder air (associated with the nocturnal boundary layer on land and with convective precipitation) is brought over the ocean by the wind blowing from land (Figure 8e), while the occurrence of convection starts to generate cold pools along the convergence line. As a result, NSTD(Δ ) increases from 0.08 to 0.43 (Table 3). NSTD(U) reaches its maximum, mostly because of the convergence line pattern (weak contribution of the cold pool circulation at this time, not shown). It is concomitant with the largest meso-scale contribution to the GCM-scale fluxes (e.g., 92% for the momentum flux; Table 3). Finally, at 0600 (Figures 8c, f, i), the convergence line has moved further offshore and the squall-line organisation of the convection has been enhanced. The associated evaporation of precipitation in the low levels generates large pools of cold air (Figures 8f), with strong winds at their front (Figures 8c). The superimposition between the convergence line and cold-pool gusts keeps large values of NSTD(U) and of the meso-scale contribution to the surface fluxes, though slightly weaker than that at 0300 (47, 10 and 16% for the momentum, sensible heat and latent heat fluxes, respectively; Table 3). Then, the convergence line and the squall-line organized convection further propagate offshore. Therefore they influence a larger region, reaching the Comoros Islands, and thus explain the regional pattern of the meso-scale contribution to GCM-scale fluxes (Figure 3).

Other meso-scale mechanisms enhancing GCM-scale surface fluxes
The previous paragraphs highlighted two meso-scale mechanisms that enhance surface fluxes, namely the interaction between a mean flow and the orography of an island, and the occurrence of a convergence line associated with the land-sea contrast which generates a meso-scale organisation of convection. Such mechanisms are not restricted to the two specific examples detailed above. The first one is probably relevant for many tropical islands or coastlines with significant orography, while the second is presumably important for any coastline with strong land-sea breezes or where convection exhibits similar organisation properties (e.g., Bergemann et al., 2015). Besides, squall lines are often observed over the open ocean (e.g. Jorgensen et al., 1997).
A systematic analysis of the present dataset (not shown) also indicates a significant role of isolated or less organized convective cells (e.g., in the southern Mozambique Channel) as well as of tropical depressions (e.g., northeast quarter of the domain and open ocean east of Madagascar). An example of each these two mechanisms is provided in the Supporting Information. A more objective classification of the processes that lead to a meso-scale enhancement of GCM-scale surface fluxes would probably help in establishing a process-based parametrization for large-scale models. Although it is beyond the scope of the present work, we think such a classification is not straightforward as we have often observed a superimposition of processes in our dataset (e.g., dynamical and convective features). Such difficulties might arguably justify the use of stochastic approaches (e.g., Bessac et al., 2019). Moreover, the observed process superimposition is consistent with the results of Bessac et al. (2019), when they regress the meso-scale flux enhancement on both precipitation (proxy for convection) and the "resolved" (or GCM) flux. In their study, the latter is simplified by  (23) the resolved velocity, which is associated with dynamical processes.

RELEVANCE OF THE GUSTINESS APPROACH
The gustiness approach has been intuitively invoked to deal with subgrid-scale surface heterogeneities associated with the convective activity and which generally increase the surface fluxes. In this framework, this convective contribution to the surface fluxes is parametrized as a correction of the wind speed used in the surface flux bulk formulae (Equation (12)). However, to the authors' knowledge, there is no rigorous or quantitative justification that this approach is able to account for most of the convective effects on surface fluxes. Besides, it is not obvious that such an approach is relevant for other meso-scale processes impacting surface fluxes (see previous section). The present numerical framework provides the opportunity to quantitatively assess the relevance of the gustiness approach to parametrize meso-scale effects on surface fluxes.

Flux decomposition
The previous sections assessed the total meso-scale contribution to surface fluxes. Here, this contribution is further

General case
One should recall that the present study does not address the sea surface current, which is thus supposed to be equal to 0. It leads to the simplification ΔU = U a = U. Let first reconsider the generic form of the surface flux bulk formula from Equation (4). This formula is supposed to be applicable locally when the state variables are horizontally homogeneous above the targeted area. In the presence of subgrid meso-scale heterogeneities, and following the work of Hanawa and Toba (1985) for time averaging, Equation (4) is first decomposed using the spatial Reynolds TA B L E 3 Mean and standard deviation (raw and normalized) of a few variables and meso-scale contribution to fluxes (raw and normalized), computed for the central GCM-like grid cell in Figure 8, for the three considered snapshots decomposition X = X + X ′ for all variables except A ( a , c pa and L E have negligible variations at the considered scales, not shown) and then averaged. Averaged quantities X are computed from Equation (6), and X ′ is the local deviation from the averaged value. This yields: where the l.h.s term represents the flux averaged over a large-scale (or GCM-scale) area, while the r.h.s terms introduce a decomposition of the flux with several field covariances. At this point, the decomposition is similar to the one used, for instance, in the study of Jabouille et al. (1996) or in that of Williams (2001) (discussion in Section 6). In the homogeneous specific case, the mean flux F is simply represented by the mean quantities T 1 A . For the following steps, it is useful to introduce the exchange coefficient that would be computed using mean quantities:Ĉ C differs from the GCM transfer coefficientC (Equation (8)) only by the use of U instead of the GCM wind magnitudeŨ. Hence, the difference between the averaged exchange coefficient and the exchange coefficient computed from the averaged quantities is defined as: The two first r.h.s. terms of Equation (15) can be further decomposed: and

Application to a GCM
As previously mentioned, the GCM is supposed to represent the surface fluxes integrated over each grid cell, which is the SM flux F . However, in general, GCM-computed fluxes are calculated from the resolved quantitiesŨ, Δ and Δq, and therefore correspond to the formulation of Equation (7). Hence, in order to evidence the GCM-computed flux contribution, Equation (18) can be decomposed as: T 1a represents the flux computed by the GCM (F ). T 1b is the contribution from the difference between SAM and VAM wind magnitude calculation within the GCM bulk formula, which impacts both the wind speed and the exchange coefficient. T 1c is the effect of using an exchange coefficient computed from average quantities instead of computing the averaged transfer coefficients.
This decomposition provides evidence that, in the case of subgrid-scale heterogeneities, the GCM-computed fluxF represents only one component of the true GCM-scale flux (i.e., the SM flux). Finally, the meso-scale contributions to surface fluxes F MS (or the opposite of the GCM-computed flux error) read as: The details of the complete decomposition for momentum, sensible heat and latent heat fluxes read as: All the meso-scale contributions introduced by Equations (22) to (24) have been explicitly computed using the numerical framework presented in Section 3. Only T 1 terms have a similar formulation as for the initial bulk formula. The terms T 2 to T 4 correspond to the contributions of the covariance of each variable pair, weighted by the average remaining variable. If it is still possible to physically interpret the terms within T 2 , it is more difficult to do so for the covariances which include C ′ (T 3 and T 4 terms) and for the third-order moment (T 5 term).
The complete decomposition for the momentum flux components x and y is presented in Equations S1 and S2.

Governing terms in case-studies
The relative contributions to surface fluxes (T i ( )∕(−F ∕A ) for i = 1 to 5) are presented for the case induced by La Réunion orography in Figures 5e-g. GCM-computed fluxes (F = T 1a , black curves) deviate from the total fluxes when the meso-scale contributes to fluxes, and for a few punctual events, the meso-scale contributions are so large that the GCM-computed flux (T 1a ) is no longer the dominant term. In most cases, T 1b is the main meso-scale contribution to the surface fluxes. It reaches almost 60% of the total momentum flux in the present example. The other terms are generally weaker and tend to balance each other. The meso-scale relative contributions (T 1b to T 5 terms) are similar for the squall line case offshore Madagascar (Figure 7d).

Governing terms of the meso-scale contribution
Even though the meso-scale contribution can be large when it occurs, and frequent over some regions (Section 4), it remains rather frequently weak (Figures 3b, d, f). In order to better understand the meso-scale contribution to a given flux, the following analysis focuses on subdatasets for which the meso-scale relative contribution to the flux is larger than 10% of its value (F MS ∕F ≥ 10%). This threshold is rather arbitrary, but the main conclusions of this study are not sensitive to this threshold. It leads to a subdataset of 26, 9 and 8% of the entire dataset for momentum, sensible heat and latent heat fluxes, respectively (around 9.3×10 4 , 3.2×10 4 and 2.9×10 4 points). The meso-scale contributions to the momentum flux meeting the previous condition occur mainly at moderate wind speed, around 4.2 m⋅s −1 on average (Figure 2). The meso-scale contributions to both the sensible and latent heat fluxes occur at comparable wind speed, around 3.6 and 3.4 m⋅s −1 on average, respectively, with a shorter distribution tail toward high wind speed.
Based on these subdatasets, the difference between SAM and VAM wind speed calculations within the bulk formula (T 1b term) is found to be the major contribution to the meso-scale variability (Figure 9). It reaches 20, 22 and 24% of the total momentum, sensible heat and latent heat fluxes, respectively. For the momentum flux, the term associated with wind velocity variance (T 2a ) has also a significant contribution of 9%. Then, the sum of T 1c , T 2b , T 3 , T 4 and T 5 leads, on average, to a residual contribution of 1%. For the sensible heat flux, the contribution of the exchange coefficient C H variability (T 1c ) reaches 4%. In contrast, the covariance between wind velocity and temperature (T 2a and T 2b ) only contributes up to 2%. The residual contribution of the remaining terms T 1c to T 5 is less than 1%. Regarding the latent heat flux, even if the contributions of the terms T 1c and T 4 reach 5 and -6%, respectively, the residual contribution of T 1c to T 5 remains below 2%.
As a summary, to reconstruct at least 99% of the true GCM-scale flux (SM-flux) on average, it is necessary to consider, in addition to the GCM-computed contribu-tionF ∕A = T 1a ( ), the terms T 1b ( ) and T 2a ( ) for the momentum flux (the term T 1b ( ) only allows for reconstructing 90% of the GCM-scale flux), the term T 1b (H) for the sensible heat flux and the term T 1b (LE) for the latent heat flux. The GCM-scale fluxes can thus be simplified from Equations (22) to (24) to the following approximations: As a reminder, the term T 1b results from the difference (U −Ũ) between the SAM and the VAM wind magnitude calculations in the bulk formula, both in the transfer coefficient computation and in the wind speed factor in Equation (7). Note that the difference U −Ũ can be associated with the variability of both the wind magnitude and direction whereas T 2a ( ) is associated only with the wind magnitude variability.
The relative contributions of each term to the momentum flux components x and y ( Figure S4) provide evidence of the same governing terms (Equations S3 and S4) as for the momentum flux.
Surprisingly, the subgrid variability of temperature and humidity does not play a significant role in the reconstruction of the sensible and latent heat fluxes. In their analysis of a convective case, Jabouille et al. (1996) found a contribution of 5% for the terms T 2a + T 2b to the sensible heat flux. Here, the mean contribution of the term T 2 reaches only 2%. A careful analysis (quadrant analysis and visual inspections, not shown) of the wind and temperature covariances reveals that, under non-zero horizontal wind conditions, the perturbation due to downdraughts and induced cold pools leads to the systematic vanishing of the covariance. In fact, the wind perturbation U ′ induced by a cold pool is divided into positive and negative contributions depending on the position within the cold pool: (i) leeward, the wind perturbation is downstream and hence positive (U ′ > 0), and (ii) windward, the wind perturbation is upstream and hence negative (U ′ < 0).
Once integrated on the whole cold pool in which the temperature perturbation can be supposed uniform to the first order (Δ ′ < 0), the U ′ Δ ′ covariance generally vanishes. Although observed on average, this mechanism is a simplification of the process and specific events may lead to a residual covariance, as shown by the interquartile range of the term T 2a distribution (Figure 9). The convection case-study of Jabouille et al. (1996) is thus arguably included in the distribution of the term T 2a .
To conclude, the main meso-scale contributions are mostly caused by heterogeneities in the horizontal wind field. Temperature and humidity horizontal heterogeneities do not have, on average, a significant meso-scale contribution to GCM-scale fluxes.

A priori reconstruction of GCM-scale flux
The flux reconstruction of Equations (25) to (27) is further assessed by comparing the true GCM-scale fluxes (or the SM fluxes F ) to the proposed approximations, in the context of the present dataset (a priori flux reconstruction from the CPM data). If we were able to parametrize these terms in a GCM, this a priori reconstruction represents the maximum improvement that could be achieved. Various intermediate flux reconstructions are compared to the GCM-scale (SM) fluxes for both the F MS ∕F ≥ 10% subdatasets and the entire dataset ( Figure 10). The addition to the GCM-computed flux of the difference between the SAM and VAM wind magnitudes in the complete bulk formulae (T 1a + T 1b , red) significantly reduces the initial GCM-computed flux error. For the 10% meso-scale contributions datasets, the residual error in the momentum, sensible heat and latent heat fluxes is reduced to -9.9, -1.6 and 0.9%, respectively (NRMSEs to 14.4, 7.0 and 3.4%, respectively). If the reconstruction T 1a + T 1b (red points) already significantly reduces the momentum flux initial error, the addition of the term T 2a (black points) clearly further improves the flux reconstruction: the mean relative residual is reduced to -2.9%, and the NRMSE to

Considerations for the GCM-scale momentum flux components
At the GCM scale, Equation (2) is valid only in case of horizontally homogeneous field, which means the momentum flux components are only equal to their respective term T 1a . In case of meso-scale contributions (that is replacing x and y by their governing terms), Equations (2) do not hold a priori any more because of covariance terms (Equations (25), S3 and S4). When usingŨ as denominator in the r.h.s. of Equations (2), and thus projecting the corrected momentum flux module onto the GCM mean wind components as usual, the mean biases on x and y reach 2.9 and 2.4% and NRMSEs are of 14.3 and 20.5%, respectively ( Figure S6). The errors are significantly reduced when using the reference wind modulus U as denominator in the r.h.s. of Equation (2) instead ofŨ (mean biases magnitude lower than 1% and NRMSEs of about 5%; Figure  S5). Therefore, the parametrization of meso-scale effect in a GCM should also consider a correction when projecting the momentum flux module in the zonal and meridional directions. This point is further detailed in the Supporting Information.

Wind range dependency
The quality of the flux reconstruction of Equations (25)- (27) is further analysed as a function of the wind intensity ( Figure 11). The largest meso-scale relative contributions are found for the lowest wind subrange (entire dataset; Figure 11b). Based in the 10% meso-scale relative contributions datasets (Figure 11a), the wind magnitude difference between SAM and VAM in the bulk formula (T 1b term) is the main meso-scale contribution for the three fluxes over the entire wind range. For the heat fluxes, this term is mainly due to the U −Ũ wind difference, which decreases as the wind speed increases (Figure 2). For the momentum, the addition of the term T 2a yields an overestimate of the flux at low wind speed (≤ 2.5 m⋅s −1 ), but it remains relevant for higher wind ranges. At low wind speed, the term T 2a ( ) is balanced by the negative terms T 3 ( ) + T 4 ( ) (not shown) in the complete reconstruction (Equation (22)). In fact, the significant negative covariance C ′ D U ′ results from the combined effect of the high wind speed sensitivity of the transfer coefficient at low wind speed and the relative importance of U ′ compared to the mean wind speed (e.g. Figure 2). Below 2.5 m⋅s −1 , U 2 −Ũ 2 = U 2 g (the main variable of T 1b ( ) term) thus corresponds to the main meso-scale contribution. It suggests that a wind dependence should be considered when correcting flux errors by compensating terms. Above 3 m⋅s −1 , the reconstruction has good performance (Figure 2). The GCM-computed flux (T 1a ) is sensitive to the resolution. For the three fluxes, the corresponding correlation with the reference SM-flux moderately increases as the grid gets finer, while the relative standard deviation (RSTD) remains similar across resolutions (Figures 12a, c, e). The mean residuals and NRMSEs are also reduced with finer GCM resolutions (Figures 12b, d, f). The flux reconstructions (T 1a + T 1b for the three fluxes plus T 1a + T 1b + T 2a for the momentum flux only) are not sensitive to the resolution. For the three fluxes, the reconstruction T 1a + T 1b results in correlations with the reference SM-fluxes higher than 0.99 and RSTDs between 0.9 and 1. For the momentum flux, the reconstruction T 1a + T 1b + T 2a further improves both the correlation and the normalized standard deviation.
Overall, the flux reconstruction based on Equations (25)-(27) is thus relevant for a wide range of GCM grid resolutions.

Sensitivity to the surface flux parameterization
Most surface flux parametrizations are based on the bulk approach, and thus the flux decomposition proposed above applies. An additional simulation set is performed following the protocol described in Section 3.3, using the ECUME surface flux parameterization and a 100 km grid resolution for the large grid-cell scale. The governing terms are the same as for the COARE 3.0 algorithm and the flux reconstructions behave similarly (Figure 12, empty thick markers). Only the sensible heat flux reconstruction shows a wider residual of 6%. It is due to a higher sensitivity of the heat transfer coefficient to wind variations (mainly at low wind speed), which leads to a small increase of the term T 1c (H). The GCM resolution sensitivity with the ECUME parametrization is also very weak (not shown), as for COARE 3.0.

Geographical dependence
The regional dependence of the meso-scale contributions to fluxes is analysed using another domain, located in the Caribbean Sea, on which the AROME model is run operationally. The AROME configuration is identical to the one used over the Indian Ocean and the same numerical framework (Section 3) has been set up. The month of August 2017 is chosen for the occurrence of a large convective activity over the southern part of the domain (Figure 13a). Easterly trade winds dominate over the domain (Figure 13b), with a mean wind speed of 7.12 m⋅s −1 (Figure S7), higher than for the Indian Ocean dataset. The subdataset corresponding to meso-scale relative contributions to the momentum flux greater than 10% has also a lower mean wind speed (5.6 m⋅s −1 ) than the entire dataset.
The proportion of the dataset that presents meso-scale relative contributions greater than 10% (9.5, 1.3 and 1.0% of the total of 108720 GCM points for the momentum, sensible heat and latent heat fluxes, respectively) is lower than that for the Indian Ocean configuration. However, as in the previous configuration, the meso-scale contribution has strong regional features and can reach locally large values and frequencies of occurrence (up to 30% of time for the momentum flux; Figure 13c). Once again, the sources of the meso-scale motions are multiple. The main cases that can be identified correspond to isolated convective cells, organized convection, orographic perturbation on the lee side of the Caribbean archipelago and a few tropical depressions.
The analysis of the flux decompositions of Equations (22)-(24) leads to the selection of the same governing terms as for the previous case-study ; Figure S8). It confirms that temperature and humidity meso-scale heterogeneities do not contribute significantly to the surface fluxes and that the variability of the wind magnitude and direction is the dominant meso-scale contribution to the surface fluxes.
A priori reconstruction based on the approximations of Equations (25)-(27), leads to performance similar to those over the Indian Ocean.

Subdataset selection dependence
As mentioned in Section 4.2.3, it is not possible to classify the different situations that lead to meso-scale contributions in a category of identified processes, as they often result from a combination of mechanisms. It is thus hardly possible to assess the sensitivity of the flux reconstruction to the process type. However, the analysis has been conducted on several locations that contain a strong redundancy of an identified process (the ones detailed in Section 4.2) as well as on other locations and specific time periods that contain isolated convective cells ( Figure S2) or tropical depressions ( Figure S3) and no variation on the governing terms has been shown significant between those examples. Meso-scale contributions to surface flux are often solely associated with convective activity Williams, 2001) and thus the precipitation rate becomes a proxy. In the current study, we have also selected a subdataset based on a precipitation criteria. The selected governing terms remain the same as for the minimum meso-scale contribution criteria and thus support the robustness of the former analysis.

DISCUSSION
The question of the meso-scale contribution to larger-scale flux has been addressed in several studies. However, no consensus has been reached yet either in the selection of the governing contributions, or in the way of representing them. A review of previous studies that addressed these questions is summarised in Table 4, where the F approximations proposed by each study is formulated in the present flux decomposition framework (Equations (22)-(24)). As suggested by Redelsperger et al. (2000), it might be relevant to separate by scale category the atmospheric processes that affect the surface field heterogeneity. They suggest a scale separation based on two convection scales, namely the boundary-layer convection scale (of the order of a kilometre) and the deep convection scale (of the order of ten kilometres). This scale separation has also been used by Williams (2001) and, indirectly, in most of the studies since Miller et al. (1992). The current study, based on a 2.5 km resolution CPM, addresses deep convection scales as all other studies except those of Mondon and Redelsperger (1998) and Redelsperger et al. (2000) (boundary-layer part). Note, however, that the present analysis reveals that other meso-scale processes need to be considered in the context of surface turbulent fluxes. The steadiness assumption in the context of the surface flux estimation was first questioned by Esbensen and Reynolds (1981), to assess the time-average surface turbulent fluxes from time-average wind speed, temperature and humidity. Through the analysis of time series, other authors also analyse the representativeness of the fluxes computed from averaged quantities, or the impact of meso-scale motions on the flux estimation (Hanawa and Toba, 1985;Ledvina et al., 1993). The Taylor's frozen turbulence approximation links two equivalent questions: the impact of temporal and spatial meso-scale variabilities on the flux computation (or the representation of the fluxes using time-and space-averaged quantities). Williams (2001) addressed this question in a horizontal 1D framework, using airborne measurements along transects, while Jabouille et al. (1996) followed by Redelsperger et al. (2000) addressed it in a 2D horizontal space of convection-permitting or -resolving simulations. Even if these studies highlight different potentially important meso-scale contributions to the fluxes, all of them point to a pragmatic solution: to simplify these contributions only Jabouille et al.
The F approximations proposed by each study is formulated in the present flux decomposition framework (Equations (23)- (25)). "n.c." stands for not communicated.
by the effect of the wind magnitude difference between the SAM and the VAM in the bulk formulae. Moreover, the addition of this wind difference is usually done through a quadratic decomposition U 2 =Ũ 2 + U 2 g . As recalled by Beljaars (1995), this wind magnitude correction was introduced intuitively for the free convection impact on the mean wind magnitude (firstly by Godfrey and Beljaars, 1991;Miller et al., 1992). In order to address the deep convection scales, Jabouille et al. (1996) "suggest" (in their own words) to use the same wind decomposition, that has later been accepted and used by Redelsperger et al. (2000), Williams (2001) and Zeng et al. (2002). This approach is known as the gustiness correction.
Even though it is barely discussed or mentioned in many studies, and probably neglected in Williams (2001), Mondon and Redelsperger (1998) and Redelsperger et al. (2000), the gustiness velocity can be used also to update the bulk transfer coefficients. Not applying such a correction leads to the modification of the term T 1b = (Ĉ U −CŨ)Δ intoC (U −Ũ)Δ , which thus still uses the GCM transfer coefficientC . This simplification is referred to below as the partial gustiness approach. If it reduces the mean residual compared to the initial GCM-computed fluxes, it generally overestimates the SM-fluxes (Figure 11, brown points), and yields large NRMSEs. This overestimation is particularly large at low wind speed (not shown), where transfer coefficients are the most sensitive to wind speed.
The studies addressing the impact of meso-scale motions on surface turbulent fluxes propose different governing terms (Table 4), except for the latent heat flux for which several studies similarly proposed the gustiness correction (or the partial gustiness correction). In the latter case, the meso-scale contribution to the latent heat flux is represented only by a wind correction in the bulk formula (that is, adding the T 1b (LE) term). Regarding the sensible heat flux, additional terms were shown non-negligible. For instance, Jabouille et al. (1996) considered T 2a + T 2b , which corresponds to an enhancement of the surface fluxes by the temperature meso-scale variability induced by convective downdraughts. Ledvina et al. (1993) found that T 2 + T 3 account for around 40% of the total sensible heat flux in their longest averaging time of 72 hr. Our findings indicate that these terms are on average negligible and probably irrelevant for the parametrization of the meso-scale impact on surface sensible heat flux. In contrast, our results support the gustiness approach of Esbensen and Reynolds (1981) and, in the limit of the partial gustiness simplification, that of Redelsperger et al. (2000) and Williams (2001).
Even though our results indicate that the momentum flux is the most sensitive to horizontal heterogeneities (frequency of occurrence both in space and time, magnitude of the meso-scale contribution), only a few studies have analysed it. Hanawa and Toba (1985) proposed a parametric modification of the term T 1a , while Ledvina et al. (1993) conserved all terms except T 4 . Williams (2001) found the partial gustiness approach to be valid. Hence, no study has addressed the momentum flux decomposition in a context of 2D horizontal heterogeneities. The present study shows that, on average, the gustiness approach fails to represent all the meso-scale contributions to the momentum flux, and that the wind magnitude meso-scale variance (T 2a term) contributes to around 10% of the reference SM-flux. However, at low wind speed (<2 m⋅s −1 ), the term T 2a tends to over-correct the momentum flux, suggesting that a wind dependency may be relevant.
The boundary-layer convection scale is not addressed in the present study. The scale separation proposed by Redelsperger et al. (2000) suggests that this would not influence our present results. The relevance of this hypothesis is beyond the scope of the present study, but we expect that datasets similar to the present one, but resolving explicitly the large eddies of the boundary layer, are becoming available and will enable us to assess it. The work of Ledvina et al. (1993) was the first to analyse the flux decomposition for these scales. In fact, by varying the averaging time from 2 to 72 hr of their time series, they considered scales from the boundary-layer free convection to the deep convection. At the shorter averaging times, they found out that the gustiness approach (equivalent to T 1a + T 1b terms) does approximate correctly both mean heat fluxes LE and H. Regarding the momentum flux, they showed that only the term T 4 has a negligible contribution. The impact of boundary-layer convection scales on the surface fluxes has also been studied by Mondon and Redelsperger (1998) through a large-eddy simulation, but the flux decomposition corresponding to these scales was not detailed. However, the gustiness approach was considered as valid and thus was used. Hence, the gustiness approach is commonly accepted for those scales, but no thorough analysis of the flux decomposition behaviour for boundary-layer scales has been conducted yet in 2D horizontal space.

CONCLUSIONS AND PERSPECTIVES
The effect of meso-scale motions on the air-sea fluxes at the scale of a GCM grid cell is quantified by setting up a numerical framework that uses a CPM as reference. Two different domains, covering the tropical Indian Ocean and the Caribbean Sea are considered, as the forecasts from the operational CPM AROME model are available over long periods. The relatively high resolution (2.5 km) allows for the explicit representation of a wide range of meso-scale motions, such as orographic perturbations and deep convection. Surface and atmospheric variables at the first atmospheric level (5 m) are used to run an offline surface model (SURFEX, which computes surface turbulent fluxes). Surface fluxes at different GCM-scale resolutions (from 20 km to 200 km) are computed following two coarse-graining methods: (a) the GCM-scale flux is the average of the reference fluxes computed at the CPM resolution, and (b) the GCM-computed flux is calculated from the forcing which is averaged at the GCM-scale. The difference between the GCM-scale flux and the GCM-computed flux defines the meso-scale contributions to the flux, which is also the GCM-computed flux error.
At the basin scale, monthly mean patterns of the meso-scale relative contribution to the fluxes are identified, highlighting areas of high frequency of occurrence of relative contribution exceeding a 10% threshold (up to 76% for the momentum flux, and more than 40% for the heat fluxes) and important amplitude (up to 25, 15 and 10% of the GCM-scale fluxes of momentum, sensible heat and latent heat, respectively). Due to few single events, the 5th highest percentile of meso-scale relative contributions reaches more than 70% for the momentum flux, and more than 50% for the heat fluxes. The processes generating these meso-scale contributions patterns are analysed and it is found that orographically induced perturbations, such as under island wakes or coastal interactions, are predominant and generate the largest meso-scale contribution magnitude and contribution occurrence. Organised deep convection, tropical depressions and convergence lines also show significant meso-scale contributions to surface fluxes. Finally, even though the isolated deep convection is a source of meso-scale flux enhancement, this process has an impact that is more spread at the basin and seasonal scales. It has been possible to identify and analyse the atmospheric processes responsible for the main regional patterns of meso-scale contributions on selected examples. However, a systematic attribution of meso-scale contributions to fluxes to specific mechanisms is difficult, given that they often result from a combination of different meso-scale processes.
A common way of representing the meso-scale contributions (the subgrid heterogeneities) to large-scale fluxes is the gustiness approach, which consists in adding to the GCM-computed flux, only the meso-scale contribution due to a wind magnitude correction in the bulk formulae. This wind correction corresponds to the difference between the average of the meso-scale wind magnitude (scalar average method) and the wind magnitude calculated from the average wind components (vector average method). However, no consensus on the relevance of this approach exists in the literature. Using a Reynolds decomposition of the different variables in the bulk formulae, all the contributions, both of the average (large-scale) and fluctuating (subgrid) terms, are quantified. One-month-long simulations are run for two numerical domains, providing around 5 × 10 5 GCM grid cell samples, which constitutes a unique dataset, robust enough for generalizing conclusions on the governing terms representing the meso-scale contributions. The gustiness approach is found to work for representing meso-scale contributions to heat fluxes. On the contrary, this study further revealed that the gustiness approach fails by up to 10% in representing the meso-scale contribution to momentum flux on average, and that a term involving the wind magnitude variance has a significant meso-scale contribution. However, at low wind speed (U < 2 m⋅s −1 ), this term contribution leads to a weak overestimation of the momentum flux. Furthermore, it is proven that the gustiness approach is relevant only if the wind speed correction is considered in the bulk formulae, both in the wind factor and in the transfer coefficients. The partial gustiness approach, which does not consider the contribution of the transfer coefficient sensitivity to the wind, leads to significant flux residuals and to an important spread (mainly at low wind speed). In conclusion, to the first order, and within a good approximation, the meso-scale contribution to surface turbulent fluxes can be parametrized only by representing the wind magnitude and wind direction heterogeneities (included in U −Ũ and U ′2 ). The temperature and humidity meso-scale heterogeneities do not play a significant role, at least on average over the regions analysed in the present study. These conclusions hold for all tested GCM-scale resolutions and for both COARE 3.0 and ECUME surface flux parametrizations.
These conclusions are based on a numerical framework that uses operational forecasts from the 2.5 km resolution CPM AROME, over two tropical domains. Our results thus possibly depend on the model and domains used here and therefore should be further tested with other CPMs and over contrasted regions. The use of relevant observation datasets could also provide further robustness of our results. In further work, it will be of interest to extend the present results to the smaller boundary-layer scale, using large datasets based on large-eddy simulations. The a priori approach also provides a benchmark to evaluate the existing gustiness parametrizations, often based on convective activity, but also to develop more comprehensive parametrizations, which can account for the various sources of meso-scale variability.
Meso-scale heterogeneities over land are even more likely to occur due to, for instance, surface type or soil moisture heterogeneities. The framework proposed in the present work can easily be extended to such contexts.

ACKNOWLEDGMENTS
This work 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, Knowledge, Dissemination and Outreach" funded under the programm SC5-01-2014: Advanced Earth-system models; grant agreement No 641816. 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 would like to thank both referees for their careful reading of this work and their relevant comments. Figure S1: Monthly average SM or true GCM-scale surface fluxes and meso-scale relative contribution to GCM-scale fluxes. Figure S2: Illustration of an isolated convective cells case. Figure S3: Illustration of a tropical low case. Figure S4: Box and whisker plot of the relative contributions of each term of Equations S1 and S2, for the Indian Ocean subdataset. Figure S5 and S6: Joint probability distribution of the reconstruction of the momentum flux components and their relative approximations. Figure S7: Wind distribution for the Caribbean Sea domain. Figure S8: Box and whisker plot of the relative contributions of each term of Equations 22 to 24, for the Caribbean Sea subdataset.