Sensitivity to changes in the surface‐layer turbulence parameterization for stable conditions in winter: A case study with a regional climate model over the Arctic

The modeling of the atmospheric boundary layer over sea ice is still challenging because of the complex interaction between clouds, radiation and turbulence over the often inhomogeneous sea ice cover. There is still much uncertainty concerning sea ice roughness, near‐surface thermal stability and related processes, and their accurate parameterization. Here, a regional Arctic climate model forced by ERA‐Interim data was used to test the sensitivity of climate simulations to a modified surface flux parameterization for wintertime conditions over the Arctic. The reference parameterization as well as the modified one is based on Monin–Obukhov similarity theory, but different roughness lengths were prescribed and the stability dependence of the transfer coefficients for momentum, heat and moisture differed from each other. The modified parameterization accounts for the most comprehensive observations that are presently available over sea ice in the inner Arctic. Independent of the parameterization used, the model was able to reproduce the two observed dominant winter states with respect to cloud cover and longwave radiation. A stepwise use of the different parameterization assumptions showed that modifications of both surface roughness and stability dependence had a considerable impact on quantities such as air pressure, wind and near‐surface turbulent fluxes. However, the reduction of surface roughness to values agreeing with those observed during the Surface Heat Budget of the Arctic Ocean campaign led to an improvement in the western Arctic, while the modified stability parameterization had only a minor impact. The latter could be traced back to the model's underestimation of the strength of stability over sea ice. Future work should concentrate on possible reasons for this underestimation and on the question of generality of the results for other climate models.


| INTRODUCTION
In climate and weather prediction models, the nearsurface turbulent fluxes in the atmospheric boundary layer (ABL) and related transfer coefficients are usually parametrized on the basis of the Monin-Obukhov similarity theory (MOST, Monin and Obukhov, 1954). It relates the mean profiles of meteorological quantities to the respective surface fluxes of momentum, heat and moisture and to surface parameters. MOST provides transfer coefficients, which are given as a product of coefficients for neutral stratification and of stability correction functions (SCF, see Section 2.2). The scatter between results of available SCFs is large and increases with stability. Gryanik and Lüpkes (2018) (in the following referred to as GL18) as well as Gryanik et al. (2020) showed that for strong stability the often-applied SCFs of Louis (1979) and Louis et al. (1982) (LTG82 in the following) overestimate turbulent fluxes by far. The increased mixing is one of the reasons why climate and weather prediction models tend to overestimate ABL thickness as compared to measurements (Savijärvi, 2009).
One possible strategy to overcome this difficulty was suggested by GL18. They propose new SCFs for stable stratification to improve the model's representation of surface turbulent fluxes of momentum and heat. The new SCFs are based on MOST stability functions of Grachev et al. (2007), which were established using data from the Surface Heat Budget of the Arctic Ocean (SHEBA) campaign (Uttal et al., 2002;Andreas et al., 2007). They represent still the most comprehensive data set available for Arctic conditions. Transfer coefficients and the new SCFs depend on the surface roughness. GL18 formulated their functions with roughness lengths for the inner Arctic sea ice based on well accepted measurements (see e.g., Lüpkes et al., 2012;Vihma et al., 2014;Lüpkes and Gryanik, 2015;Elvidge et al., 2016). These values are smaller than those often used in climate models. As found, for example, by Renfrew et al. (2019) the surface roughness of sea ice can strongly affect the results of atmospheric models in polar regions.
The purpose of this study is to undertake a first assessment of an implementation of the GL18 SCFs in a regional climate model over the Arctic, with the focus on winter, the time of predominantly stable atmospheric conditions. We consider the year 2007, a year of low Arctic summer sea ice extent, so that the near-surface stratification might still have been less stable in early winter than in other years. This choice is motivated first of all by the data availability from the North Pole-35 (NP-35) sea ice drift station in this year Makhotina et al., 2019).
The guiding questions are: What is the impact of changes from default LTG82 SCFs to GL18 ones on the surface turbulent fluxes and ABL structure, and consequently on regional atmospheric circulation? How large is the sensitivity of the simulation results when the GL18 SCFs are implemented in different levels of complexity? Can we arrive at an improved skill in the simulations by the most complete GL18 implementation? 2 | MODEL, PARAMETERIZATION, AND SIMULATIONS 2.1 | Regional model and simulations Simulations were performed with HIRHAM5 (Christensen et al., 2007), a regional atmospheric climate model that was run over the circum-Arctic domain approximately north of 65 N (see Figure S1). The model was applied with a horizontal resolution of 0.25 (ca. 27 km) and 40 vertical levels (10 levels within the lowest 1 km, lowest model level at ca. 10 m). The physical parameterizations are based on ECHAM5 (Roeckner et al., 2003). The model configuration corresponds to the HIRHAM5 description given by Dorn et al. (2019).
Previous studies have shown that although the model can reproduce the basic observed spatiotemporal patterns of atmospheric circulation, the simulations show limited skill in ABL processes (Tjernström et al., 2005;Sedlar et al., 2020). This calls for improvement of the representation of surface turbulent fluxes of momentum and heat. Accordingly, sensitivity simulations were set up using parameterization options explained in Section 2.2. They represent a step-by-step transition from the default scheme to the modified one. The procedure provides more transparent diagnostics of the sensitivity of the model results to modifications. The approach introduced can thus guide potential model improvement in the future.
The simulations, initialized and forced at the boundaries by ERA-Interim (Dee et al., 2011), were run for the year 2007, a year of low summer Arctic sea ice extent. To account for internal variability, ensemble simulations with 10 members were carried out for each of the parameterization options. The individual ensemble members differ only in their initialization date, which was shifted successively by 6 h, that is, the first member was initialized with data from 0000 UTC on January 1, 2007 and each following member of the ensemble was initialized with data 6 h before the previous member. All simulations run until 0000 UTC on January 1, 2008.

| Parameterization of surface fluxes for stable atmospheric conditions
This subsection provides some basic background information, which is needed to follow the numerical experiments, while we refer to the given references for details of the parameterizations. Like in most climate and numerical weather prediction models, the near-surface turbulent fluxes of momentum, heat and moisture and related transfer coefficients are parameterized in HIRHAM5 by default on the basis of MOST applying a noniterative scheme by Louis et al. (1982). This means that the turbulent fluxes are calculated via bulk formulae between the surface and the first model level. Following MOST the surface momentum flux components τ x and τ y as well as the surface sensible heat flux H can be parameterized according to where ρ is the air density and U ! (z) is the horizontal wind vector at height z of the lowest model level with components U and V. c p is the specific heat capacity, θ(z) is the potential temperature at z and index 0 refers to the surface values. C m and C h are the momentum transfer (or drag) and heat transfer coefficients, respectively. They are calculated as with the neutral transfer coefficients C m,n and C h,n , which solely depend on roughness parameters z/z 0 , z/z t with length scales z 0 for momentum and z t for heat. The nondimensional functions f m/h represent normalized transfer coefficients and depend on the bulk Richardson number Ri b and also on z/z 0 and z/z t (e.g., GL18). We will call them stability correction functions (SCFs) in the following to stress their dependence on stability. We stress that this notation should not be mixed up with the MOST stability correction functions (e.g., Gryanik et al., 2020). Equations (1)-(5) as well as similar formulations for humidity form the basic equations for all parameterizations of the turbulent surface fluxes used here. Differences between the schemes refer to the functional form of SCFs and to the prescribed values of z 0 and z t .
The default HIRHAM5 setup uses the reference parameterization, in the following called P0, which uses SCFs as proposed by Louis et al. (1982). P0 depends on Ri b only and applies a constant value of z 0 = 10 À3 m, and the same value for z t and humidity z q over sea ice (in a corresponding equation for humidity transport), over the entire sea ice covered region.
This parameterization has several limitations for Arctic applications. Particularly, in the range Ri b > 0.1 the SCFs show large deviations from those obtained from the SHEBA campaign measurements (Gryanik et al., 2020). Specifically, C m,h values obtained by the original Louis et al. (1982) SCFs are too high for Ri b > 0.05. Also, the values used for z 0 and z t overestimate average values obtained from SHEBA (Lüpkes et al., 2012;Castellani et al., 2014;Lüpkes and Gryanik, 2015). Moreover, the ratio z 0 /z t = 1 contradicts findings by Andreas et al. (2010) who showed that this ratio is mostly much smaller than 1. For this reason, we apply modified parameterizations in HIRHAM5, hereafter referred to as P1 À P3 (Table 1).
P1 refers to the new GL18 scheme with SCFs that represent an essentially improved fit to the SHEBA data. This scheme can be run with variable surface roughness, but as a first step and to reduce complexity, constant values (z 0 = 3.3 Â 10 À4 m and z t = z q = 0.66 Â 10 À4 m) were used. These values are representative for Arctic sea ice including SHEBA conditions (Castellani et al., 2014;Lüpkes and Gryanik, 2015). One should keep in mind, however, that sea ice roughness varies a lot, roughly between 7 Â 10 À6 m and 5 Â 10 À2 m (e.g., Guest and Davidson, 1987;Andreas et al., 2010;Lüpkes et al., 2012;Elvidge et al., 2016) and it is justified to test the sensitivity in a wide range of values, including the values in P0. The specific choice in this work aims to clarify the type of feedbacks when the parameters mentioned are varied rather than to define a final optimum setup.
While P1 differs from P0 in both the SCFs and the roughness, P2 differs from P0 only in the new representative roughness parameters and P3 only in the new SCFs (Table 1). The modifications referring to SCFs were applied over sea ice for stable stratification, while the roughness modifications were applied independent of stratification. All experiments applied z t = z q , no changes of z 0 , z t , and z q were made over water and land, and the original equations for the stability functions remained unchanged for unstable conditions.

| RESULTS
The analysis is focused on winter, the time of predominantly stable conditions in the Arctic atmosphere. In Section 3.1, we present changes in the monthly mean patterns of turbulent fluxes and atmospheric circulation for December mean fields as a representative example for the other winter months (Schneider, 2020). Section 3.2 represents an analysis of changes by considering selected ABL process relationships based on 6-hourly data.

| Monthly mean turbulent fluxes and atmospheric circulation
The parameterization changes by GL18 aim for significantly smaller mean values of the transfer coefficients C h and C m over sea ice. Figure 1 shows that the changes applied to the transfer coefficients in P1 and P2 relative to P0 are dominated by changes in roughness because the largest changes appear for near-neutral conditions (0 ≤ Ri b < 0.02), for which z 0 and z t have the dominant impact on the transfer coefficients. However, specifically for P1, the largest absolute change in C m is at Ri b = 0.15-0.2, indicating additional changes from stability effects. Changes in SCFs influence the transfer coefficients only for Ri b ≳ 0.05, which is clearly visible in the difference curves P0 minus P3 (Figure 1a, Figure S1). The probability density distribution of modeled Ri b (Figure 1b) resulting from the different runs shows that near-neutral conditions (0 ≤ Ri b < 0.02) are dominating (48% of the time in December for P0) (Figure 1, Table S1). 24% of all cases are weakly or strongly stable and about 28% of all cases are slightly unstable. The terms weakly and strongly stable refer to the stability regimes with continuous turbulence and intermittent turbulence, respectively, according to Grachev et al. (2013) and Gryanik et al. (2021). Observations from SHEBA show, however, a much larger percentage of cases in the stable regimes (55%), a considerably smaller percentage (28%) in the near-neutral regime, and about 17% in the unstable regime (Table S1).
As shown in Figure S1 the transfer coefficients largely change for P1 and P2 compared to P0. The domain averaged C m and C h values for P1 (P2) are 23% and 32% (21% and 32%) smaller than for P0. Corresponding values for P3 amount to only 3% (C m ) and 2% (C h ). The decrease for P3 would probably be much stronger if very stable conditions were modeled as often as observed (Figure 1, Table S1).
The changes in sensible heat flux, where negative values correspond to upward fluxes, over sea ice in the T A B L E 1 Model experiment setups for the different levels of surface layer parameterization change P0 ("Control") Original 1.0 Â 10 À3 1.0 Â 10 À3 Variable P1 ("New") New 3.3 Â 10 À4 6.6 Â 10 À5 Variable P2 Original 3.3 Â 10 À4 6.6 Â 10 À5 Variable P3 New 1.0 Â 10 À3 1.0 Â 10 À3 Variable Note: Given are the applied values for the stability correction functions f m/h of momentum and heat transfer for stable conditions over sea ice, surface roughness length over sea ice z 0 , temperature and moisture roughness lengths over sea ice z t and z q , and bulk Richardson number over sea ice Ri b for stable conditions. The P0 setting represents the standard HIRHAM5 ("Control"), P1 setting represents the new GL18 scheme, while P2, P3 represent stepwise changes in between. See Section 2 for details.  Table S1 for the frequency of occurrence of the stability regimes). The PDF is based on 6-hourly data and grid points with a sea ice concentration of 95% and higher. Ri b was calculated over sea ice using data at the lowest model level at around 9-10 m and at the surface. The dashed black line shows the PDF of the SHEBA data for December 1997 (valid for measurements at 8.9 m) central Arctic relative to P0 are small and mostly insignificant ( Figure 2), but not necessarily insignificant in their impact. A 1 WÁm À2 flux change equates to a 10 cm ice melt in a year, a significant fraction of the ice budget (Bourassa et al., 2013). With absolute differences of about 3 WÁm À2 they correspond to relative differences of up to 90% in some regions ( Figure S2), so that they might amplify the modification of other climate relevant processes like cloud formation. Surprisingly, the largest changes occur over open water in the Barents and Kara seas with a stabilizing effect (i.e., less upward flux) for all levels of the new parameterizations. This indicates a nonlocal effect linked to circulation changes. In the Greenland Sea, both stabilizing and destabilizing effects (larger upward fluxes) are seen for P1 and P2. The decrease in the transfer coefficients for momentum C m leads to a significantly decreased friction velocity u Ã as a representative of the turbulent momentum flux (of up to 15%) over most of the sea ice areas (Figure 2). Interestingly, the changes in u Ã are of the same order of magnitude for all parameterization levels even though C m is significantly lower for P1 and P2 than for P3. This can be explained by Equations (1) and (2), showing that the effect of changes in the transfer coefficient and wind speed can cancel each other. This is supported by the modeled increase in wind speed (Figures 2 and S2) even for areas where the friction velocity decreases. The F I G U R E 2 December 2007 ensemble mean difference fields of the friction velocity over sea ice (in mÁs À1 ; top), the surface sensible heat flux (in WÁm À2 ; upward/downward fluxes are negative/positive; middle), and the wind speed (in mÁs À1 ; bottom) at the lowest model level (height of ca. 9-10 m) for the differences P1 À P0, P2 À P0, and P3 À P0 with p-values ≤0.05 in black shading. The marginal sea ice zone (15% < ice concentration < 80%) is indicated by the dotted and dashed pink lines. The associated relative changes (relative to P0) are shown in Figure S2 changes in wind speed suggest that circulation changes play a role. Figure 3 shows significant changes in mean sea level pressure (MSLP) for P2, namely a decrease over the Barents and Kara seas, where heat fluxes also changed (see above). This result indicates an amplification and eastward extension of the northern North Atlantic storm track into this region, which also agrees with the significantly increased transport of warm moist air from the northern North Atlantic into the Kara Sea region (Schneider, 2020). The circulation changes are quasibarotropic reaching from the surface to the mid-troposphere. As a representative for the latter, the changes in the 500 hPa geopotential height (Z500) are also shown in Figure 3 revealing a similar spatial pattern as the MSLP changes. These results fit with Renfrew et al. (2019), who found that reducing roughness over sea ice leads to lower pressure over the sea ice, increase of wind speed and reduced bias of their model (see their Figure 8 and description). The combined effect of both stability and roughness changes (P1) in Figure 3 shows weaker and statistically insignificant changes in the atmospheric circulation, indicating that the considered changes in the SCFs (P3) mostly counteracts the influence of the changes in roughness in a nonlinear way. In general, the modifications do not operate in a linear way, neither for the fluxes nor for the circulation response. (P2 À P0) + (P3 À P0) does not simply sum up and can be even reverse to (P1 À P0) ( Figure S3).
Compared to the driving ERA-Interim reanalysis, the standard HIRHAM5 run (P0) has biases in MSLP of about ±4 hPa ( Figure 4) and in Z500 of ±50 m ( Figure S4). Significantly higher MSLP and Z500 are seen over the Atlantic side of the Arctic Ocean. Contrary, MSLP and Z500 are significantly underestimated over the Pacific side and the central Arctic Ocean. From all the parameterization changes, only P2 tends to reduce MSLP and Z500 biases mainly over the Beaufort and Barents Seas. The same differences show up compared to the ERA5 reanalysis ( Figure S5).

| ABL process diagnostics
One approach to investigate changes in the state of the wintertime ABL is the study of the two atmospheric states, the so-called opaquely cloudy and radiatively clear states, associated with a weakly stable and strongly stable ABL, respectively (Stramler et al., 2011;Morrison et al., 2012). These states are apparent in the Probability Density Functions (PDFs) of the net longwave radiation (LWnet) at the surface ( Figure 5). Based on observations F I G U R E 3 December 2007 ensemble mean difference fields of the mean sea level pressure (in hPa; top) and the 500 hPa geopotential height (in m; bottom) for the differences P1 À P0, P2 À P0, and P3 À P0 with p-values ≤0.05 in black shading. The marginal sea ice zone (15% < ice concentration < 80%) is indicated by the dotted and dashed pink lines from the North Pole-35 sea ice drift station (NP-35) ( Figure S6) in 2007, the year of this study, the bimodal distribution of LWnet shows two distinctive peaks with a cloudy-state peak close to 0 WÁm À2 and a radiatively clear-state peak around À40 WÁm À2 . This distinct feature fits with observations from SHEBA, N-ICE2015 and other North Pole drift stations (Graham et al., 2017;Makhotina et al., 2019). Figure 5 also demonstrates the limited skill of this aspect in ERA-Interim and ERA5.
HIRHAM5 shows a well-positioned bimodal distribution for all parameterization levels and higher skill than other climate models for the reproduction of the two states (Pithan et al., 2014). HIRHAM5 even outperforms the two reanalyses in this regard. However, compared to NP-35, the ensemble simulations show a lower probability for the opaquely cloudy state. This agrees with the model evaluation by Sedlar et al. (2020), who showed additionally that the radiatively clear conditions in HIRHAM5 are associated with only modest lower tropospheric stability (difference in equivalent potential temperature between the 950-hPa level and the surface of 0-10 K). In contrast, the observed radiatively clear F I G U R E 4 December 2007 ensemble mean difference fields of the mean sea level pressure (in hPa) compared to ERA-Interim (ERAI) for the differences P0-ERAI, P1-ERAI, P2-ERAI, and P3-ERAI with pvalues ≤0.05 in black shading. The marginal sea ice zone (15% < ice concentration < 80%) is indicated by the dotted and dashed pink lines F I G U R E 5 Probability density function of net longwave radiation at the surface (in WÁm À2 ), December 2007, based on 6-hourly data. Data are shown from observations of the North Pole drifting station NP-35, from ERA5 and ERA-Interim reanalyses, and from the four simulations P0 À P3. For the latter, the solid lines show the complete ensemble, while the shading shows the individual ensemble member range. For the reanalyses and simulations only grid points within the NP-35 region (81-87 N, 45-115 E; Figure S6) and with sea ice concentration >95% are included conditions are associated with stronger lower tropospheric stability of 10-20 K, based on the Arctic Clouds in Summer Experiment (ACSE) 2014 campaign data. A possible reason might be related to an over/ underestimation of decoupling processes between surface and atmosphere in the model, which should be considered in future research.
In the NP-35 region, the considered modifications of the surface-layer turbulence parameterization for stable conditions do not significantly affect the representation of the two ABL winter states in the simulations. The PDF differences are within the ensemble member spread although there are hints that P2 and P0 simulate more cloudy states and P1 and P3 more clear-sky states. In accordance, in the monthly mean spatial distribution, there are several regions with changes in cloud cover and net longwave radiation in the range of 5-10% relative to P0 (not shown). Thus, roughness and stability can have significant local impacts on both quantities. Also, on short time scales, changes in roughness and transfer coefficients clearly affect boundary layer characteristics, especially in case of strong stability. In the model runs, these situations apparently do not occur often enough with respect to time and space to significantly alter the modeled winter states. Alternatively, events with positive and negative deviations from the mean state are canceling each other.

| SUMMARY AND CONCLUSIONS
Based on this case study, the modifications of the surface-layer turbulence parameterization proposed by GL18 lead to significant changes in simulated mean surface-layer turbulent fluxes and atmospheric circulation from the surface to the mid-troposphere over the entire Arctic Ocean. This study demonstrates the significant impact of changes in both the surface roughness of sea ice and the modified stability correction function (SCF) for stable conditions. By stepwise implementation of the GL18 features, we found that an overall decrease of the surface roughness to values that agree better with values from the SHEBA campaign had a much stronger effect on the monthly averaged fields than the modification of SCFs. This is, however, not surprising because the model simulated slightly unstable conditions and nearneutral stability (Ri b < 0.02) over sea ice at about 76% of all sea ice grid points, so that the stability impact in the model is relatively small as compared to measurements during SHEBA, where Ri b ≥ 0.02 occurs more often. A possible reason for this drawback might be the underrepresentation of stability, but also the simplifying assumption of a constant roughness all over the huge Arctic sea ice region. One might speculate that the inclusion of variable roughness accounting for the form drag caused by sea ice ridges and floe edges might lead to an improvement also when the physically grounded changes in the SCFs are applied. Despite the positive impact of the decreased roughness, the constant value might be problematic in regions far from the SHEBA station location in the Beaufort Sea. Future work combining the modified change in SCFs and allowing more variability in sea ice roughness, for example, by the scheme proposed by Lüpkes and Gryanik (2015) or with parameter values proposed by Elvidge et al. (2016) should be a next step.
Changes in the turbulent fluxes affect the regional circulation and thus cloud cover and radiation. However, the wintertime ABL state, represented by the two opaquely cloudy and radiatively clear states, is not significantly affected by any of the surface-layer turbulence parameterization levels when analyzed specifically over the NP-35 region. The changes are within the ensemble member spread. This might indicate that improvements in the model's cloud representation (e.g., cloud liquid), which interact with boundary layer turbulence also in the lowest levels, might lead to a more significant improvement. However, the next necessary steps towards further improved surface-layer turbulence parameterization with an expansion of the parameterization to a scheme with spatially varying roughness length scales and with a modified stability dependence not only in the surface layer but also in the whole ABL might lead to more significant impacts.
We stress also that a change of the transfer coefficients alone was not sufficient to overcome the model's inability to reproduce the stability as strong as observed during SHEBA. Other reasons than the turbulence parameterization must play a role. A possible candidate is the snow cover on sea ice and its isolating effect (e.g., Batrak and Müller, 2019), which might be inadequately represented in HIRHAM, but affects surface temperature and eventually the stability regimes. This is another direction future work should focus on.
Research Center "ArctiC Amplification: Climate Relevant Atmospheric and SurfaCe Processes, and Feedback Mechanisms (AC) 3 ". SHEBA observations were provided by NCAR/EOL under the sponsorship of the U.S. National Science Foundation (NSF). ERA-Interim and ERA5 reanalyses were provided by ECMWF. We are grateful to Dimitry Sidorenko for useful discussions. We thank Ian Brooks and two anonymous reviewers for their helpful comments and suggestions.
All individual authors listed in this manuscript are qualified to be authors and have given their permission to submit this manuscript to Atmospheric Science Letters. The work contained within the manuscript is original and has not been submitted to any other peer-reviewed journals. The manuscript does not include any previously published material that would necessitate reproduction permission.
AUTHOR CONTRIBUTIONS Thea Schneider: Formal analysis; investigation; software; visualization; writing original draft; writing review and editing. Christof Lüpkes: Conceptualization; funding acquisition; investigation; methodology; writing original draft; writing review and editing. Wolfgang Dorn: Conceptualization; data curation; investigation; methodology; software; writing review and editing. Dmitry Chechin: Data curation; resources; writing review and editing. Doerthe Handorf: Conceptualization; funding acquisition; methodology; writing review and editing. Sara Khosravi: Writing review and editing. Vladimir Gryanik: Methodology; writing review and editing. Irina Makhotina: Data curation; resources; writing review and editing. Annette Rinke: Conceptualization; funding acquisition; methodology; writing original draft; writing review and editing.