A Synergistic Use of a High-Resolution Numerical Weather Prediction Model and High-Resolution Earth Observation Products to Improve Precipitation Forecast

: The Mediterranean region is frequently struck by severe rainfall events causing numerous casualties and several million euros of damages every year. Thus, improving the forecast accuracy is a fundamental goal to limit social and economic damages. Numerical Weather Prediction (NWP) models are currently able to produce forecasts at the km scale grid spacing but unreliable surface information and a poor knowledge of the initial state of the atmosphere may produce inaccurate simulations of weather phenomena. The STEAM (SaTellite Earth observation for Atmospheric Modelling) project aims to investigate whether Sentinel satellites constellation weather observation data, in combination with Global Navigation Satellite System (GNSS) observations, can be used to better understand and predict with a higher spatio-temporal resolution the atmospheric phenomena resulting in severe weather events. Two heavy rainfall events that occurred in Italy in the autumn of 2017 are studied—a localized and short-lived event and a long-lived one. By assimilating a wide range of Sentinel and GNSS observations in a state-of-the-art NWP model, it is found that the forecasts beneﬁt the most when the model is provided with information on the wind ﬁeld and/or the water vapor content. and M.L.; writing—review and editing, A.P., L.P., M.L., A.N.M., B.R., E.R., A.G., G.T., F.S.M., N.P., G.B, A.M.G. and S.B.; visualization, M.L., A.N.M., S.B.; supervision, A.P. and B.R.; funding acquisition, A.P., L.P., F.S.M., L.L., B.R.


Abstract:
The Mediterranean region is frequently struck by severe rainfall events causing numerous casualties and several million euros of damages every year. Thus, improving the forecast accuracy is a fundamental goal to limit social and economic damages. Numerical Weather Prediction (NWP) models are currently able to produce forecasts at the km scale grid spacing but unreliable surface information and a poor knowledge of the initial state of the atmosphere may produce inaccurate simulations of weather phenomena. The STEAM (SaTellite Earth observation for Atmospheric Modelling) project aims to investigate whether Sentinel satellites constellation weather observation data, in combination with Global Navigation Satellite System (GNSS) observations, can be used to better understand and predict with a higher spatio-temporal resolution the atmospheric phenomena resulting in severe weather events. Two heavy rainfall events that occurred in Italy in the autumn of 2017 are studied-a localized and short-lived event and a long-lived one. By assimilating a wide range of Sentinel and GNSS observations in a state-of-the-art NWP model, it is found that the forecasts benefit the most when the model is provided with information on the wind field and/or the water vapor content.

Introduction
The advance of numerical weather prediction (NWP) models to increasingly higher grid spacings (km and sub-km) has lead to new grounds being reached for what concerns their potential synergy with spaceborne systems. On the one hand, to feed high-resolution (HR) NWP models, HR input data and boundary conditions potentially updatable on, at least, a weekly basis are needed. On the other hand, the present state of the art HR NWP models coincides with the availability of several Earth Observation (EO) data characterized by high spatial and/or temporal resolution, such as the Sentinel missions developed by the European Space Agency (ESA) in the framework of the Copernicus programme (EU Regulation No 377/2014). Copernicus is a space infrastructure conceived to deliver worldwide and timely EO data. The free data policy of Copernicus and the fact that it is a long term programme (all investments are secured for 20 years of lifetime and the second generation of Sentinel satellites is in preparation) makes the use of Copernicus data very attractive, not only for the Remote Sensing community. At present, services like emergency management, land surface and marine applications make use of Sentinel data (https://www.copernicus.eu/en/services), while Copernicus atmospheric services focus on air quality and atmospheric composition, ozone layer and ultra-violet radiation, emissions and surface fluxes, solar radiation and climate forcing but not on NWP applications.
Taking advantage of Copernicus Sentinel data even for HR NWP applications is a challenging but very attractive problem, because Sentinel data have a clear potential to provide additional information for this kind of applications. In fact, Sentinel missions are able to provide high spatio-temporal resolution information over the surface boundary. Moreover, the Synthetic Aperture Radar (SAR) Interfeormetry (InSAR) technique [1,2] applied to Sentinel-1 SAR data enables to retrieve HR spatial-albeit low-temporal-information on the spatially and temporally highly variable atmospheric water vapour fields [3][4][5]. To increase the temporal frequency of EO derived atmospheric water vapour data, Sentinel-1 estimates can be complemented by GNSS (Global Navigation Satellite Systems) data. These are point measurements, thus having a coarser resolution but they reach a temporal resolution of 30 s. NWP models' boundary conditions observations, such as soil moisture, land and sea surface temperature, are currently provided at coarse spatial and/or temporal resolutions, strongly sub-optimal for this application. Poor representation of these boundary conditions is thought to be a significant source of weather forecast uncertainties [6]. It can be expected that ingesting products derived from Sentinel data and/or GNSS measurements into NWP models might significantly reduce severe weather forecast uncertainties.
While some investigations on the assimilation in NWP models of low resolution (tens of km) EO derived products (e.g., soil moisture extracted from the Soil Moisture and Ocean Salinity mission data), are available in the literature [7], few studies were conducted on the ingestion of HR EO products. Only a few papers on the assimilation of water vapour maps derived from InSAR data [8][9][10] and one paper on the ingestion of SAR derived soil moisture maps [11] are available in the literature. These works consider a small spatial domain [8,10] or grid spacing at the cloud-permitting limit (3 km [9]). Hence, to the best of the authors knowledge, there is currently no investigation that combines a big domain with a very high grid spacing (of the order of 1 km, in the cloud resolving range) and considers more than only one EO derived variable. Thus, the ingestion of high spatial and/or temporal resolution EO products (in particular Sentinel ones, that have been used only in few works, such as Reference [9]) in HR NWP models appears to be a significant contribution to the scientific literature.
The scientific questions at the base of the "SaTellite Earth observation for Atmospheric Modelling" (STEAM) project (awarded in 2017 by the European Space Agency (ESA)) concern the synergistic use of high-resolution numerical atmosphere models with space-borne systems. A subset of experiments were presented in Reference [12] with a focus on operational issues, while in this study the complete outcomes of the STEAM project are presented. This paper introduces and discusses the results of a set of experiments on the assimilation in the Weather Research and Forecasting (WRF) model [13], of surface information derived from Sentinel products and information about the atmospheric water vapour content derived from both Sentinel-1 data and GNSS measurements. The objective is to verify whether and to which extent this ingestion produces benefits in terms of the uncertainty in severe weather forecast. In the framework of the STEAM project, we decided to focus on the Mediterranean region because it is frequently struck by severe floods and flash floods causing impressive losses of lives and several million euros of damages every year [14,15]. This kind of extreme events are generally caused by intense rains that result from the complex interactions between the humid and warm air flows coming from the sea and the steep orography that characterizes the coastlines of the basin, for example, Reference [16] and references therein. These intense meteorological events are still challenging from a forecasting point of view because they are often controlled by a combination of relatively small-scale processes (such as surface convergence lines, interaction with the orography or the sea surface temperature, low level moist jets, to cite some), which can be very hard to forecast [17][18][19]. In addition, the Mediterranean region (in particular the Italian territory) is covered by a dense network of rain gauges, allowing for a reliable evaluation of the quantitative precipitation forecast (QPF) quality.
Several simulations are performed to assess the possible benefit of assimilating either a single EO Sentinel derived variable or a combination of them. In the STEAM project framework two test cases are considered, a flash-flood producing storm occurred in the Livorno town (Tuscany, Central Italy) in September 2017 and a precipitation event generating a flood in the Silvi Marina village (Abruzzo, South Italy) in November 2017. The set of experiments includes the open loop (i.e., no ingestion of any variable in the WRF model), used as benchmark and some runs with the ingestion of one variable at a time and the ingestion of a combination of variables.
The paper is organized as follows. Section 2 describes the spatio-temporal evolution of the two selected case studies and the WRF model setup. There is a special emphasis on the cloud-resolving grid spacing nature (1.5 km) of the performed simulations, over 1200 × 1200 km 2 wide area to accommodate large portions of ESA Sentinel snapshots. Section 3 describes the data assimilated and their respective data assimilation procedure. Section 4 presents validation of model outputs and the results discussion. Finally, in Section 5 some conclusions and future works are highlighted.

Case Studies Description
According to the criterion in Reference [20], severe rainfall events in the Mediterranean area can be classified in two main types: type I-long-lived (duration ≥ 12 h) and spatially distributed (more than 50 × 50 km 2 ) and type II-brief and localized, having a shorter duration (≤ 12 h) and a spatial extent smaller than 50 × 50 km 2 . The experiments of this study were conducted by considering one event of type I and one of type II. Events of type I correspond to the equilibrium convection, where it is assumed that CAPE (Convective Available Potential Energy [20]) production by large-scale processes is nearly balanced by its consumption by convective phenomena and, thus, CAPE values stay small. In events of type II, a larger amount of CAPE is available, as a result of building up from large-scale processes over long time-scales. However, the extent to which it produces convection and precipitation is restricted by the need for a trigger to overcome the convective inhibition energy (CIN).
As underlined in the Introduction, Sentinel missions have guaranteed a significant improvement of the temporal resolution of EO data with respect to previous satellite missions (e.g., ENVISAT). However, finding two severe weather events, corresponding to type I and II categories and characterized by suitable Sentinel observations (i.e., performed in a time range compliant with the reliability of the model forecasts ∼48 h, see Reference [11]) turned out to be a difficult task.

The Livorno Event
The first case study (type II) is the Livorno extreme weather event that took place in the night between 9 and 10 September 2017. Starting from the morning of 9 September a large trough deepened on the western Mediterranean, generating an intense flow of mild and extremely humid air from the South, over all the Tyrrhenian sectors and the eastern part of the Ligurian Sea. From the evening of 9 September, the cold airflow associated with vorticity at 500 hPa was supportive of instability conditions on Tuscany region (not shown).
The environment was also conducive to the development of intense persistent local convective precipitation systems, not only because of the slow evolution of the depression area but also because of the wind shear (not shown).
The resulting observed Quantitative Precipitation Estimate (QPE) for the time period 18UTC 9 September 2017-06UTC 10 September 2017 is shown in Figure 1. The map is obtained interpolating the raingauge data using a semivariogram directly estimated from data themselves [21] (courtesy of the Italian Civil Protection Department). The lower panel of Figure 1 highlight the very timely concentrate rainfall event (Type II) with more than 200 mm of rain felt down in three hours. It was also found that an intense lightning activity, indicative of the ongoing convective dynamics, took place (not shown).

The Silvi Marina Event
The second case study (type I) is the Silvi Marina rainfall event that took place on 14-15 November 2017 in central Italy affecting mainly Abruzzo, Marche and Umbria regions. Around 00UTC on 14 November 2017 an isolated trough at 500 hPa appeared over the Italian peninsula, shrinking and slowly moving towards southern Italy 24 h later. The synoptic scale system advected moist and warm air in the lower troposphere (850 hPa) from Africa towards Apennines mountains range in Central Italy.
The resulting observed QPE map for the 24 h of 15 November 2017 is shown in Figure 2 (courtesy of the Italian Civil Protection Department). The raingauge observation (lower panel in Figure 2) confirm that it is a type I event with a constant rainfall falling for more than 12 h during the day with a 15 min rain rate less than 10 mm. In this respect, also the number of lightning strokes observed during the two days of the event was negligible (not shown).

WRF Model Setup
The Weather Research and Forecasting (WRF) model [13] was selected as numerical weather model to accomplish all the experiments carried out within the framework of the STEAM project. It is a compressible non-hydrostatic model with mass-based terrain-following coordinates that was developed at the National Center for Atmospheric Research (NCAR) in collaboration with several institutes and universities for operational weather forecasting and atmospheric science research. In this work the adopted version for WRF is 3.8.1, while for the data assimilation package (WRFDA [22]) is 3.9.1. For this study the three domains setup shown in Figure 3 has been adopted. The domains have, respectively, a horizontal grid spacing of 13.5 (250 × 250 grid points), 4.5 (451 × 450) and 1.5 km (943 × 883) with 50 vertical levels ( Figure 3, all domains top reach 50 hPa). All the simulations are performed with the same set of physical parameterizations, described hereafter, that have already been successfully tested in the study of similar events [19,23,24]. Concerning the surface layer, the MM5 scheme has been adopted [25][26][27]. A convective velocity following [28] is used to enhance surface fluxes of heat and moisture. For the land surface processes, the Rapid Update Cycle (RUC) scheme has a multi-level soil model (6 levels) with higher resolution in the upper soil layer (0, 5, 20, 40, 160, 300 cm) [29,30]. The planetary boundary layer (PBL) dynamics is parameterized with the diagnostic non-local Yonsei University PBL scheme [31]. The WSM6 microphysics six-class scheme is adopted to include graupel and its associated processes [32]. Finally, the radiative processes are parameterized by means of the longwave and shortwave RRTMG schemes [33].
The initial and lateral boundary conditions are derived from NCEP-GFS (National Centers for Environmental Prediction Global Forecast System) analysis and forecast data available at a horizontal grid spacing of 0.25 • × 0.25 • and a time resolution of 3 h, as well as the ECMWF-IFS (European Center for Medium-Range Weather Forecasts Integrated Forecasting System) at a horizontal resolution of 0.125 • × 0.125 • and a time resolution of 3 h.

Data Assimilation Techniques
The data ingestion was performed according to three different methodologies: direct insertion, a nudging-like technique and 3DVAR assimilation [22,34]. With direct insertion we mean that a given variable in the NWP model is substituted with the corresponding one retrieved by EO sensors. This technique was used with the Sea Surface Temperature (SST) and the Land Surface Temperature (LST). The nudging-like technique was applied for the Soil Moisture (SM) observations and included different steps. Firstly, a difference map between the Sentinel-1-derived SM (re-projected over the WRF grid) and the WRF SM from the first level of soil model (in correspondence with the points observed by Sentinel-1) was produced; then the resulting map was interpolated in the horizontal plane to fill the gaps in the Sentinel-1-derived maps and added to the original WRF SM field in the first (superficial) level. At this point, the new SM field was propagated in the underlying vertical levels by linearly interpolating the difference between the observed and simulated SM assuming a difference equal to 0 in the deepest level. This is done because the Sentinel-1 observation is only superficial and in the deepest layers the value of the model is hypothesized to be more reliable.
Wind speed (WS) and wind direction (WD) from Sentinel-1 and Zenith Total Delay (ZTD) from InSAR and GNSS were assimilated using the 3DVAR technique. In particular, for the GNSS a 3 h cycle technique was implemented to take advantage of the high temporal resolution.The 3 h cycling 3DVAR implies each time the update of the background for the current analysis using the 3 h forecast of the previous cycle. The 3DVAR main purpose is to provide an optimal estimate of the atmospheric state, by minimising the cost function of Equation (1) [35] where x is the analysis to be found that minimizes the cost function J(x), x b is the first guess of the NWP model, y 0 is the assimilated observation and y = H(x) is the model-derived observation transformed from the analysis x by the observation operator H for comparison against y 0 . The solution of (1) represents an a posteriori maximum likelihood estimate of the true state given two sources of data: the first guess x b and the observation y 0 . The analysis fit to this data is weighted by the estimates of their errors, contained in B and R, which are the background error covariance matrix and the observation error covariance matrix, respectively. The R matrix is the sum of two distinct error covariance matrices: the observation (instrumental) matrix and the representativity error matrix. This matrix is assumed to be diagonal, as done in most of the models [34] assuming the correlations between different instruments and between different observations made by the same instruments equal to zero. This work adopts the Control Variable option 5 (CV5) of the WRFDA package for the B matrix calculation using the National Meteorological Center (NMC) method [36] (for more details refer to WRFDA User Guide ). The NMC method was applied over the entire month of October 2013 with a 24-hour lead time for the forecasts starting at 00:00 UTC and a 12-hour lead time for the ones initialized at 12:00 UTC of the same day. The differences between the two forecasts (t + 24 and t + 12) valid for the same reference time were used to calculate the domains specific error statistics.

The EO Variables of Interest
As mentioned in the Introduction, this research focuses on Sentinel-derived and GNSS-derived data. The family of missions called Sentinels were specifically developed for the operational needs of the Copernicus programme. Each mission is based on a constellation of two satellites to ensure that the revisit time (less than 6 days) and the spatial coverage of the data are compliant with the requirements of the programme. In this study, Sentinel-1 and Sentinel-3 data were used. Sentinel-1 is a polar-orbiting, all-weather, day-and-night C-band (5.4 GHz) SAR mission for land and ocean services. Sentinel-3 is a multi-instrument mission to measure sea and land surface temperature, ocean colour and land colour with high-end accuracy and reliability.
For a detailed description of the reasoning behind the selection of the EO variables to be ingested into the WRF model, that was based on the Observing Systems Capability Analysis and Review (OSCAR) tool, the reader is referred to Reference [12]. Through OSCAR, which is a resource developed by World Meteorological Organization (WMO) in support of Earth Observation applications, studies and global coordination (https://www.wmo-sat.info/oscar/), a set of requirements for observation of physical variables in support of WMO programs is available.
By giving priority to the availability of data at high spatial resolution among those indicated by OSCAR as the most relevant for NWP applications, we selected the following variables: ZTD (that can be retrieved using SAR interferometric technique from Sentinel-1 or GNSS data), surface soil moisture and surface wind vector over the sea (that can be retrieved using Sentinel-1 data), sea and land surface temperature (that can be retrieved using Sentinel-3 data). Some of these variables (wind vector, sea and land surface temperature) are directly available as Sentinel products, while others (soil moisture and integrated water vapour) have to be retrieved. A description of the methods used to estimate these variables is provided in the Appendices A and B.

GNSS Observations for the Livorno and Silvi Marina Test Cases
ZTD time series at 30 s resolution were estimated for all the geodetic permanent stations within the WRF domain to perform data assimilation experiments. For the selected events the data of 375 permanent stations were available ( Figure 4). The iono-free phase range combination of dual-frequency GPS-only observations for three days including the selected events were processed in PPP mode. An accuracy assessment [37] of the GNSS-derived ZTD time series was performed by comparison with 43 radiosonde launches from eight Italian radiosounding datasets. This comparison yielded a mean difference in ZTD of 0.3 cm and a standard deviation of 1.46 cm (not shown).
GNSS ZTDs are used also to validate the retrieved maps of InSAR ZTD in the Silvi Marina test case. The statistics of the differences for this case are reported in Table 1. A good agreement is reached between GNSS and InSAR after orbital error removal, with a standard deviation of the residuals in the order of few cm.

Sentinel Observations for the Livorno Test Case
The following data were available directly through the Sentinel catalogues: (1) land and sea surface temperatures derived from the Sea and Land Surface Temperature Radiometer (SLSTR) on board of Sentinel-3, generated on a 1 × 1 km 2 grid; (2) wind vector (speed and direction at 10 m a.s.l.) included in the Ocean Wind Field (OWI) component of the Sentinel-1 Level-2 Ocean (OCN) product (again generated on a 1 × 1 km 2 grid). In particular, for the SST, the Sentinel-3 data acquired on 9 September 2017 at 20:36 UTC were gathered, while for LST, the Sentinel-3 data acquired on 9 September 2017 at 09:50 UTC were used. For what concerns SM, WD and WS, Sentinel-1 observed the Livorno area on 8 September 2017 at 17:14 UTC. To apply the multi-temporal SM retrieval algorithm, even the Sentinel-1 data acquired on 2 September, 27 August, 21 August and 15 August 2017 were used. Sentinel-1 ground range detected products were collected to estimate SM; these data were multi-looked (10 × 10), calibrated and geocoded using the 30 m SRTM DEM. Figure 5 shows the maps of soil moisture (a), wind field (b), sea surface temperature (c) and land surface temperature (d). The SM map shows that the area corresponding to the observed QPE maxima is well captured by the Sentinel-1 data and corresponds to quite high SM values, around 0.3 m 3 /m 3 . As for the wind field, the map shows that an intense Scirocco wind was blowing along the southern Tuscany and northern Lazio coastline, while the pattern of the wind over the sea in front of northern Tuscany was more disorganized and chaotic. Additionally, an intense wind jet was blowing near the Strait of Bonifacio, between Corsica and Sardinia. The Sentinel-3 derived SST map shows valid values (without cloud coverage) mainly on the southern part of the area of interest around Sicily Island, while not valid observations are available in front of the Tuscany coastlines. The problem of the scarcity of valid observations when severe weather events occur is common for optical data. For what concerns the Sentinel-3 derived LST data, it can be seen that the northern part of Tuscany, near the area mostly affected by the observed torrential rainfall phenomenon, shows a significantly lower LST than the surrounding areas, especially along the Adriatic coastlines and the eastern portion of Pianura Padana. An inverse square distance weighting interpolation of the Sentinel data on the WRF domains computational grids allowed obtaining the Sentinel maps into to the resolution of the innermost domain (1.5 km) before the assimilation. As mentioned in Section 2.3, the SM maps are also propagated in all the model soil layers (refer to Figure 3 in Reference [12]). Concerning ZTD measurements, for the Livorno test case the ZTD from InSAR was not available, thus only GNSS data were assimilated.

Sentinel Observations for the Silvi Marina Test Case
For the SST, the Sentinel-3 data acquired on 14 November 2017 at 20:24 UTC were used, while cloud free optical data to derive LST were not found for the Silvi Marina event. For what concerns SM, WD and WS, Sentinel-1 observed the area of interest area on 14 November 2017 at 05:10 UTC (descending orbit) and at 17:06 UTC (ascending orbit). To apply the multi-temporal SM retrieval algorithm, even the Sentinel-1 data acquired on 8 November, 2 November, 27 October and 21 October 2017 were gathered (both ascending and descending orbits). Figure 6 shows the available EO observational data.
Looking at the SM maps, the area corresponding to the observed QPE maxima, over Marche and Abruzzo Apennines is well captured by the Sentinel data and it corresponds to SM values around 0.3 m 3 /m 3 . The SST map shows valid values only two days ahead the observed event, mainly on the swat southern part on the Ionian Sea, as well as in front of the Abruzzo and Marche coastline.
Also in this case an inverse square distance weighting interpolation of the Sentinel data on the WRF domains computational grids allowed to obtain the Sentinel maps (SM, SST, WIND) at the resolution of the innermost domain (1.5 km) before the assimilation.
For the Silvi Marina case, also ZTD map retrieved from Sentinel-1 is available. In this work 49 Sentinel-1 images have been downloaded and processed by SqueeSAR [38]. Details of the images processed are reported in Table 2.
Atmospheric products have been subsampled, roughly at 100 × 100 [m 2 ]. In Figure 7 an example of the derived atmospheric product is given, depicting the Area of Interest (AOI). Along with the Atmospheric Phase Screen (APS) products a digital terrain model of the AOI used in data processing and an incidence angle map have been generated, to allow the data projection along the zenith local direction.   This variable has been further under sampled for the assimilation due to the fact that the original resolution at 100 m increases the computational cost of the simulation and, most importantly, violates the assumption of spatially independent observation errors made in Section 2.3 for the R matrix [9,10]. For this reason, Reference [10] applied a thinning procedure choosing the nearest available InSAR measurement to each model grid point to obtain a 1 km observed map (same resolution as the inner domain of their setup). With the same method, Reference [9] used a thinning distance of 15 km (for an inner domain resolution of 3 km). To choose the best thinning resolution we performed few sensitivity tests investigating the different observed and modeled ZTD in terms of BIAS and RMSE before and after the assimilation using various InSAR ZTD map resolutions. The thinning resolutions adopted for the sensitivity are 1.5 km, 4.5 km and 13.5 km (the three WRF nested domain grid spacings). For each experiment, the observed point nearest to the center of the corresponding grid cell is kept and the rest is discarded.
Because of the differences between the WRF model orography and the InSAR Digital Elevation Model (DEM), a correction of the ZTD has to be introduced while comparing the WRF ZTD field and the InSAR-derived one. In particular, since the ZTD is the vertical integral of the refractivity N [39], which is a function of the pressure of dry air, the partial pressure of water vapour and the temperature [40][41][42], a simple ZTD correction based on the model variables can be implemented.
Let Z obs be the height of the DEM associated to the ZTD InSAR observation in a specific position. When it is lower than the WRF DEM, Z obs1 < Z 0 , a positive correction using the first model level refractivity N 1 , is introduced. In case the InSAR DEM is higher than the WRF DEM, a negative correction to the model ZTD has to be applied, namely where k 0 corresponds to the last model level that is below the InSAR DEM, Z obs2 . To find the best resolution of the InSAR ZTD to be assimilated, the modeled ZTD calculated on the innermost domain is compared with the InSAR ZTD map resampled on its grid (1.5 km resolution), after the correction due to the height differences between the DEMs described above. The GFS-driven OL (taken as an example) simulation has a BIAS of −0.15 cm and a RMSE value of 2.4 cm with respect to the InSAR ZTD. The simulations assimilating InSAR ZTD at 1.5, 4.5 and 13.5 km resolution have a BIAS of −0.69, −0.27 and −0.35 cm and a RMSE of 1.70, 1.64 and 1.92 cm respectively. Thus, all the simulation assimilating the InSAR ZTD at different resolutions (1.5, 4.5 and 13.5 km) reduce the RMSE but increase the BIAS. Since the RMSE is reduced at most with the assimilation of the InSAR ZTD at 4.5 km resolution, this thinning resolution is chosen for the InSAR assimilation in this case. Future work is planned to estimate the effects on the rain forecast of assimilating the InSAR ZTD maps at different resolutions and with different uncertainties.

Experimental Design
To comply with the objective of the STEAM project (discussed in the Introduction), for extreme events investigation, it was crucial to choose test cases with Sentinel 1 and 3 observations in a reasonable, considering the satellite revisiting times (6 days for Sentinel 1 (every 6 days over Italy there are two passages, one in the moring and another in the evening) and 1-2 days for Sentinel 3), temporal interval before the event. The available observations that follow this requirement are introduced in Sections 3.3 and 3.4. Furthermore, to assess the impact on the forecast of each satellites product it was essential to perform a set of experiment assimilating one variable at a time and afterwards one or more experiments assimilating a combination of the most influencing variables. This produced a set of 7 to 9 experiments (for Silvi Marina and Livorno test case respectively). Table 3 summarizes the experiments methodology, data assimilated (alone or their combinations) and timing for the two use cases. To further consider also the assimilation response to different initialization each set of experiments is performed with two different Global Circulation Models (GCM) namely IFS and GFS. Consequently, for each use case 14 to 18 (for Silvi Marina and Livorno use case respectively) numerical weather modeling experiments were investigated for a total of 32. Table 3 third column summarize the assimilation technique used for each experiment. The main reason of using different assimilation technique deal with the STEAM project aim. For each EO derived variable the most suitable way of assimilation is chosen to assess the impact of Sentinel satellites products assimilation. It was necessary above all to take into account their revisiting time. Furthermore, another important aspect to deal with is the observation type available from Sentinel, how it is modeled and evolved by the WRF model and what was already available in the WRFDA suite for each type of observation assimilation. Thus, concerning wind field and Zenith Total Delay, the 3dvar technique was chosen among the available assimilation technique because it allows both an assimilation to update initial condition when the observation is available only at a given time (i.e., SM for Livorno use case) and a cycling assimilation technique when observations are available with higher temporal resolution (i.e., ZTD from GNSS). Concerning sea surface temperature and land surface temperature, they are boundary observations, so they were ingested through a direct insertion technique. Finally, the soil moisture is again a boundary condition evolved in time by the model with a 1D (vertical) soil model discretized in 6 vertical levels. Thus, it was assimilated with a nudging-like (it doesn't take into account model and observation weights, refer to Section 2.3 for the details).
As it is possible to see from Table 3 the two use cases have slightly different experiments. One of the reason is linked to the absence of LST observation for the Silvi Marina use case and of InSAR ZTD for the Livorno one. While, another different is in the ZTD from GNSS assimilation: for the Silvi Marina use cases only one of the two metodologies applied in the Livorno case is implemented because it was the best experiment setup the Livorno use case and it is also the most compliant one with the R matrix assumption of uncorrelated observation (see Section 2.3).
The Livorno case study is initialized at 18 UTC on 8 September 2017, while the Silvi Marina one is initialized at 00 UTC on 14 November 2017. Both use cases have a 3 h update of the boundary conditions and a 48 h forecasting time interval.

Validation and Results
To validate all the modeling experiments and identify the best performing WRF setup for each case study, the Method for Object-Based Evaluation [43,44] [MODE] is applied by comparing the Quantitative Precipitation Forecast (QPF) of WRF with the raingauges QPE. MODE identifies precipitation structures in both forecast and observed fields and performs a spatial evaluation of the model capability of reproducing the identified observed objects. Such validation method overcomes the so-called "double-penalty" issue [45], that traditional verification methods suffer from. This is particularly true when comparing high-resolution observational data analysis and cloud-resolving meteorological forecasts in case of deep moist convective and highly localized phenomena. Since traditional methods cannot provide a measure of spatial and temporal match between the forecast and the observed rainfall patterns, it is preferable to use feature-based verification techniques, such as MODE. In this paper 13 different indices provided by MODE validation are considered. They include both pairs of objects attributes and classical statistical scores: centroid distance, angle difference, area ratio, symmetric difference and percentile intensity (in this work above the 90th percentile threshold), Frequency BIAS (FBIAS), Probability of Detection Yes (PODY), False Alarm Ratio (FAR), Critical Success Index (CSI), Hanssen and Kuipers discriminant (HK), Heidke Skill Score (HSS). For a complete description of the indices refer to References [23,24]. The setup of all MODE parameters (except the rainfall threshold) is taken as it is by default. This means that in this work it is used a convolution radius of 1 grid square and a merging threshold above or equal to 1.25.
The main goal of this meteorological validation, from a QPF standpoint, is to select the best meteorological forecast out of the whole set of the sensitivity experiments performed to investigate the impact of assimilating Sentinel and GNSS derived products. The outputs of the MODE validation tool are calculated for three different rainfall thresholds, namely 24, 48 and 72 mm. The evaluation metrics to assess the most reliable meteorological forecast in each set of experiments is performed as in Reference [23]: all the indices and statistical scores described above are calculated for each sensitivity experiment, then the times in which a simulation has been the best for each score is counted. Finally, the simulation ranking as the best for the highest number of times is identified.
Additionally, the behavior of the experiments performed is discussed in the following sections both in terms of spatial indices and in terms of cumulated rainfall thresholds, to infer as many physical insights as possible. For a given experiment or a given index, the MODE spatial indices are compared to the OL ones to see if the assimilation has a positive impact on the forecast. The most relevant results are then highlighted and explained on a physical basis. This twofold approach (number of best scores and discussion on the indices with respect to the OL) is suggested as an alternative to the standard MODE approach [43,44]. This is done because, in its original formulation, MODE makes use of an interest function that depends on a certain number of parameters that are not simple to define. Thus, to avoid extracting partial or wrong information, we decided to follow our simple and physically-oriented approach, as in References [12,24].  Figure 8 refers to GFS-driven cases, while Figure 9 refers to IFS-driven cases. The MODE analysis results for the three different thresholds are reported in Table A2 (24 mm), Table A3 (48 mm) and Table A4 (72 mm). In the case of the GFS-driven experiments, the best performing run is WIND+SM+ZTD_GNSS3h_1ist (3DVAR wind at the analysis time, SM assimilation at the analysis time, 3DVAR 3h cycling with the 1-minute ZTD observations temporally nearest to the analysis time) with a total of 15 best scores (Tables A5 and A6). Conversely, in the case of the IFS-driven experiments, the best performing one is WIND (3DVAR at the analysis time, namely 18 UTC on 08 September 2018), with a total of 16 best scores.

Livorno Test Case
Looking at the different rain depth threshold values, the behaviour of the assimilation of different variables can be discussed. For example, the assimilation of SST performs well at the lowest threshold, which is in line with the findings that the average SST controls the larger scale precipitation [46][47][48]. In fact, these works show that a good representation of the SST is important for a good precipitation forecast, indicating, in particular, that the average SST in the area near the convective event controls the total amount of rain (the warmer the sea, the higher the volume of rain). Other recent works, such as Reference [49], show that the small scales in the SST field may affect the surface wind response and, thus, the formation of convergence lines, that can drive upward motion and precipitation [50]. One of the reason why the SST experiments do not show good performances at the high rain intensity might be the fact that SST is retrieved from Sentinel-3 observations in the optical bands. This means that the presence of clouds limits the extent of the observed area, especially in the vicinity of the convective activity. For the higher thresholds, thus, the assimilation of the SST does not produce positive results, probably because of the lack of fine spatial structure in the SST field near the convective activity.
Assimilating the wind field shows very good values of the maximum rainfall values and gets better and better the higher value of the threshold, meaning that the dynamical features are of crucial importance to model the rain peaks. Also the location of the rainfall area gets more realistic, both in terms of mean position and in terms of extent, with respect to the OL when assimilating the wind field.
Assimilating the ZTD only seems to improve the orientation of the rainfall area consistently through the cases analyzed. This is probably related to a more correct spatial distribution of the water vapor field with respect to the OL case. Moreover, an improvement in the maximum rainfall values and in the extent of the rainfall area is observed, also at the higher thresholds.
Among the runs that assimilate only one variable, SST has the highest number of best scores for the low threshold, WIND has the highest number of best scores for the two higher thresholds with GFS; instead, with IFS, the OL has the highest number of best scores for 24 mm and WIND dominates for 48 mm and 72 mm: this confirms the importance of the wind structure in the simulation of an intense rain event.
The rather heterogeneous response to the assimilated Sentinel and GNSS derived variables by the GFS and IFS-driven experiments is most probably due to the difference in their respective initial conditions fields. For example, Figure 10 represents the difference map between the IFS-driven SM experiment and the Sentinel-derived SM (SM IFS -SM S1 , on the right) compared to the difference between the GFS-driven SM experiment and the Sentinel-derived SM (SM GFS -SM S1 , on the left). Furthermore, the SM mean value over the Area Of Interest (AOI) hit by the event has been calculated for the IFS and GFS driven experiments and the Sentinel observation. It is possible to see that the average SM observed by Sentinel over the AOI is closer to the IFS-driven case values with respect to the GFS ones. In fact, Figure 10 reveals bigger differences in the SM GFS -SM S1 picture (on the left) with respect to the SM IFS -SM S1 one (on the right). This means that the assimilation in the GFS experiment produced a larger change with respect to the assimilation in the IFS one, leading to a bigger impact on the forecast. It is also worth noticing that the assimilation of wind, soil moisture and ZTD all together ranks well for both driving Global Circulation Models (GCMs). Thus, for operational purposes, this can be a recommended assimilation setup possibly with a 3-h cycling 3DVAR because this setup allowed to obtain an improvement both in the timing and the intensity of the rainfall peaks, as highlighted in Reference [12] for the GFS-driven experiment.   To gain a deeper insight of the modelling experiments, the best performing IFS-driven case, namely the WIND one, is compared with the OL. In particular, Figure 11 shows the atmospheric state of the two experiments at 02UTC on 10 September 2017, corresponding to the most intense phase of the observed event, displayed with the VAPOR (Visualization and Analysis Platform for Ocean, Atmosphere and Solar Researchers, www.vapor.ucar.edu) software. Over the same scene (Figure 11), 3D isosurfaces (5 × 10 −5 kg/kg) for the rainwater, snow and graupel variables have been rendered in combination with the wind field at 10 m in case of the OL (panel a) and IFS 3DVAR-WIND case (panel c). The results show that over the area interested by the most intense observed rainfall phenomena (blue circle) the OL run is not able to produce any significant convective phenomena, while the IFS 3DVAR-WIND does so. This is also visible in the vertical cross sections of the reflectivity field for the two modelling experiments (green circles panels b and d): while reflectivity is nearly absent nearby the Tuscany coastlines (Livorno area) for the OL experiment (panel b), conversely strong activity is apparent in panel d over the Livorno area and downshear the main convective system, with values peaking up to 50-55 dBz. Similar improvement is achieved also in the GFS-driven case, as shown in Reference [12].
To better understand the impacts of the wind assimilation, the instantaneous 10 m wind field during the intense phase of the event is shown in Figure 12. It is important to notice that one of the most important ingredients responsible of the triggering of this kind of highly-precipitating back-building Mesoscale Convective System (MCS) is the persistence and intensity of the convergence line over the sea in the area in front of the coast hit by the event [19]. Figure 12 shows that IFS OL is not able to produce a significant convective phenomenon over the area hit by the event, because the convergence line over the sea is shifted to the South with respect to the real location for the entire event duration. Conversely, IFS WIND assimilation is able to produce a stronger convergence line with better localization and persistence. Concerning the comparison between IFS and GFS OL simulations, both OL are not able to produce a significant convergence line in the correct position. Furthermore, GFS OL simulation started at 01 UTC to develop a convergence line in front of the Livorno coastline but it rapidly disappeared just one hour later. Instead, the GFS WIND simulation achieved better performance with the presence of the convergence line during the whole event.    Table A7 (24 mm), Table A8 (48 mm) and Table A9 (72 mm). The best performances are achieved with the ZTD assimilation from GNSS observations (ZTD_GNSS3h_1ist), both for GFS and IFS-driven cases (Tables A10 and A11). It is important to highlight that, since this is a type-I event, the large scales are more important in controlling the dynamics with respect to the Livorno case. This is the reason why the OL runs perform quite well, as discussed in what follows.

Silvi Marina Test Case
More specifically, comparing the MODE spatial indices with respect to the OL ones, one can say that, in general, the assimilation of the wind field and/or the columnar water vapor content improve the forecast at all levels of rain intensity. This is interpreted with the fact that all the Sentinel derived fields concur to provide a better localization of the most intense core of the observed event, thanks to a better representation of the spatial distribution of the water vapor and/or the structure of the moisture flux.
As in case of the Livorno experiment, additional important comments can be drafted, when analyzing overall the rain depth threshold values. The distance of the simulated center of mass of the rainfall field with respect to the observed one seems to benefit at most from Sentinel fields when IFS-driven experiments are compared with IFS OL, while the orientation of the rainfall area gets more realistic in general from Sentinel products for lower threshold (24 mm) both in case of GFS and IFS. This statement holds true also for higher threshold in case of GFS-driven experiments, while IFS OL behaves very well for this metric when the most intense core of the rainfall field is targeted, thus making harder to outperform it. Assimilating Sentinel products has in general a positive impacts on the spatial extent both for GFS and IFS-driven experiments, especially at low and intermediate thresholds.
Concerning the overlap between the simulated and the observed rainfall areas, at low (24 mm) and intermediate (48 mm) thresholds, both GFS and IFS-driven experiments exhibit an improvement because of ingestion of Sentinel products. The maximum rainfall value benefits at most from Sentinel products at lower threshold (24 mm) both for GFS and IFS experiments. This is not true for 72 mm threshold but it is worth mentioning how in that case both GFS and IFS OL achieves values well above 0.8 and rather close to the maximum, equal to 1. Finally, among the runs that assimilate only one variable, WIND and ZTD_GNSS3h_1ist have the highest number of best scores at low threshold for both IFS and GFS while, as already argued, ZTD_GNSS3h_1ist guarantees the best performances for higher thresholds. Also the ZTD_INSAR experiment improves the forecast, as indicated by the high number of improved indices with respect to the OL (Tables A7-A9), in particular concerning the rainfall pattern location. This confirms the importance of water vapour data for rainfall events corresponding to convective equilibrium conditions (type I). The substantial difference of the InSAR and GNSS derived ZTD observations is in the spatio-temporal resolution, namely the InSAR ZTD is characterized by high spatial and low temporal resolution, while the opposite holds for the GNSS ZTD signals. Despite the assimilation of both InSAR and GNSS derived ZTD has an added value, the fact that this event is long-lived can be the reason why a constant update of the ZTD every 3 h from GNSS during the entire event allowed to obtain the best result with respect to a single update performed with the Sentinel observations. Furthermore, in this test case both the OL simulations have already a very good performance, thus it is more difficult to obtain a significant improvement from data assimilation.

Conclusions
Two sets of experiments have been presented in order to attempt to answer to the scientific question on whether Sentinel satellite constellation data can be used to better understand and predict severe weather events. For this purpose two different events have been chosen. The first event is the Livorno flash-flood producing storm that, according to the criterion in Reference [20] can be classified as type II-short-lived (duration < 12 h) and very localized (less than 50 × 50 km 2 ). The second one is the Silvi Marina event that can be classified as type I-long-lived (≥12 h) and spatially distributed (more than 50 × 50 km 2 ). For each event two groups of experiments have been performed, the first one driven by IFS global model and the second one by GFS. For each global model, a set of sensitivity experiments has been executed through the assimilation (with different techniques) of a single observed variable, either from Sentinel or from the GNSS network. Then, building on the results of the single-variable assimilation experiments, few simulations assimilating all the most influencing observations together have been performed.
As expected, given their different nature, the two case studies showed different responses to the same type of assimilated observations. In particular, the 3-h cycling 3DVAR of GNSS ZTD offers good results in both Livorno and Silvi-Marina cases but it is not fundamental in the type-II Livorno event that significantly benefits also from the adjustment of the wind field with high-resolution Sentinel observation. Conversely, the type-I Silvi Marina case takes mostly advantage from a 3-h cycling update of GNSS ZTD and the single time Sentinel observations are more penalized because they allow to adjust the model only the day before the event. Furthermore, IFS and GFS-driven experiments revealed sometimes significantly different responses to the same assimilated variable: this is most probably due to the difference in their respective initial conditions fields. For example, comparing the SM analysis fields for the Livorno case, there is a quite prominent difference between GFS and IFS values, namely IFS soil layers look definitely drier than GFS ones. This can lead to different response to the same SM observation. When initialized with GFS initial and boundary conditions, SAR derived SM seems to have an added value with respect to IFS.
In general the assimilation of Sentinel and GNSS derived variables always resulted into an improvement in the weather forecast, even if sometimes relatively small for some variables (like SST or LST). This, as discussed in the previous section, could be linked to the partial coverage of the data that are assimilated in the test cases or the already good initial conditions provided by the global models. Thus, it is worth further exploring the synergy between HR Sentinel observations and weather modelling, maybe in different test cases and geographical contexts.
Concerning the use of water vapor estimates, the assimilation of both InSAR and GNSS derived ZTD observations has positive impacts on the forecast. These results also pave the way to further explorations of the impact on hydro-meteorological predictive capability of data provided by the potential next-generation Earth Explorer geostationary InSAR satellite [51], currently under evaluation by ESA. This satellite would provide ZTD observations at high resolution both in time and space, combining the advantages of both InSAR and GNSS products, highlighted in this work. The algorithm developed to estimate SM was conceived to combine a multi-temporal approach with a good computational efficiency. The latter characteristic is particularly suitable for an operational application, such as the use of SM maps in NWP models. A multi-temporal approach to estimate SM from SAR assumes that the temporal scale of variation of soil roughness is considerably slower than that of SM. Hence, if a dense time-series of SAR data is available, as expected using Sentinel-1, short term changes in the backscattering coefficient σ 0 (that represents the SAR measurement) are basically related to SM variations [52][53][54]. The algorithm is based on a multi-temporal maximum likelihood approach, for example that in Reference [53], used to invert a direct scattering model of backscattering from a bare soil, namely that proposed by Oh [55]. The Oh Model relates the state of the soil, described by the pair (SM, s), where s denotes the height standard deviation of the rough surface, to the σ 0 measured at an incidence angle θ. This relationship is represented in form of lookup table that was generated by applying the Oh model as described in Reference [56]. The retrieval algorithm assumes that a time series of M + 1 measurements (at the current time t and at M previous times t − 1, . . . , t − M) is available. It performs a least square search minimizing the square difference between measured and modelled values of σ 0 : where d is the cost function that has to be minimized and the symbol | dB indicates that the backscatter values are expressed in logarithmic units. In (2), σ 0 VV,soil refers to the soil contribution to the Sentinel-1 measurements, so that, if vegetation is present, its effects on the Sentinel-1 backscattering measurements must be corrected before determining the cost function d. The choice of M is a compromise between SM estimation accuracy and computational time: M = 4 was chosen, as done by Reference [54]. The Water Cloud Model [57], [WCM], with the parameters proposed by Reference [58] for all land uses, was used for the correction of the vegetation effects on σ 0 . More details on the SM retrieval algorithm can be found in Reference [56].
Appendix A.1.1. Zenith Total Delay Retrieval from GNSS The meteorological applications of GNSS are well known [39,41]. A GNSS receiver is able to determine its coordinates in a global reference frame by exploiting the simultaneous observation of its distance from a number of satellites of known position. The distances are derived from the time it takes to the GNSS signals to cover the satellite-receiver paths, by assuming that they travel with the constant velocity of light in vacuum. This assumption results in an observed distance different from the geometrical one by an amount called atmospheric delay. The delay contains an ionospheric component, which can be removed by a proper (iono-free) combination of the dual frequency GNSS signals and a tropospheric component, which is partly due to tropospheric water vapor partly to the dry gases, called wet and dry tropospheric delays respectively. The total tropospheric delay (wet+dry) has to be accounted for in the GNSS data processing to get an accurate positioning. To this aim, the troposheric delays affecting the signals received by a GNSS station from all the satellites simultaneously in view are expressed as a common delay in the zenith direction above the receiver. That is, each slant delay is projected onto the zenith direction by a known mapping function. The common zenith projection is then estimated in the general adjustment of GNSS observations, for instance by modelling it as a random walk process, resulting in a high temporal resolution time series of zenith delays for each considered station. In the STEAM experiments the free and open source software goGPS [59,60] was used to estimate ZTD time series from the GNSS raw observations. goGPS is developed and maintained by GReD srl, it employs a Precise Point Positioning (PPP) strategy to deal with the systematic errors that affect GNSS observations [61], applied trough a multi-epoch joint least squares adjustment.
Appendix A.1.2. Zenith Total Delay Retrieval from InSAR Synthetic Aperture Radar data can be profitably used to monitor the tropospheric water vapor with unprecedented spatial resolution by means of interferometric techniques (e.g., References [9,10] ).
The basic idea behind SAR interferometry is simple. Satellite Radars transmit Electro Magnetic (EM) signals toward the earth, which are scattered by different targets, such as man made structures, bare soil, vegetation, bouncing back toward the satellite. These "back scattered" signals are captured by satellite sensors, at different times, according to their distance from the satellite and are used to generate radar images of the earth's surface.
The key information for ZTD retrieval is the phase of the returned signal, measured for each single target in the scene. This phase is the number of cycles of the EM wave from the satellite to the target and then back to the satellite, plus a contribution that depend on the specific electric and geometric properties of the target.
The phase is proportional to the Optical Path Length (OPL) R from the sensor, located in S(t) (t being the sensing time), to the target, located in P: where v(l; t is the actual velocity of the EM signal through the atmosphere, c is the velocity of light in vacuum, n(l; t) is the refraction index along the path R 0 (P) = |S(P; t) − P| is the sensor-target geometric distance and R N = P S (n(l; t) − 1)dl = 10 −6 P S N(l; t)dl is the extra path due to the variation of the EM signal velocity with in the atmosphere with respect to vacuum, expressed as the integral of the refractivity index N = 10 6 (n − 1). The phase is then: where f 0 the Radar frequency, φ 0 (P) is the target own phase. The phase φ is not usable to estimate the optical path length, R(P; t), since φ 0 (P) is unknown. However, this contribution is removed in the difference between the phases of the same target observed in two different epochs: This is the interferometric phase, that is modelled as the sum of the contributions due difference in the geometrical distances between the target and the satellite (at the two epochs) and the difference in the integrals of the refractivity indexes along the path [62], henceforth leading to the defined Atmospheric Phase Screen (APS). In deriving (A4) from (A2) and (A3) we have assumed that the sensor-target path does not change appreciably, to give a variation of the refractive index.
The retrieval of the APS from (A4) implies: 1-the precise knowledge of ∆R 0 (P) the differential position of the satellite, the orbit and of each target in the scene, with an error much lower than the wavelength (5.7 cm for Sentinel-1), 2-the stability of the target dielectric and geometric properties with time, that is assuming φ 0 (P) constant, 3-the separation of the ionospheric contribution, that is also affecting the propagation delay, 4-the reconstruction of the absolute phase field, from that measured by the data, that is known in modulus 2π.
All these issues has been addressed in these past decades and good solutions have been found, though leading to some limitations. The most challenging is indeed the second, the requirement of the target dielectric and geometrical stability in a time interval of many days. It has been shown (e.g., References [38,63]) the existence of natural targets, named Permanent Scatteres, that are stable over very long periods, encompassing years, so stable that can be used to track sub-millimetric crustal deformations.
The Permanent Scatterer Interferometry (PSI) makes use of long stacks of SAR images, repeating the same geometry, both to identify the stable targets and to separate the two contributions due to APS and long term deformations. This separation is carried out by a spatio-temporal filtering that assume APS low pass in space and incorrelated in time, from crustal deformations, correlated over long term or fast local deformations, like infrastuctures [63]. The APS maps are then generated by converting phases into delays by inverting (A3).
The reconstruction of the absolute phase field needs a further processing step, phase unwrapping, to solve for phase ambiguity [64]. These methods are based on the assumption that the spatial variation of the phases, from target to target, is smooth respect to the ambiguity, that, from A3 is φ a = 2∆R/λ approximately 2.8 cm in the case of Sentinel-1.
As PS are man-made structures but also exposed rocks and stones, their density is pretty large in urbanized areas, ranging from tens to thousands PS per square kilometer, preventing wrapping errors is most cases. However, no reliable measure can be achieved in forest and highly vegetated areas, whereas all water bodies are pure noise, lacking temporal stability.
PSI solves at one time the problem of target stability, by definition of PS, of precise location, by space time filtering and of satellite precise location, by using PS as ground control points, that is items 1,2 and 4 in the list above. The ionospheric contribution, is solved by using frequency diversity, similarly to GNSS, while angular diversity can be used to estimate and remove scintillations [65,66].
In the end, a set of time-differentiated maps can be generated, representing the delay variation, due to troposhere, projected along the satellite line of sight, that corresponds to incidence angles range from 20°to 45°, depending on satellite and acquisition mode.
The characteristics of the APS maps that can be generated by SqueeSAR, a technique that provides consistent improvements respect to PSI over vegetated areas, are summarized in Table A1. Despite the SAR temporal resolution and timeliness do not comply with the OSCAR requirements for Integrated Water Vapor (IWV), it is the only instrument whose horizontal resolution meets the OSCAR requirements at the goal level. Hence, it is worth using the InSAR technique to produce APS maps.

. ZTD Maps Derivation
The differential procedure proper of the interferometric process, see (A4), suppresses any constant term between the two observations. Moreover, the wrapping of the phases in that time interval, generates ambiguities at multiple of 2.8 cm, that limits the knowledge of the differential delays up to an unknown constant.
We can model the Atmospheric Phase Screen measured between the generic i-th image of the stack and one image assumed as reference or "master", as follows: where: -R N (P; t i ), is the troposheric delay in the SAR Line of Sight, LOS, at target P and acquisition epoch t i that we wish to estimate; -t M is the acquisition epoch of the reference master image; -d i the unknown phase constant generated by temporal wrapping discussed above; -OE(P; t i , t M ) is the contribution of the differential orbital error between the master and the i-th image -ν(P; t i , t M ) is an observation noise that we will assume to be independent both in time and space. In order to derive InSAR ZTD maps out of the corresponding R N (P; t i ), i = 1, . . . , N one has to estimate and remove the contribution due to the orbital error and the constant value d i , to project the residual maps onto the zenith direction and to estimate and add back the master delay map. This can be done by using independent estimates of ZTD for each APS epoch as those provided by GACOS (Generic Atmospheric Correction Online Service), an online tool to download ZTD maps based on the ECMWF model [67].        Table A10. Summary of the sensitivity performances (GFS-driven cases). The times in which each forecast has the best result for each score is counted for each threshold and summarized in a total count (summing Tables A7-A9) that is used to find the best simulation.  Table A11. Summary of the sensitivity performances (IFS-driven cases). The times in which each forecast has the best result for each score is counted for each threshold and summarized in a total count (summing Tables A7-A9) that is used to find the best simulation.