Numerical simulation of rainfall with assimilation of conventional and GPS observations over north of Iran

In this work, the effect of assimilation of synoptic, radiosonde and ground-based GPS precipitable water vapor (PWV) data has been investigated on the short-term prediction of precipitation, vertical relative humidity and PWV fields over north of Iran. We selected two rainfall events (i.e. February 1, 2014, and September 17, 2014) caused by synoptic systems affecting the southern coasts of the Caspian Sea. These systems are often associated with a shallow and cold high pressure located over Russia that extends towards the southern Caspian Sea. The three dimensional variational (3DVAR) data assimilation system of the weather research and forecasting (WRF) model is used in two rainfall cases. In each case, three numerical experiments, namely CTRL, CONVDA and GPSCONVDA, are performed. The CTRL experiment uses the global analysis as the initial and boundary conditions of the model. In the second experiment, surface and radiosonde observations are inserted into the model. Finally, the GPSCONVDA experiment uses the GPS PWV data in the assimilation process in addition to the conventional observations. It is found that in CONVDA experiment, the mean absolute error (MAE) of the accumulated precipitation is reduced about 5 and 13 percent in 24h model simulation of February and September cases, respectively, when compared to CTRL. Also, the results in both cases suggest that the assimilation of GPS data has the greatest impact on model PWV simulations, with maximum root mean squares error (RMSE) reduction of 0.7 mm. In the GPSCONVDA experiment, comparison of the vertical profiles of 12h simulated relative humidity with the corresponding radiosonde observations shows a slight improvement in the lower levels.


Introduction
The ground-based global positioning system (GPS) is used as a new technique of remote sensing of the atmospheric water vapor since more than two decades.The GPS signals passing through the atmosphere encounter delay due to the atmospheric contents especially water vapor.The total delay is divided into two parts: wet and dry.Using surface temperature and pressure measurements the wet delay is mapped into the precipitable water vapor (PWV) [Bevis et al. 1992, Rocken et al. 1993].
Sub hourly measurements of the GPS PWV are accurate for values of less than 2 mm [Duan et al. 1996, Rocken et al. 1997].The ground-based GPS observations compared with conventional methods of measuring water vapor in the atmosphere such as radiosonde observations, have a high temporal resolution, long-term stability and are independent of all-weather condition [Rocken et al. 1997, Yunck et al. 2000].
The formation of clouds and precipitations are directly related to the distribution of the atmospheric water vapor.Therefore, high-resolution measurements of this quantity can play an important role in numerical weather prediction (NWP) models [Mazany et al. 2002, Zhang et al. 2007].
Accurate initial conditions are necessary to have a better numerical prediction of the atmosphere.Assimilation of conventional (such as synoptic and radiosonde) and non-conventional (such as satellite and groundbased GPS) data can improve the initial conditions [Govindankutty andChandrasekar 2011, Srinivas et al. 2012].The variational data assimilation methods have the ability of direct assimilation of the observed variables such as satellite radiance, radar reflectivity and GPS PWV that are different from the model prognostic variables.Using this type of data assimilation approaches, several studies showed a positive impact in the simulation of the rainfall events [Lipton et al. 1995, Kalnay 2003, Zapotocny et al. 2007].
Many attempts have been made to evaluate the impact of the GPS water vapor observations assimilation on the NWP model outputs [Kuo et al. 1993, Kuo et al. 1996, Falvey and Beavan 2002, Zhang et al. 2007, Zhang et al. 2008, Bauer et al. 2011, Leiming et al. 2012].Kuo et al. [1993Kuo et al. [ , 1996] ] showed that the assimila-
tion of PWV measurements into a NWP model improve the vertical structure of the predicted moisture field.Falvey and Beavan [2002] examined the potential of GPS PWV in the advanced regional prediction system (ARPS) model retrievals of precipitation.Based on their results, assimilation of hourly GPS PWV improved the predicted total precipitation of the model up to 3 percent.Also, the use of GPS PWV data assimilation can eliminate the wrongly predicted rainfall systems and increase the quality of rainfall events predictions in timing, location and intensity [Zhang et al. 2007, Zhang et al. 2008].
In a study, Leiming et al. [2012] assimilated the GPS PWV measurements together with the automatic weather station (AWS) observations into a short-range numerical model system using 3-dimensional variational (3DVAR) data assimilation scheme.They found that the data assimilation of these observations can improve the rainfall and temperature predictions.The Iranian permanent GPS network (IPGN) has been established to study the crustal deformation and also for navigation and surveying purposes.At present, this network consists of about more than 90 active stations which are mostly distributed in the northern part of Iran.So far, no effort has been made in using the observations of IPGN stations to improve the performance of NWP models.The main objective of this study is to present preliminary the results on the impact of 3DVAR assimilation of GPS PWV measurements together with the conventional synoptic and radiosonde observations in the simulation of rainfall over the northern part of Iran using the weather research and forecasting (WRF) model.
In Section 2, a brief description of conventional data, GPS PWV measurements and their estimation process are provided.The description of the numerical experiments, the WRF model settings used in this work and 3DVAR assimilation methodology are presented in Section 3. The result of the three data assimilation experiments (with different selection of observations) in the prediction of rainfall and moisture fields are compared in Section 4. The summary and conclusions are given in the last section.

Conventional data
Most of the heavy precipitations with more than 100 mm/day, occurred in the southern coasts of Caspian Sea are associated with a thermal shallow and cold high pressure system located over Russia and the northern Caspian Sea.Due to the temperature difference between the warm water of Caspian Sea and relatively cold air in the shallow high pressure, intense evaporation occurs and thus air becomes saturated in the boundary layer.The low pressure aloft causes vertical motions and precipitations take place in the southern Caspian Sea.Accurate predictions of precipitations associated with the above mentioned mechanism is a challenging task and has much importance in operational forecasting.Thus in this research, two of such synoptic systems are selected to investigate the impact of the assimilation of radiosonde and synoptic station observations as conventional data and GPS water vapor measurements in numerical prediction of two rainfall cases (February 1, 2014, andSeptember 17, 2014) over the northern part of Iran.The study area (inner domain of simulation) is located between geodetic latitudes of 33.06 to 39.44 and longitude 42.63 to 61.82 degrees.
Conventional observations of temperature, pressure, humidity and wind for 50 synoptic stations in the inner domain and 150 stations in the outer domain considered were used in the assimilation process.Figure 1 shows the spread of the Zagros and Alborz mountain ranges over the western and northern parts of Iran, respectively.Distribution of the synoptic stations over the country is more or less uniform.Three hourly synoptic data collected at surface stations along with the vertical temperature and humidity profiles observed at four radiosonde stations (yellow stars in Figure 1) were obtained from the I. R. of Iran Meteorological Organization (IRIMO) network.The radiosondes are lunched twice a day at 00 and 12 UTC.Total daily precipitation observations from rain gauges in the study area were used to validate the NWP model simulation of precipitation.

GPS PWV data
The GPS PWV values derived from the processing of the ground-based observations of the IPGN network with sub hourly temporal resolution are used to improve the numerical simulation of rainfall.The inner domain of the model is located in the northern part of Iran, mainly because heavy convective rainfall events occur in the northern part of Iran along the coast of Caspian Sea and in many cases NWP models fail to accurately forecast the precipitation over this region.Moreover, most of the ground-based GPS stations are located in the north and north-west parts of the country that lie inside our inner domain.Accurate precipitation forecast over north of Iran is thus a challenging task for professional forecasters.Distribution of GPS stations in the IPGN network over Iran topographic map along with the boundaries of the inner domain in the simulation is depicted in Figure 1.
The total atmospheric delay in the zenith direction and the components of position are estimated for each station.Total tropospheric zenith path delay of GPS signal is divided into two parts, zenith hydrostatic delay (ZHD) and zenith wet delay (ZWD) [Iwabuchi et al. 2000].
Using the Sasstamoinen model [Davis et al. 1985] and surface pressure measurements (P 0 ) at GPS station, the ZHD can be calculated accurately through: (1) where, { is the latitude and H is the height from geoid in kilometers.The total zenith delay of the GPS signals is calculated during the processing of IPGN observations.By subtracting ZHD from the estimated total zenith delay, the ZWD is easily derived.
The ZWD is a function of the water vapor and atmospheric temperature.So, the PWV could be estimated when weighted average of the atmospheric temperature, T m is calculated.The PWV values are related to the ZWD as follows [Bevis et al. 1994]: where, M w (18.0152 gr / mol ) and M d (28.9644 gr / mol ) are the molar masses of water vapor and dry air, respectively.The physical constants k 1 (77.689K / hPa ), k 2 (71.295K / hPa ) and k 3 (375463 K 2 / hPa ) belong to the formula for atmospheric refractivity [Rüeger 2002].

Description of the model
For each rainfall case, three numerical experiments are conducted.The initial and boundary conditions of the first experiment (CTRL) are prepared from the National Center for Environmental Prediction (NCEP) global forecast system (GFS) analysis data.These data are used at a spatial horizontal resolution of 0.5°× 0.5°a nd a time resolution of 3 h.In the second experiment (CONVDA), using 3DVAR data assimilation system, the conventional surface and radiosonde observations are assimilated to generate model initial and boundary conditions.The third experiment (GPSCONVDA) is similar to the second one, but the GPS PWV measurements from the IPGN network are also assimilated into the WRF model.In each case, the model is integrated 54 hours ahead starting from the initial time at 1200 UTC one day before the rainfall event chosen.
The configuration of model includes 41 vertical layers and two nested domains.The grids spacing of . . . .H m.hPa P 0 cos ZHD HD 0 0022768 0 0000015 1 0 00265 2 0 000285 (3) the inner and of the outer domains are 7 and 21 km, respectively.The other details of domains and physical settings of the model are presented in Table 1.

Assimilation methodology
To ingest and assimilate the observational data into the WRF model, we utilized 3DVAR data assimilation technique [Barker et al. 2004].This method provides the optimal estimate of the true atmospheric state via the minimization of the following cost function [Ide et al. 1997]: (4) where, x, x b and y o are the vectors of analysis variable, background forecast and observation, respectively.B is the background error covariance matrix and O is the covariance matrix of the observational error.
In the WRF data assimilation system, the background error covariance matrix B is calculated using the National Meteorological Center (NMC) method [Parish and Derber 1992], which uses the WRF forecasts differences with different ranges and valid at the same time.In this work, the differences of 12h and 24h WRF model forecasts were used and the background error statistics were made with one month forecasts in January 2014 for the domains.
The analysis control variables in WRF 3DVAR system are the amplitudes of empirical orthogonal functions (EOFs) of stream function, pseudo relative humidity, unbalanced part of temperature, velocity potential and surface pressure [Barker et al. 2004].PWV is not a state parameter and the observation operator transforms the model control variables into the PWV counter-parts at the time and location of the observations.The observation operator of the GPS PWV which was used in the WRF 3DVAR system is expressed as follows [Iwabuchi et al. 2005]: (5) where, t ( kg / m 3) is the air density, q( kg / kg ) is specific humidity and dz(m) is the height of the vertical model grid.
The two modes of performing the 3DVAR data assimilation are cold-start and cycling.In cold-start mode, 3DVAR generates the new background analysis only at the WRF starting time of forecast.In cycling mode, the forecasts of the WRF after a specific time window is used as a new background for assimilation instead of the background produced from GFS [Guo et al. 2004, Routray et al. 2010, Hsiao et al. 2012, Eiserloh 2014].In this study, we use WRF 3DVAR data assimilation in cycling mode to account for the spin-up time of the forecasts.So, the forecasts of the WRF are used as a first guess for subsequent cycles with the 3h time window.The structure of the 3 hourly cycling update is depicted in Figure 2. As it can be seen in Figure 2, before the start of the forecast (at 12 UTC) the assimilation is done three times.In other words, the assimilation spin-up time takes 6 hours.

Numerical results
In this section, the 3DVAR data assimilation technique in cycling mode is used to evaluate the impact of the observations on the WRF model simulation in precipitation, PWV and vertical relative humidity fields.

Domains of integration
WSM 3-class simple ice scheme [Hong et al. 2004] Microphysic RRTM scheme [Mlawer et al. 1997 Using the three numerical experiments introduced in Section 3.1, the simulations are conducted for the two rainfall cases previously introduced.

Rainfall
To evaluate the WRF model forecasts of precipitation in each experiment, we used observations of rain gauges in the study area.The mean absolute errors (MAEs) of the accumulated precipitation in 24h simulations (day-1 and day-2) initialized at 12UTC, January 31, 2014 (Case-1) are given in Table 2.As it can be seen, in day-1 the inclusion of conventional data into the WRF model in the CONVDA experiment caused a 5 percent reduction in the accumulated precipitation forecast error as compared with the CTRL.The positive impact of data assimilation on the precipitation forecasts declined in day-2.The MAE values in Case-1 show that there is no significant difference between CONVDA and GPSCONVDA simulated precipitations.So, the assimilation of GPS PWV measurements into the model has almost no effect on the accumulated precipitation forecast in Case-1.Similarly, the MAE values of accumulated rainfall in 24h simulations (day-1 and day-2) initialized at 12UTC, September 16, 2014 (Case-2) are given in Table 3.According to the error values in this rainfall case, the data assimilation reduces the error on the first day of precipitation forecast.The comparison of the three experiments forecast errors shows that assimilation of observational data in CONVDA and GPSCONVDA experiments reduced the accumulative precipitation errors in day-1 up to 13 and 17 percent respectively as compared to the CTRL.Therefore, the inclusion of PWV measurements obtained from IPGN stations into assimilation cycles decreases the rainfall simulation error of day-1 up to 4 percent.Also, the positive impact of conventional and non-conventional data assimilation is declined on the second day of the forecast.

Precipitable water vapor
The simulated PWV values from different experiments during the two rainfall cases are evaluated in this section.For this purpose, the estimated PWV values prepared from the IPGN GPS (using Equation ( 2)) stations located in the inner domain of the model are considered as reference and used to verify the corresponding simulated values of PWV in CTRL, CONVDA and GP-SCONVDA experiments.Figure 3 shows the six hourly root mean squares error (RMSE) of the PWV obtained from the experiments in Case-1 (Figure 3a) and Case-2 (Figure 3a).The RMSE values have been calculated for all of the GPS stations located in the study area and are shown for different forecast lengths in this figure .For each case study, GPSCONVDA and CTRL experiment simulated PWV have the minimum and maximum RMSE over the forecast length, respectively.The assimilation of conventional surface observations in comparison with the CTRL has almost no positive effect on PWV simulation in Case-1 while in Case-2 after 12 h from the start of WRF forecast, the CONVDA reduced the forecast error up to 0.2 mm.The superiority of GPSCONVDA with respect to the CONVDA is clear in the first 24h and 12h forecast lengths of Case-1 and Case-2, respectively.As it can be seen in Figure 3, it is found that the GPSCONVDA has reduced the PWV simulation error from 0.2 to 0.7 mm.In other words, during the first 12h and 24h forecast lengths the assimilation of GPS observations improved the PWV prediction RMSE between 8 to 22 percent.

Vertical relative humidity
It is also important to validate the simulated vertical structure of relative humidity in both rainfall cases.The relative humidity prediction of three different experiments in vertical pressure levels is compared with the available corresponding radiosonde observations at Tehran Mehrabad airport (latitude=35.68°N,longi-tude=51.35°E).The 12h and 24h simulated relative humidity profiles together with the verifying observations (black graphs) are depicted in Figures 4 and 5 for the two cases, respectively.
The results of the model run starting at 12UTC January 31, for relative humidity profiles valid at 00UTC (Figure 4a) and 12UTC (Figure 4b) February 1 indicate that in the first 12 h the simulated vertical structure of moisture from the CONVDA (red graph) and GP-SCONVDA (blue graph) experiments are slightly closer to the verifying observations as compared to the experiment without assimilation.In Case-2, the predicted 12h relative humidity profiles from GPSCONVDA and CTRL experiments (blue and green graphs in Figure 5a) are fairly closer to the observations.Moreover, for the 24h simulation of vertical relative humidity in both cases, there is no significant difference in the results of the three experiments.
From Figures 4a and 5a, it is seen that in general, GPS PWV assimilation may slightly modify the humidity values only in the lower levels without really touching the vertical structure.Even though there is a change in the structure, it is mostly ascribable to the vertical structure of background error covariances.Errors of the background are necessary to give an appropriate weight of contribution to the observations and background during analysis estimation [Bölöni and Horvath 2010].Background error covariance (BE) of the WRF data assimilation system is in Eigen space but this information could be used to extract the BE in physical space [Barker et al. 2004].We use single observation test on the model levels and extract the vertical profile of humidity background error standard deviation at radiosonde location.Figure 6 shows the structure of specific humidity background error standard deviation in the vertical model levels.As seen from Figure 6, the background error standard deviation of humidity decreases from surface to upper levels.It is thus expected that the incorporation of the observations should have the largest impact in the lower levels.
According to Figures 4a and 5a, the three numerical experiments have resulted similar results in terms of relative humidity in the upper levels.The inclusion of the observations of conventional and GPS PWV data is more effective in the lower levels below 500 hPa.This could be probably attributed to the higher values of background standard deviation in lower levels (Figure 6).

Conclusions
In this study, 3DVAR data assimilation scheme for the WRF mesoscale model with two domains at 21 and 7 km horizontal resolutions were used to examine the impact of conventional (surface and radiosonde) observations and GPS PWV measurements as non-conventional observations in the simulation of precipitation, relative humidity and PWV over north of Iran.Using three sets of experiments, the simulations were performed for two rainfall events (February 1, 2014, andSeptember 17, 2014) in the study area.The first experiment named as CTRL run without data assimilation.The other two experiments CONVDA and GPSCONVDA use 3DVAR data assimilation of conventional and GPS PWV along with the conventional observations, respectively.
We investigate the impact of inclusion of observational data using the 3DVAR scheme with respect to the GFS analysis data in generating the initial and boundary conditions used in WRF forecasts.Generally, comparison of the model results with the IRIMO observations indicate a better agreement of the simulations with CONVDA and GPSCONVDA experiments compared with the CTRL simulation as many studies reported earlier [Falvey andBeavan 2002, Leiming et al. 2012].In both rainfall case studies, inclusion of observation into the model improved the predicted cumulative precipitation on the first day of the forecast while the positive impact of data assimilation declines on the second day.Based on the results in Case-1 and Case-2 using assimilation of conventional observations, the MAE of the accumulated precipitation forecast in day-1 reduced up to 5 and 13 percent, respectively.Also, the use of GPS PWV measurements together with conventional data reduced the rainfall simulation error up to 17 percent in the first forecast day of Case-2.

NUMERICAL SIMULATION OF RAINFALL
A comparison of simulated PWV values using different numerical experiments with corresponding values that were estimated from GPS observations showed that forecasts without data assimilation in CTRL experiment led to more errors in prediction compared with the CONVDA and GPSCONVDA experiments.Based on the result, the assimilation of GPS data reduced the RMSE of simulated PWV between 8 to 22 percent during the first 24 h of the forecast length.
To evaluate the effect of assimilation of observations on the prediction of vertical structure of relative humidity, the profiles of 12h and 24h simulated relative humidity were compared with the corresponding radiosonde observations at Tehran Mehrabad station.The 12h simulations of relative humidity from the GP-SCONVDA experiment which used 3DVAR system were slightly closer to the observations only in the lower levels compared to the CTRL experiment.In general, the results showed slightly more improvement at lower levels compared to higher levels.The reason for this is in agreement with higher values of background error variances for humidity at lower levels.
Overall, this study showed a positive impact on local simulations of the precipitation and PWV using surface, radiosonde and GPS observations in the 3DVAR data assimilation system.It should be noted that using a more sophisticated method of 4DVAR may improve the result further.The main reason for this is that the GPS data have high temporal resolutions and have to be used in their own time.

Figure 1 .
Figure 1.Distribution of IPGN stations together with topographic information over Iran.The red rectangle shows the inner domain of model.The yellow stars show the location of radiosondes in the domain.

Figure 3 .
Figure 3. Six hourly RMSEs in predicted PWV using different experiments: CTRL, CONVDA and GPSCONVDA from the WRF forecast initialization in both case studies (a) February 1, 2014 and (b) September 17, 2014.

Figure 6 .
Figure 6.Vertical profile of specific humidity background error standard deviation at Tehran Mehrabad airport station.The profile is obtained from WRF 3DVAR single observation test in model levels.

Table 1 .
Details of domains and physics used in the model.

Table 3 .
The structure of cycling 3DVAR data assimilation with 3h time windows.MAE of the 24h accumulated precipitation (in mm) simulated in different experiments.The WRF initialized at 12UTC, September 16, 2014.