Two-sided turbulent surface layer parameterizations for computing air-sea fluxes

Standard methods for determining air – sea ﬂuxes typically

rely on bulk algorithms set in the frame of Monin-Obukhov stability theory (MOST), using ocean surface fields and atmosphere near-surface fields.In the context of coupled ocean -atmosphere simulations, the shallowest ocean vertical level is usually used as bulk input and by default, the turbulent closure is one-sided: it extrapolates atmosphere near-surface solution profiles (for wind speed, temperature and humidity) to the prescribed ocean surface values.Using near-surface ocean fields as surface ones is equivalent to considering that in the ocean surface layer, solution profiles are constant instead of also being determined by a turbulent closure.Here we introduce a method for extending existing turbulent parameterizations to a two-sided framework by explicitely including the ocean surface layer within the aforementioned parameterizations.The formalism we use for this method is derived from that of classical turbulent closures, so that our novelties can easily be implemented within existing formulations.Special care is taken to ensure the smoothness of resulting solution profiles.Other physical phenomena, such as the penetration of radiative fluxes in the ocean and the formation of waves, are then included within our formalism, and their effects are assessed.
We also investigate the impact of such two-sided bulk formulations on air -sea fluxes evaluated from a setting similar to those of coupled ocean -atmosphere simulations.

K E Y W O R D S
Turbulent parameterizations, air-sea fluxes, bulk formulae, ocean surface layer, numerical methods, ocean-atmosphere coupling

| INTRODUCTION
Air-sea interactions play a crucial role for the evolution of the Earth system at both meteorological (e.g.Emanuel, 1986) and climatological (e.g.Neelin et al., 1992) scales.In climate models, the interactions between these Earth system components are primarily conveyed through the exchange of momentum, mass and heat fluxes.A significant part of these fluxes is linked to turbulent processes in the surface layer (SL) of the atmosphere and ocean, roughly defined as between 1 m below and 10 m above their common interface.The specific physical processes of the SL are central in determining turbulent air-sea fluxes, which are then used as boundary conditions for ocean and atmosphere models.The procedures for obtaining turbulent air-sea fluxes from large-scale quantities (extractable from climate models) are referred to as bulk closures (e.g.Fairall et al., 2002;Large, 2006).Consequently, establishing an adequate parameterization of the SL between the atmosphere and ocean has been a steady point of interest for the development of numerical weather and climate prediction models.
The overwhelming majority of air-sea SL parameterization schemes are expressed in the framework of Monin-Obukhov similarity theory (MOST, Monin and Obukhov, 1954) applied to the near-surface atmosphere.While additional physical effects have gradually been included within SL parameterizations, the fundamental hypotheses at their core have persisted: the atmosphere SL is described by the constant flux layer approximation, resulting from a combination of "law of the wall" and buoyancy effects, which calls for characterizing the surface roughness and column stability.Airsea fluxes parameterizations have mostly been designed for atmosphere circulation models, assuming ocean surface properties to be known.However, over the last decades, several ocean-specific processes have been progressively integrated within flux computations.Saunders (1967) first distinguished subsurface temperatures (at depth O (1 m)) from skin temperature (at depth O (1cm)), the latter being relevant for assessing the upward longwave, sensible and latent heat air-sea fluxes.Donlon et al. (2002) established a classification of near-surface ocean temperature measurements, insisting that measured temperatures are instrument-dependent, e.g. on the depth of the probe used for field measurements or on the wavelength used in radiometry.Ward et al. (2004) and Ward (2006) performed field measurements of the "skin" temperature layer, showing that assimilating it to the subsurface one could yield errors on the air-sea heat fluxes in the order of 10 to 50 W/m 2 .In parallel, additional parameterizations were included within bulk closures to account for such effects (e.g., Fairall et al., 1996;Zeng and Beljaars, 2005;Bellenger et al., 2017).
The wind stress dependency to surface currents has originally been neglected in bulk closures.However, both numerical experiments (Pacanowski, 1987;Duhaut and Straub, 2006;Dawe and Thompson, 2006;Renault et al., 2016) and flux measurements (Kelly et al., 2001) have shown that surface currents bear an impact on air-sea fluxes.More recently, the effects of this wind stress modulation on coupled ocean -atmosphere have been investigated (Renault et al., 2019a), and turned out to have a non-negligible impact on the energetics in coupled simulations (Renault et al., 2019b).
The parameterizations listed above are all part of a community effort to include the influence of ocean-specific processes within surface layer parameterizations.In this paper, our objective is to develop a general method for adding a simple representation of the ocean near-surface layer within existing bulk closures.The main idea is to extend a given one-sided bulk formulation to account for the evolution of currents and temperature, by extrapolating them from the depth at which the ocean information is available up to the surface, as shown in Fig. 1.At first, our approach is built in a very idealized framework.However, the formalism is flexible enough to seamlessly include additional or alternative parameterizations, such as the effects of waves of air-sea momentum exchanges and the potential wind stress deflection resulting from it.
The paper is organized as follows.Section 2 briefly introduces existing bulk closures and underlines their limitations in a coupled ocean-atmosphere context.Section 3 introduces a slight modification on atmosphere bulk closures, allowing them to explicitly treat the air-sea interface.In Sec. 4, an idealized parameterization for the ocean SL is introduced, solely accounting for shear turbulence.In Sec. 5, the inclusion of some specific non-turbulent phenomena within this framework (effects of waves and of radiation penetration) are discussed, thus illustrating its flexibility.
Section 6 investigates the effects of our novelty of offline turbulent flux computations.Finally, concluding remarks and perspectives are given in Sec. 7.

| EXISTING BULK CLOSURES AND THEIR LIMITATIONS IN A COUPLED CONTEXT
Throughout this paper, local horizontal homogeneity is assumed.Our study is therefore set in a 1D air-sea column configuration.z , the vertical coordinate, is orientated upwards, with the mean sea level assumed flat and defined by {z = 0}.Therefore, z < 0 in the ocean and z > 0 in the atmosphere.The index α ∈ {o, a } distinguishes ocean and atmosphere variables.The physical variables of interest are the horizontal velocity u ∈ and the potential temperature θ ∈ * + .The atmosphere is assumed to be dry, and latent heat resulting from mass exchanges is neglected.
This physically limiting assumption is made for easing the formulation of our framework.Moisture effects can be implicitly included by considering virtual potential temperatures, and latent heat can be explicitly included by treating the moisture profiles in the same manner as temperature ones.While u is assumed always aligned in the same direction, Sec.5.1 investigates the effects of 2D horizontal velocities on the potential deflection of surface stress from nearsurface winds.The letter x will be used as a general notation for either u or θ.Our focus is on the surface layer (SL), defined for each variable as the to the heights at which the information is extracted from the vertical finite-difference scheme of each model.These values of z 1 o and z 1 a are used in all numerical examples below.The SL is nested within the ocean-atmosphere boundary layer, being roughly 10 times thinner.In forced or coupled numerical simulations, it corresponds to the layer which is not covered by the vertical discretization of the considered model (ocean or atmosphere).Table 1 contains a nonexhaustive list of mathematical symbols introduced in this paper.

| Air-sea turbulent fluxes and their relationship to SL solution profiles
This section is intended for general readers to recall the basic aspects about the derivation of bulk formulations necessary for the proper understanding of our approach.The boundary conditions enforced to the atmosphere at z = z 1 a are: where K a,u is the momentum diffusivity; K a,θ is the thermal diffusivity; τ is the wind stress; ρ a is the atmosphere density; Q H is the sensible heat flux; c p a is the dry atmosphere heat capacity.In Equation 1, ρ a and c p a are assumed constant and known.τ and Q H are turbulent fluxes to be determined.Both diffusivities K a,x also depend on turbulent scales, and thus need to be parameterized.In the atmosphere SL, the relevant turbulent scales for velocity and temperature are u * a > 0 and θ * a , respectively.Formally, MOST states that the atmosphere SL contains a "purey turbulent" sublayer (grey shading in Fig. 1), above the direct influence of surface roughness and below the heights at which non-turbulent processes (e.g.pressure gradients, Coriolis effect) become important.In this layer of MOST validity, (u * a , θ * a ) can be linked to (∂ z u, ∂ z θ) by building dimensionless groups and applying the π-theorem (Buckingham, 1914).Throughout this manuscript, we assume z 1 a to be within the layer of MOST validity.Practically, this implies that Equation 1 holds on this part of the SL, with both diffusivities given by: where κ ≈ 0.4 is the von Kármán constant.L a = u * a 2 θ(z 1 a ) g κθ * a , where g ≈ 9.81 m/s 2 , is a signed characteristic length (Obukhov, 1971) rendering the fluid stratification's effect on K M O a,x through the two stability functions (φ m a , φ h a ) (e.g.Businger et al., 1971).From dimensional analysis, τ and Q H can be scaled as τ = ρ a (u * a ) 2 and Q H = ρ a c p a u * a θ * a .
From now on, we also assume that the SL is a constant flux layer, hence u * a and θ * a are constant as well.Combining the dimensionless groups with the constant flux assumption, by injecting Equation 2into Equation 1 and integrating with respect to z yields: where ) and (ψ m a , ψ h a ) are the integrated forms of the stability functions (Paulson, 1970).In the following, for the sake of conciseness, we assume that the atmosphere information is available at the same height z 1 a for both u and θ.Equation 3 most notably introduces (z r a,u , z r a,θ ), a set of two roughness heights, which are used as lower integration boundaries of the invariant dimensionless groups constituted by Equation 1. Indeed, physically, phenomena unrelated to MOST such as surface roughness are expected to dominate in the direct vicinity of the surface; mathematically, Equation 1 with K a,x given by Equation 2 is not integrable down to z = 0. Hence the introduction of roughness heights, which can be defined as the heights at which MOST-derived profiles (of u or θ) reach their respective surface values.In that regard, Equation 3 can be understood as an integration of idealized MOST profiles rather than the actual physical profiles, which are not known as z → 0. In Equation 3, the surface is assimilated to the roughness heights, so its left members are defined as . Equation 3 corresponds to classical "law of the wall" profiles (the logarithm term) perturbated by a stability-rendering term (the ψ x a term).Over the ocean, the stability at roughness heights can be neglected as z r a,x z 1 a , and Air-sea fluxes can be determined from Equation 3 by parameterizing z r a,u and z r a,θ as functions of u * a (e.g.Smith, 1988;Fairall et al., 2002).Equation 3 can then be exploited as a set of two nonlinear equations on (u * a , θ * a ).Solving it, usually through a carefully initialized fixed-point algorithm, leads to τ and Q H , as they are defined from u * a and θ * a .
This procedure is usually referred to as a bulk algorithm.Outside of a few exceptions (e.g. Louis, 1979;Dubrulle et al., 2002), in their vast majority, the theoretical basis of bulk closures is the MOST formalism described above, and their practical implementations arise from parameterizing the roughness heights and stability functions.

| Limitation on the use of classical bulk closures in a coupled context
Bulk closures described as in Sec.2.1 have been developed in the framework of atmosphere-only simulations, with ocean surface properties given as external forcings.Such simulations are usually carried with prescribed sea surface temperature (SST) and neglected sea surface currents since in most cases |u(z = 0) | |u(z 1 a ) |.Using bulk closures with these assumptions is consistent with field turbulent flux measurements, which are assessed at heights z ≈ z 1 a above the ocean, from "skin" surface temperature (at depth z ≈ −1 cm), unless an explicit parameterization is used (e.g. a cool skin parameterization such as Fairall et al., 1996).Since the closures include a representation of a vertical coordinate z 1 a , they can be used for both measuring turbulent fluxes (matching the measurement height), and for computing them in numerical models, with adapting z 1 a from context.As a consequence, when using transfer-coefficient based bulk closures (e.g.Large, 2006), atmosphere-model extracted values for u and θ are shifted to reference height levels, through MOST-derived invariant groups, in order to match parameterizations calibrated from observations.In other words, classical bulk closures as described in Sec.2.1 are expressed with a dependency on a reference height, so that the resulting fluxes are independent from it.In practice, this ensures that in the limits of MOST, air-sea fluxes resulting from classical bulk closures do not depend on the atmosphere model's vertical discretization.
To our understanding, directly extending (i.e.without near-surface parameterizations such as warm layers) classical bulk closures to a forced-ocean or coupled context yields inconsistencies, mostly related to the method (or lack thereof) used for incorporating near-surface ocean fields as bulk closure inputs.Unlike atmosphere-only simulations, in most coupled ones, the shallowest ocean informations available, usually located at a depth of z 1 o ≈ −1 m, is directly used as ocean surface information.This is equivalent to assuming that velocity (current) and temperature profiles within the ocean SL are constant with respect to z .Therefore, the depth at which the ocean information used as bulk input is taken does not have any influence on the formulation of the bulk closure.As a consequence, in forced-ocean or coupled experiments, turbulent fluxes arising from classical bulk closures are tributary to the vertical discretization of the ocean model.For example, carrying two "perfect ocean model" experiments with different near-surface vertical discretizations would yield distinct turbulent air-sea fluxes, which can be problematic.
In the following, we aim at building a formulational framework within which cross-medium bulk closures, more adapted to the context of coupled air-sea simulations, could be expressed.In particular, our formalism allows for vertical shifts to be performed within the ocean SL, so that the resulting bulk closures are depth-input dependent, in the same way classical ones are atmosphere height-input dependent.By design, the obtained air-sea fluxes would then be more robust and independent from the discretization of the ocean model.

| A SLIGHT ADAPTATION OF THE ATMOSPHERE BULK FORMULATION
In this section, we prepare our framework by revisiting atmosphere-only bulk closures.Below it is assumed that the surface currents are zero (i.e., u(0) = 0) and that the potential temperatures at z = z 1 a and z = 0 are known and used as inputs.Assuming surface properties to be known is an idealization, as such information cannot be measured (Donlon et al., 2002;Kent et al., 2017).A discussion on circumventing this issue is proposed in appendix B. Our objective here is to build a bulk closure which results from integrating solution profiles within the atmosphere surface layer (ASL) from the reference z 1 a height, down to the ocean -atmosphere interface (z = 0), instead of the traditional nonzero "roughness heights", so to get a direct connection to the underlying ocean.As in typical MOST applications, we assume the ASL to be a constant flux layer.Our bulk closure is thus based on integrating dimensionless groups with the following diffusivities, slightly adapted from Equation 2: where K v a,u , K v a,θ > 0 is a set of two constant diffusivities, which are to be determined.Their objective is to represent the molecular effects which dominate over the turbulent ones in the viscous sublayers, as z → 0. In Equation 4, they are divided by stability functions for convenience.The dimensionless groups arising from Equations 1 and 4 are: Mathematically, the inclusion of K v a,u , K v a,θ > 0 makes the dimensionless groups Equation 5 integrable on 0; z1 a .
In other words, while classical dimensionless groups are only relevant for describing the purely turbulent part of the ASL, Equation 5 can describe its entirety, down to z = 0. Using state-of-the-art stability functions (e.g., Högström, 1988), analytical integrated forms of Equations 1 and 4 can be obtained (see appendix A.1).They are however hard to manipulate.Obtaining simpler forms is possible by assuming: which is physically justified, as molecular effects are expected to be negligible compared to: (i) turbulent ones at z = z 1 a ; (ii) stability-induced ones.Assuming once again that the SL is a constant flux layer, integrating Equation 5 downwards from z 1 a and using Equation 6 yields (see appendix A.2 for proof): In particular, assessing Equation 7 at z = 0 leads to: where . Equation 8 can be used as a bulk closure, i.e. a set of equation on (u * a , θ * a ), depending on input fields at z = 0 and z = z 1 a .Practically speaking, K v a,u and K v a,θ can be parameterized so that the Equation 8closure is strictly equivalent 1 to the classical one Equation 3, using roughness heights as lower integration boundaries.
This can be done by setting: In other words, the new ASL closure Equation 8 is transparent in that it only differs from classical ones in terms of formalism, when assuming surface ocean fields to be known.More adequate tuning for K v a,x requires the ocean surface layer to be treated, and are discussed in appendix B. When perturbating diffusivities with positive constant factors as in Equation 4, the resulting viscous profiles' asymptotes are logarithmic functions, stopping at an equivalent roughness height of K v a,x / κu * a .It should be underlined that any other type of profiles (as long as they are monotonous with respect to z ) can be generated by relaxing Equation 4and perturbating K v a,x with a carefully built z -dependent function.However, this endeavor is not pursued here as our objective is to obtain a formalism that simply treats the full [0; z 1 a ] interval, rather than describing the viscous sublayers as accurately as possible.Hence, the viscous parameterization is kept simple, with a minimal (null here) impact on bulk closure outputs.
Using the adapted closure Equation 8 instead of the classical one Equation 3 has three assets.First, it is directly derived from integrating a dimensionless group down to z = 0, instead of using the roughness heights as an arbitrary lower integration bound.In our opinion, including K v a,x (and tuning it through Equation 9) within the effective surface layer diffusivities is more intuitive than stopping the dimensionless groups' integration at the z r a,x nonzero heights.Second, the new closure includes a smooth transition between the turbulent (z z r a,x ) and viscous (z z r a,x ) sublayers of the ASL.Integrating Equations 1 and 4 from z 1 a down to z ∈ 0; z 1 a allows for a full coverage of the ASL, including the viscous sublayers which cannot be represented by classical dimensionless groups.As Fig. 2 shows, using the modified bulk closure allows for smooth profiles to be integrated down to z = 0, and leads to negligible differences far from the viscous sublayers, as soon as z 10 −3 m.Third, at the expense of physically crude hypotheses on the viscous sublayers, it permits unambiguously expressing the solution profiles and their derivatives at z = 0 + , which will useful in the next section.

| IDEALIZED SYMMETRICAL OCEAN-ATMOSPHERE BULK FORMULATION
In this section, our objective is to extend the classical MOST framework so that two-sided ocean -atmosphere closures can be expressed within it.Strictly speaking, this should be distinguished from deriving new closures: below, we are not establishing new parameterizations for roughness and/or stability.We are simply proposing a method for vertically extending existing closures into the ocean.Hence, our novelties are more to be understood in terms of framework rather than closure in itself.In the following, we assume that the surface current velocity and potential temperature are known at a reference depth z 1 o ≈ −1 m instead of the surface.This setting is similar to that of forced orcean or coupled ocean -atmosphere simulations, where the shallowest ocean information available is located at a nonzero depth.Before including more realistic physics in Sec. 5, here we simply extend the formulations of Sec. 3 to build idealized two-sided ocean-atmosphere bulk closures, aiming at determining turbulent scales from u(z and θ(z 1 a ).We refer to the "ocean near-surface layer" (ONSL) as the thin layer located between z 1 o and z = 0, above the shallowest ocean vertical level.The ONSL is much thinner than the ocean boundary layer (OBL), whose depth can reach up to a few hundreds of meters.

| MOST-derived ocean surface layer
In this section, the ONSL is idealized and its description is based on the following assumptions (which will all be gradually relaxed in Sec.5): (iONSL-1) The ONSL is a constant flux layer.
(iONSL-2) As in the ASL, only turbulence and molecular effects play a role on a the ONSL, i.e. the wave effects and radiative fluxes are not explicitly represented.
(iONSL-3) The wind stress and sensible heat fluxes are conserved across the ocean-atmosphere interface.
The relevancy of such an idealized description has been validated by direct numerical simulations (Tsai et al., 2005), as well as both laboratory (Wu, 1975(Wu, , 1984;;Mcleish and Putland, 1975) and field (Churchill and Csanady, 1983;Csanady, 1984) experiments.However, under moderate to strong winds, other sources of air-sea exchanges (such as wave-induced stress) develop and interplay with purely turbulent and viscous effects.In particular, as soon as u(z 1 a ) 3 m/s, (iONSL-2) becomes physically invalid as waves start playing a crucial role in the momentum transfer to the ocean.Therefore, the assumptions above are expressed to build our framework, which will be extended in Sec. 5 to account for more realistic parameterizations.
As for z 1 a in Sec. 2, we also assume that z 1 o , the shallowest ocean level, is located within the domain of MOST validity.(iONSL-1) and (iONSL-2) thus lead to modelling the ONSL in the same way as the ASL was, through ocean dimensionless groups analogous to the atmosphere ones Equation 5, for z 1 o ≤ z ≤ 0: where The integrated forms of Equation 11b and Equation 11c are given by (for x ∈ {m, h }): where

| Turbulent closure for the idealized ONSL
In this section, we rely on the hypotheses above to constrain the four new unknown quantities introduced by Equation 10.First, (iONSL-3) yields: where Second, unlike classical bulk closures, the adapted ones allow assessing solution profiles at the interface z = 0. Across the interface, the gradients of the solution profiles are assumed to satisfy the following molecular constraint: where K m α are the kinematic viscosities for α ∈ {o, a } medium: K m a = 1.5 × 10 −5 m 2 /s and K m o = 10 −6 m 2 /s.Equation 14 implies that at the interface, the intrinsic properties of each medium determine the slope break of inbetween fluxes.Imposing Equation 14with Equations 5 and 10 injected yields: where µ m = K m o /K m a ≈ 6.7 × 10 −2 .Equations 13 and 15 introduce four new constraints which bind the four ocean turbulent and molecular quantities to their atmosphere counterparts.Yet, achieving turbulent closure cannot directly be done by transposing Equation 8, as in this section, surface currents and temperature are unknown.This can be overcome by integrating Equation 10 on z 1 o ; 0 , with assuming as in Equation 6that and by injecting Equations 9 and 15: Combining Equations 8 and 16, and assuming that u and θ are continuous at the interface (i.e., u(0 + ) = u(0 − ) and , which is only relevant with revised closures (as classical ones cannot treat z = 0), lead to: where the terms depending on either u * a or θ * a are in bold font, z r a,x are assessed from existing roughness parameterizations (e.g., Smith, 1988;Fairall et al., 2002).The λ x Unlike classical bulk closures, which are limited to [z r a,x ; z 1 a ], the revised atmosphere closure introduced in Sec. 3 permits unambiguously describing the surface layer arbitrarily close to the z = 0 interface.This asset has been used for enforcing the continuity of solution profiles, and the surface gradient condition Equation 14 at z = 0. Hence, although the revised atmosphere SL closure is transparent in terms of bulk outputs, it paved the way for obtaining the two-sided closure Equation 17.
Figure 3 represents solutions profiles derived from classical, revised one-sided (derived from Sec. 3) and cross-interface two-sided bulk closures.In this idealized case, cross-interface profiles are expectedly smooth, with sharper gradients very close to the ocean-atmosphere interface, and slope break at z = 0, as specified by Equation 14.This is physically relevant as z = 0 corresponds to a physical interface.While the ocean contribution to the surface layer variations of θ are barely noticeable (in Fig. 3, θ 0 ), those to the variations of u are more prevalent ( u 0 ).This can be explained by the fact that λ u ≈ 3.75 × λ θ .Shear turbulence rendering within the SL has a remarkably weaker effect on SST compared to diurnal heating (±3 K, e.g.Halpern and Reed, 1976;Stuart-Menteth et al., 2003) or even cool-skin effect (−0.2 K, e. g.Saunders, 1967;Fairall et al., 1996).However, twosided closures lead to u(z = 0) being closer to u(z 1 a ) than what one-sided closures would predict.Since the relevant large-scale shear for assessing turbulent fluxes is u z 1 a 0 , two-sided closures will then lead to distinct turbulent fluxes compared to one-sided ones.

| Impact on turbulent fluxes
By default, traditional bulk closures neglect ∂ z u and ∂ z θ in the ONSL.Unlike the revised atmosphere closure, where K v a,x had been tuned from z r a,x so that bulk outputs are unchanged, using two-sided cross-interface closures has an impact of the resulting turbulent scales (u * a , θ * a ), and thus on the air-sea turbulent fluxes.This is represented in Fig. 5, which shows that using our two-sided bulk versions leads to dampened turbulent fluxes.The effect is stronger on wind stress, which was to be expected, since Fig. 3 suggested that two-sided closures affect velocities more than potential temperatures.
The differences in fluxes displayed in Fig. 5 should be understood as the potential error made when using such closures with ocean inputs at nonzero depth.This does not correspond to an error in classical bulk closures.Since two-sided bulk closures depend on z 1 o , they harbor an extra degree of freedom compared to classical closures.Hence the results obtained from two-sided closures cannot be fully reproduced by one-sided closures, even through retuning.The results presented above are based on two-sided bulk closures which have been built as extensions of classical ones with the minimal hypotheses of Sec.4.2.In particular, the surface roughness of two-sided bulk closures has been extended from their one-sided counterpart.A longer-term perspective is to recalibrate surface roughness in the context of two-sided closures.This could be achieved from colocated air-sea turbulent flux measurements, relying on ocean inputs at nonzero depths, which is well beyond the scope of this paper.
Our idealized study has shown that accounting for shear turbulence within the ONSL may have a non-negligible impact on the representation of surface currents, and a very limited impact on surface temperature.These effects lead to perceivable changes on wind stress and sensible heat computations.In our idealized context, one way of transparently representing the ONSL is to rely on two-sided bulk closures, since they account for variations of currents and temperature due to shear-generated turbulence within the ONSL, and include a dependency to the depth from which the ocean information is extracted.In that regard, using two-sided closures is equivalent to extrapolating ocean currents and temperatures from z 1 o up to the surface, so that the fields considered as bulk formula inputs match what these formulations have been calibrated from.

| TOWARDS INCREMENTALLY MORE REALISTIC TWO-SIDED SURFACE LAYER
In this section, more elaborated SL physics, rendering processes other than turbulent shear, are included in our twosided framework developed in Sec. 4. The objective is to show the flexibility of our framework.Section 5.1 focuses on the wind deflection with strong currents; Sec.5.2 on representing the impact of surface waves; Sec.5.3 on including radiative fluxes within the ONSL.

| Velocity profile deflection under misaligned winds and currents
Throughout Sec. 4 near-surface winds and surface currents were assumed aligned with the i -axis and u z 1 a z 1 o was a scalar quantity.Relaxing this hypothesis can be carried out by representing 2D horizontal vectors as complex numbers, i.e.
Unlike other studies (e.g.Bressan and Constantin, 2019), here we assume the deflection between wind and near-surface currents to be known and given as an input.Our objective is then to investigate the velocity's rotation in the (i , j ) plane, focus on the direction of the surface currents and its influence on wind stress.
In a 2D context, assuming that the SL is a constant flux layer implies conservation of wind stress in both amplitude and direction.∂ z u is always co-aligned with τ in either side of the interface.Let us call ϕ τ ∈ [−π; π] this direction.
Integrating ∂ z u between any pair Equation 18 means that if ∂ z u is always aligned with τ, then so is the velocity shear between any z -increasing pair of vertical levels located within the SL, regardless of the chosen pair.In particular, Equation 18 implies that the stress directions obtained from relative-winds one-sided and two-sided closures are identical, as Arg u In other words, within a constant flux layer, the sampling heights of velocities (currents or winds) have no impact on the stress direction, as long as shear is considered.However, including the shear direction (whatever sampling heights it comes from) does have an impact on wind stress norm and direction, in comparison with bulk closures relying on absolute winds.Both velocity subgroups Equation 5a and Equation 10a can be rewritten as: for z ∈ 0; z 1 α and α ∈ {o, a }.Equation 19 can then be integrated to obtain a two-sided closure similar to Equation 17a, with the left member substituted by u  (Wu, 1983), where the impact of currents direction on |τ | is of the order of a few percents.

| Impact of waves on adapted bulk closures
In this section, the impact of waves within our adapted bulk closures is discussed.Under moderate to strong winds (u(z 1 a ) ≥ 3 m/s), waves develop on the sea surface and Langmuir turbulence is generated within the upper ocean.Both processes affect the air -sea momentum exchanges.As a consequence, a wave boundary layer (WBL) is generated, where the velocity profiles are dependent on both turbulence and wave-induced stresses.In this paper, our focus is on the air-sea surface layer, defined as roughly the top 1 m of the ocean and the bottom 10 m of the atmosphere.
While the atmosphere WBL is nested within our region of interest, as soon as significant waves develop, the ocean WBL spans outside our region of interest.Hence, an investigation on the effect of waves on our adapted bulk closures is needed.
Prior to carrying on, it should be clarified that the scientific question we want to address here is not as broad as that of the effects of waves on the full ocean boundary layer, whose depth is typically O (10 − 100 m), and the adequate parameterization for rendering them (see Esters et al., 2018;Li et al., 2019, for reviews).Our focus is on the extrapolation of ocean values from z 1 o ≈ −1 m up to the surface from a MOST-derived formalism.The scientific question we want to address is: assuming u(z 1 o ) known (which potentially includes Stokes drift contributions), how can the MOST formalism be further adapted to account for an ONSL perturbated by waves, and what is the subsequent impact on the wind-induced stress τ?In other words, we are focusing on the interplay between wave-generated momentum and shear-driven turbulence within the direct vicinity of the ocean surface.To answer this question, we first make a comment the implicit rendition of waves by existing parameterizations (Sec.5.2.1); then, we investigate the impact of wave-induced momentum stress on our adapted closures (Sec.5.2.2).

| On the default and implicit inclusion of wave effects within bulk closures
By default, all bulk closures include the effects of waves, at least to a minimal extent.Indeed, surface layer parameterizations used within bulk closures have been calibrated from field measurements, which may already partly incorporate the effects of waves.For example, the roughness height z r a,u is directly affected by the presence of waves; this can be tuned via the Charnock parameter (Kitaigorodskii, 1965).Its piecewise linear definition in the COARE bulk formula (Fairall et al., 2002) is one example: without prior knowledge of current wave state, it aims at representing the impact of wind-generated waves.Donelan (1982) Accounting for waves within SL parameterization can be further carried out by adapting the effective viscosities K o,x .
McWilliams and Sullivan (2000) proposed: where c w = 0.08.La = u * o /u st k is the Langmuir number (its typical range is 0.2 ≤ La 0.7), with u st k the Stokes drift at the surface.While more elaborate diagnoses have also been implemented, derived from Equation 20(e.g.Smyth et al., 2002) or resulting from large-eddy simulations (Large et al., 2019), in the following we will retain Equation 20as our framework aims at being applicable for climate models run at coarse resolutions, where Stokes drift is not necessarily available.u profiles in the ONSL with the diffusivities defined by Equation 20 are shown in Fig. 9.The middle panels of Fig. 9 show that the surface velocities derived from Equation 20("Langmuir", dotted lines) are slightly closer to u(z 1 a ) than the classical ones ("no wave", full lines).As a result, using Equation 20 comparatively decreases the shear u z 1 a 0 and both friction velocities u * α , in agreement with other results from large eddy simulations (e.g.et al., 1997) or simplified models (e.g.Teixeira, 2018).

| Effects of additional wave-induced ocean momentum input on bulk closures
In this section, we investigate the impact of an additional surface stress linked to the presence of waves, denoted τ w , considered here as an external momentum source term.We now consider that the effective stress in the ONSL is: In numerical simulations, τ w could be assessed by an external wave model and sent to the ocean as an additional boundary condition.To our understanding, the impact of injecting τ w on τ depends on the choice of bulk closure: • One-sided absolute-winds bulk closures: no impact.Currents are completely neglected, which implies that regardless of τ w and u(z ≤ 0), the resulting τ stress remains the same.
• One-sided relative-winds bulk closures: indirect impact.The integration of τ w in the ocean momentum boundary condition has an impact on u(z 1 o ), which is assumed equal to u(z = 0), which is itself used as bulk input, hence τ is indirectly impacted.
• Two-sided bulk closures: both direct and indirect impacts.In addition to the impact presented above, the ocean velocity dimensionless group Equation 10a is affected by τ w , as the relevant momentum flux describing the ONSL becomes τ + τ w .This in turn leads to a different velocity closure, hence an additional direct impact of waves on τ, which emerges because the ONSL is treated with a dimensionless momentum group.
In the following, the direct impact referred in the last point above is investigated.Using Equation 21implies that in the ONSL, ∂ z u is aligned with τ o,ef f , which can be misaligned with τ.As a consequence, wind stress is not necessarily aligned with u z 1 a z 1 o : it can be deflected when crossing the air-sea interface.Such phenomena have already been observed in the field (Geernaert, 1988;Geernaert et al., 1993;Grachev et al., 2003;Chen et al., 2018) and obtained from large-eddy simulations (Large et al., 2019;Patton et al., 2019).We define (ϕ a , ϕ o , ϕ w ) ∈ [−π; + π[ 3 the directions of the atmosphere, ocean and wave wind stresses.Since τ w is considered as a known source term, ϕ w is assumed constant and known.ϕ a and ϕ o are unknown constants, aligned with the potentially distinct shear directions in either the atmosphere and ocean part of the SL.Using Equation 21 instead of wind stress conservation yields substituting Equation 13a with: which is a two-fold constraint on both the norm and direction of the τ and τ o,ef f .As a result, the closure equation for momentum is also bidimensional and set in ¼, which has an impact on the bulk closure and its algorithmic implementation (see appendix C for more details).
Wind stresses τ obtained from bulk closures including τ w are shown in Figs. .The wave-dominated stress balance represented in Fig. 7d is analogous to low-wind conditions, when the wind stress is small compared to the swell-induced one, hence the momentum flux can be transferred upward (from the ocean to the atmosphere, e.g.Jiang et al., 2016;Högström et al., 2018), with negative drag transfer coefficient (Smedman et al., 1994).As previously mentioned, with our 2D framework derived from Equation 22, the effective drag coefficient is a complex number.Hence it can account for upward momentum transfer (when Arg τ/ u .However, Fig. 8b also shows that Arg τ, the wind stress norm, is also lightly impacted.This can also be perceived upon careful examination of the red arrows in Figs.7c and 7d.While the first effect described above is attributable to τ w increasing the effective stress transmitted to the ocean, this second effect derives from the impact of τ w on the τ closure itself.Indeed, by including τ w , the stress exerced on the ONSL is enhanced by τ w ; as a reaction to this, the ocean velocity dimensionless group is changed.Consequently, the bulk closure, which includes an integration of the velocity group in the ONSL (see Equation 10a), is altered.The turbulent closure adapts itself so that the resulting u * o,ef f (resp.u * a ) properly connects u(z 1 o ) (resp.u(z 1 a )) at the bottom (resp.top) of the surface layer, with ∂ z u being driven by τ + τ w (instead of τ) in the ONSL.As for the first effect, this second effect of τ w on momentum closure cannot be represented by one-sided closure, since they do not integrate velocity dimensionless groups on the ONSL.
Solution profiles arising from two-sided closures, with or without waves, are displayed on Fig. 9.The profiles without waves, or with effective diffusivities taking into account the Langmuir number Equation 20 are very close to each other, hence, the impact of Equation 20 is negligible in our framework, focusing on the ONSL.This does not mean that Equation 20 has a small impact on the whole ocean, as this parameterization is to be used on the whole OBL instead of, for example, a standard "K-profile" parameterization (Troen and Mahrt, 1986).Solution profiles with τ o,ef f = τ + τ w are also represented in colored dashed plots.Injecting τ w has an impact in both the atmosphere and the ocean SL, with the changes being more prevalent in the ocean, as expected from Equation 21.In Fig. 9, the near-surface winds and currents are assumed aligned, hence u j (z ) = 0 for the waveless (black lines) and Langmuir profiles (dotted black).Closer inspection of the u i profiles in the direct vicinity of the ocean surface (see Figs. 9c, 9i and 9o) illustrates the second effect described above on τ w perturbating the turbulent closure.If τ w is in the same direction as u here), then including τ w decreases u * a .Indeed, in that case, τ w contributes jointly with τ to make u(z 1 a ) connect u(z 1 o ).
Hence, u * a , which scales the shear-induced stress, is dampened, because the connection between u(z 1 a ) and u(z 1 o ) is also partly sustained by τ w .On the contrary, if τ w is opposing u z 1 a z 1 o (Arg τ w = π here), then including τ w increases u * a : the shear-induced stress needs to be stronger to connect u(z 1 o ) with u(z 1 a ), because it also has to counteract τ w in the ONSL.The intermediate cases in terms of direction (Arg τ w = π/4, π/2, 3π/4 here) cover the spectrum between both extreme cases presented above.In such intermediate cases, the solution profiles are deflected from the i -direction (hence u j 0), so that when the velocity groups are closed in the ONSL with (τ w ) j 0, u(z 1 o ), which is aligned in the i -direction, is properly reached, and the j -component of τ w is cancelled out.
The methods and results presented in this section attempt at representing wave impact while staying within the framework of MOST-derived bulk closure algorithms.It should however be reminded that since MOST does not hold in the presence of a significant WBL, an accurate and more legitimate representation of atmosphere -wind -ocean coupling cannot be formulated in this framework.Coupled wave boundary layer models are the designated tool for tackling this problem (see Chalikov and Rainchik, 2011, for a review).

| Radiative flux including ONSL parameterization
Radiative fluxes have been neglected from classical atmosphere-only bulk closures, since they are assumed independent of z in the lowest few meters of the atmosphere (see footnote 2 in Monin and Obukhov, 1954), and the radiative budget of the atmosphere SL is in equilibrium.Both these hypotheses are reasonably accurate in the context of atmosphere-only closures.However, in our two-sided framework, accounting for radiative fluxes (and thus lifting the (iONSL-2) hypothesis) is required for two reasons.First, the radiative budget on the full SL is not in equilibrium: the net radiative flux at the air-sea interface is not zero (see Sect. 5.3.1).Second, shortwave radiative fluxes can display perceivable vertical gradients over the ONSL (see Sect. 5.3.2).
Let us call Q 0 l w the net longwave flux at the ocean surface (i.e., accounting for surface blackbody radiation) and Q 0 sw the net surface solar radiation (i.e., accounting for surface albedo), both fluxes being positive downwards.The boundary condition on θ at the ocean surface (z = 0 − ) now reads: where ) can be derived from Equation 10b.The minus sign in front of radiative fluxes is due to the z -axis being orientated upwards.It should also be mentioned that models typically do not inject Q 0 sw as a boundary condition and instead prescribe it as an additional volumetric heat source term in the few first vertical levels (Jerlov, 1976).Since here we only focus on the SL, which is not treated by models, we will retain the shortwave-including form Equation 23.Theoretically, Q 0 l w should be unknown since it includes an upward flux ∝ (θ(0)) 4 and only θ(z 1 o ) is known.We will however assume Q 0 l w to be known accurately enough, since θ(z = 0) ≈ 290 K and previous results have shown that within our hypotheses,

| Constant ONSL radiative flux approximation
As a first step, let us consider the simplified case where both radiative fluxes are deemed constant on the ONSL, i.e.
Equation 23 is valid for z 1 o ≤ z ≤ 0. In that case, the ONSL is still a constant flux layer, which can be described where the first term arises from sensible heat conservation through the ocean -atmosphere interface.Since here, Equation 23 is assumed true on z 1 o ; 0 , integrating it on this interval with Equation 24 injected leads to: which is a radiation-including version of Equation 17b, in the idealized constant flux case, and can therefore be used jointly with Equation 17a as a bulk closure.

| Ocean surface layer with depth-varying solar flux
Longwave radiation is absorbed in (and emitted from) the first few millimeters of the ocean, hence we limit ourselves to directly injecting it as a boundary condition, and consider that it does not play a significant role on the dimensionless groups defined on ]z 1 o ; 0[.On the other hand, the shortwave (solar) flux can display perceivable gradients on the [z 1 o ; 0] interval, with z 1 o ≈ −1 m.In the upper ocean, its penetration can be parameterized as a combination of exponential modes depending on various factors such as incident angle and ocean biochemistry (e.g., Soloviev and Vershinsky, 1982;Morel and Antoine, 1994;Ohlmann and Siegel, 2000).Here we use a parameterization established by Paulson and Simpson (1981): where k i > 0 m −1 , 1 ≤ i ≤ 9 characterize the typical damping depths, and A i their relative intensities ( A i = 1).
Values of A i and k i can be extracted from Table 1 of Paulson and Simpson (1981).Including Equation 26 in the ONSL parameterization breaks the constant flux layer hypothesis (iONSL-1) and thus requires further adaptation of the bulk closure, compared to Sec. 5.3.1.Rendering depth-varying solar fluxes can be carried out by relaxing θ * ,0 o,r ad and letting it vary with z : In particular, θ * o,r ad (z = 0 − ) = θ * ,0 o,r ad .Integrating Equation 10b on z 1 o ; 0 with θ * o substituted by θ * ,0 o,r ad from Equation 27 (this can no longer be considered an as "invariant group", since the flux depends on z ), injecting Equation 26and rearranging terms yields the following two-sided closure on θ: Under stable stratification, Figs.10d to 10f suggest that similar conclusions hold, with θ profiles decreasing with z in the ONSL during daytime (see Fig. 10f), in contrast with shear-driven θ profile are then increasing with z .While the effects are expectedly more perceivable with constant fluxes than with fluxes progressively dampened with depth, the explanations given above hold in both cases.

| AN OFFLINE NUMERICAL ASSESSMENT OF TWO-SIDED CLOSURES
In both coupled or forced-ocean simulations, the implementation of two-sided SL parameterizations can swiftly be carried out by patching the existing bulk fixed-point algorithm to make it solve Equation 17 instead of Equation 3.Here we investigate the impact of our novelties by reassessing and comparing turbulent fluxes on offline computations, using ocean and atmosphere reanalyses as input data.For doing so, we compare fluxes obtained from global largescale ocean and atmosphere reanalyses.The reference bulk formula is COARE 3.0 (Fairall et al., 2002) as per the aerobulk package (Brodeau et al., 2017).We perform standard absolute-winds, relative-winds and two-sided bulk closures for the year 2006, which has been close to the current era climatology.Ocean inputs are given by the GLORYS2V4 dataset (Mercator Ocean, 2019), extracted at z 1 o ≈ −0.5 m, and atmosphere forcings are given by ERA-Interim (Dee et al., 2011), extracted from z 1 a = 10 m.Since the two-sided bulk formula has only been introduced for the open ocean, all grid cells where the sea ice concentration exceeds 15% have been excluded from computations.
While atmosphere moisture is included, in order to fit with our hypothesis, turbulent-scale moisture effects have been screened out by assuming q (z = 0) = q (z 1 a ), thus resulting in null latent heat flux.Numerical results from these offline tests are shown in Figs.11 and 12.
Fig. 11 shows the yearly mean discrepancy of turbulent fluxes when using two-sided bulk closures in comparison with standard relative-winds ones.A negative bias is observed on |τ |, which was to be expected since including previouslyneglected layers leads to velocity shear damping.The positive bias in Q H is explained by the fact that on average, stratification is unstable, i.e., Q H is "less negative" with two-sided closures than with standard ones.The larger mitigation values are reached where ocean currents are strong, in the Antarctic circumpolar current, the Kuroshio current and the Gulf stream.Fig. 12 shows 2006 daily time series of turbulent fluxes and surface fields at one location in the Kuroshio current ((149.25 • E, 36.75 • S), marked in Fig. 11).Fig. 12 suggests that punctually, turbulent flux mitigation from using two-sided closures can significantly alter the resulting fluxes.The biases between classical closures and two-sided ones are marked with a very high temporal variability, which was to be expected: our framework assumes that the ONSL directly responds to surface forcing from the atmosphere, hence the strong atmosphere variability is transferred into the upper layer of the ocean.Figure 12a suggests that while including relative winds yields stronger wind stress, using our two-sided lightly weakens it.This can be explained by the shear, which is typically less important in our two-sided framework, since across the ONSL, the currents will progressively adapt to the near-surface winds (see Fig. 12c).This observation is coherent with results presented in Sec.5.1.
The impact on Q H can be at times quite large, although it is negligible most of the year, for both relative-winds and twosided closures (see Fig. 12b).On Fig. 11, sensible heat mitigation is relatively low at high latitudes because grid cells covered in sea ice have been screened out.As already pointed out in Sec.4.1, the effects of our two-sided closures on the SST are negligible on average (in the order of 0.05 K).However, Fig. 12d features a few extreme events where two-sided closures significantly cool down the SST (early August, late September).These events are concomitant with low-wind conditions under unstable stratification, which is consistent with results already presented in Fig. 5b.The relatively weak reaction of SSTs can be explained by a combination of three factors.First, temperature (or equivalently, heat) mitigation is globally weaker than the velocity (stress) one since λ θ λ u .Second, this idealized test does not include radiative fluxes and their subsequent effects, which would generate stronger SST variability.Third, in this online test, the input temperature are not allowed to drift in time: only their instant reaction to SL parameterization is shown.

| DISCUSSION
We have introduced a formalism for extending air-sea turbulent flux parameterizations in order to make them account for shear turbulence driven effects on both sides of the interface, including the ocean.Special care has been taken to ensure by design the smoothness of solution profiles within the SL.The impact of our novelties in the SL treatment has been investigated on both idealized and more realistic cases.In general, our findings affect near-surface velocity profiles more than temperature ones.Occurences with significant differences on turbulent air-sea flux determination (up to 20%) have been underlined.The effects are mostly concentrated on the representation of surface velocity, which then impacts the wind stress and sensible heat fluxes through their dependency on z r α,u and u * α .Such results may have implications for describing debris transport in the upper ocean.Recent results on that field suggest strong variations in the top meter of the ocean, albeit linearly varying with z (e.g.Laxague et al., 2018), instead of the logarithmic profiles used in Sec. 4.An interesting perspective would then to build sturdier dimensionless groups for velocity profiles within the ONSL to match the results found in such studies.
The main message of this paper is that typical bulk closures have not been designed for being directly used in forcedocean or coupled ocean-atmosphere settings, when ocean surface properties are unknown.One central inconsistency we have lifted supporting this idea is that classical bulk closures do not depend on the ocean vertical discretization.
The two-sided formalism introduced above is adapted for using state-of-the-art bulk closures with nonzero depth ocean information as input.We do admit relying on a crude surface layer representation and first neglecting physically determining phenomena, such as wave-induced enhanced momentum transfer or radiative penetration within the ocean SL.Yet, classical bulk closures fully neglect this part of the ONSL, hence, they implicitly rely on even cruder assumptions.Stripping down the SL parameterization to the simpler, mathematically more ergonomic formalism that we have relied on was a necessary step for developing our framework.Moreover, historical bulk closures have first been developed within a similarly idealized setting (i.e., shear-driven SL).Our approach is to propose a relevant framework within which incorporating new parameterizations could be carried out, without altering the global consistency of the SL scheme.Section 5 proposed a few examples of such extensions.We stand by the idea that explicitly parameterizing the ONSL, albeit in a crude way, is more legitimate than implicitly neglecting it.Indeed, since the ONSL is usually assumed passive in classical closures, its impact of the surface physics is "hidden".Consequently, we believe that explicitly acknowledging the ONSL, by formulating bulk closures including it, may attract attention towards developing physically more realistic two-sided closures.
Four specific further development perspectives retain our attention.First, adapting our framework to two-sided wave-permitting boundary layers effects on turbulent fluxes.Results obtained using simplified wave formulations are briefly established and discussed in Sec.5.2.However, our formalism cannot be used per se for representing the ocean surface at the viscous sublayers under conditions of strong winds, where the problem geometry is changed and waveinduced micro instabilities overshadow viscous stress.Accurately simulating fluxes under heterogeneous surfaces, such as wave-deformated oceans, has been a considerable research challenge for decades, even from the broader perspective of boundary layer meteorology (LeMone et al., 2019).Literature on this matter includes both simplified models (Troitskaya and Rybushkina, 2008) and three-component ocean -wave -atmosphere coupling (Hristov et al., 2003;Chen et al., 2007Chen et al., , 2013)), which is well beyond the scope of this paper.Second, including moisture and salinity influence on our two-sided algorithms ought to be carried out, since such effects are already present in one-sided turbulent closures.This could be done by adapting the study of Bellenger et al. (2017), which proposes an extension of the Zeng and Beljaars (2005) warm layer model for enhancing the representation of saline stratification in the upper ocean layer.Third, the Equation 14 molecular constraint imposed at the z = 0 is a serious simplification compared to the very rich physics of the viscous sublayers.In particular, the stress ratio at the interface has already been investigated in the literature, with different values or dependencies underlined (Saunders, 1967;Robinson et al., 1984;Ward and Donelan, 2006).Since no clear consensus has arisen from the aforementioned studies, Equation 14 has been used as a minimal representation of surface constraint.However, our framework could swiftly incorporate any explicit parameterization by reformulating Equation 14 and integrating it within the two-sided closure.More generally speaking, since our new formalism includes new physics, it also calls for calibrating anew bulk closures (roughness and stability representations) from two-sided turbulent measurements.This could potentially limit the spread between observations and parameterizations.Fourth, our study assumes that the ONSL immediately responds to above-surface fast changes, which is usually not the case as the ocean kinematic viscosity is ≈ 30 times greater than the atmosphere's.Soloviev et al. (2001) suggests rendering the nonstationarity of the ONSL by using a gradient Richardson number and linking it to the Obukhov length.Such an endeavor is beyond the scope of the current paper, where the stationarity of the ONSL is used as a working hypothesis.
On the longer-term perspective, we believe that parameterizing the SL in full at a continuous level clarifies the mathematical nature of the boundary conditions enforced between the ocean and the atmosphere.Due to their explicit form (see Equation 1), classical air-sea boundary conditions are subject to being erroneously assimilated to Neumann conditions.We however argue that the turbulent air-sea boundary conditions are, at the continuous level, equivalent to a combination of Dirichlet and Neumann boundary conditions.If we consider that the surface layer solution profiles obtained from the parameterization scheme are correct, then imposing Equation 1 is mathematically equivalent to imposing the continuity of the solutions and a constraint on their gradients at z = 0.In our opinion, transcripting the air-sea coupling problem into such a simpler yet sturdier mathematical formalism would ease the further theoretical development of turbulence-including coupling algorithms.At the practical, discrete level, connecting the surface layer with computational domains could be implemented by using specifically designed splines built from the chosen parameterization set.
which directly yields: Doing a zeroth-order approximation of Equation 31 in the K v a,u /(κu * a ) z 1 a , |L a | limit yields a form compatible with Equation 7. In the unstable case (L a , ζ < 0), integrating Equation 5a with Equation 29b injected and using the z → η = (1 − 16 z La ) 1/4 change of variable yields: where ν 0 = −κu * a L a /16 > 0. Depending on the sign of ξ = 1 − K v a,u /ν 0 , integrating Equation 32 yields: z 1 a , |L a | in Equation 33 yields ξ ≈ 1, and thus only the ξ > 0 case of Equation 33 is relevant.
Assessing it at z = 0 yields:

A.2 | Asymptotic development on the molecular effect including ASL closure
In this section, we prove that integrating Equation 5 in the K v a,u , K v a,θ κu * a z 1 a , κu * a |L a | limit leads to Equation 7.
For doing so, we go back to the φ m a , φ h a stability functions being generic, and assume that they satisfy the following constraints: (h1) φ m a , φ h a are smooth over (i.e., they are continuous); Note that both these constraints are satisfied by most classical stability functions.Integrating either subequations of Equation 5 leads to: where x ∈ {u, θ } and x v a ∈ {K v a,u , K v a,θ }, chosen accordingly.Performing the change of variable z → z = z + x v a /(κu * a ) and rearranging terms in Equation 35 leads to: Since x v a κu * a z 1 a , the numerator in the logarithm of Equation 36 can be reduced to z 1 a .For the second, stabilityrendering term of Equation 36, we argue that it is equivalent to −ψ x a (z 1 a ) + ψ x a (z ).Indeed: 1, and (h2) guarantees that the lower integration boundary can be approximated as z = 0; B | ADAPTING TWO-SIDED BULK CLOSURES TO TURBULENT FLUX MEA-

SUREMENTS CALIBRATED WITH NONZERO DEPTH OCEAN FIELDS
In Sec. 3, an estimate of K v a,x directly deduced from z r a,x has been given (see Equation 9).This corresponds to our understanding of the roughness heights in classical bulk closures, which match the heights at which MOST-derived solution profiles reach their expected "surface values".As a consequence, the adapted atmosphere bulk closure evaluated from the velocity and temperature jumps between z = 0 and z 1 a are equivalent to classical bulk closures.However, ocean surface properties cannot be measured: only near-surface properties can, at a depth of a few millimeters (Donlon et al., 2002) at least.Below K v a,x , x ∈ {u, θ } are evaluated anew so that the resulting two-sided bulk closure matches the velocity and temperature jumps from an arbitrary reference depth z ms o ≈ −1 mm to z 1 a .As a consequence, the computed fluxes will match the experimental setting measurements z r a,u and z r a,θ have actually been tuned for.The method described below can also be generalized for cases in which the reference currents and temperatures measurement depths are distinct.
Integrating Equation 10 on z ms o ; 0 , with Equations 13 and 15 injected, leads to, for x ∈ {u, θ }: where stratification has been neglected, since Summing Equations 8 and 37, rearranging terms, and identifying the resulting denominator to z r a,x , yields, for x ∈ {u, θ }: which is a condition for two-sided bulk closures to match classical bulk closures, assuming they have been tuned from ocean measurements at z = z ms o .In order words, evaluating z r a,x from a given classical bulk formula, and then solving where: A light adaptation of fixed-point bulk algorithms is needed in order to make them solve Equation 40 instead of Equation 17a.Indeed, non-wave including bulk algorithms rely on iterating over u * a with the following procedure: where k ∈ denotes bulk algorithm iterations.Solving Equation 40 can no longer be done by simply injecting Equation 42 as in the presence of τ w , u * o,ef f is not proportional to u * a anymore.Instead, we propose the following four-step procedure, deduced from Equations 22 and 40: where G a > 0 is implicitly assumed, which is reasonable since |z 1 α | |z r α,u | and ψ m a ≤ 1 over the ocean.F I G U R E 2 Solution profiles for (a) u and (b) θ in the ASL, under unstable stratification, arising from classical dimensionless groups (full plots) and the modified ones (dashed plots) relying on Equations 7 and 9.The roughness height parameterizations and the stability functions are taken from the COARE bulk formula (Fairall et al., 2002).b) and (e) represent solution profiles in the direct vicinity of the ocean-atmosphere interface (note the linear, zoomed-in z axis).Here, as in Fig. 2, the roughness parameterizations and stability functions are taken from the COARE3.0algorithm (Fairall et al., 2002) o > 0 ocean molecular viscosity and thermal diffusivity; φ m o , φ h o a set of two stability functions; and L o the ocean Obukhov length defined as in Large (1998), i.e.L o = u * o 2 κg α eos θ * o , with α eos ≈ 1, 8 × 10 −4 K −1 the ocean thermal expansion coefficient.The φ x o stability functions can differ from their atmosphere counterparts.In the following, we will use the ocean stability functions from Large et al. (2019):

Fig
Fig.4suggests that by using classical bulk closures in coupled simulations, the u and θ differentials considered as bulk closure inputs are systematically slightly overestimated.

Figure 6
Figure 6 illustrates the impact of shear direction in an extreme case, under neutral stratification ( θ z 1 a z 1 o ;Geernaert et al. (1987);Johnson et al. (1998);Taylor and Yelland (2001);Oost et al. (2002);Drennan et al. (2003) all propose more sophisticated examples of such parameterizations, with the Charnock parameter usually depending on the wave age c w p /u * a , where c w p is the wave phase speed at its peak frequency.
5 m/s, aligned with the i direction (Arg = 0).Figures 7 and 8 also feature τ 0 , the wind stress that would be obtained in the absence of external wave stress (τ w = 0), i.e. using the two-sided bulk closures discussed in Sec.4.3 (black lines, same on all subfigures).With the parameters described above, τ 0 ≈ 0.23 N/m 2 , aligned in the i -direction (Arg τ 0 = 0).On Fig. 7, each subfigure corresponds to a different value for |τ w | (green lines), and covers the full interval −π < Arg τ w ≤ +π; on the other hand, Fig. 8 represents τ depending on |τ w | for a limited set of values for Arg τ w (legend box).A first effect, observed in Fig.7only, arises from comparing τ o,ef f (blue lines) and τ 0 (black lines), the stresses transmitted to the ocean with and without the external τ w wave stress.Expectedly, including τ w drastically impacts τ o,ef f (blue lines) in both norm and direction.As |τ w | increases, τ o,ef f is progressively deviated from τ 0 towards τ w .This can be easily inferred from Equation21: τ w is an external ocean stress contribution, in addition to the wind-stress τ (red lines).For example, for Arg τ w = π (wave stress opposite to shear): in Fig.7c(|τ w | ≈ |τ 0 |), the wind and wave stresses counterbalance each other, so that in τ o,ef f is nearly zero.In Fig.7d(|τ w | ≈ 3|τ 0 |), the wave stress dominates over the wind stress, and thus the resulting stress transmitted to the ocean can be contrary to the direction u but also any stress deflection due to a given external wave-induced stress τ w , which can be fully independent from near-surface velocities (e.g., in the presence of swell).Since classical bulk closures neglect the velocity profile evolution in the ONSL, they imply C D ∈ * + , and thus can only transfer momentum in the direction of u z 1 a 0 (resp.u(z 1 a )) for relative-winds (resp.absolute-winds) closures.In that regard, our ONSL-including formalism allows more flexibility for representing wave-induced deflections of the wind stress, as their mathematical structure allows Arg τ o,ef f to be decorrelated , observed in both Figs.7 and 8, arises from comparing τ (red lines) and τ 0 (black lines), the wind stresses with and without τ w .This effect is less obvious, but it becomes more and more perceivable with increasing |τ w |: the presence of an external wave-induced stress also perturbates τ, the wind stress resulting from the bulk closure.As Figs. 7c, 7d and 8a suggest, the impact becomes prevalent on |τ | with strong |τ w | in opposite by an invariant group similar to Equation 10b.The right member of Equation 23 can thus be described as u * o θ * ,0 o,r ad with a constant θ * ,0 o,r ad scale to be determined.The presence of radiative fluxes in Equation 23 yields substituting Equation 13b with: = µ m z r a,u /λ u , e = exp(1) and E 1 (x ) = ∫ ∞ x e −t /t dt is the exponential integral with index 1.Expectedly, Equation 28 encompasses Equation 25; its last integral is a radiative-inclusive variant of the ψ h o integrated stability function, and cannot be determined analytically, in general.Since the radiative contributions of Equations 25 and 28 are part of ONSL-specific terms of the closure, they can only be rendered by two-sided closures.Indeed, default classical closures (cool-skin including ones being a notable exception) neglect all processes potentially impacting ∂ z θ in the ONSL, including radiation.θ solution profiles arising from different bulk closures, including radiative-inclusive ones, are illustrated in Fig. 10.As expected, including radiative fluxes has an impact of the θ solution profiles.Under unstable stratification Figs.10a to 10c suggest that during daytime (e.g. at Q 0 sw = 300 W/m 2 ), compared to radiation-neglecting bulks, the perceived surface temperature can be diminished by a few 0.1 K.This behavior is due to the bulk closure aiming at connecting the same θ(z 1 a ) and θ(z 1 o ) couple: if a positive radiative flux is enforced, then the bulk closure leads to a slightly colder surface temperature, to compensate for the additional heating in the ONSL, and still connect θ(z 1 o ).It should be clear that in Fig. 10, θ(z 1 o ) is taken as input, hence the results discussed only hold for the evolution of θ in the do not mean that the ocean surface is cooled down by incoming shortwave radiation.In nonstationary simulations, the solar flux would warm θ(z 1 o ) through its impact on the θ bounday condition, and dominate over the relative cooling observed in Figs.10a to 10c.At nighttime (Q 0 sw = 0 W/m 2 ), the total net radiation flux is negative since Q 0 l w < 0, and the solution profiles are increasing with z (see Fig. 10c): the surface temperature is slightly warmer than θ(z 1 o ), since the ONSL θ profile is dominated by a longwave-induced flux, cooling it down anew to θ(z 1 o ).
b) Two-sided scheme.F I G U R E 1 Surface layer parameterized schemes for a given transverse variable x ∈ {u, θ } (velocity or temperature) with arbitrary and nonuniform vertical scales.Grey shadings indicate the layers of MOST validity.(a) is the standard methodology: only a subset of the atmosphere surface layer, the Monin-Obukhov atmosphere surface layer (MO-ASL), is parameterized.Profiles are assumed constant on the atmosphere viscous sublayer (AVSL) and the ocean near-surface layer (ONSL), leading to gradient discontinuities at z = z r a,x and z = z 1 o .(b) is a two-sided parameterization scheme: the full surface layer is parameterized and thus the solution profile is mathematically regular, except at z = 0 where the solution gradient can be discontinuous.
Ocean near-surface layer contribution to full surface layer jump for (a) u and (b) θ arising from the idealized two-sided closure Equation 17. Mind the distinct x -axes and legend boxes.Here z 1 o = −1 m and z 1 a = 10 m.Turbulent fluxes ((a): wind stress, (b): sensible heat) arising from classical and two-sided bulk closures (solid and dashed lines, respectively).The x -axis is the wind-current surface layer shear u regima.Here the bulk formula used is COARE without gustiness nor warm layers parameterizations, with z 1 o = −1 m and z 1 a = 10 m.Influence of shear direction on surface layer properties.Here a two-sided 2D version of the COARE bulk formula is used with |u(z 1 a ) | = 2 m/s, Arg u(z 1 a ) = 0 (i.e.alignment with the i -direction), |u(z 1 o ) | = 0.5 m/s, and Arg u(z 1 o ) varying over [0; π].(a) displays: |τ | (black line) and Arg τ (colored arrows) for five specific values of Arg u(z 1 o ).With the same color code, (b) and (c) display (u i , u j ) in the ASL for the same five values of Arg u(z 1 o) (lines), with (u i (z 1 o ), u j (z 1 o ))specified by the thin vertical lines.Impact of an external wave-induced stress τ w under various wave directions (angle axis) and norm (radius axis and specified |τ w | values) on the wind stress.Plotted here are: wind stresses in the absence of wave (τ 0 , black); wave-induced stress (τ w , green); wind stress with a τ w -including closure (τ, red); ocean stress including wave effect (τ o,ef f blue).The norms are represented as lines (in polar coordinates) and the arrows indicate directions.Note that the radii of the polar axis and the scales of the arrow lengths vary inbetween subfigures.Obtained using a wave-including adapted version of the COARE bulk formula as described in appendix C. Here,z 1 o = −1 m, z 1 a = 10 m, u(z 1 a ) = 12 m/s, u(z 1 o ) = 0.5 m/s, Arg u(z 1 a ) = Arg u(z 1 o ) = 0, θ(z 1 a ) = 293 K, θ(z 1 o ) =295 K. Colored lines: wind stress τ ((a): norm, (b) angle) obtained from a wave-including bulk closure, with varying τ w (norm: x -axis; angle: color code).The black line represents |τ 0 |, the wind stress without τ w .The green dashed line represents |τ w |, for comparison with |τ |.Here, as in Fig. 7, u(z 1 a ) = 12 m/s, u(z 1 o ) = 0.5 m/s, Arg u(z 1 a ) = Arg u(z 1 o ) = 0. Note that for readability, the y -axis of (b) is defined in degrees (radians are the angle unit elsewhere).
Yearly 2006 mean of the differences on turbulent fluxes (a) τ and (b) Q H between two-sided bulk closures and classical ones, with relative winds (accounting for surface currents).The light blue dot locates the grid cell of the time series shown in Fig. 12. Symbol Description Unit α Atmosphere (α = a) or ocean (α = o) u Winds (α = a) or currents (α = o) Momentum roughness height (α = a) or depth (α = o) m z r α,θ Temperature roughness height (α = a) or depth (α = o) m Non-exhaustive list of symbols used in this paper.Symbols with the α subscript are twofolds: they are defined in both the atmosphere (α = a) and the ocean (α = o), and are distinct from one medium to the other.
[• • • ] terms in Equation 17, which encapsulate the novelties of two-sided closures, are neglected in the standard ones.Equation 17 can be used as a cross-interface turbulent closure: in other words, it provides a set of two nonlinear equations on u * a and θ * a for determining them from four large-scale, near-surface quantities: