A simplified atmospheric boundary layer model for an improved representation of air–sea interactions in eddying oceanic models: implementation and first evaluation in NEMO (4.0)

A simplified model of the atmospheric boundary layer (ABL) of intermediate complexity between a bulk parameterization and a three-dimensional atmospheric model is developed and integrated to the Nucleus for European Modelling of the Ocean (NEMO) general circulation model. An objective in the derivation of such a simplified model, called ABL1d, is to reach an apt representation in ocean-only numerical simulations of some of the key processes associated with air–sea interactions at the characteristic scales of the oceanic mesoscale. In this paper we describe the formulation of the ABL1d model and the strategy to constrain this model with large-scale atmospheric data available from reanalysis or real-time forecasts. A particular emphasis is on the appropriate choice and calibration of a turbulent closure scheme for the atmospheric boundary layer. This is a key ingredient to properly represent the air–sea interaction processes of interest. We also provide a detailed description of the NEMO-ABL1d coupling infrastructure and its computational efficiency. The resulting simplified model is then tested for several boundary-layer regimes relevant to either ocean–atmosphere or sea-ice–atmosphere coupling. The coupled system is also tested with a realistic 0.25^∘ resolution global configuration. The numerical results are evaluated using standard metrics from the literature to quantify the wind–sea-surface-temperature (a.k.a. thermal feedback effect), wind–current (a.k.a. current feedback effect), and ABL–sea-ice couplings. With respect to these metrics, our results show very good agreement with observations and fully coupled ocean–atmosphere models for a computational overhead of about 9 % in terms of elapsed time compared to standard uncoupled simulations. This moderate overhead, largely due to I/O operations, leaves room for further improvement to relax the assumption of horizontal homogeneity behind ABL1d and thus to further improve the realism of the coupling while keeping the flexibility of ocean-only modeling.

1 Introduction

Owing to advances in computational power, global oceanic models used for research or operational purposes are now configured with increasingly higher horizontal and vertical resolution, thus resolving the baroclinic deformation radius in the tropics (e.g., Deshayes et al., 2013; Metzger et al., 2014; von Schuckmann et al., 2018). Meanwhile fine-scale local models are routinely used to simulate submesoscales, which occur on scales on the order of 0.1–20 km horizontally, and their impact on larger scales (e.g., Marchesiello et al., 2011; McWilliams et al., 2019). By increasing the oceanic model resolution, small-scale features are explicitly resolved, but an apt representation of the associated processes also requires the relevant scales to be present in the surface forcings including the proper interaction with the low-level atmosphere.

1.1 Historical context

Historically, oceanic general circulation models (OGCMs) were forced by specified wind stress and thermal boundary conditions (from observations or reanalysis) independent from the oceanic state, thus often leading to important drifts in model sea surface properties. To minimize such drifts, a flux correction in the form of a restoration of sea-surface temperature and salinity toward climatological values can be added (e.g., Haney, 1971; Barnier et al., 1995). To overcome the shortcomings of the forcing with specified flux, Takano et al. (1973) proposed to use a parameterization of the atmospheric surface layer (ASL) constrained by large-scale meteorological data and by the sea state (essentially the sea-surface temperature and sometimes the surface currents) to compute the turbulent components of air–sea fluxes. Currently, whatever the target applications, such a technique is widely used in the absence of a concurrently running atmospheric model. Such parameterization of the ASL (known as bulk parameterization, e.g., Beljaars, 1995; Large, 2006), which corresponds to a generalization of the classical neutral wall law to stratified conditions (Monin and Obukhov, 1954), is expected to be valid in the first few tens of meters in the atmosphere. In practice, unless a fully coupled ocean–atmosphere model is used, atmospheric quantities at 10 m, either from existing numerical simulations of the atmosphere or from observations, are prescribed as input to the bulk parameterization. Throughout the paper, this approach will be referred to as “ASL forcing strategy”. A problem with such methodology is that the fast component of the system (the atmosphere) is specified to force the slow component (the ocean), whereas the inertia is in the latter. Indeed, a change in wind stress or heat flux will affect 10 m winds and temperature more strongly than sea surface currents and temperature. In the “ASL forcing strategy”, the key marine atmospheric boundary layer (MABL) processes are not taken into account, and thus feedback loops between the MABL and the upper ocean are not represented.

1.2 Air–sea interactions at oceanic mesoscales

An increasing number of studies based either on observational studies and/or on air–sea coupled simulations have unambiguously shown the existence of air–sea interactions at oceanic mesoscales (e.g., Giordani et al., 1998; Bourras et al., 2004; Chelton and Xie, 2010; Frenger et al., 2013; Schneider and Qiu, 2015; Oerder et al., 2016). Those interactions affect the mass, heat, and momentum exchange between the atmosphere and the ocean. We focus in this work on the dynamical response of the surface wind stress to the sea-surface properties (sea surface temperature (SST) and currents) which directly affects low-level winds, temperature, and humidity. Several mechanisms responsible for the surface wind-stress response to SST and oceanic currents can be invoked:

• i.

  Downward momentum mixing. SST-induced changes in the stratification produce significant changes of wind speed and turbulent fluxes throughout the MABL with an increase (decrease) in wind speed over warm water (cold water). As the wind blows over warm water, the MABL becomes more unstable, which leads to an increased vertical mixing, resulting in a downward mixing of momentum from the upper atmosphere to the surface strengthening surface winds on the warm side of an SST front (e.g., Wallace et al., 1989). This mechanism results in a proportional relationship between wind-stress intensity and SST mesoscale anomalies which has been identified in observations and coupled simulations (e.g., O'Neill et al., 2010; Oerder et al., 2016). Considering spatial derivatives of this proportional relationship leads to a correlation between wind-stress divergence (curl) and downwind (crosswind) SST–gradient (e.g., Chelton and Xie, 2010; Schneider and Qiu, 2015).

• ii.

  Atmospheric pressure adjustment. This mechanism corresponds to an adjustment of the atmospheric pressure gradient to the underlying SST, which manifests itself as a linear relation between horizontal wind divergence and the Laplacian of SST (Lindzen and Nigam, 1987; Minobe et al., 2008; Lambaerts et al., 2013).

• iii.

  Oceanic current feedback. The momentum exchange between the ocean and the atmosphere is also largely affected by a dynamical coupling through the dependence of surface wind stress on oceanic surface currents (e.g., Dewar and Flierl, 1987). This coupling results in a drag exerted by the air–sea interface on the ocean which leads to a systematic reduction of the wind power input to the oceanic circulation.

Even if these three mechanisms are mainly active at oceanic eddy scales, they can induce significant effects at larger scales in regions with large SST gradients and/or surface currents (Hogg et al., 2009; Bryan et al., 2010; Renault et al., 2016a). They jointly leave their imprint on the wind divergence, and identifying the relative importance of each mechanism on the momentum balance is difficult because it depends on the dynamical regime and on the spatial and temporal scales of interest (Schneider and Qiu, 2015; Ayet and Redelsperger, 2019).

In the ASL coupling strategy the pressure adjustment mechanism is absent, and only a small fraction of the downward momentum mixing mechanism is accounted for through the modification of the surface drag coefficient depending on the ASL stability (Businger and Shaw, 1984; Chelton and Xie, 2010). As far as the current feedback is concerned, Renault et al. (2016b) showed that the reduction of wind power input to the ocean is systematically overestimated in oceanic simulations based on an ASL forcing strategy compared to air–sea coupled simulations. A simulation that neglects the MABL adjustment to the current feedback cannot represent the partial re-energization of the ocean by the atmosphere and hence overestimates the drag effect by more than 30 % (e.g., Renault et al., 2016b, 2019a). The ASL forcing strategy used in most oceanic models will thus overestimate the current feedback effect and underestimate the downward momentum mixing.

1.3 The proposed approach and focus for this paper

The various aspects discussed so far suggest that a relevant coupling at the characteristic scales of the oceanic mesoscales requires nearly the same horizontal resolution in the ocean and the atmosphere (since the atmosphere must “see” oceanic eddies and fronts) as well as an atmospheric component more complete than a simple ASL parameterization to estimate air–sea fluxes. This assessment raises numerous questions on current practices to force oceanic models across all scales¹ in the absence of an interactive atmospheric model. The computational cost associated with the systematic use of fully coupled ocean–atmosphere models of similar horizontal resolution is generally unaffordable and comes with practical issues like the proper definition of initial conditions via data assimilation techniques (e.g., Mulholland et al., 2015) and the proper choice of a parameterization set. Moreover, in the fully coupled case at basin or global interannual scales the temporal consistency with the observed variability is generally lost unless a nudging toward observations or reanalysis is done in the atmosphere above the MABL (e.g., Bielli et al., 2009). There is thus clearly room for improvement in the methodology to compute the surface boundary conditions for an ocean model. Alternatives to the ASL forcing strategy have already been suggested by Kleeman and Power (1995), Seager et al. (1995), and Deremble et al. (2013). They proposed a vertically integrated thermodynamically active and dynamically passive MABL model where the wind and the MABL height are specified as in the current practices. Such a model allows a better feedback between SSTs and low-level air temperature and humidity because the latter two are prognostic (Abel, 2018). However, by construction, such models do not reproduce the various aforementioned coupling mechanisms affecting the surface wind stress. Their focus is on the improvement of the large-scale thermodynamics while ours is on the improvement of the eddy-scale momentum exchanges. In the present study we propose an alternative methodology to improve the representation of the downward momentum mixing and of the current feedback effect in ocean-only simulations, leaving aside the pressure gradient adjustment from now on. Our aim is to account for the modulation of atmospheric turbulence by anomalies in sea-surface properties in the air–sea flux computation, which is thought to be the main coupling mechanism at the characteristic scales of the oceanic mesoscales. As a step forward beyond the ASL forcing strategy we propose to complement the ASL parameterization with an ABL parameterization while keeping a single-column frame. By construction our approach excludes horizontal advection whose effect can be important in the vicinity of strong SST fronts (e.g., Kilpatrick et al., 2014; Ayet and Redelsperger, 2019). However, we considered that finding a simple and efficient MABL parameterization is the top priority to start investigating the viability of our approach in terms of practical implementation and computational cost. Indeed there exists a large variety of parameterization schemes to represent the effects of subgrid-scale turbulent mixing in the ABL (see LeMone et al., 2019, and references therein). The schemes based on a diagnostic or prognostic turbulent kinetic energy (TKE) are very popular for operational and research purposes despite well-identified shortcomings (e.g., Baklanov et al., 2011). For our purposes we do not need the full complexity of the schemes used in practice in atmospheric models because aspects like cloud processes and complex terrains are outside our scope. For this reason, the guideline in this paper is the development and the testing of a simplified version of the TKE-based scheme proposed by Cuxart et al. (2000) for over-water and over-sea-ice conditions. Note that the single-column approximation for our simplified model selected in this study is only a temporary choice to provide evidence on the viability of the whole approach. More advanced formulations allowing a more realistic momentum balance (i.e., including advection) to be recovered will be studied in future work.

1.4 Content

The objective of the present study is to introduce a simplified model of the MABL of intermediate complexity between a bulk parameterization and a full three-dimensional atmospheric model and to describe its integration to the Nucleus for European Modelling of the Ocean (NEMO) general circulation model (Madec, 2012). This approach will be referred to as the “ABL coupling strategy”. A constraint in the conception of such a simplified model is to allow an apt representation of the downward momentum mixing mechanism and partial re-energization of the ocean by the atmosphere while keeping the computational efficiency and flexibility inherent to ocean-only modeling. The paper is organized as follows. In Sect. 2, we describe the continuous formulation of the simplified model called ABL1d, including the parameterization scheme used to represent vertical turbulent mixing in the MABL and the strategy to constrain this model with large-scale atmospheric conditions. Section 3 provides the description of the dicretization and of the practical implementation of the ABL1d model in the NEMO framework. In Sects. 4 and 5 numerical results obtained for some atmosphere-only simplified test cases available in the literature and for a coupled NEMO-ABL1d simulation in a global configuration are shown. Finally, our conclusions and perspectives are summarized and discussed in Sect. 6.

2 Model equations

In this section we first provide some basic elements on model reduction to motivate our approach and mention possible alternatives (Sect. 2.1). Then we detail the continuous formulation of the ABL1d model and discuss the assumptions made. In particular the governing equations and necessary boundary conditions are given in Sect. 2.2 and the turbulence closure scheme for the MABL in Sect. 2.3. Finally in Sect. 2.4 we discuss the methodology to relax the ABL1d prognostic quantities toward large-scale data.

2.1 Motivations and proposed approach

Global oceanic models can be run at higher resolution than global atmospheric models because of their affordable computational cost. From an oceanic perspective, we generally simulate at high resolution (in space and time) ocean fields ${\mathit{\varphi }}_{\mathrm{HR}}^{\mathrm{oce}}(x,y,z,t)$ over a time interval $t\in [\mathrm{0},T]$ over which only large-scale atmospheric data ${\mathit{\varphi }}_{\mathrm{LS}}^{\mathrm{atm}}(x,y,z,t)$ are known from the integration of a model ℳ_atm using lower-resolution surface oceanic data ${\mathit{\varphi }}_{\mathrm{LS}}^{\mathrm{oce}}(x,y,z=\mathrm{0},t)$ to compute its surface boundary conditions, namely

$${\displaystyle }{\mathit{\varphi }}_{\mathrm{LS}}^{\mathrm{atm}}(x,y,z,t)={\mathcal{M}}_{\mathrm{atm}}\left({\mathit{\varphi }}_{\mathrm{LS}}^{\mathrm{oce}}\right(x,y,z=\mathrm{0},t\left)\right),t\in [\mathrm{0},T].$$

Instead of directly using ${\mathit{\varphi }}_{\mathrm{LS}}^{\mathrm{atm}}(x,y,z=\mathrm{10}\mathrm{m},t)$ to constrain the oceanic model as in the ASL forcing strategy, our objective is to estimate (without running the full atmospheric model again) the correction to the 10 m large-scale atmospheric data associated both with the fine resolution in the oceanic surface fields and with the two-way air–sea coupling. Somehow we aim to find a methodology to get a cheap estimate ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{HR}}^{\mathrm{atm}}(x,y,z=\mathrm{10}\mathrm{m},t)$ of the solution that would have been obtained using a coupling of ℳ_atm and the oceanic model at high resolution. To do so we could imagine several approaches: (i) estimate $\frac{\partial {\mathcal{M}}_{\mathrm{atm}}}{\partial {\mathit{\varphi }}_{\mathrm{LS}}^{\mathrm{oce}}}$ (i.e., the derivatives of the atmospheric solution with respect to the oceanic parameters) via sensitivity analysis which would require to have the possibility to operate ℳ_atm; (ii) build a surrogate model via learning strategies which would require a huge amount of data and computing time; (iii) select the feedback loops of interest and define a simplified model to mimic the underlying physical mechanisms. Following the terminology of Razavi et al. (2012), the first two approaches enter the class of statistical or empirical data-driven models emulating the original model responses while the third one enters the class of low-fidelity physically based surrogates which are built on a simplified version of the original system of equations. In the present study we consider this latter approach, in the spirit of Giordani et al. (2005), who derived a simplified oceanic model by degenerating the primitive equation system and prescribing geostrophic currents into the momentum equation in substitution of the horizontal pressure gradient. In this model, a simple 1D oceanic mixed layer is three-dimensionalized via advective terms to couple the vertical columns with each other. The idea here is to translate this idea to the MABL context. In the rest of this section we describe the continuous formulation of our simplified MABL model which will be referred to as ABL1d.

2.2 Formulation of a single-column approach

The formulation of the ABL1d model is derived under the following assumptions: (i) horizontal homogeneity (i.e., ${\partial }_x\cdot ={\partial }_y\cdot =\mathrm{0}$); (ii) the atmosphere in the computational domain being transparent (i.e., ${\partial }_z\mathcal{I}=\mathrm{0}$ with ℐ the radiative flux) meaning that cloud physics is ignored and solar radiation and precipitations at the air–sea interface are specified as usual from observations (e.g., Large and Yeager, 2009); (iii) vertical advection being neglected. Such assumptions prevent the model from prognostically accounting for the SST-induced adjustment of the atmospheric horizontal pressure gradient and for horizontal advective processes associated with a higher resolution boundary condition at the air–sea interface. The focus here is on the proper representation of the modulation of the MABL turbulent mixing by the air–sea feedback, which is thought to be the main coupling mechanism at the characteristic scales of the oceanic mesoscales impacting ${\mathit{\varphi }}^{\mathrm{atm}}(z=\mathrm{10}\mathrm{m},t)$ and hence air–sea fluxes. This mechanism is expected to explain most of the eddy-scale wind–SST and wind–current interactions and is key to properly downscaling large-scale atmospheric data produced by a coarse-resolution GCM to the oceanic resolution. At a given location in space, the ABL1d model for the Reynolds-averaged profiles of horizontal velocities u_h(z,t), potential temperature θ(z,t), and specific humidity q(z,t), given a suitable initial condition, is

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{rcl}{\partial }_t{\mathbf{u}}_{\mathrm{h}}& =& -f\mathit{k}\times {\mathbf{u}}_{\mathrm{h}}+{\partial }_z\left(K_{\mathrm{m}}{\partial }_z{\mathbf{u}}_{\mathrm{h}}\right)+{\mathbf{R}}_{\mathrm{LS}}\\ {\partial }_t\mathit{\theta }& =& {\partial }_z\left(K_{\mathrm{s}}{\partial }_z\mathit{\theta }\right)+{\mathit{\lambda }}_{\mathrm{s}}({\mathit{\theta }}_{\mathrm{LS}}-\mathit{\theta })\\ {\partial }_{t}q& =& {\partial }_z\left(K_{\mathrm{s}}{\partial }_{z}q\right)+{\mathit{\lambda }}_{\mathrm{s}}(q_{\mathrm{LS}}-q)\end{array}\right.\end{array}$$

for the height z between a lower boundary z_sfc and an upper boundary z_top, which will be considered horizontally constant because only the ocean and sea-ice-covered areas are of interest. In Eq. (1), $\mathit{k}=(\mathrm{0},\mathrm{0},\mathrm{1}{)}^t$ is a vertical unit vector, f is the Coriolis parameter, K_m and K_s are the eddy diffusivity for momentum and scalars respectively, the subscript LS is used to characterize large-scale quantities known a priori, λ_s(z,t) is the inverse of a relaxation timescale, and R_LS denotes a large-scale forcing for the momentum equation. R_LS can either represent a forcing by geostrophic winds u_G (i.e., ${\mathbf{R}}_{\mathrm{LS}}=f\mathit{k}\times {\mathbf{u}}_{\mathrm{G}}$) or equivalently by a horizontal pressure gradient (i.e., ${\mathbf{R}}_{\mathrm{LS}}={\left(\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{a}}}{\mathbf{\nabla }}_{\mathrm{h}}p\right)}_{\mathrm{LS}}$) combined with a standard Newtonian relaxation (i.e., ${\mathbf{R}}_{\mathrm{LS}}={\mathit{\lambda }}_{\mathrm{m}}({\mathbf{u}}_{\mathrm{LS}}-{\mathbf{u}}_{\mathrm{h}})$). Because of the simplifications made to derive the ABL1d model the R_LS term and a nonzero λ_s are necessary to prevent the prognostic variables from drifting very far away from the large-scale values used to “guide” the model. By itself a relaxation term does not directly represent any real physical process, but the rationale is that it accounts for the influence of large-scale three-dimensional circulation processes not explicitly represented in a simple 1D model. Note that this methodology is currently used to evaluate GCM parameterizations in 1D column models. Once the turbulent mixing and the Coriolis term have been computed to provide a provisional prediction ${\mathit{\varphi }}^{n+\mathrm{1},\star }$ at time n+1 for any ABL1d prognostic variable ϕ, the relaxation term provides a weighting between this prediction and the large-scale quantities:

$$\begin{array}{}\text{(2)}& {\mathit{\varphi }}^{n+\mathrm{1}}=\mathrm{\Delta }t\mathit{\lambda }{\mathit{\varphi }}_{\mathrm{LS}}+(\mathrm{1}-\mathrm{\Delta }t\mathit{\lambda }){\mathit{\varphi }}^{n+\mathrm{1},\star },\end{array}$$

with Δt the increment of the temporal discretization. Above the boundary layer, the ABL1d formulation is unable to properly represent the physics; therefore the λ parameter should be large, while in the first tens of meters near the surface we expect the ABL1d model to accurately represent the interaction with the fine-resolution oceanic state, and thus the relaxation toward ϕ_LS should be small. The exact form of the λ_s and λ_m coefficients is discussed later in Sect. 2.4. Note that because of the relaxation term, three-dimensional atmospheric data for u_LS, θ_LS, q_LS, and possibly ${\left(\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{a}}}{\mathbf{\nabla }}_{\mathrm{h}}p\right)}_{\mathrm{LS}}$ sampled between z_sfc and z_top must be provided to the oceanic model instead of the two-dimensional data (usually at 10 m) necessary for an ASL forcing strategy. Since the ABL1d model does not include any representation of radiative processes and microphysics, the radiative fluxes and precipitation at the air–sea interface are similar to the one provided for a standard uncoupled oceanic simulation. The model requires boundary conditions for the vertical mixing terms which are computed via a standard bulk formulation:

$$\begin{array}{ll}\text{(3)}& {\displaystyle }& {\displaystyle }{\left.K_{\mathrm{m}}{\partial }_z{\mathbf{u}}_{\mathrm{h}}\right|}_{z=z_{\mathrm{sfc}}}=C_D\Vert {\mathbf{u}}_{\mathrm{h}}\left(z_{\mathrm{sfc}}\right)-{\mathbf{u}}_{\mathrm{oce}}\Vert \left({\mathbf{u}}_{\mathrm{h}}\right(z_{\mathrm{sfc}})-{\mathbf{u}}_{\mathrm{oce}}),\text{(4)}& {\displaystyle }& {\displaystyle }{\left.K_{\mathrm{s}}{\partial }_z\mathit{\varphi }\right|}_{z=z_{\mathrm{sfc}}}=C_{\mathit{\varphi }}\Vert {\mathbf{u}}_{\mathrm{h}}\left(z_{\mathrm{sfc}}\right)-{\mathbf{u}}_{\mathrm{oce}}\Vert \left(\mathit{\varphi }\right(z_{\mathrm{sfc}})-{\mathit{\varphi }}_{\mathrm{oce}}),{\displaystyle }& {\displaystyle }\mathrm{with}\mathit{\varphi }=\mathit{\theta },q.\end{array}$$

For the sake of consistency, it is preferable to use a bulk formulation as close as possible to the one used to compute the three-dimensional large-scale atmospheric data ${\mathit{\varphi }}_{\mathrm{LS}}^{\mathrm{atm}}$. Because in the present study the plan is to use a large-scale forcing from ECMWF reanalysis products, we use the IFS (Integrated Forecasting System: https://www.ecmwf.int/en/forecasts/documentation-and-support/changes-ecmwf-model/ifs-documentation, last access: 20 January 2021) bulk formulation such as implemented in the AeroBulk (https://brodeau.github.io/aerobulk/, last access: 20 January 2021) package (Brodeau et al., 2017) to compute C_D, C_θ, and C_q in realistic simulations (see Sect. 5). Note that for an ASL forcing strategy u_h(z_sfc) and ϕ(z_sfc) in Eq. (3) would be equal to u_LS(z=10 m) and ϕ_LS(z=10 m) respectively, while in the ABL coupling strategy those variables are provided prognostically by an ABL1d model. As far as the boundary conditions at z=z_top are concerned, Dirichlet boundary conditions u_h(z_top)=u_LS(z_top) and ϕ(z_top)=ϕ_LS(z_top) are prescribed. Model (1) is a first step before evolving toward a more advanced surrogate model including horizontal advection and fine-scale pressure gradients in the future. A particular focus of the present study is on the appropriate choice of a closure scheme to diagnose the eddy diffusivities K_m and K_s. This is a key step to properly represent the downward mixing process.

2.3 Turbulence closure scheme

This subsection describes the turbulence scheme used to compute the eddy diffusivity for momentum and scalars. Those eddy diffusivities are responsible for a vertical mixing of atmospheric variables due to turbulent processes. The turbulence scheme we have implemented in our ABL1d model is very similar to the so-called CBR-1d scheme of Cuxart et al. (2000), which is used operationally at Meteo France (Bazile et al., 2012). We chose to recode the parameterization from scratch for several reasons: computational efficiency, consistency with the NEMO coding rules, use of a geopotential vertical coordinate, and flexibility to add elements specific to the marine atmospheric boundary layer. CBR-1d is a one-equation turbulence closure model based on a prognostic TKE and a diagnostic computation of appropriate length scales. The prognostic equation for the TKE $e=\frac{\mathrm{1}}{\mathrm{2}}\left(〈u^{\prime }{u}^{\prime }〉+〈v^{\prime }{v}^{\prime }〉+〈w^{\prime }{w}^{\prime }〉\right)$ (with 〈⋅〉 the Reynolds averaging operator) is

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\partial }_{t}e=-〈{\mathbf{u}}_{\mathrm{h}}^{\prime }{w}^{\prime }〉\cdot {\partial }_z〈{\mathbf{u}}_{\mathrm{h}}〉+{\displaystyle \frac{g}{{\mathit{\theta }}_{\mathrm{v}}^{\mathrm{ref}}}}〈w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }〉\\ & -{\partial }_z\left(〈e^{\prime }{w}^{\prime }〉+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{a}}}}〈p^{\prime }{w}^{\prime }〉\right)-\mathit{\epsilon },\end{array}\end{array}$$

where horizontal terms and vertical advection are neglected, as usually done in mesoscale atmospheric models. Here θ_v is the virtual potential temperature, ρ_a is atmospheric density, and ε is a dissipation term. In order to express the evolution of e in terms of Reynolds-averaged atmospheric variables we consider the standard closure assumptions for the first order turbulent fluxes (Cuxart et al., 2000) to obtain the classical TKE prognostic equation

$$\begin{array}{}\text{(6)}& {\partial }_{t}e=K_{\mathrm{m}}\Vert {\partial }_z〈{\mathbf{u}}_{\mathrm{h}}〉{\Vert }^{\mathrm{2}}-K_{\mathrm{s}}{N}^{\mathrm{2}}+{\partial }_z\left(K_{\mathrm{e}}{\partial }_{z}e\right)-{\displaystyle \frac{c_{\mathit{\epsilon }}}{l_{\mathit{\epsilon }}}}{e}^{\mathrm{3}/\mathrm{2}},\end{array}$$

where l_ε is a dissipative length scale, c_ϵ a constant, and N² is the moist Brunt–Väisälä frequency computed as $N^{\mathrm{2}}=(g/{\mathit{\theta }}_{\mathrm{v}}^{\mathrm{ref}})\left({\partial }_z〈\mathit{\theta }〉+\mathrm{0.608}{\partial }_z\left(〈\mathit{\theta }〉〈q〉\right)\right)$ with ${\mathit{\theta }}_{\mathrm{v}}^{\mathrm{ref}}=\mathrm{288}\mathrm{K}$. The eddy diffusivities for momentum K_m, TKE K_e, and scalars K_s all depend on e and on a mixing length scale l_m:

$$(K_{\mathrm{m}},K_{\mathrm{s}},K_{\mathrm{e}})=(C_{\mathrm{m}},C_{\mathrm{s}}{\mathit{\varphi }}_z,C_{\mathrm{e}})l_{\mathrm{m}}\sqrt{e},$$

with $(C_{\mathrm{m}},C_{\mathrm{s}},C_{\mathrm{e}})$ a triplet of constants and ϕ_z a stability function proportional to the inverse of a turbulent Prandtl number, given by

$${\displaystyle }{\mathit{\varphi }}_z\left(z\right)={\left(\mathrm{1}+max\left\{C_{\mathrm{1}}{l}_{\mathrm{m}}{l}_{\mathit{\epsilon }}{N}^{\mathrm{2}}/e,-\mathrm{0.5455}\right\}\right)}^{-\mathrm{1}}.$$

The ϕ_z function is bounded not to exceed ${\mathit{\varphi }}_z^{max}=\mathrm{2.2}$ as done in the Arpege model of Meteo France (e.g., Bazile et al., 2012). Assuming that the minimum of ϕ_z is attained in the linearly stratified limit (i.e., for $l_{\mathrm{m}}=l_{\mathit{\epsilon }}=\sqrt{\mathrm{2}e/N^{\mathrm{2}}}$), values of the maximum Prandtl number ${\mathrm{Pr}}_{\mathrm{t}}=C_{\mathrm{m}}/\left(C_{\mathrm{s}}{\mathit{\varphi }}_z\right)$ are given in Table 1. Constant values for C_m, C_s, C_e, c_ε, and C₁ can be determined from different methods, leading to nearly similar values. The traditional way is to use the inertial–convective subrange theory of locally isotropic turbulence (Lilly, 1967; Deardorff, 1974). Another way relies on a theoretical turbulence model partly based on renormalization group methods (see Cheng et al., 2002). For the present study, the sets proposed by Cuxart et al. (2000) and Cheng et al. (2002) will be considered (Table 1). A major difference between the two sets concerns the value of C_m. This difference is explained by a reevaluation of the energy redistribution among velocity components by pressure fluctuations, whose magnitude is assumed to be proportional to the degree of energy anisotropy as initially introduced by Rotta (1951). Note that the constant set of Cheng et al. (2002) is now used by default in both research and operational Meteo France models.

Cuxart et al. (2000)Cheng et al. (2002)

Table 1 Set of turbulence scheme constants from Cuxart et al. (2000) and Cheng et al. (2002). $P{r}_{\mathrm{t}}=C_{\mathrm{m}}/\left(C_{\mathrm{s}}{\mathit{\varphi }}_z\right)$ is the turbulent Prandtl number.

The Dirichlet boundary condition for TKE applied at the top z=z_top is $e(z=z_{\mathrm{top}})=e_{min}={\mathrm{10}}^{-\mathrm{6}}$ m² s⁻² and at the bottom z=z_sfc we have

$$\begin{array}{}\text{(7)}& e(z=z_{\mathrm{sfc}})=e_{\mathrm{sfc}}={\displaystyle \frac{u_{\star }^{\mathrm{2}}}{\sqrt{C_{\mathrm{m}}{c}_{\mathit{\epsilon }}}}}+\mathrm{0.2}{w}_{\star }^{\mathrm{2}},\end{array}$$

with u_⋆ and w_⋆ the friction and convective velocities given by the bulk formulation. The value for e_min has been chosen empirically as well as background values $K_{\mathrm{m}}^{min}={\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹ and $K_{\mathrm{s}}^{min}={\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹ for eddy diffusivities.

The minimum value for l_m is simply set as $l_{\mathrm{m}}^{min}=\frac{K_{\mathrm{m}}^{min}}{C_{\mathrm{m}}\sqrt{e_{min}}}$. There are multiple options to compute the mixing lengths l_m and l_ε (this point will be discussed later in Sect. 3.2.2), but all options have identical boundary conditions $l_{\mathrm{m}}(z=z_{\mathrm{top}})=l_{\mathrm{m}}^{min}$ and

$$\begin{array}{}\text{(8)}& l_{\mathrm{m}}(z=z_{\mathrm{sfc}})=L_{\mathrm{sfc}}=\mathit{\kappa }{\displaystyle \frac{(C_{\mathrm{m}}{c}_{\mathit{\epsilon }}{)}^{\mathrm{1}/\mathrm{4}}}{C_{\mathrm{m}}}}(z_{\mathrm{sfc}}+z_{\mathrm{0}}).\end{array}$$

The value of L_sfc results from the similarity theory in the neutrally stratified surface layer (Sect. 4.1 in Redelsperger et al., 2001, and Appendix A). In Eq. (8), κ is the von Karman constant and z₀ a roughness length computed within the bulk algorithm. The way e_sfc and L_sfc are obtained is detailed in Appendix A.

Our current implementation of boundary layer subgrid processes is an eddy-diffusivity approach which does not include any explicit representation of boundary-layer convective structures. This could be done via a mass-flux representation (e.g., Hourdin et al., 2002; Soares et al., 2004) or the introduction of a countergradient term (e.g., Troen and Mahrt, 1986). This point is left for future developments of the ABL1d model.

2.4 Processing of large-scale forcing and Newtonian relaxation

As mentioned earlier, the ABL1d model (1) requires three-dimensional $(x,y,z)$ large-scale atmospheric variables ${\mathit{\varphi }}_{\mathrm{LS}}^{\mathrm{atm}}$, while existing uncoupled oceanic forcing strategies require only two-dimensional (x,y) atmospheric variables. This is a difficulty for efficiency reasons since it substantially increases the number of I/Os but also for practical reasons because it requires the development of a dedicated tool to extract large-scale atmospheric data and interpolate them on prescribed geopotential heights from their native vertical grid, which can be either pressure based or arbitrary Lagrangian Eulerian. Such tools have been developed specifically to work with ERA-Interim, ERA5, and operational IFS datasets and are described in Appendix B.

Beyond the particular values of ${\mathit{\varphi }}_{\mathrm{LS}}^{\mathrm{atm}}$, the form of the relaxation timescale has a great impact on model solutions. The vertical profile for the λ_m and λ_s coefficients in Eq. (1) is chosen to nudge strongly above the MABL and moderately in the MABL with a smooth transition between its minimum and maximum value to avoid large vertical gradients in λ_m and λ_s, which would result in artificially large vertical gradients in atmospheric variables. In practice the λ_m(z) and λ_s(z) functions depend on the following parameters:

• $({\mathit{\lambda }}_{\mathrm{m}}^{\max },{\mathit{\lambda }}_{\mathrm{m}}^{\min })$ and $({\mathit{\lambda }}_{\mathrm{s}}^{\max },{\mathit{\lambda }}_{\mathrm{s}}^{\min })$, which define the maximum and minimum of the nudging coefficient for momentum and scalars respectively. Following Eq. (2), a guideline to set reasonable values for those parameter values would be to make sure that $\mathrm{\Delta }t{\mathit{\lambda }}_{\mathrm{s}}^{max}\approx \mathrm{1}$ (i.e., the large-scale value is imposed above the boundary layer) and choose ${\mathit{\lambda }}_{\mathrm{s}}^{min}$ based on the typical adjustment timescale of the ABL to surface perturbations. Broadly speaking the ABL can be defined as the region that responds to surface forcings with a timescale of about an hour (e.g., LeMone et al., 2019). In the realistic numerical experiments shown in Sect. 5, we used ${\mathit{\lambda }}_{\mathrm{s}}^{min}=\frac{\mathrm{1}}{\mathrm{90}\left[\min \right]}$, which, for an oceanic dynamical time step Δt=1080 s, would lead to $\mathrm{\Delta }t{\mathit{\lambda }}_{\mathrm{s}}^{min}=\mathrm{0.2}$. (i.e., the boundary layer values are the result of a weighting with a weight 0.8 for the ABL1d prediction and 0.2 for the large-scale value).

• $({\mathit{\beta }}_{min},{\mathit{\beta }}_{max})$, which defines the extent of the transition zone separating the maximum and minimum of the nudging coefficient

We considered the following general form for λ_s(z) and λ_m(z), with h_bl the boundary layer height whose value is diagnosed using an integral Richardson number criteria (Sect. 3.2 and 3.3 in Lemarié et al., 2012) with a critical value equal to C₁:

$$\begin{array}{}\text{(9)}& {\mathit{\lambda }}_{\mathrm{s}}\left(z\right)=\left\{\begin{array}{ll}{\mathit{\lambda }}_{\mathrm{s}}^{\min },& z\le {\mathit{\beta }}_{min}{h}_{\mathrm{bl}},\\ {\sum }_{m=\mathrm{0}}^{\mathrm{3}}{\mathit{\alpha }}_{\mathrm{m}}{\left(\frac{z}{h_{\mathrm{bl}}}\right)}^m,& z\in ]{\mathit{\beta }}_{min}{h}_{\mathrm{bl}};{\mathit{\beta }}_{max}{h}_{\mathrm{bl}}[,\\ {\mathit{\lambda }}_{\mathrm{s}}^{\max },& z\ge {\mathit{\beta }}_{max}{h}_{\mathrm{bl}},\end{array}\right.\end{array}$$

where four α_m coefficients are necessary to guarantee the continuity of λ_s(z) and its derivative ∂_zλ_s at $z={\mathit{\beta }}_{min}{h}_{\mathrm{bl}}$ and $z={\mathit{\beta }}_{max}{h}_{\mathrm{bl}}$. We easily find

$${\displaystyle }\begin{array}{l}{\mathit{\alpha }}_{\mathrm{0}}=\frac{(\mathrm{3}{\mathit{\beta }}_{max}-{\mathit{\beta }}_{min}){\mathit{\beta }}_{min}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{s}}^{\max }+({\mathit{\beta }}_{max}-\mathrm{3}{\mathit{\beta }}_{min}){\mathit{\beta }}_{max}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{s}}^{\min }}{({\mathit{\beta }}_{max}-{\mathit{\beta }}_{min}{)}^{\mathrm{3}}},\\ {\mathit{\alpha }}_{\mathrm{1}}=-\frac{\mathrm{6}{\mathit{\beta }}_{max}{\mathit{\beta }}_{min}({\mathit{\lambda }}_{\mathrm{s}}^{\max }-{\mathit{\lambda }}_{\mathrm{s}}^{\min })}{({\mathit{\beta }}_{max}-{\mathit{\beta }}_{min}{)}^{\mathrm{3}}},\\ {\mathit{\alpha }}_{\mathrm{2}}=\mathrm{3}\frac{({\mathit{\beta }}_{max}+{\mathit{\beta }}_{min})({\mathit{\lambda }}_{\mathrm{s}}^{\max }-{\mathit{\lambda }}_{\mathrm{s}}^{\min })}{({\mathit{\beta }}_{max}-{\mathit{\beta }}_{min}{)}^{\mathrm{3}}},\\ {\mathit{\alpha }}_{\mathrm{3}}=-\frac{\mathrm{2}({\mathit{\lambda }}_{\mathrm{s}}^{\max }-{\mathit{\lambda }}_{\mathrm{s}}^{\min })}{({\mathit{\beta }}_{max}-{\mathit{\beta }}_{min}{)}^{\mathrm{3}}}.\end{array}$$

The value of h_bl is bounded beforehand to guarantee that at least 3 grid points are such that $z\le {\mathit{\beta }}_{min}{h}_{\mathrm{bl}}$ and $z\ge {\mathit{\beta }}_{max}{h}_{\mathrm{bl}}$. A typical profile of the λ_s(z) is shown in Fig. 1a.

Figure 1 (a) Typical profile of the nudging coefficient λ_s(z) with respect to the parameters ${\mathit{\lambda }}_{\mathrm{s}}^{\max },{\mathit{\lambda }}_{\mathrm{s}}^{\min },{\mathit{\beta }}_{min},{\mathit{\beta }}_{max},h_{\mathrm{bl}}$. (b) Equatorial restoring function r_eq with respect to the Coriolis frequency f.

When the model is forced by the large-scale pressure gradient (or the geostrophic winds), the parameter λ_m(z) should be theoretically zero at high and middle latitudes. However, for the equatorial region, a Newtonian relaxation toward the large-scale winds should be maintained. To do so, the coefficient λ_m(z) is multiplied by a coefficient r_eq, which is a function of the Coriolis parameter f. The r_eq coefficient is equal to zero for large values of $\left|f\right|$ and increases to 1 when approaching the Equator. The following form satisfies those constraints (see also Fig. 1b):

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& r_{\mathrm{eq}}\left(f\right)=\sin {\left({\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\left[{\displaystyle \frac{f-f_{max}}{f_{max}}}\right]\right)}^{\mathrm{6}},\\ & f_{max}={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{\mathrm{12}\times \mathrm{3600}}}{\mathrm{s}}^{-\mathrm{1}}.\end{array}\end{array}$$

3 Numerical discretization and implementation within NEMO

We have introduced so far the continuous formulation of the ABL1d model. In this section we describe the discretization methods used and how this model is included in the NEMO modeling framework. In particular, the discretization of the Coriolis term and of the TKE Eq. (6) and associated mixing lengths are described in Sect. 3.1 and 3.2 respectively. Details about the practical implementation in NEMO are given in Sect. 3.3 for the coupling aspects and 3.4 for the computational aspects. The ABL1d model (1) is discretized in time with an Euler backward scheme for the vertical diffusion terms, semi-implicitly for the Coriolis term, and explicitly for the relaxation term, which means that the model is stable as long as λ_sΔt≤1. The variables are defined on a non-staggered grid in the horizontal (a.k.a Arakawa A-grid). Because we consider a computational domain exclusively over water or sea ice, topography is not considered and vertical levels are flat and fixed in time which, among other things, allows the large-scale data ϕ_LS to be interpolated on the vertical grid offline. The position of the various quantities introduced so far on the computational grid is given in Fig. 2.

Figure 2 Vertical grid variable arrangements and important notations.

3.1 Coriolis term treatment

Since in our implementation the horizontal velocity components are collocated, the discretization of the Coriolis term is straightforward and is energetically neutral. In the event the ABL1d is integrated with a time step much larger than the oceanic time step, specific care must be given to the stability of the Coriolis term time stepping. A semi-implicit scheme with weighting parameter γ reads ${\mathbf{u}}_{\mathrm{h}}^{n+\mathrm{1},\star }=-\left(f\mathrm{\Delta }t\right)\mathit{k}\times \left((\mathrm{1}-\mathit{\gamma }){\mathbf{u}}_{\mathrm{h}}^n+\mathit{\gamma }{\mathbf{u}}_{\mathrm{h}}^{n+\mathrm{1},\star }\right)$, where the exponent ⋆ is used here to emphasize that ${\mathbf{u}}_{\mathrm{h}}^{n+\mathrm{1},\star }$ is a temporary value at time n+1 before vertical diffusion and Newtonian relaxation are applied. For a given grid cell with index (i,j), the semi-implicit scheme can be written in a more compact way as

$$\begin{array}{l}{\displaystyle }{u}_{i,j}^{n+\mathrm{1},\star }={\displaystyle \frac{(\mathrm{1}-\mathit{\gamma }(\mathrm{1}-\mathit{\gamma }\left)\right(f_{i,j}\mathrm{\Delta }t{)}^{\mathrm{2}})u_{i,j}^n+(f_{i,j}\mathrm{\Delta }t)v_{i,j}^n}{\mathrm{1}+(f_{i,j}\mathrm{\Delta }t{)}^{\mathrm{2}}{\mathit{\gamma }}^{\mathrm{2}}}},\\ {\displaystyle }{v}_{i,j}^{n+\mathrm{1},\star }={\displaystyle \frac{(\mathrm{1}-\mathit{\gamma }(\mathrm{1}-\mathit{\gamma }\left)\right(f_{i,j}\mathrm{\Delta }t{)}^{\mathrm{2}})v_{i,j}^n-(f_{i,j}\mathrm{\Delta }t)u_{i,j}^n}{\mathrm{1}+(f_{i,j}\mathrm{\Delta }t{)}^{\mathrm{2}}{\mathit{\gamma }}^{\mathrm{2}}}}.\end{array}$$

The associated amplification factor modulus is $\left|{\mathcal{A}}_{\mathrm{cor}}\right|=\sqrt{\frac{\mathrm{1}+(\mathrm{1}-\mathit{\gamma }{)}^{\mathrm{2}}(f\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{1}+{\mathit{\gamma }}^{\mathrm{2}}(f\mathrm{\Delta }t{)}^{\mathrm{2}}}}$ meaning that unconditional stability is obtained as long as $\mathit{\gamma }\ge \mathrm{1}/\mathrm{2}$. For the numerical results obtained below in Sects. 4 and 5 we used γ=0.55, which is deliberately slightly dissipative.

3.2 Discretization of TKE equation

In Sect. 2.3 we have presented the continuous formulation of the TKE-based turbulence closure of the ABL1d model. In the following we describe how the positivity of TKE can be preserved and how the mixing lengths l_m and l_ε are computed. We provide a substantial discussion on the latter aspect because numerical results are very sensitive to the choices made.

3.2.1 TKE positivity preservation

The TKE equation is discretized using a backward Euler scheme in time with a linearization of the dissipation term $\frac{c_{\mathit{\epsilon }}}{l_{\mathit{\epsilon }}}{e}^{\mathrm{3}/\mathrm{2}}$, which is discretized as $\frac{c_{\mathit{\epsilon }}}{l_{\mathit{\epsilon }}}\sqrt{e^n}{e}^{n+\mathrm{1}}$. However, such discretization is not unconditionally positivity-preserving for TKE, which could give rise to unphysical solutions (e.g., Burchard, 2002b). Ignoring the diffusion term, the TKE prognostic Eq. (6) can be written as an ordinary differential equation (ODE) of the form

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& {\partial }_{t}e=S({\mathbf{u}}_{\mathrm{h}},N^{\mathrm{2}})-D(e,t)e,\text{with}S({\mathbf{u}}_{\mathrm{h}},N^{\mathrm{2}})=\\ & K_{\mathrm{m}}\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}{\Vert }^{\mathrm{2}}-K_t{N}^{\mathrm{2}},D(e,t)={\displaystyle \frac{c_{\mathit{\epsilon }}}{l_{\mathit{\epsilon }}}}\sqrt{e^n},\end{array}\end{array}$$

where the last term can be seen as a damping term. For ODEs like Eq. (11) it can be shown that for an initial condition e(0)≥0 and $S({\mathbf{u}}_{\mathrm{h}},N^{\mathrm{2}})\ge \mathrm{0}$, the solution e(t) keeps the same sign as e(0) whatever the sign of the damping coefficient D(e,t). Assuming that S(u_h,N²) and D(e,t) are positive, a backward Euler discretization of the damping term in Eq. (11) would lead to $e^{n+\mathrm{1}}=\frac{e^n+\mathrm{\Delta }tS({\mathbf{u}}_{\mathrm{h}},N^{\mathrm{2}})}{\mathrm{1}+\mathrm{\Delta }tD(e,t)}$, which preserves positivity since for eⁿ≥0 we obtain $e^{n+\mathrm{1}}\ge \mathrm{0}$. However, there is no guarantee that the forcing term S(u_h,N²) is positive, in particular when the shear is weak and the stratification is large. When S(u_h,N²) is negative a specific treatment (known as the “Patankar trick”; see Deleersnijder et al., 1997; Burchard, 2002b) is required. In the event of a negative S(u_h,N²), the idea is to move the buoyancy term from S to D after dividing it by eⁿ, such that $S({\mathbf{u}}_{\mathrm{h}},N^{\mathrm{2}})=K_{\mathrm{m}}\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}{\Vert }^{\mathrm{2}}$ is now strictly positive and $D(e,t)=\frac{c_{\mathit{\epsilon }}}{l_{\mathit{\epsilon }}}\sqrt{e^n}+K_{\mathrm{s}}\frac{N^{\mathrm{2}}}{e^n}$. Such a procedure is a sufficient condition to preserve the positivity of the TKE without ad hoc clipping of negative values. Moreover our discretization of the shear and buoyancy terms in the TKE equation is done in an energetically consistent way following Burchard (2002a).

3.2.2 Mixing length computation

Another challenging task when implementing a TKE scheme is the discretization of the mixing lengths. As mentioned earlier, four different discretizations of l_m (l_ε) have been coded. All discretizations consider the boundary conditions given in Eq. (8). The values of l_m and l_ε are traditionally computed from two intermediate length scales l_up and l_dwn, which respectively correspond to the maximum upward and downward displacement of a parcel of air with a given initial kinetic energy. Once l_up and l_dwn have been estimated by one of the methods described below, the dissipative and mixing length scales l_m and l_ε are computed as

$$\begin{array}{}\text{(12a)}& {\displaystyle }{l}_{\mathrm{m}}& {\displaystyle }={\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left\{l_{\mathrm{up}}^{\frac{\mathrm{1}}{a}}+l_{\mathrm{dwn}}^{\frac{\mathrm{1}}{a}}\right\}\right)}^a,\text{(12b)}& {\displaystyle }{l}_{\mathit{\epsilon }}& {\displaystyle }=min\left(l_{\mathrm{up}},l_{\mathrm{dwn}}\right),\end{array}$$

where $a\approx -\frac{\mathrm{3}}{\mathrm{2}}$ for CBR00 and $a\approx -\frac{\mathrm{6}}{\mathrm{7}}$ for CCH02 (see Appendix A). The impact of the weighting between l_up and l_dwn to compute l_m can be significant for idealized experiments like the ones presented in Sect. 4.2 but for more realistic cases results are weakly sensitive and equivalent to the ones obtained with the simpler weighting $l_{\mathrm{m}}=\sqrt{l_{\mathrm{up}}{l}_{\mathrm{dwn}}}$.

In the following we provide the continuous form of the various ways to compute l_up and l_dwn implemented in the ABL1d model. The discretization aspects are detailed in Appendix C.

1. Bougeault and Lacarrère (1989) length scale. A classical approach in atmospheric models is the use of the Bougeault and Lacarrère (1989) mixing length (see also Bougeault and André, 1986) which defines l_up and l_dwn as

   $$\begin{array}{}\text{(13)}& \begin{array}{rl}& \underset{z}{\overset{z+l_{\mathrm{up}}}{\int }}{N}^{\mathrm{2}}\left(s\right)(s-z)\mathrm{d}s=e\left(z\right),\\ & \underset{z-l_{\mathrm{dwn}}}{\overset{z}{\int }}{N}^{\mathrm{2}}\left(s\right)(z-s)\mathrm{d}s=e\left(z\right).\end{array}\end{array}$$

   By construction such mixing lengths are bounded by the distance to the bottom and the top of the computational domain. It is worth noting that for constant values of N², Eq. (13) gives $\frac{l_{\mathrm{up}}^{\mathrm{2}}{N}^{\mathrm{2}}}{\mathrm{2}}=e\left(z\right)$ and $\frac{l_{\mathrm{dwn}}^{\mathrm{2}}{N}^{\mathrm{2}}}{\mathrm{2}}=e\left(z\right)$ respectively, which is equivalent to the Deardorff (1980) length scale. In the remainder we will note l_BL89, the mixing length obtained from Eq. (13).

2. Adaptation of NEMO's length scale. The standard NEMO algorithm (Sect. 10.1.3 in Madec, 2012) is simple and efficient compared to Eq. (13). This algorithm is based on the Deardorff (1980) length scale $l_{\mathrm{D}\mathrm{80}}=\sqrt{\mathrm{2}e\left(z\right)/N^{\mathrm{2}}}$. l_up and l_dwn are first initialized to $l_{\mathrm{up}}=l_{\mathrm{dwn}}=l_{\mathrm{D}\mathrm{80}}$. The resulting length scales are then limited not only by the distance to the surface and to the top but also by the distance to a strongly stratified portion of the air column. This limitation amounts to control of the vertical gradients of l_up(z) and l_dwn(z) such that they are not larger than the variations of altitude. The resulting mixing length will be simply referred to as l_D80. Note that the Taylor expansion of the integral in Eq. (13) is

   $${\displaystyle }\underset{z}{\overset{z+l_{\mathrm{up}}}{\int }}{N}^{\mathrm{2}}\left(s\right)(s-z)\mathrm{d}s\approx {\displaystyle \frac{N^{\mathrm{2}}\left(z\right)l_{\mathrm{up}}^{\mathrm{2}}}{\mathrm{2}}}+{\displaystyle \frac{\frac{d{N}^{\mathrm{2}}}{\mathrm{d}z}{l}_{\mathrm{up}}^{\mathrm{3}}}{\mathrm{3}}}+\mathcal{O}\left(l_{\mathrm{up}}^{\mathrm{4}}\right),$$

   which shows that the l_D80 mixing length is an approximation of l_BL89, which is obtained by retaining only the leading order term in the Taylor expansion.

3. Rodier et al. (2017) length scale. Recently, Rodier et al. (2017) proposed a modification of the Bougeault and Lacarrère (1989) mixing length. This modification turns out to improve results for stably stratified boundary layers typical of areas covered by ice. They propose to add a shear-related term to Eq. (13) such that the definition of l_up and l_dwn becomes

   $$\begin{array}{}\text{(14)}& \begin{array}{rl}& \underset{z}{\overset{z+l_{\mathrm{up}}}{\int }}\left[N^{\mathrm{2}}\left(s\right)(s-z)+c_{\mathrm{0}}\sqrt{e\left(s\right)}\Vert {\partial }_{\mathrm{s}}{\mathbf{u}}_{\mathrm{h}}\Vert \right]\mathrm{d}s=e\left(z\right),\\ & \underset{z-l_{\mathrm{dwn}}}{\overset{z}{\int }}\left[N^{\mathrm{2}}\left(s\right)(z-s)+c_{\mathrm{0}}\sqrt{e\left(s\right)}\Vert {\partial }_{\mathrm{s}}{\mathbf{u}}_{\mathrm{h}}\Vert \right]\mathrm{d}s=e\left(z\right),\end{array}\end{array}$$

   where c₀ is a parameter whose value should be smaller than $\sqrt{C_{\mathrm{m}}/c_{\mathit{\epsilon }}}$. The value of c₀ will be chosen based on numerical experiments presented in Sect. 4. In the following this mixing length will be referred to as l_R17.

4. A local buoyancy- and shear-based length scale. For the sake of computational efficiency, we have derived a local version of the Rodier et al. (2017) length scale which is original to the present paper. Under the assumption that l_up (l_dwn) is small compared to the spatial variations of N², e, and $\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}\Vert$, we end up with the following second-order equation for l_up:

   $${\displaystyle \frac{N^{\mathrm{2}}\left(z\right)}{\mathrm{2}}}{l}_{\mathrm{up}}^{\mathrm{2}}+c_{\mathrm{0}}\sqrt{e\left(z\right)}\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}\Vert l_{\mathrm{up}}=e\left(z\right),$$

   whose unique positive solution is

   $${\displaystyle }{l}_{\mathrm{D}\mathrm{80}}^{\star }\left(z\right)={\displaystyle \frac{\mathrm{2}\sqrt{e\left(z\right)}}{c_{\mathrm{0}}\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}\Vert +\sqrt{c_{\mathrm{0}}^{\mathrm{2}}\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}{\Vert }^{\mathrm{2}}+\mathrm{2}{N}^{\mathrm{2}}\left(z\right)}}}.$$

   We easily find that $l_{\mathrm{D}\mathrm{80}}^{\star }=l_{\mathrm{D}\mathrm{80}}$ for $\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}\Vert =\mathrm{0}$, and $l_{\mathrm{D}\mathrm{80}}^{\star }=\frac{\sqrt{e\left(z\right)}}{c_{\mathrm{0}}\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}\Vert }$ for N²=0, which is consistent with the shear-based length scale of Wilson and Venayagamoorthy (2015). Once $l_{\mathrm{D}\mathrm{80}}^{\star }$ has been computed we apply the same algorithmic approach as in the l_D80 case.

The performance of those four length scales for various physical flows is discussed in Sect. 4.

3.3 Coupling with ocean and sea ice

For the practical implementation of the ABL coupling strategy within a global oceanic model, a proper coupling method is required for stability and consistency purposes (e.g., Beljaars et al., 2017; Renault et al., 2019a), and the ABL1d must have the ability to handle grid cells partially covered by sea ice. For the coupling strategy, a so-called implicit flux coupling which is unconditionally stable (Appendix B in Beljaars et al., 2017) and asymptotically consistent for Δt→0 (Renault et al., 2019a) is used. Because vertical diffusion in ABL1d is handled implicitly in time, the boundary conditions (Eqs. 3 and 4) should be provided at time n+1. The implicit flux coupling amounts to discretize the boundary conditions Eqs. (3) and (4) as

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }{\left.K_{\mathrm{m}}{\partial }_z{\mathbf{u}}_{\mathrm{h}}\right|}_{z=z_{\mathrm{sfc}}}^{n+\mathrm{1}}=C_D\Vert {\mathbf{u}}_{\mathrm{h}}^n\left(z_{\mathrm{sfc}}\right)-{\stackrel{\mathrm{̃}}{\mathbf{u}}}_{\mathrm{oce}}\Vert \left({\mathbf{u}}_{\mathrm{h}}^{n+\mathrm{1}}\right(z_{\mathrm{sfc}})-{\stackrel{\mathrm{̃}}{\mathbf{u}}}_{\mathrm{oce}}),\text{(16)}& {\displaystyle }& {\displaystyle }{\left.K_{\mathrm{s}}{\partial }_z\mathit{\varphi }\right|}_{z=z_{\mathrm{sfc}}}^{n+\mathrm{1}}=C_{\mathit{\varphi }}\Vert {\mathbf{u}}_{\mathrm{h}}^n\left(z_{\mathrm{sfc}}\right)-{\stackrel{\mathrm{̃}}{\mathbf{u}}}_{\mathrm{oce}}\Vert \left({\mathit{\varphi }}^{n+\mathrm{1}}\right(z_{\mathrm{sfc}})-{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{oce}}),\end{array}$$

where ${\stackrel{\mathrm{̃}}{\mathbf{u}}}_{\mathrm{oce}}$ and ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{oce}}$ are either the instantaneous values at time n if NEMO and ABL1d have the same time step or an average over the successive oceanic substeps otherwise.

Particular care has also been given to the compatibility between the ABL1d model and SI³ (Sea Ice model Integrated Initiative) the sea-ice component of NEMO. SI³ is a multi-category model whose state variables relevant for our study are the ice surface temperature $T_l^{\mathrm{ice}}$ with associated fractional area a_l (for the lth category), and the ice velocity u^ice (same for all categories). Note that the values of the exchange coefficients over sea ice $C_D^{\mathrm{ice}}$, $C_{\mathit{\theta }}^{\mathrm{ice}}$, and $C_q^{\mathrm{ice}}$ are different from their oceanic counterparts but are the same over all sea-ice categories. At this point there are several strategies for the ABL1d∕SI³ coupling:

1. Run the ABL1d model over the whole ABL for each category l and then average atmospheric variables weighted by a_l.

2. Run a single ABL1d model with a category-averaged surface flux. In the current version of NEMO $C_{\mathit{\theta }}^{\mathrm{ice}}$ is a function of the averaged temperature T^ice which means that it is equivalent to compute a flux over each category before averaging them and to compute a single flux using the averaged surface temperature, indeed

   $$\begin{array}{ll}{\displaystyle }& {\displaystyle }\sum _l{a}_l\left[C_{\mathit{\theta }}^{\mathrm{ice}}\Vert {\mathbf{u}}_{\mathrm{h}}\left(z_{\mathrm{sfc}}\right)-{\mathbf{u}}^{\mathrm{ice}}\Vert \left(\mathit{\theta }\right(z_{\mathrm{sfc}})-T_l^{\mathrm{ice}})\right]=\\ {\displaystyle }& {\displaystyle }{C}_{\mathit{\theta }}^{\mathrm{ice}}\Vert {\mathbf{u}}_{\mathrm{h}}\left(z_{\mathrm{sfc}}\right)-{\mathbf{u}}^{\mathrm{ice}}\Vert \left(\mathit{\theta }\left(z_{\mathrm{sfc}}\right)-\sum _l{a}_l{T}_l^{\mathrm{ice}}\right).\end{array}$$

The second option has been preferred because it is much easier to implement and more computationally efficient. It amounts to consider an ice surface temperature averaged over all categories $T^{\mathrm{ice}}={\sum }_{l=\mathrm{1}}^{n_{\mathrm{cat}}}{a}_l{T}_l^{\mathrm{ice}}$ for the computation of ice–atmosphere turbulent fluxes (T^ice also enters in the computation of q_ice). Noting F_oce the fraction of open water (lead), the boundary condition (15) and (16) are modified in

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\left.K_{\mathrm{m}}{\partial }_z{\mathbf{u}}_{\mathrm{h}}\right|}_{z=z_{\mathrm{sfc}}}^{n+\mathrm{1}}=F_{\mathrm{oce}}{C}_D\Vert {\mathbf{u}}_{\mathrm{h}}^n\left(z_{\mathrm{sfc}}\right)-{\stackrel{\mathrm{̃}}{\mathbf{u}}}_{\mathrm{oce}}\Vert \left({\mathbf{u}}_{\mathrm{h}}^{n+\mathrm{1}}\right(z_{\mathrm{sfc}})-{\stackrel{\mathrm{̃}}{\mathbf{u}}}_{\mathrm{oce}})\\ {\displaystyle }& {\displaystyle }+(\mathrm{1}-F_{\mathrm{oce}})C_D^{\mathrm{ice}}\Vert {\mathbf{u}}_{\mathrm{h}}^n\left(z_{\mathrm{sfc}}\right)-{\stackrel{\mathrm{̃}}{\mathbf{u}}}^{\mathrm{ice}}\Vert \left({\mathbf{u}}_{\mathrm{h}}^{n+\mathrm{1}}\right(z_{\mathrm{sfc}})-{\stackrel{\mathrm{̃}}{\mathbf{u}}}^{\mathrm{ice}}),\\ {\displaystyle }& {\displaystyle }{\left.K_{\mathrm{s}}{\partial }_z\mathit{\varphi }\right|}_{z=z_{\mathrm{sfc}}}^{n+\mathrm{1}}=F_{\mathrm{oce}}{C}_{\mathit{\varphi }}\Vert {\mathbf{u}}_{\mathrm{h}}^n\left(z_{\mathrm{sfc}}\right)-{\stackrel{\mathrm{̃}}{\mathbf{u}}}_{\mathrm{oce}}\Vert \left({\mathit{\varphi }}^{n+\mathrm{1}}\right(z_{\mathrm{sfc}})-{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{oce}})\\ {\displaystyle }& {\displaystyle }+(\mathrm{1}-F_{\mathrm{oce}})C_{\mathit{\varphi }}^{\mathrm{ice}}\Vert {\mathbf{u}}_{\mathrm{h}}^n\left(z_{\mathrm{sfc}}\right)-{\stackrel{\mathrm{̃}}{\mathbf{u}}}^{\mathrm{ice}}\Vert \left({\mathit{\varphi }}^{n+\mathrm{1}}\right(z_{\mathrm{sfc}})-{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}^{\mathrm{ice}}).\end{array}$$

Because the dynamics of sea ice is computed before the thermodynamics (see Fig. 1 in Rousset et al., 2015), the ABL1d–SI³ coupling follows these different steps:

1. compute surface fluxes over ice and ocean and integrate the ABL1d model for given values $F_{\mathrm{oce}}^n$ and $a_l^n$,

2. compute the dynamics of sea ice,

3. update $F_{\mathrm{oce}}^n$ and $a_l^n$ in $F_{\mathrm{oce}}^{\star }$ and $a_l^{\star }$ because of step 2,

4. distribute the fluxes over each ice category considering the updated values $a_l^{\star }$ (Sect. 3.6 in Rousset et al., 2015)

5. compute the thermo-dynamics of sea ice.

3.4 Computational aspects

As described in Maisonnave and Masson (2015), the NEMO source code is organized to separate the ocean routines on one side and the routines responsible for the surface boundary conditions computation (including sea ice and the coupling interfaces) on the other side. This makes a clear separation between the standard ocean model (OCE component) and the so-called surface module (SAS component). As schematically described in Fig. 3, the ABL1d model has been implemented within the SAS component, which allows the following useful features:

• The ABL1d model can be run in standalone mode (coupled or not with sea ice) with prescribed oceanic surface fields.

• The ABL1d model can be run in detached mode; i.e., the OCE and SAS components run on potentially separate processors and computational grids communicating via the OASIS3-MCT coupling library (Craig et al., 2017).

Figure 3 Schematic representation of the ASL forcing strategy (left) and ABL coupling strategy (right) in terms of code organization and required external data. The OCE and SI³ components represent the oceanic and sea-ice dynamics and thermodynamics respectively while the ASL component is in charge of providing boundary conditions related to atmospheric conditions. In the NEMO computational framework the so-called surface module (SAS component), delineated by dashed line polygons, is virtually separated from the OCE component, which allows SAS to be run in standalone or detached mode (see Sect. 3.4).

An other capability implemented within the NEMO modeling framework is the possibility to interpolate forcing fields on the fly. This is particularly useful for the ABL coupling strategy since three-dimensional atmospheric data must be interpolated on the ABL1d computational grid. As the current implementation of the on-the-fly interpolation only works in the horizontal, the vertical interpolation of large-scale atmospheric data on the ABL1d vertical grid is done offline. Nevertheless it means that the size of input data compared to an ASL forcing strategy is N times larger with N the number of vertical levels in the ABL. A possibility to improve the efficiency for the reading of input data would be to take advantage of the parallel I/O capabilities provided by the XIOS library (XML-IO-Server; Meurdesoif et al., 2016) which is currently used in NEMO only for writing output data. This technical development is left for future work. This is a key aspect because, as discussed later in Appendix D, the main source of computational overhead associated with the ABL coupling strategy is due to the time spent waiting for input files to be read.

Table 2 Description of the idealized experiments performed in Sect. 4. L_MO is the Monin–Obukhov length.

4 Atmosphere-only numerical experiments

4.1 Sensitivity experiments and objectives

To check the relevance of our ABL1d model for idealized atmospheric situations typical of the atmospheric boundary layer over water or sea ice, we performed a set of single-column experiments. Each of those experiments are evaluated with benchmark large eddy simulations (LESs). Moreover, we use standardized test cases from the literature to allow our results to be cross-compared with other well-established ABL schemes. In the following we consider a neutrally stratified (Sect. 4.2) and a stably stratified (Sect. 4.3) case as well as a case with a transition from stable to unstable stratification representative of an atmospheric flow over an SST front (Sect. 4.4). All ABL1d simulations presented here have been performed directly within the SAS component of the NEMO modeling framework and can be reproduced using the code available at https://zenodo.org/record/3904518 (last access: 20 January 2021) (Lemarié and Samson, 2020), which also includes the scripts to generate the figures. An objective of the present section is to illustrate the type of sensitivity we can expect from the ABL1d model and discriminate between the various options available in the code. The experiments showed in Sect. 4.2 and 4.3 are meant to investigate the impact of (i) the set of constant coefficients (CBR00 vs. CCH02), (ii) the various formulations of l_m and l_ε among the algorithms described in Sect. 3.2.2, and (iii) the parameter value c₀ in the l_R17 and $l_{\mathrm{D}\mathrm{80}}^{\star }$ mixing length computation. Those experiments will allow several options to be discarded. The ability of the remaining options to represent the downward mixing mechanism, discussed in Sect. 1.2, is then evaluated in Sect. 4.4. The robustness of the results to the bulk formulation and to the nudging coefficient is also checked. For each experiment we explicitly provide the initial and boundary conditions as well as all the necessary parameter values (see Table 2) so that the experiments can be reproduced easily by other modeling groups.

Figure 4 Results obtained for the neutral boundary-layer case of Andren et al. (1994) with the CBR00 model constants (left panels, a to d) and CCH02 model constants (right panels, d to e) for different parameter values for c₀ and different mixing-length formulations ($l_{\mathrm{D}\mathrm{80}}^{\star }$ for black lines or l_R17 for gray lines). Results are shown for u (a, e), v (b, f), e (c, g), and l_m (d, h). In the top four panels results are compared with LESs from Cuxart et al. (2000) (their Fig. 16). As in Andren et al. (1994), simulations were run over a period of 10∕f and results are averaged over the last 3∕f period.

4.2 Neutral turbulent Ekman layer

We first propose to investigate the simulation of a neutrally stratified atmosphere analogous to a classical turbulent Ekman layer. The selected case is based on the setup described in Andren et al. (1994). The initial conditions for this experiment are not defined analytically; they are given by Table A1 in Andren et al. (1994)². This test case is mainly used to check the adequacy of our surface boundary conditions with similarity theory and the proper calibration of the parameter c₀ in the $l_{\mathrm{D}\mathrm{80}}^{\star }$ and l_R17 formulations of the mixing lengths. In theory, the l_D80 and l_BL89 mixing lengths do not support the asymptotic limit N²=0 but for the integrity of numerical results a minimum threshold $N_{\mathit{\epsilon }}^{\mathrm{2}}$ on the stratification is imposed in the code. In this case the procedure to compute those mixing lengths as described in Appendix C will provide identical results, namely $l_{\mathrm{up}}=z_{\mathrm{top}}-z$ and $l_{\mathrm{dwn}}=z-z_{\mathrm{sfc}}$ (i.e., the distance from the top and from the bottom of the computational domain). We test here the $l_{\mathrm{D}\mathrm{80}}^{\star }$ and l_R17 introduced in Sect. 3.2.2. The reference solution is taken from Cuxart et al. (2000) (panels a and b in their Fig. 16). Results are obtained using the ABL1d model with either the CBR00 (Fig. 4a–d) or the CCH02 (Fig. 4e–h) set of parameters. All experiments have been done with c₀=0.15 and c₀=0.2. All simulations are able to reproduce the overall behavior of the LES case. The main outcomes are as follows:

• The best agreement is obtained when using the CCH02 constants along with $l_{\mathrm{D}\mathrm{80}}^{\star }$ mixing length and c₀=0.2.

• The results obtained for l_D80 and l_BL89 are identical and close to the l_R17 results with c₀=0.15 (not shown).

• All simulations with the CCH02 set of parameters show reasonable results.

Figure 5 Results obtained for the stably stratified boundary-layer case of Cuxart et al. (2006) for the parameter values CBR00 (left panels, a to c) and CCH02 (right panels, d to f) with different mixing-length formulations: l_D80 for black solid lines, $l_{\mathrm{D}\mathrm{80}}^{\star }$ with c₀=0.15 for dashed black lines (c₀=0.2 for dotted black lines), l_BL89 for solid gray lines, and l_R17 with c₀=0.15 for dashed gray lines (c₀=0.2 for dotted gray lines). Results are shown for potential temperature θ (a, d), wind speed  (b, e), and l_m (c, f). Dotted red lines represent LES results from Rodier et al. (2017). Instantaneous profiles after 9 h are shown.

4.3 Stably stratified boundary layer (GABLS1)

Within the Global Energy and Water Exchanges (GEWEX) atmospheric boundary layer study (GABLS), idealized cases for stable surface boundary layers have been investigated (e.g., Cuxart et al., 2006). Such conditions are typical of areas covered with sea ice. Here we consider the GABLS1 case, whose technical description is available at http://turbulencia.uib.es/gabls/gabls1d_desc.pdf (last access: 20 January 2021). This experiment is particularly interesting as significant differences generally exist between solutions obtained from LES and single-column simulations, for example when the Bougeault and Lacarrère (1989) length scale is used (e.g., Cuxart et al., 2006; Rodier et al., 2017). A large-scale geostrophic wind is imposed as well as a cooling of the surface temperature θ_s(t) given by ${\mathit{\theta }}_{\mathrm{s}}\left(t\right)=\mathrm{263.5}-\mathrm{0.25}(t/\mathrm{3600}\mathrm{s})$. The parameter values for this test are reported in Table 2 and the initial conditions are ${\mathbf{u}}_{\mathrm{h}}(z,t=\mathrm{0})={\mathbf{u}}_{\mathrm{G}}$, and

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\theta }(z,t=\mathrm{0})=\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}\mathrm{265}& z\le \mathrm{100}\mathrm{m}\\ \mathrm{265}+\mathrm{0.01}(z-\mathrm{100})& \mathrm{otherwise}\end{array}\right.,\\ {\displaystyle }& {\displaystyle }e(z,t=\mathrm{0})=\left\{\begin{array}{ll}{e}_{min}+\mathrm{0.4}(\mathrm{1}-z/\mathrm{250}{)}^{\mathrm{3}}& z\le \mathrm{250}\mathrm{m}\\ e_{min}& \mathrm{otherwise}\end{array}.\right.\end{array}$$

The solutions after 9 h of simulation are shown in Fig. 5 (left panels) for CBR00 parameter values and in Fig. 5 (right panels) for CCH02 parameter values. The reference solution is taken from Rodier et al. (2017) LESs. As expected, solutions based on a mixing length ignoring the contribution from the vertical shear exhibit a boundary layer that is too thick and a wind speed maximum located too high in altitude. Using a buoyancy- and shear-based mixing length mitigates the issue and provides very good agreement with reference solutions when the CCH02 model constants are used. The best results are obtained for $l_{\mathrm{D}\mathrm{80}}^{\star }$ with c₀=0.2 and l_R17 with c₀=0.15. Solutions obtained with the CBR00 model constants systematically predict larger turbulent kinetic energy and mixing lengths, resulting in large values of K_s in the first 100 m near the surface (not shown). The mismatch in terms of TKE is partially explained by the difference in boundary conditions since with CBR00 constants we have $e_{\mathrm{sfc}}=\mathrm{4.628}{u}_{\star }^{\mathrm{2}}$ while with CCH02 constants we get $e_{\mathrm{sfc}}=\mathrm{3.065}{u}_{\star }^{\mathrm{2}}$ from Eq. (7). Note that the proper calibration of the c₀ constant jointly with the c_ε is the subject of several ongoing studies. Since our simulations reproduce the known sensitivity to those parameters, the ABL1d model could directly benefit from new findings on that topic.

The main outcomes are as follows:

• The CCH02 set of parameters provides results of better quality than the CBR00 constants. For the sake of simplicity, we will retain only the CCH02 parameters for the numerical results shown in the remainder.

• The buoyancy- and vertical-shear-based mixing lengths l_R17 and $l_{\mathrm{D}\mathrm{80}}^{\star }$ are superior to the buoyancy-based mixing lengths l_D80 and l_BL89 for stable boundary layers.

4.4 Winds across a midlatitude SST front

4.4.1 Setup and reference solutions

An idealized experiment particularly relevant for the coupling of the MABL with mesoscale oceanic eddies (and potentially submesoscale fronts) was initially suggested by Spall (2007) and then revised by Kilpatrick et al. (2014). More recently Ayet and Redelsperger (2019) derived an analytical model based on a similar setup. The geometry of the problem is two-dimensional x–z with an SST front along the x axis:

$${\displaystyle }{\mathit{\theta }}_{\mathrm{s}}\left(x\right)=\mathrm{288.95}+{\displaystyle \frac{\mathrm{\Delta }\mathit{\theta }}{\mathrm{2}}}\tanh \left({\displaystyle \frac{x}{L_{\mathit{\theta }}}}\right),$$

where Δθ=3 K, L_θ=100 km, and $x\in [-\mathrm{1800}\mathrm{km},\mathrm{1800}\mathrm{km}]$. As indicated in Table 2, a zonal geostrophic wind of 15 m s⁻¹ is prescribed, balanced by a vertically homogeneous meridional pressure gradient. The wind thus flows over cold water before reaching a warm SST anomaly, which is 3 K warmer. We consider a dry case, in which the model is initialized ∀x with

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\theta }(z,t=\mathrm{0})=\mathrm{288.95}+\left(N^{\mathrm{2}}{\mathit{\theta }}^{\mathrm{ref}}/g\right)z,\\ {\displaystyle }& {\displaystyle }q(z,t=\mathrm{0})=\mathrm{0},\end{array}$$

where $N^{\mathrm{2}}={\mathrm{10}}^{-\mathrm{4}}$ s⁻² and θ^ref=288 K.

The velocities are systematically initialized with geostrophic winds. All simulations are run for 36 h when the flow reaches a quasi-equilibrium state. For this configuration the reference solution is obtained from the mesoscale non-hydrostatic model (MesoNH) v5.3.0 (Lafore et al., 1998; Lac et al., 2018), where microphysics and radiation packages have not been activated. The horizontal resolution is Δx=1 km and the model is discretized with 91 vertical levels from the surface to 20 km height. The vertical resolution near the surface is Δz=10 m and around Δz=100 m at 2000 m height. The turbulence scheme is the 1.5-order closure of Cuxart et al. (2000) in its one-dimensional form with the l_BL89 mixing length and CCH02 set of parameters. Sea surface fluxes are computed using the bulk parameterization COARE3.0, which is also available in NEMO from the AeroBulk package. As far as the ABL1d model is concerned, the top of the computational domain is z_top=2000 m and the vertical grid is stretched with a typical resolution of 20 m near the surface and 100 m near z=z_top with a first grid point located at z=10 m. In the horizontal, the resolution is Δx=6 km.

Figure 6 Zonal (a, b, c) and meridional (d, e, f) components of atmospheric winds for the reference MesoNH simulation (c, f), for ABL1d simulations with l_D80 mixing length and CCH02 model constants (a, d), and with l_BL89 mixing length and CCH02 model constants (b, e). Temperature contours are shown in white with a contour every 0.5 ^∘C between 15 and 17.5 ^∘C. The SST front is centered at x=0 km.

Figure 7 Zonal (a) and meridional (b) components of 10 m winds and 10 m temperature (c) for the reference MesoNH simulation (dashed red) and for ABL1d simulations with different mixing-length formulations for the winds across a midlatitude SST front experiment.

4.4.2 Numerical results

For this configuration, results will be mostly evaluated in terms of 10 m winds u₁₀ and temperature θ₁₀. As an illustration of the type of result we get, we first compare the MesoNH solution and the ABL1d solution obtained with the l_D80 and l_BL89 mixing lengths in Fig. 6. It is worth noting that the MesoNH solution closely compares with the solution of Kilpatrick et al. (2014) (their Fig. 2) with a shallow boundary layer height (around 400 m) before the front and a thicker one (around 800 m) after the front where momentum mixing is enhanced. Over the front, as noted by Ayet and Redelsperger (2019) with a similar setup, the effect of advection is predominant for meridional winds, thus explaining the differences seen with the ABL1d simulations. Indeed with ABL1d, whatever the numerical options, the atmospheric column will locally adjust to the underlying oceanic conditions since horizontal advection is neglected. This explains the absence of horizontal lag when passing over the front in the ABL1d solution compared to the MesoNH solution. However, away from the SST front the solutions are very similar in terms of boundary layer height and vertical wind structure. In anticipation of a coupling with an oceanic model, the most important quantities to look at are the 10 m atmospheric variables rather than the full 3D vertical structure of the MABL. In Fig. 7, the 10 m wind components and temperature when the ABL1d model reaches a quasi-equilibrium state are shown for different mixing length options, as well as the MesoNH results. First, the results obtained with the l_R17 are very different from the expected behavior, and we will focus the discussion on other options. In terms of zonal 10 m wind the buoyancy-based l_BL89 and l_D80 mixing lengths provide a good agreement with the MesoNH solution, which could be expected as the MesoNH solution has been generated using the l_BL89 mixing length. As soon as the mixing length is a function of buoyancy and vertical shear (as is the case for $l_{\mathrm{D}\mathrm{80}}^{\star }$) the simulated winds are weaker because the boundary layer is thinner. This leads to improved results in the stably stratified case shown earlier, but in the present case, which is more representative of realistic configuration in the MABL, it leads to a mixing that is too weak. However, compared to the l_R17 mixing length the $l_{\mathrm{D}\mathrm{80}}^{\star }$ still performs reasonably well but the winds on the warm side of the front are about 1 m s⁻¹ weaker than the MesoNH winds for c₀=0.15 and become weaker and weaker as c₀ increases.

The main outcomes are as follows:

• In the frontal region the effect of horizontal advection is predominant and the ABL1d model cannot reproduce the horizontal lag seen in the reference solution when passing over the front.

• The ABL1d model reproduces the downward momentum mixing mechanism correctly. The best results are obtained with the buoyancy-based l_BL89 and l_D80 mixing lengths.

• The l_R17 mixing length will be discarded from the comparison.

Although relevant for the present study this 2D x–z setup is not fully representative of realistic conditions because the air column has time to adjust to the underlying oceanic state, which is kept frozen in time.

Figure 8 2D time vs. height sections representing the temporal evolution of the zonal (a, b, c) and meridional (d, e, f) components of atmospheric winds for ABL1d simulations of an air column crossing an SST front with COARE bulk formulation (a, d) and IFS bulk formulation (middle and right panels). For the case presented in the right panels, a Newtonian relaxation toward the initial temperature profile was added with ${\mathit{\lambda }}_{\mathrm{s}}^{min}=\frac{\mathrm{1}}{\mathrm{48}\left[\mathrm{h}\right]}$ and ${\mathit{\lambda }}_{\mathrm{s}}^{max}=\frac{\mathrm{1}}{\mathrm{6}\left[\mathrm{h}\right]}$. The simulations were performed with l_D80 mixing length and CCH02 model constants.

Figure 9 Temporal evolution of the zonal (a, c, e) and meridional (b, d, f) components of 10 m atmospheric winds for ABL1d simulations of an air column crossing an SST front. The temporal evolution of SST (solid red lines) is also shown. For each panel the results from four different simulations are shown: with COARE bulk formulation (solid lines) or IFS bulk formulation (dashed lines), with Newtonian relaxation on temperature such that ${\mathit{\lambda }}_{\mathrm{s}}^{min}=\frac{\mathrm{1}}{\mathrm{48}\left[\mathrm{h}\right]}$ and ${\mathit{\lambda }}_{\mathrm{s}}^{max}=\frac{\mathrm{1}}{\mathrm{6}\left[\mathrm{h}\right]}$ (gray lines) or no relaxation (black lines). Panels (a, b) are obtained from simulations performed with l_D80 mixing length, (c, d) with $l_{\mathrm{D}\mathrm{80}}^{\star }$ (c₀=0.15), and (e, f) l_BL89.

4.4.3 A single-column version

An alternative to the x–z setup would be to formulate the test case as a Lagrangian advection of an air column over an SST front by prescribing a temporal evolution of sea surface temperature θ_s(t) as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\mathrm{s}}\left(t\right)=\mathrm{288.95}+\mathrm{3}\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\tanh \left({\displaystyle \frac{\mathrm{3}(t-\mathrm{144}\times {\mathrm{10}}^{\mathrm{3}}\mathrm{s})}{\mathrm{20}\mathrm{000}}}\right)\right],\\ {\displaystyle }& {\displaystyle }t\in [\mathrm{0},\mathrm{80}\times \mathrm{3600}\mathrm{s}].\end{array}$$

In this case the air column does not necessarily have time to adjust to the underlying oceanic conditions. Initial conditions are the same as the ones of the 2D x–z case. For this test case we do not have a reference solution, but it is expected that the temporal evolution of the solution should be relatively similar to the spatial evolution in the MesoNH 2D x–z case studied in previous subsection. This can be seen from Fig. 8, where there is a clear similarity between the time vs. height sections obtained with the ABL1d simulations and the x vs. height sections shown for MesoNH in Fig. 6. The ABL1d solution shows a temporal lag analogous to the horizontal lag in the reference solution for the 2D x–z case. In addition to that, we also use this test case to investigate the sensitivity of the solutions to the bulk formulation and to the Newtonian relaxation which was absent in simulations discussed so far. We consider the COARE and IFS bulk formulations, which are relatively close to each other, to check the robustness of the results to small perturbations in surface fluxes. We also consider simulations with a relaxation of the temperature variable toward the initial condition with a fast relaxation timescale ${\mathit{\lambda }}_{\mathrm{s}}=\frac{\mathrm{1}}{\mathrm{6}\left[\mathrm{h}\right]}$ above the ABL and a slower one ${\mathit{\lambda }}_{\mathrm{s}}=\frac{\mathrm{1}}{\mathrm{48}\left[\mathrm{h}\right]}$ in the ABL. This is meant to check that the relaxation does not completely overwrite the physics of the coupling we aim to represent with the ABL1d model. Results from those sensitivity experiments are shown in Figs. 8 and 9. In particular in Fig. 9 the evolution of the 10 m winds across the SST front closely resembles the one shown in Fig. 7 (dashed red lines) for MesoNH. Moreover, the results in Fig. 9 are robust to a change of bulk formulation to compute the surface fluxes. Reassuringly, adding a relaxation toward large-scale data which did not see the SST front does not deteriorate the realism of the solutions, as can be seen from Fig. 8c and f and 9 (gray lines).

The main outcomes are as follows:

• The response of the ABL1d model to evolving oceanic conditions is not local in time (it shows a temporal lag).

• The good representation of the downward momentum mixing process is not sensitive to the bulk formulation.

• Adding a relaxation term toward large-scale data does not deteriorate the realism of the solutions significantly.

Based on the results reported in this section, the best balance between efficiency and physical relevance is obtained when using the parameter values from CCH02 and the modified Deardorff (1980) mixing-length formulation l_D80 or $l_{\mathrm{D}\mathrm{80}}^{\star }$. In particular we could imagine using the l_D80 formulation over water and the $l_{\mathrm{D}\mathrm{80}}^{\star }$ formulation over sea ice.

5 Coupled numerical experiments

Using atmosphere-only experiments, we have been focused so far on the good representation of the downward momentum mixing mechanism and of the stable boundary layers typical of areas covered with sea ice. In the following we check that those two aspects are still adequately represented in a realistic coupled NEMO-ABL1d simulation. This simulation will also be used to look at the wind–current interaction, which was left aside so far. We performed a 5-year global simulation using the ORCA025 configuration. Details and illustrations are given hereafter.

5.1 Coupled NEMO-ABL1d configuration

We use here a global ORCA025 configuration at a 0.25^∘ horizontal resolution (Barnier et al., 2006) with 75 vertical z levels forced by the ECWMF ERA-Interim 6 h analysis (Dee et al., 2011). This configuration is identical to the one described in Couvelard et al. (2020) (see their Sect. 4.1.1). The ABL1d-NEMO coupled simulation is carried out with the same numerical options as in a standard ASL forcing strategy. However, in the ABL coupling strategy, the two-dimensional near-surface air temperature, humidity, and winds used in the usual ASL forcing are replaced by three-dimensional atmospheric variables sampled between the surface and 2000 m preprocessed following the different steps described in Appendix B. The large-scale pressure gradient computed during the preprocessing is used as a geostrophic forcing for the ABL1d model dynamics. Three-dimensional atmospheric variables are generated over the 2014–2018 period and vertically interpolated on 50 levels between 10 and 2000 m with a vertical resolution increasing with height. Grid resolution is about 20 m near the air–sea interface and reaches 70 m at the top of the ABL1d domain. The choice of a vertical extent of 2000 m and 50 vertical levels in the ABL1d model is somewhat arbitrary and the robustness of numerical results to these choices will be investigated in future studies. For the simulations presented here, the same horizontal grid and time step (Δt=1200 s) are chosen in the ABL1d and NEMO models. The options associated with the ABL coupling available through the NEMO standard name list are reported in Table 3, and a detailed profiling of the code is presented in Appendix D in order to assess the overhead associated with the ABL coupling strategy vs. ASL coupling strategy. This profiling shows that the overhead associated with the ABL1d (when using the l_D80 mixing length) is on the order of 4 % and the one associated with the input part of the I/O operations is 5 %. Overall there is an increase of 9 % in elapsed time compared to the standard ASL forcing strategy.

Table 3 Name list parameters in the NEMO(v4.0) to set in the name list section namsbc_abl before running a simulation coupled with ABL1d.

5.2 Numerical results

In this section, we evaluate the ABL coupling strategy in a realistic context for a set of relevant metrics. The objective is not to conduct a thorough physical analysis of the numerical results but to illustrate the potential of the ABL coupling strategy and its proper implementation in NEMO. To evaluate our numerical results, we use standard metrics from the literature to quantify the wind–SST (a.k.a. thermal feedback effect), wind–currents (a.k.a. current feedback effect), and MABL–sea-ice couplings (e.g., Bryan et al., 2010; Renault et al., 2019b).

5.2.1 Thermal feedback effect

To quantify the surface wind response to SST, we show in Fig. 10a a global map of the temporal correlation between the high-pass-filtered 10 m wind speed from the first vertical level in the ABL1d model and the SST. The same correlation is shown in Bryan et al. (2010) from satellite observations (their Fig. 1d) and from coupled numerical experiments between a 0.1^∘ ocean and a 0.25^∘ atmospheric model (their Fig. 1c). Consistent with observations and fully coupled models, the correlation obtained from the coupled NEMO-ABL1d simulation shows large positive correlations over regions like the Southern Ocean, Kuroshio, and Gulf Stream extensions as well as in the Gulf of Guinea. Correlations are, however, weaker than observations in the northern and equatorial Pacific between 90 and 180^∘ W. As the thermal feedback strength is related to the ocean model resolution (Bryan et al., 2010), we can expect a better agreement with observations using a higher resolution configuration such as ORCA12 ($\mathrm{1}/\mathrm{12}{}^{\circ }$ resolution). This coupling sensitivity to the oceanic resolution will be addressed in a future study.

Figure 10 (a) Global map of temporal correlation of high-pass-filtered wind speed at the first vertical level of the ABL1d model with SST from NEMO. Both NEMO and ABL1d are configured at 0.25^∘ resolution. (b, c) Global maps of the coupling coefficient between the surface current vorticity and the wind curl (s_w, b) and between the surface current vorticity and the wind-stress curl (s_τ, c) estimated from a 0.25^∘ resolution coupled NEMO-ABL1d global simulation. The fields are first temporally averaged using a 29 d running mean and spatially high-pass filtered.

5.2.2 Current feedback effect

Other processes of interest are those related to the coupling between oceanic surface currents, wind stress, and wind. Such coupling is responsible for a dampening of the eddy mesoscale activity in the ocean. In Renault et al. (2019b), two coupling coefficients called s_w and s_τ are defined to quantify this effect. s_τ is a measure of the sink of energy from the eddies and fronts to the ABL and s_w quantifies the partial re-energization of the ocean by the wind response to the wind–current coupling. This re-energization is absent in the ASL forcing strategy, which results in an excessive dampening of the oceanic eddy mesoscale activity. In practice, s_τ (s_w) corresponds to the slope of the linear relationship between high-pass-filtered surface current vorticity and surface wind-stress (wind) curl. Global maps of s_τ and s_w computed from our coupled NEMO-ABL1d global simulation are shown in Fig. 10. Large negative values of s_τ indicate an efficient dampening of the eddy mesoscale activity by the current feedback (i.e., a large sink of energy from the ocean to the atmosphere), and the large positive values of s_w indicate an efficient wind response and re-energization of the mesoscale currents. Our numerical experiment provides results very consistent with the results obtained from coupled simulations between NEMO and the WRF shown in Renault et al. (2019b) (their Fig. 1b for s_τ and 2c for s_w). As mentioned earlier, with an ASL forcing strategy we would systematically have s_w=0 and stronger s_τ values.

Figure 11 Yearly average of sea-ice cover (contours) and atmospheric boundary layer height (shaded) over the antarctic (b) and arctic (a) regions.

5.2.3 MABL and sea-ice coupling

The last illustration of our implementation presented in this section is the coupling of ABL1d with sea ice. As described in Sect. 4.3, sea ice generally induces a shallow stably stratified boundary layer due to the near-surface air cooling. This increased vertical stability tends to reduce atmospheric turbulence, producing shallower ABL heights over sea ice. This relationship between sea-ice concentration and ABL height is clearly visible from Fig. 11 on both Arctic and Antarctic domains, where the ABL height follows a progressive decrease from about 800 to 200 m in the transition zone between the open ocean and fully ice-covered regions. This coupling between the ABL and sea ice have important effects on near-surface wind, temperature, and humidity, and consequently on sea-ice concentration evolution, which will need to be specifically assessed in future ABL-based studies.

6 Conclusions

6.1 Summary

A simplified atmospheric boundary layer (ABL) model has been developed and integrated to an oceanic model. This development is made with the objective to improve the representation of air–sea interactions in eddying oceanic models compared to the standard forcing strategy where the 10 m height atmospheric quantities are prescribed. For this preliminary study, the simplified ABL model takes the form of a single-column model including a turbulence scheme coupled to each oceanic grid point. A crucial hypothesis is that the dominant process at the characteristic scale of the oceanic mesoscale is the so-called downward mixing process which stems from a modulation of atmospheric turbulence by sea surface temperature (SST) anomalies. Our approach can be seen as an extended bulk approach: instead of prescribing atmospheric quantities at 10 m to compute air–sea fluxes via an atmospheric surface layer (ASL) parameterization, atmospheric quantities in the first few hundred meters are used to constrain an ABL model which provides 10 m atmospheric values to the ASL parameterization. An important point is that our modeling strategy keeps the computational efficiency and flexibility inherent to ocean-only modeling. Indeed, the overhead generally observed in terms of computational cost compared to the usual ASL forcing strategy is roughly 10 %, and half of this overhead is due to I/O operations since the ABL model is constrained by 3D atmospheric data. In this paper the key components of such an approach have been described. This includes the large-scale forcing strategy, the coupling with the ocean and sea ice and last but not least the ABL turbulence closure scheme based on a prognostic equation for the turbulence kinetic energy. The resulting simplified model, called ABL1d, has been tested for several boundary-layer regimes relevant to either ocean–atmosphere or sea-ice–atmosphere coupling. Results have systematically been evaluated against large eddy simulations (LESs). Furthermore we have investigated the behavior of the model to several parameters including the formulation of the mixing length and the turbulence model constants. First results from a global ABL1d-NEMO configuration show an excellent behavior in terms of wind–SST two-way coupling. A first analysis of the impact of the coupling with ABL1d from a physical viewpoint is presented in Brivoal et al. (2020).

6.2 Future work

Now that an adequate computational framework and an efficient turbulent scheme that can be operated for a reasonable computational overhead have been developed, the next step is to investigate the relevance of the single-column representation of the ABL selected for the present paper. Indeed, several studies have already shown that momentum vertical turbulent mixing, pressure gradient, Coriolis, and nonlinear advection are all important to the momentum balance in the marine atmospheric boundary layer at the vicinity of oceanic fronts (see for example Spall, 2007; Small et al., 2008; O'Neill et al., 2010). It is well known that the relative importance of those terms depends on the wind regime: for strong winds the vertical mixing is the dominant mechanism while for weak winds the pressure adjustment mechanism dominates. The current single-column approximation is based on the assumption that the balance is dominated by vertical turbulent mixing, and the effect of other terms is roughly represented by the geostrophic guide and/or the nudging term. The test case presented in Sect. 4.4 clearly illustrates the limitations of a single-column approach ignoring advective effects. However, before moving to more advanced formulations, our rationale was that the two main bottlenecks in terms of computational cost inherent to the ABL coupling strategy are the reading of 3D atmospheric data and the choice of ABL turbulent scheme. As a first step, we focused on those two aspects to assess whether or not our approach can be a viable option. Even if the justification of our model is not beyond reproach, it already brings an improvement compared to the ASL forcing strategy.

Several ways to improve the methodology presented here are currently under investigation. At a practical level, ways to lower the computational overhead due to I/O operations will be investigated using the parallel I/O capabilities provided by the XIOS library which is currently used in NEMO only for outputs. At a more fundamental level, the continuous formulation of the ABL1d model will be completed to improve the representation of the momentum balance by integrating the effect of horizontal advection and fine-scale pressure gradients. Increasing the complexity of the model should allow the impact of the nudging term on the ABL solutions to be lowered. In the event our approach turns out to be physically sound for a reasonable complexity it could be useful not only for offline oceanic simulations but also in coupled simulations to downscale the information from a low-resolution atmospheric component to a high-resolution oceanic component. A standalone ABL model of intermediate complexity could also play a role in coupled data assimilation where the current practice is generally to assimilate data separately in the ocean and the atmosphere, ignoring the air–sea interactions, which results in inconsistencies at the air–sea interface in the initial conditions, causing initial shocks in the coupled forecasts (e.g., Mulholland et al., 2015). We wish to conclude this study by clarifying that the framework we have developed within NEMO is general enough to allow alternative approaches (e.g., via model-driven empirical models) to be seamlessly tested and confronted with the ABL coupling strategy.

Appendix A: Surface boundary conditions for TKE and mixing lengths

In this appendix, following the methodology of Redelsperger et al. (2001) re-expressed with our notations, we quickly recall how the surface boundary conditions for the turbulent kinetic energy (TKE) and for the mixing lengths l_m and l_ε are determined via a matching between the subgrid turbulence scheme and the surface-layer theory. The simplest form of atmospheric surface layer (ASL) theory, namely for neutral stratification, is considered. Under the quasi-equilibrium hypothesis the evolution Eq. (6) for the TKE reduces to the equilibrium between turbulence production and dissipation $K_{\mathrm{m}}\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}{\Vert }^{\mathrm{2}}=\frac{c_{\mathit{\epsilon }}}{l_{\mathit{\epsilon }}}{e}^{\mathrm{3}/\mathrm{2}}$ which, combined with $K_{\mathrm{m}}=C_{\mathrm{m}}{l}_{\mathrm{m}}\sqrt{e}$, leads to

$$\begin{array}{}\text{(A1)}& e={\displaystyle \frac{C_{\mathrm{m}}}{c_{\mathit{\epsilon }}}}\left(l_{\mathrm{m}}{l}_{\mathit{\epsilon }}\right)\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}{\Vert }^{\mathrm{2}}.\end{array}$$

The similarity theory for the ASL in the neutral case is such that $\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}\Vert =\frac{u_{\star }}{\mathit{\kappa }(z+z_{\mathrm{0}})}$ with κ the Von Karman constant and z₀ the roughness length which can be combined with Eq. (A1) to get that in the surface layer:

$$\begin{array}{}\text{(A2)}& e={\displaystyle \frac{C_{\mathrm{m}}}{c_{\mathit{\epsilon }}}}\left({\displaystyle \frac{l_{\mathrm{m}}{l}_{\mathit{\epsilon }}}{{\mathit{\kappa }}^{\mathrm{2}}(z+z_{\mathrm{0}}{)}^{\mathrm{2}}}}\right)u_{\star }^{\mathrm{2}},\end{array}$$

with $u_{\star }^{\mathrm{2}}=\Vert 〈{\mathbf{u}}_{\mathrm{h}}^{\prime }{w}^{\prime }〉\Vert$ the friction velocity. Moreover, enforcing the consistency between the eddy diffusivity for momentum given by the ASL theory ($K_{\mathrm{m}}=\mathit{\kappa }{u}_{\star }(z+z_{\mathrm{0}})$) and the one given by the TKE closure ($K_{\mathrm{m}}=C_{\mathrm{m}}{l}_{\mathrm{m}}\sqrt{e}$) leads to

$$\begin{array}{}\text{(A3)}& l_{\mathrm{m}}={\displaystyle \frac{\mathit{\kappa }(z+z_{\mathrm{0}})}{C_{\mathrm{m}}}}\left({\displaystyle \frac{u_{\star }}{\sqrt{e}}}\right)\to l_{\mathrm{m}}^{\mathrm{3}}{l}_{\mathit{\epsilon }}=c_{\mathit{\epsilon }}{C}_{\mathrm{m}}^{-\mathrm{3}}{\left(\mathit{\kappa }(z+z_{\mathrm{0}})\right)}^{\mathrm{4}}.\end{array}$$

We thus have two relations (A2) and (A3) for three unknowns (e, l_m, and l_ε). At this point our derivation will differ from Redelsperger et al. (2001) as we will assume that $l_{\mathrm{m}}=l_{\mathit{\epsilon }}=L_{\mathrm{sfc}}$ in the ASL.

Under this assumption, combining Eqs. (A2) and (A3) we easily obtain

$${\displaystyle }e={\displaystyle \frac{u_{\star }^{\mathrm{2}}}{\sqrt{C_{\mathrm{m}}{c}_{\mathit{\epsilon }}}}},L_{\mathrm{sfc}}\left(z\right)=\mathit{\kappa }{\displaystyle \frac{(C_{\mathrm{m}}{c}_{\mathit{\epsilon }}{)}^{\mathrm{1}/\mathrm{4}}}{C_{\mathrm{m}}}}(z+z_{\mathrm{0}}),$$

for z≤δ_asl, where δ_asl corresponds to the extent of the ASL.

The expression of L_sfc(z) for z≤δ_asl is also used as a constraint to define the weighted average needed to determine l_m from l_up and l_dwn:

$${\displaystyle }{l}_{\mathrm{m}}={\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(l_{\mathrm{up}}^{\mathrm{1}/a}+l_{\mathrm{dwn}}^{\mathrm{1}/a}\right)\right)}^a$$

(equivalent to Eq. 12a). In the ASL we further assume that l_dwn≈0 and l_up≈δ_asl; for l_m(z=δ_asl) to be consistent with L_sfc(z=δ_asl) we should have

$$\begin{array}{}\text{(A4)}& {\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\delta }}_{\mathrm{asl}}^{\mathrm{1}/a}\right)\right)}^a=\mathit{\kappa }{\displaystyle \frac{(C_{\mathrm{m}}{c}_{\mathit{\epsilon }}{)}^{\mathrm{1}/\mathrm{4}}}{C_{\mathrm{m}}}}({\mathit{\delta }}_{\mathrm{asl}}+z_{\mathrm{0}}).\end{array}$$

Considering that δ_asl≫z₀, Eq. (A4) is satisfied for

$${\displaystyle }a=-\left({\displaystyle \frac{\log \left(c_{\mathit{\epsilon }}\right)-\mathrm{3}\log \left(C_{\mathrm{m}}\right)+\mathrm{4}\log \left(\mathit{\kappa }\right)}{\log \left(\mathrm{16}\right)}}\right).$$

Considering the CBR00 model constants we obtain $a=-\mathrm{1.4796}\approx -\mathrm{3}/\mathrm{2}$ and $a=-\mathrm{0.860834}\approx -\mathrm{6}/\mathrm{7}$ for the CCH02 constants (see Table 1).

Appendix B: Preprocessing of atmospheric data from IFS

B1 Altitude of IFS vertical levels

The ABL1d model is discretized on fixed in time and space geopotential levels while the IFS model uses a pressure-based sigma coordinate. A first step is to recover the altitude associated with each sigma level. The pressure p_k+½ defined at cell interfaces between two successive vertical layers is given by

$${\displaystyle }{p}_{k+\frac{\mathrm{1}}{\mathrm{2}}}=A_{k+\frac{\mathrm{1}}{\mathrm{2}}}+B_{k+\frac{\mathrm{1}}{\mathrm{2}}}{p}_{\mathrm{s}},k\in \left[\right[\mathrm{1},N_{\mathrm{ifs}}\left]\right],$$

where A_k+½ (Pa) and B_k+½ (dimensionless) are constants given by a smooth analytical function defining the vertical grid stretching. Typical values of the altitude of grid points in the vertical for a standard 60-level grid (L60) and a surface pressure of 1013 hPa are given in Table B1. Once the values of p_k+½ and p_s are known, the altitude of cell interfaces can be computed by integrating the hydrostatic equilibrium

$$\begin{array}{}\text{(B1)}& {\partial }_z\mathit{\varphi }=-{\displaystyle \frac{R_{\mathrm{d}}{T}_{\mathrm{v}}}{p}}{\partial }_{z}p\end{array}$$

vertically. In Eq. (B1), ϕ is the geopotential, T_v the virtual temperature, and R_d the specific gas constant for dry air. At a discrete level we get

$${\displaystyle }\underset{z_{k-\frac{\mathrm{1}}{\mathrm{2}}}}{\overset{z_{k+\frac{\mathrm{1}}{\mathrm{2}}}}{\int }}{\partial }_{\mathrm{s}}\mathit{\varphi }\mathrm{d}s=-R_{\mathrm{d}}{T}_{\mathrm{v}}\left(z_k\right)\underset{z_{k-\frac{\mathrm{1}}{\mathrm{2}}}}{\overset{z_{k+\frac{\mathrm{1}}{\mathrm{2}}}}{\int }}\left({\displaystyle \frac{{\partial }_{\mathrm{s}}p}{p}}\right)\mathrm{d}s,$$

which gives

$${\displaystyle }\mathrm{e}\mathrm{3}{\mathrm{t}}_k^{\mathrm{ifs}}=-{\displaystyle \frac{R_{\mathrm{d}}{T}_{\mathrm{v}}\left(z_k\right)}{g}}\ln \left({\displaystyle \frac{p_{k+\frac{\mathrm{1}}{\mathrm{2}}}}{p_{k-\frac{\mathrm{1}}{\mathrm{2}}}}}\right).$$

Once the layer thicknesses $\mathrm{e}\mathrm{3}{\mathrm{t}}_k^{\mathrm{ifs}}$ are known, horizontal wind components, potential temperature, and specific humidity can be interpolated on the ABL1d vertical levels. Under the constraint that ${\int }_{z_{\mathrm{sfc}}}^{z_{\mathrm{top}}}{\mathit{\psi }}^{\mathrm{ifs}}\left(z\right)\mathrm{d}z={\int }_{z_{\mathrm{sfc}}}^{z_{\mathrm{top}}}\mathit{\psi }\left(z\right)\mathrm{d}z$ for any IFS quantity ψ^ifs to be interpolated. Wind components are interpolated using a fourth-order compact scheme while tracers are interpolated using a WENO-like PPM scheme (Alexander Shchepetkin, personal communication, 2001) which is monotonic.

Table B1 Altitude z_k and layer thickness e3t_k of the IFS L60 vertical grid in the first 2000 m with respect to the parameter values A_k and B_k of a surface pressure p_s=1013 hPa.

B2 Filtering in the presence of boundaries

Because of the IFS numerical formulation and of the post-processing of output data, the solutions sometimes contain high-frequency oscillations at the vicinity of the land-sea interface. This problem is further compounded when the nearshore topography is steep. The atmospheric fields over water thus need to be smoothed horizontally to specifically remove the 2Δx noise. We use a standard two-dimensional Shapiro filter which, in the absence of lateral boundaries, can be formulated as

$$\begin{array}{l}{\displaystyle }{\mathit{\psi }}_{i,j}^{\star }={\mathit{\psi }}_{i,j}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left({\mathit{\delta }}_{i+\mathrm{1}/\mathrm{2},j}^{\left(x\right)}-{\mathit{\delta }}_{i-\mathrm{1}/\mathrm{2},j}^{\left(x\right)}\right),\\ {\displaystyle }{\mathit{\psi }}_{i,j}^{\mathrm{f}}={\mathit{\psi }}_{i,j}^{\star }+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left({\mathit{\delta }}_{i,j+\mathrm{1}/\mathrm{2}}^{(y,\star )}-{\mathit{\delta }}_{i,j-\mathrm{1}/\mathrm{2}}^{(y,\star )}\right),\end{array}$$

where ${\mathit{\delta }}_{i+\mathrm{1}/\mathrm{2},j}^{\left(x\right)}={\mathit{\psi }}_{i+\mathrm{1},j}-{\mathit{\psi }}_{i,j}$ and ${\mathit{\delta }}_{i,j+\mathrm{1}/\mathrm{2}}^{(y,\star )}={\mathit{\psi }}_{i,j+\mathrm{1}}^{\star }-{\mathit{\psi }}_{i,j}^{\star }$. The amplification factor associated with this filter is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathcal{A}}_{\mathrm{shap}}({\mathit{\theta }}_x,{\mathit{\theta }}_y)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left(\mathrm{1}+\cos {\mathit{\theta }}_x\right)\left(\mathrm{1}+\cos {\mathit{\theta }}_y\right),{\mathit{\theta }}_x=k_x\mathrm{\Delta }x,\\ {\displaystyle }& {\displaystyle }{\mathit{\theta }}_y=k_y\mathrm{\Delta }y,\end{array}$$

which guarantees that one iteration of the filter is sufficient to remove the grid-scale noise since ${\mathcal{A}}_{\mathrm{shap}}(\mathit{\pi },\mathit{\pi })={\mathcal{A}}_{\mathrm{shap}}(\mathit{\pi },{\mathit{\theta }}_y)={\mathcal{A}}_{\mathrm{shap}}({\mathit{\theta }}_x,\mathit{\pi })=\mathrm{0}$ and that 𝒜_shap≤1 (i.e., no waves are amplified). In the presence of solid boundaries we would like to retain those properties as much as possible. A straightforward approach would be to impose a no-gradient condition at the coast, i.e., ${\mathit{\delta }}_{i+\mathrm{1}/\mathrm{2},j}^{\left(x\right)}=\mathrm{0}$ as soon as ${\mathrm{tmask}}_{i+\mathrm{1},j}\times {\mathrm{tmask}}_{i,j}=\mathrm{0}$ (${\mathit{\delta }}_{i,j+\mathrm{1}/\mathrm{2}}^{(y,\star )}=\mathrm{0}$ as soon as ${\mathrm{tmask}}_{i,j+\mathrm{1}}\times {\mathrm{tmask}}_{i,j}=\mathrm{0}$), with tmask the indicator function equal to 1 over water and 0 over land. Let us also consider the following alternative boundary conditions

$$\begin{array}{}\text{(B2)}& \left\{\begin{array}{r}{\mathit{\delta }}_{i+\mathrm{1}/\mathrm{2},j}^{\left(x\right)}=-{\mathit{\delta }}_{i-\mathrm{1}/\mathrm{2},j}^{\left(x\right)},\mathrm{if}{\mathrm{tmask}}_{i+\mathrm{1},j}=\mathrm{0},\\ {\mathit{\delta }}_{i-\mathrm{1}/\mathrm{2},j}^{\left(x\right)}=-{\mathit{\delta }}_{i+\mathrm{1}/\mathrm{2},j}^{\left(x\right)},\mathrm{if}{\mathrm{tmask}}_{i-\mathrm{1},j}=\mathrm{0},\end{array}\right.\end{array}$$

and similar in the y direction. We do not elaborate on this choice but it can be shown theoretically that boundary conditions (B2) provide a better control of grid-scale noise near the coast. To illustrate this point, in Fig. B1 the surface pressure gradients are shown for different boundary conditions. In particular it can be seen near the coast that the no-gradient boundary condition (panels b and e) leaves some artificial patterns in gradients, especially in the Peru–Chile current system, while the boundary condition (B2) efficiently mitigates this issue. Note that it is particularly essential to make sure that the surface pressure field is sufficiently smooth because gradients of this field are used to compute geostrophic winds which are important for the large-scale forcing of the ABL1d model.

Figure B1 Atmospheric surface pressure horizontal gradients in x (a, b, c) and y (d, e, f) directions obtained from the original IFS data (a, d), after a Shapiro filtering with no-gradient boundary conditions (b, e), and after a Shapiro filtering with boundary conditions (B2) (c, f). The area in red is covered by land.

B3 Large-scale pressure gradient computation

The last aspect of the pre-processing of atmospheric data we would like to discuss is the computation of the large-scale pressure gradient (or equivalently of the geostrophic wind components) The objective is to estimate the following terms:

$${\displaystyle }{R}_{\mathrm{LS}}^u={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{a}}}}({\partial }_{x}p{)}_z,R_{\mathrm{LS}}^v={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{a}}}}({\partial }_{y}p{)}_z,$$

where (⋅)_z denotes a gradient along constant geopotential height. Using the hydrostatic balance we have $\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{a}}}=-g({\partial }_{z}p{)}^{-\mathrm{1}}$, which leads to

$$\begin{array}{}\text{(B3)}& R_{\mathrm{LS}}^u=-g\left({\partial }_{z}p{)}^{-\mathrm{1}}\right({\partial }_{x}p{)}_z,R_{\mathrm{LS}}^v=-g\left({\partial }_{z}p{)}^{-\mathrm{1}}\right({\partial }_{y}p{)}_z.\end{array}$$

Assuming a generalized vertical coordinate $s=s(x,y)$ the computation of gradients along constant height is not straightforward since $({\partial }_{x}p{)}_z=({\partial }_{x}p{)}_{\mathrm{s}}-\left({\partial }_{z}p\right)({\partial }_{x}z{)}_{\mathrm{s}}$ leading to

$${\displaystyle }\left({\partial }_{x}p{)}_z\right({\partial }_{z}p{)}^{-\mathrm{1}}=\left({\partial }_{z}p{)}^{-\mathrm{1}}\right({\partial }_{x}p{)}_{\mathrm{s}}-({\partial }_{x}z{)}_{\mathrm{s}}.$$

In the particular case of the IFS coordinate s we have

$$\begin{array}{}\text{(B4)}& \left({\partial }_{z}p{)}^{-\mathrm{1}}\right({\partial }_{x}p{)}_{\mathrm{s}}={\displaystyle \frac{B\left(z\right){\partial }_x{p}_{\mathrm{s}}}{\left({\partial }_{z}A\right)+\left({\partial }_{z}B\right)p_{\mathrm{s}}}},\end{array}$$

and (∂_xz)_s can be estimated after integrating the hydrostatic balance. Starting from the layer interface height $z_{i,j,k+\mathrm{1}/\mathrm{2}}^{\mathrm{ifs}}$, surface pressure (p_s)_i,j, and parameter values $A_k,B_k,A_{k+\mathrm{1}/\mathrm{2}},B_{k+\mathrm{1}/\mathrm{2}}$ the different steps are the following:

1. compute Δx_i,j and Δy_i,j from latitudes and longitudes;

2. compute horizontal gradients ∂_xp_s and ∂_yp_s for surface pressure

   $$\begin{array}{l}{\displaystyle }{\mathrm{FX}}_{i+\mathrm{1}/\mathrm{2},j}={\displaystyle \frac{\mathrm{2}\left\{(p_{\mathrm{s}}{)}_{i+\mathrm{1},j}-(p_{\mathrm{s}}{)}_{i,j}\right\}}{\mathrm{\Delta }{x}_{i,j}+\mathrm{\Delta }{x}_{i+\mathrm{1},j}}},\\ {\displaystyle }{\mathrm{FY}}_{i,j+\mathrm{1}/\mathrm{2}}={\displaystyle \frac{\mathrm{2}\left\{(p_{\mathrm{s}}{)}_{i,j+\mathrm{1}}-(p_{\mathrm{s}}{)}_{i,j}\right\}}{\mathrm{\Delta }{y}_{i,j}+\mathrm{\Delta }{y}_{i,j+\mathrm{1}}}};\end{array}$$

3. compute horizontal gradients (∂_xz)_s and (∂_yz)_s

   $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{dZx}}_{i+\mathrm{1}/\mathrm{2},j,k}=\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{z_{i+\mathrm{1},j,k+\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}-z_{i,j,k+\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}+z_{i+\mathrm{1},j,k-\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}-z_{i,j,k-\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}}{\mathrm{\Delta }{x}_{i,j}+\mathrm{\Delta }{x}_{i+\mathrm{1},j}}},\\ {\displaystyle }& {\displaystyle }{\mathrm{dZy}}_{i,j+\mathrm{1}/\mathrm{2},k}=\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{z_{i,j+\mathrm{1},k+\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}-z_{i,j,k+\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}+z_{i,j+\mathrm{1},k-\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}-z_{i,j,k-\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}}{\mathrm{\Delta }{y}_{i,j}+\mathrm{\Delta }{y}_{i,j+\mathrm{1}}}};\end{array}$$

4. compute $\left({\partial }_{z}p{)}^{-\mathrm{1}}\right({\partial }_{x}p{)}_{\mathrm{s}}$ via (B4)

   $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{wrkX}}_{i,j,k}=\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{B_k\left(z_{i,j,k+\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}-z_{i,j,k-\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}\right)\left({\mathrm{FX}}_{i+\mathrm{1}/\mathrm{2},j}+{\mathrm{FX}}_{i-\mathrm{1}/\mathrm{2},j}\right)}{\left(p_s{)}_{i,j}\right(B_{k+\mathrm{1}/\mathrm{2}}-B_{k-\mathrm{1}/\mathrm{2}})+(A_{k+\mathrm{1}/\mathrm{2}}-A_{k-\mathrm{1}/\mathrm{2}})}},\\ {\displaystyle }& {\displaystyle }{\mathrm{wrkY}}_{i,j,k}=\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{B_k\left(z_{i,j,k+\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}-z_{i,j,k-\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}\right)\left({\mathrm{FY}}_{i+\mathrm{1}/\mathrm{2},j}+{\mathrm{FY}}_{i-\mathrm{1}/\mathrm{2},j}\right)}{\left(p_s{)}_{i,j}\right(B_{k+\mathrm{1}/\mathrm{2}}-B_{k-\mathrm{1}/\mathrm{2}})+(A_{k+\mathrm{1}/\mathrm{2}}-A_{k-\mathrm{1}/\mathrm{2}})}};\end{array}$$

5. finalize (we get a minus sign in $R_{\mathrm{LS}}^u$ because the grid in the y direction is flipped in the raw data)

   $$\begin{array}{ll}{\displaystyle }& {\displaystyle }(R_{\mathrm{LS}}^u{)}_{i,j,k}=\\ {\displaystyle }& {\displaystyle }-g\left({\mathrm{wrkY}}_{i,j,k}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathrm{dZy}}_{i,j+\mathrm{1}/\mathrm{2},k}+{\mathrm{dZy}}_{i,j-\mathrm{1}/\mathrm{2},k})\right),\\ {\displaystyle }& {\displaystyle }(R_{\mathrm{LS}}^v{)}_{i,j,k}=\\ {\displaystyle }& {\displaystyle }-g\left({\mathrm{wrkX}}_{i,j,k}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathrm{dZx}}_{i+\mathrm{1}/\mathrm{2},j,k}+{\mathrm{dZx}}_{i-\mathrm{1}/\mathrm{2},j,k})\right).\end{array}$$

Appendix C: Discrete algorithms to compute l_up and l_dwn

In the following we describe the discrete algorithms used to provide the mixing lengths l_up and l_dwn given in Sect. 3.2.2. Four different ways to compute those quantities have been implemented in the ABL1d model.

C1  Bougeault and Lacarrère (1989) length scale

The Bougeault and Lacarrère (1989) mixing length defines l_up and l_dwn as

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}& \underset{z}{\overset{z+l_{\mathrm{up}}}{\int }}{N}^{\mathrm{2}}\left(s\right)(s-z)\mathrm{d}s=e\left(z\right),\\ & \underset{z-l_{\mathrm{dwn}}}{\overset{z}{\int }}{N}^{\mathrm{2}}\left(s\right)(z-s)\mathrm{d}s=e\left(z\right).\end{array}\end{array}$$

By construction such mixing lengths are bounded by the distance to the bottom and the top of the computational domain and revert to the Deardorff (1980) length scale (i.e., $l_{\mathrm{up}}=l_{\mathrm{dwn}}=\sqrt{\mathrm{2}e\left(z\right)/N^{\mathrm{2}}}$) for N²=cste. An objective is to also satisfy this last property at a discrete level. Considering a simple trapezoidal rule to approximate the integral in Eq. (C1) over each grid cell, the procedure for the computation of $l_{\mathrm{up}}\left(z_{k+\mathrm{1}/\mathrm{2}}\right)$ is given in Algorithm 1. In the case $N^{\mathrm{2}}\left(z_{p+\mathrm{1}/\mathrm{2}}\right)=N^{\mathrm{2}}\left(z_{p-\mathrm{1}/\mathrm{2}}\right)=N_{\mathrm{cst}}^{\mathrm{2}}(\forall p)$, Algorithm 1 gives the following sequence:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{FC}\left(z_{k+\mathrm{1}/\mathrm{2}}\right)=-e\left(z_{k+\mathrm{1}/\mathrm{2}}\right),\\ {\displaystyle }& {\displaystyle }\mathrm{FC}\left(z_{k+\mathrm{3}/\mathrm{2}}\right)=-e\left(z_{k+\mathrm{1}/\mathrm{2}}\right)+N_{\mathrm{cst}}^{\mathrm{2}}{\displaystyle \frac{\mathrm{e}\mathrm{3}\mathrm{t}(z_{k+\mathrm{1}}{)}^{\mathrm{2}}}{\mathrm{2}}},\\ {\displaystyle }& {\displaystyle }\mathrm{FC}\left(z_{k+\mathrm{5}/\mathrm{2}}\right)=-e\left(z_{k+\mathrm{1}/\mathrm{2}}\right)+N_{\mathrm{cst}}^{\mathrm{2}}{\displaystyle \frac{{\left[\mathrm{e}\mathrm{3}\mathrm{t}\left(z_{k+\mathrm{1}}\right)+\mathrm{e}\mathrm{3}\mathrm{t}\left(z_{k+\mathrm{2}}\right)\right]}^{\mathrm{2}}}{\mathrm{2}}}\\ {\displaystyle }& {\displaystyle }\mathrm{\dots }\end{array}$$

As soon as $\mathrm{FC}\left(z_{p+\mathrm{1}/\mathrm{2}}\right)$ changes sign we stop the procedure because l_up such that $-e\left(z_{k+\mathrm{1}/\mathrm{2}}\right)+N_{\mathrm{cst}}^{\mathrm{2}}{l}_{\mathrm{up}}^{\mathrm{2}}=\mathrm{0}$, which corresponds to the Deardorff (1980) length scale, has been found. We note l_BL89 the mixing length corresponding to the Bougeault and Lacarrère (1989) algorithm.

C2 Adaptation of NEMO's length scale

The standard NEMO algorithm (Sect. 10.1.3 in Madec, 2012) is much easier to discretize. As a first step the Deardorff (1980) length scale l_D80 is computed at cell interfaces, such that

$${\displaystyle }{\left(l_{\mathrm{D}\mathrm{80}}\right)}_{k+\mathrm{1}/\mathrm{2}}=max\left(\sqrt{{\displaystyle \frac{\mathrm{2}{e}_{k+\mathrm{1}/\mathrm{2}}}{max\left(N^{\mathrm{2}},N_{\mathit{\epsilon }}^{\mathrm{2}}\right)}}},l_{min}\right),$$

with $N_{\mathit{\epsilon }}^{\mathrm{2}}$ the minimum stratification allowed whose value is set to the smallest positive real computer value. The vertical gradients of l_D80 are then limited such that they stay smaller than the variations of height. This amounts to compute l_up and l_dwn as

$$\begin{array}{}\text{(C2)}& {\displaystyle }& {\displaystyle }{\left(l_{\mathrm{up}}\right)}_{k-\mathrm{1}/\mathrm{2}}=min\left({\left(l_{\mathrm{up}}\right)}_{k+\mathrm{1}/\mathrm{2}}+\mathrm{e}\mathrm{3}{\mathrm{t}}_k,{\left(l_{\mathrm{D}\mathrm{80}}\right)}_{k-\mathrm{1}/\mathrm{2}}\right),\text{(C3)}& {\displaystyle }& {\displaystyle }{\left(l_{\mathrm{dwn}}\right)}_{k+\mathrm{1}/\mathrm{2}}=min\left({\left(l_{\mathrm{dwn}}\right)}_{k-\mathrm{1}/\mathrm{2}}+\mathrm{e}\mathrm{3}{\mathrm{t}}_k,{\left(l_{\mathrm{D}\mathrm{80}}\right)}_{k+\mathrm{1}/\mathrm{2}}\right),\end{array}$$

with e3t_k the thickness of vertical layer k (Fig. 2). The resulting mixing length is simply referred to as l_D80.

C3 Rodier et al. (2017) length scale

Recently, Rodier et al. (2017) proposed to add a shear-related term to Eq. (C1) such that the definition of l_up and l_dwn becomes

$$\begin{array}{}\text{(C4)}& \begin{array}{rl}& \underset{z}{\overset{z+l_{\mathrm{up}}}{\int }}\left[N^{\mathrm{2}}\left(s\right)(s-z)+c_{\mathrm{0}}\sqrt{e\left(s\right)}\Vert {\partial }_{\mathrm{s}}{\mathbf{u}}_{\mathrm{h}}\Vert \right]\mathrm{d}s=e\left(z\right),\\ & \underset{z-l_{\mathrm{dwn}}}{\overset{z}{\int }}\left[N^{\mathrm{2}}\left(s\right)(z-s)+c_{\mathrm{0}}\sqrt{e\left(s\right)}\Vert {\partial }_{\mathrm{s}}{\mathbf{u}}_{\mathrm{h}}\Vert \right]\mathrm{d}s=e\left(z\right),\end{array}\end{array}$$

where c₀ is a parameter whose value should be smaller than $\sqrt{C_{\mathrm{m}}/c_{\mathit{\epsilon }}}$. At a discrete level, the FC function in Algorithm 1 is replaced by

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{FC}\left(z_{p+\mathrm{1}/\mathrm{2}}\right)=\mathrm{FC}\left(z_{p-\mathrm{1}/\mathrm{2}}\right)+{\displaystyle \frac{\mathrm{e}\mathrm{3}\mathrm{t}\left(z_p\right)}{\mathrm{2}}}\left(N^{\mathrm{2}}\left(z_{p+\mathrm{1}/\mathrm{2}}\right)(z_{p+\mathrm{1}/\mathrm{2}}-\right.\\ {\displaystyle }& {\displaystyle }\left.z_{k+\mathrm{1}/\mathrm{2}})+N^{\mathrm{2}}(z_{p-\mathrm{1}/\mathrm{2}}\left)\right(z_{p-\mathrm{1}/\mathrm{2}}-z_{k+\mathrm{1}/\mathrm{2}})\right)\\ {\displaystyle }& {\displaystyle }+c_{\mathrm{0}}{\displaystyle \frac{\mathrm{e}\mathrm{3}\mathrm{t}\left(z_p\right)}{\mathrm{2}}}\left(\sqrt{e\left(z_{p+\mathrm{1}/\mathrm{2}}\right)}\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}\left(z_{p+\mathrm{1}/\mathrm{2}}\right)\Vert \right.\\ {\displaystyle }& {\displaystyle }\left.+\sqrt{e\left(z_{p-\mathrm{1}/\mathrm{2}}\right)}\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}\left(z_{p-\mathrm{1}/\mathrm{2}}\right)\Vert \right).\end{array}$$

This mixing length will be referred to as l_R17.

C4 A local buoyancy- and shear-based length scale

For the sake of computational efficiency, we have derived a local version of the Rodier et al. (2017) length scale (denoted as $l_{\mathrm{D}\mathrm{80}}^{\star }$) which is original to the present paper:

$${\displaystyle }{l}_{\mathrm{D}\mathrm{80}}^{\star }\left(z\right)={\displaystyle \frac{\mathrm{2}\sqrt{e\left(z\right)}}{c_{\mathrm{0}}\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}\Vert +\sqrt{c_{\mathrm{0}}^{\mathrm{2}}\Vert {\partial }_z{\mathbf{u}}_{\mathrm{h}}{\Vert }^{\mathrm{2}}+\mathrm{2}{N}^{\mathrm{2}}\left(z\right)}}}.$$

Once $l_{\mathrm{D}\mathrm{80}}^{\star }$ has been computed at cell interfaces $z=z_k+\mathrm{1}/\mathrm{2}$ we apply the limitations (C2) and (C3) as in the NEMO algorithm.

Appendix D: Code performance and profiling

To finalize our description of the implementation of the simplified atmospheric boundary layer model in NEMO, we assess in this appendix the computational efficiency of our approach. We compare the performance of two simulations: one with a coupling with the ABL1d model (with 50 vertical levels) which requires reading 3D atmospheric data in input files, and one with a standard ASL forcing strategy which necessitates reading only 2D atmospheric data. For the coupling with ABL1d, we consider the l_D80 mixing length which gave robustly good results across the different numerical tests investigated earlier in the paper. The simulations are performed with NEMO version 4.0 for the ORCA025 configuration previously described on 128 cores (Intel(R) Xeon(R) E5 processors 2.6 GHz) compiled with ifort (v13.0.1) using the “$-\mathrm{i}\mathrm{4}-\mathrm{r}\mathrm{8}-\mathrm{O}\mathrm{3}-\mathrm{fp}-\mathrm{model}\mathrm{precise}-\mathrm{fno}-\mathrm{alias}$” options. The I/Os are handled via the Lustre file system. Each MPI subdomain has 80×130 points in the horizontal and 75 points in the vertical. The various reports given below have been obtained from a built-in NEMO code instrumentation dedicated to calculation measurement (e.g., Maisonnave and Masson, 2019). As mentioned earlier, the outputs are done using the parallel I/O capabilities provided by the XIOS library. Thanks to XIOS, we do not expect any significant difference between the two simulations regarding the cost of output operations. However, the use of XIOS to handle input operations is still under development, and because of the significant amount of data to read in the ABL coupling strategy it makes sense to assess the associated overhead. We ran the ASL forced and ABL coupled NEMO simulations for 20 d such that the cost of the initialization step is no longer visible in the averaged cost per time step. Moreover, for the two simulations, the atmospheric data necessary for the computation of the turbulent components of air–sea fluxes are provided every 6 h.

We first show in Fig. D1 the elapsed time for each time step over the first 48 h of the simulations with different ways to specify the surface fluxes. For most time steps, the overhead associated with the ABL1d when using the l_D80 mixing length is very small (on the order of 4 %); however, every 18 time steps (i.e., every 6 h), there is a larger overhead due to the input part of the I/O operations. To further refine our assessment, we report in Table D1 the elapsed and CPU time spent on average over all the processors in the 11 most expansive sections of the code. As expected, the CPU time is not significantly affected by the ABL1d model (increase of 4 %), but the elapsed time is increased by about 9 % because of the time spent in waiting for I/O operations. The overhead associated with input operations could be mitigated by reducing the number of vertical levels in the ABL1d model (we used 50 levels here to get an upper bound on the computational overhead) and either by using XIOS to handle input operations or by running ABL1d in detached mode as explained in Sect. 3.4 such that the time spent reading input files is covered by actual computations. Nonetheless the small increase in CPU time leaves room for further improvements of the ABL model to relax the horizontal homogeneity assumption.

Figure D1 Elapsed time for each time step of a 48 h simulation with standard ASL forcing strategy (black circles) and ABL forcing strategy using the l_D80 mixing length (gray diamonds). For the two simulations the time step size is Δt=1200 s.

Table D1 Report of the elapsed time and CPU time in different sections of the NEMO (v4.0) code for the ASL forcing strategy (left portion of the table) and the ABL coupling strategy (right portion of the table). The timing is averaged on all processors. The right-most column provides a quick description of the task handled by the corresponding section. On top of the timing in seconds the percentage of the total CPU and elapsed time associated with each section is reported in parentheses. The computational overhead associated with the ABL coupling strategy can be estimated from the sbc section and the elapsed and CPU time.