The impact of global navigation satellite system (GNSS) zenith total delay data assimilation on the short-term precipitable water vapor and precipitation forecast over Italy using the Weather Research and Forecasting (WRF) model

The impact of assimilating GNSS-ZTD (global navigation satellite system–zenith total delay) on the precipitable water vapor and precipitation forecast over Italy is studied for the month of October 2019, which was characterized by several moderate to intense precipitation events, especially over northwestern Italy. The WRF (Weather Research and Forecasting) model, version 4.1.3, is used with its 3D-Var data assimilation system to assimilate ZTD observations from 388 GNSS receivers distributed over the country. The dataset was built collecting data from all the major national and regional GNSS permanent networks, achieving dense coverage over the whole area. The water vapor forecast is verified for the forecast hours of 1–6 h after the last data assimilation time. Results show that WRF underestimates the atmospheric water vapor content for the period, and GNSS-ZTD data assimilation improves this underestimation. The precipitation forecast is verified in the phases of 0–3 and 3–6 h after the last data assimilation time using more than 3000 rain gauges spread over Italy. The application of GNSS-ZTD data assimilation to a case study improved the precipitation forecast by increasing the rainfall maximum and by better focusing the precipitation pattern over northeastern Italy, with the main drawback being the prediction of false alarms. Considering the study over the whole period, GNSS-ZTD data assimilation had a positive impact on rainfall forecast, with an improvement in the performance up to 6 h and with statistically significant results for moderate to intense rainfall thresholds (25–30 mm (3 h)⁻¹).

1 Introduction

The Mediterranean area is often struck by severe weather and deep convective events because of the presence of the warm sea, the complex orography of the area and the specific synoptic-scale environment. This scenario is worsened by climate change, which is affecting many weather and climate extremes, and the frequency and intensity of heavy-precipitation events have increased in most of the world (IPCC, 2023). Numerical weather prediction (NWP) models are useful tools to predict adverse weather conditions and to guide responsive actions for mitigating the impact of severe weather. Over the past years, the use of NWP models, along with an increasing availability of computing power, has led to an improvement in the forecast accuracy. However, NWPs have well-known difficulties in representing the physical processes at small spatial and temporal scales, which are involved in convective or severe weather events (Stensrud et al., 2009).

One of the common problems of NWP-based nowcasting and short-term forecasting is the spin-up time because the model needs a few hours to balance the inconsistencies between the initial and boundary conditions to properly reproduce the small-scale dynamic (Lagasio et al., 2019). Data assimilation of local observations in NWP has been reported as a key factor in reducing this issue and improving the prediction of high-impact weather events (Federico et al., 2021). Among local observations, water vapor plays a key role because of its importance in humid and energetic exchanges in the atmosphere. Therefore, good knowledge of water vapor distribution in space and time is a fundamental requirement for improving NWP forecasts of convective and severe weather events.

Global navigation satellite system (GNSS – a collective term used to address all global and regional satellite navigation systems, including the Global Positioning System (GPS), Galileo, GLONASS and BeiDou) routine observation processing for geodetic and geophysical purposes can provide estimates of the tropospheric delays (generally ZTD – zenith total delay), directly connected to the water vapor content in the atmosphere, which can be very useful to improve NWP forecasts. Relevant, albeit non-exhaustive, experiences are shortly summarized hereafter:

• Vedel and Huang (2004) assimilated GNSS-ZTD into the High Resolution Limited Area Model using three-dimensional variational data assimilation (3D-Var) and found improvements for the forecast of geopotential height and high precipitation levels.

• Poli et al. (2007), assimilated GNSS-ZTD using four-dimensional variational data assimilation (4D-Var) and the ARPEGE (Action de Researche Petite Echelle Grande Echelle) global model. Results show the positive impact of the GNSS-ZTD data assimilation on the forecast of synoptic-scale circulations and precipitation in spring and summer. Other studies followed in France (Boniface et al., 2009; Yan et al., 2009) and found a positive impact of GNSS-ZTD data assimilation on the NWP forecast.

• Bennitt and Jupp (2012) assimilated GNSS-ZTD observations by both 3D-Var and 4D-Var using the Met Office North Atlantic and European model. The assimilation of GNSS-ZTD resulted in an improvement of the cloud forecast.

• Lindskog et al. (2017) performed GNSS-ZTD data assimilation, assimilating data into the HARMONIE-AROME model at 2.5 km horizontal resolution. The assimilation was performed by 3D-Var and improved the forecast by up to 1.5 d.

• Rohm et al. (2019) assimilated GNSS (both ZTD and precipitable water) in the Weather Research and Forecasting (WRF) model at 4 km horizontal resolution over Poland for 2 months. They found an improvement in predicting both water vapor and precipitation. In the same direction, Trzcina et al. (2020) showed the positive impact of assimilating the GNSS tomography wet refractivity field on temperature and precipitation forecasting at the short range (0–6 h) over Poland.

• Giannaros et al. (2020) showed the positive impact of GNSS-ZTD data assimilation on both precipitation and water vapor forecast over Greece, and Caldas-Alvarez and Khodayar (2020) showed similar results with the climatic setting of the COSMO model over a large portion of the Mediterranean and central Europe.

• Risanto et al. (2021) assessed the impact of GPS precipitable water vapor (PWV) data assimilation on short-range North American monsoon precipitation forecasts over northwest Mexico. They showed that GPS-PWV data assimilation created more favorable conditions for nocturnal convection organization, leading to a better precipitation forecast.

• Singh et al. (2019) studied the impact of GNSS-ZTD data assimilation over the Indian region with the WRF model. They showed a positive impact on the precipitation forecast and on the lower- to middle-tropospheric moisture, upper-air temperature, and middle- and upper-tropospheric wind.

• Yang et al. (2020) assimilated both radar reflectivity and GNSS-ZTD for a heavy-rainfall event over Taiwan. They showed the complementary role of both observations in improving the precipitation forecast.

The assimilation of GNSS data is also used in an operational context. The Rapid Refresh model over the USA assimilates GPS-derived integrated water vapor (IWV) every hour from 300 stations across the USA (Benjamin et al., 2016). The study shows that there is a clear benefit in using GNSS observations in rapid refresh weather forecasting.

Considering the GNSS data assimilation over Italy, Faccani et al. (2005), using the Mesoscale Model 5 (MM5) at 9 km horizontal resolution and 3D-Var to assimilate GNSS-ZTD over Italy, found improvements in the precipitation forecast during the transition from winter to spring.

By assimilating a wide range of Sentinel-1 and GNSS-ZTD observations over Italy into the WRF model, Lagasio et al. (2019) found that the forecasts benefit the most when the model is provided with information on the wind field and/or the water vapor content. Mascitelli et al. (2019, 2021) reported two successful experiments of GNSS-ZTD and PWV data assimilation with the RAMS@ISAC model in Italy. The 3D-Var system was used to assimilate GNSS-ZTD, while nudging was used to assimilate the PWV. In both cases the assimilation showed a significant improvement in the short-term prediction of water vapor with smaller impact on the precipitation forecast.

This paper enriches the numerical experiments made over Italy to improve the precipitation and water vapor forecast through GNSS-ZTD data assimilation because it refers to a different period compared to previous studies (October 2019) and uses the data of 388 GNSS receivers widely spread over the country for the whole period, giving a robust assessment of the impact that GNSS-ZTD data assimilation can have on forecasting at the local scale. In addition, it considers two issues that are important in the operational context: the optimal spacing of GNSS receivers for data assimilation and bias removal. The first experiment estimates the observation-error decorrelation length scale and applies it to optimize the GNSS-ZTD data assimilation by data thinning, while the second experiment quantifies the impact of the bias removal on the forecast performance because the bias cannot be completely removed in an operational context.

It is also noted that a similar configuration of the WRF model used in this paper is already operational at the National Research Council of Italy, Institute of Atmospheric Sciences and Climate (CNR-ISAC), in the framework of the agreement between the Civil Protection Department (DPC) and CNR-ISAC to improve the NWP forecast in different time ranges. So the results of this paper are of practical importance as GNSS-ZTD data could also be assimilated in the near future in the operational run.

This paper is organized as follows: Sect. 2 shows the data and methods used, including the model setting, the GNSS-ZTD dataset and the statistics for the verification. Section 3 discusses the results, showing the impact of GNSS-ZTD data assimilation on the WRF analyses (both ZTD and PWV), on the precipitation forecast for a case study, and on the PWV and precipitation forecast for the whole period. Conclusions are given is Sect. 4. Finally, more details on the statistics used in the paper as well as the resampling test are reported in the Appendix.

2 Data and methods

2.1 WRF model configuration and assimilation method

The model used in this study is the Weather Research and Forecasting model with advanced WRF dynamics (WRF-ARW, Version 4.1.3; Skamarock et al., 2019). The model was run over one domain covering the whole Italian territory and central Europe (Fig. 1) with a spatial horizontal resolution of 3 km. The model grid has 635 × 635 grid points and 50 vertical levels, with the model top at 50 hPa. The main physical parameterizations used are the following: the Thompson scheme (Thompson et al., 2008) for microphysics; the Mellor–Yamada–Janjić (Eta) turbulent kinetic energy (TKE) scheme for the boundary layer (Janjic, 1994); and the Dudhia scheme (Dudhia, 1989) and the Rapid Radiative Transfer Model (RRTM; Mlawer et al., 1997) for shortwave and longwave radiation, respectively.

Figure 1 WRF model domain and GNSS receivers' height [m] above sea level.

The experiment aims at evaluating the impact of GNSS-ZTD data assimilation on the precipitation prediction in the short term (up to 6 h in this paper). For this purpose, we considered a 1-month period, from 2 to 31 October 2019. The choice of the period is due to two main reasons: on the one hand, this month was characterized by both widespread and localized precipitation events over Italy; on the other hand, a dense network of GNSS receivers (about 500) is available for the country. It is important to note that the number of receivers actually used for data assimilation was reduced to 388 by applying data thinning as discussed in Sect. 2.2.

GNSS-ZTD data assimilation was performed using the 3D-Var tool distributed with the WRF model, which is one of the components of the WRF Data Assimilation (WRFDA) system and also includes 4D-Var and ensemble data assimilation systems (Barker et al., 2004, 2012; Huang et al., 2009).

We consider two kinds of simulations: control simulations, without GNSS-ZTD data assimilation, hereafter also CTRL, and simulations assimilating GNSS-ZTD, hereafter also GNSS. The European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS) 3-hourly operational analysis–forecast cycle at 0.25^∘ starting at 12:00 UTC on the day before the actual day to forecast is used for initial and boundary conditions to simulate a real forecasting context. The temporal scheme used for the simulations uses a very short-term forecast (VSF) approach, with a 6 h update (Fig. 2).

Figure 2 Rapid update cycle at 6 h. Red dots denote analysis times.

In this scheme, for each day, we run four simulations starting from a cold start. Each simulation lasts 12 h. The first 6 h of each run is used for the model spin-up and for data assimilation in GNSS simulations, while the last 6 h is used as forecast. Therefore, four runs are necessary to cover a whole day. For GNSS simulations, in the assimilation phase, we considered an analysis every 1 h (red points of Fig. 2), starting from the beginning of the simulation and reaching the sixth hour, so a total of seven analyses are performed for each run. For CTRL simulations we use only initial and boundary conditions from ECMWF IFS and no other data are assimilated.

GNSS-ZTD observations are assimilated by 3D-Var, which is a variational approach involving the calculation of the analysis that minimizes a cost function by measuring its distances from the background and observations. The cost function is given by

$$\begin{array}{rl}J\left(x\right)=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{x}-{\mathit{x}}_{\mathrm{b}}\right)}^T{\mathbf{B}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_{\mathrm{b}}\right)\\ & +{\left[{\mathit{y}}_{\mathrm{o}}-H\left(\mathit{x}\right)\right]}^T{\mathbf{R}}^{-\mathrm{1}}\left[{\mathit{y}}_{\mathrm{o}}-H\left(\mathit{x}\right)\right],\end{array}$$

where x is the state vector; x_b is the background field; H is the forward observational operator which transforms the x vector in the observational space; y_o is the observation vector; and B and R are the error covariance matrices, respectively, for the background and observations.

The R matrix is given by the sum of instrumental and representation errors. To prepare observational data for 3D-Var, the obsproc tool was employed for GNSS-ZTD observations. Data for each analysis were considered on a 1 h time range, centered at the hour of the analysis. The R matrix is diagonal, which requires data thinning (see next section), and the ZTD error is set to 5 mm. Different regional networks are considered to reach the considerable number of GNSS-ZTD receivers used in this paper; however, we used a constant value for the errors in all receivers because the software and the processing method are the same for all the receivers. Also, the GNSS-ZTD time series were visually checked, and no specific differences among networks arose. The value of 5 mm was not specifically computed for this experiment but comes from previous comparisons that, in any case, do not extend to the whole of Italy (Tagliaferro, 2021; Krietemeyer et al., 2018; Mascitelli et al., 2019, 2021). In these works, the GNSS-ZTD retrieved with the method used in our paper was compared with other methods and with radiosoundings. In general, comparison with radiosondes shows differences in the range of 1.0–1.5 cm (i.e., larger than the error used in this paper), while comparison with other methods shows differences between 0.1 and 0.8 mm. Now, the comparison with radiosondes is less representative of the GNSS-ZTD error because radiosondes can move far from the GNSS receiver, and the 0.5 cm used in this paper comes from the comparison with other methods to estimate ZTD. However, future experiments considering different errors for different networks should be done to assess this point in more detail.

As regards the B matrix, the calculation for this study is performed via the GEN_BE tool, employing the National Meteorological Center (NMC) method (Parrish and Derber, 1992) for the month of October and the option CV5, which accounts for five control variables.

The GNSS-ZTD observations are considered to have unbiased errors compared to the WRF model. To achieve this goal, a statistical bias correction was applied to the ZTD data for the whole period. First the raw GNSS data are assimilated into 3D-Var to calculate the corrections that come from the background. The difference between the observation and the background is saved for each receiver and for each time giving the quantity $(O-B{)}_{k,t}$, where k is the receiver index and t is the time. The quantity $(O-B{)}_{k,t}$ takes into account the difference between the model orography and receiver height that, in our case, is never larger than 300 m. For each receiver we compute the background bias by averaging $(O-B{)}_{k,t}$ over all times:

$$\stackrel{\mathrm{‾}}{{\left(O-B\right)}_k}=\sum _{t=\mathrm{1}}^N{\displaystyle \frac{{\left(O-B\right)}_{k,t}}{N}},$$

where N is the total number of times (i.e., observations) available for each GNSS receiver. Then we use the corrected observation $O_{k,t}^{\prime }$ in 3D-Var:

$$O_{k,t}^{\prime }=O_{k,t}-\stackrel{\mathrm{‾}}{(O-B{)}_k}.$$

This method is like that applied in many papers including some cited in the Introduction.

2.2 Observational datasets and verification procedure

2.2.1 GNSS data

To achieve a good spatial density of the ZTD measurements, the GNSS data have been collected from all major Italian national and regional GNSS networks. To derive ZTDs, GPS L1 and L2 pseudorange and phase measurements were adjusted in precise point positioning (PPP) mode (Bevis et al., 1992; Zumberge et al., 1997) by using the GNSS data processing suite Breva, which is a proprietary, closed-source version of the free and open-source goGPS software (Herrera et al., 2016). Both Breva and goGPS are developed by GReD.

Several of the used stations tracked only GPS satellites; thus we decided to process only GPS observations, disregarding data from other GNSS constellations. Table 1 describes the corrections applied to the observations, as well as the stochastic model.

Table 1 Corrections applied to the observations and stochastic model.

Figure 1 presents the locations and the heights of the 388 GNSS stations used over the Italian Peninsula and its surroundings after a reduction in their number. Although the initial number of GNSS receivers was about 500, more than 100 GNSS receivers were discarded considering the following two requirements: (a) the difference between the receiver and WRF model height at the closest grid point to the GNSS receiver must be less than 300 m (similarly to Bennitt and Jupp, 2012; Rohm et al., 2019; and Mascitelli et al., 2019) and (b) in the case that two or more receivers fell in the same WRF grid cell, the one whose height was closer to the model orographic height was retained. At the end of the process, 388 GNSS receivers were used for data assimilation and the closest distance between any two receivers was larger than 10 km. No further data thinning was applied.

2.2.2 Forecast verification

As regards the evaluation of the precipitation forecast performances, five precipitation scores are calculated, i.e., the frequency bias (FBIAS), probability of detection (POD), threat score (TS), equitable threat score (ETS) and false alarm rate (FAR). A detailed description of these scores is provided in Appendix A.

In the following, we will show the results of the scores, summarizing them through a performance diagram (Roebber, 2009), with the exception of the ETS. In this diagram the perfect score is the one that reaches the upper-right corner of the diagram. The x axis represents the success ratio (SR), which is defined as 1 − FAR, while the POD is on the y axis. The straight lines from the origin represent the FBIAS, while the hyperboles branches are the TS.

We considered the forecast for both types of run, CTRL and GNSS, taking into account two main periods after the last analysis time in order to evaluate the impact of GNSS-ZTD data assimilation on the precipitation forecast at different ranges: the first 3 h (from the 6th to the 9th hour of run) and the second 3 h of forecast (from the 9th to the 12th hour of run). The model forecast in correspondence to a given rain gauge is computed using the nearest-neighbor method. By this method, we consider all WRF precipitation values at grid points in a radius of $\mathrm{2}\mathrm{\Delta }x\sqrt{\mathrm{2}}$ from the rain gauge, with Δx being the WRF model grid spacing, and we select the grid point with the rainfall value closest to the rain gauge observation.

Precipitation data come from the Italian rain gauge network, with more than 4000 rain gauges over Italy. This network belongs to the Italian regional administrations, and data are collected nationwide by the Civil Protection Department (Davolio et al., 2015).

Furthermore, a test to assess the statistical significance of the differences between CTRL and GNSS precipitation forecasts is performed. The resampling test of Hamill (1999) is used (see Appendix B for details).

For precipitable water vapor verification, we focused on both the assimilation and forecast phases. Verification is made calculating two main scores.

The first is BIAS, which measures the mean differences between forecast (F) and observation (O):

$$\begin{array}{}\text{(1)}& \mathrm{BIAS}={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N\left(F_i-O_i\right).\end{array}$$

The second is root mean square error (RMSE), which measures the mean of the squared differences between forecast and observation:

$$\begin{array}{}\text{(2)}& \mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N{\left(F_i-O_i\right)}^{\mathrm{2}}},\end{array}$$

with N number of forecast–observation pairs considered in the statistics.

3 Results and discussion

3.1 ZTD and precipitable water vapor results for the analysis phase

In this section we consider the differences between the first guess (FG) and the observations and between the analyses (ANL) and the observations. Statistics are presented considering all simulations (30 × 4). For each simulation, six times are considered for the first guess and analysis, one for each time (red dots of Fig. 2), except for the analysis at the initial time of each simulation because it is coincident with the last analysis time of the previous simulation, the latter being considered in the results.

Two types of statistics are shown. In the first statistics we present the time series of the bias and RMSEs aggregating all receivers together for each analysis time; in this case, referring to Eqs. (1) and (2), N is the number of GNSS receivers considered at the analysis time. In the second statistics, we show the PWV BIAS and RMSE computed for each receiver; in this case, referring to Eqs. (1) and (2), N is the number of analyses performed for the whole period (30 × 4 × 6 = 720 first guess and analysis pairs).

Figure 3 shows the time series for ZTD BIAS (Fig. 3a) and ZTD RMSE (Fig. 3b) calculated for the whole period. From Fig. 3a it is apparent that the ANL BIAS (blue curve) has lower absolute values compared to that of FG (red curve) for all times considered. The better results obtained by ANL is confirmed also by the BIAS values calculated over all times together, shown in the upper-right part of the figure, which is halved after the analysis (BIAS_A, −0.1 cm) compared to the first guess (BIAS_F, −0.2 cm). Interestingly, the WRF first-guess bias is mainly negative, showing an underestimation of the water vapor content in the atmosphere. This underestimation is reduced by GNSS-ZTD data assimilation that increases the water vapor in the model. The positive impact of the GNSS-ZTD data assimilations is confirmed by the RMSE calculation (Fig. 3b). The RMSE for the whole period for the analyses (RMSE_A) is 0.7 cm, compared to the value of 1.3 cm for the first guess (RMSE_F).

Figure 3 Time series of ZTD BIAS (a) and of ZTD RMSE (b) calculated for all sensors together for first guess (red curve) and analysis (blue curve) during the analysis phase. Dates are shown along the x axis.

It is important to discuss the impact of GNSS-ZTD data assimilation over the Italian territory. This is shown considering the PWV estimated by the ZTD. The PWV in millimeters is given by

$$\mathrm{PWV}=Q({\mathrm{ZTD}}_{\mathrm{GNSS}}-{\mathrm{ZHD}}_{\mathrm{WRF}}),$$

where ZHD_WRF is the hydrostatic delay calculated using the Saastamoinen (1972) equation given the WRF surface pressure (p_sfc), latitude (φ) and height (h):

$${\mathrm{ZHD}}_{\mathrm{WRF}}={\displaystyle \frac{\mathrm{0.0022767}{p}_{\mathrm{sfc}}}{\mathrm{1.0}-\mathrm{0.0266}\cos \left(\mathrm{2}\mathit{\varphi }\right)-\mathrm{0.00000029}h}}.$$

Similarly to other studies (for example, Rohm et al., 2019), we estimate the ZHD from the WRF surface pressure because no pressure observations were available in correspondence to the GNSS-ZTD receivers.

The proportionality factor Q is computed as follows:

$$Q={\displaystyle \frac{{\mathrm{10}}^{\mathrm{6}}}{R_{\mathrm{w}}(\frac{k_{\mathrm{3}}}{T_{\mathrm{m}}}+k_{\mathrm{2}}^{\prime })}},$$

where R_w= 461 J (kg K)⁻¹ is the gas constant for water vapor, $k_{\mathrm{2}}^{\prime }=$ 22.9726 and k₃= 375 463 K hPa⁻¹ are the refractivity constants from Rueger (2002), and T_m is the mean temperature given by

$$T_{\mathrm{m}}=\mathrm{70.2}+\mathrm{0.72}{T}_{\mathrm{WRF}},$$

where T_WRF is the WRF surface temperature.

Figure 4 shows the PWV RMSE (panels a and b for first guess and analysis, respectively) and the BIAS (panels c and d for first guess and analysis, respectively) for the whole period at each GNSS receiver. Both statistics are improved by data assimilation. The impact of GNSS-ZTD data assimilation on PWV RMSE is apparent as most of the green–orange–red colors (RMSE between 1.0 and 2.5 mm, Fig. 4a) are reduced to dark-green–green colors (RMSE between 0.5 and 1.5 cm, Fig. 4b). For the BIAS the improvement is less apparent, but it is clearly shown by the increase in the orange, red and yellow dots in panel (d) compared to panel (c) of Fig. 4.

Figure 4 PWV RMSE for FG (a) and ANL (b) and PWV BIAS for FG (c) and ANL (d) calculated for each sensor over all times of analysis.

All in all, the results of this section show a positive impact of GNSS-ZTD data assimilation on the ZTD and PWV over the period, and the net result is the increase in the model water vapor content over the area, because the first-guess underestimation of the ZTD (and PWV) is reduced by data assimilation.

3.2 Rainfall forecast for a case study

The selected date is the 15 October 2019, in the 3 h between 18 and 21:00 UTC. This phase was selected because it is representative of the improvement that the assimilation of GNSS-ZTD had for several cases on the 3 h precipitation forecast. During this time period, three main thunderstorms developed: the first two over the northeast and northwest of Italy, respectively, with a precipitation maximum between 60 and 70 mm (3 h)⁻¹, and the third in southern Italy, with precipitation ranging between 40 and 50 mm (3 h)⁻¹. We start examining the innovations, i.e., the analysis minus first-guess fields, at 18:00 UTC on 15 October (Fig. 5a), which is the last analysis before the 3 h forecast considered in this section and has an important role in the 3 h rainfall forecast. Indeed, as shown below, the precipitation between 18:00 and 21:00 UTC has several correspondences with the innovations at 18:00 UTC.

Figure 5 Innovations of water vapor at 18:00 UTC on 15 October 2019; (a) horizontal map around 1800 m above the surface. Cross section of (b) is taken along the red line; (b) latitude–height cross section of water vapor innovations with the first-guess temperature (black contours).

Figure 5a shows the innovations at 18:00 UTC on 15 October 2019 at about 1800 m above the ground surface. A complex pattern of positive and negative innovations is shown thanks to the many observations available, which add information on the water vapor field at the subregional scale (the distance between two closest maxima or minima of Fig. 5a can be roughly estimated in 50–70 km).

Figure 5b shows the latitude–height cross section of the innovations at the same time of Fig. 5a. The cross section is for the longitude 12.6^∘ E. Water vapor values in the range 0.4–1.0 g kg⁻¹ are shown in several parts of the cross section at about 1600 m height a.g.l. in western Sicily (37.5^∘ N) revealing a considerable impact of GNSS-ZTD data assimilation (> 10 % of the first-guess value) on the water vapor field.

Between 18:00 and 21:00 UTC, moderate to intense precipitation are shown in NE and NW of Italy by rain gauge observations, with maximum intensities larger than 60 mm (3 h)⁻¹ (Fig. 6a). Some areas of moderate precipitation are apparent in central and southern Italy with maxima between 20 and 50 mm (3 h)⁻¹.

Figure 6 Precipitation between 18:00 and 21:00 UTC on 15 October 2019: (a) rain gauge observations; (b) CTRL forecast; (c) GNSS forecast.

The CTRL, Fig. 6b, has a good forecast because it represents the two main precipitation areas over NW and NE of Italy. There are, however, less satisfactory points in the CTRL rainfall forecast: first, the maximum over NE of Italy is composed by two branches and one of them is displaced close to the sea, starting from Venice and going toward NE. This branch is a false alarm, at least in its southernmost part. Second, the maximum over central Italy, starting from Tuscany and displaced in the SW–NE direction is displaced to the north, compared to observations. Third, the precipitation over southern Italy is missed.

The assimilation of GNSS-ZTD improves the precipitation forecast. The pattern of the precipitation over NE is more in agreement with observations because precipitation does not appear separated in two branches as in the CTRL and the precipitation amount for GNSS simulations is higher, better catching the observation maximum. This agrees with the increase in the water vapor over NE of Italy given by the last assimilation (Fig. 5a). A similar improvement is apparent in the south of Italy, where the precipitation was missed by CTRL. In this case the observed precipitation is well forecast by GNSS and the increase in precipitation is determined by the increase in water vapor over the area (Fig. 5a). The maximum over central Italy is not much improved by the assimilation of GNSS-ZTD because, even if there is better superposition with observations compared to the control forecast, there are more false alarms in the GNSS forecast.

Section S3 of the Supplement shows the rainfall scores for CTRL and GNSS for this case study. GNSS forecast has a better performance compared to CTRL, especially for thresholds larger than 40 mm (3 h)⁻¹.

3.3 Results for the water vapor forecast for the whole period

In this section we show the impact of GNSS-ZTD data assimilation on the PWV forecast. Three different statistics are considered: the time series of the differences between the PWV RMSE of CTRL and GNSS simulations for different forecasting times; the time series of the difference between the PWV absolute value of the BIAS of CTRL and GNSS simulations for different forecasting times; the RMSE maps of PWV for CTRL and GNSS on two different forecasting times. The PWV is computed as shown in Sect. 3.1, replacing ZTD_GNSS with ZTD_WRF

The first two statistics are shown in Fig. 7 for the first (red curve), third (blue curve) and sixth hour of forecast (green curve). Considering the differences between CTRL and GNSS PWV RMSE, positive values show better performance of the forecasts assimilating GNSS-ZTD, and, from Fig. 7a, it is apparent that the GNSS-ZTD data assimilation improves the PWV forecast at different forecast ranges, with the exceptions of a few cases. The RMSE of the GNSS simulations can have an RMSE lower than CTRL simulations up to 2.5 mm, and, from Fig. 7a, a decrease in the improvement of GNSS-ZTD data assimilation for longer forecasting times is also notable as the red curve shows, with some exceptions, larger values than the blue curve, which, in turn, has larger values that the green curve. However, it is important to note that the positive impact of GNSS-ZTD data assimilation on the PWV forecast is still apparent after 6 h of forecast.

Figure 7 Time series of the differences between PWV-RMSE for CTRL and GNSS (a) and of the difference in the absolute value of the BIAS for the CTRL and GNSS simulations (b). The red curve is for the first forecast hour, the blue curve is for the third forecast hour and the green curve is for the sixth forecast hour. The date format is the month followed by the day.

The bias is also improved by the GNSS-ZTD data assimilation as the difference in the absolute value of the bias for CTRL and GNSS simulations shows positive values. As for RMSE, the bias improvement decreases with forecasting time.

It is interesting to evaluate the spatial distribution of the improvement of GNSS-ZTD data assimilation on the PWV forecast over Italy. This is shown in Fig. 8 for the first hour of forecast (panels a and b for CTRL and GNSS, respectively) and for the sixth hour of forecast (panels c and d for CTRL and GNSS, respectively). Figure 8 shows that the improvement of RMSE is not limited to a specific area, thanks to the good coverage of the GNSS receivers used in this work, but it is widespread with RMSE more than halved in correspondence to most receivers.

Figure 8 Maps of the PWV RMSE for the first (a, b) and sixth (c, d) forecasting hour. Panels (a) and (c) are for CTRL simulations; panels (b) and (d) are for GNSS simulations.

The RMSE for the sixth forecast hour increases for both CTRL and GNSS simulations compared to the first forecast hour, as expected because the forecast error increases with forecasting time. However, the increase in the RMSE is larger for the GNSS simulations (comparison between panels b and d) than for the CTRL simulation (panels a and c), and the impact of GNSS-ZTD data assimilation decreases with the forecasting time over the whole Italian territory.

3.4 Results for the precipitation forecast over the whole period

In this section we show the impact of GNSS-ZTD data assimilation on the precipitation forecast statistics for the whole period (2–31 October 2019). Two precipitation phases are considered: the first 3 h of each simulation after the last assimilation time (also 0–3 h) and the following 3 h, i.e., from the third to the sixth hour of each simulation after the last assimilation time (hereafter 3–6 h), to evaluate for how long the GNSS-ZTD data assimilation impacts the precipitation forecast.

Performance diagrams are shown in Fig. 9 for both CTRL and GNSS and for three different precipitation thresholds, namely 1, 10 and 30 mm (3 h)⁻¹. Starting from the analysis of the 0–3 h forecast, the following four points can be noticed: (a) the GNSS-ZTD assimilation improves the precipitation forecast for all precipitation thresholds as the GNSS symbols are always closer to the upper-right corner compared to the CTRL; (b) the FBIAS of the CTRL is underestimated for all precipitation thresholds; (c) the FBIAS, POD, FAR and TS of the GNSS are always larger than the corresponding CTRL; (d) the performance of both GNSS and CTRL decreases with the increasing precipitation thresholds. Point (a) shows that the precipitation forecast is improved by GNSS-ZTD data assimilation; points (b) and (c) show that the assimilation of GNSS-ZTD improves the CTRL underestimation of the precipitation events for all thresholds. Point (c) shows that, while the performance of GNSS is improved compared to CTRL, the added value of GNSS-ZTD data assimilation is reduced by the FAR increase. This point will be confirmed by the results of the statistical test shown later in this section. Point (d) shows the well-known difficulties of correctly predicting the correct number, location and timing of precipitation events as their intensity increases.

Figure 9 Performance diagram showing precipitation scores over the whole period for CTRL (squares) and GNSS (circles) for the 0–3 h period (a) and for the 3–6 h period (b) after the assimilation phase for three different precipitation thresholds: 1 mm (3 h)⁻¹ (in red), 10 mm (3 h)⁻¹ (in green) and 30 mm (3 h)⁻¹ (in magenta).

Results for the 3–6 h phase are like those of the 0–3 h period, highlighted in the four points above, with the exception that the FAR for the 30 mm (3 h)⁻¹ threshold is reduced by GNSS-ZTD data assimilation. In addition, the comparison of Fig. 9a and b shows that the improvement of the model performance decreases with increasing forecasting time because the GNSS and CTRL symbols are closer for the 3–6 h phase compared to the 0–3 h forecast. This is in agreement with the analysis of the PWV presented in the previous section.

All in all, there are two important points to remark upon for the 0–3 and 3–6 h performance diagrams: (a) the performance of CTRL is improved by GNSS-ZTD data assimilation; (b) the impact of assimilating GNSS-ZTD on the precipitation forecast lasts at least 6 h.

To better examine the difference in the statistical scores of the CTRL and GNSS forecasts shown in Fig. 9, the results of a bootstrap statistical test (see Appendix B for details) are shown in Table 2 for the 0–3 and 3–6 h forecast phases. The level of statistical significance is shown in plain text (90 %, 95 %, 99 %) if the GNSS simulations perform better than CTRL and in bold for the opposite case. The statistical test was computed for 1 and 5 mm (3 h)⁻¹ and then every 5 mm (3 h)⁻¹ up to 50 mm (3 h)⁻¹ thresholds. Results are shown up to 30 mm (3 h)⁻¹ thresholds because differences were not statistically significant for larger precipitation events. From Table 2, it is apparent that the assimilation of GNSS-ZTD has a good and significant impact on the precipitation forecast because the FBIAS and POD are improved for several thresholds. The FAR is worsened by the data assimilation because of the larger number of false alarms predicted by GNSS compared to CTRL. This limits the improvement of the ETS, which accounts for both hits and false alarms; however, the ETS has statistically significant improvements at 25 and 30 mm (3 h)⁻¹, showing a significant and positive impact of GNSS-ZTD data assimilation for moderate–intense precipitation forecast. Finally the improvement of the FBIAS, POD and ETS scores for both the 0–3 and the 3–6 h phases is noted, showing that the impact of GNSS-ZTD data assimilation lasts at least 6 h into the forecast.

Table 2 Results of the resampling statistical test. Statistical significance of 90 %, 95 % and 99 % is shown in plain-text when GNSS-ZTD data assimilation improves the statistics, while it is shown in bold when the GNSS-ZTD data assimilation has a negative impact on the statistics. Significant values less than 90 % are not shown. The first number in each cell refers to the first 3 h of forecast (0–3 h), while the second number refers to the second 3 h of forecasts (3–6 h). Significant values less than 90 % are represented by the forward slash symbol.

4 Sensitivity tests

4.1 Data-thinning experiment

As discussed in Sect. 2, the ZTD observation-error covariances are assumed to be uncorrelated in space and the R matrix is diagonal. However, the ZTD observation errors are correlated to some degree and the assimilation of GNSS-ZTD data from receivers that are too close is sub-optimal. To consider in more detail this issue, we performed a sensitivity test for 16 d of October 2019 and data thinning. The days selected are those from 14 to 23 October that, as better discussed in the Supplement of this paper, were characterized by heavy and widespread rainfall, especially in the NW of Italy, and those from 5 to 10 October, which were also characterized by moderate to intense rainfall. In this experiment (THIN), we first derive a distance where the observation errors become uncorrelated (d) and then GNSS-ZTD observations are assimilated only if their distance is larger than the d value. To estimate this distance, we use the method of Desroziers et al. (2005; see also Bennit et al., 2017).

In this method, observation minus background (innovation) and analysis minus observations (residual) statistics are used. To compute the observation-error decorrelation distance, GNSS receivers are binned every 12.5 km and the covariance is computed as a function of the distance:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{cov}=& {\displaystyle \frac{\mathrm{1}}{{\sum }_{i=\mathrm{1}}^n{m}_i}}{\sum }_{i=\mathrm{1}}^n{\sum }_{i={\mathrm{1}}_i^m}\\ & \left[{\mathit{y}}_i-H_i\left({\mathit{x}}_{\mathrm{a}}\right)\right]\left[{\mathit{y}}_j-H_j\left({\mathit{x}}_{\mathrm{b}}\right)\right],\end{array}\end{array}$$

where m_i is the number of GNSS receivers that fall into a specific bin at a certain distance from the ith receiver, n is the number of GNSS receivers and H_i is the observation operator. By Eq. (3) the errors covariance is computed as a function of the distance. Another parameter needed to assess the observation-error correlation distance d is the observation standard error variance, which is the covariance for the zero distance. It is computed as

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n\left[{\mathit{y}}_i-H_i\left({\mathit{x}}_{\mathrm{a}}\right)\right]\left[{\mathit{y}}_i-H_i\left({\mathit{x}}_{\mathrm{b}}\right)\right].\end{array}$$

Figure 10a shows the behavior of the error covariance as a function of the distance (6.25 km, which is the central distance of the first bin; 18.75 km, which is the distance of the second bin; etc.). The values of covariances have been normalized by the variance s² (Eq. 4). In statistics, correlation values below 0.2 are usually assumed to be negligible, and, from Fig. 10a, the 20 km distance can be reasonably assumed to be the distance for observation-error decorrelation. The assimilation of two receivers whose distance is less than 20 km is considered sub-optimal and is avoided in the THIN experiment.

Figure 10 (a) Error covariance as a function of the distance (red curve). Covariances are normalized with the variance s². The y= 0.2 curve is plotted for reference. (b) GPS receivers used for the forecast verification; (c) GPS receivers used for data assimilation.

To obtain GNSS receivers with about 20 km spacing, we divided the domain (35–50^∘ N and 5–20^∘ E) with a regular latitude–longitude lattice with 0.25^∘ spacing and we selected one GNSS receiver for each grid box. If more receivers are present in the grid box, the one with the lowest topographical difference from the WRF model is selected. The number of GNSS receivers used in the data-thinning experiment is 259 (Fig. 10b).

Figure 11 shows the FBIAS, ETS, POD and FAR for the CTRL, GNSS, and THIN and BIAS experiments (the last is considered in the following section) for the 0–3 h forecast phase. Here we limit the discussion to this phase as the differences between THIN and BIAS are lower for the 3–6 h phase. Considering the four scores, it is apparent that the difference between THIN and GNSS is small, especially for precipitation thresholds larger than 40 mm (3 h)⁻¹. For lower rainfall thresholds, the THIN experiment is slightly worse. All in all, these results do not show an improvement for the data-thinning experiment. This could be related to several factors, and three of them seem more relevant: (a) the limited period considered, (b) the importance of the local scale in water vapor distribution and (c) the fact that the Desroziers method is an estimate of the optimal distance between two receivers to be assimilated. So, considering the results of the data-thinning experiment, further research must be done to obtain a definitive answer.

Figure 11 (a) The FBIAS for the CTRL, GNSS, THIN and BIAS experiments; (b) as in (a) for the ETS; (c) as in (a) for the POD; (d) as in (a) for the FAR. Scores are computed for the 1 mm (3 h)⁻¹ threshold and every 2 mm (3 h)⁻¹ from 2 to 60 mm (3 h)⁻¹. Scores are computed for the 10 d between 14 and 23 October 2019.

4.2 Impact of the bias removal

It is interesting to study the impact of the bias removal on the precipitation forecast. The 3D-Var system assumes that observations are unbiased with respect to the WRF model, and the bias was removed before data assimilation as discussed in Sect. 2.1. This approach assumes that we know the observations for the forecasted period, which is not possible in an operational context. While other methods are possible (Bennit and Jupp, 2017; Giannaros et al., 2020; Benjamin et al., 2016, among others), it is expected that in our case the impact of the bias removal is small because we limited the data assimilation to GNSS receivers with a maximum height difference from the model orography of 300 m and because WRFDA takes into account the difference between the GNSS receiver and model height. To better assess this point, we compared the results of simulations with or without bias removal for the 10 d between 14 and 23 October 2019. The experiment without bias removal is named BIAS (Fig. 11). As expected, results show small differences between BIAS and GNSS experiments with the latter performing slightly better.

5 Conclusions

In this paper we studied the impact of GNSS-ZTD data assimilation on the short-term (up to 6 h) forecast over Italy for the month of October 2019, using the WRF model. A dense dataset of 388 GNSS receivers was used.

The comparison between first guess and ZTD observations showed that the forecast without data assimilation underestimates the water vapor content for the period. The GNSS-ZTD data assimilation partially compensates for this underestimation, increasing the water vapor content in the atmosphere. The data assimilation roughly halves both the BIAS and the RMSE statistics for the ZTD. The analysis over the Italian territory shows a general reduction in BIAS and RMSE of the PWV thanks to the rather homogeneous and widespread coverage of GNSS receivers.

A case study on 15 October 2019 was chosen to show the impact of assimilating GNSS-ZTD on the precipitation forecast. The analysis at 18:00 UTC shows that 3D-Var spreads the water vapor innovations both horizontally and vertically, and the improvement of the short-term precipitation forecast is notable, the main drawback being the increase in false alarms over central Italy.

Considering the statistics over the whole period, the analysis of the PWV forecast shows a clear improvement for the simulations with data assimilation. This improvement is widespread over the Italian territory. The PWV RMSE is almost halved for the first forecasting time. As expected, the improvement of the PWV RMSE decreases with forecasting time as the effects caused by data assimilation are partially advected out of Italy.

Assimilating GNSS-ZTD increases the precipitation of the short-term forecast compared to CTRL, which shows a systematic underestimation of the FBIAS. Two 3 h periods were considered after the last analysis time: 0–3 and 3–6 h. For both periods, the FBIAS and the POD are increased by data assimilation and the control forecast is improved. As a drawback, the number of false alarms was increased by GNSS-ZTD data assimilation. The ETS is also improved by data assimilation, even if the increase in the number of false alarms limits the impact of data assimilation on this score.

The results of the resampling statistical test show that the improvement for the FBIAS and POD is statistically significant for several thresholds up to 30 mm (3 h)⁻¹. The FAR increase is statistically significant up to 20 mm (3 h)⁻¹, while the ETS has significant improvements for 1, 25 and 30 mm (3 h)⁻¹.

We started to evaluate the impact of data thinning and of the bias removal in this work. The impact of bias removal is expected to be small because we assimilated data of the GNSS receiver whose height difference from the model does not exceed 300 m. This was confirmed by the results of a numerical experiment spanning 10 heavy-precipitation days. Also the data-thinning experiment did not show consistent differences between the GNSS and THIN experiments. We believe that a longer numerical experiment must be considered to better assess this point.

All in all, the statistical analysis reveals a positive impact of GNSS-ZTD and the improvement is apparent up to 6 h.

While the results of this paper are encouraging, there are several points that need future studies. First, the Mediterranean climate has important seasonal variability and the impact of assimilating GNSS-ZTD must be studied in different seasons. Second, this study refers to only 1 month: longer periods must be considered to give a more robust assessment of the impact of GNSS-ZTD data assimilation on the forecast over Italy. Third, other techniques to assimilate GNSS data, such as assimilating the precipitable water vapor (PWV), and the sensitivity of the results to the background error matrix should be considered in future studies. Fourth, assimilating GNSS-ZTD in real time over Italy and assimilating the ZTD gradients are two subjects that deserve detailed future research.

Appendix A: Precipitation scores

Scores are computed starting from contingency tables (Table A1) for dichotomous events, i.e., events that can have only two values. In this case the two values are “yes” or “no” and the event is “precipitation is above or equal to a certain threshold”.

Table A1 Contingency table for dichotomous events.

In Table A1, $a,b,c$ and d have the following meaning:

• a represents the hits. A hit occurs when both the precipitation forecast and the corresponding rain gauge observation are above or equal to a rainfall threshold.

• b represents the false alarms. A false alarm occurs when the precipitation forecast is above or equal to a rainfall threshold, while the corresponding rain gauge observation is below the same threshold.

• c represents the misses, i.e., when the precipitation forecast is below a rainfall threshold, while the corresponding rain gauge observation is above or equal to the same threshold.

• d represents the correct no forecasts, i.e., when both the precipitation forecast and the corresponding observation are below a rainfall threshold.

The rain gauges used for computing the elements of the contingency tables are shown in Fig. A1. Not all rain gauges report data for each 3 h time interval considered, and the average number of reports is 2700.

Figure A1 Rain gauge elevations [m].

Starting from the values in Table A1, the following scores can be defined:

• Frequency bias (FBIAS).

  $$\begin{array}{}\text{(A1)}& \mathrm{FBIAS}={\displaystyle \frac{a+b}{a+c}},\end{array}$$

  which represents the frequency of the predicted events above a rainfall threshold with respect to observed events and can assume values between 0 and ∞, with 1 being the best value possible.

• Probability of detection (POD).

  $$\begin{array}{}\text{(A2)}& \mathrm{POD}={\displaystyle \frac{a}{a+c}},\end{array}$$

  which is the ratio between the number of correctly predicted events and the number of observed events and can assume values between 0 and 1, with 1 being the best value possible.

• Threat score (TS).

  $$\begin{array}{}\text{(A3)}& \mathrm{TS}={\displaystyle \frac{a}{a+b+c}},\end{array}$$

  which is given by the ratio between the number of events correctly predicted and the sum of observed and predicted events. The TS assumes values between 0 and 1, and 1 is the best value that can be obtained.

• Equitable threat score (ETS).

  $$\begin{array}{}\text{(A4)}& \mathrm{ETS}={\displaystyle \frac{a-a_{\mathrm{r}}}{a+b+c-a_{\mathrm{r}}}},\mathrm{where}{a}_{\mathrm{r}}={\displaystyle \frac{\left(a+b\right)\left(a+c\right)}{a+b+c+d}}.\end{array}$$

  The ETS is similar to the TS but takes into account the possibility of correctly forecasting an event by chance. It varies in the range between $-\mathrm{1}/\mathrm{3}$ and 1, with 1 being the best value. A value of zero represents a useless forecast in which the probability of correctly predicting an event is random.

• False alarm rate (FAR).

  $$\begin{array}{}\text{(A5)}& \mathrm{FAR}={\displaystyle \frac{b}{a+b}},\end{array}$$

  which is the ratio between false alarms and the number of predicted events and can assume values between 0 and 1. In this case the best value possible is 0. Scores were computed after summing contingency table elements for each model over all simulations (4 × 30 = 120 times in this paper, for both the 0–3 and the 3–6 h forecast intervals).

Appendix B: The resampling method

The resampling method from Hamill (1999) is used for assessing if score differences are statistically significant in a confidence interval. For this purpose, a hypothesis test is performed. The null hypothesis is that the scores of the considered models do not differ.

In our case, the null hypothesis is that the differences between the two model scores, i.e., CTRL and GNSS, are zero. Let S₁ and S₂ be a generic score, namely the FBIAS, ETS, POD and FAR (see Appendix A for details), for the two model types, CTRL and GNSS, respectively. The null hypothesis can be then written as

$$\begin{array}{}\text{(B1)}& H_{\mathrm{0}}:S_{\mathrm{1}}-S_{\mathrm{2}}=\mathrm{0}.\end{array}$$

The test statistic is calculated after summing contingency table elements for each model over all simulations (4 × 30 = 120 times in this paper, for both the 0–3 and the 3–6 h forecast intervals). The contingency tables can be written as a vector:

$$\begin{array}{}\text{(B2)}& {\mathit{x}}_{i,j}={\left(a,b,c,d\right)}_{i,j},\end{array}$$

where i is the model type (CTRL or GNSS) and j= 1, …, 120 is the contingency table for each time interval.

The contingency table elements are then summed over all times for both model forecasts CTRL and GNSS,

$$\begin{array}{}\text{(B3)}& S_i=f\left(\sum _{j=\mathrm{1}}^n{x}_{i,j}\right),\end{array}$$

and the test statistic is given by the difference between S1 and S2.

Resampled sums can be calculated after randomly choosing either one or the other model for each time. Let I_j be a random number which can assume values of 1 and 2 for CTRL and GNSS, respectively, with j= 1, …, 120. We can then sum the shuffled vectors over all times,

$$\begin{array}{}\text{(B4)}& S_{\mathrm{1}}^{\ast }=f\left(\sum _{j=\mathrm{1}}^n{x}_{Ij,j}\right),\end{array}$$

and sum separately data not selected for the first sum, given by the index (3 – I_j):

$$\begin{array}{}\text{(B5)}& S_{\mathrm{2}}^{\ast }=f\left({\sum }_{j=\mathrm{1}}^n{x}_{\left(\mathrm{3}-Ij\right),j}\right).\end{array}$$

Each sample is produced to be consistent with a null distribution; i.e., score differences ($S_{\mathrm{1}}^{\ast }-$ $S_{\mathrm{2}}^{\ast }$) are zero.

This random sampling is performed many times (10 000).

We consider the significance levels of α= 0.01, α= 0.05 and α= 0.1, and we test the null hypothesis H₀ with the percentiles. Let t_L be the $\frac{\mathit{\alpha }}{\mathrm{2}}$ percentile of S₁−S₂ distribution and t_U be the $\frac{\left(\mathrm{1}-\mathit{\alpha }\right)}{\mathrm{2}}$ percentile of S₁−S₂ distribution; the null hypothesis H₀ is rejected at the level of 90 % (α= 0.1), 95 % (α= 0.05) or 99 % (α= 0.01) if

$$\begin{array}{}\text{(B6)}& S_{\mathrm{1}}-S_{\mathrm{2}}<t_L\end{array}$$

or

$$\begin{array}{}\text{(B7)}& S_{\mathrm{1}}-S_{\mathrm{2}}>t_U,\end{array}$$

where S₁ and S₂ are the non-resampled scores.