A dynamic extension of the pragmatic blending scheme for scale‐dependent sub‐grid mixing

A recent pragmatic blending approach treats sub‐grid turbulent mixing using a weighted average of a 1D mesoscale model and a 3D Smagorinsky formulation. Here the approach is modified and extended to incorporate a scale‐dependent dynamic Smagorinsky scheme instead of a static Smagorinsky scheme. Results from simulating an evolving convective boundary layer show that the new scheme is able to improve the representation of turbulence statistics and potential temperature profiles at grey‐zone resolutions during the transition from the shallow morning to the deep afternoon boundary layer. This is achieved mainly because the new scheme enables and controls an improved spin‐up of resolved turbulence. The dynamic blending scheme is shown to be more adaptive to the evolving flow and somewhat less sensitive to the blending parameters. The new approach appears to offer a more robust and more flexible formulation of blending and the results strongly encourage further assessment and development.


INTRODUCTION
Operational numerical weather prediction (NWP) is now moving towards sub-kilometric horizontal grid spacings as more computer power becomes available from massively parallelized high-performance computing systems. For example, the UK Met Office is running its Unified Model (UM) at 333 m horizontal resolution over the broader London area and Heathrow airport (Boutle et al., 2015), providing operational forecasts especially for aviation purposes. Moreover, Hanley et al. (2015) and Stein et al. (2015) have simulated a number of convective days in the UK using the UM for a range of grid spacings between 100 m and 1.5 km and compared model results against radar data. The comparisons revealed marked impacts of sub-grid turbulent mixing on cloud representation and morphology. Although enhanced resolutions have exhibited some significant benefits, such as for fog prediction in complex terrain (Boutle et al., 2015), there are significant challenges in the "uncharted waters" of sub-kilometric meteorological modelling (e.g. Wyngaard, 2004;Honnert et al., 2011;Efstathiou and Beare, 2015;Ito et al., 2017).
The traditional approach in NWP is to consider the entire spectrum of turbulent motions to be unresolved, with the vertical turbulent transfers requiring parametrization. In contrast, when the numerical grid is sufficiently fine that models are able to resolve the most energetic turbulent eddies, then smaller-scale motions can be parametrized with a local 3D sub-grid diffusion scheme, designed to remove energy from the resolved scales. This approach is common in performing large-eddy simulations (LESs), and such simulations have proved valuable for our understanding of boundary-layer (BL) turbulence and for developing BL parametrizations suitable for NWP. However, if the dominant turbulence length-scales are similar to the grid scale, then models will partially resolve the dominant modes of BL turbulence. These resolutions comprise the terra incognita or the grey zone of BL turbu-lence (Wyngaard, 2004) where the fundamental assumptions behind the traditional parametrizations are no longer valid (Wyngaard, 2004;Honnert et al., 2011;Beare, 2014;Efstathiou and Beare, 2015). Boutle et al. (2014) proposed a new BL parametrization for use across the entire range of resolutions, which blends between the standard Smagorinsky 3D formulation (Lilly, 1967) and the 1D NWP scheme of Lock et al. (2000). A transition function was used to weigh between the two schemes depending upon the resolution and boundary-layer regime, following the work of Honnert et al. (2011). This pragmatic blending (PGB) approach was shown to be effective in simulating a stratocumulus deck across grey-zone resolutions (Boutle et al., 2014) and is now used operationally by the Met Office. Efstathiou et al. (2016) tested a simple version of the PGB scheme in the Met Office large-eddy model (LEM) and compared it with the bounding approach (Efstathiou and Beare, 2015) for idealized simulations of an evolving convective BL (CBL) in the grey zone. In the LEM framework, PGB was able to reproduce firstand second-order quantities at both the LES and mesoscale limits. However, in the transition regime there was a notable deficit of resolved turbulence which adversely affected the representation of potential temperature profiles.
In a similar way, Shin and Hong (2015) used a partitioning function to make the relative contribution of local and non-local turbulent transport terms scale-dependent in a non-local BL parametrization. Their results showed improvements over the conventional non-local scheme for the representation of mean first-order and vertical transport profiles. Honnert et al. (2016) used similarity relationships proposed by Honnert et al. (2011) to modify a mass-flux scheme for use at grey-zone resolutions. Moreover, Zhang et al. (2018) tried to extend the use of a 3D turbulent kinetic energy (TKE) model across the grey zone by blending the master mixing and dissipation length-scales from the LES to the mesoscale limit, using the transition function of Boutle et al. (2014) (see also the work of Ito et al., 2015).
In this contribution, we extend the PGB approach by testing a blended formulation in which we replace the standard static Smagorinsky scheme (SMAG) with the Lagrangian-averaged scale-dependent dynamic Smagorinsky scheme (LASD) of Bou-Zeid et al. (2005). Efstathiou et al. (2018) modified and implemented the LASD scheme in the LEM and used it to simulate an evolving CBL at near-grey-zone resolutions at and beyond the LES regime. The dynamic model was shown to significantly improve the initial BL development due to the faster spin-up of appropriate resolved motions. Nonetheless, the LASD model exhibited a usability limit when applied at coarser grid spacings further into the grey zone (Efstathiou et al., 2018).
Here we demonstrate that the combination of a blending approach with the LASD model is an attractive one, seemingly able to combine their benefits while offseting their defects. LASD is blended with a non-local 1D BL scheme in the LEM and tested for an evolving CBL similarly to Efstathiou et al. (2016;. In particular, the ability of the new blending scheme to provide a smooth transition of the first-and second-order quantities across the grey zone is examined in comparison with the standard blending scheme and with coarse-grained LES results. The study therefore provides a proof of concept. It is reasonable to suppose that the full benefits of the method might not be immediately obvious in these idealized, dry CBL LEM simulations. Surface heterogeneity, complex topography and the presence of clouds are some of the situations where the inclusion of a dynamic element in the sub-grid scheme could improve the representation of turbulent mixing (e.g. Kirkpatrick et al., 2006;Huang et al., 2008) across a wide range of scales.

Pragmatic blending implementation
The implementation of the PGB scheme as outlined here follows Efstathiou et al. (2016) which may be consulted for full details. It is based on the approach of Boutle et al. (2014). Specifically, SMAG is blended here with the Yonsei University (YSU) BL scheme of Hong et al. (2006) for the sub-grid sensible heat flux while momentum fluxes (u ′ i u ′ j ) are treated locally: . (1) The sub-grid sensible heat flux (u ′ j ′ ) includes a non-local contribution given by: where w ′ ′ z h denotes the heat flux at the BL height z h and is the counter-gradient term (calculation of these quantities is given in Hong et al., 2006;Efstathiou et al., 2016). K M and K H are the eddy diffusivities for heat and momentum respectively and they will be specified below. Throughout this paper an overbar denotes the averaging that produces the model's grid-scale quantities, while primes denote deviations from the mean. The blending function (W 1D ) which appears in Equation 2 is taken from Boutle et al. (2014) as with the parameter b having a default value of b = 0.15. If a simulation is in the mesoscale limit (z h ∕Δx ≥ 4, W 1D = 1) Equation 2 results in a non-local 1D model. As a simulation tends to the LES limit (z h ∕Δx → ∞, W 1D → 0) then the non-local flux contributions (counter-gradient and entrainment heat flux) to Equation 2 vanish because of the W 1D weighting and Equation 2 tends to a local down-gradient scheme where the non-local fluxes are expected to be explicitly resolved.
The eddy diffusivities for momentum and heat are determined by comparing the standard Smagorinsky model diffusion with a weighted 1D K diffusion profile (K M,H(1D) ) according to (see also Lock et al., 2000): ] . (4) The K M,H(1D) profiles are chosen as in Efstathiou et al. (2016), S is the modulus of the strain rate tensor S ij and f M,H are stability functions for momentum and heat (Brown et al., 1994) which depend upon the gradient Richardson number (Ri) given by: with 0 expressing a reference temperature and g the acceleration due to gravity. Eddy diffusivities are treated in the same way in the vertical and horizontal, according to Equation 4. The blending mixing length (l blend ) is formulated following Boutle et al. (2014), Mixing lengths follow the functional form of Mason and Thomson (1992) with the Smagorinsky mixing length l SMAG based on Brown et al. (1994) as: where is the von Kármán constant and z is the height variable. l 1D represents a "mesoscale" 1D mixing length, The asymptotic values of l 1D account for the weak mixing in shallow, unresolved BLs with a minimum sub-grid mixing length scale of 40 m (Boutle et al., 2014). The Smagorinsky mixing is formulated using its default representation in both the LEM and UM, with Smagorinksy coefficient C S = 0.23.

Dynamic blending scheme
The concept of a dynamic Smagorinsky scheme is to compute the Smagorinsky coefficient as a function of the evolving flow, rather than prescribing it in advance. This is achieved by comparing the flow properties at the grid scale and after filtering onto a test scale (Germano et al., 1991). The concept can be extended to account for scale dependence of the Smagorinsky coefficient by considering a second test scale. The dynamic blending scheme (DNB) considered here incorporates the dynamic Smagorinsky model of Bou-Zeid et al. (2005) as implemented for the LEM by Efstathiou et al. (2018). The dynamic blending uses Equations 1-4 as above, but Equation 6 is replaced by while the dynamic model mixing length (l LASD ) is given by It may be noted that the dynamic calculation of the Smagorinsky coefficient is capable of capturing the effect of the solid boundary on turbulent mixing and so obviates the need to impose a wall-damping function (compare Equations 7 and 10). As in the PGB implementation, eddy diffusivities do not differ between the horizontal and the vertical in DNB. Full details of the procedure for the calculation of C SΔ can be found in Efstathiou et al. (2018).

SIMULATIONS
The PGB and DNB approaches were used in the LEM to simulate the evolving Wangara day 33 convective boundary layer (Clarke et al., 1971) following the set-up of Zhou et al. (2014), similar to Efstathiou et al. (2016;. 3D dry LEM simulations are presented for two different horizontal resolutions of Δx = 400 and 800 m which are considered representative of the grey zone (also Efstathiou and Beare, 2015;Efstathiou et al., 2016). We shall show later that both the 1D and 3D treatments of mixing play important roles in the behaviour of these simulations (e.g. Figure 2 below). Efstathiou et al. (2018) showed that at Δx = 200 m, the stand-alone LASD scheme can well reproduce the filtered LES fields, while at 400 m it does not produce sufficient sub-grid mixing, especially when the BL is shallow. Moreover, for the Wangara case-study, simulations at horizontal resolutions of more than 800 m would be dominated by the 1D non-local scheme. The vertical resolution was set to Δz = 40 m, which is similar to that used within the BL in operational high-resolution NWP models. The domain size was 48 × 48 grid points in the horizontal. The model top was set to 2,500 m, with a damping layer above 2,000 m to avoid any reflection of gravity waves from model lid. Results from the grey-zone simulations are compared with LES data for the same case, which were produced using the standard static Smagorinsky scheme with Δx = 25 m and Δz = 10 m for a 384 × 384 point horizontal domain.
The LES fields were coarse-grained at 400 and 800 m in order to derive the reference fields for the comparison of fluxes at the different resolutions, following the procedure of Honnert et al. (2011) and Honnert and Masson (2014). The resolved (w ′ ′ res ) fluxes at the filtered (coarse-grained) scale Δ f are given by: where w and are resolved variables, hence w ≡ w and ≡ . The average over the filter scale ( Δ f ) expresses box averaging of the reference LES fields over Δ f following Honnert et al. (2011). The angle brackets denote horizontal averaging over the domain. The LEM simulations were integrated using the Boussinesq approximation with the Piacsek-Williams centred-difference advection scheme (Piacsek and Williams, 1970) for the momentum and the Total Variation Diminishing scheme  (Leonard et al., 1993) for the perturbation potential temperature equations. Simulations started from 0900 LST (local time) and ran until 1800 LST. Figure 1 shows horizontally-averaged potential temperature profiles for several times in the 400 and 800 m simulations, along with the corresponding LES profiles. At 1000 LST, when the BL is still shallow, both PGB and DNB are able to reproduce the LES profile. As shown in Figure 2, the weighting function (W 1D ) is close to 1 at this time so that both blending schemes behave like a non-local 1D scheme (see Equation 2). As the BL starts to deepen, some differences between the PGB and DNB simulations emerge, despite the fact that the diagnosed z h are very similar so that the W 1D values in the 400 m simulations are almost identical ( Figure 2). For the 400 m simulation (Figure 1a), PGB produces a superadiabatic temperature profile throughout the BL at 1140 LST whereas DNB and LES are well-mixed above the surface layer. The PGB profile becomes well-mixed an hour later (1240 LST). The same issue is observed in the 800 m simulations (Figure 1b), although in this case the superadiabatic profile produced by PGB is evident at 1240 LST and has become well-mixed by 1340 LST. The behaviour of the PGB temperature profiles indicates insufficient non-local mixing during the PGB simulations. One notable difference between the LES and grey-zone simulations is the stronger inversion strength in the latter, especially at earlier times (1140 LST) during the BL development. One should recall that our LES has four times finer vertical resolution than the grey-zone simulations, so some discrepancies between them in the surface layer and especially in the inversion layer are expected. The delayed onset of resolved turbulence with the PGB scheme is obvious in Figures 3 and 4 where the horizontally averaged non-dimensional resolved vertical velocity variance (w ′2 res ∕w 2 * ) and heat flux (w ′ ′ res ∕w ′ ′ 0 ) are shown together with the corresponding quantities computed from the coarse-grained LES fields following Equation 11. Results are presented for two times for each resolution, before and after the formation of well-mixed potential temperature profiles. For the 400 m simulation with the PGB scheme, there is almost no resolved motion ( Figure 3a) and therefore very little resolved heat flux (Figure 3c) at 1140 LST. This is despite the fact that the mixing from the 1D non-local mixing has been down-weighted such that the parametrized part of the flux is insufficient to produce a well-mixed profile. By the same time, the DNB scheme has already enabled considerable resolved upward motion, matching the intensity of the coarse-grained LES fields. The onset of resolved w ′2 res by 1240 LST in the PGB simulation (Figure 3b) gives rise to non-local heat flux (Figure 3d) which does result in well-mixed potential temperature profiles (Figure 1a) and reasonable agreement in resolved turbulence statistics between the grey-zone simulations and the coarse-grained LES fields. In the 800 m PGB simulation, there is little resolved turbulence at 1240 LST (Figure 4a,c)  slightly superadiabatic profiles, even though most of the turbulent heat transfer is being done through the 1D scheme (Figure 2). At the same time, the DNB run produces substantially stronger resolved fluxes (especially heat fluxes) than the coarse-grained fields. Nevertheless, as the BL evolves further, both of the 800 m grey-zone simulations are able to produce resolved motions similar to the coarse-grained LES (Figure 4b,d). The DNB simulations at both 400 and 800 m remain slightly more energetic than PGB (Figures 3b,d and  4b,d). In additional runs with a 600 m grid length, we found that the benefits of the DNB approach are still evident, with the results being intermediate between those of the 400 and 800 m runs in terms of the impact of the dynamic Smagorinsky scheme (not shown).

Results
To further explore the effects of sub-grid mixing on the resolved turbulence, Figure 5a presents the time evolution of the horizontally averaged vertical profiles of heat eddy diffusivity (K H ) for the 400 m PGB and DNB simulations, while Figure 5b,c compare the two components of the heat eddy diffusivity according to Equation 4 for two different times. The eddy diffusivity profiles are also averaged in time over a 0.5 h period before the output time, to capture the effects of sub-grid diffusion in the temperature and flux profiles, prior to the spin-up of resolved turbulence. The behaviour of K H at 400 m is similar to that in the 800 m simulations; however the impact of the blending approach is more pronounced at 400 m (cf. W 1D evolution in Figure 2) as is evident in the response of the potential temperature profiles to sub-grid diffusion in Figure 1. In accordance with Figures 1 and 2, K H is almost identical for both PGB and DNB runs in the morning (1000 LST). However, near the time of resolved turbulence onset (1140 LST), DNB produces significantly less sub-grid mixing than PGB especially in the lower part of the BL (albeit with slightly more mixing in the upper BL). Examining Figure 5b, it becomes obvious that this is due to the reduced mixing from the LASD in the DNB simulation within the lower part of the BL (as 1D mixing is almost identical between the two runs). This results in the faster onset of resolved turbulence in DNB. The standard Smagorinsky scheme has been shown to be excessively diffusive in the grey zone (Efstathiou and Beare, 2015;de Roode et al., 2017). When the BL deepens (1240 LST), K H profiles from both schemes are similar with DNB remaining slightly more diffusive in agreement with its more energetic fluxes  Figure 5c is dominated by the Smagorinsky diffusion in the lower part of the BL and by the W 1D K H(1D) contribution in the upper part for both PGB and DNB schemes. The strong impact of the different Smagorinsky formulations in the resolved turbulence field demonstrates the importance of the lower BL representation in simulating the CBL at grey-zone resolutions (see also Zhou et al., 2017).

Sensitivity to the blending parameter
The sensitivity of the PGB and DNB schemes to the blending parameter b was investigated by performing additional simulations at Δx = 400 m, setting b = 0.10 and 0.20 as opposed to the default value of 0.15. Time series of the weighting function are shown in Figure 2 and results for l blend , K H and the resolved vertical velocity variances are shown in Figure 6. In the PGB simulations, a change to b does not change the imposed l SMAG function (Equation 7) with its wall function and fixed C S value. l 1D also remains fairly constant between simulations as z h does not change significantly (Equation 8). Thus, changing the b parameter has a tightly-constrained impact on l blend (Figure 6c), with relatively modest impacts in the surface layer. The vertical profile of the K H is shown in In contrast, the DNB scheme has much more scope in responding to a changed b by adjusting the l blend through the dynamically derived C SΔ . As b increases (decreases) the weighting function (Equation 3) allows for more (less) resolved turbulence (Figure 2) and the LASD scheme is found to adapt to the resolved flow field in such a way as to increase (decrease) l blend (Figure 6c). The impact of l LASD can be seen in the DNB eddy diffusivity profiles (Figure 6d) where K H values are almost doubled in the surface layer when b increases from 0.10 to 0.20, to account for the changes in the resolved flow field. The outcome is a more moderate increase of w ′2 res than in the PGB run. For both the PGB and DNB simulations, above 0.3z h , K H is controlled by the 1D scheme and its magnitude is determined by b through the blending function (Equation 3).
The adjustment of the blending mixing length to changes in the blending function is a manifestation of the dynamic interaction of the resolved fields with the adaptive sub-grid model which can in turn modify the partitioning between resolved and sub-grid fluxes. It is notable that in the upper part of the BL the sub-grid mixing length scales in the DNB simulations reduce in the proximity of the inversion layer and tend to zero in the stable free troposphere, while the mixing length in PGB remains almost constant with height above ∼ 0.3z h . However, this does not have any significant effect on the K H profiles as these are dominated by the 1D scheme (Equation 6). The outcome of these dynamic behaviours is that the resolved turbulence is somewhat less sensitive to the weighting parameter b in the DNB (Figure 6b) than in the PGB simulations (Figure 6a) especially for the higher b value. It should be remarked that increasing b does not improve the problems with the spin-up of appropriate resolved motions in the PGB simulation (not shown).

DISCUSSION AND CONCLUDING REMARKS
A new blending approach (DNB) has been demonstrated, based on the pragmatic blending (PGB) of Boutle et al. (2014) but using a dynamic rather than a static Smagorinsky model. There was little difference found in the simulations when the BL was shallow or when the resolution was too coarse, since the blending function then ensures that heat transfer is handled by the 1D non-local scheme. However, the dynamic approach improves the representation of CBL potential temperature profiles and turbulence statistics during the handover from a predominantly 1D non-local scheme towards a more LES formulation of diffusion. This is due to an earlier onset of an appropriate level of resolved turbulence in the DNB grey-zone simulations in comparison with the PGB scheme. The benefits for spin-up of a dynamic Smagorinsky model were shown in Efstathiou et al. (2018). A late onset of resolved TKE in PGB runs and adverse effects on potential temperature profiles has also been found in Efstathiou et al. (2016). Even though the stand-alone LASD scheme exhibits a usability limit when the BL is very shallow compared to the grid length (Efstathiou et al., 2018), the blending of LASD with a non-local mesoscale scheme ensures a seamless representation of sub-grid mixing across the scales. This is due to the fact that diffusion is handled by a more appropriate non-local scheme when Δx∕z h is in the mesoscale regime. In any case, we have shown evidence that the dynamic blending approach has more adaptability, and that this is useful partly in terms of the improved spin-up and partly also in terms of a reduced sensitivity to the parameter controlling the blending function. The dynamic approach alleviates the need for a specified functional form of the Smagorinsky mixing length or for well-chosen C S values. The use of the LASD can add up to a factor of 2 in computational time relative to Smagorinsky scheme calculations. However, we note that the Smagorinsky calculations are a relatively small expense within the context of a full grey-zone NWP model.
The present study is a proof of concept for the practical applicability and significance of the dynamic blending approach. The approach shows promising results for the simulation of an evolving dry CBL. This motivates further work on applying the new method in more complex and realistic simulations. A scale-dependent and flow-adaptive mixing length has the potential to improve the representation of turbulent transfers where surface variability is pronounced (e.g. in the urban environment), in strong stratification (e.g. the nocturnal BL; Basu et al., 2008) or where turbulent length-scales might vary between the BL and the cloud layer . A caveat is that the implementation of the blending procedure presented here is not the same as in the operational Met Office Unified Model where a more sophisticated treatment of the inversion layer is used. For these reasons, an implementation in the Met Office UM would provide an ideal test bed for exploring the potential benefits of dynamic blending in simulating a wide range of atmospheric phenomena.