# Spatial Variability of Nocturnal Stability Regimes in an Operational Weather Prediction Model

Marvin Kähnert, Harald Sodemann, Teresa M Remes, Carl Fortelius, Eric

Bazile, Igor Esau

► To cite this version:

Marvin Kähnert, Harald Sodemann, Teresa M Remes, Carl Fortelius, Eric Bazile, et al. Spatial Variability of Nocturnal Stability Regimes in an Operational Weather Prediction Model. Boundary-Layer Meteorology, 2022, 186 (2), pp.373 - 397. 10.1007/s10546-022-00762-1. meteo-04444315

HAL Id: meteo-04444315

https://meteofrance.hal.science/meteo-04444315

Submitted on 16 May 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
# Spatial Variability of Nocturnal Stability Regimes in an Operational Weather Prediction Model

Marvin Kähnert, Harald Sodemann, Teresa M Remes, Carl Fortelius, Eric

Bazile, Igor Esau

► To cite this version:

Marvin Kähnert, Harald Sodemann, Teresa M Remes, Carl Fortelius, Eric Bazile, et al. Spatial Variability of Nocturnal Stability Regimes in an Operational Weather Prediction Model. Boundary-Layer Meteorology, 2022, 186 (2), pp.373 - 397. 10.1007/s10546-022-00762-1. meteo-04444315

HAL Id: meteo-04444315

https://meteofrance.hal.science/meteo-04444315v1

Submitted on 16 May 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
# Spatial Variability of Nocturnal Stability Regimes in an Operational Weather Prediction Model

Marvin Kähnert<fnref n="1" /> Ⓒ · Harald Sodemann<fnref n="1" /> · Teresa M. Remes<fnref n="2" /> · Carl Fortelius<fnref n="3" /> · Eric Bazile<fnref n="4" /> · Igor Esau<fnref n="5" />

Received: 26 March 2022 / Accepted: 26 October 2022 / Published online: 15 November 2022
© The Author(s) 2022

## Abstract

Forecast errors in near-surface temperatures are a persistent issue for numerical weather prediction models. A prominent example is warm biases during cloud-free, snow-covered nights. Many studies attribute these biases to parametrized processes such as turbulence or radiation. Here, we focus on the contribution of physical processes to the nocturnal temperature development. We compare model timestep output of individual tendencies from parametrized processes in the weather prediction model AROME-Arctic to measurements from Sodankylä, Finland. Thereby, we differentiate between the weakly stable boundary layer (wSBL) and the very stable boundary layer (vSBL) regimes. The wSBL is characterized by continuous turbulent exchange within the near-surface atmosphere, causing near-neutral temperature profiles. The vSBL is characterized by a decoupling of the lowermost model level, low turbulent exchange, and very stable temperature profiles. In our case study, both regimes occur simultaneously on small spatial scales of about 5 km. In addition, we demonstrate the model's sensitivity towards an updated surface treatment, allowing for faster surface cooling. The updated surface parametrization has profound impacts on parametrized processes in both regimes. However, only modelled temperatures in the vSBL are impacted substantially, whereas more efficient surface cooling in the wSBL is compensated by an increased turbulent heat transport within the boundary layer. This study demonstrates the utility of individual tendencies for understanding process-related differences between model configurations and emphasizes the need for model studies to distinguish between the wSBL and vSBL for reliable model verification.

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

<footnote>Marvin Kähnert
marvin.kahnert@uib.no</footnote>

<footnote n="1">¹ Geophysical Institute and Bjerknes Centre for Climate Research, University of Bergen, Postbox 7803, 5020 Bergen, Norway</footnote>

<footnote n="2">² Development Centre for Weather Forecasting, Norwegian Meteorological Institute, Oslo, Norway</footnote>

<footnote n="3">³ Finnish Meteorological Institute, Helsinki, Finland</footnote>

<footnote n="4">⁴ CNRM - UMR3589, Météo-France and CNRS, Toulouse, France</footnote>

<footnote n="5">⁵ Nansen Environmental and Remote Sensing Center,UiT - The Arctic University of Norway, Tromso, Norway</footnote>
# Spatial Variability of Nocturnal Stability Regimes in an Operational Weather Prediction Model

Marvin Kähnert<fnref n="1" /> Ⓒ · Harald Sodemann<fnref n="1" /> · Teresa M. Remes<fnref n="2" /> · Carl Fortelius<fnref n="3" /> · Eric Bazile<fnref n="4" /> · Igor Esau<fnref n="5" />

Received: 26 March 2022 / Accepted: 26 October 2022 / Published online: 15 November 2022
© The Author(s) 2022

## Abstract

Forecast errors in near-surface temperatures are a persistent issue for numerical weather prediction models. A prominent example is warm biases during cloud-free, snow-covered nights. Many studies attribute these biases to parametrized processes such as turbulence or radiation. Here, we focus on the contribution of physical processes to the nocturnal temperature development. We compare model timestep output of individual tendencies from parametrized processes in the weather prediction model AROME-Arctic to measurements from Sodankylä, Finland. Thereby, we differentiate between the weakly stable boundary layer (wSBL) and the very stable boundary layer (vSBL) regimes. The wSBL is characterized by continuous turbulent exchange within the near-surface atmosphere, causing near-neutral temperature profiles. The vSBL is characterized by a decoupling of the lowermost model level, low turbulent exchange, and very stable temperature profiles. In our case study, both regimes occur simultaneously on small spatial scales of about 5 km. In addition, we demonstrate the model's sensitivity towards an updated surface treatment, allowing for faster surface cooling. The updated surface parametrization has profound impacts on parametrized processes in both regimes. However, only modelled temperatures in the vSBL are impacted substantially, whereas more efficient surface cooling in the wSBL is compensated by an increased turbulent heat transport within the boundary layer. This study demonstrates the utility of individual tendencies for understanding process-related differences between model configurations and emphasizes the need for model studies to distinguish between the wSBL and vSBL for reliable model verification.

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

<footnote>Marvin Kähnert
marvin.kahnert@uib.no</footnote>

<footnote n="1">¹ Geophysical Institute and Bjerknes Centre for Climate Research, University of Bergen, Postbox 7803, 5020 Bergen, Norway</footnote>

<footnote n="2">² Development Centre for Weather Forecasting, Norwegian Meteorological Institute, Oslo, Norway</footnote>

<footnote n="3">³ Finnish Meteorological Institute, Helsinki, Finland</footnote>

<footnote n="4">⁴ CNRM - UMR3589, Météo-France and CNRS, Toulouse, France</footnote>

<footnote n="5">⁵ Nansen Environmental and Remote Sensing Center,UiT - The Arctic University of Norway, Tromso, Norway</footnote>
**Keywords** Numerical weather prediction · Physical tendencies · Stability regimes · Stable boundary layer

# 1 Introduction

The accurate simulation of near-surface temperature variations in polar regions poses a long-standing challenge for numerical weather prediction (NWP) models (Viterbo et al. 1999; Holtslag et al. 2013; Sandu et al. 2013). During nocturnal, cloud-free periods over snow-covered surfaces, where the surface energy balance is negative, many NWP models underestimate observed drops in near-surface temperature. This underestimation becomes apparent as pronounced warm biases during model verification (Haiden et al. 2018; Køltzow et al. 2019). Previous studies point to shortcomings in the calculation of turbulent heat fluxes, radiative fluxes, and ground heat fluxes as possible causes (Sodemann and Foken 2004; Tjernström et al. 2005; Vihma et al. 2014). Esau et al. (2018) proposed that the heat capacity of a too deep boundary layer could cause too much thermal inertia, preventing a sufficiently fast reaction of near-surface temperatures. Furthermore, there may be contributions from an inappropriate treatment of reflected radiative fluxes in the radiation parametrization (Edwards 2009a, b), and, during the presence of snow, an erroneous representation of the strong thermal gradients within the snow pack (Arduini et al. 2019; Day et al. 2020). The complex interplay of these many different processes at the interface between the lower atmosphere and the surface has so far prevented an identification of the leading causes behind such NWP model deficiencies.

The problem is further aggravated by the possible occurrence of different stability regimes in the SBL that exhibit diverse characteristics between radiative cooling, stratification, and turbulence (Vogelezang and Holtslag 1996; Derbyshire 1999). Zilitinkevich et al. (2008) distinguishes between the so-called weakly stable boundary layer (wSBL) regime, and a very stable boundary layer (vSBL) regime:

1. The wSBL is characterized by a negative feedback between turbulence and radiative cooling. In this regime, the radiative surface cooling is balanced by turbulent downward transport of sensible heat from higher up in the boundary layer towards the surface.

2. The vSBL is characterized by a positive feedback between turbulence and radiative cooling. In this regime, the radiative surface cooling is not sufficiently counteracted by turbulent transport. Instead, increased atmospheric stability leads to a reduction or cessation of turbulence, resulting in further rapid cooling of the surface layer.

Both regimes exhibit considerably different observed temperature profiles. While the wSBL shows near-neutral temperature profiles, the strong stratification in the vSBL creates pronounced temperature inversions (André and Mahrt 1982; Vignon et al. 2017). The presence or absence of mechanical forcing appears to be a key factor in developing either regime. While the wSBL is often associated with a minimal wind speed, the vSBL is common in low-wind conditions (Ulden and Holtslag 1985; Vogelezang and Holtslag 1996). From idealized model simulations, McNider et al. (1995) found a bifurcation in the phase space of the SBL. Depending on small changes in initial conditions, feedback mechanisms can lead to a nearly discontinuous transition from the wSBL to the vSBL regime (Derbyshire 1999). In another idealized study, de Wiel et al. (2017) not only found sudden regime transitions in the nocturnal boundary layer, but also a coexistence of the vSBL and wSBL at certain wind speed ranges. In the context of the above-mentioned challenges in representing the SBL in NWP models, such chaotic behaviour of SBL regimes adds an additional level of complexity.
**Keywords** Numerical weather prediction · Physical tendencies · Stability regimes · Stable boundary layer

# 1 Introduction

The accurate simulation of near-surface temperature variations in polar regions poses a long-standing challenge for numerical weather prediction (NWP) models (Viterbo et al. 1999; Holtslag et al. 2013; Sandu et al. 2013). During nocturnal, cloud-free periods over snow-covered surfaces, where the surface energy balance is negative, many NWP models underestimate observed drops in near-surface temperature. This underestimation becomes apparent as pronounced warm biases during model verification (Haiden et al. 2018; Køltzow et al. 2019). Previous studies point to shortcomings in the calculation of turbulent heat fluxes, radiative fluxes, and ground heat fluxes as possible causes (Sodemann and Foken 2004; Tjernström et al. 2005; Vihma et al. 2014). Esau et al. (2018) proposed that the heat capacity of a too deep boundary layer could cause too much thermal inertia, preventing a sufficiently fast reaction of near-surface temperatures. Furthermore, there may be contributions from an inappropriate treatment of reflected radiative fluxes in the radiation parametrization (Edwards 2009a, b), and, during the presence of snow, an erroneous representation of the strong thermal gradients within the snow pack (Arduini et al. 2019; Day et al. 2020). The complex interplay of these many different processes at the interface between the lower atmosphere and the surface has so far prevented an identification of the leading causes behind such NWP model deficiencies.

The problem is further aggravated by the possible occurrence of different stability regimes in the SBL that exhibit diverse characteristics between radiative cooling, stratification, and turbulence (Vogelezang and Holtslag 1996; Derbyshire 1999). Zilitinkevich et al. (2008) distinguishes between the so-called weakly stable boundary layer (wSBL) regime, and a very stable boundary layer (vSBL) regime:

1. The wSBL is characterized by a negative feedback between turbulence and radiative cooling. In this regime, the radiative surface cooling is balanced by turbulent downward transport of sensible heat from higher up in the boundary layer towards the surface.

2. The vSBL is characterized by a positive feedback between turbulence and radiative cooling. In this regime, the radiative surface cooling is not sufficiently counteracted by turbulent transport. Instead, increased atmospheric stability leads to a reduction or cessation of turbulence, resulting in further rapid cooling of the surface layer.

Both regimes exhibit considerably different observed temperature profiles. While the wSBL shows near-neutral temperature profiles, the strong stratification in the vSBL creates pronounced temperature inversions (André and Mahrt 1982; Vignon et al. 2017). The presence or absence of mechanical forcing appears to be a key factor in developing either regime. While the wSBL is often associated with a minimal wind speed, the vSBL is common in low-wind conditions (Ulden and Holtslag 1985; Vogelezang and Holtslag 1996). From idealized model simulations, McNider et al. (1995) found a bifurcation in the phase space of the SBL. Depending on small changes in initial conditions, feedback mechanisms can lead to a nearly discontinuous transition from the wSBL to the vSBL regime (Derbyshire 1999). In another idealized study, de Wiel et al. (2017) not only found sudden regime transitions in the nocturnal boundary layer, but also a coexistence of the vSBL and wSBL at certain wind speed ranges. In the context of the above-mentioned challenges in representing the SBL in NWP models, such chaotic behaviour of SBL regimes adds an additional level of complexity.
Representing SBL regimes at the correct time and location can be decisive for meaningful model forecasts and verification studies.

So far, process representations in both the wSBL and vSBL have been addressed within idealized frameworks and in single-column model (SCM) experiments. SCM studies show that the combination of parametrizations employed in operational NWP models is principally capable of representing both the wSBL and the vSBL (Krishna et al. 2003; Steeneveld et al. 2006). Using 11 years of daily simulations over the Cabauw observatory, Baas et al. (2018) demonstrate a realistic response of the RACMO SCM to different mechanical forcing. Thereby, the turbulence scheme predominately balances the radiative loss at the surface in the wSBL, whereas in the vSBL, that role is taken over by the soil heat flux, preventing excessive surface cooling. With SCMs, it is however not possible to assess the impact of spatial variation due to surface heterogeneity, topography, and advection between neighbouring model columns on simulated SBL regimes.

Here, we attempt to disentangle the contributions to nocturnal near-surface temperature variations from different model components in the operational, three-dimensional limited-area NWP model AROME-Arctic (Müller et al. 2017). The application of a three-dimensional model distinguishes our study from previous SCM work, as the effects of dynamics, surface heterogeneity, and topography are included. In our analysis, we focus on a case study for a 3-day period from 14 to 17 March 2018 during the Year of Polar Prediction (YOPP) Special Observing Period Northern Hemisphere 1 (SOP-NH1, Køltzow et al. 2019). During this period, near-surface temperatures over the sub-Arctic supersite Sodankylä, Finland, showed pronounced differences between simulated and observed values. To account for the rapid transitions between the wSBL and vSBL in the model, we require model output at a higher time resolution than the usual hourly intervals. Such output is enabled here using the tool diagnostics par domaines horizontaux (DDH, Météo-France, 2019), which provides the output of selected variables in a limited domain at every model integration timestep. To identify the contribution of different parametrized processes to the temperature evolution in the NWP model, we in addition employ the output of individual physics parametrization tendencies from AROME-Arctic (Kähnert et al. 2021). This unique combination of two diagnostics provides us unprecedented insight into the SBL development in a three-dimensional NWP model. Using a sensitivity study with a new scheme for surface processes considered for a future operational cycle, we study the interplay of resolved and parametrized processes during SBL development in time and space.

# 2 Data and Method

## 2.1 AROME-Arctic

AROME-Arctic is an operational, convection-permitting forecasting system covering the European Arctic with a horizontal grid spacing of 2.5 km (Müller et al. 2017). In the vertical, the model uses 65 sigma-hybrid pressure levels, with the lowermost two mass levels at approximately 11 m and 32 m. The model employs a semi-Lagrangian advection scheme and several advanced physical parametrizations (Bengtsson et al. 2017). The most important parametrizations for a clear-night SBL are described in more detail below.

Radiation is calculated by the rapid radiative transfer model (RRTM; Fouquart and Bonnel 1980; Mlawer et al. 1997), providing radiative fluxes to the atmospheric column and to the surface. Six spectral bands are used for shortwave radiation, while climatological distribu-
Representing SBL regimes at the correct time and location can be decisive for meaningful model forecasts and verification studies.

So far, process representations in both the wSBL and vSBL have been addressed within idealized frameworks and in single-column model (SCM) experiments. SCM studies show that the combination of parametrizations employed in operational NWP models is principally capable of representing both the wSBL and the vSBL (Krishna et al. 2003; Steeneveld et al. 2006). Using 11 years of daily simulations over the Cabauw observatory, Baas et al. (2018) demonstrate a realistic response of the RACMO SCM to different mechanical forcing. Thereby, the turbulence scheme predominantly balances the radiative loss at the surface in the wSBL, whereas in the vSBL, that role is taken over by the soil heat flux, preventing excessive surface cooling. With SCMs, it is however not possible to assess the impact of spatial variation due to surface heterogeneity, topography, and advection between neighbouring model columns on simulated SBL regimes.

Here, we attempt to disentangle the contributions to nocturnal near-surface temperature variations from different model components in the operational, three-dimensional limited-area NWP model AROME-Arctic (Müller et al. 2017). The application of a three-dimensional model distinguishes our study from previous SCM work, as the effects of dynamics, surface heterogeneity, and topography are included. In our analysis, we focus on a case study for a 3-day period from 14 to 17 March 2018 during the Year of Polar Prediction (YOPP) Special Observing Period Northern Hemisphere 1 (SOP-NH1, Køltzow et al. 2019). During this period, near-surface temperatures over the sub-Arctic supersite Sodankylä, Finland, showed pronounced differences between simulated and observed values. To account for the rapid transitions between the wSBL and vSBL in the model, we require model output at a higher time resolution than the usual hourly intervals. Such output is enabled here using the tool diagnostics par domaines horizontaux (DDH, Météo-France, 2019), which provides the output of selected variables in a limited domain at every model integration timestep. To identify the contribution of different parametrized processes to the temperature evolution in the NWP model, we in addition employ the output of individual physics parametrization tendencies from AROME-Arctic (Kähnert et al. 2021). This unique combination of two diagnostics provides us unprecedented insight into the SBL development in a three-dimensional NWP model. Using a sensitivity study with a new scheme for surface processes considered for a future operational cycle, we study the interplay of resolved and parametrized processes during SBL development in time and space.

# 2 Data and Method

## 2.1 AROME-Arctic

AROME-Arctic is an operational, convection-permitting forecasting system covering the European Arctic with a horizontal grid spacing of 2.5 km (Müller et al. 2017). In the vertical, the model uses 65 sigma-hybrid pressure levels, with the lowermost two mass levels at approximately 11 m and 32 m. The model employs a semi-Lagrangian advection scheme and several advanced physical parametrizations (Bengtsson et al. 2017). The most important parametrizations for a clear-night SBL are described in more detail below.

Radiation is calculated by the rapid radiative transfer model (RRTM; Fouquart and Bonnel 1980; Mlawer et al. 1997), providing radiative fluxes to the atmospheric column and to the surface. Six spectral bands are used for shortwave radiation, while climatological distribu-
tions of aerosols and ozone are utilized for the longwave radiation. RRTM's contribution to modelled atmospheric temperatures stems from the vertical divergence of radiative fluxes between the model levels (see Sect. 2.2, Eq. 4).

The turbulence parametrization HARMONIE-AROME with RACMO Turbulence (HARATU; Lenderink and Holtslag 2004; Bengtsson et al. 2017) employs an eddy diffusivity approach. Here, the vertical turbulent heat flux $\overline{w'\theta'}$ is taken proportionally to an eddy diffusivity coefficient for heat $K_h$, and the vertical gradient of potential temperature $\theta$:

$$ \overline{w'\theta'} = -K_h \frac{\partial\theta}{\partial z}. \qquad (1) $$

To determine $K_h$, HARATU employs a 1.5th-order closure, combining a prognostic equation for turbulence kinetic energy (TKE, $E_k$) with a diagnostic length scale $l$:

$$ K_h = \sqrt{E_k} \cdot l. \qquad (2) $$

For stable conditions, this length scale is a function of local stability, expressed as the sum of inverses of a neutral $l_n$ and a stable length scale $l_s$ (Deardorff 1980; Baas et al. 2008):

$$ \frac{1}{l} = \frac{1}{l_n} + \frac{1}{l_s} = \frac{1}{c_n \kappa z} + \frac{N}{c_s \sqrt{E_k}}. \qquad (3) $$

Here, $c_n$ and $c_s$ are model internal parameters, $\kappa$ is the von Kármán constant, $z$ is height above ground, and $N$ is the Brunt-Väisälä frequency. The turbulence parametrization affects atmospheric temperatures by reducing the vertical temperature gradient between adjacent model levels. Hereby, stability and TKE affect the speed at which the gradient can be reduced. The HARATU parametrization has been compared to results from the first and second GEWEX Atmospheric Boundary Layer Studies (GABLS, Holtslag 2006; Svensson and Holtslag 2006) and yielded good agreements (Baas et al. 2008; de Rooy et al. 2022). Further, the SCM experiments by Baas et al. (2018) demonstrate that the parametrization is capable of representing both the wSBL and vSBL.

Regarding the vSBL regime, the relatively coarse vertical resolution of the lowermost model level can be in conflict with the metre-scale of vertical eddies, the anisotropy of turbulence, and the occasional very shallow depth of the SBL (Zilitinkevich and Esau 2005; Vihma et al. 2014). As AROME-Arctic does not have a dedicated setting for very stable situations, the same set of schemes needs to represent the feedback mechanism connected to the wSBL and vSBL, respectively (Sect. 1). Even though configurations of the model with higher vertical resolution exist, we want our results to relate to the operational forecasts and thus kept the set-up with 65 vertical levels.

The surface parametrization (SURFEX; Moigne 2009) calculates, among others, surface temperatures, surface sensible, and latent heat fluxes, as well as wind stress. The parametrization applies specific formulations depending on the underlying surface type and distinguishes between sea, nature, urban areas, and water. These calculations also include the treatment of snow cover and its effects on, e.g. albedo and emissivity. A more detailed depiction of important process representations in regard to this study, namely the treatment of the vegetation and snow cover, is provided in Sect. 2.3.1.

All model simulations are initialized at 1200 UTC 14 March 2018 and are integrated for 72 h. To allow for some adjustment of the soil temperatures, a 4-day spin-up period is chosen. 3D-VAR data assimilation of conventional observations is applied to upper-air variables, whereas observed screen level temperature, relative humidity, and snow depth are assimilated by optimal interpolation and used to update surface variables. Boundary
tions of aerosols and ozone are utilized for the longwave radiation. RRTM's contribution to modelled atmospheric temperatures stems from the vertical divergence of radiative fluxes between the model levels (see Sect. 2.2, Eq. 4).

The turbulence parametrization HARMONIE-AROME with RACMO Turbulence (HARATU; Lenderink and Holtslag 2004; Bengtsson et al. 2017) employs an eddy diffusivity approach. Here, the vertical turbulent heat flux $\overline{w'\theta'}$ is taken proportionally to an eddy diffusivity coefficient for heat $K_h$, and the vertical gradient of potential temperature $\theta$:

$$ \overline{w'\theta'} = -K_h \frac{\partial\theta}{\partial z}. \qquad (1) $$

To determine $K_h$, HARATU employs a 1.5th-order closure, combining a prognostic equation for turbulence kinetic energy (TKE, $E_k$) with a diagnostic length scale $l$:

$$ K_h = \sqrt{E_k} \cdot l. \qquad (2) $$

For stable conditions, this length scale is a function of local stability, expressed as the sum of inverses of a neutral $l_n$ and a stable length scale $l_s$ (Deardorff 1980; Baas et al. 2008):

$$ \frac{1}{l} = \frac{1}{l_n} + \frac{1}{l_s} = \frac{1}{c_n \kappa z} + \frac{N}{c_s \sqrt{E_k}}. \qquad (3) $$

Here, $c_n$ and $c_s$ are model internal parameters, $\kappa$ is the von Kármán constant, $z$ is height above ground, and $N$ is the Brunt-Väisälä frequency. The turbulence parametrization affects atmospheric temperatures by reducing the vertical temperature gradient between adjacent model levels. Hereby, stability and TKE affect the speed at which the gradient can be reduced. The HARATU parametrization has been compared to results from the first and second GEWEX Atmospheric Boundary Layer Studies (GABLS, Holtslag 2006; Svensson and Holtslag 2006) and yielded good agreements (Baas et al. 2008; de Rooy et al. 2022). Further, the SCM experiments by Baas et al. (2018) demonstrate that the parametrization is capable of representing both the wSBL and vSBL.

Regarding the vSBL regime, the relatively coarse vertical resolution of the lowermost model level can be in conflict with the metre-scale of vertical eddies, the anisotropy of turbulence, and the occasional very shallow depth of the SBL (Zilitinkevich and Esau 2005; Vihma et al. 2014). As AROME-Arctic does not have a dedicated setting for very stable situations, the same set of schemes needs to represent the feedback mechanism connected to the wSBL and vSBL, respectively (Sect. 1). Even though configurations of the model with higher vertical resolution exist, we want our results to relate to the operational forecasts and thus kept the set-up with 65 vertical levels.

The surface parametrization (SURFEX; Moigne 2009) calculates, among others, surface temperatures, surface sensible, and latent heat fluxes, as well as wind stress. The parametrization applies specific formulations depending on the underlying surface type and distinguishes between sea, nature, urban areas, and water. These calculations also include the treatment of snow cover and its effects on, e.g. albedo and emissivity. A more detailed depiction of important process representations in regard to this study, namely the treatment of the vegetation and snow cover, is provided in Sect. 2.3.1.

All model simulations are initialized at 1200 UTC 14 March 2018 and are integrated for 72 h. To allow for some adjustment of the soil temperatures, a 4-day spin-up period is chosen. 3D-VAR data assimilation of conventional observations is applied to upper-air variables, whereas observed screen level temperature, relative humidity, and snow depth are assimilated by optimal interpolation and used to update surface variables. Boundary
conditions are provided by the high-resolution version of the global ECMWF Integrated Forecasting System (IFS-HRES).

## 2.2 Diagnostic Methods

We utilize a combination of two diagnostics to assess the contribution of physical processes in AROME-Arctic during the SBL, namely the individual tendency output and DDH. The individual tendency output provides the contributions of each parametrization to the prognostic variables such as temperature, specific humidity, and momentum. A detailed description of this diagnostic can be found in Kähnert et al. (2021). To better frame this output, let us look at the evolution of temperature in the (dry) stable boundary layer in AROME-Arctic:

$$ \frac{\partial T}{\partial t} = -\underbrace{u_H \Delta H T}_{\text{advection}} - w \underbrace{\left(\frac{\partial T}{\partial z} + \frac{g}{c_p}\right)}_{\text{subsidence}} - \underbrace{\frac{1}{\rho c_p} \frac{\partial}{\partial z} (F^+ - F^-)}_{\text{radiation}} + \underbrace{\frac{1}{\rho} \frac{\partial}{\partial z} \left(\rho K_h \frac{\partial \theta}{\partial z}\right)}_{\text{turbulence}} \quad (4) $$

Here, $u_H$ depicts the wind speed, $w$—the vertical velocity component, $g$—the acceleration due to gravity, $c_p$—the heat capacity of air, $\partial(F^+ - F^-)/\partial z$—the divergence of the net upward radiative flux, and $\rho$—the density of air. The first two terms on the right-hand-side of Eq. 4 are resolved by the model dynamics and are incorporated in the dynamical tendency, including horizontal diffusion (Kähnert et al. 2021). The latter two terms represent the model physics and are calculated by the respective parametrizations. The surface parametrization is not represented by an individual term in Eq. 4, instead its contribution in form of turbulent fluxes and radiative fluxes influences the turbulence and the radiation parametrizations (Seity et al. 2011).

Typically, the individual tendency output provides accumulated values over the (hourly) output interval of the model. However, we require high-temporal resolution to account for the rapid transitions in the SBL (see Sect. 1). Therefore, we combine this diagnostic with the high-frequency output provided by DDH.

DDH enables the output of selected variables for every model timestep (here 75 s) during a full model run. This output can be constrained to a user-defined sub-domain of the model. In our case, we chose an area of 35 km × 35 km centred around the supersite of Sodankylä (67.36 N, 26.64 E), a total of 196 grid boxes. DDH provides output for each of these grid boxes, enabling us to investigate the contribution and interplay of resolved and parametrized processes on a very detailed temporal and spatial resolution.

## 2.3 Experiments, Measurement Site, and Case Study

### 2.3.1 Sensitivity Experiments

The representation of the surface is of particular importance for modelling the SBL (de Wiel et al. 2017; Baas et al. 2018), especially in the presence of snow (Sterk et al. 2013; Day et al. 2020). We address this importance by two experiments, named REF and NEW. A reference run (REF) utilizes the model configuration described in Sect. 2.1. Sensitivity experiment NEW, on the other hand, contains a set of sophisticated updates for the surface parametrization. These updates are detailed below, contrasted against the configuration in REF. A schematic depiction of both configurations is given in Fig. 1.
conditions are provided by the high-resolution version of the global ECMWF Integrated Forecasting System (IFS-HRES).

## 2.2 Diagnostic Methods

We utilize a combination of two diagnostics to assess the contribution of physical processes in AROME-Arctic during the SBL, namely the individual tendency output and DDH. The individual tendency output provides the contributions of each parametrization to the prognostic variables such as temperature, specific humidity, and momentum. A detailed description of this diagnostic can be found in Kähnert et al. (2021). To better frame this output, let us look at the evolution of temperature in the (dry) stable boundary layer in AROME-Arctic:

$$ \frac{\partial T}{\partial t} = -\underbrace{u_H \Delta H T}_{\text{advection}} - \underbrace{w \left( \frac{\partial T}{\partial z} + \frac{g}{c_p} \right)}_{\text{subsidence}} - \underbrace{\frac{1}{\rho c_p} \frac{\partial}{\partial z} (F^+ - F^-)}_{\text{radiation}} + \underbrace{\frac{1}{\rho} \frac{\partial}{\partial z} \left( \rho K_h \frac{\partial \theta}{\partial z} \right)}_{\text{turbulence}}. \quad (4) $$

Here, $u_H$ depicts the wind speed, $w$—the vertical velocity component, $g$—the acceleration due to gravity, $c_p$—the heat capacity of air, $\partial(F^+ - F^-)/\partial z$—the divergence of the net upward radiative flux, and $\rho$—the density of air. The first two terms on the right-hand-side of Eq. 4 are resolved by the model dynamics and are incorporated in the dynamical tendency, including horizontal diffusion (Kähnert et al. 2021). The latter two terms represent the model physics and are calculated by the respective parametrizations. The surface parametrization is not represented by an individual term in Eq. 4, instead its contribution in form of turbulent fluxes and radiative fluxes influences the turbulence and the radiation parametrizations (Seity et al. 2011).

Typically, the individual tendency output provides accumulated values over the (hourly) output interval of the model. However, we require high-temporal resolution to account for the rapid transitions in the SBL (see Sect. 1). Therefore, we combine this diagnostic with the high-frequency output provided by DDH.

DDH enables the output of selected variables for every model timestep (here 75 s) during a full model run. This output can be constrained to a user-defined sub-domain of the model. In our case, we chose an area of 35 km × 35 km centred around the supersite of Sodankylä (67.36 N, 26.64 E), a total of 196 grid boxes. DDH provides output for each of these grid boxes, enabling us to investigate the contribution and interplay of resolved and parametrized processes on a very detailed temporal and spatial resolution.

## 2.3 Experiments, Measurement Site, and Case Study

### 2.3.1 Sensitivity Experiments

The representation of the surface is of particular importance for modelling the SBL (de Wiel et al. 2017; Baas et al. 2018), especially in the presence of snow (Sterk et al. 2013; Day et al. 2020). We address this importance by two experiments, named REF and NEW. A reference run (REF) utilizes the model configuration described in Sect. 2.1. Sensitivity experiment NEW, on the other hand, contains a set of sophisticated updates for the surface parametrization. These updates are detailed below, contrasted against the configuration in REF. A schematic depiction of both configurations is given in Fig. 1.
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Idealized schematic of the surface parametrization in the two experiments REF and NEW. See text for details</figcaption></figure>

In REF, the soil is represented by two levels and its temperature is calculated according to a force-restore method (Noilhan and Planton 1989). Hereby, the upper level describes a composite layer consisting of soil, vegetation, and snow (Boone and Etchevers 2001). In NEW, ground heat fluxes and soil temperatures are explicitly calculated by a multi-layer soil parametrization (Decharme et al. 2011). In contrast to the force-restore method, such diffusive approaches are able to capture the near-term effects of, i.e. sudden changes in cloudiness or day-night transitions on soil temperatures or ground heat fluxes (van der Linden et al. 2022). This modification, however, is not of primary importance during our case study due to the presence of a deep snow layer that insulates the ground (Sect. 2.3.3).

The treatment of the energy balance between the surface and the atmosphere is also modified. In REF, only one energy balance is considered for the entire ground-vegetation-snow system. The composite soil layer provides heat fluxes that are weighted, based on the respective fractional coverage of bare ground, snow, and vegetation. Hereby, the model distinguishes between high and low vegetation, expressed by, e.g. different roughness lengths, indicated as a difference in elevation of the respective patches in Fig. 1a. These roughness lengths influence, among others, the estimation of total snow cover fraction by the snow scheme (see below), which in turn influences the surface emissivity.

NEW, in contrast, employs a multi-energy-budget explicit vegetation parametrization (ISBA-MEB; Boone et al. 2017; Napoly et al. 2017) that considers the fully coupled energy budgets between the soil, the snow, and the vegetation. Hereby, high vegetation is represented by a canopy layer parametrization (Fig. 1b) to better capture the resistance between the ground and the canopy air space (Boone et al. 2017). For grid boxes covered by forest (see Sect. 2.3.2), this enhances the estimated roughness length by about 0.16 m on average. Furthermore, ISBA-MEB improves the calculation of total snow cover fraction. Now the vegetation canopy can (partly) be buried by snow, influencing the exchange between snow-covered surfaces and the atmosphere. Among others, ISBA-MEB updates the calculation of surface emissivity by taking the absorption of longwave radiation within the canopy and the effect of partially buried vegetation into account. For the domain considered here, ISBA-MEB, together with the new snow parametrization (see below), reduces the effective surface emissivity by about 0.01 on average and enhances the snow cover fraction by 0.4.

Finally, the treatment of snow-covered surfaces is changed. In REF, snow is calculated by a single-layer snow parametrization (D95; Douville et al. 1995a, b) which estimates the
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Idealized schematic of the surface parametrization in the two experiments REF and NEW. See text for details</figcaption></figure>

In REF, the soil is represented by two levels and its temperature is calculated according to a force-restore method (Noilhan and Planton 1989). Hereby, the upper level describes a composite layer consisting of soil, vegetation, and snow (Boone and Etchevers 2001). In NEW, ground heat fluxes and soil temperatures are explicitly calculated by a multi-layer soil parametrization (Decharme et al. 2011). In contrast to the force-restore method, such diffusive approaches are able to capture the near-term effects of, i.e. sudden changes in cloudiness or day-night transitions on soil temperatures or ground heat fluxes (van der Linden et al. 2022). This modification, however, is not of primary importance during our case study due to the presence of a deep snow layer that insulates the ground (Sect. 2.3.3).

The treatment of the energy balance between the surface and the atmosphere is also modified. In REF, only one energy balance is considered for the entire ground-vegetation-snow system. The composite soil layer provides heat fluxes that are weighted, based on the respective fractional coverage of bare ground, snow, and vegetation. Hereby, the model distinguishes between high and low vegetation, expressed by, e.g. different roughness lengths, indicated as a difference in elevation of the respective patches in Fig. 1a. These roughness lengths influence, among others, the estimation of total snow cover fraction by the snow scheme (see below), which in turn influences the surface emissivity.

NEW, in contrast, employs a multi-energy-budget explicit vegetation parametrization (ISBA-MEB; Boone et al. 2017; Napoly et al. 2017) that considers the fully coupled energy budgets between the soil, the snow, and the vegetation. Hereby, high vegetation is represented by a canopy layer parametrization (Fig. 1b) to better capture the resistance between the ground and the canopy air space (Boone et al. 2017). For grid boxes covered by forest (see Sect. 2.3.2), this enhances the estimated roughness length by about 0.16 m on average. Furthermore, ISBA-MEB improves the calculation of total snow cover fraction. Now the vegetation canopy can (partly) be buried by snow, influencing the exchange between snow-covered surfaces and the atmosphere. Among others, ISBA-MEB updates the calculation of surface emissivity by taking the absorption of longwave radiation within the canopy and the effect of partially buried vegetation into account. For the domain considered here, ISBA-MEB, together with the new snow parametrization (see below), reduces the effective surface emissivity by about 0.01 on average and enhances the snow cover fraction by 0.4.

Finally, the treatment of snow-covered surfaces is changed. In REF, snow is calculated by a single-layer snow parametrization (D95; Douville et al. 1995a, b) which estimates the
evolution of the equivalent water content and temperature of a bulk layer of snow. In NEW, a multi-layer snow parametrization (3-L; Boone and Etchevers 2001) enables an explicit treatment of the large thermal gradients that exist in the snow pack. The 3-L parametrization employs 12 layers, whose depths are typically a few cm, with finer layers towards the respective interfaces between air/snow and snow/ground (Decharme et al. 2016).<fnref n="1" />

The thermal properties of snow take snow density and the latent heat transfer during metamorphism into account (Yen 1981; Sun et al. 1999; Barrere et al. 2017). For the case studied here, this yields an approximate thermal diffusivity of $3 \times 10^{-7} \text{ m}^2 \text{ s}^{-1}$. Further, effects of storing capabilities of liquid water in the snow are included, as well as the absorption of incident radiation within the snow pack. Compared to the single-layer parametrization D95 employed in REF, the 3-L scheme improves the representation of snow depth, snow water equivalent, and surface temperature.

All of these changes in NEW influence the coupling between the surface and the atmosphere. Especially, the turbulence and radiation parametrizations are expected to be impacted, which in turn can affect the representation of the wSBL or vSBL (Sect. 1).

### 2.3.2 Sodankylä Supersite

The supersite of Sodankylä is located in the boreal and sub-Arctic zone (67.37°N, 26.63°E, Fig. 3a, red dot). The site provides a comprehensive set of measurements, well-suited for atmospheric and environmental research. The available observations cover conventional meteorological data, heat and carbon fluxes, remote sensing observations, and snow pack properties. For our study, we mostly focus on the observations taken at the 48 m-high micrometeorological mast. An extensive instrumentation provides temperature, wind, humidity, and radiation measurements at various levels between 8 and 48 m, post-processed to a consistent temporal resolution of 10 min. For a detailed description of the mast, see Kangas et al. (2016).

The Sodankylä supersite is located within predominately flat terrain, reflected by the topography in AROME-Arctic with elevations of 175–350 m above sea level (Fig. 2a). The micrometeorological mast (black cross) is located within a slight surface depression in that region. During the investigated period, a snow-layer depth of 74–82 cm is present over the entire domain (Fig. 2b). The most common land cover type in the region is the boreal needle-leaf forest (Fig. 2c), influencing the estimated roughness length (Fig. 2d).

### 2.3.3 Case Study

A selected 3-day period from 14 to 17 March 2018 during the YOPP Special Observing Period Northern Hemisphere 1 (SOP-NH1, Køltzow et al. 2019) forms our case study. The weather situation in all three nights during the case study was dominated by a high-pressure system, leading to subsiding air masses that dissolved clouds over Scandinavia (Fig. 3a). In the first two nights, the observations show a decrease in near-surface temperatures by up to 15 K, reaching values of –20 °C at 3 m (Fig. 3b). Wind speeds are comparably low during the nocturnal periods, dropping from 3 to 1 m s⁻<fnref n="1" /> (Fig. 3c). During the third night (17 March 2018), a warm front passed over Sodankylä, causing warmer temperatures of up to 5 K at 3 m and 8 m compared to the previous nights (Fig. 3b) and even led to a halt of the nocturnal temperature decreases at 18 m, 32 m, and 48 m. The front can be identified by the brief period of low wind speed, as the wind direction changed. Due to the difference in

<footnote n="1">¹ For a snow-layer depth of 78 cm (see Sect. 2.3.2), the layer depths from top to bottom are 1 cm, 5 cm, 6.5 cm, 6.5 cm, 6.5 cm, 9.3 cm, 12.4 cm, 9.3 cm, 6.5 cm, 6.5 cm, 6.5 cm, and 2 cm.</footnote>
evolution of the equivalent water content and temperature of a bulk layer of snow. In NEW, a multi-layer snow parametrization (3-L; Boone and Etchevers 2001) enables an explicit treatment of the large thermal gradients that exist in the snow pack. The 3-L parametrization employs 12 layers, whose depths are typically a few cm, with finer layers towards the respective interfaces between air/snow and snow/ground (Decharme et al. 2016).<fnref n="1" />

The thermal properties of snow take snow density and the latent heat transfer during metamorphism into account (Yen 1981; Sun et al. 1999; Barrere et al. 2017). For the case studied here, this yields an approximate thermal diffusivity of $3 \times 10^{-7} \text{ m}^2 \text{ s}^{-1}$. Further, effects of storing capabilities of liquid water in the snow are included, as well as the absorption of incident radiation within the snow pack. Compared to the single-layer parametrization D95 employed in REF, the 3-L scheme improves the representation of snow depth, snow water equivalent, and surface temperature.

All of these changes in NEW influence the coupling between the surface and the atmosphere. Especially, the turbulence and radiation parametrizations are expected to be impacted, which in turn can affect the representation of the wSBL or vSBL (Sect. 1).

### 2.3.2 Sodankylä Supersite

The supersite of Sodankylä is located in the boreal and sub-Arctic zone (67.37°N, 26.63°E, Fig. 3a, red dot). The site provides a comprehensive set of measurements, well-suited for atmospheric and environmental research. The available observations cover conventional meteorological data, heat and carbon fluxes, remote sensing observations, and snow pack properties. For our study, we mostly focus on the observations taken at the 48 m-high micrometeorological mast. An extensive instrumentation provides temperature, wind, humidity, and radiation measurements at various levels between 8 and 48 m, post-processed to a consistent temporal resolution of 10 min. For a detailed description of the mast, see Kangas et al. (2016).

The Sodankylä supersite is located within predominately flat terrain, reflected by the topography in AROME-Arctic with elevations of 175–350 m above sea level (Fig. 2a). The micrometeorological mast (black cross) is located within a slight surface depression in that region. During the investigated period, a snow-layer depth of 74–82 cm is present over the entire domain (Fig. 2b). The most common land cover type in the region is the boreal needle-leaf forest (Fig. 2c), influencing the estimated roughness length (Fig. 2d).

### 2.3.3 Case Study

A selected 3-day period from 14 to 17 March 2018 during the YOPP Special Observing Period Northern Hemisphere 1 (SOP-NH1, Køltzow et al. 2019) forms our case study. The weather situation in all three nights during the case study was dominated by a high-pressure system, leading to subsiding air masses that dissolved clouds over Scandinavia (Fig. 3a). In the first two nights, the observations show a decrease in near-surface temperatures by up to 15 K, reaching values of −20 °C at 3 m (Fig. 3b). Wind speeds are comparably low during the nocturnal periods, dropping from 3 to 1 m s⁻<fnref n="1" /> (Fig. 3c). During the third night (17 March 2018), a warm front passed over Sodankylä, causing warmer temperatures of up to 5 K at 3 m and 8 m compared to the previous nights (Fig. 3b) and even led to a halt of the nocturnal temperature decreases at 18 m, 32 m, and 48 m. The front can be identified by the brief period of low wind speed, as the wind direction changed. Due to the difference in

<footnote n="1">¹ For a snow-layer depth of 78 cm (see Sect. 2.3.2), the layer depths from top to bottom are 1 cm, 5 cm, 6.5 cm, 6.5 cm, 6.5 cm, 9.3 cm, 12.4 cm, 9.3 cm, 6.5 cm, 6.5 cm, 6.5 cm, and 2 cm.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Physiography of AROME-Arctic for the extracted 35 × 35 km² model domain investigated in this study. The black cross indicates the location of the micrometeorological mast. Shown are **a** topography, **b** snow depth, **c** area fraction from the most abundant land cover type, in this case boreal needleleaf forest, and **d** roughness length. In sensitivity experiment NEW, the roughness length is increased on average by 0.16 m in the depicted domain, yet the spatial pattern is not affected. The location of the micrometeorological mast is marked by the black cross</figcaption></figure>

forcing associated with such a frontal passing, the remaining analysis focuses solely on the first two nights of 15–16 March 2018, dominated by the radiative surface cooling conducive to SBL development.

# 3 Results

## 3.1 The Nocturnal Boundary Layer in the Reference Run

We now utilize the combination of individual tendencies and DDH (Sect. 2.2) to inspect the contributions of parametrized and resolved processes during nocturnal, stable periods in AROME-Arctic. First, we compare the temperature development, observed at the micrometeorological mast to the closest grid point in our reference run REF (Fig. 4). We focus on the periods when the model internal boundary layer type (Kähnert et al. 2021) is stable, i.e. when a downward surface-buoyancy flux is present.
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Physiography of AROME-Arctic for the extracted 35 × 35 km² model domain investigated in this study. The black cross indicates the location of the micrometeorological mast. Shown are **a** topography, **b** snow depth, **c** area fraction from the most abundant land cover type, in this case boreal needleleaf forest, and **d** roughness length. In sensitivity experiment NEW, the roughness length is increased on average by 0.16 m in the depicted domain, yet the spatial pattern is not affected. The location of the micrometeorological mast is marked by the black cross</figcaption></figure>

forcing associated with such a frontal passing, the remaining analysis focuses solely on the first two nights of 15–16 March 2018, dominated by the radiative surface cooling conducive to SBL development.

# 3 Results

## 3.1 The Nocturnal Boundary Layer in the Reference Run

We now utilize the combination of individual tendencies and DDH (Sect. 2.2) to inspect the contributions of parametrized and resolved processes during nocturnal, stable periods in AROME-Arctic. First, we compare the temperature development, observed at the micrometeorological mast to the closest grid point in our reference run REF (Fig. 4). We focus on the periods when the model internal boundary layer type (Kähnert et al. 2021) is stable, i.e. when a downward surface-buoyancy flux is present.
<figure><img src="image_5.png" /><figcaption>**Fig. 3** a Synoptic situation at 0000 UTC 16 March 2018 depicting total cloud cover (shading) and mean sea level pressure in hPa (orange contours). Values depict 12 h forecasts by AROME-Arctic (upper domain) and the mesoscale ensemble prediction system (MEPS, lower domain) extracted from the operational archive of the Norwegian Meteorological Institute. The red dot marks the Sodankylä supersite. b Temperature observed at the micrometeorological mast at the Sodankylä supersite. c Wind speed observed at the micrometeorological mast at the Sodankylä supersite. Unfortunately, no wind speed estimates at lower heights are available in the vicinity of the micrometeorological mast</figcaption></figure>

We integrate the respective model tendencies coming from the individual parametrizations of turbulence and radiation, as well as the model dynamics (Eq. 4). This provides us with the change in temperature due to different processes in the SBL throughout each night, labelled $\partial T_{SBL}$. We further depict the contribution from the model physics as the sum of the contributions from the turbulence and radiation parametrizations.

First, we examine the development of modelled and observed temperatures (red and black, solid lines). During the first two nights, the nocturnal cooling in the model (7.5 K first night/10 K second night) is lower compared to the cooling in the observations (13 K in both nights). This development would yield a considerable warm bias of 5.5 K at the end of the first night and a warm bias of about 3 K at the end of the second night. As shown in the following paragraphs, the reduced bias between the first and second night coincides with a change in stability regime at the investigated grid point, from (mostly) weakly stable during the first night to very stable during the second night. This change is apparent by the contribution of the turbulence parametrization.

On both nights, the cooling due to turbulence starts similarly as the cold surface triggers an overall downward transport of heat in the PBL. However, on the first night, this cooling by turbulence comes to a halt after approximately 3 h (1800 UTC). This halt describes the compensation of the turbulent heat loss to the surface by an efficient transport of heat from upper levels in the SBL. Such an efficient heat transport is indicative of a strongly coupled SBL and the formation of a wSBL in the model. Later during the first night (2100 UTC),
<figure><img src="image_5.png" /><figcaption>**Fig. 3** a Synoptic situation at 0000 UTC 16 March 2018 depicting total cloud cover (shading) and mean sea level pressure in hPa (orange contours). Values depict 12 h forecasts by AROME-Arctic (upper domain) and the mesoscale ensemble prediction system (MEPS, lower domain) extracted from the operational archive of the Norwegian Meteorological Institute. The red dot marks the Sodankylä supersite. b Temperature observed at the micrometeorological mast at the Sodankylä supersite. c Wind speed observed at the micrometeorological mast at the Sodankylä supersite. Unfortunately, no wind speed estimates at lower heights are available in the vicinity of the micrometeorological mast</figcaption></figure>

We integrate the respective model tendencies coming from the individual parametrizations of turbulence and radiation, as well as the model dynamics (Eq. 4). This provides us with the change in temperature due to different processes in the SBL throughout each night, labelled $\partial T_{SBL}$. We further depict the contribution from the model physics as the sum of the contributions from the turbulence and radiation parametrizations.

First, we examine the development of modelled and observed temperatures (red and black, solid lines). During the first two nights, the nocturnal cooling in the model (7.5 K first night/10 K second night) is lower compared to the cooling in the observations (13 K in both nights). This development would yield a considerable warm bias of 5.5 K at the end of the first night and a warm bias of about 3 K at the end of the second night. As shown in the following paragraphs, the reduced bias between the first and second night coincides with a change in stability regime at the investigated grid point, from (mostly) weakly stable during the first night to very stable during the second night. This change is apparent by the contribution of the turbulence parametrization.

On both nights, the cooling due to turbulence starts similarly as the cold surface triggers an overall downward transport of heat in the PBL. However, on the first night, this cooling by turbulence comes to a halt after approximately 3 h (1800 UTC). This halt describes the compensation of the turbulent heat loss to the surface by an efficient transport of heat from upper levels in the SBL. Such an efficient heat transport is indicative of a strongly coupled SBL and the formation of a wSBL in the model. Later during the first night (2100 UTC),
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Temperature change $\partial T_{SBL}$ at 8 m observed at the micrometeorological mast at Sodankylä (red line) and modelled from the lowermost model level (11 m) of the closest model grid point in REF (black solid line) during stably stratified, nocturnal periods. The simulated $\partial T_{SBL}$ is decomposed into the contributions from the turbulence scheme (orange line), and radiation scheme (blue line), the model physics (black dashed line), and the model dynamics (black dotted line)</figcaption></figure>

this coupling reduces and the turbulent heat loss to the surface dominates again, leading to cooling by turbulence.

The second night is different. Here, turbulence cools consistently as the lowermost model level decouples from the rest of the PBL, hindering the turbulent exchange of heat with higher model levels, and a vSBL forms. The continued heat loss to the surface causes a strong vertical flux divergence leading to a pronounced cooling of 9 K by turbulence throughout the night (Fig. 4, orange line). As only a relatively shallow layer is cooled, temperatures drop stronger compared to the wSBL in the previous night.

The two nights also exhibit changes in the contribution of the radiation tendency. In the first night, the continuous radiative cooling amounts to about 2 K. During the second night, the radiative cooling stops around 2300 UTC, leading to a reduced cooling of 1 K throughout the night. This reduction is due to the formation of a strong inversion atop the lowest model level in the vSBL and the associated increase in downward longwave radiation (see Sect. 3.3, Fig. 10d).

Our diagnostic yields a detailed insight into the changed properties of the local, physical processes at the investigated grid point. However, no clear statement can be drawn from the differences in the dynamical tendency (Fig. 4, black dotted line), as both small- and large-scale advection (or subsidence), as well as horizontal diffusion, can play an important role. Therefore, spatial information is required.
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Temperature change $\partial T_{SBL}$ at 8 m observed at the micrometeorological mast at Sodankylä (red line) and modelled from the lowermost model level (11 m) of the closest model grid point in REF (black solid line) during stably stratified, nocturnal periods. The simulated $\partial T_{SBL}$ is decomposed into the contributions from the turbulence scheme (orange line), and radiation scheme (blue line), the model physics (black dashed line), and the model dynamics (black dotted line)</figcaption></figure>

this coupling reduces and the turbulent heat loss to the surface dominates again, leading to cooling by turbulence.

The second night is different. Here, turbulence cools consistently as the lowermost model level decouples from the rest of the PBL, hindering the turbulent exchange of heat with higher model levels, and a vSBL forms. The continued heat loss to the surface causes a strong vertical flux divergence leading to a pronounced cooling of 9 K by turbulence throughout the night (Fig. 4, orange line). As only a relatively shallow layer is cooled, temperatures drop stronger compared to the wSBL in the previous night.

The two nights also exhibit changes in the contribution of the radiation tendency. In the first night, the continuous radiative cooling amounts to about 2 K. During the second night, the radiative cooling stops around 2300 UTC, leading to a reduced cooling of 1 K throughout the night. This reduction is due to the formation of a strong inversion atop the lowest model level in the vSBL and the associated increase in downward longwave radiation (see Sect. 3.3, Fig. 10d).

Our diagnostic yields a detailed insight into the changed properties of the local, physical processes at the investigated grid point. However, no clear statement can be drawn from the differences in the dynamical tendency (Fig. 4, black dotted line), as both small- and large-scale advection (or subsidence), as well as horizontal diffusion, can play an important role. Therefore, spatial information is required.
### 3.1.1 Spatial Pattern

We now broaden our analysis to the extracted sub-domain around Sodankylä (Sect. 2.2) and analyse the mean contribution of the individual processes. To avoid diluting the different signatures of turbulence and radiation within the wSBL and vSBL, respectively, we first examine the spatial occurrence of these two SBL regimes (Fig. 5). We distinguish between the wSBL and vSBL, using a threshold of 1.5 K for the mean temperature difference between the two lowermost model levels during each respective night. Even though the depicted mean values and standard deviations in Fig. 6 are sensitive to this threshold, partly due to the limited amount of grid boxes considered here, the displayed characteristics of each curve are robust. The final value of 1.5 K approximately minimizes the standard deviation.

The spatial occurrence of the two regimes is very similar in both nights, allowing us to depict it by using a single timestep here (0335 UTC 16 March 2018, Fig. 5). All shown fields exhibit a small-scale coexistence of the wSBL and vSBL. The wSBL is characterized by comparably warmer surface and 11 m-temperatures (Fig. 5a, b), comparably higher wind speeds ($> 1.5 \text{ m s}^{-1}$, Fig. 5c), and, most notably, little to no temperature inversions between the two lowermost model levels (Fig. 5d). The vSBL on the other hand is characterized by colder surface and 11 m-temperatures, low wind speeds, and pronounced temperature inversions of up to 4 K between the two lowermost model levels. It is noteworthy that the grid point closest to the micrometeorological mast is located at the border between these two SBL regimes. Therefore, a shift by one grid point would lead to pronounced differences in modelled temperatures (see Sect. 3.3).

A comparison with Fig. 2 shows that the occurrence of the SBL regimes does not follow any evident pattern in the snow cover or model physiography. However, some resemblance can be seen with the model topography, even though the area is relatively flat. For example, the “warm” region covered by the wSBL at the bottom left corner (Fig. 5a, b) is located at the highest elevation of the domain (Fig. 2a) and also displays the strongest wind speeds (Fig. 5c). Such wind speeds favour turbulent mixing and the formation and upkeep of the wSBL. The region covered by the vSBL, on the other hand, partly aligns with the surface depression found in the region.

We again depict the temperature change due to individual processes (see Fig. 4), but instead of a grid point perspective, we now broaden our analysis to the two stability regimes (Fig. 6). Thereby, we gain insights into how individual processes contribute within the wSBL and vSBL, respectively.

In the wSBL (Fig. 6a), the turbulence tendency strongly cools the lowermost model level by up to 4 K for the first 3–4 h in both nights (orange line). After 1800 UTC in both nights, the cooling rate due to turbulence reduces and even transitions into a warming in the first night at around 0200 UTC. Here, the downward, turbulent heat transport from the higher model levels in the wSBL (over)compensates the turbulent heat loss to the surface. Radiation, on the other hand, cools consistently, and amounts to a total cooling of about 2 K per night (blue dashed line). The model dynamics also contribute with 1–3 K to the nocturnal cooling each night (black dotted line).

In the vSBL (Fig. 6b), turbulence contributes continuously to the cooling of the lowest model level (8–13 K per night, orange line), as the turbulent heat loss to the surface cannot be efficiently compensated by any heat transport from upper model levels. The cooling by radiation also changes compared to the wSBL. Throughout both nights, the radiative cooling is not continuous anymore but reduces during the first night and even transitions to
### 3.1.1 Spatial Pattern

We now broaden our analysis to the extracted sub-domain around Sodankylä (Sect. 2.2) and analyse the mean contribution of the individual processes. To avoid diluting the different signatures of turbulence and radiation within the wSBL and vSBL, respectively, we first examine the spatial occurrence of these two SBL regimes (Fig. 5). We distinguish between the wSBL and vSBL, using a threshold of 1.5 K for the mean temperature difference between the two lowermost model levels during each respective night. Even though the depicted mean values and standard deviations in Fig. 6 are sensitive to this threshold, partly due to the limited amount of grid boxes considered here, the displayed characteristics of each curve are robust. The final value of 1.5 K approximately minimizes the standard deviation.

The spatial occurrence of the two regimes is very similar in both nights, allowing us to depict it by using a single timestep here (0335 UTC 16 March 2018, Fig. 5). All shown fields exhibit a small-scale coexistence of the wSBL and vSBL. The wSBL is characterized by comparably warmer surface and 11 m-temperatures (Fig. 5a, b), comparably higher wind speeds (>$1.5 \text{ m s}^{-1}$, Fig. 5c), and, most notably, little to no temperature inversions between the two lowermost model levels (Fig. 5d). The vSBL on the other hand is characterized by colder surface and 11 m-temperatures, low wind speeds, and pronounced temperature inversions of up to 4 K between the two lowermost model levels. It is noteworthy that the grid point closest to the micrometeorological mast is located at the border between these two SBL regimes. Therefore, a shift by one grid point would lead to pronounced differences in modelled temperatures (see Sect. 3.3).

A comparison with Fig. 2 shows that the occurrence of the SBL regimes does not follow any evident pattern in the snow cover or model physiography. However, some resemblance can be seen with the model topography, even though the area is relatively flat. For example, the “warm” region covered by the wSBL at the bottom left corner (Fig. 5a, b) is located at the highest elevation of the domain (Fig. 2a) and also displays the strongest wind speeds (Fig. 5c). Such wind speeds favour turbulent mixing and the formation and upkeep of the wSBL. The region covered by the vSBL, on the other hand, partly aligns with the surface depression found in the region.

We again depict the temperature change due to individual processes (see Fig. 4), but instead of a grid point perspective, we now broaden our analysis to the two stability regimes (Fig. 6). Thereby, we gain insights into how individual processes contribute within the wSBL and vSBL, respectively.

In the wSBL (Fig. 6a), the turbulence tendency strongly cools the lowermost model level by up to 4 K for the first 3–4 h in both nights (orange line). After 1800 UTC in both nights, the cooling rate due to turbulence reduces and even transitions into a warming in the first night at around 0200 UTC. Here, the downward, turbulent heat transport from the higher model levels in the wSBL (over)compensates the turbulent heat loss to the surface. Radiation, on the other hand, cools consistently, and amounts to a total cooling of about 2 K per night (blue dashed line). The model dynamics also contribute with 1–3 K to the nocturnal cooling each night (black dotted line).

In the vSBL (Fig. 6b), turbulence contributes continuously to the cooling of the lowest model level (8–13 K per night, orange line), as the turbulent heat loss to the surface cannot be efficiently compensated by any heat transport from upper model levels. The cooling by radiation also changes compared to the wSBL. Throughout both nights, the radiative cooling is not continuous anymore but reduces during the first night and even transitions to
<figure><img src="image_7.png" /><figcaption>Fig. 5 Domain of 14 × 14 grid boxes around the Sodankylä supersite extracted from REF at 0335 UTC 16 March 2018 depicting a effective surface temperature $T_{0m}$, b temperature at the first model level $T_{11m}$, c wind speed at the first model level $U_{11m}$, and d temperature difference between the first two model levels $T_{32m} - T_{11m}$. The black contour lines depict the differentiation between wSBL and vSBL, based on the threshold of 1.5 K for the mean, nocturnal temperature difference between the first two model levels. The location of the micrometeorological mast is marked by the red cross. The two black dots mark the grid points used in Sect. 3.3</figcaption></figure>

slight warming in the second night (see also Sect. 3.3). The model dynamics contribute with pronounced warming of up to 5 K per night, which counteracts the cooling by turbulence.

In summary, AROME-Arctic introduces both the wSBL and vSBL during our case study on small spatial scales of 5 km (two grid boxes) which has strong implications for validation studies of the SBL against point observations (see Sect. 3.3). Both SBL regimes exhibit different contributions by the model dynamics and physics. The turbulence tendency illustrates this difference best. In the wSBL, the atmospheric levels are well coupled leading to an efficient heat transport by turbulence. In the vSBL, on the other hand, such a coupling is absent and the lowermost model level strongly cools due to the turbulent heat loss to the surface. We now investigate how the distribution and representation of the wSBL and vSBL change when simulating the same case with the updated surface parametrization NEW (see Sect. 2.3.1).
<figure><img src="image_7.png" /><figcaption>Fig. 5 Domain of 14 × 14 grid boxes around the Sodankylä supersite extracted from REF at 0335 UTC 16 March 2018 depicting a effective surface temperature $T_{0m}$, b temperature at the first model level $T_{11m}$, c wind speed at the first model level $U_{11m}$, and d temperature difference between the first two model levels $T_{32m} - T_{11m}$. The black contour lines depict the differentiation between wSBL and vSBL, based on the threshold of 1.5 K for the mean, nocturnal temperature difference between the first two model levels. The location of the micrometeorological mast is marked by the red cross. The two black dots mark the grid points used in Sect. 3.3</figcaption></figure>

slight warming in the second night (see also Sect. 3.3). The model dynamics contribute with pronounced warming of up to 5 K per night, which counteracts the cooling by turbulence.

In summary, AROME-Arctic introduces both the wSBL and vSBL during our case study on small spatial scales of 5 km (two grid boxes) which has strong implications for validation studies of the SBL against point observations (see Sect. 3.3). Both SBL regimes exhibit different contributions by the model dynamics and physics. The turbulence tendency illustrates this difference best. In the wSBL, the atmospheric levels are well coupled leading to an efficient heat transport by turbulence. In the vSBL, on the other hand, such a coupling is absent and the lowermost model level strongly cools due to the turbulent heat loss to the surface. We now investigate how the distribution and representation of the wSBL and vSBL change when simulating the same case with the updated surface parametrization NEW (see Sect. 2.3.1).
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Mean temperature change $\overline{\partial T_{SBL}}$ in AROME-Arctic for the two consecutive nights in REF, distinguished between all grid points allocated with **a** the wSBL and **b** the vSBL. Line colouring as in Fig. 4. Error bars at the end of each respective line indicate the mean standard deviation throughout each night. *N* indicates the number of grid points allocated to the respective SBL regime per night</figcaption></figure>

## 3.2 The Nocturnal Boundary Layer in the Sensitivity Experiment

In the entire region, we notice reduced wind speeds in NEW (Fig. 7c) compared to REF (Fig. 5c), as a result of the updated roughness lengths due to the explicit treatment of the canopy layer. Interestingly, the wSBL and vSBL occur approximately at the same locations in REF and NEW (Figs. 5 and 7), even though the vSBL now covers a larger area of the
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Mean temperature change $\overline{\partial T_{SBL}}$ in AROME-Arctic for the two consecutive nights in REF, distinguished between all grid points allocated with **a** the wSBL and **b** the vSBL. Line colouring as in Fig. 4. Error bars at the end of each respective line indicate the mean standard deviation throughout each night. *N* indicates the number of grid points allocated to the respective SBL regime per night</figcaption></figure>

## 3.2 The Nocturnal Boundary Layer in the Sensitivity Experiment

In the entire region, we notice reduced wind speeds in NEW (Fig. 7c) compared to REF (Fig. 5c), as a result of the updated roughness lengths due to the explicit treatment of the canopy layer. Interestingly, the wSBL and vSBL occur approximately at the same locations in REF and NEW (Figs. 5 and 7), even though the vSBL now covers a larger area of the
<figure><img src="image_9.png" /><figcaption>Fig. 7 As in Fig. 5 but for experiment NEW</figcaption></figure>

extracted sub-domain. Again we find an agreement between the occurrence of the stability regimes, the ambient wind speed, and partially the model topography (Fig. 2a). Identifying the underlying reason for these nearly unchanged spatial patterns is outside the scope of this study, but we provide some recommendations in Sect. 4.4.

A comparison of the depicted temperatures in the wSBL between REF (Fig. 5) and NEW (Fig. 7) reveals only minor differences between the sensitivity experiments, even for surface temperatures. For the vSBL, on the other hand, the same comparison yields substantially colder grid boxes in NEW (Fig. 7a, b) and displays a stronger inversion between the two lowermost model levels (Fig. 7d). Initially, we expected overall colder surface and near-surface temperatures in NEW due to the multi-layer snow parametrization 3-L (Sect. 2.3.1). This parametrization captures the strong cooling of the upper snow layers due to radiative forcing better (Arduini et al. 2019), compared to the single-layer scheme employed in REF. However, only in the vSBL, this impact is clearly visible. To understand why NEW impacts the vSBL more than the wSBL, we again decompose the temperature change in these two SBL regimes with the high-resolution tendency output (Fig. 8).

In the wSBL, the mean nocturnal cooling amounts to about 7–8 K per night (Fig. 8a, black solid line), similar to citation issue “REF”. However, compared to REF, the underlying processes are distinctly different (see Fig. 6a). In NEW, the turbulence parametrization warms the lowermost model level by up to 3 K per night (orange line). This heating reflects an imbalance between the downward turbulent transport of heat and the turbulent heat loss to
<figure><img src="image_9.png" /><figcaption>Fig. 7 As in Fig. 5 but for experiment NEW</figcaption></figure>

extracted sub-domain. Again we find an agreement between the occurrence of the stability regimes, the ambient wind speed, and partially the model topography (Fig. 2a). Identifying the underlying reason for these nearly unchanged spatial patterns is outside the scope of this study, but we provide some recommendations in Sect. 4.4.

A comparison of the depicted temperatures in the wSBL between REF (Fig. 5) and NEW (Fig. 7) reveals only minor differences between the sensitivity experiments, even for surface temperatures. For the vSBL, on the other hand, the same comparison yields substantially colder grid boxes in NEW (Fig. 7a, b) and displays a stronger inversion between the two lowermost model levels (Fig. 7d). Initially, we expected overall colder surface and near-surface temperatures in NEW due to the multi-layer snow parametrization 3-L (Sect. 2.3.1). This parametrization captures the strong cooling of the upper snow layers due to radiative forcing better (Arduini et al. 2019), compared to the single-layer scheme employed in REF. However, only in the vSBL, this impact is clearly visible. To understand why NEW impacts the vSBL more than the wSBL, we again decompose the temperature change in these two SBL regimes with the high-resolution tendency output (Fig. 8).

In the wSBL, the mean nocturnal cooling amounts to about 7–8 K per night (Fig. 8a, black solid line), similar to citation issue “REF”. However, compared to REF, the underlying processes are distinctly different (see Fig. 6a). In NEW, the turbulence parametrization warms the lowermost model level by up to 3 K per night (orange line). This heating reflects an imbalance between the downward turbulent transport of heat and the turbulent heat loss to
<figure><img src="image_10.png" /><figcaption>Fig. 8 As in Fig. 6 but for experiment NEW</figcaption></figure>

the surface. Meanwhile, radiation strongly cools the lowermost model level by up to 7K per night (blue line), almost quadrupling its contribution in the wSBL compared to REF, despite similar surface temperatures (Figs. 5a and 7a). This increase in radiative cooling is related to the enhanced snow cover fraction and updated calculation of surface emissivity by the ISBA-MEB parametrization (Sect. 2.3.1), yielding a much lower radiative surface temperature. However, the enhanced cooling by radiation is compensated by the downward, turbulent transport of heat at 0000 UTC during the first night and 2100 UTC during the second night, leading to only little temperature changes by the model physics (Fig. 8a, black
<figure><img src="image_10.png" /><figcaption>Fig. 8 As in Fig. 6 but for experiment NEW</figcaption></figure>

the surface. Meanwhile, radiation strongly cools the lowermost model level by up to 7K per night (blue line), almost quadrupling its contribution in the wSBL compared to REF, despite similar surface temperatures (Figs. 5a and 7a). This increase in radiative cooling is related to the enhanced snow cover fraction and updated calculation of surface emissivity by the ISBA-MEB parametrization (Sect. 2.3.1), yielding a much lower radiative surface temperature. However, the enhanced cooling by radiation is compensated by the downward, turbulent transport of heat at 0000 UTC during the first night and 2100 UTC during the second night, leading to only little temperature changes by the model physics (Fig. 8a, black
dashed line). The model dynamics nearly continuously cool the lowest model level in the wSBL (Fig. 8a, black dotted line) with 3–4 K per night.

In the vSBL, the total amount of mean cooling reaches values of 10–12 K and is solely driven by the model physics (Fig. 8b). Both the radiation and turbulence parametrizations cool the lowermost model level. In the first night, the mean cooling by the turbulence parametrization is not continuous, but transitions to a warming between 1800 to 2100 UTC, correlating with the occurrence of an ice cloud. This cloud led to snowfall over parts of the extracted domain (not shown) and to warming of the lowermost model levels by the turbulence parametrization. On the second night, the cooling by turbulence dominates again, reaching mean values of about 12 K. Furthermore, the mean cooling by radiation does not come to a halt throughout both nights, in contrast to REF. The model dynamics warm the lowermost model level in the vSBL, again counteracting the cooling, as in REF.

In summary, the updated surface treatment in NEW strongly increases the radiative cooling of the near-surface atmosphere in AROME-Arctic during the night. Depending on the stability regime, this increase stems from the updated surface emissivity (wSBL) or from a combination of the updated surface emissivity and colder surface temperatures (vSBL). In the wSBL, the surface temperatures between REF and NEW are nearly unchanged. However, the updated surface emissivity in NEW reduces the emitted longwave radiation from the surface, causing a more negative net upward radiative flux at the lowermost model level (and above). Since the atmospheric temperature change due to the radiation parametrization depends on the vertical radiative flux divergence at each level (Eq. 4), this decrease in emitted radiation from the surface leads to pronounced radiative cooling in the near-surface atmosphere. In the wSBL, this enhanced cooling by radiation is compensated by an increased downward turbulent transport of heat, whereas in the vSBL, it contributes to the stronger cooling of near-surface temperatures.

### 3.3 Vertical Depiction of the Two Stability Regimes

So far our analysis focussed on the lowermost model level. We now extend this analysis and investigate the contributions of individual processes in the entire boundary layer, again separated by stability regime. Therefore, we will use two grid points that represent the wSBL and vSBL, respectively.

In our case study, the micrometeorological mast is located at the border between the wSBL and vSBL. We choose two grid boxes closest to the mast, which are distinctly located in either the wSBL or vSBL throughout the second night in both experiments (Figs. 5 and 7, black dots). These two grid boxes are approximately 7 km apart from each other, emphasizing the small-scale heterogeneity within the representation of the nocturnal SBL in AROME-Arctic. We chose the second night here as the first night in NEW featured a brief snowfall event, not present in REF, which had an impact on the modelled temperatures (Sect. 3.2).

The modelled temperature profiles throughout the second night (Fig. 9) exhibit concave (wSBL) and convex (vSBL) shaped forms, respectively, forming near-neutral or strongly stratified near-surface atmospheres (Ulden and Holtslag 1985; Vogelezang and Holtslag 1996; Vignon et al. 2017). The temperature profiles of the wSBL (Fig. 9a) do not distinctly differ between REF and NEW, agreeing with our findings in Sect. 3.2. The temperature profiles of the vSBL (Fig. 9b), on the other hand, show colder near-surface temperatures in NEW throughout the night, also in agreement with our previous analysis. The observed temperature profiles at the mast (Fig. 9c) coincide well with the temperature development in the vSBL and would yield large temperature biases if being compared to profiles of the wSBL (compare
dashed line). The model dynamics nearly continuously cool the lowest model level in the wSBL (Fig. 8a, black dotted line) with 3–4 K per night.

In the vSBL, the total amount of mean cooling reaches values of 10–12 K and is solely driven by the model physics (Fig. 8b). Both the radiation and turbulence parametrizations cool the lowermost model level. In the first night, the mean cooling by the turbulence parametrization is not continuous, but transitions to a warming between 1800 to 2100 UTC, correlating with the occurrence of an ice cloud. This cloud led to snowfall over parts of the extracted domain (not shown) and to warming of the lowermost model levels by the turbulence parametrization. On the second night, the cooling by turbulence dominates again, reaching mean values of about 12 K. Furthermore, the mean cooling by radiation does not come to a halt throughout both nights, in contrast to REF. The model dynamics warm the lowermost model level in the vSBL, again counteracting the cooling, as in REF.

In summary, the updated surface treatment in NEW strongly increases the radiative cooling of the near-surface atmosphere in AROME-Arctic during the night. Depending on the stability regime, this increase stems from the updated surface emissivity (wSBL) or from a combination of the updated surface emissivity and colder surface temperatures (vSBL). In the wSBL, the surface temperatures between REF and NEW are nearly unchanged. However, the updated surface emissivity in NEW reduces the emitted longwave radiation from the surface, causing a more negative net upward radiative flux at the lowermost model level (and above). Since the atmospheric temperature change due to the radiation parametrization depends on the vertical radiative flux divergence at each level (Eq. 4), this decrease in emitted radiation from the surface leads to pronounced radiative cooling in the near-surface atmosphere. In the wSBL, this enhanced cooling by radiation is compensated by an increased downward turbulent transport of heat, whereas in the vSBL, it contributes to the stronger cooling of near-surface temperatures.

### 3.3 Vertical Depiction of the Two Stability Regimes

So far our analysis focussed on the lowermost model level. We now extend this analysis and investigate the contributions of individual processes in the entire boundary layer, again separated by stability regime. Therefore, we will use two grid points that represent the wSBL and vSBL, respectively.

In our case study, the micrometeorological mast is located at the border between the wSBL and vSBL. We choose two grid boxes closest to the mast, which are distinctly located in either the wSBL or vSBL throughout the second night in both experiments (Figs. 5 and 7, black dots). These two grid boxes are approximately 7 km apart from each other, emphasizing the small-scale heterogeneity within the representation of the nocturnal SBL in AROME-Arctic. We chose the second night here as the first night in NEW featured a brief snowfall event, not present in REF, which had an impact on the modelled temperatures (Sect. 3.2).

The modelled temperature profiles throughout the second night (Fig. 9) exhibit concave (wSBL) and convex (vSBL) shaped forms, respectively, forming near-neutral or strongly stratified near-surface atmospheres (Ulden and Holtslag 1985; Vogelezang and Holtslag 1996; Vignon et al. 2017). The temperature profiles of the wSBL (Fig. 9a) do not distinctly differ between REF and NEW, agreeing with our findings in Sect. 3.2. The temperature profiles of the vSBL (Fig. 9b), on the other hand, show colder near-surface temperatures in NEW throughout the night, also in agreement with our previous analysis. The observed temperature profiles at the mast (Fig. 9c) coincide well with the temperature development in the vSBL and would yield large temperature biases if being compared to profiles of the wSBL (compare
<figure><img src="image_11.png" /><figcaption>**Fig. 9** Air temperature profiles throughout the second night of **a** a grid box exhibiting the wSBL, **b** a grid box exhibiting the vSBL, and **c** the profiles observed at the 48 m-micrometeorological mast. The two grid boxes depicted in (**a**) and (**b**) are about 3 km away from the mast and about 7 km apart from one another. The markers indicate the time stamps of the profiles in UTC</figcaption></figure>

Fig. 4). Furthermore, the observed profiles agree better with modelled profiles in NEW than in REF.

We now investigate the individual processes at these two grid points in both REF and NEW. Starting with REF, radiation cools most of the SBL in both the wSBL and vSBL (Fig. 10a, d). The cooling is more pronounced in the vSBL throughout most of the PBL depth, evidencing a larger radiative imbalance. The most notable difference between the two regimes appears at the lowermost level. In the wSBL, a continuous cooling by radiation is evident, while in the vSBL, a transition from cooling to warming occurs. This warming originates from the strong inversion that forms close to the surface in the vSBL (Fig. 9a). The warmer temperatures at the top of the inversion increase the downward longwave radiation and lead to a positive radiative imbalance at the lowermost model level. Thereby, the radiation parametrization partly counteracts the strong cooling of near-surface temperatures in the vSBL in REF.

<figure><img src="image_12.png" /><figcaption>**Fig. 10** Contributions of the radiation parametrization (**a**, **d**), the turbulence parametrization (**b**, **e**), and the model dynamics (**c**, **f**) throughout the second night in the lowest 220 m in for the reference run REF for the two grid points representing the wSBL (top row) and vSBL (bottom row), respectively</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 9** Air temperature profiles throughout the second night of **a** a grid box exhibiting the wSBL, **b** a grid box exhibiting the vSBL, and **c** the profiles observed at the 48 m-micrometeorological mast. The two grid boxes depicted in (**a**) and (**b**) are about 3 km away from the mast and about 7 km apart from one another. The markers indicate the time stamps of the profiles in UTC</figcaption></figure>

Fig. 4). Furthermore, the observed profiles agree better with modelled profiles in NEW than in REF.

We now investigate the individual processes at these two grid points in both REF and NEW. Starting with REF, radiation cools most of the SBL in both the wSBL and vSBL (Fig. 10a, d). The cooling is more pronounced in the vSBL throughout most of the PBL depth, evidencing a larger radiative imbalance. The most notable difference between the two regimes appears at the lowermost level. In the wSBL, a continuous cooling by radiation is evident, while in the vSBL, a transition from cooling to warming occurs. This warming originates from the strong inversion that forms close to the surface in the vSBL (Fig. 9a). The warmer temperatures at the top of the inversion increase the downward longwave radiation and lead to a positive radiative imbalance at the lowermost model level. Thereby, the radiation parametrization partly counteracts the strong cooling of near-surface temperatures in the vSBL in REF.

<figure><img src="image_12.png" /><figcaption>**Fig. 10** Contributions of the radiation parametrization (**a**, **d**), the turbulence parametrization (**b**, **e**), and the model dynamics (**c**, **f**) throughout the second night in the lowest 220 m in for the reference run REF for the two grid points representing the wSBL (top row) and vSBL (bottom row), respectively</figcaption></figure>
The differences between the wSBL and vSBL can best be seen by the turbulence tendency (Fig. 10b, e). At the beginning of the night, the entire PBL is cooled by turbulence due to the surface cooling. However, from 1800 UTC onward, differences emerge. The wSBL demonstrates an efficient heat transport within the SBL by cooling the higher model levels (150–50 m) by up to 3 $Kh^{-1}$ and warming the near-surface levels (< 50 m) with 1–2 $Kh^{-1}$. The reduced magnitude of the warming at the lowermost model levels can be explained by the turbulent heat loss to the surface. In the vSBL, the lowest model levels cool by 1–2 $Kh^{-1}$ due to turbulent heat loss to the surface, and the SBL collapses within the first 3 h of the night. Thereby, the lowermost model level decouples from the other levels, forming a very shallow layer close to the ground that is continuously cooled by turbulence. The model levels between 175 and 25 m in the vSBL exhibit downward transport of heat from higher to lower levels, similar to a wSBL with values of ±0.5 $Kh^{-1}$. The model dynamics counteract the cooling or warming by the turbulent heat transport in both the wSBL and vSBL (Fig. 10c, f).

In NEW, the radiation parametrization exhibits a much stronger cooling compared to REF throughout the entire SBL, in both regimes (Fig. 11a, d). In the wSBL, this enhanced cooling is compensated by an increase in warming by the turbulence parametrization (Fig. 11b), yielding similar near-surface temperatures in REF and NEW (Fig. 9b). In the vSBL, radiation in NEW does not warm the lowermost model level despite the strong inversion in near-surface temperatures, but instead exhibits a decreasing cooling rate throughout the night. The absence of the warming by radiation partly explains the stronger decrease in near-surface temperatures in NEW within the vSBL compared to REF (compare Figs. 5a and 7a) in addition to the more pronounced surface cooling due to the multi-layer snow scheme in NEW.

<figure><img src="image_13.png" /><figcaption>**Fig. 11** As Fig. 10, but for sensitivity experiment NEW</figcaption></figure>

## 4 Discussion

The presented results demonstrate that operational NWP models such as AROME-Arctic can simulate both the wSBL and vSBL (agreeing with Baas et al. 2018). In the wSBL, the model levels are strongly coupled by turbulence, leading to efficient downward transport of sensible heat towards the surface, cooling upper levels, and warming lower ones. Pronounced negative
The differences between the wSBL and vSBL can best be seen by the turbulence tendency (Fig. 10b, e). At the beginning of the night, the entire PBL is cooled by turbulence due to the surface cooling. However, from 1800 UTC onward, differences emerge. The wSBL demonstrates an efficient heat transport within the SBL by cooling the higher model levels (150–50 m) by up to 3 K h⁻¹ and warming the near-surface levels (< 50 m) with 1–2 K h⁻¹. The reduced magnitude of the warming at the lowermost model levels can be explained by the turbulent heat loss to the surface. In the vSBL, the lowest model levels cool by 1–2 K h⁻¹ due to turbulent heat loss to the surface, and the SBL collapses within the first 3 h of the night. Thereby, the lowermost model level decouples from the other levels, forming a very shallow layer close to the ground that is continuously cooled by turbulence. The model levels between 175 and 25 m in the vSBL exhibit downward transport of heat from higher to lower levels, similar to a wSBL with values of ±0.5 K h⁻¹. The model dynamics counteract the cooling or warming by the turbulent heat transport in both the wSBL and vSBL (Fig. 10c, f).

In NEW, the radiation parametrization exhibits a much stronger cooling compared to REF throughout the entire SBL, in both regimes (Fig. 11a, d). In the wSBL, this enhanced cooling is compensated by an increase in warming by the turbulence parametrization (Fig. 11b), yielding similar near-surface temperatures in REF and NEW (Fig. 9b). In the vSBL, radiation in NEW does not warm the lowermost model level despite the strong inversion in near-surface temperatures, but instead exhibits a decreasing cooling rate throughout the night. The absence of the warming by radiation partly explains the stronger decrease in near-surface temperatures in NEW within the vSBL compared to REF (compare Figs. 5a and 7a) in addition to the more pronounced surface cooling due to the multi-layer snow scheme in NEW.

<figure><img src="image_13.png" /><figcaption>**Fig. 11** As Fig. 10, but for sensitivity experiment NEW</figcaption></figure>

# 4 Discussion

The presented results demonstrate that operational NWP models such as AROME-Arctic can simulate both the wSBL and vSBL (agreeing with Baas et al. 2018). In the wSBL, the model levels are strongly coupled by turbulence, leading to efficient downward transport of sensible heat towards the surface, cooling upper levels, and warming lower ones. Pronounced negative
surface sensible heat fluxes of − 40 W m⁻² occur in the model (not shown). However, due to the large heat capacity of the coupled SBL, the boundary layer cools comparably slow.

The vSBL exhibits strongly reduced turbulence, and the lowest model level becomes decoupled. Here, the surface sensible heat flux in the model becomes very small (about −5 W m⁻²). However, as only a shallow layer needs to be cooled, strongly negative turbulent heating rates occur and temperatures drop fast.

In both cases, the turbulence tendency is indicative of the respective SBL regime. In that regard, it should be noted that the halt in turbulent cooling used for characterizing the wSBL (Figs. 4 and 6a) could also stem from a decoupling between the surface and the lowermost model level, as that would cause a cessation of the turbulent heat loss towards the surface. However, AROME-Arctic prevents such a decoupling from happening. The maximum Richardson number at the surface is set to zero, implying neutral stratification, which always allows for an exchange between the surface and the atmosphere. The influence of this operational setting on the SBL in AROME-Arctic is the subject of ongoing research (personal comm., M. Homleid, E. Samuelsen, MET-Norway 2022) and outside the scope of this study.

Both SBL regimes coexist on small spatial scales of about 5 km. This scale is most likely set by the model as it equals two grid lengths and does not necessarily have a direct correspondence to nature. However, the pronounced grid-scale and subgrid-scale differences between the wSBL and vSBL raise questions on how such a small-scale heterogeneity might impact model verification or model development. These questions deserve a more detailed discussion, following below.

## 4.1 Impacts on Verification

The small scales, at which the wSBL and vSBL coexist, highlight the importance of local effects for the SBL (Atlaskin and Vihma 2012) and have consequences for verification studies. As evident in Fig. 9, verification based on a single grid point can yield large model biases that would not adequately reflect the model’s performance. The estimated error would thus be a representativeness error instead of a model error (Kanamitsu and DeHaan 2011). For a number of regional forecasting systems including AROME-Arctic, Køltzow et al. (2019) find that the differences between observed and forecasted variables, such as 2-m temperature or 10-m wind speed, are substantially impacted by representativeness errors. The overall lack of observations in polar regions (Casati et al. 2017) combined with the comparably high-spatial resolution of regional forecasting models can further aggravate this problem (Kanamitsu and DeHaan 2011).

In order to mitigate such problems for model verification in the SBL, we recommend that a single-point observation should not be compared to the closest model grid point. Instead, the stability regimes of surrounding grid points should be classified and compared to the conditions at the observation site. In this study, we used a threshold value for the inversion strength between the two lowermost atmospheric model levels. Other possibilities are provided by the scaled curvature parameter for the temperature profiles (André and Mahrt 1982) or a wind speed threshold successfully applied in several studies such as Vignon et al. (2017) and Baas et al. (2018). This way, a more tailored verification of the SBL can be achieved.
surface sensible heat fluxes of − 40 W m⁻² occur in the model (not shown). However, due to the large heat capacity of the coupled SBL, the boundary layer cools comparably slow.

The vSBL exhibits strongly reduced turbulence, and the lowest model level becomes decoupled. Here, the surface sensible heat flux in the model becomes very small (about −5 W m⁻²). However, as only a shallow layer needs to be cooled, strongly negative turbulent heating rates occur and temperatures drop fast.

In both cases, the turbulence tendency is indicative of the respective SBL regime. In that regard, it should be noted that the halt in turbulent cooling used for characterizing the wSBL (Figs. 4 and 6a) could also stem from a decoupling between the surface and the lowermost model level, as that would cause a cessation of the turbulent heat loss towards the surface. However, AROME-Arctic prevents such a decoupling from happening. The maximum Richardson number at the surface is set to zero, implying neutral stratification, which always allows for an exchange between the surface and the atmosphere. The influence of this operational setting on the SBL in AROME-Arctic is the subject of ongoing research (personal comm., M. Homleid, E. Samuelsen, MET-Norway 2022) and outside the scope of this study.

Both SBL regimes coexist on small spatial scales of about 5 km. This scale is most likely set by the model as it equals two grid lengths and does not necessarily have a direct correspondence to nature. However, the pronounced grid-scale and subgrid-scale differences between the wSBL and vSBL raise questions on how such a small-scale heterogeneity might impact model verification or model development. These questions deserve a more detailed discussion, following below.

## 4.1 Impacts on Verification

The small scales, at which the wSBL and vSBL coexist, highlight the importance of local effects for the SBL (Atlaskin and Vihma 2012) and have consequences for verification studies. As evident in Fig. 9, verification based on a single grid point can yield large model biases that would not adequately reflect the model’s performance. The estimated error would thus be a representativeness error instead of a model error (Kanamitsu and DeHaan 2011). For a number of regional forecasting systems including AROME-Arctic, Køltzow et al. (2019) find that the differences between observed and forecasted variables, such as 2-m temperature or 10-m wind speed, are substantially impacted by representativeness errors. The overall lack of observations in polar regions (Casati et al. 2017) combined with the comparably high-spatial resolution of regional forecasting models can further aggravate this problem (Kanamitsu and DeHaan 2011).

In order to mitigate such problems for model verification in the SBL, we recommend that a single-point observation should not be compared to the closest model grid point. Instead, the stability regimes of surrounding grid points should be classified and compared to the conditions at the observation site. In this study, we used a threshold value for the inversion strength between the two lowermost atmospheric model levels. Other possibilities are provided by the scaled curvature parameter for the temperature profiles (André and Mahrt 1982) or a wind speed threshold successfully applied in several studies such as Vignon et al. (2017) and Baas et al. (2018). This way, a more tailored verification of the SBL can be achieved.
## 4.2 Impacts on Model Performance and Development

We do not believe that the narrow misrepresentation of the vSBL and wSBL is responsible for the frequently documented warm bias of NWP models in polar regions during stably stratified conditions. It appears more likely that the SBL in NWP models is too well coupled, i.e. by too much internal turbulent mixing. Such mixing would result in a too deep boundary layer (Svensson and Holtslag 2009; Holtslag et al. 2013), which in turn yields too large heat capacities, resulting in too much thermal inertia. Such inertia prevents a fast reaction of atmospheric temperatures to surface forcing in the model, causing the formation of a warm bias (Beesley et al. 2000; Esau et al. 2018). The common practice of enhancing turbulent mixing in global NWP models such as the Integrated Forecasting System, albeit for well-founded reasons (Sandu et al. 2013), aggravates this problem.

Another important factor for the temperature development in the near-surface atmosphere is the representation of snow cover. A single-layer snow scheme often underestimates the magnitude and speed of temperature drops in the snow pack as it distributes the surface cooling over the entire snow depth. A multi-layer snow scheme on the other hand is capable of capturing the strong temperature gradients in the snow pack and yields stronger drops of near-surface temperatures (Haiden et al. 2018; Arduini et al. 2019; Day et al. 2020). However, the analysis of the sensitivity experiment NEW, with a multi-layer snow scheme, provides a more differentiated picture. Grid boxes covered by the wSBL exhibit nearly unchanged effective surface temperatures, as well as nearly unchanged temperatures, at the lowest model level. Despite the ability of the parametrized snow in NEW to react faster and more adequately to the radiative forcing, the continuous heat supply by turbulence from the boundary layer balances a more efficient cooling of the surface. In the vSBL, however, where this heat supply does not exist, the modifications to the surface are not counteracted and can have a pronounced impact on near-surface temperatures and stability.

Model optimization can be susceptible to such case-depending responses of the model, especially in regards to the common practice of model tuning. A possible consequence can be an artificial enhancement of a specific process or parameter by tuning, for example the heat conductivity of the snow pack, due to the absence of a desired effect, for example a more pronounced cooling of the surface. Such an artificial enhancement might deteriorate the model's performance for different atmospheric conditions and impede model development. Thus, such a varying and sometimes compensating interplay of subgrid-scale parameterization schemes, in this case between the snow scheme and the turbulence scheme, needs to be accounted for when testing and optimizing updates to the model configuration.

## 4.3 Role of Model Dynamics

Both SBL regimes exhibit a distinct temperature contribution from the model dynamics at the lowermost model level: cooling in the wSBL and warming in the vSBL. The model dynamics appear to balance the contributions by the turbulence scheme to a certain degree (Figs. 10 and 11), but not fully (Figs. 6 and 8). The dynamical tendency is comprised of advection, subsidence, and horizontal diffusion (Kähnert et al. 2021). All of these three processes can contribute to the compensating interplay observed here.

The horizontal advection of air will most likely be important at the “inflow” borders of the vSBL regime, when air is advected from the wSBL to the vSBL. In that case, warmer air flows over a colder surface in the vSBL, supplying heat to the turbulent heat flux and having a warming effect on the atmosphere. Likewise, the wSBL loses heat due to the flux divergence
## 4.2 Impacts on Model Performance and Development

We do not believe that the narrow misrepresentation of the vSBL and wSBL is responsible for the frequently documented warm bias of NWP models in polar regions during stably stratified conditions. It appears more likely that the SBL in NWP models is too well coupled, i.e. by too much internal turbulent mixing. Such mixing would result in a too deep boundary layer (Svensson and Holtslag 2009; Holtslag et al. 2013), which in turn yields too large heat capacities, resulting in too much thermal inertia. Such inertia prevents a fast reaction of atmospheric temperatures to surface forcing in the model, causing the formation of a warm bias (Beesley et al. 2000; Esau et al. 2018). The common practice of enhancing turbulent mixing in global NWP models such as the Integrated Forecasting System, albeit for well-founded reasons (Sandu et al. 2013), aggravates this problem.

Another important factor for the temperature development in the near-surface atmosphere is the representation of snow cover. A single-layer snow scheme often underestimates the magnitude and speed of temperature drops in the snow pack as it distributes the surface cooling over the entire snow depth. A multi-layer snow scheme on the other hand is capable of capturing the strong temperature gradients in the snow pack and yields stronger drops of near-surface temperatures (Haiden et al. 2018; Arduini et al. 2019; Day et al. 2020). However, the analysis of the sensitivity experiment NEW, with a multi-layer snow scheme, provides a more differentiated picture. Grid boxes covered by the wSBL exhibit nearly unchanged effective surface temperatures, as well as nearly unchanged temperatures, at the lowest model level. Despite the ability of the parametrized snow in NEW to react faster and more adequately to the radiative forcing, the continuous heat supply by turbulence from the boundary layer balances a more efficient cooling of the surface. In the vSBL, however, where this heat supply does not exist, the modifications to the surface are not counteracted and can have a pronounced impact on near-surface temperatures and stability.

Model optimization can be susceptible to such case-depending responses of the model, especially in regards to the common practice of model tuning. A possible consequence can be an artificial enhancement of a specific process or parameter by tuning, for example the heat conductivity of the snow pack, due to the absence of a desired effect, for example a more pronounced cooling of the surface. Such an artificial enhancement might deteriorate the model's performance for different atmospheric conditions and impede model development. Thus, such a varying and sometimes compensating interplay of subgrid-scale parameterization schemes, in this case between the snow scheme and the turbulence scheme, needs to be accounted for when testing and optimizing updates to the model configuration.

## 4.3 Role of Model Dynamics

Both SBL regimes exhibit a distinct temperature contribution from the model dynamics at the lowermost model level: cooling in the wSBL and warming in the vSBL. The model dynamics appear to balance the contributions by the turbulence scheme to a certain degree (Figs. 10 and 11), but not fully (Figs. 6 and 8). The dynamical tendency is comprised of advection, subsidence, and horizontal diffusion (Kähnert et al. 2021). All of these three processes can contribute to the compensating interplay observed here.

The horizontal advection of air will most likely be important at the “inflow” borders of the vSBL regime, when air is advected from the wSBL to the vSBL. In that case, warmer air flows over a colder surface in the vSBL, supplying heat to the turbulent heat flux and having a warming effect on the atmosphere. Likewise, the wSBL loses heat due to the flux divergence
and cools. However, this process is less relevant at the interior and “outflow” borders of the vSBL domain. Here, and elsewhere, the horizontal diffusion might be responsible for the observed behaviour. AROME-Arctic employs a semi-Lagrangian horizontal diffusion scheme (SLHD, Váňa et al. 2008) which locally modulates the interpolation done by the semi-Lagrangian advection scheme. Thereby, the signals of neighbouring grid cells are mixed which can lead to a cooling of warmer grid cells and vice versa. In the case of the vSBL, a warming of the lowermost model level by subsidence is also plausible due to the presence of the pronounced near-surface inversion.

To investigate this interplay and the role of individual processes, one could modify the strength of the SLHD as in Bengtsson et al. (2012) and evaluate the impacts with the individual tendency output.

## 4.4 Role of the Bifurcation Between the Stability Regimes

The bifurcation between the wSBL and vSBL poses a potential limit of predictability of the SBL for NWP models (McNider et al. 1995). The analysis presented in this study is not suited to identify the underlying causes that lead to the development of either the wSBL or vSBL in AROME-Arctic. While wind speed and to some extend topography plays an important role (agreeing with Derbyshire 1999; de Wiel et al. 2017; Vignon et al. 2017; Baas et al. 2018), we cannot directly link physiography fields such as snow cover and vegetation to the occurrence of the stability regimes.

To investigate the bifurcation, one can make use of the framework of demand and supply, introduced by de Wiel et al. (2017). Hereby, the demand describes the (radiative) heat loss of a well-mixed SBL, which needs to be met by the supply. The supply is represented by turbulent heat fluxes and ground heat fluxes. If the demand is met, a wSBL is maintained, whereas if the demand is not met, a vSBL forms.

For a certain critical wind speed window, two regimes exist (with weak and with strong surface inversion) that are in balance with the demand (de Wiel et al. 2017). It is within this critical wind speed window that the bifurcation occurs. Recasting the framework of demand and supply by the means of individual tendencies will allow to investigate the bifurcation in an operational NWP model such as AROME-Arctic.

## 5 Conclusion

The present work utilizes high-resolution output of individual physical tendencies from the operational NWP model AROME-Arctic to investigate the temperature development of the stable boundary layer during cloud-free, snow-covered nights. Hereby, the analysis distinguishes between the weakly stable and very stable boundary layer regimes that can develop during these conditions. The employed combination of DDH and tendency output enables insight into the contributions from parametrized and resolved processes at every model timestep, helping us to understand the evolution of temperature in the SBL in AROME-Arctic to great detail.

AROME-Arctic is capable of simulating both wSBL and vSBL. In the wSBL, turbulence responds to the radiative cooling of the surface by efficiently mixing heat downwards, warming the lower layers and cooling the upper layers of the SBL. In the vSBL, the radiative cooling leads to an increase in atmospheric stability that reduces turbulent mixing and leads to a decoupling and rapid cooling of the lowermost model level. The subsequent develop-
and cools. However, this process is less relevant at the interior and “outflow” borders of the vSBL domain. Here, and elsewhere, the horizontal diffusion might be responsible for the observed behaviour. AROME-Arctic employs a semi-Lagrangian horizontal diffusion scheme (SLHD, Váňa et al. 2008) which locally modulates the interpolation done by the semi-Lagrangian advection scheme. Thereby, the signals of neighbouring grid cells are mixed which can lead to a cooling of warmer grid cells and vice versa. In the case of the vSBL, a warming of the lowermost model level by subsidence is also plausible due to the presence of the pronounced near-surface inversion.

To investigate this interplay and the role of individual processes, one could modify the strength of the SLHD as in Bengtsson et al. (2012) and evaluate the impacts with the individual tendency output.

## 4.4 Role of the Bifurcation Between the Stability Regimes

The bifurcation between the wSBL and vSBL poses a potential limit of predictability of the SBL for NWP models (McNider et al. 1995). The analysis presented in this study is not suited to identify the underlying causes that lead to the development of either the wSBL or vSBL in AROME-Arctic. While wind speed and to some extend topography plays an important role (agreeing with Derbyshire 1999; de Wiel et al. 2017; Vignon et al. 2017; Baas et al. 2018), we cannot directly link physiography fields such as snow cover and vegetation to the occurrence of the stability regimes.

To investigate the bifurcation, one can make use of the framework of demand and supply, introduced by de Wiel et al. (2017). Hereby, the demand describes the (radiative) heat loss of a well-mixed SBL, which needs to be met by the supply. The supply is represented by turbulent heat fluxes and ground heat fluxes. If the demand is met, a wSBL is maintained, whereas if the demand is not met, a vSBL forms.

For a certain critical wind speed window, two regimes exist (with weak and with strong surface inversion) that are in balance with the demand (de Wiel et al. 2017). It is within this critical wind speed window that the bifurcation occurs. Recasting the framework of demand and supply by the means of individual tendencies will allow to investigate the bifurcation in an operational NWP model such as AROME-Arctic.

## 5 Conclusion

The present work utilizes high-resolution output of individual physical tendencies from the operational NWP model AROME-Arctic to investigate the temperature development of the stable boundary layer during cloud-free, snow-covered nights. Hereby, the analysis distinguishes between the weakly stable and very stable boundary layer regimes that can develop during these conditions. The employed combination of DDH and tendency output enables insight into the contributions from parametrized and resolved processes at every model timestep, helping us to understand the evolution of temperature in the SBL in AROME-Arctic to great detail.

AROME-Arctic is capable of simulating both wSBL and vSBL. In the wSBL, turbulence responds to the radiative cooling of the surface by efficiently mixing heat downwards, warming the lower layers and cooling the upper layers of the SBL. In the vSBL, the radiative cooling leads to an increase in atmospheric stability that reduces turbulent mixing and leads to a decoupling and rapid cooling of the lowermost model level. The subsequent develop-
ment of a strong near-surface inversion in the vSBL can lead to a warming by radiation at the lowermost model level. Follow-up studies should test the impact of an increased vertical resolution on the process representation.

Both stability regimes coexist on small spatial scales of about 5 km during our study, most likely set by the model grid, and yield substantially different atmospheric temperature profiles. Such small-scale heterogeneity in the modelled SBL needs to be considered for validation studies, i.e. by investigating a sufficiently large sub-domain around the observational site.

The updated surface scheme, NEW, leads to a substantial increase in radiative cooling of the lower atmosphere due to a reduced effective emissivity of the surface which again leads to a larger radiative flux divergence. Yet, in terms of the investigated surface and near-surface temperature fields, NEW only exhibits pronounced changes to the vSBL. Here, the multi-layer snow scheme better captures the strong cooling of the upper layers in the snow pack compared to the single-layer snow scheme in REF. In the wSBL, however, this more efficient cooling of the surface is counteracted by the continuous heat supply of the turbulence scheme, yielding nearly unchanged surface and near-surface temperatures in both experiments. Such a varying and potentially compensating interplay between parametrization schemes is important to consider during model development, especially in regard to the common practice of model tuning.

A natural next step now is to validate the representation of these processes against obser-
vations. A suitable measurement set-up required for such a study is described by Sun et al.
(2003), who used observations of sensible heat flux and radiative heat flux at two different
heights above the ground to estimate the heating rates of turbulence and radiation within
this atmospheric layer. Importantly, this would enable a comparison with model tendencies.
Expanding such a measurement network over a certain area would further allow compar-
ing the spatial distribution of different stability regimes between nature and model. Such a
processed-based validation would provide novel insights into the shortcomings of the repre-
sentation of the SBL and benefit model development for high latitudes.

**Acknowledgements** We firmly thank Bas van de Wiel and one anonymous reviewer for their valuable and constructive comments on the manuscript. We thank Roberta Pirazzini from the Finnish Meteorological Institute for preparing and providing the Sodankylä observations during YOPP SOP-NH1. We also thank Patrick Samuelsson for helpful discussions about the SURFEX scheme. We further thank Christiane Duscha for designing the schematic displayed in Fig. 1. This work was funded by the Norwegian Research Council through the ALERTNESS project in the POLARPROG programme (Project 280573). The model data that were created and used for this study are openly available from the THREDDS data server of the Norwegian Meteorological Institute (https://thredds.met.no/thredds/catalog/alertness/users/marvink/YOPP/catalog.html). The DDH data can be made available upon request. The Sodankylä observations can be accessed at https://litdb.fmi.fi/.

**Funding** Open access funding provided by University of Bergen (incl Haukeland University Hospital)

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
ment of a strong near-surface inversion in the vSBL can lead to a warming by radiation at the lowermost model level. Follow-up studies should test the impact of an increased vertical resolution on the process representation.

Both stability regimes coexist on small spatial scales of about 5 km during our study, most likely set by the model grid, and yield substantially different atmospheric temperature profiles. Such small-scale heterogeneity in the modelled SBL needs to be considered for validation studies, i.e. by investigating a sufficiently large sub-domain around the observational site.

The updated surface scheme, NEW, leads to a substantial increase in radiative cooling of the lower atmosphere due to a reduced effective emissivity of the surface which again leads to a larger radiative flux divergence. Yet, in terms of the investigated surface and near-surface temperature fields, NEW only exhibits pronounced changes to the vSBL. Here, the multi-layer snow scheme better captures the strong cooling of the upper layers in the snow pack compared to the single-layer snow scheme in REF. In the wSBL, however, this more efficient cooling of the surface is counteracted by the continuous heat supply of the turbulence scheme, yielding nearly unchanged surface and near-surface temperatures in both experiments. Such a varying and potentially compensating interplay between parametrization schemes is important to consider during model development, especially in regard to the common practice of model tuning.

A natural next step now is to validate the representation of these processes against observations. A suitable measurement set-up required for such a study is described by Sun et al. (2003), who used observations of sensible heat flux and radiative heat flux at two different heights above the ground to estimate the heating rates of turbulence and radiation within this atmospheric layer. Importantly, this would enable a comparison with model tendencies. Expanding such a measurement network over a certain area would further allow comparing the spatial distribution of different stability regimes between nature and model. Such a processed-based validation would provide novel insights into the shortcomings of the representation of the SBL and benefit model development for high latitudes.

**Acknowledgements** We firmly thank Bas van de Wiel and one anonymous reviewer for their valuable and constructive comments on the manuscript. We thank Roberta Pirazzini from the Finnish Meteorological Institute for preparing and providing the Sodankylä observations during YOPP SOP-NH1. We also thank Patrick Samuelsson for helpful discussions about the SURFEX scheme. We further thank Christiane Duscha for designing the schematic displayed in Fig. 1. This work was funded by the Norwegian Research Council through the ALERTNESS project in the POLARPROG programme (Project 280573). The model data that were created and used for this study are openly available from the THREDDS data server of the Norwegian Meteorological Institute (https://thredds.met.no/thredds/catalog/alertness/users/marvink/YOPP/catalog.html). The DDH data can be made available upon request. The Sodankylä observations can be accessed at https://litdb.fmi.fi/.

**Funding** Open access funding provided by University of Bergen (incl Haukeland University Hospital)

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
## References

André JC, Mahrt L (1982) The nocturnal surface inversion and influence of clear-air radiative cooling. J Atmos Sci 39(4):864–878

Arduini G, Balsamo G, Dutra E, Day JJ, Sandu I, Boussetta S, Haiden T (2019) Impact of a multi-layer snow scheme on near-surface weather forecasts. J Adv Model Earth Syst 11(12):4687–4710

Atlainkin E, Vihma T (2012) Evaluation of NWP results for wintertime nocturnal boundary-layer temperatures over Europe and Finland. Q J R Meteorol Soc 138(667):1440–1451

Baas P, de Roode SR, Lenderink G (2008) The scaling behaviour of a turbulent kinetic energy closure model for stably stratified conditions. Boundary-Layer Meteorol 127(1):17–36

Baas P, van de Wiel BJH, van der Linden SJA, Bosveld FC (2018) From near-neutral to strongly stratified: adequately modelling the clear-sky nocturnal boundary layer at Cabauw. Boundary-Layer Meteorol 166(2):217–238

Barrere M, Domine F, Decharme B, Morin S, Vionnet V, Lafaysse M (2017) Evaluating the performance of coupled snow-soil models in SURFEXv8 to simulate the permafrost thermal regime at a high Arctic site. Geosci Model Dev 10(9):3461–3479

Beesley JA, Bretherton CS, Jakob C, Andreas EL, Intrieri JM, Utal TA (2000) A comparison of cloud and boundary layer variables in the ECMWF forecast model with observations at surface heat budget of the arctic ocean (SHEBA) ice camp. J Geophys Res Atmos 105(D10):12337–12349

Bengtsson L, Tijm S, Váňa F, Svensson G (2012) Impact of flow-dependent horizontal diffusion on resolved convection in AROME. J Appl Meteorol Climatol 51(1):54–67

Bengtsson L, Andrae U, Aspelien T, Batrak Y, Calvo J, de Rooy W, Gleeson E, Hansen-Sass B, Homleid M, Hortal M, Ivarsson KI, Lenderink G, Niemelä S, Nielsen KP, Onvlee J, Rontu L, Samuelsson P, Muñoz DS, Subias A, Tijm S, Toll V, Yang X, Koltzow MO (2017) The HARMONIE–AROME model configuration in the ALADIN–HIRLAM NWP system. Mon Weather Rev 145(5):1919–1935

Boone A, Etchevers P (2001) An intercomparison of three snow schemes of varying complexity coupled to the same land surface model: local-scale evaluation at an Alpine site. J Hydrometeorol 2(4):374–394

Boone A, Samuelsson P, Gollvik S, Napoli A, Jarlan L, Brun E, Decharme B (2017) The interactions between soil–biosphere–atmosphere land surface model with a multi-energy balance (ISBA-MEB) option in SURFEXv8—Part 1: model description. Geosci Model Dev 10(2):843–872

Casati B, Haiden T, Brown B, Nurmi P, Lemieux JF (2017) Verification of environmental prediction in polar regions: recommendations for the year of polar prediction. WWRP 2017-1

Day JJ, Arduini G, Sandu I, Magnusson L, Beljaars A, Balsamo G, Rodwell M, Richardson D (2020) Measuring the impact of a new snow model using surface energy budget process relationships. J Adv Model Earth Syst 12(12):e2020MS002,144

de Rooy WC, Siebesma P, Baas P, Lenderink G, de Roode SR, de Vries H, van Meijgaard E, Meirink JF, Tijm S, van't Veen B (2022) Model development in practice: a comprehensive update to the boundary layer schemes in HARMONIE–AROME cycle 40. Geosci Model Dev 15(4):1513–1543

de Wiel BJHV, Vignon E, Baas P, van Hooijdonk IGS, van der Linden SJA, van Hooft JA, Bosveld FC, de Roode SR, Moene AF, Genthon C (2017) Regime transitions in near-surface temperature inversions: a conceptual model. J Atmos Sci 74(4):1057–1073

Deardorff JW (1980) Stratocumulus-capped mixed layers derived from a three-dimensional model. Boundary-Layer Meteorol 18(4):495–527

Decharme B, Boone A, Delire C, Noilhan J (2011) Local evaluation of the interaction between soil biosphere atmosphere soil multilayer diffusion scheme using four pedotransfer functions. J Geophys Res Atmos. https://doi.org/10.1029/2011JD016002

Decharme B, Brun E, Boone A, Delire C, Le Moigne P, Morin S (2016) Impacts of snow and organic soils parameterization on northern Eurasian soil temperature profiles simulated by the ISBA land surface model. Cryosphere 10(2):853–877

Derbyshire SH (1999) Boundary-layer decoupling over cold surfaces as a physical boundary-instability. Boundary-Layer Meteorol 90(2):297–325

Douville H, Royer JF, Mahfouf JF (1995a) A new snow parameterization for the Météo-France climate model. Clim Dyn 12(1):21–35

Douville H, Royer JF, Mahfouf JF (1995b) A new snow parameterization for the météo-France climate model. Clim Dyn 12(1):37–52

Edwards JM (2009a) Radiative processes in the stable boundary layer: Part I. Radiative aspects. Boundary-Layer Meteorol 131(2):105

Edwards JM (2009b) Radiative processes in the stable boundary layer: Part II. The development of the nocturnal boundary layer. Boundary-Layer Meteorol 131(2):127–146
References

André JC, Mahrt L (1982) The nocturnal surface inversion and influence of clear-air radiative cooling. J Atmos Sci 39(4):864–878

Arduini G, Balsamo G, Dutra E, Day JJ, Sandu I, Boussetta S, Haiden T (2019) Impact of a multi-layer snow scheme on near-surface weather forecasts. J Adv Model Earth Syst 11(12):4687–4710

Atlainkin E, Vihma T (2012) Evaluation of NWP results for wintertime nocturnal boundary-layer temperatures over Europe and Finland. Q J R Meteorol Soc 138(667):1440–1451

Baas P, de Roode SR, Lenderink G (2008) The scaling behaviour of a turbulent kinetic energy closure model for stably stratified conditions. Boundary-Layer Meteorol 127(1):17–36

Baas P, van de Wiel BJH, van der Linden SJA, Bosveld FC (2018) From near-neutral to strongly stratified: adequately modelling the clear-sky nocturnal boundary layer at Cabauw. Boundary-Layer Meteorol 166(2):217–238

Barrere M, Domine F, Decharme B, Morin S, Vionnet V, Lafaysse M (2017) Evaluating the performance of coupled snow-soil models in SURFEXv8 to simulate the permafrost thermal regime at a high Arctic site. Geosci Model Dev 10(9):3461–3479

Beesley JA, Bretherton CS, Jakob C, Andreas EL, Intrieri JM, Utal TA (2000) A comparison of cloud and boundary layer variables in the ECMWF forecast model with observations at surface heat budget of the arctic ocean (SHEBA) ice camp. J Geophys Res Atmos 105(D10):12337–12349

Bengtsson L, Tijm S, Váňa F, Svensson G (2012) Impact of flow-dependent horizontal diffusion on resolved convection in AROME. J Appl Meteorol Climatol 51(1):54–67

Bengtsson L, Andrae U, Aspelien T, Batrak Y, Calvo J, de Rooy W, Gleeson E, Hansen-Sass B, Homleid M, Hortal M, Ivarsson KI, Lenderink G, Niemelä S, Nielsen KP, Onvlee J, Rontu L, Samuelsson P, Muñoz DS, Subias A, Tijm S, Toll V, Yang X, Koltzow MO (2017) The HARMONIE–AROME model configuration in the ALADIN–HIRLAM NWP system. Mon Weather Rev 145(5):1919–1935

Boone A, Etchevers P (2001) An intercomparison of three snow schemes of varying complexity coupled to the same land surface model: local-scale evaluation at an Alpine site. J Hydrometeorol 2(4):374–394

Boone A, Samuelsson P, Gollvik S, Napoli A, Jarlan L, Brun E, Decharme B (2017) The interactions between soil–biosphere–atmosphere land surface model with a multi-energy balance (ISBA-MEB) option in SURFEXv8—Part 1: model description. Geosci Model Dev 10(2):843–872

Casati B, Haiden T, Brown B, Nurmi P, Lemieux JF (2017) Verification of environmental prediction in polar regions: recommendations for the year of polar prediction. WWRP 2017-1

Day JJ, Arduini G, Sandu I, Magnusson L, Beljaars A, Balsamo G, Rodwell M, Richardson D (2020) Measuring the impact of a new snow model using surface energy budget process relationships. J Adv Model Earth Syst 12(12):e2020MS002,144

de Rooy WC, Siebesma P, Baas P, Lenderink G, de Roode SR, de Vries H, van Meijgaard E, Meirink JF, Tijm S, van't Veen B (2022) Model development in practice: a comprehensive update to the boundary layer schemes in HARMONIE–AROME cycle 40. Geosci Model Dev 15(4):1513–1543

de Wiel BJHV, Vignon E, Baas P, van Hooijdonk IGS, van der Linden SJA, van Hooft JA, Bosveld FC, de Roode SR, Moene AF, Genthon C (2017) Regime transitions in near-surface temperature inversions: a conceptual model. J Atmos Sci 74(4):1057–1073

Deardorff JW (1980) Stratocumulus-capped mixed layers derived from a three-dimensional model. Boundary-Layer Meteorol 18(4):495–527

Decharme B, Boone A, Delire C, Noilhan J (2011) Local evaluation of the interaction between soil biosphere atmosphere soil multilayer diffusion scheme using four pedotransfer functions. J Geophys Res Atmos. https://doi.org/10.1029/2011JD016002

Decharme B, Brun E, Boone A, Delire C, Le Moigne P, Morin S (2016) Impacts of snow and organic soils parameterization on northern Eurasian soil temperature profiles simulated by the ISBA land surface model. Cryosphere 10(2):853–877

Derbyshire SH (1999) Boundary-layer decoupling over cold surfaces as a physical boundary-instability. Boundary-Layer Meteorol 90(2):297–325

Douville H, Royer JF, Mahfouf JF (1995a) A new snow parameterization for the Météo-France climate model. Clim Dyn 12(1):21–35

Douville H, Royer JF, Mahfouf JF (1995b) A new snow parameterization for the météo-France climate model. Clim Dyn 12(1):37–52

Edwards JM (2009a) Radiative processes in the stable boundary layer: Part I. Radiative aspects. Boundary-Layer Meteorol 131(2):105

Edwards JM (2009b) Radiative processes in the stable boundary layer: Part II. The development of the nocturnal boundary layer. Boundary-Layer Meteorol 131(2):127–146
Esau I, Tolstykh M, Fadeev R, Shashkin V, Makhnorylova S, Miles V, Melnikov V (2018) Systematic errors in northern Eurasian short-term weather forecasts induced by atmospheric boundary layer thickness. *Environ Res Lett* 13(12):125009

Fouquart Y, Bonnel B (1980) Computations of solar heating of the earth's atmosphere—a new parameterization. *Beitr Phys Atmos* 53:35–62

Haiden T, Sandu I, Gianpaolo B, Gabriele A, Anton B (2018) Addressing biases in near-surface forecasts. *ECMWF Newsl* 157:20–25

Holtslag B (2006) Preface: GEWEX atmospheric boundary-layer study (GABLS) on stable boundary layers. *Boundary-Layer Meteorol* 118(2):243–246

Holtslag AAM, Svensson G, Baas P, Basu S, Beare B, Beljaars ACM, Bosveld FC, Cuxart J, Lindvall J, Steeneveld GJ, Tjernström M, Wiel BJHVD (2013) Stable atmospheric boundary layers and diurnal cycles: challenges for weather and climate models. *Bull Am Meteorol Soc* 94(11):1691–1706

Kähnert M, Sodemann H, de Rooy WC, Valkonen TM (2021) On the utility of individual tendency output: revealing interactions between parameterized processes during a marine cold air outbreak. *Weather Forecast* 36(6):1985–2000

Kanamitsu M, DeHaan L (2011) The added value index: a new metric to quantify the added value of regional models. *J Geophys Res Atmos*. https://doi.org/10.1029/2011JD015597

Kangas M, Rontu L, Fortelius C, Aurela M, Poikonen A (2016) Weather model verification using Sodankylä mast measurements. *Methods Data Syst* 5(1):75–84

Køltzow M, Casati B, Bazile E, Haiden T, Valkonen T (2019) A NWP model inter-comparison of surface weather parameters in the European Arctic during the year of polar prediction special observing period northern hemisphere 1. *Weather Forecast* 34:959–83

Krishna TBPSRV, Sharan M, Gopalakrishnan SG, Aditi (2003) Mean structure of the nocturnal boundary layer under strong and weak wind conditions: EPRI case study. *J Appl Meteorol Climatol* 42(7):952–969

Le Moigne P (2009) SURFEX scientific documentation. *CNRM Tech Rep p* 211 pp

Lenderink G, Holtslag AAM (2004) An updated length-scale formulation for turbulent mixing in clear and cloudy boundary layers. *Q J R Meteorol Soc* 130(604):3405–3427

McNider RT, England DE, Friedman MJ, Shi X (1995) Predictability of the stable atmospheric boundary layer. *J Atmos Sci* 52(10):1602–1614

Météo-France (2019) Diagnostics in horizontal domains (DDH)—variables and budget equations, in horizontal mean ARPEGE, ALADIN and AROME models. *Météo-France, Specific Documentation*

Mlawer EJ, Taubman SJ, Brown PD, Iacono MJ, Clough SA (1997) Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlated-k model for the longwave. *J Geophys Res Atmos* 102(D14):16,663–16,682

Müller M, Batrak Y, Kristiansen J, Køltzow M, Noer G, Korosov A (2017) Characteristics of a convective-scale weather forecasting system for the European arctic. *Mon Weather Rev* 145(12):4771–4787

Napoly A, Boone A, Samuelsson P, Gollvik S, Martin E, Seferian R, Carrer D, Decharme B, Jarlan L (2017) The interactions between soil-biosphere-atmosphere (ISBA) land surface model multi-energy balance (MEB) option in SURFEXv8—Part 2: introduction of a litter formulation and model evaluation for local-scale forest sites. *Geosci Model Dev* 10(4):1621–1644

Noilhan J, Planton S (1989) A simple parameterization of land surface processes for meteorological models. *Mon Weather Rev* 117(3):536–549

Sandu I, Beljaars A, Bechtold P, Mauritsen T, Balsamo G (2013) Why is it so difficult to represent stably stratified conditions in numerical weather prediction (NWP) models? *J Adv Model Earth Syst* 5(2):117–133

Seity Y, Brousseau P, Malardel S, Hello G, Bénard P, Bouttier F, Lac C, Masson V (2011) The AROME-France convective-scale operational model. *Mon Weather Rev* 139(3):976–991

Sodemann H, Foken T (2004) Empirical evaluation of an extended similarity theory for the stably stratified atmospheric surface layer. *Q J R Meteorol Soc* 130(602):2665–2671

Steeneveld GJ, van de Wiel BJH, Holtslag AAM (2006) Modeling the evolution of the atmospheric boundary layer coupled to the land surface for three contrasting nights in CASES-99. *J Atmos Sci* 63(3):920–935

Sterk HAM, Steeneveld GJ, Holtslag AAM (2013) The role of snow-surface coupling, radiation, and turbulent mixing in modeling a stable boundary layer over Arctic sea ice. *J Geophys Res Atmos* 118(3):1199–1217

Sun S, Jin J, Xue Y (1999) A simple snow–atmosphere–soil transfer model. *J Geophys Res Atmos* 104(D16):19,587–19,597

Sun J, Burns SP, Delany AC, Oncley SP, Horst TW, Lenschow DH (2003) Heat balance in the nocturnal boundary layer during CASES-99. *J Appl Meteorol Climatol* 42(11):1649–1666

Svensson G, Holtslag B (2006) Single column modeling of the diurnal cycle based on CASES99 data—GABLS second intercomparison project. In: 17th symposium on boundary layers and turbulence, San Diego, CA, USA, 22–25 May 2006 *Amer Meteorol Soc*, Boston p 8.1
Esau I, Tolstykh M, Fadeev R, Shashkin V, Makhnorylova S, Miles V, Melnikov V (2018) Systematic errors in northern Eurasian short-term weather forecasts induced by atmospheric boundary layer thickness. *Environ Res Lett* 13(12):125009

Fouquart Y, Bonnel B (1980) Computations of solar heating of the earth's atmosphere—a new parameterization. *Beitr Phys Atmos* 53:35–62

Haiden T, Sandu I, Gianpaolo B, Gabriele A, Anton B (2018) Addressing biases in near-surface forecasts. *ECMWF Newsl* 157:20–25

Holtslag B (2006) Preface: GEWEX atmospheric boundary-layer study (GABLS) on stable boundary layers. *Boundary-Layer Meteorol* 118(2):243–246

Holtslag AAM, Svensson G, Baas P, Basu S, Beare B, Beljaars ACM, Bosveld FC, Cuxart J, Lindvall J, Steeneveld GJ, Tjernström M, Wiel BJHVD (2013) Stable atmospheric boundary layers and diurnal cycles: challenges for weather and climate models. *Bull Am Meteorol Soc* 94(11):1691–1706

Kähnert M, Sodemann H, de Rooy WC, Valkonen TM (2021) On the utility of individual tendency output: revealing interactions between parameterized processes during a marine cold air outbreak. *Weather Forecast* 36(6):1985–2000

Kanamitsu M, DeHaan L (2011) The added value index: a new metric to quantify the added value of regional models. *J Geophys Res Atmos.* https://doi.org/10.1029/2011JD015597

Kangas M, Rontu L, Fortelius C, Aurela M, Poikonen A (2016) Weather model verification using Sodankylä mast measurements. *Methods Data Syst* 5(1):75–84

Køltzow M, Casati B, Bazile E, Haiden T, Valkonen T (2019) A NWP model inter-comparison of surface weather parameters in the European Arctic during the year of polar prediction special observing period northern hemisphere 1. *Weather Forecast* 34:959–83

Krishna TBPSRV, Sharan M, Gopalakrishnan SG, Aditi (2003) Mean structure of the nocturnal boundary layer under strong and weak wind conditions: EPRI case study. *J Appl Meteorol Climatol* 42(7):952–969

Le Moigne P (2009) SURFEX scientific documentation. *CNRM Tech Rep p* 211 pp

Lenderink G, Holtslag AAM (2004) An updated length-scale formulation for turbulent mixing in clear and cloudy boundary layers. *Q J R Meteorol Soc* 130(604):3405–3427

McNider RT, England DE, Friedman MJ, Shi X (1995) Predictability of the stable atmospheric boundary layer. *J Atmos Sci* 52(10):1602–1614

Météo-France (2019) Diagnostics in horizontal domains (DDH)—variables and budget equations, in horizontal mean ARPEGE, ALADIN and AROME models. *Météo-France, Specific Documentation*

Mlawer EJ, Taubman SJ, Brown PD, Iacono MJ, Clough SA (1997) Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlated-k model for the longwave. *J Geophys Res Atmos* 102(D14):16,663–16,682

Müller M, Batrak Y, Kristiansen J, Køltzow M, Noer G, Korosov A (2017) Characteristics of a convective-scale weather forecasting system for the European arctic. *Mon Weather Rev* 145(12):4771–4787

Napoly A, Boone A, Samuelsson P, Gollvik S, Martin E, Seferian R, Carrer D, Decharme B, Jarlan L (2017) The interactions between soil-biosphere-atmosphere (ISBA) land surface model multi-energy balance (MEB) option in SURFEXv8—Part 2: introduction of a litter formulation and model evaluation for local-scale forest sites. *Geosci Model Dev* 10(4):1621–1644

Noilhan J, Planton S (1989) A simple parameterization of land surface processes for meteorological models. *Mon Weather Rev* 117(3):536–549

Sandu I, Beljaars A, Bechtold P, Mauritsen T, Balsamo G (2013) Why is it so difficult to represent stably stratified conditions in numerical weather prediction (NWP) models? *J Adv Model Earth Syst* 5(2):117–133

Seity Y, Brousseau P, Malardel S, Hello G, Bénard P, Bouttier F, Lac C, Masson V (2011) The AROME-France convective-scale operational model. *Mon Weather Rev* 139(3):976–991

Sodemann H, Foken T (2004) Empirical evaluation of an extended similarity theory for the stably stratified atmospheric surface layer. *Q J R Meteorol Soc* 130(602):2665–2671

Steeneveld GJ, van de Wiel BJH, Holtslag AAM (2006) Modeling the evolution of the atmospheric boundary layer coupled to the land surface for three contrasting nights in CASES-99. *J Atmos Sci* 63(3):920–935

Sterk HAM, Steeneveld GJ, Holtslag AAM (2013) The role of snow-surface coupling, radiation, and turbulent mixing in modeling a stable boundary layer over Arctic sea ice. *J Geophys Res Atmos* 118(3):1199–1217

Sun S, Jin J, Xue Y (1999) A simple snow–atmosphere–soil transfer model. *J Geophys Res Atmos* 104(D16):19,587–19,597

Sun J, Burns SP, Delany AC, Oncley SP, Horst TW, Lenschow DH (2003) Heat balance in the nocturnal boundary layer during CASES-99. *J Appl Meteorol Climatol* 42(11):1649–1666

Svensson G, Holtslag B (2006) Single column modeling of the diurnal cycle based on CASES99 data—GABLS second intercomparison project. In: 17th symposium on boundary layers and turbulence, San Diego, CA, USA, 22–25 May 2006 *Amer Meteorol Soc*, Boston p 8.1
Svensson G, Holtslag AAM (2009) Analysis of model results for the turning of the wind and related momentum fluxes in the stable boundary layer. Boundary-Layer Meteorol 132(2):261–277

Tjernström M, Žagar M, Svensson G, Cassano JJ, Pfeifer S, Rinke A, Wyser K, Dethloff K, Jones C, Semmler T, Shaw M (2005) Modelling the arctic boundary layer: an evaluation of six Arcmip regional-scale models using data from the Sheba project. Boundary-Layer Meteorol 117(2):337–381

Ulden APV, Holtslag AAM (1985) Estimation of atmospheric boundary layer parameters for diffusion applications. J Appl Meteorol Climatol 24(11):1196–1207

van der Linden SJA, Kruis MT, Hartogensis OK, Moene AF, Bosveld FC, van de Wiel BJH (2022) Heat transfer through grass: a diffusive approach. Boundary-Layer Meteorol 184:251–276

Váňa F, Bénard P, Geleyn JF, Simon A, Seity Y (2008) Semi-Lagrangian advection scheme with controlled damping: an alternative to nonlinear horizontal diffusion in a numerical weather prediction model. Q J R Meteorol Soc 134(631):523–537

Vignon E, van de Wiel BJH, van Hooijdonk IGS, Genthon C, van der Linden SJA, van Hooft JA, Baas P, Maurel W, Traullé O, Casasanta G (2017) Stable boundary-layer regimes at Dome C, Antarctica: observation and analysis. Q J R Meteorol Soc 143(704):1241–1253

Vihma T, Pirazzini R, Fer I, Renfrew IA, Sedlar J, Tjernström M, Lüpkes C, Nygård T, Notz D, Weiss J, Marsan D, Cheng B, Birnbaum G, Gerland S, Chechin D, Gascard JC (2014) Advances in understanding and parameterization of small-scale physical processes in the marine Arctic climate system: a review. Atmos Chem Phys 14(17):9403–9450

Viterbo P, Beljaars A, Mahfouf JF, Teixeira J (1999) The representation of soil moisture freezing and its impact on the stable boundary layer. Q J R Meteorol Soc 125(559):2401–2426

Vogelezang DHP, Holtslag AAM (1996) Evaluation and model impacts of alternative boundary-layer height formulations. Boundary-Layer Meteorol 81(3):245–269

Yen YC (1981) Review of the thermal properties of snow, ice and sea ice. Tech Rep, Cold Reg Res Eng Lab pp 81–10

Zilitinkevich SS, Esau IN (2005) Resistance and heat-transfer laws for stable and neutral planetary boundary layers: old theory advanced and re-evaluated. Q J R Meteorol Soc 131(609):1863–1892

Zilitinkevich SS, Elperin T, Kleeorin N, Rogachevskii I, Esau I, Mauritsen T, Miles MW (2008) Turbulence energetics in stably stratified geophysical flows: strong and weak mixing regimes. Q J R Meteorol Soc 134(633):793–799

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Svensson G, Holtslag AAM (2009) Analysis of model results for the turning of the wind and related momentum fluxes in the stable boundary layer. Boundary-Layer Meteorol 132(2):261–277

Tjernström M, Žagar M, Svensson G, Cassano JJ, Pfeifer S, Rinke A, Wyser K, Dethloff K, Jones C, Semmler T, Shaw M (2005) Modelling the arctic boundary layer: an evaluation of six Arcmip regional-scale models using data from the Sheba project. Boundary-Layer Meteorol 117(2):337–381

Ulden APV, Holtslag AAM (1985) Estimation of atmospheric boundary layer parameters for diffusion applications. J Appl Meteorol Climatol 24(11):1196–1207

van der Linden SJA, Kruis MT, Hartogensis OK, Moene AF, Bosveld FC, van de Wiel BJH (2022) Heat transfer through grass: a diffusive approach. Boundary-Layer Meteorol 184:251–276

Váňa F, Bénard P, Geleyn JF, Simon A, Seity Y (2008) Semi-Lagrangian advection scheme with controlled damping: an alternative to nonlinear horizontal diffusion in a numerical weather prediction model. Q J R Meteorol Soc 134(631):523–537

Vignon E, van de Wiel BJH, van Hooijdonk IGS, Genthon C, van der Linden SJA, van Hooft JA, Baas P, Maurel W, Traullé O, Casasanta G (2017) Stable boundary-layer regimes at Dome C, Antarctica: observation and analysis. Q J R Meteorol Soc 143(704):1241–1253

Vihma T, Pirazzini R, Fer I, Renfrew IA, Sedlar J, Tjernström M, Lüpkes C, Nygård T, Notz D, Weiss J, Marsan D, Cheng B, Birnbaum G, Gerland S, Chechin D, Gascard JC (2014) Advances in understanding and parameterization of small-scale physical processes in the marine Arctic climate system: a review. Atmos Chem Phys 14(17):9403–9450

Viterbo P, Beljaars A, Mahfouf JF, Teixeira J (1999) The representation of soil moisture freezing and its impact on the stable boundary layer. Q J R Meteorol Soc 125(559):2401–2426

Vogelezang DHP, Holtslag AAM (1996) Evaluation and model impacts of alternative boundary-layer height formulations. Boundary-Layer Meteorol 81(3):245–269

Yen YC (1981) Review of the thermal properties of snow, ice and sea ice. Tech Rep, Cold Reg Res Eng Lab pp 81–10

Zilitinkevich SS, Esau IN (2005) Resistance and heat-transfer laws for stable and neutral planetary boundary layers: old theory advanced and re-evaluated. Q J R Meteorol Soc 131(609):1863–1892

Zilitinkevich SS, Elperin T, Kleeorin N, Rogachevskii I, Esau I, Mauritsen T, Miles MW (2008) Turbulence energetics in stably stratified geophysical flows: strong and weak mixing regimes. Q J R Meteorol Soc 134(633):793–799

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.