Energetically consistent scale‐adaptive stochastic and deterministic energy backscatter schemes for an atmospheric model

Current climate models still suffer from many biases which are partly due to excessive subgrid‐scale dissipation. Here we systematically develop energetically consistent stochastic energy backscatter (SEB) and deterministic energy backscatter (DEB) parametrization schemes. We implement our schemes in a simplified spectral atmospheric General Circulation Model (GCM). It is shown that the SEB scheme performs better than the DEB scheme at low horizontal resolutions (T21 and T31), whereas the performance of both schemes becomes comparable as the resolution increases to T42 when comparing with our reference simulation at T127 resolution. The energy backscatter parametrization schemes improve eddy variability in low‐resolution models and correctly capture the dominant mode of zonal‐mean zonal wind variability. The autocorrelation time‐scale of low‐resolution models is also found to be more consistent with the reference simulation when applying the SEB and DEB parametrizations. Our schemes are scale‐adaptive and computationally efficient.


INTRODUCTION
With the manifold increase in computing power over the last few decades, it has become possible for weather and climate forecast centres to run numerical prediction models at very high spatio-temporal resolutions. For example, modelling centres have been able to run global forecast models at O(1 km) horizontal resolution leading to considerably improved forecasting accuracy (Hólm et al., 2016). However, running General Circulation Models (GCMs) at such high resolutions requires tremendous computational resources and associated data storage capabilities. Even after such heavy investments, the models still fail to accurately capture many of the important small-scale features (such as cloud dynamics, tropical convection, gravity wave drag, and other microphysical processes). Since these small-scale unresolved (subgrid-scale) processes interact with and influence large-scale resolved processes, their inaccurate representation in the forecast models leads to, among others, under-dispersive ensemble forecasts, and biases in e.g. 500 hPa geopotential height and precipitation (Franzke et al., 2015;Berner et al., 2017;Frederiksen et al., 2017). The use of stochastic modelling helps to overcome some of these problems. Berner et al. (2017) argue that stochastic parametrizations are essential for estimating uncertainty in weather and climate predictions, reducing systematic model errors arising from unrepresented subgrid-scale fluctuations, triggering noise-induced regime transitions, and capturing the response to changes in the external forcing. There have been many studies in recent years in which stochastic parametrization such as the Stochastically Perturbed Parametrization Tendency (SPPT) scheme or a Stochastic Kinetic Energy Backscatter Scheme (SKEBS) has been applied to atmospheric model components in order to improve predictability of numerical forecast models (Weisheimer et al., 2014;Batté and Doblas-Reyes, 2015;Dawson and Palmer, 2015;Sanchez et al., 2016;Berner et al., 2009;2017 and references therein).
The construction of robust and efficient low-resolution models is of practical interest for the modelling community (Zurita-Gotor et al., 2015), especially for understanding long-term climatic processes (e.g. Claussen et al., 2002), including palaeoclimates, or extreme events (e.g. Franzke, 2017). Relatively low-resolution models will also allow the generation of large ensembles for uncertainty estimation. One of the main advantages of stochastic parametrizations is that these can be used to improve low-resolution atmospheric and oceanic models so that they achieve simulation statistics comparable to high-resolution models (Jansen and Held, 2014;Zurita-Gotor et al., 2015;Berner et al., 2017 and references therein). The main aim of this article is to exploit this potential of stochastic parametrization for improving coarser-resolution atmospheric models. For this purpose, we are using the spectral Portable University Model of the Atmosphere (PUMA: Frisius et al., 1998;Franzke et al., 2000;2001;Franzke, 2002;Fraedrich et al., 2005). Here we propose an energetically consistent modified spectral SKEBS parametrization to overcome the problem of over-dissipation at smaller scales and the neglect of kinetic energy backscatter from unresolved subgrid-scale processes.
In addition to the development of stochastic parametrizations, studies have also focused on developing deterministic energy backscatter schemes for improving low-resolution models. Jansen and Held (2014) proposed a novel energetically consistent deterministic (as well as stochastic) backscatter parametrization scheme to compensate for the loss of energy due to hyperdiffusion at smaller scales in an idealized two-layer quasi-geostrophic model. Their scheme combines a hyperviscous closure with a deterministic forcing term and represents the backscatter of kinetic energy from the unresolved subgrid-scales into the resolved flow. The deterministic forcing term is chosen such that it cancels the spurious energy dissipation associated with the hyperviscosity, while maintaining a net dissipation of enstrophy . Later, Jansen et al. (2015) and Zurita-Gotor et al. (2015) showed the advantage of their deterministic backscatter scheme in more realistic ocean and atmospheric models, respectively. Another aim of our study is to use a modified Jansen and Held (2014) deterministic energy backscatter parametrization (described in section 2 below) for improving PUMA at coarse resolutions. The results of the stochastic and deterministic energy backscatter parametrizations will also be compared and contrasted.
We organize the article as follows. Section 2 briefly describes PUMA and the stochastic and deterministic kinetic energy backscatter parametrization schemes. The results of implementing these schemes and their quality in PUMA at different horizontal resolutions are given in section 3 and conclusions are given in section 4.
For the purpose of this article, we only discuss the vorticity equation below. The non-dimensional prognostic equation for vorticity ( ) in PUMA may be written as: where T ′ is the temperature deviation from a constant reference temperature T 0 . U and V are the zonal and meridional wind components multiplied by cos . The horizontal coordinates are = sin with the latitude and the longitude . The vertical coordinate is = p/p s with p and p s denoting pressure and surface pressure, respectively, and f is the Coriolis parameter. See Hoskins and Simmons (1975) and Frisius et al. (1998) for more details on the model formulation.
The model is driven by Newtonian cooling for the diabatic heating, i.e. the model temperature is relaxed toward a prescribed reference temperature. The dissipative processes in the atmosphere are parametrized using Rayleigh friction, which describes the effects of surface drag and vertical transport of the horizontal momentum due to small-scale turbulence in the boundary layer, and damps vorticity and divergence towards a state of rest with the time-scale F = 1 day. The term H denotes the hyperdiffusion in the vorticity equation. The hyperdiffusion parametrizes both the subgrid-scale horizontal mixing and the energy cascade into these scales and its subsequent dissipation, because the dissipative range of the wave-number energy spectrum is not included with the relatively coarse model resolution (Seiffert et al., 2006). For a prognostic variable , the subgrid-scale hyperdiffusion in the model is given by: where the diffusion coefficient K is expressed as ) ℎ with a being the radius of the earth. We use h = 4 and the spectral modes with total wave number n = n T are damped with a time-scale H = 1 ∕ 4 day. The hyperdiffusion term is applied in the model so that we can ensure that it reaches its highest (lowest) impact at the smallest (largest) resolved scales. We run the model at varying horizontal resolutions (from low to moderately high) represented by triangular truncation of spherical harmonics with the T21, T31, T42 and T127 resolutions. On a latitude × longitude grid, these model resolutions correspond to 32 × 64, 48 × 96, 64 × 128 and 192 × 384 grids, respectively, which is approximately equivalent to 625, 415, 310 and 104 km horizontal resolutions, respectively, at the Equator. We take the model run with T127 horizontal resolution as the reference simulation in this study. Our simulations are with 10 vertical -levels. We use an aqua-planet set-up without topography. The model uses time-steps of 60, 40, 30 and 4 min for carrying out the simulations at T21, T31, T42 and T127 resolutions, respectively. The model runs are carried out for a period of 10 years out of which the first 2 years have been discarded as spin-up period. We use daily data of 8 years for the analysis presented in this article.
At low horizontal resolutions, the excessive subgrid-scale kinetic energy dissipation (due to hyperdiffusion) affects the quality of the simulations. It not only weakens the eddy variability, it also affects the propagation of waves and their interaction with the mean flow (Zurita-Gotor et al., 2015). To address this problem in an effort to make the low-resolution models more realistic, we use stochastic energy backscatter (SEB) and deterministic energy backscatter (DEB) parametrizations. For the flow-dependent SEB parametrization, we combine the schemes proposed by Berner et al. (2009) with that of Jansen and Held (2014). The key idea is to use the first-order flow-dependent auto-regressive stochastic forcing function (Berner et al., 2009) in which the subgrid-scale kinetic energy due to hyperdiffusion is injected back at each time-step (Jansen and Held, 2014) to make the scheme energetically consistent. Following Berner et al. (2009), we introduce a stochastic perturbation in the vorticity obtained from Equation 1 in the form of a first-order autoregressive (AR1) process: where n is the total wave number and (0, 1) is Gaussian white noise with zero mean and unit variance. We choose = 0.2 and p = −1.27 following Berner et al. (2009). The 1 and 2 are appropriately chosen constants at different horizontal resolutions (discussed later). The 1 is the fractional coefficient of the linear autoregressive parameter (1 − ), and 2 is the fraction of the wave-number-dependent noise amplitude. The perturbation kinetic energy input per unit mass (ΔE) (due to subgrid-scale hyperdiffusion) is computed at each time step and injected back into the resolved model scales at run time. The ΔE at each time step is computed with the help of the vorticity tendency term due to hyperdiffusion as: where n, m and N are total, zonal and truncation wave numbers, respectively. In the SEB parametrization scheme, we use horizontally uniform ΔE at each vertical model level in contrast to Jansen and Held (2014) who used the same temporally varying ΔE for all vertical levels and Berner et al. (2009) who used a spatio-temporally constant ΔE.
For the energetically consistent DEB parametrization, we add a negative biharmonic hyperdiffusion backscatter term H DEB in Equation 1 (with h = 2 in Equation 2) following Jansen and Held (2014). In other words, the net hyperdiffusion H net in the vorticity Equation 1 modifies to: where determines the fraction of kinetic energy (ΔE) dissipated by the hyperdiffusion and available for backscatter into the resolved flow. The goal is to bring the low-resolution models as close to the reference simulation as possible. We experienced no issues with numerical stability of our two backscatter schemes.

RESULTS AND DISCUSSION
We begin by showing the kinetic energy spectrum of PUMA at horizontal resolutions of T21, T31, T42 and T127 at 500 hPa ( Figure 1, dashed curves). It is clear from the figure that the impact of the hyperdiffusion is quite strong at high wave numbers in low-resolution simulations. To investigate as to whether the energetically consistent SEB and DEB parametrizations can help us alleviate this problem, we apply these schemes in PUMA at low horizontal resolutions. The kinetic energy (KE) wave-number spectrum of the stochastically forced model with the T21, T31 and T42 horizontal resolution is shown in Figure 1a. We see from Figure 1a that our scheme substantially reduces the excessive subgrid-scale KE dissipation at smaller scales. For example, the effect of hyperdiffusion is almost absent in the low-resolution model simulations with the SEB parametrization. The KE spectrum of the T21+SEB, T31+SEB and T42+SKB simulations follows a k −3 power law and almost exactly overlaps with the corresponding spectrum of our T127 reference simulation at large wave numbers, thus demonstrating the quality of our  Figure 1b. We notice that the KE spectrum of the T21+DEB, T31+DEB and T42+DEB simulations are much improved when compared to the simulations without backscatter. The KE spectrum in this case also closely follows the k −3 power law. However, the kinetic energy spectrum in the DEB experiments is not as good as the corresponding spectrum of the SEB experiments. Any further increase in the fraction of energy re-injected into the system in the DEB experiments distorts the KE spectra at smaller scales. The values of 1 and 2 used for SEB (in Equation 3) and for DEB (in Equation 4) simulations are given in Table 1. The values in Table 1 include the non-dimensional time step. Interestingly, 1 shows an exponential decay, whereas 2 and show a power-law scaling with truncation wave number N as follows: 1 = 28exp(− 1 N), 2 = 145.5 − 2 and = 9.23N − ; with N = 21, 31, 42 and the exponents 1 = 0.3, 2 = 0.88, = 0.73. This suggests that our parametrization schemes are scale-adaptive; thus, no re-tuning of the parameters is needed when changing the horizontal resolution. A PUMA simulation at T85 confirms that our power-law relationship is able to successfully predict the parameter values at different resolutions (not shown).
To demonstrate the advantages of our backscatter schemes, we show in Figure 2 the eddy variability produced by these models at different horizontal resolutions. The zonally averaged eddy kinetic energy (EKE), eddy heat flux (EHF) and eddy momentum flux (EMF) respectively. Since the model is hemispherically symmetric, we are showing the eddy variability over the Northern Hemisphere only. We see from the figures that the EKE, EHF and EMF systematically become weaker as the resolution is reduced compared to our reference simulation of T127. The eddy variables are especially weak in the T21 model simulation, suggesting that baroclinic instability processes are not well resolved. When comparing with the T127 simulation, we find that not only the magnitude of the jet is weaker, but also its position is not correct in the T21 simulation. When applying the SEB scheme, we see a noticeable improvement in the EKE and eddy fluxes in all the low-resolution models with maximum improvement seen in the T21+SEB simulation. Not only the magnitude of EKE and eddy fluxes have improved significantly, the extratropical jets are now located at the correct position in the SEB simulations. Similarly, we also notice improvements in eddy variability exhibited by low-resolution models when applying the DEB scheme, but the improvement is not as much as is obtained with the SEB scheme. To further diagnose the reduction in biases with subgrid-scale energy backscatter, we show in Figure 3 the difference map of EKE, EHF and EMF. The differences of these variables for all the low-resolution model simulations (with and without energy backscatter) are computed with respect to the corresponding values from our T127 reference simulation. It is clear from Figure 3 that the low-resolution simulations (T21 and T31) with the SEB scheme have the highest reduction in biases. The application of the DEB scheme in low-resolution models also reduces the biases, but not as effectively as with the SEB scheme. This is especially the case for very low-resolution models at T21 and T31. With increasing resolution, the model runs with SEB and DEB schemes move closer to each other. However, it is worth noting from Figure 3, that the application of the SEB scheme introduces extra energy at those grid points also where it is not desired (especially in the T21+SEB simulation). This drawback is less pronounced in the model runs with the DEB scheme. For example, for [T127-(T21+SEB)], even though we see a substantial reduction in positive EKE differences as compared to the [T127-T21] figure in the extratropical jet location, however, we also see high negative EKE differences elsewhere. On the other hand, [T127-(T21+DEB)] with higher positive differences at the jet location as compared to [T127-(T21+SEB)] shows lesser negative differences at other grid points. The middle and bottom panel of Figure 3 clearly shows that the representation of eddy fluxes is greatly improved in low-resolution models when applying the SEB and DEB schemes. We also quantify the improvements due to the implementation of stochastic and deterministic parametrizations in terms of the root-mean-square error (r.m.s.e.). The results are summarized in To further diagnose the improvement in the quality of model simulations on the application of the KE backscatter parametrization, we analyse in Figure 4a   (i.e. very long persistence) as compared to the reference simulation, which is improved in the T21+SEB and T21+DEB simulations, with the T21+SEB simulation giving more realistic autocorrelation time-scales when compared to the reference simulation. With an increase in horizontal resolution, the autocorrelation time-scale of the perturbed simulations moves closer to the reference simulation. We show in Figure 4b the difference between the autocorrelation time-scale of the reference simulation and the coarse simulations with and without backscatter. The differences clearly suggest that the parametrization schemes improve the autocorrelation time-scale (and hence predictability) of coarse-resolution models. We notice from Figure 4 that the SEB scheme works better than the DEB scheme for the T21 and T31 resolutions, whereas, at T42 resolution, the results from both the simulations are comparable. We also compute the first Empirical Orthogonal Function (EOF1) of the zonally averaged zonal wind to see the effect of energy backscatter on the dominant mode of variability. We show in Figure 5 EOF1 of different simulations with and without energy backscatter. The EOF1 of the reference simulation (T127) shows an annular mode of variability in the midlatitude upper troposphere with a dipole structure (of opposite signs) centred at approximately 35 • N and 55 • N. We see that the low-resolution models (T21 and T31) also capture the dipole structure though with weaker amplitudes as compared to the reference simulation. Furthermore, the position of the dipoles is not correctly simulated in the T21 simulation and the annular mode is shifted equatorwards from its mean position. When applying the SEB and DEB parametrizations, we see that not only the position of the annular mode is now correct, but that also the magnitude is improved (for example, in T21+SEB and T21+DEB simulations).
We also analyse the autocorrelation function (ACF) of the first principal component (PC1) of zonally averaged zonal wind ( Figure 6). In Figure 6a we show the ACF up to a lag of 100 days for the T21, T42 and T127 simulations. The e-folding time-scale of the T127 reference simulation is 25 days. We see that the ACF values of the T21 simulation are higher than the reference simulation up to approximately the e-folding time-scale. These values become much closer to the reference simulation for smaller lags by the use of the SEB parametrization (T21+SEB). However, at higher lags, the ACF values of the T21+SEB simulation become higher than in the reference simulation. The effect of SEB is less pronounced at T42 resolution. The T42+SEB simulation shows slightly higher (lower) ACF values for smaller (larger) lags. The application of the DEB parametrization on the T21 and T42 resolutions is depicted in Figure 6b. We see that in the T21+DEB simulation, the ACF values at all lags remain much higher compared to the reference simulation and represent an unrealistic persistence time-scale. The ACF values of the T42+DEB simulation become as good as the reference simulation up to a lag of approximately 15 days, after which they decrease before becoming higher than in the reference simulation after a lag of 60 days. It is important to note from Figure 6a,b that the minima of the T42 ACF curve shift towards lower lags when applying the backscatter schemes, with an even stronger shift for the DEB scheme. For example, in the ACF curves of the T42, T42+SEB and T42+DEB simulations, the minimum value is obtained approximately at a lag of 70, 60 and 50 days, respectively. This suggests that the "return of skill" or "rebound in predictability" (Guo et al., 2012) occurs much earlier in the T42+DEB simulation as compared to the simulations without backscatter. To examine the effect of our backscatter schemes more closely, we show in Figure 6c,d the ACF up to the lag of 30 days only. The e-folding time-scale of 26 days at the T42 horizontal resolution becomes 25 days and 20 days in the T42+SEB and T42+DEB simulations, respectively. On the other hand, at the T21 horizontal resolution, the e-folding time-scale, which was 29 days without backscatter, decreases to 23 days in the T21+SEB simulation, whereas it increases to 45 days in the T21+DEB simulation. On comparing Figure 6a,c with Figure 6b,d, we find that the ACF values with the SEB parametrization are much closer to each other for smaller lags than at higher lags.
To further highlight the practical advantage of using energy backscatter schemes, we collected the information regarding CPU time for running PUMA at different horizontal resolutions with and without energy backscatter schemes (summarized in Table 3). We see from the table that the additional computational expense of using DEB parametrization on the CPU time is very small. For example, the T21 simulation without energy backscatter takes 8 s as against 9 s with DEB, to complete 1 year of simulation on an Intel Xeon E7-8837 Processor at 2.67 GHz running PUMA in parallel on four cores. Similarly, for the T31 (27 s) and  T42 (69 s) resolutions, we notice an increase of just 2 and 3 s, respectively, when using the DEB scheme. Even though on using the SEB scheme, the CPU time to complete the same simulation increases as compared to the corresponding values of the DEB scheme, however, the increase is not too high. On the other hand, with the same computational resources, a T127 reference simulation takes 8,278 s to complete 1 year of simulation (higher by three orders to T21 and two orders of magnitude to T42). Considering the robust improvements obtained in the coarse-resolution simulations using the SEB/DEB parametrization schemes, their use in place of higher-resolution models can help save a lot of computational time without compromising much of the quality of the simulations.

CONCLUSIONS
Low-resolution atmosphere-ocean models suffer from the problem of excessive subgrid-scale dissipation. The artificial loss of kinetic energy at smaller scales affects the model's ability to correctly capture the mean flow and eddy variability, thus making them of limited practical use. On the other hand, the construction of robust and realistic low-resolution models is of practical importance, especially for long-term climate, palaeoclimatic studies and climatic extremes. In this article, we demonstrated that the energetically consistent subgrid-scale kinetic energy backscatter parametrization schemes are able to alleviate this problem to a significant extent. We apply scale-adaptive stochastic energy backscatter (SEB) and deterministic energy backscatter (DEB) parametrization schemes in the low horizontal resolution spectral model PUMA for this purpose. We have shown that the SEB scheme performs better as compared to the DEB scheme at the low T21 and T31 horizontal resolutions. We also find that with increasing resolution (T42), the performance of these schemes become comparable. It is important to mention here that for the backscatter schemes presented here, the purpose is not so much the nonlinear interactions of unresolved scales with the resolved flow, but the reduction of spurious dissipation by the viscosity scheme. This is a crucial difference, since the backscatter is not intended as a scheme that highly correlates with the unresolved subgrid turbulence, which is highlighted by the fact that both the stochastic and the deterministic schemes perform nearly equally well although the structure of the respective backscatter operators is quite different.
The quality of our schemes is comparable to the deterministic backscatter scheme in Zurita-Gotor et al. (2015) for an atmospheric model. In contrast to the European Centre for Medium-range Weather Forecasts (ECMWF) model (Berner et al., 2009), our schemes are energy consistent, which will be of importance when performing long-term climate simulations. The main purpose of the ECMWF model, however, is extended-range forecasts where energy consistency might be of lesser importance.
Our analysis clearly suggests that the application of these parametrization schemes in low-resolution models greatly improves the eddy variability. The application of the SEB scheme generates better eddy variability as compared to the DEB scheme. However, this is achieved at the cost of producing extra energy in other areas where it is not desirable. We also show that the autocorrelation time-scale of low-resolution models becomes more consistent with the reference simulation (T127) when they are run with the SEB or DEB parametrization schemes. Furthermore, it is demonstrated that the dominant mode of variability is realistically reproduced in model simulations using the SEB and DEB schemes. These schemes not only improve the magnitude of the dominant mode of variability, but its position is also correctly captured. It is also shown that the SEB and DEB schemes have little impact on the predictability (e-folding) time-scale of low-resolution models.
Our future work will focus on combining the stochastic and deterministic energy backscatter schemes. It will also be worthwhile to see the advantage of these schemes in more realistic models and set-ups for practical applications. While our schemes can be straightforwardly implemented in spectral models, the implementation of stochastic backscatter (e.g. Gugole and Franzke, 2019) and deterministic backscatter (e.g. Juricke et al., 2019) schemes is more involved in grid-point models. For this, novel efficient ways of deriving the noise covariance matrix are needed. In spectral models like PUMA, the noise covariance matrix can be diagonal which is not possible in grid-point models.