<figure><img src="image_1.png" /></figure>

HAL
open science

# Formulation of a three-dimensional atmospheric boundary layer model for an improved representation of air-sea interactions in eddying oceanic models

Florian Lemarié

► To cite this version:

Florian Lemarié. Formulation of a three-dimensional atmospheric boundary layer model for an improved representation of air-sea interactions in eddying oceanic models. 2024. hal-04563051

HAL Id: hal-04563051

https://inria.hal.science/hal-04563051

Preprint submitted on 29 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Formulation of a three-dimensional atmospheric boundary layer model for an improved representation of air-sea interactions in eddying oceanic models

Florian Lemarié

► To cite this version:

Florian Lemarié. Formulation of a three-dimensional atmospheric boundary layer model for an improved representation of air-sea interactions in eddying oceanic models. 2024. hal-04563051

HAL Id: hal-04563051

https://inria.hal.science/hal-04563051v1

Preprint submitted on 29 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Formulation of a three-dimensional atmospheric boundary layer model for an improved representation of air-sea interactions in eddying oceanic models

F. Lemarié¹

¹Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France

**Abstract**

In this document, we derive a three-dimensional extension of the single-column atmospheric boundary layer (ABL1d) model derived in Lemarié et al. (2021). The idea is to formulate the ABL3d model in terms of perturbation around a large-scale ambient state obtained from analyses or reanalyses. The proposed formulation is both simple and efficient from the computational point of view. Using numerical experiments and comparison with LES simulations, we illustrate the improvement brought by ABL3d compared with ABL1d. The ABL1d approach was designed to represent the modulation of the atmospheric boundary layer by sea-surface anomalies. The ABL3d model has the same ability, but it can also represent horizontal advection and the adjustment of the atmospheric horizontal pressure gradient to sea surface temperature anomalies.

## 1 Introduction

To improve the consistency between the scales resolved by oceanic models and the scales present in surface fluxes at the air-sea interface (often calculated using low-resolution atmospheric data), an aggregation (a.k.a. downscaling) of atmospheric data available at low resolution to the resolution of the oceanic model has been developed in Lemarié et al. (2021). In their work, these authors have designed a simplified atmospheric boundary layer model via a so-called physically based low-fidelity approach. The objective was to find an approach alternative to usual forcing strategies that makes the atmospheric boundary layer interactive while maintaining synoptic forcing (see Lemarié et al. (2021) for a detailed description of the motivations). For this approach, the focus was on the representation of the modulation of atmospheric turbulence by sea surface temperature and current anomalies. The atmospheric boundary layer model (referred to as ABL1d) thus defined was implemented in a general circulation model and was tested on a series of idealized and realistic benchmarks. The objective here is to investigate the possibility of relaxing the single-column assumption behind the ABL1d model as a first step toward a more realistic ABL3d model which would be able to represent horizontal advection and the adjustment of the atmospheric horizontal pressure gradient to the underlying sea surface temperature.

Our approach throughout this document is to consider the problem as a multiscale problem representing the interaction between a large-scale, a priori known, state (provided by analysis or reanalysis data) and the so-called micro-scales (i.e. the scales of the high-resolution oceanic model). In Sec. 2 different ways to represent such multiscale interactions are presented and illustrated via a simple tracer equation with source term. In Sec. 3 the full atmospheric primitive equations are considered and a perturbation method is applied as a technique to derive a simplified three-dimensional ABL model. The resulting model is then discretized and applied to an idealized experiment for strong winds crossing a persistent mid-latitude sea surface temperature (SST) gradient from cold to warm SST (Sec. 4.2) and from warm to cold SST (Sec. 4.3).

Corresponding author: Florian Lemarié, florian.lemarie@inria.fr
# Formulation of a three-dimensional atmospheric boundary layer model for an improved representation of air-sea interactions in eddying oceanic models

F. Lemarié¹

¹Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France

## Abstract

In this document, we derive a three-dimensional extension of the single-column atmospheric boundary layer (ABL1d) model derived in Lemarié et al. (2021). The idea is to formulate the ABL3d model in terms of perturbation around a large-scale ambient state obtained from analyses or reanalyses. The proposed formulation is both simple and efficient from the computational point of view. Using numerical experiments and comparison with LES simulations, we illustrate the improvement brought by ABL3d compared with ABL1d. The ABL1d approach was designed to represent the modulation of the atmospheric boundary layer by sea-surface anomalies. The ABL3d model has the same ability, but it can also represent horizontal advection and the adjustment of the atmospheric horizontal pressure gradient to sea surface temperature anomalies.

## 1 Introduction

To improve the consistency between the scales resolved by oceanic models and the scales present in surface fluxes at the air-sea interface (often calculated using low-resolution atmospheric data), an aggregation (a.k.a. downscaling) of atmospheric data available at low resolution to the resolution of the oceanic model has been developed in Lemarié et al. (2021). In their work, these authors have designed a simplified atmospheric boundary layer model via a so-called physically based low-fidelity approach. The objective was to find an approach alternative to usual forcing strategies that makes the atmospheric boundary layer interactive while maintaining synoptic forcing (see Lemarié et al. (2021) for a detailed description of the motivations). For this approach, the focus was on the representation of the modulation of atmospheric turbulence by sea surface temperature and current anomalies. The atmospheric boundary layer model (referred to as ABL1d) thus defined was implemented in a general circulation model and was tested on a series of idealized and realistic benchmarks. The objective here is to investigate the possibility of relaxing the single-column assumption behind the ABL1d model as a first step toward a more realistic ABL3d model which would be able to represent horizontal advection and the adjustment of the atmospheric horizontal pressure gradient to the underlying sea surface temperature.

Our approach throughout this document is to consider the problem as a multiscale problem representing the interaction between a large-scale, a priori known, state (provided by analysis or reanalysis data) and the so-called micro-scales (i.e. the scales of the high-resolution oceanic model). In Sec. 2 different ways to represent such multiscale interactions are presented and illustrated via a simple tracer equation with source term. In Sec. 3 the full atmospheric primitive equations are considered and a perturbation method is applied as a technique to derive a simplified three-dimensional ABL model. The resulting model is then discretized and applied to an idealized experiment for strong winds crossing a persistent mid-latitude sea surface temperature (SST) gradient from cold to warm SST (Sec. 4.2) and from warm to cold SST (Sec. 4.3).

Corresponding author: Florian Lemarié, florian.lemarie@inria.fr
## 2 Techniques to represent multiscale interactions: the tracer transport equation with source term

To establish the theoretical framework of the study and for pedagogical reasons, we start by considering a generic transport equation with source term

$$ \partial_{\tau} \phi + u \partial_{\xi} \phi + w \partial_{z} \phi = S_{\phi} \quad (1a) $$

$$ \partial_{\xi} u + \partial_{z} w = 0 \quad (1b) $$

with $\phi$ a given quantity transported by the model, $\tau$ the time variable, $\xi$ the horizontal coordinate and $z$ the vertical coordinate.

### 2.1 Perturbation method around an ambient state

A first approach consists in defining a perturbation form of the original PDE around an ambient state. In this case, the ambient state and the perturbations are assumed to represent the same spatial and temporal scales. Such approach is quite common in the atmospheric community to improve the computational efficiency by filtering fast waves (e.g. Smolarkiewicz et al., 2019). Considering the system (3), an arbitrary ambient state $(\phi_a, u_a, w_a)$ that satisfies

$$ \partial_{\tau} \phi_a + u_a \partial_{\xi} \phi_a + w_a \partial_{z} \phi_a = S_{\phi_a} \quad (2a) $$

$$ \partial_{\xi} u_a + \partial_{z} w_a = 0 \quad (2b) $$

is assumed. Solutions of the original PDEs (3) are looked at in the form

$$ \psi' := \psi - \psi_a, \quad \psi = (\phi, u, v) $$

The system (3) in perturbation form reads

$$ \partial_{\tau} \phi' + u \partial_{\xi} \phi' + w \partial_{z} \phi' = -u' \partial_{\xi} \phi_a - w' \partial_{z} \phi_a + (S_{\phi} - S_{\phi_a}) \quad (3a) $$

$$ \partial_{\xi} u' + \partial_{z} w' = 0 \quad (3b) $$

The basic idea here would be to consider that the ambient state is provided by the, a priori known, large-scale atmospheric data. However, with this approach it is implicitly assumed that the ambient state and the perturbations live at the same scales in space and time which is not the case in the downscaling type of problem we are tackling here. In the following, we present alternative approaches assuming a difference in scale between the prescribed large-scale flow and the small-scale perturbations.

### 2.2 Multiple scales asymptotics

A more appropriate approach would be to adapt the derivation of the multi-scale model proposed by Majda (2007). To do so, we consider that the system of equations (3) is relevant to the typical scales of interest (the oceanic mesoscales) roughly characterized by

$$ L_m = 10 \text{ km} \quad \text{and} \quad T_m = 10 \text{ min} \quad (4) $$

Those scales are referred to as micro-scales. The aim is to exhibit multi-space multi-time scale solutions representing the interaction between the micro scales (4) and the atmospheric mesoscales characterized by

$$ L_M = \mathcal{O}(100 \text{ km}) \sim \varepsilon^{-1} L_m \quad \text{and} \quad T_M = \mathcal{O}(2 \text{ h}) \sim \varepsilon^{-1} T_m \quad (5) $$

know as meso-$\beta$ scales. We consider a regime with meso-$\beta$ scales interacting with microscales by looking for general asymptotic solutions of (3) in the form

$$ \phi = \bar{\phi}(X, z, T, t) + \phi'(X, x, z, T, t) + \mathcal{O}(\varepsilon^2) \quad (6a) $$

$$ u = \bar{u}(X, z, T, t) + u'(X, x, z, T, t) + \mathcal{O}(\varepsilon^2) \quad (6b) $$

$$ w = \varepsilon \bar{w}(X, z, T, t) + w'(X, x, z, T, t) + \mathcal{O}(\varepsilon^2) \quad (6c) $$
## 2 Techniques to represent multiscale interactions: the tracer transport equation with source term

To establish the theoretical framework of the study and for pedagogical reasons, we start by considering a generic transport equation with source term

$$ \partial_{\tau} \phi + u \partial_{\xi} \phi + w \partial_{z} \phi = S_{\phi} \quad (1a) $$

$$ \partial_{\xi} u + \partial_{z} w = 0 \quad (1b) $$

with $\phi$ a given quantity transported by the model, $\tau$ the time variable, $\xi$ the horizontal coordinate and $z$ the vertical coordinate.

### 2.1 Perturbation method around an ambient state

A first approach consists in defining a perturbation form of the original PDE around an ambient state. In this case, the ambient state and the perturbations are assumed to represent the same spatial and temporal scales. Such approach is quite common in the atmospheric community to improve the computational efficiency by filtering fast waves (e.g. Smolarkiewicz et al., 2019). Considering the system (3), an arbitrary ambient state $(\phi_a, u_a, w_a)$ that satisfies

$$ \partial_{\tau} \phi_a + u_a \partial_{\xi} \phi_a + w_a \partial_{z} \phi_a = S_{\phi_a} \quad (2a) $$

$$ \partial_{\xi} u_a + \partial_{z} w_a = 0 \quad (2b) $$

is assumed. Solutions of the original PDEs (3) are looked at in the form

$$ \psi' := \psi - \psi_a, \quad \psi = (\phi, u, v) $$

The system (3) in perturbation form reads

$$ \partial_{\tau} \phi' + u \partial_{\xi} \phi' + w \partial_{z} \phi' = -u' \partial_{\xi} \phi_a - w' \partial_{z} \phi_a + (S_{\phi} - S_{\phi_a}) \quad (3a) $$

$$ \partial_{\xi} u' + \partial_{z} w' = 0 \quad (3b) $$

The basic idea here would be to consider that the ambient state is provided by the, a priori known, large-scale atmospheric data. However, with this approach it is implicitly assumed that the ambient state and the perturbations live at the same scales in space and time which is not the case in the downscaling type of problem we are tackling here. In the following, we present alternative approaches assuming a difference in scale between the prescribed large-scale flow and the small-scale perturbations.

### 2.2 Multiple scales asymptotics

A more appropriate approach would be to adapt the derivation of the multi-scale model proposed by Majda (2007). To do so, we consider that the system of equations (3) is relevant to the typical scales of interest (the oceanic mesoscales) roughly characterized by

$$ L_m = 10 \text{ km} \quad \text{and} \quad T_m = 10 \text{ min} \quad (4) $$

Those scales are referred to as micro-scales. The aim is to exhibit multi-space multi-time scale solutions representing the interaction between the micro scales (4) and the atmospheric mesoscales characterized by

$$ L_M = \mathcal{O}(100 \text{ km}) \sim \varepsilon^{-1} L_m \quad \text{and} \quad T_M = \mathcal{O}(2 \text{ h}) \sim \varepsilon^{-1} T_m \quad (5) $$

know as meso-$\beta$ scales. We consider a regime with meso-$\beta$ scales interacting with microscales by looking for general asymptotic solutions of (3) in the form

$$ \phi = \bar{\phi}(X, z, T, t) + \phi'(X, x, z, T, t) + \mathcal{O}(\varepsilon^2) \quad (6a) $$

$$ u = \bar{u}(X, z, T, t) + u'(X, x, z, T, t) + \mathcal{O}(\varepsilon^2) \quad (6b) $$

$$ w = \varepsilon \bar{w}(X, z, T, t) + w'(X, x, z, T, t) + \mathcal{O}(\varepsilon^2) \quad (6c) $$
where $X = \varepsilon x$ and $T = \varepsilon t$ are respectively the large spatial and long temporal scales. The overline corresponds to a horizontal spatial average of $\phi$ over the microscales (thus removing the dependency on $x$):

$$ \bar{\phi}(X, z, T, t) = \lim_{L \to \infty} \int_{-L}^{L} \phi(X, x, z, T, t) dx $$

Furthermore, a temporal average over the microscales (thus removing the dependency on $t$) is also introduced

$$ \langle \phi \rangle (X, x, z, T) = \lim_{T \to \infty} \int_{-T}^{T} \phi(X, x, z, T, t) dt $$

In this multiple-scale framework, the space and time gradients are redefined as

$$ \partial_{\xi} \phi \to (\epsilon \partial_{X} + \partial_{x}) \phi, \quad \partial_{\tau} \phi \to (\epsilon \partial_{T} + \partial_{t}) \phi $$

Substituting the ansatz (6) into (3), we obtain:

$$ \partial_t \phi + u \partial_x \phi' + w' \partial_z \phi + \varepsilon [\partial_T \phi + u \partial_X \phi + \bar{w} \partial_z \phi] = S_{\phi} \quad (7a) $$

$$ \partial_x u + \partial_z w' + \varepsilon [\partial_x u + \partial_z \bar{w}] = 0 \quad (7b) $$

The derivation of the multiple scales model follows the 3 steps:

### 1. Mesoscale spatial mean on the microtime scales

Considering the order 1 terms in (7) we obtain

$$ \partial_t \phi + u \partial_x \phi' + w' \partial_z \phi = S_{\phi} \quad (8a) $$

$$ \partial_x u + \partial_z w' = 0 \quad (8b) $$

which can be written in a conservative way

$$ \partial_t \phi + \partial_x (u \phi') + \partial_z (w' \phi') + w' \partial_z \bar{\phi} = S_{\phi} \quad (9) $$

A mesoscale spatial averaging of this equation (assuming that $a'\bar{b} = 0$) leads to

$$ \partial_t \bar{\phi} + \partial_z \overline{w'\phi'} = \bar{S}_{\phi} \quad (10) $$

It is interesting to note here that we recover a single-column model similar to the ABL1d model proposed in Lemarié et al. (2021). This means that the ABL1d model is relevant to represent the fast temporal evolution at the micro time scales but not to represent micro scales in space as its formulation is obtained through mesoscale spatial averaging.

### 2. Space and time mesoscale-averaged equations

Now considering the order $\varepsilon$ terms in (7) we obtain

$$ \partial_T \phi + u \partial_X \phi + \bar{w} \partial_z \phi = 0 \quad (11a) $$

$$ \partial_x u + \partial_z \bar{\phi} = 0 \quad (11b) $$

which can be written in a conservative way as

$$ \partial_T \phi + \partial_x (u' \phi') + \partial_x (\bar{u} \phi) + \partial_x (\bar{u}' \phi') + \partial_x (u' \bar{\phi}) + \partial_z (\bar{w} \phi) + \partial_z (\bar{w}' \phi) = 0 $$

After a spatiotemporal averaging, we obtain

$$ \partial_T \langle \bar{\phi} \rangle + \partial_x \langle u' \phi' \rangle + \partial_x \langle \bar{u} \phi \rangle + \partial_z \langle \bar{w} \phi \rangle = \langle \bar{S}_{\phi} \rangle $$

which corresponds to the meso-$\beta$ scale averaged equations.
where $X = \varepsilon x$ and $T = \varepsilon t$ are respectively the large spatial and long temporal scales.
The overline corresponds to a horizontal spatial average of $\phi$ over the microscales (thus
removing the dependency on $x$):

$$
\bar{\phi}(X, z, T, t) = \lim_{L \to \infty} \int_{-L}^{L} \phi(X, x, z, T, t) \, dx
$$

Furthermore, a temporal average over the microscales (thus removing the dependency
on $t$) is also introduced

$$
\langle \phi \rangle (X, x, z, T) = \lim_{T \to \infty} \int_{-T}^{T} \phi(X, x, z, T, t) \, dt
$$

In this multiple-scale framework, the space and time gradients are redefined as

$$
\partial_{\xi} \phi \rightarrow (\epsilon \partial_{X} + \partial_{x}) \phi, \quad \partial_{\tau} \phi \rightarrow (\epsilon \partial_{T} + \partial_{t}) \phi
$$

Substituting the ansatz (6) into (3), we obtain:

$$
\partial_t \phi + u \partial_x \phi' + w' \partial_z \phi + \epsilon [\partial_T \phi + u \partial_X \phi + \bar{w} \partial_z \phi] = S_{\phi} \quad (7a)
$$

$$
\partial_x u + \partial_z w' + \epsilon [\partial_x u + \partial_z \bar{w}] = 0 \quad (7b)
$$

The derivation of the multiple scales model follows the 3 steps:

1. Mesoscale spatial mean on the microtime scales

Considering the order 1 terms in (7) we obtain

$$
\partial_t \phi + u \partial_x \phi' + w' \partial_z \phi = S_{\phi} \quad (8a)
$$

$$
\partial_x u + \partial_z w' = 0 \qquad (8b)
$$

which can be written in a conservative way

$$
\partial_t \phi + \partial_x (u\phi') + \partial_z (w'\phi') + w' \partial_z \bar{\phi} = S_{\phi} \quad (9)
$$

A mesoscale spatial averaging of this equation (assuming that $a'b = 0$) leads to

$$
\partial_t \bar{\phi} + \partial_z \overline{w'\phi'} = \bar{S}_\phi \quad (10)
$$

It is interesting to note here that we recover a single-column model similar to the
ABL1d model proposed in Lemarié et al. (2021). This means that the ABL1d model
is relevant to represent the fast temporal evolution at the micro time scales but
not to represent micro scales in space as its formulation is obtained through mesoscale
spatial averaging.

2. Space and time mesoscale-averaged equations

Now considering the order $\varepsilon$ terms in (7) we obtain

$$
\partial_T \phi + u \partial_X \phi + \bar{w} \partial_z \phi = 0 \quad (11a)
$$

$$
\partial_x u + \partial_z \bar{\phi} = 0 \qquad (11b)
$$

which can be written in a conservative way as

$$
\partial_T \phi + \partial_X (u' \phi') + \partial_X (\bar{u} \bar{\phi}) + \partial_X (\bar{u} \phi') + \partial_X (u' \bar{\phi}) + \partial_z (\bar{w} \bar{\phi}) + \partial_z (\bar{w} \phi') = 0
$$

After a spatiotemporal averaging, we obtain

$$
\partial_T \langle \bar{\phi} \rangle + \partial_X \langle u' \phi' \rangle + \partial_X \langle \bar{u} \bar{\phi} \rangle + \partial_z \langle \bar{w} \bar{\phi} \rangle = \langle \bar{S}_\phi \rangle
$$

which corresponds to the meso-$\beta$ scale averaged equations.
### 3. Micro-scale equations

A slight reformulation of (9) leads to

$$ \partial_t \bar{\phi} + \partial_t \phi' + \partial_x((\bar{u} + u')\phi') + \partial_z(w'\phi') + w' \partial_z \bar{\phi} = S_{\phi} $$

Now expressing $\partial_t \bar{\phi}$ using (10) we finally get

$$ \partial_t \phi' + \partial_x((\bar{u} + u')\phi') + w' \partial_z \bar{\phi} + \partial_z(w'\phi') = (S_{\phi} - \bar{S}_{\phi}) + \partial_z(\bar{w}'\phi') \quad (12) $$

In our context the mesoscale averaged quantities $\bar{\phi}$, $\bar{u}$ and $\bar{w}$ are known but additional assumptions are necessary to express the right-hand-side in a closed form. First, we consider that the source term is not significantly affected by the mesoscale spatial average, meaning that $S_{\phi} \approx \bar{S}_{\phi}$ (this implicitly assumes that disturbances have no significant influence on micro-physics and radiation.). The closure for the $\bar{w}'\phi'$ term will be discussed later.

## 2.3 Perturbation method around a time-dependent reference state

We now propose to apply the perturbation method described in Brereton et al. (2019) to derive a simplified ABL model. We again start from the generic transport equation with a source term $S_{\phi}$ meant to represent sub-grid scales (SGS)

$$ \partial_t \phi + u \partial_x \phi + w \partial_z \phi = S_{\phi} \quad (13a) $$

$$ \partial_x u + \partial_z w = 0 \quad (13b) $$

with $t$ the time coordinate and $x$ the horizontal coordinate. We assume that a horizontally averaged mean flow $(\tilde{\phi}, \tilde{u}, \tilde{w})$ is known and are interested in the fluctuating part $(\delta\phi, \delta u, \delta w)$ around this state. A general solution is sought in the form

$$ \phi(x, y, z, t) = \tilde{\phi}(z, t) + \delta\phi(x, y, z, t) \quad (14a) $$

$$ u(x, y, z, t) = \tilde{u}(z, t) + \delta u(x, y, z, t) \quad (14b) $$

$$ w(x, y, z, t) = \tilde{w}(z, t) + \delta w(x, y, z, t) \quad (14c) $$

The partitioning between the ambient state and the fluctuations is achieved using the horizontal mean $(\bar{\phi} = \frac{1}{A} \int_A \phi \, dx \, dy)$ such that

$$ \bar{\phi} = \tilde{\phi}(z, t), \quad \bar{\delta\phi} = 0 $$

The equations for the horizontally averaged flow are

$$ \partial_t \tilde{\phi} + (\delta u) \partial_x (\delta \phi) + \tilde{w} \partial_z \tilde{\phi} + (\delta w) \partial_z (\delta \phi) = \bar{S}_{\phi} \quad (15a) $$

$$ \partial_z \tilde{w} = 0 \quad (15b) $$

Subtracting this equation from (13) we obtain

$$ \partial_t(\delta\phi) + \partial_x((\tilde{u} + \delta u)\delta\phi) + (\delta w)\partial_z\tilde{\phi} + \partial_z(\delta w\delta\phi) \\ = (S_{\phi} - \bar{S}_{\phi}) - \partial_z(\tilde{w}\delta\phi) + (\overline{\delta w})\partial_z(\delta\phi) + (\overline{\delta u})\partial_x(\delta\phi) \quad (16) $$

with the continuity equation

$$ \partial_x(\delta u) + \partial_z(\delta w) = 0 $$

The left-hand side of (16) is formally equivalent to that of (12). The additional terms on the right-hand side are necessary to ensure that $\partial_t(\delta\phi) = 0$ (this can be interpreted as nudging to prevent the solution from drifting too far from its large-scale reference state). As mentioned in Brereton et al. (2019), instead of computing $(\delta u)\partial_x(\delta\phi)$ explicitly it is possible to simply remove $\partial_t(\delta\phi)$ from (16).
### 3. Micro-scale equations

A slight reformulation of (9) leads to

$$ \partial_t \bar{\phi} + \partial_t \phi' + \partial_x ((\bar{u} + u')\phi') + \partial_z (w'\phi') + w' \partial_z \bar{\phi} = S_{\phi} $$

Now expressing $\partial_t \bar{\phi}$ using (10) we finally get

$$ \partial_t \phi' + \partial_x ((\bar{u} + u')\phi') + w' \partial_z \bar{\phi} + \partial_z (w'\phi') = (S_{\phi} - \bar{S}_{\phi}) + \partial_z (w'\phi') \quad (12) $$

In our context the mesoscale averaged quantities $\bar{\phi}$, $\bar{u}$ and $\bar{w}$ are known but additional assumptions are necessary to express the right-hand-side in a closed form. First, we consider that the source term is not significantly affected by the mesoscale spatial average, meaning that $S_{\phi} \approx \bar{S}_{\phi}$ (this implicitly assumes that disturbances have no significant influence on micro-physics and radiation.). The closure for the $w'\phi'$ term will be discussed later.

## 2.3 Perturbation method around a time-dependent reference state

We now propose to apply the perturbation method described in Brereton et al. (2019) to derive a simplified ABL model. We again start from the generic transport equation with a source term $S_{\phi}$ meant to represent sub-grid scales (SGS)

$$ \partial_t \phi + u \partial_x \phi + w \partial_z \phi = S_{\phi} \quad (13a) $$

$$ \partial_x u + \partial_z w = 0 \quad (13b) $$

with $t$ the time coordinate and $x$ the horizontal coordinate. We assume that a horizontally averaged mean flow $(\tilde{\phi}, \tilde{u}, \tilde{w})$ is known and are interested in the fluctuating part $(\delta\phi, \delta u, \delta w)$ around this state. A general solution is sought in the form

$$ \phi(x, y, z, t) = \tilde{\phi}(z, t) + \delta\phi(x, y, z, t) \quad (14a) $$

$$ u(x, y, z, t) = \tilde{u}(z, t) + \delta u(x, y, z, t) \quad (14b) $$

$$ w(x, y, z, t) = \tilde{w}(z, t) + \delta w(x, y, z, t) \quad (14c) $$

The partitioning between the ambient state and the fluctuations is achieved using the horizontal mean $(\bar{\phi} = \frac{1}{A} \int_A \phi \, dx \, dy)$ such that

$$ \bar{\phi} = \tilde{\phi}(z, t), \quad \bar{\delta\phi} = 0 $$

The equations for the horizontally averaged flow are

$$ \partial_t \tilde{\phi} + (\delta u) \partial_x (\delta \phi) + \tilde{w} \partial_z \tilde{\phi} + (\delta w) \partial_z (\delta \phi) = \bar{S}_{\phi} \quad (15a) $$

$$ \partial_z \tilde{w} = 0 \quad (15b) $$

Subtracting this equation from (13) we obtain

$$ \partial_t (\delta \phi) + \partial_x ((\tilde{u} + \delta u)\delta \phi) + (\delta w) \partial_z \tilde{\phi} + \partial_z (\delta w \delta \phi) = (S_{\phi} - \bar{S}_{\phi}) - \partial_z (\tilde{w} \delta \phi) + (\overline{\delta w}) \partial_z (\delta \phi) + (\overline{\delta u}) \partial_x (\delta \phi) \quad (16) $$

with the continuity equation

$$ \partial_x (\delta u) + \partial_z (\delta w) = 0 $$

The left-hand side of (16) is formally equivalent to that of (12). The additional terms on the right-hand side are necessary to ensure that $\partial_t (\delta \phi) = 0$ (this can be interpreted as nudging to prevent the solution from drifting too far from its large-scale reference state). As mentioned in Brereton et al. (2019), instead of computing $(\delta u)\partial_x (\delta \phi)$ explicitly it is possible to simply remove $\partial_t (\delta \phi)$ from (16).
We have seen so far two possible approaches to derive a simplified model for the microscales forced by the mesoscales. In the two approaches the microscale model is formulated in such a way that the spatial and/or temporal average of the microscale perturbations is zero. It means that we should not need any nudging (a.k.a. restoring) toward the mesoscale data as in the ABL1d.

We now extend the philosophy behind those approaches to the full system of equations representative of atmospheric circulation.

# 3 A three-dimensional atmospheric boundary layer model

## 3.1 Continuous equations (dry case)

The starting point for the derivation is to select an adequate system of equations capable of representing the multi-scale phenomenon of interest. While most large-scale models rely on the compressible quasi-hydrostatic system of equations (the "primitive equations"), the non-hydrostatic anelastic system is used for smaller-scale motions. By prioritizing computational efficiency over accuracy we will consider the primitive equations as formulated by (Arakawa & Konor, 2009) and (Konor, 2014). In the 2D x-z case without rotation and moisture effects:

$$ \partial_t u + u \partial_x u + w \partial_z u = -c_p \theta \partial_x \pi + P_u(\theta, u) \quad (17a) $$

$$ c_p \theta \partial_z \pi = -g \quad (17b) $$

$$ \partial_x(\rho u) + \partial_z(\rho w) = -\partial_t \rho \quad (17c) $$

$$ \partial_t \theta + u \partial_x \theta + w \partial_z \theta = \frac{Q(\theta, u)}{c_p \pi} \quad (17d) $$

$$ \rho = \frac{p_0}{R\theta} \pi^{\gamma}, \quad \gamma = \frac{c_p - R}{R} \quad (17e) $$

with $\pi = T/\theta$, $Q$ a term representing the heating rate per unit mass, $P_u$ the parameterized physics, and $p_0$ a constant reference pressure.

System (17) is written in a classical z-coordinate which will be the coordinate used in ABL3d for the sake of simplicity and for reasons related to the computational cost associated with interpolations of large-scale data in the case of a vertical grid changing over time. Note that in system (17) the vertical velocity is computed in a diagnostic way from the continuity equation (17c) where $\rho$ is determined through the equation of state using $\theta$ and $\pi$.

Our derivation follows from the successive ansatzes:

* Continuity equation

Ansatz: $\rho = \bar{\rho}(X, z, T, t) + \epsilon \rho'(X, x, z, T, t)$, $u = \bar{u} + u'(X, x, z, T, t)$, $w = \epsilon \bar{w} + w'(X, x, z, T, t)$
The order 1 equation is

$$ \partial_z(\bar{\rho}w') = -\partial_t \bar{\rho} - \bar{\rho} \partial_x u' $$

where $\bar{\rho}$ is provided by the large-scale data.

* Hydrostatic equation

Ansatz: $\theta = \bar{\theta}(X, z, T, t) + \theta'(X, x, z, T, t)$, $\pi = \bar{\pi} + \pi'(X, x, z, T, t)$

$$ \begin{aligned} c_p \bar{\theta} \partial_z \bar{\pi} + c_p \bar{\theta} \partial_z \pi' + c_p \theta' \partial_z \pi' + c_p \theta' \partial_z \bar{\pi} &= -g \\ c_p \bar{\theta} \partial_z \pi' &\approx -c_p \theta' \partial_z \bar{\pi} \end{aligned} $$
We have seen so far two possible approaches to derive a simplified model for the microscales forced by the mesoscales. In the two approaches the microscale model is formulated in such a way that the spatial and/or temporal average of the microscale perturbations is zero. It means that we should not need any nudging (a.k.a. restoring) toward the mesoscale data as in the ABL1d.

We now extend the philosophy behind those approaches to the full system of equations representative of atmospheric circulation.

# 3 A three-dimensional atmospheric boundary layer model

## 3.1 Continuous equations (dry case)

The starting point for the derivation is to select an adequate system of equations capable of representing the multi-scale phenomenon of interest. While most large-scale models rely on the compressible quasi-hydrostatic system of equations (the "primitive equations"), the non-hydrostatic anelastic system is used for smaller-scale motions. By prioritizing computational efficiency over accuracy we will consider the primitive equations as formulated by (Arakawa & Konor, 2009) and (Konor, 2014). In the 2D x-z case without rotation and moisture effects:

$$ \partial_t u + u \partial_x u + w \partial_z u = -c_p \theta \partial_x \pi + P_u(\theta, u) \quad (17a) $$

$$ c_p \theta \partial_z \pi = -g \quad (17b) $$

$$ \partial_x(\rho u) + \partial_z(\rho w) = -\partial_t \rho \quad (17c) $$

$$ \partial_t \theta + u \partial_x \theta + w \partial_z \theta = \frac{Q(\theta, u)}{c_p \pi} \quad (17d) $$

$$ \rho = \frac{p_0}{R\theta} \pi^{\gamma}, \quad \gamma = \frac{c_p - R}{R} \quad (17e) $$

with $\pi = T/\theta$, $Q$ a term representing the heating rate per unit mass, $P_u$ the parameterized physics, and $p_0$ a constant reference pressure.

System (17) is written in a classical z-coordinate which will be the coordinate used in ABL3d for the sake of simplicity and for reasons related to the computational cost associated with interpolations of large-scale data in the case of a vertical grid changing over time. Note that in system (17) the vertical velocity is computed in a diagnostic way from the continuity equation (17c) where $\rho$ is determined through the equation of state using $\theta$ and $\pi$.

Our derivation follows from the successive ansatzes:

* Continuity equation

Ansatz: $\rho = \bar{\rho}(X, z, T, t) + \epsilon \rho'(X, x, z, T, t)$, $u = \bar{u} + u'(X, x, z, T, t)$, $w = \epsilon \bar{w} + w'(X, x, z, T, t)$
The order 1 equation is

$$ \partial_z(\bar{\rho}w') = -\partial_t \bar{\rho} - \bar{\rho} \partial_x u' $$

where $\bar{\rho}$ is provided by the large-scale data.

* Hydrostatic equation

Ansatz: $\theta = \bar{\theta}(X, z, T, t) + \theta'(X, x, z, T, t)$, $\pi = \bar{\pi} + \pi'(X, x, z, T, t)$

$$ \begin{aligned} c_p \bar{\theta} \partial_z \bar{\pi} + c_p \bar{\theta} \partial_z \pi' + c_p \theta' \partial_z \pi' + c_p \theta' \partial_z \bar{\pi} &= -g \\ c_p \bar{\theta} \partial_z \pi' &\approx -c_p \theta' \partial_z \bar{\pi} \end{aligned} $$
and thus

$$c_p \bar{\theta} \partial_z \pi' = g \frac{\theta'}{\bar{\theta}}$$

## Horizontal velocity equation

The original equation is

$$\partial_t u + u \partial_x u + w \partial_z u = -c_p \theta \partial_x \pi + P_u(\theta, u)$$

The order 1 equation is

$$\partial_t u' + (\bar{u} + u') \partial_x u' + w' \partial_z (\bar{u} + u') = -c_p (\bar{\theta} + \theta') \partial_x \pi' + (P_u(\bar{\theta} + \theta', \bar{u} + u') - \partial_t \bar{u})$$

## Thermodynamic equation

Considering the following thermodynamics equation

$$\partial_t \theta + u \partial_x \theta + w \partial_z \theta = \frac{Q(\theta, u)}{c_p \pi}$$

we simply obtain

$$\partial_t \theta' + (\bar{u} + u') \partial_x \theta' + w' \partial_z (\bar{\theta} + \theta') = \left( \frac{Q(\bar{\theta} + \theta', \bar{u} + u')}{c_p \bar{\pi}} - \partial_t \bar{\theta} \right)$$

We obtain the closed system of equations in a conservative form

$$\partial_t (\bar{\rho} u') + \partial_x (\bar{\rho} (\bar{u} + u') u') + \partial_z (\bar{\rho} w') u' = -\bar{\rho} w' \partial_z \bar{u} - \bar{\rho} c_p (\bar{\theta} + \theta') \partial_x \pi' + \bar{\rho} (P_u - \partial_t \bar{u}) \quad (18a)$$

$$c_p \bar{\theta} \partial_z \pi' = g \frac{\theta'}{\bar{\theta}} \quad (18b)$$

$$\partial_z (\bar{\rho} w') = -\partial_t \bar{\rho} - \bar{\rho} \partial_x u' \quad (18c)$$

$$\partial_t (\bar{\rho} \theta') + \partial_x (\bar{\rho} (\bar{u} + u') \theta') + \partial_z (\bar{\rho} w' \theta') = -\bar{\rho} w' \partial_z \bar{\theta} + \bar{\rho} \left( \frac{Q}{c_p \bar{\pi}} - \partial_t \bar{\theta} \right) \quad (18d)$$

The boundary conditions are

$$w'(z = z_{\text{sfc}}) = w'(z = z_{\text{top}}) = 0, \quad \theta'(z = z_{\text{top}}) = 0, \quad \theta'(z = z_{\text{sfc}}) = \bar{\theta}_s - \theta_s$$

with $\theta_s$ the sea surface temperature (SST) from the high-resolution oceanic model and $\bar{\theta}_s$ the SST used to compute the ambient state. For $\pi'$ we have $\pi'(z = z_{\text{top}}) = 0$ (since $\theta'(z = z_{\text{top}}) = 0$). To enforce this condition we can simply consider that

$$\pi'(z_{\text{sfc}}) = -\frac{g}{c_p} \int_{z_{\text{sfc}}}^{z_{\text{top}}} \frac{\theta'}{\bar{\theta}^2} \, dz$$

A simple procedure to enforce the condition $w'(z = z_{\text{sfc}}) = w'(z = z_{\text{top}}) = 0$ is to use the so-called O'Brien adjustment of vertical velocity (see below in the discretization aspects section).

## 3.2 Virtual temperature effect

Considering only water vapor mixing ratio $q_v$, the virtual potential temperature is defined by

$$\theta_v = \theta (1 + \beta q_v)$$

where $\beta = \frac{R_{\text{vap}}}{R_{\text{gas}}} - 1 \approx 0.6078$. $\theta$ must be replaced by $\theta_v$ in the hydrostatic and horizontal velocity equations and an additional equation for the water vapor mixing ratio is required

$$\partial_t (\bar{\rho} q_v') + \partial_x (\bar{\rho} (\bar{u} + u') q_v') + \partial_z (\bar{\rho} w' q_v') = -\bar{\rho} w' \partial_z \bar{q}_v + \bar{\rho} (P_q - \partial_t \bar{q}_v)$$
and thus

$$c_p \bar{\theta} \partial_z \pi' = g \frac{\theta'}{\bar{\theta}}$$

### Horizontal velocity equation

The original equation is

$$\partial_t u + u \partial_x u + w \partial_z u = -c_p \theta \partial_x \pi + P_u(\theta, u)$$

The order 1 equation is

$$\partial_t u' + (\bar{u} + u') \partial_x u' + w' \partial_z (\bar{u} + u') = -c_p (\bar{\theta} + \theta') \partial_x \pi' + (P_u(\bar{\theta} + \theta', \bar{u} + u') - \partial_t \bar{u})$$

### Thermodynamic equation

Considering the following thermodynamics equation

$$\partial_t \theta + u \partial_x \theta + w \partial_z \theta = \frac{Q(\theta, u)}{c_p \pi}$$

we simply obtain

$$\partial_t \theta' + (\bar{u} + u') \partial_x \theta' + w' \partial_z (\bar{\theta} + \theta') = \left( \frac{Q(\bar{\theta} + \theta', \bar{u} + u')}{c_p \bar{\pi}} - \partial_t \bar{\theta} \right)$$

We obtain the closed system of equations in a conservative form

$$\partial_t (\bar{\rho} u') + \partial_x (\bar{\rho}(\bar{u} + u') u') + \partial_z (\bar{\rho} w') u' = -\bar{\rho} w' \partial_z \bar{u} - \bar{\rho} c_p (\bar{\theta} + \theta') \partial_x \pi' + \bar{\rho} (P_u - \partial_t \bar{u}) \quad (18a)$$

$$c_p \bar{\theta} \partial_z \pi' = g \frac{\theta'}{\bar{\theta}} \quad (18b)$$

$$\partial_z (\bar{\rho} w') = -\partial_t \bar{\rho} - \bar{\rho} \partial_x u' \quad (18c)$$

$$\partial_t (\bar{\rho} \theta') + \partial_x (\bar{\rho}(\bar{u} + u') \theta') + \partial_z (\bar{\rho} w' \theta') = -\bar{\rho} w' \partial_z \bar{\theta} + \bar{\rho} \left( \frac{Q}{c_p \bar{\pi}} - \partial_t \bar{\theta} \right) \quad (18d)$$

The boundary conditions are

$$w'(z = z_{\text{sfc}}) = w'(z = z_{\text{top}}) = 0, \quad \theta'(z = z_{\text{top}}) = 0, \quad \theta'(z = z_{\text{sfc}}) = \bar{\theta}_s - \theta_s$$

with $\theta_s$ the sea surface temperature (SST) from the high-resolution oceanic model and $\bar{\theta}_s$ the SST used to compute the ambient state. For $\pi'$ we have $\pi'(z = z_{\text{top}}) = 0$ (since $\theta'(z = z_{\text{top}}) = 0$). To enforce this condition we can simply consider that

$$\pi'(z_{\text{sfc}}) = -\frac{g}{c_p} \int_{z_{\text{sfc}}}^{z_{\text{top}}} \frac{\theta'}{\bar{\theta}^2} \, dz$$

A simple procedure to enforce the condition $w'(z = z_{\text{sfc}}) = w'(z = z_{\text{top}}) = 0$ is to use the so-called O'Brien adjustment of vertical velocity (see below in the discretization aspects section).

## 3.2 Virtual temperature effect

Considering only water vapor mixing ratio $q_v$, the virtual potential temperature is defined by

$$\theta_v = \theta (1 + \beta q_v)$$

where $\beta = \frac{R_{\text{vap}}}{R_{\text{gas}}} - 1 \approx 0.6078$. $\theta$ must be replaced by $\theta_v$ in the hydrostatic and horizontal velocity equations and an additional equation for the water vapor mixing ratio is required

$$\partial_t (\bar{\rho} q_v') + \partial_x (\bar{\rho}(\bar{u} + u') q_v') + \partial_z (\bar{\rho} w' q_v') = -\bar{\rho} w' \partial_z \bar{q}_v + \bar{\rho} (P_q - \partial_t \bar{q}_v)$$
### 3.3 Discretization aspects

We assume a standard grid arrangement (density and velocities are located in the center of the cells on the vertical) and we consider N grid cells in the vertical with thickness $h_k = z_{k+1/2} - z_{k-1/2}$ ($z_{1/2} = 0$ and $z_{N+1/2} = z_{top}$). The turbulent quantities like turbulent kinetic energy k and eddy diffusivities $K_X$ are naturally located on the interfaces at $z_{k+1/2}$. To enforce the condition $(w')_{1/2} = (w')_{N+1/2} = 0$ we use the following procedure:

1. Compute $w_1$ by integration of the continuity equation

$$ (w_1)_{1/2} = 0 \\ (w_1)_{k+1/2} = (w_1)_{k-1/2} - h_k(\partial_x u)_k $$

2. redistribute linearly the excess $(w_1)_{N+1/2}$

$$ w_{k+1/2} = (w_1)_{k+1/2} - \left( \frac{z_{k+1/2}}{z_{N+1/2}} \right) (w_1)_{N+1/2} $$

For the first numerical tests, a very simple discretization strategy has been implemented:

1. Compute surface fluxes using variables at time n

2. Compute turbulent viscosity/diffusivity to compute $P_u$, $Q$, and $P_q$

3. Compute vertical velocity $(w')^n$ from continuity equation using $(u')^n$

4. Advance tracers to time $n+1$ via horizontal and vertical advection (first-order upwind scheme), $(w')^n \partial_z \bar{\theta}$, and vertical diffusion

5. Compute open boundary conditions for tracers

6. Compute Exner function via the hydrostatic balance using $(\theta_v')^{n+1}$

7. Advance horizontal velocities to time $n+1$ via horizontal and vertical advection (first-order upwind scheme), $(w')^n \partial_z \bar{u}$, horizontal pressure gradient, vertical velocity

8. Compute open boundary conditions for velocities

Advection is computed using a first-order upwind scheme while other terms are computed using simple second-order schemes.

## 4 Open ocean test: Winds across a midlatitude SST front

An idealized experiment particularly relevant for the coupling of the marine atmospheric boundary layer (MABL) with mesoscale oceanic eddies (and potentially submesoscale fronts) has been initially suggested by Spall (2007) and then revised by Kilpatrick et al. (2014). In Lemarié et al. (2021), limitations of the ABL1d model were illustrated using this experiment, in particular it was shown that in the frontal region, the effect of horizontal advection is predominant and the ABL1d model can not reproduce the horizontal lag seen in the reference solution when passing over the front. In the following, we will be interested again in this experiment to illustrate the improvements brought by the ABL3d model (18). For each configuration, a reference solution is obtained from the mesoscale non-hydrostatic model (MesoNH) v5.5.0 (Lac et al., 2018), where microphysics and radiation packages have not been activated.

### 4.1 Setup

The geometry of the problem is two-dimensional x-z with an SST front along the x-axis

$$ \theta_s(x) = 288.95 + \frac{\Delta\theta}{2} \tanh\left(\frac{x}{L_\theta}\right) $$
### 3.3 Discretization aspects

We assume a standard grid arrangement (density and velocities are located in the center of the cells on the vertical) and we consider N grid cells in the vertical with thickness $h_k = z_{k+1/2} - z_{k-1/2}$ ($z_{1/2} = 0$ and $z_{N+1/2} = z_{top}$). The turbulent quantities like turbulent kinetic energy k and eddy diffusivities $K_X$ are naturally located on the interfaces at $z_{k+1/2}$. To enforce the condition $(w')_{1/2} = (w')_{N+1/2} = 0$ we use the following procedure:

1. Compute $w_1$ by integration of the continuity equation

$$ (w_1)_{1/2} = 0 \\ (w_1)_{k+1/2} = (w_1)_{k-1/2} - h_k(\partial_x u)_k $$

2. redistribute linearly the excess $(w_1)_{N+1/2}$

$$ w_{k+1/2} = (w_1)_{k+1/2} - \left( \frac{z_{k+1/2}}{z_{N+1/2}} \right) (w_1)_{N+1/2} $$

For the first numerical tests, a very simple discretization strategy has been implemented:

1. Compute surface fluxes using variables at time n

2. Compute turbulent viscosity/diffusivity to compute $P_u$, $Q$, and $P_q$

3. Compute vertical velocity $(w')^n$ from continuity equation using $(u')^n$

4. Advance tracers to time $n+1$ via horizontal and vertical advection (first-order upwind scheme), $(w')^n \partial_z \bar{\theta}$, and vertical diffusion

5. Compute open boundary conditions for tracers

6. Compute Exner function via the hydrostatic balance using $(\theta_v')^{n+1}$

7. Advance horizontal velocities to time $n+1$ via horizontal and vertical advection (first-order upwind scheme), $(w')^n \partial_z \bar{u}$, horizontal pressure gradient, vertical velocity

8. Compute open boundary conditions for velocities

Advection is computed using a first-order upwind scheme while other terms are computed using simple second-order schemes.

## 4 Open ocean test: Winds across a midlatitude SST front

An idealized experiment particularly relevant for the coupling of the marine atmospheric boundary layer (MABL) with mesoscale oceanic eddies (and potentially submesoscale fronts) has been initially suggested by Spall (2007) and then revised by Kilpatrick et al. (2014). In Lemarié et al. (2021), limitations of the ABL1d model were illustrated using this experiment, in particular it was shown that in the frontal region, the effect of horizontal advection is predominant and the ABL1d model can not reproduce the horizontal lag seen in the reference solution when passing over the front. In the following, we will be interested again in this experiment to illustrate the improvements brought by the ABL3d model (18). For each configuration, a reference solution is obtained from the mesoscale non-hydrostatic model (MesoNH) v5.5.0 (Lac et al., 2018), where microphysics and radiation packages have not been activated.

### 4.1 Setup

The geometry of the problem is two-dimensional x-z with an SST front along the x-axis

$$ \theta_s(x) = 288.95 + \frac{\Delta\theta}{2} \tanh\left(\frac{x}{L_\theta}\right) $$
where $\Delta\theta = 3$ K, $L_{\theta} = 100$ km, and $x \in [-1800$ km, $1800$ km]. A zonal geostrophic wind of $15$ m s$^{-1}$ is prescribed balanced by a vertically homogeneous meridional pressure gradient. The wind thus flows over cold water before reaching a warm SST anomaly $3$ K warmer. We consider a dry case, the model is initialized $\forall x$ with

$$ \theta(z, t=0) = 288.95 + (N^2\theta^{\text{ref}}/g)z \\ q(z, t=0) = 0 $$

where $N^2 = 10^{-4}$ s$^{-2}$ and $\theta^{\text{ref}} = 288$ K. The velocities are systematically initialized with geostrophic winds. All simulations are run for 36 hours. The horizontal resolution is $\Delta x = 6$ km for ABL1d and ABL3d simulations and $\Delta x = 1$ km for MesoNH simulations. MesoNH equations are discretized with 91 vertical levels from the surface to 20 km height. The vertical resolution near the surface is $\Delta z = 5$ m and around $\Delta z = 50$ m at 2000 m height. The turbulence scheme is a 1.5-order closure similar to the one used in the ABL1d and ABL3d models and all simulations use the COARE bulk formulation to compute turbulent components of air-sea fluxes. For the ABL models, the top of the computational domain is $z_{\text{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_{\text{top}}$ with a first grid point located at $z = 10$ m.

For the numerical experiments using the ABL3d simplified model we consider an ambient state characterized by

$$ \bar{\theta} = (288.95 + \frac{N^2}{g/\theta^{\text{ref}}}z) \text{ K}, \bar{u} = 15 \text{ m s}^{-1}, \theta_s = 288.95 \text{ K}, \bar{\rho} = 1.22 \text{ kg m}^{-3}, p_{\text{sfc}} = 1013 \text{ hPa} \qquad (19) $$

## 4.2 Cold-to-warm experiment

For this first experiment, we consider that

$$ \theta'(z_{\text{sfc}}) = \frac{\Delta\theta}{2} \tanh\left(\frac{x}{L_{\theta}}\right) $$

with $\Delta\theta = 3$ K, $L_{\theta} = 100$ km and we assess how this perturbation of the surface temperature affects the ambient state. In Fig. 1 the modeled winds are represented using a large-eddy simulation (LES) with MesoNH model for reference. The improvements brought by the ABL3d compared to ABL1d are obvious, in particular, it reproduces the horizontal lag when passing the front whereas the ABL1d solution shows a local adjustment to the underlying sea surface temperature. Those improvements are also visible in the 10 m quantities (Fig. 2) albeit too strong meridional winds are obtained with ABL1d and ABL3d compared to MesoNH. Since the target application for the proposed model is the improvement of the forcing of ocean models, the important quantities to look at are the 10 m variables, as these are the ones used to calculate turbulent components of air-sea fluxes.

The ABL3d model can reproduce a more realistic momentum budget as seen in Fig. 3. The budget terms with ABL3d compare well with those from (Kilpatrick et al., 2014) (their figure 7) and show that in the frontal region, the advection and pressure gradient contributions are significant to shape the winds.

## 4.3 Warm-to-cold experiment

For this second experiment we consider that

$$ \theta'(z_{\text{sfc}}) = \frac{\Delta\theta}{2} \tanh\left(-\frac{x}{L_{\theta}}\right) $$

with $\Delta\theta = 3$ K, $L_{\theta} = 100$ km. Note that, in (Kilpatrick et al., 2014) this test is done by reversing the direction of the wind (i.e. $\bar{u} = -15$ m s$^{-1}$) while keeping the same
where $\Delta\theta = 3$ K, $L_{\theta} = 100$ km, and $x \in [-1800$ km, $1800$ km]. A zonal geostrophic wind of $15$ m s$^{-1}$ is prescribed balanced by a vertically homogeneous meridional pressure gradient. The wind thus flows over cold water before reaching a warm SST anomaly $3$ K warmer. We consider a dry case, the model is initialized $\forall x$ with

$$ \theta(z, t=0) = 288.95 + (N^2\theta^{\text{ref}}/g)z \\ q(z, t=0) = 0 $$

where $N^2 = 10^{-4}$ s$^{-2}$ and $\theta^{\text{ref}} = 288$ K. The velocities are systematically initialized with geostrophic winds. All simulations are run for 36 hours. The horizontal resolution is $\Delta x = 6$ km for ABL1d and ABL3d simulations and $\Delta x = 1$ km for MesoNH simulations. MesoNH equations are discretized with 91 vertical levels from the surface to 20 km height. The vertical resolution near the surface is $\Delta z = 5$ m and around $\Delta z = 50$ m at 2000 m height. The turbulence scheme is a 1.5-order closure similar to the one used in the ABL1d and ABL3d models and all simulations use the COARE bulk formulation to compute turbulent components of air-sea fluxes. For the ABL models, the top of the computational domain is $z_{\text{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_{\text{top}}$ with a first grid point located at $z = 10$ m.

For the numerical experiments using the ABL3d simplified model we consider an ambient state characterized by

$$ \bar{\theta} = (288.95 + \frac{N^2}{g/\theta^{\text{ref}}}z) \text{ K}, \bar{u} = 15 \text{ m s}^{-1}, \theta_s = 288.95 \text{ K}, \bar{\rho} = 1.22 \text{ kg m}^{-3}, p_{\text{sfc}} = 1013 \text{ hPa} \qquad (19) $$

## 4.2 Cold-to-warm experiment

For this first experiment, we consider that

$$ \theta'(z_{\text{sfc}}) = \frac{\Delta\theta}{2} \tanh\left(\frac{x}{L_{\theta}}\right) $$

with $\Delta\theta = 3$ K, $L_{\theta} = 100$ km and we assess how this perturbation of the surface temperature affects the ambient state. In Fig. 1 the modeled winds are represented using a large-eddy simulation (LES) with MesoNH model for reference. The improvements brought by the ABL3d compared to ABL1d are obvious, in particular, it reproduces the horizontal lag when passing the front whereas the ABL1d solution shows a local adjustment to the underlying sea surface temperature. Those improvements are also visible in the 10 m quantities (Fig. 2) albeit too strong meridional winds are obtained with ABL1d and ABL3d compared to MesoNH. Since the target application for the proposed model is the improvement of the forcing of ocean models, the important quantities to look at are the 10 m variables, as these are the ones used to calculate turbulent components of air-sea fluxes.

The ABL3d model can reproduce a more realistic momentum budget as seen in Fig. 3. The budget terms with ABL3d compare well with those from (Kilpatrick et al., 2014) (their figure 7) and show that in the frontal region, the advection and pressure gradient contributions are significant to shape the winds.

## 4.3 Warm-to-cold experiment

For this second experiment we consider that

$$ \theta'(z_{\text{sfc}}) = \frac{\Delta\theta}{2} \tanh\left(-\frac{x}{L_{\theta}}\right) $$

with $\Delta\theta = 3$ K, $L_{\theta} = 100$ km. Note that, in (Kilpatrick et al., 2014) this test is done by reversing the direction of the wind (i.e. $\bar{u} = -15$ m s$^{-1}$) while keeping the same
<figure><img src="image_3.png" /><figcaption>Figure 1: Zonal (top panels) and meridional (bottom panels) components of atmospheric winds for the reference MesoNH simulation (right panels), for ABL1d simulation (left panels) and ABL3d simulation (middle panels). Temperature contours are shown in white with a contour every 0.5 °C between 15 °C and 17.5 °C. The SST front is centered at $x = 0$ km and $\bar{u} = 15$ m s$^{-1}$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: Zonal (top) and meridional (middle) components of 10 m winds and 10 m temperature (bottom) for the reference MesoNH simulation (dashed red), for ABL1d simulation (solid black) and for ABL3d simulation (dashed black) for the winds across a mid-latitude SST front experiment. The SST front is centered at $x = 0$ km and $\bar{u} = 15$ m s$^{-1}$.</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>Figure 1: Zonal (top panels) and meridional (bottom panels) components of atmospheric winds for the reference MesoNH simulation (right panels), for ABL1d simulation (left panels) and ABL3d simulation (middle panels). Temperature contours are shown in white with a contour every 0.5 °C between 15 °C and 17.5 °C. The SST front is centered at $x = 0$ km and $\bar{u} = 15$ m s$^{-1}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Zonal (top) and meridional (middle) components of 10 m winds and 10 m temperature (bottom) for the reference MesoNH simulation (dashed red), for ABL1d simulation (solid black) and for ABL3d simulation (dashed black) for the winds across a mid-latitude SST front experiment. The SST front is centered at $x = 0$ km and $\bar{u} = 15$ m s$^{-1}$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3: Sections of zonal momentum budget terms (×10⁻⁴ m s⁻²) for the cold-to-warm case (with $\bar{u}$ = 15 m s⁻¹) with ABL3d model: advection, Coriolis, turbulent stress divergence, and pressure (look out for the color bars, which differ from panel to panel).</figcaption></figure>

SST as in the cold-to-warm case. For comparison with their results, we must invert the sign of the zonal and meridional velocities and do a symmetry around the x = 0 axis.

This case is more complex than the previous one as a shallow stable internal boundary layer forms over cool SST. The ABL3d model represents properly the “inertial lee wave” (inertial oscillation) which results in a descent-ascent-descent pattern (see Fig. 4 bottom right panel) located between x = 100 km and x = 1000 km, consistent with (Kilpatrick et al., 2014) (their figure 6). Overall the results shown in Fig. 4, 5 and 6 compares well with the WRF results presented in Kilpatrick et al. (2014) and with MesoNH LES simulations.

## 4.4 Role of the pressure gradient

In the cases shown so far the pressure gradient had little impact on the wind structure. Following Ayet and Redelsperger (2019), the impact of pressure relative to Coriolis on the momentum balance is measured by the non-dimensional number $P_c$

$$P_c = \frac{gh_e \Delta\theta}{\theta_0 f U_g L_\theta}$$

with $L_\theta$ the characteristic horizontal length scale of the SST front, $\Delta\theta$ the characteristic temperature difference across the SST front, and $h_e$ the MABL height. Considering $h_e$ = 600 m, we had $P_c$ = 0.4 in the previous experiments. To increase the relative role of pressure gradient over Coriolis, we can set $\bar{u}$ = 10 m s⁻¹ and $f$ = 5 × 10⁻⁵ s⁻¹ which results in $P_c$ = 1.2. The momentum budget in this configuration is shown in Fig. 7. The model seems to respond adequately to this change in the large-scale ambient state,
<figure><img src="image_4.png" /><figcaption>Figure 3: Sections of zonal momentum budget terms (×10⁻⁴ m s⁻²) for the cold-to-warm case (with $\bar{u}$ = 15 m s⁻¹) with ABL3d model: advection, Coriolis, turbulent stress divergence, and pressure (look out for the color bars, which differ from panel to panel).</figcaption></figure>

SST as in the cold-to-warm case. For comparison with their results, we must invert the sign of the zonal and meridional velocities and do a symmetry around the $x = 0$ axis.

This case is more complex than the previous one as a shallow stable internal boundary layer forms over cool SST. The ABL3d model represents properly the “inertial lee wave” (inertial oscillation) which results in a descent-ascent-descent pattern (see Fig. 4 bottom right panel) located between $x = 100$ km and $x = 1000$ km, consistent with (Kilpatrick et al., 2014) (their figure 6). Overall the results shown in Fig. 4, 5 and 6 compares well with the WRF results presented in Kilpatrick et al. (2014) and with MesoNH LES simulations.

## 4.4 Role of the pressure gradient

In the cases shown so far the pressure gradient had little impact on the wind structure. Following Ayet and Redelsperger (2019), the impact of pressure relative to Coriolis on the momentum balance is measured by the non-dimensional number $P_c$

$$P_c = \frac{gh_e \Delta\theta}{\theta_0 f U_g L_\theta}$$

with $L_\theta$ the characteristic horizontal length scale of the SST front, $\Delta\theta$ the characteristic temperature difference across the SST front, and $h_e$ the MABL height. Considering $h_e = 600$ m, we had $P_c = 0.4$ in the previous experiments. To increase the relative role of pressure gradient over Coriolis, we can set $\bar{u} = 10$ m s⁻¹ and $f = 5 \times 10^{-5}$ s⁻¹ which results in $P_c = 1.2$. The momentum budget in this configuration is shown in Fig. 7. The model seems to respond adequately to this change in the large-scale ambient state,
<figure><img src="image_6.png" /><figcaption>Figure 4: Zonal (top panels), meridional (middle panels) and vertical (bottom panels) components of atmospheric winds for the ABL1d simulation (left panels), ABL3d simulation (middle panels) and reference MesoNH simulation (right panels). Temperature contours are shown in white with a contour every 0.5 °C between 15 °C and 17.5 °C. The SST front is centered at $x = 0$ km (warm for negative $x$ and cold for positive $x$) and $\bar{u} = 15$ m s$^{-1}$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Zonal (top panels), meridional (middle panels) and vertical (bottom panels) components of atmospheric winds for the ABL1d simulation (left panels), ABL3d simulation (middle panels) and reference MesoNH simulation (right panels). Temperature contours are shown in white with a contour every 0.5 °C between 15 °C and 17.5 °C. The SST front is centered at $x = 0$ km (warm for negative $x$ and cold for positive $x$) and $\bar{u} = 15$ m s$^{-1}$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: Zonal (top) and meridional (middle) components of 10 m winds and 10 m temperature (bottom) for the ABL1d simulation (solid black), the ABL3d simulation (dashed black), and the reference MesoNH simulation (dashed red). The SST front is centered at $x = 0$ km and $\bar{u} = 15$ m s$^{-1}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Sections of zonal momentum budget terms ($\times 10^{-4} \text{ m s}^{-2}$) for the warm-to-cold case (with $\bar{u} = 15 \text{ m s}^{-1}$) with ABL3d model: advection, Coriolis, turbulent stress divergence, and pressure (look out for the color bars, which differ from panel to panel).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Zonal (top) and meridional (middle) components of 10 m winds and 10 m temperature (bottom) for the ABL1d simulation (solid black), the ABL3d simulation (dashed black), and the reference MesoNH simulation (dashed red). The SST front is centered at $x = 0$ km and $\bar{u} = 15$ m s$^{-1}$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Sections of zonal momentum budget terms ($\times 10^{-4} \text{ m s}^{-2}$) for the warm-to-cold case (with $\bar{u} = 15 \text{ m s}^{-1}$) with ABL3d model: advection, Coriolis, turbulent stress divergence, and pressure (look out for the color bars, which differ from panel to panel).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Sections of zonal momentum budget terms (×10⁻⁴m s⁻²) for the cold-to-warm case (with $\bar{u} = 10$ m s⁻¹) and $f = 5 \times 10^{-5}$ s⁻¹ with ABL3d model: advection, Coriolis, turbulent stress divergence, and pressure (look out for the color bars, which differ from panel to panel).</figcaption></figure>

and as expected the impact of the pressure gradient term on the wind structure is much more preponderant than in previous cases.

## 5 Conclusion

A three-dimensional extension of the single-column atmospheric boundary layer (ABL1d) model developed by Lemarié et al. (2021) has been derived. This development is made 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. A crucial hypothesis in the derivation of the ABL1d model was that the dominant process at the characteristic scale of the oceanic mesoscale is the modulation of atmospheric turbulence by sea surface temperature (SST) anomalies. However, it is now well documented that momentum vertical turbulent mixing, pressure gradient, Coriolis, and nonlinear advection are all important to the momentum balance in the marine atmospheric boundary layer in the vicinity of oceanic fronts. The shortcomings of the ABL1d approach in this respect are remedied by the proposed ABL3d model. This is unambiguously confirmed by numerical experiments in an idealized setting. Before we can deploy the ABL3d model in a realistic configuration, as we did with the ABL1d model, we need to think about the treatment of coastlines (and open boundaries in general).

## Acknowledgments

The author is extremely grateful to Joris Pianezze for providing the Meso-NH LES simulations. The author appreciates the funding from the European Union Horizon 2020
<figure><img src="image_8.png" /><figcaption>Figure 7: Sections of zonal momentum budget terms (×10⁻⁴m s⁻²) for the cold-to-warm case (with $\bar{u}$ = 10 m s⁻¹) and f = 5 × 10⁻⁵ s⁻¹ with ABL3d model: advection, Coriolis, turbulent stress divergence, and pressure (look out for the color bars, which differ from panel to panel).</figcaption></figure>

and as expected the impact of the pressure gradient term on the wind structure is much
more preponderant than in previous cases.

**5 Conclusion**

A three-dimensional extension of the single-column atmospheric boundary layer (ABL1d) model developed by Lemarié et al. (2021) has been derived. This development is made 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. A crucial hypothesis in the derivation of the ABL1d model was that the dominant process at the characteristic scale of the oceanic mesoscale is the modulation of atmospheric turbulence by sea surface temperature (SST) anomalies. However, it is now well documented that momentum vertical turbulent mixing, pressure gradient, Coriolis, and nonlinear advection are all important to the momentum balance in the marine atmospheric boundary layer in the vicinity of oceanic fronts. The shortcomings of the ABL1d approach in this respect are remedied by the proposed ABL3d model. This is unambiguously confirmed by numerical experiments in an idealized setting. Before we can deploy the ABL3d model in a realistic configuration, as we did with the ABL1d model, we need to think about the treatment of coastlines (and open boundaries in general).

**Acknowledgments**

The author is extremely grateful to Joris Pianezze for providing the Meso-NH LES sim-
ulations. The author appreciates the funding from the European Union Horizon 2020
research and innovation programme under grant agreement No 821926 (IMMERSE). This
study was carried out as part of the techonological defense project PROTEVS2 under
the auspices of French Ministry of the armies / DGA.

References

Arakawa, A., & Konor, C. S. (2009). Unification of the anelastic and quasi-hydrostatic systems of equations. *Mon. Weather Rev.*, **137**(2), 710 - 726.
doi: 10.1175/2008MWR2520.1

Ayet, A., & Redelsperger, J.-L. (2019). An analytical study of the atmospheric boundary-layer flow and divergence over an sst front. *Quart. J. Roy. Meteorol. Soc.*, **145**(723), 2549-2567. doi: https://doi.org/10.1002/qj.3578

Brereton, A., Tejada-Martinez, A. E., Palmer, M. R., & Polton, J. A. (2019). The perturbation method - a novel large-eddy simulation technique to model realistic turbulence: Application to tidal flow. *Ocean Modell.*, **135**, 31-39. doi: https://doi.org/10.1016/j.ocemod.2019.01.007

Kilpatrick, T., Schneider, N., & Qiu, B. (2014). Boundary layer convergence induced by strong winds across a midlatitude sst front. *J. Climate*, **27**(4), 1698-1718.

Konor, C. S. (2014). Design of a dynamical core based on the nonhydrostatic “unified system” of equations. *Mon. Weather Rev.*, **142**(1), 364 - 385. doi: 10.1175/MWR-D-13-00187.1

Lac, C., Chaboureau, J.-P., Masson, V., Pinty, J.-P., Tulet, P., Escobar, J., ...
Wautelet, P. (2018). Overview of the meso-nh model version 5.4 and
its applications. *Geosci. Model Dev.*, **11**(5), 1929–1969. doi: 10.5194/gmd-11-1929-2018

Lemarié, F., Samson, G., Redelsperger, J.-L., Giordani, H., Brivoal, T., & Madec, G. (2021). 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). *Geosci. Model Dev.*, **14**(1), 543-572. doi: 10.5194/gmd-14-543-2021

Majda, A. J. (2007). Multiscale models with moisture and systematic strategies for superparameterization. *J. Atmos. Sci.*, **64**(7), 2726 - 2734. doi: 10.1175/JAS3976.1

Smolarkiewicz, P. K., Kühnlein, C., & Wedi, N. P. (2019). Semi-implicit integrations of perturbation equations for all-scale atmospheric dynamics. *J. Comp. Phys.*, **376**, 145-159. doi: https://doi.org/10.1016/j.jcp.2018.09.032

Spall, M. (2007, 08). Midlatitude wind stress–sea surface temperature coupling in the vicinity of oceanic fronts. *J. Climate*, 20. doi: 10.1175/JCLI4234.1
research and innovation programme under grant agreement No 821926 (IMMERSE). This
study was carried out as part of the techonological defense project PROTEVS2 under
the auspices of French Ministry of the armies / DGA.

References

Arakawa, A., & Konor, C. S. (2009). Unification of the anelastic and quasi-hydrostatic systems of equations. *Mon. Weather Rev.*, **137**(2), 710 - 726. doi: 10.1175/2008MWR2520.1

Ayet, A., & Redelsperger, J.-L. (2019). An analytical study of the atmospheric boundary-layer flow and divergence over an sst front. *Quart. J. Roy. Meteorol. Soc.*, **145**(723), 2549-2567. doi: https://doi.org/10.1002/qj.3578

Brereton, A., Tejada-Martinez, A. E., Palmer, M. R., & Polton, J. A. (2019). The perturbation method - a novel large-eddy simulation technique to model realistic turbulence: Application to tidal flow. *Ocean Modell.*, **135**, 31-39. doi: https://doi.org/10.1016/j.ocemod.2019.01.007

Kilpatrick, T., Schneider, N., & Qiu, B. (2014). Boundary layer convergence induced by strong winds across a midlatitude sst front. *J. Climate*, **27**(4), 1698-1718.

Konor, C. S. (2014). Design of a dynamical core based on the nonhydrostatic “unified system” of equations. *Mon. Weather Rev.*, **142**(1), 364 - 385. doi: 10.1175/MWR-D-13-00187.1

Lac, C., Chaboureau, J.-P., Masson, V., Pinty, J.-P., Tulet, P., Escobar, J., ...
Wautelet, P. (2018). Overview of the meso-nh model version 5.4 and
its applications. *Geosci. Model Dev.*, **11**(5), 1929–1969. doi: 10.5194/gmd-11-1929-2018

Lemarié, F., Samson, G., Redelsperger, J.-L., Giordani, H., Brivoal, T., & Madec, G. (2021). 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). *Geosci. Model Dev.*, **14**(1), 543-572. doi: 10.5194/gmd-14-543-2021

Majda, A. J. (2007). Multiscale models with moisture and systematic strategies for superparameterization. *J. Atmos. Sci.*, **64**(7), 2726 - 2734. doi: 10.1175/JAS3976.1

Smolarkiewicz, P. K., Kühnlein, C., & Wedi, N. P. (2019). Semi-implicit integrations of perturbation equations for all-scale atmospheric dynamics. *J. Comp. Phys.*, **376**, 145-159. doi: https://doi.org/10.1016/j.jcp.2018.09.032

Spall, M. (2007, 08). Midlatitude wind stress–sea surface temperature coupling in the vicinity of oceanic fronts. *J. Climate*, **20**. doi: 10.1175/JCLI4234.1