2.2.1. Data Assimilation Scheme

The WRF data assimilation package provides various assimilation techniques, ranging from empirical (e.g., nudging) and statistical (e.g., variational analysis) approaches, to advanced methods (e.g., ensemble Kalman filter). In the variational scheme, the best possible estimate of the model's initial state (analysis) is determined by the minimization of a prescribed cost function, *J*(*x*), that combines a background forecast (first-guess), observations and estimates of both modeling and observational error covariances. 3D and 4D-var options are available in WRF [48,49]. The main difference between the two implementations is that the observations in 4D-var/WRF model are integrated within an assimilation window at the exact time of the observations. Consequently, greater computing resources compared to 3D-var assimilation are necessary [50]. Taking the possible application of GNSS ZTD data assimilation in NOA's operational forecasting system into account, without significantly altering the timeliness of forecast delivery, the 3D-var option was employed in the present work.

According to Barker et al. [48], the cost function in WRF 3D-var operation can be written as:

$$J(\mathbf{x}) = \frac{1}{2} (\mathbf{x} - \mathbf{x}\_b)^T B^{-1} (\mathbf{x} - \mathbf{x}\_b) + \frac{1}{2} (H(\mathbf{x}) - \mathbf{y}^\sigma)^T R^{-1} (H(\mathbf{x}) - \mathbf{y}^\sigma)^T$$

expressing the weighted distance of analysis, *x*, to first-guess, *xb*, and observations, *yo*. The contributions of *xb* and *yo* to *x* are determined by the model background errors covariance matrix, B, and the observations errors covariance matrix, R, respectively. The observations operator *H* is used to transform the model's analysis to observational space. The WRF first guess can be either the initial analysis computed based on global fields (cold start) or a previous model forecast (cycling). In the last-mentioned case, the observations are used to improve the current model state, ultimately resulting in better analysis and prediction. The 3D-var/WRF system supports the assimilation of in-situ conventional measurements, remotely sensed observations, and satellite radiances. Ground-based GNSS ZTD data are included in the remote sensing observations category and they are handled by the GPSPW operator [28,39].

The success of the 3D-var assimilation depends heavily on the accuracy with which the observations and model background errors are specified. The R matrix is determined by instrument and representation errors. The latter includes errors introduced by the observations pre-processing and operator, as well as by the effect of unresolved scales in the model [51]. The observations preprocessor (OBSPROC) in the 3D-var/WRF system defines the R matrix, based on pre-specified observations errors. These error values had to be determined by the user for the case of GNSS ZTD data. No correlation in space and time between the individual observation errors is assumed by the OBSPROC package [39]. In the current study, the GNSS ZTD formal errors, derived during the GNSS raw observations processing, were used in the 3D-var/WRF application. These errors are related to uncertainties induced by satellite orbits, antennas, signals multipath, ionospheric delays, and the applied mapping function. Hence, they are a measure of the ZTD estimation uncertainty that considers both measurement and processing errors. The B matrix is considered static and its specification is based on "climatological" estimates, assuming that the model background errors correlations are homogeneous and isotropic. Its computation involves the application of the National Meteorological Center (NMC) method [52], in which the covariances for a set of five independent control variables are derived for a given domain by averaging the forecast differences between the 24 and 12 h prediction lead time over a period of time [39,48]. In the present work, the B matrix was calculated for each rainfall event, considering that the model background fields are different for each case study. In particular, the B matrix was calculated based on a series of WRF simulations of 24 h duration for the domain of interest (d02), which were initialized twice daily at 0000 UTC and 1200 UTC and covered a 10 day period, starting thus 10 days before each episode occurrence. Similar approaches for computing the B matrix have also been applied in previous studies [28,53].
