Application of multigrid NLS-4DVar in radar radial velocity data assimilation with WRF-ARW

ZHANG Lu,TIAN Xiangjun,c and ZHANG Hongqin

aInternational Center for Climate and Environment Sciences,Institute of Atmospheric Physics,Chinese Academy of Sciences,Beijing,China;

bCollege of Earth and Planetary Sciences,University of Chinese Academy of Sciences,Beijing,China; cCollaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters,Nanjing University of Information Science and Technology,Nanjing,China

ABSTRACT The nonlinear least-squares four-dimensional variational assimilation (NLS-4DVar) method introduced here combines the merits of the ensemble Kalman filter and 4DVar assimilation methods.The multigrid NLS-4DVar method can be implemented without adjoint models and also corrects small-to large-scale errors with greater accuracy.In this paper,the multigrid NLS-4DVar method is used in radar radial velocity data assimilations. Observing system simulation experiments were conducted to determine the capability and efficiency of multigrid NLS-4DVar for assimilating radar radial velocity with WRF-ARW(the Advanced Research Weather Research and Forecasting model).The results show significant improvement in 24-h cumulative precipitation prediction due to improved initial conditions after assimilating the radar radial velocity. Additionally, the multigrid NLS-4DVar method reduces computational cost.

KEYWORDS Heavy rainfall;multigrid scheme NLS-4DVar method;radar radial velocity data assimilation

1. Introduction

Radar data assimilation strongly impacts short-term weather forecasting;however,unconventionally observed parameters and mismatched resolutions between observations and models limit assimilation implementation(Fabry and Kilambi 2011; Thompson, Wicker, and Wang 2012).Recently,four-dimensional variational assimilation(4DVar)has been applied to radar data to improve quantitative precipitation forecasts of mesoscale convective systems(Wang et al. 2013; Sun and Wang 2013). However, traditional 4DVar needs adjoint models;thus,data assimilation is complicated by a series of challenges in coding,maintaining and updating the adjoint model.On the other hand,the ensemble Kalman filter(EnKF)method has attracted broad attention because of its simple conceptual formulation and relative ease with which it can be implemented;as such,it is widely used in radar data assimilation. Researchers have shown the EnKF method to be better at simulating heavy rainfall caused by mesoscale convection systems, such as those that occur during the pre-summer rainy season in South China(Bao et al.2017).However,the EnKF method lacks the constraints of numerical models and cannot correct errors over the window of assimilation; thus, it takes longer for the EnKF to spin-up in mesoscale weather applications(Kalnay and Yang 2010).

Based on the merits and demerits of the 4DVar and EnKF methods,a hybrid assimilation method combining the advantages of both has emerged: the so-called 4DEnVar method(Tian,Xie,and Dai 2008;Tian,Xie,and Sun 2011; Lorenc 2003). Tian, Xie, and Dai (2008) and Tian, Xie, and Sun (2011) proposed the PODEn4DVar method on the basis of proper orthogonal decomposition and ensemble forecasting. Zhang et al. (2015) created a PODEn4DVar-based radar data assimilation scheme to assimilate radar data for heavy-rain forecasting. Tian and Feng (2015) proposed a nonlinear leastsquares 4DVar (NLS-4DVar)-enhanced POD-4DVar method, and demonstrated that POD-4DVar is the first iteration of the NLS-4DVar version.

NLS-4DVar obtains flow-dependent background error covariance through ensemble perturbations without invoking an adjoint model, which significantly reduces computational complexity.Zhang et al.(2017a)used the NLS-4DVar method to assimilate radar data; their approach showed an improvement in the prediction accuracy of precipitation intensity and location due to the adjustments of the initial field’s wind and water vapor. Tian et al. (2018) explored the relationship between the En4DVar and 4DEnVar methods and further improved the algorithm to minimize the cost function while improving the accuracy of NLS-4DVar with minimal computational load.

The multigrid method is an efficient approach to accelerate iterative processes to solve linear and nonlinear optimization problems. For example, the multigrid method used in the Space and Time Mesoscale Analysis System has been applied to 2D radar radial-velocity data assimilation of typhoon structures(Li et al.2010).Fu et al.(2016)also used the Space and Time Mesoscale Analysis System to assimilate 3D Doppler radar radial velocity data,and reconstructed the 3D structure of a typhoon in an idealized experiment. Compared with NLS-4DVar, multigrid NLS-4DVar(referred to as MG_NLS-4DVar)minimizes the errors from different scales and has higher efficiency and precision in the assimilation of conventional data(Zhang et al.2017b;Zhang and Tian 2018a).

In this study,the assimilation capabilities of MG_NLS-4DVar were assessed in observing system simulation experiments (OSSEs) of radar radial velocity. Section 2 briefly reviews the MG_NLS-4DVar method.In section 3,we describe the OSSEs used to assess the capabilities of the main formulation of the MG_NLS-4DVar method for assimilating radar radial velocity and the advantages of the multigrid technique with the Advanced Research Weather Research and Forecasting (WRF-ARW) model.Section 4 summarizes the findings.

2. Material and methods

2.1 Brief review of multigrid variational data assimilation schemes

Ide et al.(1997)presented the cost function of multigrid variational data assimilation for any given grid level i:

where k is the observation time, and S+1 represents the total number of observation times in the assimilation window; i=1，2，...，n is the ith grid level; n represents the total number of levels;is the linearization of the observation operator; the superscripts T and -1 denote the matrix transposition and inverse, respectively;andare the background and observation error covariance matrices, respectively; δx（i）=x（i）-is the perturbation of the background field,andis the state variables;andis the innovation at the ith grid level.

A nonlinear optimization algorithm can be used to solve the cost function(Equation(1))at the ith grid level to calculate the analysis increment δx（i）of individual levels, as described elsewhere (Xie et al. 2005, 2011;Zhang and Tian 2018a).

Minimizing the 4DVar cost function at different scales requires repeated calls to the forecast model, observation operator, and their adjoint model, which is a challenge for 4DVar. Here, we introduce MG_NLS-4DVar,which avoids the need for an adjoint model.

2.2 Brief review of MG_NLS-4DVar

Generally, when NLS-4DVar is used to solve each grid level, it is necessary to repeatedly invoke the forecast model and observation operator for the simulation,resulting in higher computational costs. To minimize the computational load, MG_NLS-4DVar only simulates samples on the finest mesh. On the other hand, the ensemble samples of other grid levels are interpolated from the finest mesh. The implementation of MG_NLS-4DVar is divided into the following four steps:

Step one: According to the particular assimilation problem,determine the total number of grid levels n.

Step two:Prepare the observation yo，k，k =0，1，...，S,the background field,and samples on the finest gridwhere N represents the number of ensemble samples andis the linearization of the nonlinear forecast modelfrom t0to tk.

Step three:Calculate the analysis increment δx（i）at the ith mesh using the NLS-4DVar method(i=1，2，...，n):

where Γi→j（）（·）represents an interpolation operator from i to j,which are the numbers of grid levels.The forecast model（·）at each grid level simplifies Mt0→tk（·）.In addition,the efficient local correlation matrix decomposition scheme, proposed by Zhang and Tian (2018b), is used in the horizontal localization calculation. The vertical localization adopts the same scheme as Zhang et al.(2004).are matrix that related to the localization;model perturbationsare interpolated from corresponding data at the finest grid level. The simulated observation perturbations at the ith grid level are as follows:

When i=1,

and when 2 ≤i ≤n,

Step four: If i ＜n, we need to repeat step 3 until it reaches the finest grid level (i=n). Finally, we have xaas the sum of the background field and the analysis increment on all grid levels:

When l=1 in Equation(2),that is the first iteration,β0，（i）=0,

Equation (11) describes the POD-4DVar method (Tian,Xie,and Sun 2011;Zhang and Tian 2018a).The MG_NLS-4DVar method used in this study solves the variational minimization problem through Equation (11). The simplification for Equation (2) improves assimilation effi-ciency. As such, MG_NLS-4DVar provides the multiscale, flow-dependent background error covariance matrices B. Using the four-step approach described,the background error can be corrected sequentially from different scales to improve the final analysis. The study by Zhang and Tian (2018a) provides more details on MG_NLS-4DVar.

2.3 Observation operator for radar velocity

The observation operator for radial velocity Vradarwas proposed by Sun and Crook(1997,1998):

in which （u，v，w） are the zonal wind, meridional wind,and vertical velocity, respectively; rdisis the distance between the observation location （xobs，yobs，zobs） and the radar location （xradar，yradar，zradar）; and Vtmis the mass-weighted terminal velocity of the precipitation:

The correction factor a is defined as

where P0is the pressure at the ground,ˉP is the base state pressure, qris the rainwater mixing ratio in units of g kg-1,and ρ is the air density.

3. Observing system simulation experiments

3.1 Data and model

In the OSSEs, ARW-WRF version 3.9 was the forecast model. The number of vertical layers was 30, with the model top at 50 hPa,from η=0 to η=1.The domain was(27°-33°N,117°-125°E).The final level was n=3 for the MG_NLS-4DVar experiments in the horizontal direction.The finest grid level had 260×220(longitude×latitude)grid points with a 3-km horizontal resolution. In the horizontal direction, the number of grid points from the coarser grid to the finer grid increased two-fold(i.e.from 65 × 55 to 130× 110 and then to 260 ×220). The physics options of the forecast model were the same as employed in Zhang et al. (2015). In the OSSEs, the boundary conditions and first-guess field were generated from the National Centers for Environmental Prediction (NCEP)final(FNL)operational global analysis data. The assimilated data consisted of simulated radar radial velocity observations of Ningbo radar station,located in Jiangsu Province, China (30.07°N, 121.51°E),at an altitude of 458.4 m.

3.2 Experimental design

In 2012, Typhoon Haikui struck Jiangsu and Zhejiang provinces, making landfall at Xiangshan County,Ningbo City,at 1920 UTC 7 August 2012,and produced heavy rain from 7 August to 9 August 2012.The OSSEs in our study used data from the heavy rainfall event in the Ningbo area.The test design is described below.

The analysis time was 0000 UTC 8 August 2012,with an analysis window of 0000-0100 UTC 8 August 2012.Radar data were assimilated every 12 min. The analysis variables included the temperature, pressure, water vapor mixing ratio(qv),and the three wind components(u,v and w winds).The horizontal localization radius was 150 km.

The ‘true’ atmospheric state was initiated from 12 h prior to the analysis time at 0000 UTC 8 August 2012 with NCEP FNL analysis data. The ‘true’ initial field was used to initiate the 24-h simulation. Correspondingly,the background fields (without data assimilation) were generated from a 24-h forecast initialized by the NCEP FNL analysis data,12 h prior to the time of analysis.Then,we took the 24-h integration as the control run(referred to as CTRL) to prepare the basic states for assimilation and to assess the results. On this basis, the observation was constructed by the observation operator and the‘true’ atmospheric state in the assimilation window.Ningbo radar station had nine elevation scans, with elevation angles of 0.5°, 1.5°, 2.4°, 3.4°, 4.3°, 6.0°, 9.9°,14.6°, and 19.5°. Simulated observation perturbations and model perturbations for the MG_NLS-4DVar method were formed from two sampling runs. On the basis of a 4D moving-sampling strategy (Tian and Feng 2015),these perturbations were used to generate the ensemble samples. Specifically, we ran the model from 0000 UTC 6 August 2012 and 0000 UTC 7 August 2012 to 0000 UTC 8 August 2012 with NCEP FNL data.Then,two 14-h forecasts were intercepted from the previous simulations that included the analysis time.Each 14-h forecast included 61 1-h moving windows. Thus, the size of the ensemble was 122.

Three sets of comparative tests were designed in this study. MG_NLS-4DVar, with a final level of n=3(referred to as MG_NLS), was compared with NLS-4DVar with a maximum iteration number of Imax=3(referred to as NLS) in the OSSEs. In addition, we investigated the impact of multiscale errors with respect to the total number of grid levels n on the assimilation accuracy.Finally,we designed an experiment to explore the sensitivity of the assimilation results to the localization radius.

3.3 Experimental results

Root-mean-square errors (RMSEs) were calculated between CTRL, MG_NLS with n=3, NLS with Imax=3,and the truth at the analysis time. Figure 1 indicates the vertical profiles of the RMSEs of u (zonal) and v(meridional) winds, temperature, and the water vapor mixing ratio. The RMSEs from both MG_NLS and NLS assimilation results were far less than CTRL, especially for u and v winds. This indicated that both assimilation methods are capable of effectively assimilating radar radial velocity to improve the accuracy of the initial field.Moreover,the two methods improved accuracy to the same degree.To assess the effectiveness of MG_NLS and NLS for forecast capacity, Figure 2 compares the RMSEs of the 24-h forecast of accumulated rainfall(including analysis time) with the initial condition from CTRL, MG_NLS with n=3, and NLS with Imax=3. The RMSEs of MG_NLS and NLS were smaller than those of CTRL, which implies that both methods are effective at absorbing the observations, consequently improving the 24-h accumulated rainfall forecast. Further comparison showed that MG_NLS outperformed NLS.

To evaluate their efficiency,we compared the central processing unit (CPU) times of NLS with Imax=3 and MG_NLS with n=3. The experiments were implemented serially using a Dell PowerEdge R610 server at 24-GB memory and 2.40 GHz with 16 CPUs (16 × 4 [2-thread]Intel Xeon CPU E5620). The forecast model ran on 24 cores in parallel.It is important to note that the CPU time for model implementation depended on the case and the given CPU. The CPU times for MG_NLS with n=3 and NLS with Imax=3 is different. A forecast model was run for each iteration.The CPU time for NLS was the sum of 909.520 s and 120 s, where the former is the time of the assimilation process and the latter is twice the time of the forecast model operation. The CPU time for MG_NLS was 481.459 s for the assimilation process,which was almost half that of the NLS method.Therefore, the computational cost was reduced by the multigrid strategy.

Figure 1.Vertical profiles of the RMSEs of(a)u(zonal)wind(m s-1),(b)v(meridional)wind(m s-1),(c)temperature(°C),and(d)water vapor mixing ratio(g kg-1)from CTRL,nonlinear least-squares(NLS)assimilation with Imax=3,and multigrid NLS(MG_NLS)with n=3 results at the analysis time.

To further investigate the effect of the number of total grid levels n on multi-scale error in the assimilation,we compared the performances when using different total numbers of grid levels for MG_NLS (referred to as MG_NLSi,i=2，3,where i represents the total number of grid levels). MG_NLS2represents two grid levels (130 ×110 and 260 × 220), and MG_NLS3is the same as MG_NLS with n=3 referred to earlier. For comparison,NLS1(with one grid level; 260 × 220 grid points) is the same as NLS with Imax=3 referred to earlier.The RMSEs of the 24-h forecast of accumulated rainfall with CTRL and the three schemes, NLS1and MG_NLSi, had smaller RMSEs than CTRL (Figure 3), indicating that the three methods assimilated the observations correctly,creating a more accurate rainfall forecast. The forecast accuracy increased gradually with the number of total grid levels.MG-NLS3outperformed MG-NLS2and NLS1;thus,the three grid levels more effectively accounted for and minimized background error, starting from large-scale error and working towards error at a small scale, for a more accurate analysis.

The total times for the NLS1, MG_NLS, and MG_NLS3methods were 300.368, 483.958, and 601.459 s, respectively (Table 1). The computational cost increased with the number of horizontal grids.The time for the forecast model operation to run once was about 60 s.Although MG-NLS3had the highest computational cost,its higher forecast accuracy was worth the extra but almost negligible computational cost. Moreover, compared with NLS with Imax=3 referred to earlier,which is slightly lower accuracy, MG-NLS3had greatly improved computing efficiency.

Figure 2.RMSEs of the 24-h forecast of accumulated rainfall(from 0000 UTC 8 August to 0000 UTC 9 August)with the initial condition from CTRL,nonlinear least squares(NLS)assimilation with Imax=3,and multigrid NLS(MG_NLS)with n=3 methods.

Figure 3.RMSEs of 24-h forecast of accumulated rainfall(from 0000 UTC 8 August to 0000 UTC 9 August)with the initial condition from the CTRL, nonlinear least squares (NLS1), and multigrid NLS (MG_NLSi) (i=2，3 represents the total number of grid levels)methods.

Table 1.CPU time of each grid level and total time for NLS1 and MG_NLSi(i=2，3,where i represents the number of grid levels)(units:s).

To test the sensitivity of MG_NLS with n=3, the RMSEs of the 24-h forecast of accumulated rainfall of MG_NLS with five localization radii were compared(Figure 4).The RMSE was smallest when the localization radius was 150 km. Therefore, this radius was used for MG_NLS and NLS in the OSSEs. The MG_NLS performance was similar when the localization radius varied from 120 to 210 km, indicating that MG_NLS has low sensitivity to the localization radius of the interval.

4. Conclusions

Ensemble-based variational data assimilation methods have flow-dependent background error covariance matrices B and avoid the adjoint model; however, their accuracy depends on iterative corrections on a given scale.The assimilation of unconventionally observed radar data using a large number of observations and high time frequency incurs a high computational cost.The MG_NLS-4DVar method used in this paper not only has the advantages of ensemble-based variational data assimilation,but also corrects the background error from large to small scales.Ultimately,the analysis field is the sum of the background field and the analysis increment on all grid levels.In this paper,a triple-grid MG_NLS-4DVar method was used to assimilate the radar radial velocity of Ningbo station to assess the capabilities of MG_NLS-4DVar.The main conclusions can be summarized as follows.

Figure 4.RMSEs of the 24-h forecast of accumulated rainfall(from 0000 UTC 8 August to 0000 UTC 9 August)with the initial condition from CTRL and multigrid nonlinear least-squares assimilation(MG_NLS)(n=3)with different localization radii.

Both NLS-4DVar and MG_NLS-4DVar can effectively assimilate radar radial velocity to improve the initial field and enhance the 24-h cumulative precipitation forecast accuracy. Also, MG_NLS-4DVar outperforms NLS-4DVar.In terms of assimilation CPU time, the computational cost of NLS-4DVar with three iterations on the finest grid was found to be higher. The precision of MG_NLS-4DVar improved as the total number of grids increased.Therefore,the MG_NLS-4DVar method,with dual advantages in assimilation accuracy and computational effi-ciency, is a promising tool for radar radial velocity data assimilation. Only a limited area was covered in this paper because the radar radial velocity observations were from a single radar. Multiple radar observations and the radar reflectivity are expected to be assimilated with MG_NLS-4DVar in the future.

Acknowledgments

We thank Dr. FENG Chen at the Zhejiang Institute of Meteorological Sciences for providing the radar data and helpful comments and suggestions, and Dr. ZHANG Bin for data analysis. We thank NCAR, Boulder, Colorado, USA, for making the WRF-ARW model available (http://www2.mmm.ucar.edu/wrf/users/downloads.html). The initial and lateral boundary conditions are generated by NCEP FNL operational global analysis data(https://rda.ucar.edu/datasets/ds083.2/).

Disclosure statement

No potential conflict of interest was reported by the authors.

Funding

The work was partially supported by the National Key Research and Development Program of China [grant number 2016YFA0600203]; the National Natural Science Foundation of China [grant number 41575100]; and the Key Research Program of Frontier Sciences, Chinese Academy of Sciences[grant number QYZDY-SSW-DQC012].