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

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.


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 withing 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 limitedarea 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.

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 w θ is taken proportionally to an eddy diffusivity coefficient for heat K h , and the vertical gradient of potential temperature θ : 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: 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): Here, c n and c s are model internal parameters, κ 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).

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: 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, ∂(F + − F − )/∂z-the divergence of the net upward radiative flux, and ρ-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.

Sensitivity Experiments
The representation of the surface is of particular importance for modelling the SBL 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. 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-vegetationsnow 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 . 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 snowcovered 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). 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×10 −7 m 2 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).

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 needleleaf forest (Fig. 2c), influencing the estimated roughness length (Fig. 2d).

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 3m (Fig. 3b). Wind speeds are comparably low during the nocturnal periods, dropping from 3 to 1 m s −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 Fig. 2 Physiography of AROME-Arctic for the extracted 35 × 35 km 2 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 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.

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. 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 ∂ 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), 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 largescale advection (or subsidence), as well as horizontal diffusion, can play an important role. Therefore, spatial information is required.

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 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 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 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).

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 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 nearsurface 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  (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.

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  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. 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 −1 and warming the near-surface levels (< 50 m) with 1-2 K h −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 K h −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 K h −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.

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 −2 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 −2 ). 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.

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.

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 wellfounded 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.

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.

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 . 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.

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 an 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.