## Abstract

Kilometer-scale models allow for an explicit simulation of deep convective overturning but many subgrid processes that are crucial for convective initiation are still poorly represented. This leads to biases such as insufficient convection triggering and late peak of summertime convection. A physically based stochastic perturbation scheme (PSP) for subgrid processes has been proposed (Kober and Craig) that targets the coupling between subgrid turbulence and resolved convection. The first part of this study presents four modifications to this PSP scheme for subgrid turbulence: an autoregressive, continuously evolving random field; a limitation of the perturbations to the boundary layer that removes artificial convection at night; a mask that turns off perturbations in precipitating columns to retain coherent structures; and nondivergent wind perturbations that drastically increase the effectiveness of the vertical velocity perturbations. In a revised version, PSP2, the combined modifications retain the physically based coupling to the boundary layer scheme of the original scheme while removing undesirable side effects. This has the potential to improve predictions of convective initiation in kilometer-scale models while minimizing other biases. The second part of the study focuses on perturbations to account for convective initiation by subgrid orography. Here the mechanical lifting effect is modeled by introducing vertical and horizontal wind perturbations of an orographically induced gravity wave. The resulting perturbations lead to enhanced convective initiation over mountainous terrain. However, the total benefit of this scheme is unclear and we do not adopt the scheme in our revised configuration.

## 1. Introduction

Numerical weather prediction (NWP) models are now routinely run at horizontal grid spacings of a few kilometers, which permits an explicit representation of deep convection. However, while larger storm structures are sufficiently realistic at these resolutions (Weisman et al. 1997; Clark et al. 2016), the initiation of convection is often tightly linked to unresolved boundary layer processes. Turbulent eddies, for instance, are essential for overcoming the convective inhibition layer. Because boundary layer turbulence parameterizations typically only represent the mean mixing caused by subgrid eddies in a deterministic way, models can exhibit biases caused by an incorrect coupling between subgrid processes and resolved convection. In particular, this has been observed in the Consortium for Small-Scale Modeling (COSMO) model, which we are using in this study (Keil et al. 2014; Baldauf et al. 2011; Craig et al. 2012).

This inconsistency becomes particularly relevant when synoptic forcing is weak and local mechanisms are the main driver for overcoming convection inhibition (Keil et al. 2014). In these situations insufficient convective initiation has been commonly observed for kilometer-scale models (see e.g., Clark et al. 2016). In operational NWP systems such biases are often compensated by tuning other parameters, such as the turbulent length scale in the boundary layer parameterization (Hanley et al. 2015). Smaller mixing lengths allow the lowest model levels to heat up more, resulting in earlier and more deep convection (Hanley et al. 2015; Baldauf et al. 2011). For example, for the operational COSMO version, a decision was made to decrease the turbulent length parameter tur_len in order to reduce a lack of summertime convection (Baldauf et al. 2011). This, however, comes at the expense of temperature biases, which are widely observed in different studies (Baldauf et al. 2011; Schraff et al. 2016; Leutwyler et al. 2017; Necker et al. 2018). Baldauf et al. (2011) argue, however, that “the improved forecasts of deep convection outweigh those disadvantages” (p. 3902, l. 8, 9) and also climate simulations with a kilometer-scale COSMO version exploit this aspect (see e.g., Leutwyler et al. 2017). Such biases, however, can pose a severe problem for data assimilation systems (Schraff et al. 2016; Necker et al. 2018) and an improvement of deep convection without the loss of temperature skill is desirable.

One mathematically well-founded method to improve the coupling between subgrid and resolved processes are stochastic parameterizations [see Berner et al. (2017), for a review]. While theoretical studies rigorously derive stochastic representations of model uncertainty for idealized models (e.g., Chapron et al. 2018; Vissio and Lucarini 2018; Demaeyer and Vannitsem 2017, 2018), more pragmatic approaches dominate applications in complex NWP models today. Holistic approaches as in Hirt et al. (2018) or the Stochastically Perturbed Parameterization Tendencies (SPPT) (Buizza et al. 1999) and Stochastic Kinetic Energy Backscatter (SKEBS) (Shutts 2005) approaches provide a bulk representation for many types of model error. In contrast, process-level stochastic perturbations directly represent model uncertainties at its known sources: Ollinaho et al. (2017) and Jankov et al. (2017, 2019) propose stochastic perturbations of uncertain model parameters (SPP); Christensen et al. (2017) update the SPPT scheme to treat different processes independently; Craig and Cohen (2006) introduce a physically motivated approach based on statistical mechanics. Such process-level stochastic parameterizations have shown to better fit observations compared to bulk approaches (Christensen et al. 2017).

In this paper, we focus on process-level stochastic perturbations to improve the representation of convective initiation in kilometer-scale models. In the first part of the paper we revise the physically based stochastic perturbation (PSP) scheme for boundary layer turbulence originally developed by Kober and Craig (2016); in the second part we introduce a stochastic perturbation scheme for subgrid-scale orography (SSOSP). The developments introduced here will improve the consistency of the PSP scheme, and add additional physical processes. The aim of the present paper is to document these changes and additions, and to provide a preliminary evaluation of their impacts on simulations of convection to determine if the anticipated physical effects are present. Based on these results, a revised version of the PSP scheme will be defined that has the potential to improve forecast skill for convective precipitation. To demonstrate improved skill will require a systematic evaluation of ensemble forecasts over a long test period. This will be the subject of a future study.

The PSP scheme (Kober and Craig 2016) aims to better couple subgrid turbulence to convective initiation by reintroducing missing subgrid variance to the resolved scales. This is done by adding stochastic perturbations to the model tendencies of temperature, humidity, and vertical velocity. The perturbation structure is given by a horizontally correlated random field that is regenerated every 10 min, a time span that represents the average lifetime of a boundary layer eddy. Since such eddies span the entire depth of the convective boundary layer, the perturbations are kept vertically coherent by applying the same random field at all model levels. Most importantly, the amplitudes of the perturbations are scaled by physical information on the subgrid variance computed in the turbulence parameterization. This links the perturbation strength to the physical processes that cause small-scale turbulence. Furthermore, this enables the scheme to adapt to different atmospheric situations.

The PSP scheme showed promising results in a range of case studies (Kober and Craig 2016). Under weak synoptic forcing, the perturbations trigger more convective cells and consequently produce more precipitation. This improved the diurnal cycle and amplitude of precipitation. Further, Rasp et al. (2018) have successfully exploited its capability to create different realizations of the convective cloud field for the same synoptic situation by generating different random fields for each ensemble member. The magnitude of the ensemble spread created by different PSP realizations was also found to be comparable to operational ensembles (Keil et al. 2019). And in a convective-scale data assimilation framework, the PSP scheme produced a beneficial impact in ensemble spread and precipitation forecasts (Zeng et al. 2019, manuscript submitted to *Mon. Wea. Rev*.).

However, the original version of the PSP scheme also leads to undesirable side effects. In the first part of this paper, we investigate these issues and propose modifications to alleviate them. In particular, four topics will be covered:

In its current formulation, the evolution of the random field is temporally discontinuous, being regenerated every 10 min. To avoid these discontinuities, we introduce a continuously evolving random field using an autoregressive (AR) process.

Unrealistically high nighttime precipitation has been observed using the PSP scheme (Rasp et al. 2018). As we will show, this is related to perturbations in the free troposphere. These do not fit with the rationale of the PSP scheme, which is designed to represent vertically coherent structures in convective boundary layers. Hence to remove these artifacts we restrict the perturbations to the boundary layer.

Furthermore, organized convective structures were found to be broken apart by the perturbations. To prevent this, we investigate the use of a mask to switch off the perturbations in regions where precipitation is already occurring.

Finally, the vertical velocity perturbations were found to have an insignificant effect compared to the temperature and humidity perturbations (Kober and Craig 2016). As we will discuss in this paper, this is associated with the generation of acoustic modes that lead to rapid dissipation of the imposed vertical velocity perturbations. By introducing balanced, 3D wind perturbations, we obtain lasting and effective wind perturbations.

After evaluating the modifications separately, we propose a revised version of the scheme, which we will call PSP2.

In the second part of this manuscript, we target the influence of subgrid-scale orography (SSO) on triggering convection as a complementary process to turbulence. While the importance of resolved orography on convective initiation is unquestioned (Houze 2012), the contribution of subgrid-scale orography is unclear. For example, Kirshbaum et al. (2007b) find that scales smaller than few kilometers are not dominant for convective initiation when a range of orographic scales exist. On the other hand, Kirshbaum et al. (2007a,b) also show that small, single-scale orographic features can dictate the spacing of orographic bands, and Langhans et al. (2011) conclude that small-scale orography contributes significantly to convective initiation. Since orography can act as a source of predictability for convective scales (Bachmann et al. 2019a,b, manuscript submitted to *Mon. Wea. Rev*.; Anthes 1986), the accurate representation of its effects will be important for weather prediction.

Most current NWP models include a parameterization for SSO. However, these parameterizations primarily account for a deceleration of the flow by unresolved orography, either by modifying the roughness length or by including an orographic drag term. Here, we aim to account for the mechanical lifting caused by SSO and its effect on convective initiation with a newly developed stochastic perturbation scheme, called SSOSP. The new scheme closely follows the formulation of the PSP scheme: wind tendencies are randomly perturbed with an amplitude that scales with theoretical gravity waves excited by SSO. Investigating the impact of this scheme will allow us to estimate the importance of SSO on convective initiation and therefore contribute to the discussion mentioned in the previous paragraph.

In section 2, the formulation of the PSP scheme is explained in detail, as well as the four modifications mentioned above. Further, the SSOSP scheme is introduced. In section 3 we describe our experimental setup, which is used for testing both PSP and SSOSP. The results are presented in section 4 for the PSP scheme and in section 5 for the SSOSP, followed by a summary and discussion in section 6.

## 2. Conception of the stochastic perturbations

In this section we summarize the original perturbation scheme for boundary layer turbulence (PSP) by Kober and Craig (2016) and then describe our modifications to improve the physical consistency of the scheme. At the end of this section, we introduce the stochastic perturbations to account for mechanical lifting by subgrid-scale orography, the SSOSP scheme.

### a. Physically based stochastic perturbations for subgrid turbulence (PSP)

The PSP scheme for boundary layer turbulence was first described in Kober and Craig (2016). Here we present an updated formulation according to Rasp et al. (2018).

Boundary layer turbulence is characterized by eddies with scales ranging from millimeters to around 1 km, the boundary layer height. Traditional parameterizations that account for the impact of turbulence on the resolved grid, were developed for grid sizes that are large relative to the turbulent eddies. The parameterization predicts the mean effects of a large ensemble of eddies, and the variability associated with the individual eddies is averaged out. But the kilometer-scale region represented by a grid box in a convective-scale model will contain only a small number of the largest eddies and their variability will have a substantial effect on the gridbox mean. As a result, the variability on the smallest resolved scales in the model is underrepresented, which can lead to the difficulties in overcoming convective inhibition that were described in the introduction.

The PSP scheme reintroduces the influence of the lost small-scale variability by adding perturbations to the tendencies of the temperature, humidity, and vertical velocity field (*T*, *q _{υ}*,

*w*) on the smallest effectively resolved scale, Δ

*x*

_{eff}= 5Δ

*x*(see e.g., Bierdel et al. 2012). For this we assume a random, horizontal eddy field

*η*(

*x*,

*y*,

*t*) with temporal and spatial correlation scales

*τ*

_{eddy}and Δ

*x*

_{eff}. The temporal correlation is given by the typical lifetime of convective eddies (i.e.,

*τ*

_{eddy}= 10 min), and is also theoretically justified by the missing scale separation between resolved and subgrid scales that necessitates the representation of memory effects (see e.g., Berner et al. 2017). The perturbation amplitude is chosen to be proportional to the subgrid standard deviation $\Phi \u20322\xaf$ for Φ ∈ {

*T*,

*q*

_{υ},

*w*}, which is taken directly from the turbulence scheme. The COSMO model uses a 1.5-order closure (Raschendorfer 2001) based on level 2.5 of Mellor and Yamada (1982), in which the second moments are diagnostically computed based on the turbulence kinetic energy and the vertical gradients of the variables in question. In addition, the scheme adapts to the grid spacing. For uncorrelated subgrid eddies the variability scales with $1/Neddy$, where $Neddy=\Delta x2/leddy2$ is the number of eddies of scale

*l*

_{eddy}= 1000 m in a grid box (Craig and Cohen 2006). Dividing by

*τ*

_{eddy}ensures that the accumulated perturbations during this time period scale as $(leddy/\Delta xeff)\Phi \u20322\xaf$. The complete mathematical formulation is

where *α*_{tuning} is a free parameter that should be of order one, if the length, time and variance scales that appear in the equation have the correct orders of magnitude. A single choice of *α*_{tuning} should apply for all weather regimes and model resolutions, if the physical assumptions of the scheme are satisfied.

### b. Modifications to the PSP scheme

#### 1) Autoregressive process: PSP-AR

In the original version of the PSP scheme, a new random field *η* was generated every *τ*_{eddy} = 10 min and was held constant in between. To ensure a more realistic evolution, we introduce a continuously evolving random field by using an autoregressive (AR) process:

The autoregressive parameter *σ*_{τ} defines the temporal correlation scale as described below. A random field *ε*_{t} is generated at each time step *t* with same spatial correlation scale as *η*, and standard deviation 1 and mean 0. The parameter *σ*_{ε} scales *ε* as described below. To ensure a stationary process (mean and variance remain constant with time) that contains “memory” (i.e., subsequent time steps are positively correlated) we require 0 < |*σ*_{τ}| < 1 (Wilks 2005). To determine the two scaling parameter *σ*_{τ} and *σ*_{ε} we consider two constraints. First, the variance of the random field *η*_{t} is required to be 1, which leads to the relationship $\sigma \epsilon =1\u2212\sigma \tau 2$. Second, the autoregressive parameter *σ*_{τ} is defined in terms of the characteristic time scale for decay of the temporal correlation between two fields, *τ*. We take this correlation time scale *τ* as the exponential decay scale of the autocorrelation function $\rho \u2061(k)=\sigma \tau k=ekln\sigma \tau $, where *k* corresponds to the lag *t* − *t*_{0}. Then, the correlation time scale is linked the autoregressive parameter *σ*_{τ} by *τ* = −(1/ln*σ*_{τ}) [see Wilks (2005), for more details]. We set *τ* to the time scale at which the random field was modified in the original PSP scheme (i.e., *τ*_{eddy} = 10 min). The original PSP scheme with this modification will be denoted as PSP*-*AR.

#### 2) Adaptive boundary layer height restriction: PSP-HPBLcut

It will be shown in section 4a(2) that the original PSP showed unrealistic precipitation at night, which can be related to considerable perturbations in the free troposphere present also during the night. However, the PSP was developed with the convective boundary layer in mind and many of its assumptions (notably vertically correlated eddies) are not valid for free tropospheric turbulence. In order to improve the physical integrity of the scheme and simultaneously reduce unrealistic precipitation, we limit the perturbations to the planetary boundary layer height *H*_{PBL}. We determine *H*_{PBL} as the height at which the bulk Richardson number reaches the critical Richardson number as implemented in the COSMO model. This value is 0.33 under stable conditions (Wetzel 1982) and 0.22 under convective conditions (Vogelezang and Holtslag 1996). Below *H*_{PBL} the perturbations are applied as usual. Above, however, the perturbations are decreased linearly to zero at a height of 500 m above *H*_{PBL}. This modification, applied to the PSP-AR version of the scheme, will be termed PSP-HPBLcut.

#### 3) Masking of already precipitating areas: PSP-mask

We observed that in certain situations the original PSP scheme appeared to break up organized convection [see section 4a(3)]. We hypothesize that this is partly caused by the perturbations disturbing already existing updrafts. As a simple ad hoc fix, we test turning off the perturbations in grid columns where the vertically integrated precipitating hydrometeor content exceeds 10^{−3} kg m^{−2}. This adjustment to the PSP-AR version will be labeled as PSP-mask.

#### 4) 3D nondivergent perturbations: PSP-3D

Kober and Craig (2016) observed that the impact of the original PSP scheme was caused by the temperature and humidity perturbations while the vertical velocity perturbations had essentially no effect. Further analysis shows that the *w* perturbations decay rapidly in a few time steps. This can be related to pressure perturbations: Chagnon and Bannon (2005), for example, show that an injection of vertical momentum efficiently induces high-frequency acoustic waves causing a rapid removal of the imposed vertical updraft. Following their approach we extend the scheme to add perturbations to the horizontal velocity field that approximate three-dimensional nondivergence. The goal is to suppress the excitation of spurious acoustic waves so that the vertical velocity perturbations persist (see e.g., Fiedler 2002; Chagnon and Bannon 2005; Edson and Bannon 2008).

Nondivergence is given by

We neglect modifications due to the terrain following vertical coordinate system of the COSMO-model and the metrical terms due to the spherical Earth as well as compressibility. Neglecting the distortion to the grid due to these effects greatly simplifies the computation, but implies that the resulting wind perturbations will not be exactly nondivergent, and some perturbation energy will be lost to acoustic waves. However, this lost perturbation energy will be small in comparison to the lost perturbation energy of the original PSP scheme where the wind perturbation was entirely divergent and almost all the energy was rapidly dissipated. The COSMO-model employs a rotated spherical coordinate system with Earth radius *R*, rotated longitude *λ*, and latitude *ϕ*, so that ∂*x* = *R*cos*ϕ*∂*λ* and ∂*y* = *R*∂*ϕ*. Due to the rotated coordinate system and the comparably small domain, we neglect the trigonometric terms for simplicity.

This results in the following approximate nondivergence criterion for the PSP setup:

We define the horizontal perturbations by the velocity potential *χ*: **v**_{h} = ∇*χ*. Consequently, the horizontal perturbations are nonrotational, which is in accordance with the incentive of the perturbations.

Computing the divergence of **v**_{h} and substituting in Eq. (4) we obtain the following Poisson equation:

For a vertical velocity perturbation (increment) given by Eq. (1) we derive *u* and *υ* perturbations by solving Eq. (5) for the corresponding perturbation of *χ*. A finite-difference approximation is obtained by replacing the Laplace operator with a five-point-stencil Laplace matrix.

To avoid large vertical velocity gradients ∂*w*/∂*z* near the lower boundary, we linearly increase *w* perturbations from 0 at the surface to the original value at *z*_{0}. We chose *z*_{0} = 500 m to include a sufficient number of vertical model levels while ensuring that the full perturbation amplitude is active at the height of the daytime inversion, where we expect the perturbations to be most effective.

Applying these balanced wind perturbations to PSP-AR will be denoted as PSP-3D.

### c. Stochastic perturbations to account for mechanical lifting by subgrid-scale orography (SSOSP)

As described in the introduction, orography plays an important role in triggering convection. The impact of orography that is smaller than model resolution [i.e., subgrid-scale orography (SSO)], needs to be parameterized. In most models, this is done by an orographic drag parameterization (Lott and Miller 1997) and a modification of the roughness length (Doms et al. 2011). The former accounts for deceleration of flow by orographically induced gravity waves and blocking and the latter approximates effects within the boundary layer. However, the vertical motion responsible for triggering of convection by SSO is not addressed in these parameterizations.

For this purpose we propose stochastic perturbations to account for convective initiation by subgrid-scale orography, which we will name the SSOSP scheme. Here we focus on the effects of mechanical lifting and the resulting gravity waves. The flow over a mountain causes an upward displacement of the air parcels (i.e., a vertical velocity perturbation), which may decay with height, or propagate away from the source, depending on the stratification of the atmosphere. The upward displacement may cause the parcels to reach their level of free convection (LFC). As the relative importance of single mechanisms for convective initiation is not well understood, we neglect processes such as thermally induced upslope winds, upward and leeside convergence and flow blocking by mountains (see e.g., Geerts et al. 2008; Schmidli 2013; Langhans et al. 2013; Houze 2012; Hassanzadeh et al. 2016; Panosetti et al. 2016).

Inspired by the PSP scheme, we introduce stochastic vertical velocity perturbations that scale with the amplitude of the SSO lifting and the induced gravity waves. The horizontal wind is perturbed to approximate nondivergence as described in section 2b(4). Since the errors in the estimation of the nondivergent wind perturbations are likely to be largest over strong orography where the grid is most distorted, there will be an additional source of noise in these regions. However, we were unable to identify any obvious effects in the simulation results.

The *w* tendency is perturbed according to

Here *η*(*τ*_{η}, *σ*) denotes a horizontally correlated AR random field as described in section 2b(1). Instead of *τ* = *τ*_{eddy} we assume that slower processes such as the synoptic wind field and the diurnal cycle dominate the evolution of orographically induced flow structures. Craig and Selz (2018) for example show that the time scales of gravity waves (even on scales of 1km) can range from minutes to days. Since no single time scale stands out, we test several values, *τ* = 30 min, 2, 5 h. As in PSP, *α*_{tuning} is a tuning parameter.

The variable *w*′ represents the physical scaling of the perturbations. It depends on the amplitude *w*_{0} (at given height) and vertical wavenumber *m* of a gravity wave induced by the SSO. This includes both propagating and decaying gravity waves and also represents the effect of mechanical lifting. Both *w*_{0} and *m* are considered in terms of forced gravity waves as in Gill (1982). The formulation for *w*_{0} and *m* in terms of horizontal wind speed, Brunt–Väisälä frequency *N* and SSO information is presented in appendix A. The gridbox standard deviation, slope, and orientation of SSO is derived from 30 m ASTER elevation data (Schättler et al. 2015). To avoid numerical instabilities, we decrease the perturbations linearly to zero between *z*_{max} = 500 m and the surface. We further constrain the perturbations to the planetary boundary layer, where we assume triggering processes to be most relevant: we impose an exponential decay of the perturbations at *z*_{max} that scales to *γ* = 1/2 km^{−1} for nondamped waves and to $max\u2061(\gamma 0,\u2212m)$ for damped waves. Mathematically, this can be expressed as

Further details are described in appendix A.

## 3. Model simulations, observations, and simulation period

To evaluate the impact of the PSP variants and the SSOSP we compute several simulations with different setups. We investigate each modification of PSP individually on one day with weak synoptic forcing of convection, 6 June 2016. For comparison, we use 30 May 2016 as an example of a strong forcing day with synoptically driven, coherent precipitation. Based on these results, we decide on a suitable configuration for a revised PSP version, PSP2, and compute 10 subsequent days with total simulation time of 24 h to confirm our findings. For the SSOSP scheme we focus on 6 June 2016. See Table 1 for an overview of the experiments.

### a. Simulation period

We select a 10-day period from 29 May to 7 June 2016 for our experiments, which was characterized by heavy precipitation over Germany (Piper et al. 2016). The first half of the period was dominated by an upper-level trough and associated low pressure systems over Germany. The resulting synoptic lifting combined with southeasterly advection of moist air caused heavy rainfalls, particularly on 30 May (Fig. 1). In the second half of the period, the trough made way for a persistent ridge structure reminiscent of an omega-block. Then the synoptic situation was calm, allowing large instability to build up each day followed by strong convection. This behavior was particularly dominant on 6 June. This period provides a good testing ground for the stochastic schemes, offering both strongly and weakly forced convection.

### b. Model and simulation setup

We use the COSMO model (Baldauf et al. 2011), version 5.4 g, in a convection-permitting setup with Δ*x* = 0.025°, roughly 2.8 km, and a time step of 25 s. The setup uses 50 vertically stretched model layers ranging from 10 m above ground to 22 km above mean sea level. Shallow convection is parameterized using the Tiedtke scheme. Details on parameterizations can be found in Doms et al. (2011). The setup follows the operational COSMO-DE setup with 461 by 421 grid points centered over Germany at 50°N, 10°E.

The only major deviation from the current operational setup is a change of the tuning parameter tur_len in the boundary layer scheme to 500 m. As discussed in the introduction, a smaller value was applied as an ad hoc correction in the operational model to increase the rate of convective initiation. In this way we will examine whether the PSP scheme is able to achieve realistic convective initiation without this correction.

Initial conditions are taken from the COSMO-KENDA ensemble data assimilation system (Schraff et al. 2016), which is based on a local ensemble transform Kalman filter to directly assimilate conventional observations. Additionally, radar observations are assimilated through latent heat nudging. In our experiments we only use the deterministic analysis. Using initial conditions provided by a convective-scale data assimilation system reduces the model spinup compared to downscaled initial conditions. Boundary conditions are provided by global ICON forecasts. Our simulations are started at 0000 UTC on each of the simulation days and extend for 24 h.

A simulation without stochastic perturbations, also with tur_len = 500 m, is denoted as *Reference*.

### c. Observations

We compare our simulations to radar-derived precipitation fields. For this we use the Radar Online Aneichung (RADOLAN) quality-controlled radar observations (EY product) provided by the German Weather Service, which is based on European radar reflectivities and provides radar coverage for most of our domain (DWD 2018a,b). For all analyses of our simulations we restrict the geographical region to grid points where radar data is available.

For this study we focus on the systematic impact of the PSP scheme on the behavior of the model. The comparison to radar observation only serves as a guideline. Trying to improve the forecasts scores would require much longer simulation periods and tuning of several model parameters (e.g., tur_len). We anticipate, however, that the systematic impacts of the PSP scheme are independent of the specific model configuration.

## 4. Results

### PSP modifications

#### 1) Autoregressive Process in PSP-AR

To test the influence of the AR random field we focus on 6 June, a weakly forced day on which the perturbations are expected to have a large influence on the precipitation. The original PSP scheme (PSP), in which the random field is regenerated every 10 min, leads to an increase and earlier onset of precipitation (Fig. 2a) compared to the unperturbed reference simulation. We note that the diurnal cycle and amplitude of precipitation from radar observations (Fig. 2a) are already well matched by the *Reference* simulation. The PSP simulations consequently produce too much precipitation. This behavior is, however, not necessarily typical of longer periods [see section 4a(5)]. As mentioned above, due to the tuning problem of NWP, we concentrate on the systematic impacts of the PSP schemes, not their performance compared to observations. Switching to a continuously evolving random field (PSP-AR) further increases the precipitation amount. Example precipitation fields in Fig. 2b confirm that this increase in precipitation is a consequence of more convective cells and thus, more triggering.

We explain this behavior by the different meaning of *τ*_{eddy} in the PSP-AR setup: For the AR process *τ* represents the average decay scale while for single grid points and time ranges longer correlation times are possible. These might be more effective in triggering precipitation. The PSP-AR perturbation can be scaled by tuning *α*_{tuning}, which we will do after considering also the other modifications. Note also that these results are robust for different random seeds (not shown). In the remainder of this paper we will always use the AR random field, unless otherwise specified.

#### 2) Adaptive boundary layer height restriction in PSP-HPBLcut

Next, we investigate the impact of constraining the perturbations to the height of the planetary boundary layer (PSP-HPBLcut). Figure 3 shows the horizontal standard deviations of the perturbations. Without limiting the perturbation height (PSP-AR) the perturbations are largest in the boundary layer but significant perturbations in the midtroposphere are present as well. As discussed above this is conceptually inconsistent with the rationale of the PSP scheme. Limiting the height in HPBLcut removes this inconsistency. This modification leads to a decrease in the total precipitation amount and, most importantly, removes the precipitation peak at night (Fig. 2a).

#### 3) Masking of already precipitating areas in PSP-mask

In the original and the PSP-AR scheme, we observed a breaking up of larger-scale precipitation structures. An example for this can be seen on 30 May (Fig. 4c), a strongly forced convective day. A coherent, synoptically driven precipitation structure is present over western Germany and Belgium, which is well reproduced by the reference simulation. Adding the PSP-AR perturbations breaks up this structure into many smaller convective cells. Turning off the perturbations in precipitating columns PSP-mask corrects this problem but also reduces the impact on weak forcing days (see Fig. 2). A too small threshold for precipitating hydrometeors can be responsible for switching off the perturbations before deep convection has developed.

A quantitative measure for the structure of the precipitation field is the *S*(tructure) component of the SAL score (Wernli et al. 2008) displayed in Fig. 4b. Negative *S* values imply that the simulated precipitation cells are too small and peaked compared to the radar observations while *S* = 0 suggests a perfect match in terms of structure. For details on the computation, see appendix B. During the morning hours on 30 June, when the synoptic structure is most salient, the PSP-AR simulations show a clear reduction of the *S* score, indicating smaller cells. The precipitation mask reverts the structure back to its original value.

The precipitation mask solves the problem of the original scheme, where organized convection was broken up, but as noted above, also reduces the triggering effect in weak forcing situations.

#### 4) 3D nondivergent perturbations in PSP-3D

As a last modification to the PSP scheme, we investigate the impact of adding nondivergent wind perturbations. For these experiments we reduce the perturbation amplitude *α*_{tuning} to *α*_{1} = 2 because larger horizontal wind perturbations resulted in numerical instabilities. The domain averaged precipitation time series are displayed in Fig. 5a. As expected the precipitation increase is lower for the new perturbation amplitude *α*_{1} compared to the original perturbation amplitude *α*_{0} (see Table 1). We first investigate the impact of the original *w* perturbations without horizontal wind perturbations by comparing PSP-AR*, α*_{1} to PSP-AR *Tq*_{υ}, *α*_{1}, where only temperature and humidity are perturbed. Consistent with the results of Kober and Craig (2016), *w* has a negligible effect in this experiment.

Adding the *u* and *υ* components to enforce nondivergent wind perturbations results in a marked increase in precipitation even without the *T* and *q*_{υ} perturbations (PSP-3D *w, α*_{1}). Including perturbations to *T* and *q*_{υ} again (PSP-3D all*, α*_{1}) further increases the effect. This suggests that the effect of nondivergent wind perturbations is now of same order of magnitude as the *T* and *q* perturbations. We emphasize that the relative amplitudes of the imposed wind and thermodynamic perturbations are determined from the subgrid variances derived from the turbulence scheme.

We further investigate the impact of the nondivergent wind perturbations on the precipitation structure (Figs. 5b,c). PSP-3D all*, α*_{1} has fewer but larger convective cells compared to PSP*, α*_{0} at a comparable total precipitation amount. Again we look at the *S* score to quantify the convective structure. In general, the model simulations show negative *S* values. This confirms the visual impression that the cells are too small and intense and lack the stratiform precipitation regions in the observations. During the daytime the PSP-AR simulation has slightly improved scores compared to the unperturbed reference simulations. The nondivergent wind perturbations result in a further improvement, particularly during the time of maximum convection in the afternoon. Hence, balanced wind perturbations lead to larger cell structures.

In general, the nondivergent wind perturbations are very effective in triggering convection and also appear to induce mesoscale circulations that cause larger precipitation cells.

#### 5) Combining the modifications for PSP2

Finally, we choose a most suitable combination of modifications to PSP for further investigation. This includes the AR random field, HPBLcut and the nondivergent wind perturbations and will be denoted as PSP2. Since the 3D wind perturbations already improve the precipitation structure (see Fig. 4c, PSP2), we do not include the precipitation mask, which has the undesirable side effect of reducing convective initiation. Because of the many changes to the original scheme we start by determining a reasonable perturbation amplitude *α*_{tuning}. Corresponding domain averaged precipitation time series are displayed in Fig. 6 for *α*_{tuning} = 1, 1.5, 2. We select *α*_{final} = 1.5, which produces perturbation amplitudes comparable to the original PSP (see blue line in Fig. 2) and use this value henceforth.

We now simulate the entire high-impact weather period. For the analysis we divide it into a strong forcing period (29 May–2 June) and a weak forcing period (3–7 June). For both synoptic situations the modified PSP2 scheme produces an earlier peak of precipitation (Fig. 7a). This effect is desirable, since it addresses the known problem of late initiation of convection in kilometer-scale models. However, it should be noted that for this particular model and this weather period, the change does not result in an overall improvement in the precipitation forecast in comparison to radar observations.

For the weak forcing days the precipitation amplitude is also increased, while for the strong forcing days the amplitude is decreased slightly. This reduction could be caused by a reduction of CAPE in the evening due to earlier triggering or the 3D wind perturbations.

We again look at the *S* component of the SAL score to investigate the impact on precipitation structure (Fig. 7b). For the weak forcing period the PSP2 scheme shows a shift toward more coherent precipitation. This is a change in the right direction but compared to observations the convective cells are still too small scale.

Finally, in Fig. 7c we look at a common metric for assessing precipitation forecast skill, the fraction skill score (FSS). Displayed is the FSS for a scale of 227 km with a precipitation threshold of 0.1 mm h^{−1}. Results have also been computed for other scales and thresholds (not shown). As expected the skill decreases for smaller scales and higher thresholds, however, the relative performance of the reference and PSP2 forecasts is similar for all values. In particular, the PSP2 perturbations lead to a slight improvement in skill for the weak forcing days, particularly during the time of maximum convective activity. For the strong forcing days, there is a skill increase during the time of convective initiation but a decrease of skill toward the evening hours. This is likely connected with the too early decrease in precipitation caused by the perturbations as discussed above. To evaluate the significance of these differences, longer periods need to be investigated.

To summarize, the new version, PSP2, includes the AR random field, the restriction of the perturbations to the boundary layer and the balanced wind perturbations. As a preliminary estimate, we retune the perturbation amplitude *α*_{tuning} to 1.5 to match the precipitation amplitude of the original PSP scheme. For operational use, this choice would have to be revisited to produce the best overall precipitation forecast over a long test period. Although a clear overall improvement in forecast skill scores was not demonstrated in the preliminary evaluation here, the specific impacts of the changes to the scheme were as expected. In particular, PSP2 maintains the desired effect of the original PSP scheme, namely an earlier peak and increased amplitude of convection precipitation in weak forcing situations, while several drawbacks have been eliminated. In addition, it improves the cell structure and the spatial forecast skill compared to the unperturbed simulation.

## 5. SSOSP scheme

We now turn our attention to the SSO perturbations. First, we investigate the impact on the precipitation amount for different values of the correlation time scale *τ* and perturbation amplitude *α*_{ref} where *α*_{ref} = (*α*_{tuning} × 2 h)/*τ*. This rescales the final perturbation amplitude for different values of *τ* with respect to a reference time scale *τ*_{ref} = 2 h. To see the effects of the orography more clearly, we analyze orographic and flat grid points separately by defining orographic/flat grid points as grid points where the SSO standard deviation is larger/lower than 20 m.^{1}^{,}^{2} Corresponding domain averaged precipitation is displayed in Fig. 8. For orographic grid points increasing *α*_{ref} from 0.18 to 0.42 leads to more precipitation but, counterintuitively, a further increase of *α*_{ref} causes a reduction in the precipitation amount. Interestingly, flat grid points only experience a reduction with increased *α*_{ref}. One hypothesis for this decrease is that the perturbations enhance mixing of boundary layer and free tropospheric air, reducing instability, and leading to a reduction in precipitation. Increasing *τ* monotonically leads to more precipitation. More slowly evolving perturbations allow stronger updrafts to develop. These stronger updrafts cause an increased probability of convective initiation.

For maximum impact we chose *α*_{ref} = 0.42 and *τ* = 5 h for further ensemble experiments with 10 members. The domain averaged precipitation time series in Fig. 9 shows reasonable spread between ensemble members for both orographic and flat grid points. This spread is comparable to spread generated by the original PSP scheme [see Kober and Craig (2016), Fig. 7d]. While spread over orographic regions is desired, further insight is necessary to understand the behavior over flat grid points. Moreover, the behavior of individual ensemble members suggests a significant increase in precipitation over orographic grid points and decrease over flat grid points. This is further illustrated in Fig. 9 where ensemble precipitation fields are displayed for two time steps: precipitation within the SSOSP scheme mainly occurs in the alpine region while most parts of Germany are precipitation free.

Overall, the SSOSP scheme shows a nonmonotonic dependence on perturbation amplitude *α*_{ref}, while increasing the correlation time scale *τ* appears to give increasing precipitation over orographic regions. For flat regions however, the scheme causes an undesired reduction of precipitation.

## 6. Summary and discussion

### a. Stochastic perturbations for subgrid turbulence (PSP)

In the first part of this paper, four concerns of the original PSP scheme were addressed and corrected by modifying different aspects of the original PSP scheme:

The autoregressive process in PSP-AR generates a physically more realistic, continuous evolution of the random field and hence refines the concept of the scheme. The performance is qualitatively comparable to the original PSP.

In PSP-HPBLcut we constrain vertically correlated perturbations to the boundary layer in accordance with the vertical extent of boundary layer eddies. This removes perturbations in the free troposphere that resulted in spurious nighttime precipitation.

By excluding perturbations in precipitating regions in PSP-mask we are able to reestablish organized convective structures. However, we decided not to include the PSP-mask adjustments in the revised PSP2 version for two reasons: first, the other modifications improved the breakup of organized structures already; second, this ad hoc fix has the undesirable side effect of reducing convective initiation overall.

In PSP-3D we included 3D nondivergent wind perturbations to prevent rapid dissipation of vertical velocity perturbations. Interestingly, these consistent wind perturbations are at least as effective as the

*T*and*q*_{υ}perturbations in triggering convection. The relative importance of buoyant and mechanical lifting in initiating convective updrafts in nature is still unknown. Both processes were found by Torri et al. (2015) to play a role in secondary initiation by cold pools. It also appears that the cell structures in the simulations with nondivergent wind perturbations were more comparable to radar observations. The relationship between triggering processes and the resulting convective structures is an important topic for future research.

We have combined the first, second and fourth modification to define the PSP2 scheme and selected a tuning parameter of *α*_{tuning} = 1.5 to match the precipitation amplitude of the original PSP. That the optimal value of the tuning parameter is of order one is an encouraging sign that our physical reasoning is appropriate. By evaluating PSP2 on a longer time period, we find its behavior well in accordance with the rationale of the original scheme: when synoptic forcing is weak we observe a small shift in the diurnal cycle and an increase in total precipitation amplitude; this could counteract long-standing biases in current NWP systems (Clark et al. 2016; Baldauf et al. 2011). During strongly forced episodes or at night the scheme has only a small impact.

In addition to the systematic impact of the PSP scheme on convective initiation, studies have also demonstrated its ability to increase ensemble spread (Rasp et al. 2018; Keil et al. 2019.). This has potential to help reduce the underdispersion in many current convective-scale ensemble forecasting and data assimilation systems (see e.g., Berner et al. 2017; Ollinaho et al. 2017). Direct radar assimilation might benefit especially from such a perturbation scheme where sufficient spread in the first guess ensemble is particularly necessary. One remaining, undesired effect of the PSP2 scheme is a decrease of precipitation toward the evening. We note that this has been identified also for the original PSP by Rasp et al. (2018) who related this to a lack of convective organization. Organization of convection is caused by processes other than the ones addressed with the PSP2 scheme, most notably cold pools. A better representation of these features and their ability to trigger new convection might be necessary to capture organized evening convection. We aim to address this in a follow-up study.

Overall, we conclude the PSP revisions by recommending PSP2 instead of the original PSP: PSP2 contains comparable advantages as the PSP but without most of its drawbacks.

The question whether the PSP/PSP2 scheme improves the forecasts is more challenging to answer. Dealing with a highly complex model with dozens of adjustable parameters inevitably conjures up the curse of model tuning (Palmer and Weisheimer 2011). Over many years the parameters in an operational model are modified to achieve the best forecast scores. This most likely invoked a multitude of compensating errors. Introducing a more physical parameterization, therefore, demands a careful retuning of the parameters, which is a considerable computational effort. Furthermore, simulations using a full ensemble prediction system over a sufficiently long time period are necessary to thoroughly validate the impact of the PSP2 scheme. We intend to address these challenges in future work to allow for a more conclusive evaluation of forecast skill. In this paper we have focused on understanding the systematic changes caused by the perturbation schemes rather than on absolute forecast scores. The acquired physical understanding will help to guide development of a potential operational forecasting system.

### b. Stochastic perturbations for subgrid orography (SSOSP)

In addition to the stochastic perturbations for turbulence, we have developed SSOSP, a scheme to account for subgrid orography and its impact on convective initiation. Doing so we consider mechanical lifting by SSO taking into account the stratification of the atmosphere. Again, reasonable values of the tuning parameter *α*_{tuning} are of order one, supporting our physical reasoning.

The impact of the SSOSP on precipitation is not clear-cut. While we observe a precipitation increase over orographic regions the scheme also leads to a reduction of precipitation over flat terrain. Furthermore, there is a nonmonotonic relationship between the perturbation amplitude and the precipitation amount that may be related to a combination of mechanical triggering and a mixing-induced reduction of CAPE. An ensemble simulation with SSOSP showed that spread over orographic terrain is generated, but also over flat regions which is conceptually not desired.

These behaviors of the SSOSP make it hard to employ. Furthermore, the impacts are not clearly complementary to those of PSP2: the precipitation impact of SSOSP is either too weak or coupled with an undesired impact on flat regions, and the generation of ensemble spread by SSOSP can also be achieved with PSP2. For these reasons, we refrained from including the SSOSP perturbations in the PSP2 configuration. Conceptually, however, schemes for different subgrid processes could easily be combined.

In the introduction we referred to the open question of whether small-scale orography matters for convective initiation. Here we tested this by modeling subgrid mechanical lifting on the smallest resolved scale. Our results suggest that the subgrid orographic triggering effect most likely is not that relevant. On the one hand, to see the desired effect in mountainous regions we had to increase the perturbation strength to a point where unwanted effects occurred over flat terrain. On the other hand, the perturbations introduced by our scheme are active on the smallest resolved scales of the kilometer-scale model. Actual small-scale orography perturbations would be even smaller, and therefore likely even less influential. In general, our results confirm in a realistic setup what has been observed in other studies (Langhans et al. 2013; Schneider et al. 2018; Tucker and Crook 2005; Kirshbaum et al. 2007a): conceptually, SSO may be relevant in regions of insignificant resolved orography; yet due to the fractal, scale-invariant structure of orography (Turcotte 1987) regions of high subgrid orographic variations are generally also regions of high grid-scale orography. This larger scale, resolved orography may dominate the triggering of convection so that SSO is virtually irrelevant for convective initiation.

## Acknowledgments

The research leading to these results has been done within the subproject “Representing forecast uncertainty using stochastic physical parameterizations” of the Transregional Collaborative Research Center SFB/TRR 165 (www.wavestoweather.de) funded by the German Research Foundation (DFG). We further acknowledge K. Kober and F. Brundke for developing a previous version of the SSOSP scheme that was used as baseline for the final SSOSP version.

### APPENDIX A

#### SSOSP: Detailed Description

##### a. Scaling orographically induced vertical velocity

The orography acts as a boundary and creates internal gravity waves, which can be evanescent or propagating. The orographic height can be idealized as a two-dimensional sine wave with amplitude *h*_{0} and horizontal wavenumbers *k* and *l* (or a composition of many such sine waves). Using a general exponential ansatz for vertical velocity *w* as a solution to the Taylor–Goldstein equation and *ω* as oscillation frequency, the following two cases can be distinguished (Gill 1982): For the propagating case (i.e., *ω*^{2} ≤ *N*^{2}), the solution for *w* reads as

When *ω*^{2} > *N*^{2}, the amplitude of the vertical displacement decays with height:

with *γ*^{2} = −*m*^{2}.

For both cases, *ω* = *ku* + *lυ* and *w*_{0} and *m* are given as

where *u*, *υ* denote *x* and *y* components of the horizontal wind. Further details can be found in Gill (1982).

##### b. Parameterization

For the parameterization, we are not interested in representing the complete wave structure. Instead we aim at representing the effect of SSO-induced gravity waves on the resolved scales in a stochastic manner by scaling the random perturbations accordingly. *w*′ of Eq. (6) is hence linked to the amplitude of gravity waves triggered by SSO [see Eqs. (A1) and (A2)]. The cosine terms in Eqs. (A1) and (A2) are consequently neglected. We formulate *w*_{0} and *m* in terms of SSO information. Following Baines and Palmer (1990), the subgrid-scale orography for one grid box of the COSMO model is represented by four parameters, namely *θ*_{sso}, *γ*_{sso}, *σ*_{sso}, and *μ*_{sso}. They describe, respectively, the orientation, shape, slope and standard deviation of height of the orography. We specify the velocities *u*, *υ* as *u*′, *υ*′ in the rotated system of the SSO:

The horizontal wavenumbers *k* and *l* (in the rotated system of the SSO) are obtained by using the semiminor and semimajor axes *a*, *b* for the elliptical mountains that result when the orography is represented by sine waves:

(Lott and Miller 1997). We set *h*_{0} to be *μ*_{sso}, implying that higher mountains will result in larger vertical displacements.

Using these approximations, we can formulate the vertical wavenumber *m* and wave amplitude *w*_{0}, given by Eqs. (A4) and (A3), in the following way:

### APPENDIX B

#### Computational Setup and Details

##### a. Computational details

Analysis and visualization was done using Python. The code and data are available from the authors upon request. For reading and processing the COSMO model output files the ens_tools package, which is based on the xarray library (Hoyer and Hamman 2017), was used which is being developed within the *Waves to Weather* project and will be available publicly soon.

##### b. SAL computation

We compute the SAL score using the threshold *R** for identifying precipitation objects by using radar observations as reference; *R** is computed based on the suggestion by Wernli et al. (2009): *R** = (1/15)*R*^{95}, where *R*^{95} corresponds to the 95th percentile of precipitation of grid points exceeding 0.1 mm h^{−1}; *R*^{95} is computed separately for radar and forecast and for each time step. Different threshold choices, such as fixed precipitation amounts, did not change the outcome qualitatively.

## REFERENCES

*Statistical Methods in the Atmospheric Sciences*. Academic Press, 648 pp.

## Footnotes

This article is included in the Waves to Weather (W2W) Special Collection.

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

This threshold approximately coincides with the median of SSO-std height.

^{2}

Separating southern and northern Germany or the Alpine region instead of single grid points showed similar results.