Reduced non-Gaussianity by 30 s rapid update in convective-scale numerical weather prediction

Non-Gaussian forecast error is a challenge for ensemble-based data assimilation (DA), particularly for more nonlinear convective dynamics. In this study, we investigate the degree of the non-Gaussianity of forecast error distributions at 1 km resolution using a 1000-member ensemble Kalman filter, and how it is affected by the DA update frequency and observation number. Regional numerical weather prediction experiments are performed with the SCALE (Scalable Computing for Advanced Library and Environment) model and the LETKF (local ensemble transform Kalman filter) assimilating phased array radar observations every 30 s. The results show that non-Gaussianity develops rapidly within convective clouds and is sensitive to the DA frequency and the number of assimilated observations. The non-Gaussianity is reduced by up to 40 % when the assimilation window is shortened from 5 min to 30 s, particularly for vertical velocity and radar reflectivity.

1 Introduction

The Kalman filter (KF) is the minimum variance linear unbiased estimator of the state of a dynamical system. The ensemble Kalman filter (EnKF; Evensen, 2009; Houtekamer and Zhang, 2016) is a Monte Carlo extension to the KF suitable for nonlinear systems with a large number of variables, so it became a viable choice for data assimilation (DA) in numerical weather prediction (NWP) and other geoscience applications. The EnKF is optimal in the sense of maximum likelihood estimation when the error distributions are Gaussian (Evensen, 2009), but it becomes suboptimal when the observational and forecast error distributions depart from the Gaussian (Lei et al., 2010). Miyoshi et al. (2014, 2015) and Kondo and Miyoshi (2019) investigated non-Gaussianity in forecast error distributions using a 10 240-member EnKF with global atmospheric models. They showed that large non-Gaussianity measured by the Kullback–Leibler divergence is found frequently in the tropics, mainly due to abundance of deep moist convection, and also in other active areas with a real-world NWP model at relatively low 112 km resolution. In those experiments, temperature and moisture show generally more non-Gaussian distributions than winds. Recently, the horizontal resolution of operational NWP systems reached the order of 1 km, which is fine enough to resolve convective phenomena explicitly. Obtaining appropriate initial conditions at such a high resolution is a challenge (Sun et al., 2014). The EnKF has been successfully applied to the mesoscale assimilation of radar and satellite data (e.g., Stensrud et al., 2013). However, previous studies (Jacques and Zawadzki, 2014; Kawabata and Ueno, 2020) revealed that the underlying assumptions, such as linear error dynamics and Gaussian error distributions, are much more questionable in mesoscale than in synoptic and larger scales.

Miyoshi et al. (2016a, b) developed a so-called big data assimilation (BDA) system assimilating observations every 30 s at 100 m resolution, taking advantage of new-generation technologies, like the phased array weather radar (PAWR), which provide observations at unprecedented high temporal and spatial resolution. With the BDA configuration under an idealized observing system simulation experiment (OSSE) framework, Maejima and Miyoshi (2020) showed that DA cycles every 1 min resulted in better analyses than cycles every 15 min. However, the impact of the DA frequency upon the forecast error distribution has not been investigated in real-case convective scale NWP.

This study investigates how the DA frequencies affect non-Gaussianity using a 1000-member and 1 km mesh EnKF. The 1000 ensemble members would be useful for detecting non-Gaussian forecast error distributions, as suggested by Kondo and Miyoshi (2019). Necker et al. (2020a, b) performed similar experiments and investigated the covariance structure and the effect of sampling noise at the mesoscale in a heavy rainfall event over Germany. Although the previous research employed data assimilation with only conventional observations at a 3 h DA frequency, this study is fundamentally different in the convection-resolving rapid DA cycles with PAWR data as frequently as every 30 s. The high-frequency data allow us to investigate the sources of non-Gaussian distributions at the kilometer scale in the presence of rapidly evolving deep moist convection. The paper is organized as follows: Sect. 2 describes the methodological aspects. Results are presented in Sect. 3, and Sect. 4 provides concluding remarks and discussion.

2 Methodology

We use observations from the PAWR at Osaka University, Suita, Japan (Yoshikawa et al., 2013, Fig. 1a, red cross). This PAWR provides a unique data set suitable for this study, with various assimilation frequencies up to every 10 s at the fastest. This study follows the case study of Miyoshi et al. (2016a), focusing on the period between 04:00 and 06:00 UTC, on 13 July 2013, when heavy rains produced flash floods in Kyoto. Individual convective cells moved from west to east within a quasi-stationary intense rainband (see Fig. 1b for a snapshot at 05:30 UTC). For this period, full volume scans of the PAWR are available every 30 s, with 98 elevation angles, an azimuthal resolution of 1.2^∘, and a range resolution of 100 m up to a maximum range of 60 km (Fig. 1a, red circle). Unambiguous Doppler velocities are available in the range −50 to 50 ms⁻¹. PAWR reflectivity data are quality-controlled, following Ruiz et al. (2015). A simple quality control algorithm has also been applied to the Doppler velocity field to remove outliers.

Figure 1 (a) Terrain height of the 1 km mesh SCALE–LETKF domain (shades; meters). The red circle indicates the 60 km radar range centered at the radar site (red cross) at Osaka University, Suita, Japan. The black box indicates the area shown in panels (b)–(d). (b) Column maximum PAWR observation (decibel relative to Z; hereafter dBZ) at 05:30 UTC, half an hour after the initialization of the data assimilation cycle. (c) 5MIN and (d) 30SEC experiments analysis ensemble mean column maximum radar reflectivity (dBZ) at 05:30 UTC. Black lines indicate the locations of the cross sections displayed in Fig. 2

In this study, the regional NWP model known as the Scalable Computing for Advanced Library and Environment model (SCALE; Nishizawa et al., 2015) is used, coupled with the local ensemble transform Kalman filter (LETKF; Hunt et al., 2007). Lien et al. (2017) and Honda et al. (2018) describe the SCALE–LETKF system in detail. The model configuration follows Lien et al. (2017), with a single moment bulk microphysics scheme (Tomita, 2008), a level 2.5 boundary layer turbulence scheme (Nakanishi and Niino, 2004), the model simulation radiation transfer radiation scheme (Sekiguchi and Nakajima, 2008), and soil processes represented by a Beljaars-type soil model (Beljaars and Holtslag, 1991).

The SCALE–LETKF system is implemented over a single domain, with a horizontal resolution of 1 km and a size of 180 km by 180 km (Fig. 1a). There are 50 vertical levels extending up to 18 km elevation, with a variable grid spacing from 140 to 790 m in a hybrid sigma z terrain-following coordinate. A 1000-member ensemble is used to assimilate the observations. Kondo and Miyoshi (2019) showed significant sampling error contaminations in non-Gaussian measures when the ensemble size is smaller than 1000. The initial conditions for the first cycle and the boundary conditions are taken from the National Centers for Environmental Prediction Global Data Assimilation System final analysis (FNL). Using FNL as the boundary conditions may be overly optimistic for the forecasting purpose, but this is not relevant to the goal of this study, which focuses on non-Gaussian distributions and the impact of DA frequency.

The initial ensemble at the first assimilation cycle and the boundary condition ensemble are created by adding random perturbations which preserve the hydrostatic and nearly geostrophic equilibrium (Necker et al., 2020a; Maldonado et al., 2021). These perturbations are generated from a sample of continuous 6 h analysis states provided by the Climate Forecast System Reanalysis (CFSR; Saha et al., 2010), $\left[X_{\mathrm{CFSR}}\right(t_{\mathrm{1}}),X_{\mathrm{CFSR}}(t_{\mathrm{2}}),\mathrm{\dots },X_{\mathrm{CFSR}}(t_N\left)\right]$, where N=5840 (4 years). The horizontal grid spacing of the CFSR data is 0.5^∘. At the beginning of the assimilation cycle (t=t_s), the initial condition perturbation of the ith member X^′(i) is computed as follows:

$$X^{\prime \left(i\right)}\left(t_s\right)=\mathit{\alpha }\left[X_{\mathrm{CFSR}}\right(t_{{n_{\mathrm{1}}}^{\left(i\right)}})-X_{\mathrm{CFSR}}(t_{{n_{\mathrm{2}}}^{\left(i\right)}}\left)\right],$$

where α is a multiplicative factor equal to 0.1 so that the amplitude of the perturbations is roughly equivalent to 10 % of the climatological variability. The two CFSR analysis states are chosen by randomly selecting two numbers $n_{\mathrm{1}}^{\left(i\right)}$ and $n_{\mathrm{2}}^{\left(i\right)}$ from the N elements satisfying the condition that $t_{{n_{\mathrm{1}}}^{\left(i\right)}}$ and $t_{{n_{\mathrm{2}}}^{\left(i\right)}}$ correspond to the same time of the year and time of the day. In the following assimilation cycles at time t>t_s, we obtain the boundary perturbations as:

$$\begin{array}{rl}{X}^{\prime \left(i\right)}\left(t\right)& =\mathit{\alpha }\left[\right(\mathrm{1}-\mathit{\beta }\left)\right(X\left(t_{{l_{\mathrm{1}}}^{\left(i\right)}}\right)-X\left(t_{{l_{\mathrm{2}}}^{\left(i\right)}}\right))\\ & +\mathit{\beta }\left(X\right(t_{{u_{\mathrm{1}}}^{\left(i\right)}})-X(t_{{u_{\mathrm{1}}}^{\left(i\right)}}\left)\right],\end{array}$$

where $l_{\mathrm{1},\mathrm{2}}^{\left(i\right)}=n_{\mathrm{1},\mathrm{2}}^{\left(i\right)}+m$ and $u_{\mathrm{1},\mathrm{2}}^{\left(i\right)}=n_{\mathrm{1},\mathrm{2}}^{\left(i\right)}+m+\mathrm{1}$, with $m=\mathrm{floor}\left[\right(t-t_s)/\mathrm{6}\mathrm{h}]$ and $\mathit{\beta }=\left[\right(t-t_s)/\mathrm{6}\mathrm{h}]-m$ being a temporal linear interpolation factor to compute perturbations at arbitrary times (not necessarily a multiple of 6 h). In this way, we obtain perturbations that are smoothly varying in time and consistent with the large-scale dynamics of the atmosphere. This procedure is applied to all atmospheric and soil state variables.

In the SCALE–LETKF system, radar data can be assimilated using different localization scales for different variables. Based on preliminary experiments with the SCALE–LETKF using smaller ensemble sizes and PAWR data every 30 s, it was found that a vertical localization scale of 2 km (with a 7.3 km or similar cut-off hereafter) produced good results. For horizontal localization, better results were obtained using 4 km localization to assimilate observations with reflectivities >10 dBZ. Observations of reflectivity values ≤10 dBZ are assimilated with a fixed value of 10 dBZ to avoid large observation minus forecast departures associated with clear air reflectivities (Aksoy et al., 2009). Also, a shorter horizontal localization scale of 2 km is used to reduce the impact of non-precipitating observations at the edge of clouds. Doppler velocity observations are assimilated with horizontal and vertical localization scales of 10 and 3 km, respectively. For covariance inflation, a relaxation to prior ensemble spread (RTPS,  Whitaker and Hamill, 2012) with a relaxation parameter of 0.9 is applied. This helps consider the inhomogeneous distribution of observations as in Lien et al. (2017).

Reflectivity and Doppler velocity observations are superobbed to a horizontal resolution of 1 km and a vertical resolution of 500 m to approximately match the model resolution. This helps reduce the errors of representativeness due to the gap between what is represented by the model and observation. This procedure can also reduce the impact of possible spatial correlations in the observation errors. The observational error standard deviations for these super-observations are set at 5.0 dBZ and 3.0 ms⁻¹ for reflectivity and Doppler velocities, respectively. The radar data are assimilated up to a maximum height of 11 km.

A spin-up DA experiment with every 5 min PAWR reflectivity and Doppler velocity data is performed for an hour from 04:00 UTC on 13 July 2013. Only a single PAWR volume scan closest to the analysis time is assimilated per analysis. The 1000-member analysis ensemble at 05:00 UTC is used as the initial conditions for the DA experiments.

Experiments are performed with different DA update frequencies to study the impacts of the DA frequency and observation number on the forecast error distributions. All experiments share the configuration described above, but the only differences are the DA frequency and the amount of data assimilated. First, four experiments with 5, 2, 1, and 0.5 min DA frequencies are performed (hereafter referred to as 5MIN, 2MIN, 1MIN, and 30SEC, respectively). Here, only a single-volume scan closest to the analysis time is used per analysis. That is, more frequent updates assimilate more data. In all cases, the time difference between the observation time (center time of the radar volume scan) and the analysis time do not differ by more than 15 s.

Next, to separate the impact of DA frequency and the amount of data assimilated, two additional experiments are performed using a 5 and 1 min DA frequency, with all radar volumes every 30 s assimilated by a 4-dimensional EnKF approach (Hunt et al., 2004). These experiments are referred to as 5MIN-4D and 1MIN-4D, respectively, assimilating the same amount of data as 30SEC but using longer assimilation windows.

To measure the degree of non-Gaussianity of the error distributions, we compute the Kullback–Leibler divergence (hereafter KLD; Kullback and Leibler, 1951), which is defined as follows:

$$\begin{array}{}\text{(1)}& \mathrm{KLD}\left(P\right|\left|Q\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}p\left(x\right)\ln {\displaystyle \frac{p\left(x\right)}{q\left(x\right)}}\mathrm{d}x,\end{array}$$

where p(x) and q(x) are the probability density functions (PDFs) of P and Q, respectively. The KLD is 0 if P and Q are the same and takes positive values if P and Q differ. In our case, p(x) is either the first guess or analysis error distribution for the state variable x, and q(x) is a Gaussian distribution whose mean and standard deviation are equal to the ones of p(x). Therefore, a low KLD value corresponds to the first guess or analysis error distribution close to a Gaussian. In the EnKF, we do not have access to the continuous PDF p(x) but to its finite, limited sample. For each state variable x (e.g., temperature, wind components, etc.), and at each model grid point, we approximate p(x) with the sample histogram from the 1000-member ensemble using 32 equally sized bins covering the range where p(x) is greater than 0. This range is defined by the minimum and maximum values of x at each model grid point and time. Hence, we can approximate the KLD as follows:

$$\begin{array}{}\text{(2)}& \mathrm{KLD}\left(P\right|\left|Q\right)\approx \sum _{j=\mathrm{1}}^{j=\mathrm{32}}{p}_j\ln {\displaystyle \frac{p_j}{q_j}},\end{array}$$

where p_j is the empirical frequency of x at the jth histogram bin. q_j is the integral over the jth histogram bin of a Gaussian PDF, whose mean and standard deviation are given by the ensemble-based sample estimates. After implementing this, we end up with an estimation of the KLD of the analysis and first-guess error distributions with respect to the Gaussian for each grid point location, vertical level, and time.

3 Results

All experiments show that the analyzed reflectivity fields are in good agreement with the observation. However, some differences can be found between the experiments that assimilate different amounts of data and with different assimilation windows. For example, Fig. 1c and d show that 30SEC captures the strong reflectivity areas (>45 dBZ; orange and red shadings) better than 5MIN. 5MIN shows noisy patterns of spurious convective cells surrounding the main convective rainband.

3.1 Impact of DA frequency on the analysis of a convective cell

In this section, we explore the impact of data assimilation frequency over a convective cell located in the cross section along the black line in Fig. 1b–d at 05:30 UTC. First, the impact of the data assimilation frequency is explored by the 5MIN, 2MIN, 1MIN, and 30SEC experiments. Here, more observations are assimilated with more frequent data assimilation. Figure 2 top row (a–d) shows that the reflectivity patterns (Z; shades) are similar among all experiments, but the vertical velocity (W; contours) is different. Stronger updrafts are found in DA experiments with shorter assimilation windows. This suggests that DA frequency have a significant impact upon quantities which are not directly observed.

Figure 2 (a–h) South–north vertical cross section along the black line indicated in Fig. 1b–d at 05:30 UTC for (a–d) first-guess ensemble mean reflectivity (Z; shades; dBZ) and vertical velocity (W; contours every 2.5 m s⁻¹). (e–h) Vertical velocity KLD (shades; 10⁻²) and ensemble spread (red contours at 1.0, 2.5, 5.0, and 10.0 m s⁻¹). (i–l) Temperature KLD (shades; 10⁻²) and ensemble spread (red contours at 0.2, 0.5, 1.0, and 2.0 K). Blacked dashed contours indicate reflectivity over 30 dBZ. The black cross in panels (i)–(l) indicates the location of the maximum KLD within the grid points at which Z>30 dBZ.

Strong non-Gaussianity is observed in the first-guess ensemble in W and temperature T in the 5MIN experiment (Fig. 2e and i, respectively). Non-Gaussianity for W is stronger at the southern edge and the highest peak of the convective cell, which is probably related to the development of a new updraft in the southern edge and the top of the strong updraft, respectively. Weaker low-level maxima south of the convective line are associated with shallow convective clouds that are not effectively corrected by radar observations. The KLD maxima for T are approximately collocated as those for W. KLD maxima in T can be associated with non-Gaussianity in W through vertical advection of scalar quantities such as T and moisture. Another KLD maximum for T is found near the surface south of the convective cell, probably associated with the gust front.

Kondo and Miyoshi (2019) found that, in synoptic scales, the ensemble spread maxima are collocated with the KLD maxima. At convective scales for W, the ensemble spread maxima (Fig. 2e; red contours) are not necessarily collocated. For example, larger departures from the Gaussian are found above the ensemble spread maximum associated with the main updraft in the 5MIN experiment. For temperature, there is also no clear relation in the distribution of the ensemble spread and the KLD, although KLD maxima seem to occur within areas of a relatively large ensemble spread. As the assimilation frequency increases, it is more difficult to find a relationship between KLD and ensemble spread for either W or T (Fig. 2; second and third rows).

KLD for W and T are consistently reduced with more frequent DA (Fig. 2e–h), although the reduction is smaller for T. Overall, KLD is reduced more from 5MIN to 2MIN than from 1MIN to 30SEC. This reduction occurs mainly within the convective clouds. Non-Gaussianity in W at low levels observed outside the cloud is not significantly affected by more frequent updates. The ensemble spread for W is also reduced with more frequent DA and indicates a narrower error distribution. This result is linked with the reduced non-Gaussianity since it is expected that smaller perturbations grow in a more linear regime and contribute to reducing the departures from the Gaussian.

To better investigate the shape of the error distributions and how they are affected by the update frequency, Fig. 3 shows the sample histograms for the first guess at the location of maximum KLD (indicated with a black cross in Fig. 2). We restrict the search of the maximum KLD to the grid points at which the ensemble mean reflectivity is over 30 dBZ, where radar data impact would be large. The forecast error distribution for W and for the 5MIN experiment shows large departures from the Gaussian with a strong positive tail (Fig. 3a). A similar situation is observed for Z (Fig. 3e). This result is consistent since ensemble members with larger W are probably associated with larger reflectivity values, so both distributions become positively skewed. As the update frequency is increased, non-Gaussianity and ensemble spread are reduced for both variables. The only exception is that Z at 30SEC update frequency shows a KLD value that is slightly larger than that in the 1MIN experiment. Note that these error distributions are taken at slightly different locations based on the simulated convection locations in each experiment, and thus, the mean of the distribution can change from one experiment to the other.

Figure 3 Sample histograms for (a–d) vertical velocity (m s⁻¹ – meters per second) and (e–h) reflectivity (dBZ) at the location of the maximum KLD for vertical velocity (black cross in Fig. 2). Thick dashed curves indicate fitted Gaussian functions, and KLD non-Gaussian measures are also indicated. Each column corresponds to 5MIN, 2MIN, 1MIN, and 30SEC from left to right, respectively.

4D EnKF experiments allow us to investigate the impact of changing the assimilation frequency while keeping the observation number unchanged. 5MIN-4D shows weaker updrafts (similar to those found in 5MIN) compared with experiments with more frequent updates (Fig. 4a, b). 5MIN-4D also shows almost the same ensemble spread for W and T as 5MIN (Fig. 4c and e; red contours). KLD for W (Fig. 4c; shading) is lower, indicating that the observation number contributes to reducing non-Gaussianity. This is not the case for T for which KLD is similar or larger (Fig. 4e; shading). 1MIN-4D is close to 1MIN and 30SEC in terms of non-Gaussianity and the shape and strength of the convective cell (Fig. 4b, d and f).

Figure 4 As in Fig. 2 but for (a, c, e) 5MIN-4D and (b, d, f) 1MIN-4D.

3.2 Spatiotemporal distribution of non-Gaussianity

We further investigate the non-Gaussianity by averaging the KLD vertically and temporally (Fig. 5). In 5MIN, the central and eastern sides of the convective area show relatively low KLD values because the impact of radar DA is generally bigger in the convective areas (Fig. 5a). The impact of DA frequency on non-Gaussianity is investigated by means of the relative KLD difference between the 5MIN and all the other experiments, which is computed as follows:

$$\begin{array}{}\text{(3)}& \stackrel{\mathrm{‾}}{{\mathrm{KLD}}_{\mathrm{diff}}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathrm{KLD}}_{\mathrm{E}}}-\stackrel{\mathrm{‾}}{{\mathrm{KLD}}_{\mathrm{5}\mathrm{MIN}}}}{\stackrel{\mathrm{‾}}{{\mathrm{KLD}}_{\mathrm{5}\mathrm{MIN}}}}},\end{array}$$

where $\stackrel{\mathrm{‾}}{{\mathrm{KLD}}_{\mathrm{diff}}}$ is the relative difference between the averaged KLD in the 5MIN experiment ($\stackrel{\mathrm{‾}}{{\mathrm{KLD}}_{\mathrm{5}\mathrm{MIN}}}$) and on each of the other experiments ($\stackrel{\mathrm{‾}}{{\mathrm{KLD}}_{\mathrm{E}}}$), where E can be either 5MIN-4D, 2MIN, 1MIN, 1MIN-4D or 30SEC).

Figure 5 (a) Column-averaged KLD for zonal wind for 5MIN, averaged for the experiment period from 05:00 to 06:00 UTC. Relative KLD difference (percent) from 5MIN for (b) 2MIN, (c) 1MIN, (d) 30SEC, (e) 5MIN-4D, and (f) 1MIN-4D. Warm colors correspond to smaller KLD values.

KLD consistently decreases with increasing DA frequency (Fig. 5b–d). KLD is reduced by up to 40 % in 30SEC with respect to the 5MIN. KLD is reduced more in the convective area, where more observations are assimilated. Increasing the DA frequency and the observation number produces a more substantial impact over the western part of the convective line where KLD maxima are found associated with convective cells entering the radar range from the west.

KLD in 1MIN-4D is as low as that in 30SEC and lower than that in 1MIN. This result suggests that both observation number and DA frequency contribute to reducing non-Gaussianity, at least for high DA frequencies. KLD in 5MIN-4D is lower than that in 5MIN, so that a larger observation number helps to reduce non-Gaussianity. However, KLD in 5MIN-4D is larger than that in 30SEC or 1MIN-4D, indicating that the DA frequency is equally important. Moreover, the impact of DA frequency can be larger in the case of variables like T and moisture. As already found in the vertical cross sections (Fig. 4), for those variables, KLD in 5MIN and 5MIN-4D is almost the same, while KLD is clearly reduced for 1MIN, 1MIN-4D, and 30SEC (not shown).

We also investigate the vertical distribution of non-Gaussianity by the spatially averaged vertical profile of KLD at “precipitating” grid points, defined by the ensemble mean column maximum reflectivity >30 dBZ, and “non-precipitating” grid points, defined by the ensemble mean column maximum reflectivity <0 dBZ. At the precipitating grid points (Fig. 6a–d), KLD, for temperature and vertical velocity, is maximum at mid levels coinciding with the maximum in latent heat release within convective clouds and with the maximum ensemble spread for these two variables (not shown). KLD or temperature, vertical velocity, and specific humidity maximizes at lower heights over the non-precipitating area, since, as stated before, at such locations non-Gaussianity is mainly associated with shallow convection. For instance, for the vertical velocity, the ensemble spread in the shallow convection is usually low, but the KLD can be larger. An upper-level maximum in KLD is found for the meridional wind (Fig. 6d and h), which also coincides with the maximum ensemble spread (not shown). Convective outflows are stronger at the top of convective clouds and can be one of the mechanisms contributing to the increase in non-Gaussianity at these levels over the precipitating area. Overall, KLD in 30SEC is lower than that in 5MIN with reductions of more than 40 %. The reduction in KLD in the non-precipitating area is smaller because the radar DA is inherently less effective in these areas (Fig. 6e–h). There are some exceptions to the general reduction in non-Gaussianity with increased update frequency. Specific humidity in non-precipitating grid points shows larger KLD in the 30SEC than in the 5MIN experiments. This is also the case for the precipitating grid points at upper levels in the second half of the experiment. Also, the KLD in W in the non-precipitating grid points at middle and upper levels is slightly larger in the 30SEC experiment.

Figure 6 Time vertical cross section of KLD for 5MIN (contours; 10⁻²) and the relative difference for 30SEC (shaded), averaged over the (a–d) precipitating (>30 dBZ) and (e–h) non-precipitating (<0 dBZ) grid points for (a, e) temperature, (b, f) specific humidity, (c, g) vertical velocity, and (d, h) meridional wind.

3.3 Non-Gaussianity evolution within the DA cycle

To investigate the effect of the analysis update on non-Gaussianity, we present the time series of the KLD of the analysis and first guess vertically and horizontally averaged over the precipitating and non-precipitating grid points (Fig. 7). At most times and variables over the precipitating and non-precipitating grid points, KLD is reduced during the assimilation step. Experiments with longer windows show more KLD growth during the forecast as expected but also a larger reduction at the analysis step, which is not as effective as the more frequent updates in reducing the analysis KLD. As noted before, the specific humidity over the non-precipitating grid points behaves differently, and KLD increases during the assimilation step for almost all times and experiments, leading to larger KLD at shorter assimilation windows (Fig. 6b and f). In this area, mostly non-precipitating observations are assimilated to suppress spurious clouds. Interestingly, in the non-precipitating grid points, 5MIN-4D is the experiment providing the lowest KLD for all variables (Fig. 7b, d, and f). This result suggests the potential benefits of treating non-precipitating observations differently.

Figure 7 Sawtooth time series of the KLD (10⁻²) of the analysis and first guess vertically and horizontally averaged over the precipitating (>30 dBZ; a, c, e) and non-precipitating (<0 dBZ; b, d, f) grid points for temperature (a, b), specific humidity (c, d), and vertical velocity (e, f) and for the 5MIN (red), 5MIN-4D (blue), 2MIN (green), 1MIN (magenta), 1MIN-4D (black), and 30SEC (cyan) experiments.

To evaluate the impact of assimilation frequency on the distance between the analysis and first guess to the observations in a more systematic way, we compute the root mean squared error (RMSE) and bias for reflectivity observations (Fig. 8). The computation of the RMSE and bias between the model and the observations is done by comparing the column maximum of the reflectivity for each horizontal grid location and time. The RMSE and bias are computed only over grid points at which the observed maximum reflectivity is over 5 dBZ. The time series of RMSE shows a better fit to the observed reflectivity for shorter assimilation windows. The impact of 4D DA is not so clear; 1MIN-4D slightly outperforms the 1MIN, but 5MIN-4D and 5MIN perform similarly (Fig. 8a). This is partially because in 4D data assimilation the analysis results from the assimilation of all the observations within the assimilation window, while, to construct this figure, only the observations at the analysis time were considered. The bias, computed as the mean difference between the model and the observations does not seem to be consistently affected by the assimilation frequency (Fig. 8b). These results are in agreement with those observed in the time series of KLD for different variables. However, we should be cautious with the interpretation of these results since increasing the observation number can lead to both a reduced KLD and a better fit to the observed quantities, which does not necessarily imply a causal link between these two effects.

Figure 8 Sawtooth time series of the root mean squared error (dBZ; a) and bias (dBZ; b) of the maximum reflectivity of the analysis and first guess for the 5MIN (red), 5MIN-4D (blue), 2MIN (green), 1MIN (magenta), 1MIN-4D (black), and 30SEC (cyan) experiments.

4 Summary and discussion

We performed 1000-member 1 km resolution ensemble DA experiments using real phased array radar observations and a mesoscale NWP model to investigate the impact of DA frequency and observation number on the non-Gaussian error distributions. We found that a DA frequency of 5 min, although it was already much faster than the typical DA frequency, resulted in strong non-Gaussianity possibly affecting the performance of the EnKF. Non-Gaussianity is stronger for vertical velocity, as has been found by Kawabata and Ueno (2020)​​​​​​​. Non-Gaussianity is also larger at mid levels within convective cells, near the level of larger latent heat release and vertical accelerations associated with convective instability. At convective scales, some of the local maxima in KLD can be related directly to advection by mesoscale circulations associated with strong convective cells, but other processes not specifically presented in this study may also possibly contribute to the generation of non-Gaussianity, such as those not directly associated with clouds, like differential heating circulations or gravity waves.

We found that increasing the analysis update frequency and observation number from 5 min to 30 s has a huge impact upon non-Gaussianity in the error distributions for all model variables but particularly for vertical velocity and reflectivity, which are the ones showing larger KLD from Gaussianity at these scales. Increasing the assimilation frequency to 30 s and assimilating more observations can reduce KLD by up to 40 %. Moreover, 4D EnKF experiments revealed that, for frequent DA of every 1 min, the observation number explained most of the reduction in non-Gaussianity; in contrast, for a longer window of 5 min, even the experiments using all 30 s frequency observations present significant departures from the Gaussian. While convective clouds are particularly favorable for nonlinear error growth, non-Gaussianity is not necessarily larger within convective clouds. This is mainly due to the convective-scale radar that DA is usually most effective within precipitating clouds.

There are two possible ways in which more frequent DA can result in error distributions closer to the Gaussian. First, more frequent DA contributes to a quasi-linear evolution of the forecast error due to forecast lengths which are shorter than the predictability limit for the resolved scales. This also helps to keep the perturbation small and can, additionally, contribute to quasi-linear perturbation dynamics. Second, our results show that the analysis step effectively contributes to reducing non-Gaussianity for different variables, although this may not be the case for non-precipitating reflectivity observations that produce an increase in KLD for specific humidity. Non-Gaussianity reduction during DA is larger with longer windows. However, it is not enough to compensate for the effect of more rapid and nonlinear error growth during the forecast step in the lower update frequency experiments.

From the point of view of KLD reduction, the largest impact is found between 5MIN and 2MIN updates. This suggests that nonlinear error growth becomes more important after the first 2 min of integration at these scales. This hypothesis is partially supported by the reduction in RMSE and ensemble spread. A 2 min update frequency seems to provide a good compromise between the computational cost and non-Gaussianity of the error distributions. However, from the point of view of the analysis accuracy, more frequent DA provides a better fit to the observed quantities. The specific role of reduced non-Gaussianity on this is not clear and should be further investigated. Gaussian error distributions may contribute to more accurate analysis updates, but in the current experimental setting, other factors like the increase in the number of assimilated observations may also lead to the reduction in the RMSE for observed quantities. Maejima and Miyoshi (2020) investigated the impact of assimilation frequency at 1 km using observing system simulation experiments. They also found a significant improvement in the forecast quality when the assimilation window is reduced from 5 to 3 min and additional improvements using 1 min windows. These results are consistent with what is found in this paper with respect to Gaussianity in the error distributions.

Moreover, as has been shown in the previous studies, more frequent assimilation can produce a larger degree of imbalance in the initial conditions which can degrade the quality of the forecasts (e.g., Lange and Craig, 2014; Bick et al., 2016). Therefore, despite the potential benefits of a more Gaussian model error distribution on the analysis accuracy, other factors may degrade the forecasts initialized from more frequent data assimilation cycles. Imbalance may also be an additional source of non-Gaussianity. Gaussian error distributions can lead to more physically meaningful assimilation updates in the context of an EnKF and, therefore, more balanced initial conditions. However, a larger imbalance in the initial conditions can contribute to faster error growth and increased departure from the Gaussian in the forecast distribution. Possible interactions of these mechanisms in a data assimilation cycle have not been investigated and are a subject for future research. Our results suggest that, despite the effect of a larger imbalance, the increase in DA frequency reduces non-Gaussianity in the sample distributions with the EnKF. This is even true of variables like vertical velocity within convective clouds which are frequently used to measure the effect of imbalance in the initial conditions.

This study is the first attempt to investigate the impact of assimilation frequency and observation number on non-Gaussianity using an EnKF employing a large 1000-member ensemble and every 30 s observations from a PAWR. In this first set of experiments, we evaluate the impact on the non-Gaussianity of the ensemble-based sample distribution. Future experiments will be performed to investigate the overall quality of the analysis obtained with different assimilation windows and the number of observations and also the impact of the assimilation window on the structure of the error covariance matrix.