Systematic errors in northern Eurasian short-term weather forecasts induced by atmospheric boundary layer thickness

This study introduces and explores a new perspective on SAT error analysis in the context of ABL modeling. This perspective emerges from an underappreciated factor, which scales forcing efficacy by the ABL thickness, $h.$ Let us look at the following simple energy-balance equation (e.g. Park et al 2005, Zilitinkevich and Esau 2009, Rypdal 2012, Davy et al 2017):

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}({\rm{\Delta }}T)={c}_{s}{h}^{-1}{\rm{\Delta }}Q.\end{eqnarray} \tag{ 1 }$$

Usually, for example in Rypdal (2012), ${\rm{\Delta }}T$ is understood to be the mean temperature of a bulk well-mixed layer of thickness $h$ and constant heat capacity $C.$ However, it is known that in a well-mixed layer, the temperatures at any layer, and also the SAT, are functionally dependent. This is a trivial fact for a CBL where $\tfrac{\partial T}{\partial z}=-\gamma \equiv -{\gamma }_{a}$ and ${\gamma }_{a}$ is an adiabatic temperature gradient. This is less trivial in a SBL where more sophisticated relationships apply (Zilitinkevich and Esau 2005). The coefficient ${c}_{s}\lt 1$ should account for the weaker mixing in a SBL where $\gamma \lt {\gamma }_{a}.$ Following the accepted methodology of turbulence studies, we understand equation (1) in a statistical sense as valid only on average. Here, ${\rm{\Delta }}T={T}_{m}-{T}_{o}$ is the SAT difference (error) in the model, i.e. the deviation of the model SAT value ${T}_{m}$ from the corresponding collocated observed value ${T}_{o}.$ This deviation is caused by a proper energy imbalance ${\rm{\Delta }}Q$ in the model. In terms of equation (1), ${\rm{\Delta }}Q$ is the surface turbulent sensible heat flux anomaly. In the NWP model, this anomaly could be caused by errors in other components of the heat budget, which modify the sensible heat flux through the surface energy-balance equation.

This SAT response, ${\rm{\Delta }}T,$ is scaled by the effective heat capacity of the system, which could be expressed through $h.$ In the ABL, the energy imbalance ${\rm{\Delta }}Q$ is changing the atmosphere through vertical turbulent exchange. The fact that the turbulent exchange is faster than the advective or radiative changes has been also used in the construction of turbulence schemes in the NWP models. We will consider two such schemes in this study. In an operational NWP model, the SAT differences can be also of a dynamic origin. We acknowledge this difficulty, but do not consider it in this research letter.

In the studied regions of northern Russia, most, if not all, land is covered with snow in January. The snow acts as insulator inhibiting heat exchange between the atmosphere and the soil. The data used in this study show that the observed daily maximum temperature slightly exceeds 0 °C at some stations on some days in the north European part of Russia; however, the daily minimum temperature is always below 0 °C except for 1–3 of January for three stations. In Siberia, there are only few occasional positive daily maximum temperatures, all daily minimum temperatures being negative. This means that soil and, in particular, soil wetness effects do not play any role in the SAT differences.

The ABL thickness, $h,$ is a scaling factor moderating the temperature response, ${\rm{\Delta }}T,$ on the proper energy imbalance, ${\rm{\Delta }}Q.$ Cuxart et al (2006) and Beare et al (2006) demonstrated that $h$ could be as much as an order of magnitude larger in the turbulence schemes than in the corresponding observations or in turbulence-resolving large-eddy simulations. A small scaling factor makes SAT sensitive even to small energy imbalances. Contrary, a large $h$ damps the SAT sensitivity even to large imbalances. The scaling factor does not determine the sign of the temperature response. Davy and Esau (2014) found that the majority of climate models maintain ABLs that are too thick, particularly over continents and at high latitudes, where the observed ABLs are naturally shallow due to the mostly negative surface heat balance (Serreze et al 2007) and stronger stability of the lower atmosphere (Medeiros et al 2011, Zhang et al 2011).

It needs to be clarified that the scaling factor $h$ in equation (1) should not be considered as yet another underrepresented component of the surface energy budget. This is a factor which scales ${\rm{\Delta }}Q$ to some degree independently of its nature. We acknowledge that some errors in the heat budget would be correlated with $h.$ For instance, radiative effects of the convective cloud cover are likely to be correlated with large $h$ in the thick convective boundary layer. Detailed investigations of such correlations require a separate study, being outside of scopes of this letter.

Our previous studies of the climate models and ERA interim reanalysis (Esau et al 2012, Davy and Esau 2014, Davy et al 2017) convinced us that $h$ scaling as such could be studied statistically with no regard to the intrinsic physical structure of the forcing. The ABL thickness varies by an order of magnitude from $h\sim O({10}^{2}\,{\rm{m}})$ in a SBL to $h\sim O({10}^{3}\,{\rm{m}})\,\,$in a CBL. Therefore, ${\rm{\Delta }}Q$ that is well correlated with the SBL events will be significantly more amplified. Thus, if the SAT difference is indeed sensitive to $h,$ some kind of a reciprocal relationship ${\rm{\Delta }}T\propto {h}^{-1}$ should surface in the NWP data. To reveal such dependences, we will transform the data from the traditional analysis of ${\rm{\Delta }}T(t)$ to an analysis of ${\rm{\Delta }}T(h).$

It is possible to demonstrate that the differences in $h$ introduce a specific statistical signature. Unfortunately, the ABL thickness is not available from routine ground-based meteorological observations. Its determination from radiosounding data provides rather uncertain values, which are sensitive to the choice of method and the stability of the lower atmosphere (e.g. Seibert et al 2000, Seidel et al 2012, Wang and Wang 2016). Therefore, the observed and modeled $h$ could be compared only in special idealized studies, for example in the GEWEX Atmospheric Boundary Layer Study (GABLS) (Cuxart et al 2006, Svensson et al 2011), where the effects of intermittency and statistics of the cooling/warming events are not considered. The operational weather forecast generates such statistics providing an opportunity to study the SAT errors in ABLs of significantly different thickness.

The scaling factor not only modifies the amplitude of the response, but also imposes low-pass filtering on it (Rypdal 2012). Equation (1) has the following solution in response to a delta-pulse forcing ${\rm{\Delta }}Q\delta (t)$

$$\begin{eqnarray}&&{\rm{\Delta }}T(t)={\rm{\Delta }}Q{c}_{s}{h}^{-1}{{\rm{e}}}^{-t/(\gamma h)}\end{eqnarray} \tag{ 2 }$$

where $\gamma$ is a combined diabatic sensitivity parameter (Park et al 2005). It follows that ${\rm{\Delta }}T$ in the shallow SBL will not only be amplified but will also show a more agile response and will respond to more high-frequency external fluctuations than a thick CBL. Considering realistic chaotic forcing, the low-pass filtering in equation (2) creates a biased memory effect on the SAT. Slower fluctuations are better captured by a thick ABL. Too large an $h$ in the model should induce larger SAT differences during periods with more rapid heating fluctuations. Thus, stable anti-cyclonic conditions with a pronounced diurnal cycle of ${\rm{\Delta }}Q$ should result in a larger ${\rm{\Delta }}T$ than the conditions with dominating synoptic variability. The low-pass filtering induces insufficient diurnal range in the models. The damped modeled temperature range was found in comparisons with satellite temperature products (Liu et al 2006) and with in situ data (Yang et al 2006). Although these kinds of errors are difficult to disentangle, a look at the probability density function (PDF) may provide a useful clue to the scale of the effects (Ashkenazy et al 2008). This effect will be observed even in the situation of a perfectly modeled surface energy budget.

The SAT itself could be used for rough estimation of $h$ (for momentum) as Esau et al (2012) have demonstrated. Indeed, a shallow SBL is created by surface cooling or warm air advection. Those conditions lead to a lower temperature of the surface than of the air above it. In contrast, a deep CBL is created by surface warming or cool air advection. Statistically, the dependence ${\rm{\Delta }}T(h)$ could be transformed to ${\rm{\Delta }}T(T).$ The latter is easier to study in comparison with the observational data. Previous studies, for example the evaluation by Atlaskin and Vihma (2012) of four NWP models, demonstrated that ${\rm{\Delta }}T$ increases with decreasing temperature and strengthening of the low-level temperature inversions. They also noted the low-pass filtering of the forcing fluctuations as the models responded too slowly in situations of rapid drops in temperature. The magnitudes of SAT tendencies in the models were only 17%–20% of the observed magnitudes. A similar slow response was noted in Beesley et al (2000). Fay and Neunhauserer (2006) demonstrated that the slow response and warm bias in their limited area model forecast did not disappear even in very high-resolution (1 km) runs.

The SL-AV global atmospheric general circulation model (Tolstykh 2003, Tolstykh et al 2015) is used for operational medium-range numerical weather prediction, and is applied as a component of the probabilistic long-range forecast system at the Hydrometeorological Centre of the Russian Federation (HMCR). The model is based on a semi-Lagrangian semi-implicit dynamical core, which has been developed and described by Tolstykh et al (2017). The model uses a stretched latitude–longitude grid with a locally refined resolution over Russian territory (Tolstykh et al 2017). The original SL-AV dynamical core uses the vertical component of the absolute vorticity and divergence as prognostic variables. It employs fourth-order finite differences at the unstaggered grid to approximate non-advective equations terms. In this study, we use the operational version of the SL-AV model for medium-range weather forecasts with the resolution in latitude varying from 0.155° to 0.178° in the Northern Hemisphere to 0.24° in the Southern Hemisphere. The resolution in longitude is 0.225°; the model has 51 vertical levels with the highest model level placed at 5 hPa, with 10 levels being located in the lowermost kilometer of the model atmosphere.

Physical processes in the SL-AV at unresolved scales are parameterized mostly using the algorithms developed by the mesoscale weather prediction consortium ALADIN/LACE (Termonia et al 2018). Shortwave and longwave radiation is calculated using CLIRAD SW (Chou and Suarez 1999, Tarasova and Fomin 2007) and RRTMG LW (Mlawer et al 1997) algorithms, respectively. Some parameterizations of important physical processes, for example a multilevel soil model by Volodin and Lykossov (1998), were developed in addition to the standard version.

The SL-AV forecast was run with two different ABL schemes, which we will refer here to as the old scheme and the new scheme. Both schemes provide a diagnostic estimate of the ABL height. The old scheme is also known as the pTKE scheme. We will briefly present the central elements of these schemes following Geleyn et al (1994), Redelsperger et al (2001), Geleyn et al (2006) and Bašták Ďurán et al (2014), hereafter referred to as BD14. A reader seeking for an exhaustive description should refer to those publications.

The old scheme is based on a version of the scheme of Louis (1979) in which advection and self-diffusion of the turbulent kinetic energy (TKE) is computed by the level 2.0 Louis-type scheme taken to be the stationary solution of the system with unresolved fluxes. The latter conditions mean that the TKE is obtained from the turbulent exchange coefficients following the procedure in Redelsperger et al (2001). The ABL thickness, $h,$ is a diagnostic parameter, and is calculated following the approach of Troen and Mahrt (1986) modified by in an iteration scheme by Ayotte et al (1996). The ABL thickness is calculated as

$$\begin{eqnarray}&&h=R{i}_{c}\displaystyle \frac{({U}_{h}^{2}+{V}_{h}^{2})}{g/{\theta }_{s}({\theta }_{h}-{\theta }_{s})},\end{eqnarray} \tag{ 3 }$$

where $R{i}_{c}$ is the critical the bulk Richardson number formulated using the difference between the surface temperature and the temperature at the top of the mixed layer, ${U}_{h},{V}_{h}$ are the components of the horizontal velocity, ${\theta }_{h},{\theta }_{s}$ are the potential temperatures at the top of the ABL and at the surface and $g$ is the acceleration due to gravity. Equation (3) is tried for different model levels starting from the lowermost one until the bulk Richardson number becomes larger than the prescribed $R{i}_{c}.$ It is used to calculate the mixing length scale, ${l}_{m},$ at different elevations above the surface. The corresponding equation is shown below to demonstrate that ${l}_{m},$ and therefore the turbulent fluxes, decay smoothly with the height, $z,$ but are not confined within the layer of thickness $h$:

$$\begin{eqnarray}&&\begin{array}{l}{l}_{m}=kz{\left(1+kz{l}_{m0}\displaystyle \frac{1+{{\rm{e}}}^{\tilde{z}}}{{c}_{m}+{{\rm{e}}}^{\tilde{z}}}\right)}^{-1},\\ \tilde{z}=-{a}_{m}\sqrt{\tfrac{z}{h}}+{b}_{m}.\end{array}\end{eqnarray} \tag{ 4 }$$

Here, $k$ is the von Karman constant and ${l}_{m0},\,{a}_{m},{b}_{m},\,{c}_{m}$ are the tuning constants. Thus, a direct modification of the ABL thickness in this scheme is hardly possible. Moreover, the effects of such modifications will not be straightforward to evaluate as the turbulent fluxes cannot be adjusted independently of $h.$

An important structural feature of turbulent mixing in a SBL as described in Zilitinkevich et al (2013), is the fact that the velocity and temperature fluctuations are differently transferred by the turbulent motions in the stably stratified fluid flow. Stability damps the vertical velocity fluctuations, $w^{\prime} ,$ but it simultaneously increases the temperature fluctuations, $\theta ^{\prime} ,$ as the temperature gradients increase. As a result, the temperature flux, $\overline{w^{\prime} \theta ^{\prime} },$ remains significant even in very stably stratified flows. There are two effective consequences of this feature: the ABL thickness is not directly proportional to the stability of the layer—it are linked through a non-constant turbulent Prandtl number, Pr; and vertical heat mixing remains significant in a shallow SBL but rapidly vanishes above it. This better understanding lead to a new energy-flux balance (EFB) family of parameterization schemes.

A Third Order moments Unified Condensation Accounting and N-dependent Solver for turbulence and diffusion (TOUCANS) scheme was proposed in BD14. A variant known as the eeEFB scheme is referred to as a new scheme in this study. The new scheme is based on a compact model for the stability dependence of the turbulent kinetic energy without a critical Richardson number. It is based on an analytical framework combining and simplifying the revisited Mellor–Yamada formulation (Canuto et al 2008) and the EFB approach (Zilitinkevich and Esau 2007, Zilitinkevich et al 2013). Significant changes were introduced in the formulations of the stability correction functions. An essential feature of the new scheme is the use of the total turbulent energy (TTE) as a prognostic variable. The TTE is a sum of the TKE and the turbulent potential energy, which is given by the squared temperature fluctuations, $\overline{\theta {{\prime} }^{2}}.$ The TOUCANS formulation is valid for all stability conditions, and accounts for turbulence anisotropy in both momentum- and heat-related terms. The TOUCANS parameterization also accounts for the third-order moments for unstable stratification. The mixing length is calculated in the same way as in the old scheme (equation 3).

The SL-AV model was set to compute 27 h global forecasts starting from the HMCR operational analyses for each day of January 2015 at 12:00 UTC. The forecast data for lead times from 6 to 27 h are stored with a time step of 3 h. The forecast data for the first 6 h are dropped to exclude possible spinup effects. Longer lead times are excluded from our consideration as 'slow' diabatic processes (e.g. radiation) play a more important role for these times. The SL-AV model was run with both old and new boundary layer schemes providing two forecast series. All other model settings were identical in both series. The forecast SAT and $h$ were interpolated into locations of selected meteorological stations using the bilinear interpolation method.

The chosen region of northern Russia (figure 1) is a rather homogeneous and flat region. Only a small fraction of it is occupied by the Khibiny and Urals mountains. To maintain homogeneity of our analysis, meteorological observations from elevated stations were not considered. We also exclude several arctic coastal stations as the current model grid does not represent their surface type correctly. We include only those stations that provided eight SAT observations every day between 1 and 31 January 2015 without missing and wrong values. Wrong values are defined as those that differing by more than 12 K from the observed values 3 h before and 3 h later. The whole region was broken in two distinct climatic sub-regions, namely the region of the northern European part of Russia (58°–70° N, 30°–55° E) with a temperate continental climate and more clouds and the western Siberian part (55°–75° N, 60°–85° E) with a severe continental climate and fewer clouds. A description of the prevailing meteorological situation in January 2015 in these sub-regions is presented in the supplementary material (available online at stacks.iop.org/ERL/13/125009/mmedia). In total, we considered 33 meteorological stations in the European part and 40 stations in Siberia.

Zoom In Zoom Out Reset image size

The considered meteorological quantities varied widely during the period of the study. This variability provides us with representative statistics of advection- and radiation-dominated synoptic conditions. The first half of the month experienced a transient synoptic activity with large (up to 30 °C) but slow SAT variations. The surface heat budget changed relatively slowly with a characteristic time scale of several days. The second half of the month experienced more persistent anti-cyclonic conditions with cold temperatures. The changes in the surface heat budget were dominated by the diurnal cycle. The diurnal temperature range was small due to low winter-time solar elevation at these latitudes.

The temperature difference, ${\rm{\Delta }}T,$ between the observed, ${T}_{o},$ and predicted, ${T}_{m},$ SAT also varied considerably. The variations are not chaotic but tend to fell in two bins. One bin is characterized by a shallow ABL (the blue bars in the $h$ panel) and negative ${\rm{\Delta }}T.$ The other bin is characterized by a thick ABL (the red bars) and positive ${\rm{\Delta }}T.$ This suggests that the model under-predicts the higher SAT in the thick ABL and over-predicts the lower SAT in the shallow ABL.

The data also reveal that ${T}_{m}$ fluctuated more slowly than ${T}_{o}.$ This frequently leads to situations when ${T}_{m}$ has no time to fully respond to the imposed cooling before this forcing vanishes. One can note that the new scheme produced a thicker ABL (e.g. from 28 to 30 January). This immediately resulted in an even smaller diurnal temperature range. Previous studies tried to attribute these under- and over-predictions separately, looking at different components of the surface heat budget and considering turbulent fluxes in different stability regimes. We argue that these differences should be attributed to a common factor such as the damped SAT response in too thick an ABL.

To advance the argument, the data for this station are recast to the SAT and ABL thickness parameter space. Figure 3 shows the same data as in figure 2 but presented with respect to the SAT. Now, it is clearly seen that the positive differences are concentrated at higher temperatures, whereas the negative ones are at lower temperatures. Furthermore, the empirical SAT histograms were obtained using temperature bins (10 in this study) which cover the entire range of the SAT variability in January. A certain lack of predictions for cold temperatures can be identified.

Zoom In Zoom Out Reset image size

Contrary to common expectations, the new scheme does not improve the histogram at the low-SAT end, while it does introduce the improvements at the high-SAT end. The slopes of the dependences ${\rm{\Delta }}T({T}_{o})$ on the upper panel are almost the same for the old and new schemes, but the new scheme produces smaller differences at higher temperatures. This behavior is surprising as the new scheme was sought to primarily improve the cold temperature prediction. Thus, the model with either scheme tends to over-predict not because it predicts systematically too high temperatures; in fact, the model under-predicts a high SAT, but because it does not reach the cold and warm temperatures during the cooling or warming events (see supplementary material). The major regular features found in the SAT errors remain in the analysis of all stations in both the European and Siberian regions (figure 4).

Zoom In Zoom Out Reset image size

The observed systematic over-prediction, ${\rm{\Delta }}T\gt 0,$ of the lowest and under-prediction, ${\rm{\Delta }}T\lt 0,$ of the highest temperatures suggests that the model produces SAT tendencies that are too small, ${T}_{t}={\rm{d}}T/{\rm{d}}t.$ As the surface heat budget over longer periods of time should be approximately correct in the model without systematic temperature drift, the damped tendencies indicate too large a scaling factor. Such a relationship between ${T}_{t}$ and $h$ has been inferred from the solution in equation (2). We have also mentioned that this solution would have a low-pass filtering effect. It creates larger SAT differences under the action of higher-frequency forcing. Having enough data, we can statistically quantify this effect on the SAT.

The empirical histogram of ${T}_{t}$ reveals that the larger tendencies are under-predicted in the model runs with both schemes (figure 5). Large positive skewness of 3.6 in the European and 2.8 in the Siberian regions indicates the presence of relatively short periods with strong warming tendencies. It is physically consistent with the advective nature of the winter-time warming events. Very large excess kurtosis (up to 11) points to a significant excess of large ${T}_{t}$ driven by persistent warming or cooling. At the same time, shorter (diurnal) heating does not produce large enough tendencies.

Zoom In Zoom Out Reset image size

The ABL thickness was not identical in the forecast series with the old and new schemes. Although the new scheme does not modify $h$ directly, it enhances the turbulent fluxes, and hence increases $h,$ which is related to the mixing length scale through equation (3). Figure 5 and others show that $h$ in the new scheme is systematically larger than in the old scheme in all temperature and temperature tendency regimes. This gives us an opportunity to test the described effects of the scaling factor in equation (1) evaluating differences between these two forecasts.

The ABL thickness is larger in the forecasts with the new scheme. We divided the differences in the ABL thickness from these two forecasts into 20 bins. The SAT differences (errors), ${\rm{\Delta }}T,$ were calculated for each bin of the $h$ differences. The analysis in figure 6 reveals that the new scheme indeed produces larger $h,$ particularly in the western Siberian part of Russia. The SAT differences in the new scheme are smaller than those in the old scheme in all $h$ bins. One should note that 70% of the cases are found in the bins with the smallest $h$ differences of ±0.1 km where both schemes shows practically identical errors.

Zoom In Zoom Out Reset image size

We further split ${\rm{\Delta }}T$ into five bins defined as $h$ < 0.25 km, 0.25 km < h < 0.5 km, 0.5 km < h < 0.75 km, 0.75 km < h < 1.0 km and h > 1.0 km. The ABL thickness is taken from the forecast with the old scheme. In each bin, the histogram analysis in figure 7 helps us to understand the internal statistical structure of the model errors. The model exhibits a robust positive difference in the runs with both the old and the new schemes. The difference is slightly less positive in the runs with the new scheme. However, the largest improvement caused by the new scheme is observed in the thickest ABL. In the Siberian region, the positive difference reduction by the new scheme is statistically significant for h > 0.75 km.

Zoom In Zoom Out Reset image size