## 2 Material and Methods
### Data Assimilation Techniques  
The data ingestion was performed according to three different methodologies: direct insertion, a nudging-like technique and 3 DVAR 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 3 DVAR 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 3 DVAR implies each time the update of the background for the current analysis using the 3 h forecast of the previous cycle. The 3 DVAR main purpose is to provide an optimal estimate of the atmospheric state, by minimising the cost function of Equation (1) [35]  
\[J(x)=\frac{1}{2}\left(\mathbf{x}-\mathbf{x}^{b}\right)^{T}\mathbf{B}^{-1} \left(\mathbf{x}-\mathbf{x}^{b}\right)+\frac{1}{2}\left(\mathbf{y}-\mathbf{ y}^{0}\right)^{T}\mathbf{R}^{-1}\left(\mathbf{y}-\mathbf{y}^{0}\right), \tag{1}\]  
where \(\mathbf{x}\) is the analysis to be found that minimizes the cost function \(J(x)\), \(\mathbf{x}^{b}\) is the first guess of the NWP model, \(\mathbf{y}^{0}\) is the assimilated observation and \(\mathbf{y}=\mathcal{H}(\mathbf{x})\) is the model-derived observation transformed from the analysis \(\mathbf{x}\) by the observation operator \(\mathcal{H}\) for comparison against \(\mathbf{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 \(\mathbf{x}^{b}\) and the observation \(\mathbf{y}^{0}\). The analysis fit to this data is weighted by the estimates of their errors, contained in \(\mathbf{B}\) and \(\mathbf{R}\), which are the background error covariance matrix and the observation error covariance matrix, respectively. The \(\mathbf{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 \(\mathbf{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.