2.1.

Test Sites

Four experiments are considered in this study, while two of them are applications of the best possible spin-up time configuration $I$. The four experiments are divided into two pairs and each pair corresponds to one test site. These test sites are now briefly introduced.

2.1.1.

Test site Wettzell in Germany

Wettzell in Germany (WTZR) was selected as a test site because a GNSS station and a corner reflector with accurately known position and a long time series of SAR acquisitions are available. The test site has a topography of medium complexity; it is about 600 m above sea level and located in the Bavarian forest. A digital elevation model is shown in Fig. 1(a).

Fig. 1

SRTM digital elevation model of test sites (a) Wettzell¹⁹ and (b) Mexico City²⁰ with 900-m resolution domain bounding box of WRF. The location, related to the spotty experiments performed in Wettzell, is illustrated by a red pin in (a). The red rectangle in (b) shows the area of the Mexico City test site under consideration for these experiments. It is a clipped area acquired by the Sentinel-1 SAR satellite, as indicated by the surrounding rectangle.

The experiments that took place at this test site were the GNSS and the absolute ranging technique experiment. The methodologies used are introduced in Secs. 2.4 and 2.5, respectively.

2.2.

Weather Research and Forecasting Setup

WRF provides hindcasts to compensate for atmospherically originated wave delays. The ERA-interim reanalysis data are used as initialization data and have a resolution of 0.75 deg. This is about 83 km for the test sites under consideration. Every initialization data set is constructed from the same observation basis and assimilation model; this guarantees equal quality over time. The target resolution of the NWP was 900 m at both test sites, and this resolution was stepwise resolved: the 900-m resolution domain was nested in the 2700-m resolution domain, which was again nested in the 13,500-m resolution domain. These settings overcame the obstacle that between 3000 and 10,000 m resolution there is no clear recommendation for the physical parameter setting. In the vertical direction, 50 layers were used, whereas the layer altitudes were automatically computed by WRF. The top-level altitude was placed at 50 hPa because that altitude is widely used by the WRF community. The domain feedback option was turned off to avoid any unwanted effects related to the different physical parameter settings between the domains. The 900 m resolution domain grid sizes were $403\times 352$ and $202\times 202$ for Mexico City and Wettzell, respectively. The WRF domain setups are shown in Figs. 2(a) and 2(b) for Wettzell and Mexico City, respectively.

Fig. 2

WRF domain setup of (a) Wettzell²¹ and (b) Mexico City²² test sites.

2.3.

Electromagnetic Wave Delay Computation

The predicted delays are derived from NWP. Therefore, the temperature as well as the water vapor data were trilinearly interpolated within the grid cells. For pressure, the horizontal data were bilinearly interpolated and finally exponentially interpolated in the vertical direction.

The distance deviation at time $t$ was estimated by

2.4.

Global Navigation Satellite System Experiment

For the GNSS experiment, two 3-day weather predictions were computed for each of the months January, April, July, and October in order to cover all four seasons. The two simulations that were selected for each month represented one simulation with high and one simulation with low water vapor concentrations in the air. Additionally, the time span during the technique comparison campaign CONT11¹⁶ was simulated. The starting dates were January 1, 2012; January 13, 2012; April 7, 2012; April 27, 2012; July 1, 2012; July 21, 2012; October 3, 2012; October 28, 2012; September 15, 2011; September 18, 2011; September 21, 2012; September 29, 2012; and September 26, 2012. All 13 simulations were initialized at 0:00 GMT and spanned 3 days. Overall, 39 days of atmospheric states were simulated. The three-dimensional state of the simulated atmosphere was sampled every 15 min. Finally, a time series with a time lag of 15 min for the ZPD could be computed—this matched the ground truth GNSS ZPD time series sampling rate. The GNSS ZPD series with a sampling rate of 15 min was provided by E-GVAP.²⁴

Spin-up times of the 900-m resolution are investigated by computing the RMSE $E$ between GNSS observations $z_{\mathrm{g}}(t)$ and spotty predictions $z_{\mathrm{a}}(t)$ depending on the spin-up time $s$ and window size $w$. In doing so, the 6-h window ($w=6\mathrm{h}$) around the spin-up time $s$ that was used selected related residuals and guaranteed a significant statistic. The time-of-day–dependent RMSE affects the spin-up time–dependent RMSE, which is shown in Fig. 6(b). Due to the diurnal temperature and humidity changes, this time-of-day–dependent RMSE was computed and is shown in Fig. 6(a) to distinguish the diurnal effect from the spin-up-time effect.

2.4.1.

Time-of-day–dependent root mean squared error

A time window related to a time of day $d$ selects a subset of the time series, for which the RMSE is computed. Three windows within each simulated time series are present, since every NWP covers 3 days. For clarity, the time windows considered are shown in Fig. 3. This is formalized by

Eq. (3)

$$E_d=\sqrt{\frac{1}{\tilde{n}}\sum _{0\le (t-d)\mathrm{mod}24\times 60\le w}{[z_{\mathrm{g}}(t)-z_{\mathrm{a}}(t)]}^2},$$

where $\mathrm{mod}24\times 60$ describes the relationship to the day time (24 h times 60 min) and $\tilde{n}$ is the amount of selected observations, which are within the 6-h windows.

Fig. 3

Considered time-windows, which are related to the same time of day.

2.4.2.

Spin-up time–dependent root mean squared error

Again, a subset of the time series is selected. The selection window is related to the spin-up time $s$ in this instance. If the investigated spin-up time and the window size equal 6 h, the used time series starts from 6 h and reaches 12 h after the initialization time. As a result, just one window is present in each simulated time series, which is shown in Fig. 4. This is formalized by

Eq. (4)

$$E_s=\sqrt{\frac{1}{\tilde{n}}\sum _{s\le t\le s+w}{[z_{\mathrm{g}}(t)-z_{\mathrm{a}}(t)]}^2},$$

where $\tilde{n}$ is again the number of selected observations. By doing so, the spin-up time accuracy was derived: ranging from 0 to 66 h.

Fig. 4

Considered time-window, which is related to the spin-up time.

2.5.

Absolute Ranging Technique Application

The accurate position of the corner reflector in Wettzell is known and is also estimated by using SAR acquisitions. Therefore, the orbit and the atmospheric range displacement are corrected with respect to a 5- and an 11-h prediction. The bias and the standard deviation with respect to the ground truth position are compared.

2.6.

Spatial Frequency Comparison

The dual-tree complex wavelet transform, which is broadly used in computer vision, is applied for the frequency investigation.²⁵ Basically, it is a wavelet transform with directional wavelets and has some advantages such as shift invariance in comparison to the straight forward discrete wavelet transform.²⁶ Different wavelet domains provide different spatial frequencies and act like a filter in our approach. Frequencies below the NWP resolution are filtered out in order to keep only those frequencies for comparison that are resolved by the interferogram as well as by the predictions. Hence, the sample variance for the wavelet domains $W_1,W_m$

2.7.

Interferogram Correction Application

Remember, APSs disturb interferograms. Correspondingly, predicted APSs from NWP are used for the APS mitigation of interferograms.^3–6 For comparison, four interferograms ${\mathrm{\Phi }}_i$, with $i\in \{1,4\}$, are corrected by 6- and 12-h APS predictions ${\mathrm{\Phi }}_i^6,{\mathrm{\Phi }}_i^{12}$ denoted by

3.1.

Global Navigation Satellite System Experiment

The GNSS experiment provides the best spin-up time configuration $I$. For this, two RMSE figures are derived to distinguish between the day time–dependent effect and the spin-up time–dependent effect.

3.1.1.

Time of day–dependent root mean squared error

All in all, 39 days were simulated. As shown in Fig. 6(a), the 6-h time window of Sec. 2.4 is big enough to generate reliable statistics and small enough to see daily trends.

The following describes the significance of the derived RMSE figure; the histogram of the residuals $z_g(t)-z_a(t)$ is shown in Fig. 5. It shows that the residuals follow a Gaussian-centered distribution. Consequently, the RMSE equals the sample standard deviation $s$ of this distribution. Commonly, the sample standard deviation $s$ is close to the real standard deviation. The confidence interval $u$ is an estimate of the interval that inherits the real standard deviation to a given confidence level ($1-\alpha$) and is given as

Eq. (8)

$$u=[s\sqrt{\frac{\widehat{n}}{{\chi }_{(1-\frac{\alpha }{2};\widehat{n})}^2}};s\sqrt{\frac{\widehat{n}}{{\chi }_{(\frac{\alpha }{2};\widehat{n})}^2}}],$$

where $\widehat{n}$ equals the number of observations and ${\chi }_{(p;\widehat{n})}^2$ is the $p$-quantile of the ${\chi }^2$ distribution with $\widehat{n}$ degrees of freedom. In our case, $\widehat{n}$ is at least 1003, because each RMSE in Fig. 6(a) was derived from at least 1003 observations. The absolute error is defined by the half width of the confidence interval $u$ and provides an estimate of the quality of the RMSE. For the accuracy estimation, a 95% confidence level is assumed in the following. In doing so, the absolute error of each RMSE estimate is better than 0.87 mm.

Fig. 5

Histogram of all residuals between GNSS and WRF ZPDs.

Fig. 6

RMSE of estimated ZPD with respect to (a) time of day and (b) the spin-up time.

3.2.

Absolute Ranging Application

For the absolute ranging technique, the position of a corner reflector is measured within an SAR image, which is influenced by the atmospheric path delay. For improved estimates of the position of the corner reflector in SAR images, the atmospheric delay is predicted from WRF and is used to correct SAR signals. A more precise prediction of the atmospheric delay is equivalent to a more precise estimate of the corner reflector position.

A total of 46 radar acquisitions of the TerraSAR-X (TSX-1)/TanDEM-X (TDX-1) satellites, covering different geometries, were used to derive the bias and the standard deviation of the estimated position. The 5- and the 11-h predictions of WRF provide an accuracy of $-24.0\pm 27.6\mathrm{mm}$ and $-24.0\pm 21.9\mathrm{mm}$, respectively. The bias remains stable, whereas the standard deviation goes down for the 11-h prediction. Due to the small accuracy gain and the relatively small sample count ($n=46$), the improvement is only significant for a confidence level of 0.7. Nevertheless, the accuracy increase of about 21% confirms the observation of the first experiment. Further, the benefit of knowledge about the intersection point $I$ is demonstrated. The main advantage of the 11-h prediction lies in an outlier correction caused by erroneous precipitation prediction. The range offsets corresponding to the outlier are highlighted by the red circles in Fig. 7. Jankov et al.¹⁵ investigated the precipitation depending on the spin-up time and came to the conclusion that too-wet predictions are characteristic for the 6-h spin-up time. A 12-h spin-up time reduces this bias, such that the 11-h prediction performs better for this technique. Accordingly, the authors assume that the region for the 12-h intersection point $I$ can be generalized if similar conditions are present, i.e., in temperate climate regions with high quality initialization data.

Fig. 7

Absolute range displacement for (a) the 5-h prediction and (b) the 11-h prediction.

3.3.

Spatial Frequency Comparison

For the spatial frequency comparison, four short-term interferograms of Mexico City derived from Sentinel-1 data are compared with the predicted APSs (a sample is presented in Fig. 8). To do so, the predicted distance deviations $z_{\mathrm{a}}(t)$ were computed for each range-azimuth position ${\overrightarrow{z}}_0$ of the interferogram. For this application of APS prediction, ${\overrightarrow{z}}_{\mathrm{s}}$ denotes the position of the SAR satellite. A limitation of the wavelet transform is that the number of pixels has to be a power of two with respect to the width and the height of an image. Therefore, only a $149\mathrm{km}\times 143\mathrm{km}$ clipped section of the interferograms are compared. For the predicted APSs, the simulation duration of about 6 and 12 h allows the model to spin-up, i.e., to generate structure while a 900-m horizontal resolution is used. The utilized wavelet is provided by Selesnick²⁷ of length 12 with three vanishing moments and with a specified degree of three, because it estimates the Hurst exponent $H$ of the synthetic data at best in preliminary tests. Sample variances of domains representing larger spatial frequencies than the NWP resolution, of the predicted APS, and of the interferogram are compared to each other. In doing so, high-frequency signals like the digital elevation model error, which could disturb this comparison, are not considered. The resulting scatter plot with respect to different scales and orientations is shown in Fig. 9. The linear relationship is close to 1 in both cases but slightly better for the 12-h predictions. The 6-h prediction shows a structural over- or underestimation of the variances because the dots lay beneath or on top of the red line. This is not the case for the 12-h prediction, which illustrates the advantage of the longer spin-up time and is consistent with the kinetic energy spectra investigation.¹³

Fig. 8

(a) The interferogram of the acquisitions December 2, 2014, and November 8, 2014, east of Mexico City and (b) the corresponding predicted APS using a 12-h spin-up time.

Fig. 9

Scatterplot of variances of wavelet domains between interferograms and simulated APS with (a) 6-h spin-up (b) and 12-h spin-up.

3.4.

Interferogram Correction Application

The same interferograms as in Sec. 3.3 are considered. The covered area is east of Mexico City, and it is known that this area has no linear subsidence.²⁸ The interferometric digital elevation model error signal is related to the spatial baseline, which varies from 8.1 to 45.3 m. Consequently, the standard deviations of the interferograms ${\mathrm{\Phi }}_i$ with small baselines are mainly related to the APS. However, the predicted APS of the NWP has a much lower spatial frequency than the digital elevation model error. Therefore, the derived signal reduction is independent of the high-frequency digital elevation model error, which is present in interferograms with large spatial baselines. The standard deviations of the residuals ${\widehat{\mathrm{\Phi }}}_i^6$ and ${\widehat{\mathrm{\Phi }}}_i^{12}$ are presented in Table 1. For completeness, the related spatial baselines are additionally presented in the last column. This shows that only one interferogram has a large spatial baseline. The 12-h APS predictions always allow a reduction of the APS disturbance. This is shown by comparing the first and third columns of Table 1. In contrast, this is not true for the residuals ${\widehat{\mathrm{\Phi }}}_i^6$, because the second row in this table indicates an additional disturbance for ${\widehat{\mathrm{\Phi }}}_i^6$. On average, the APS mitigation reduces the APS disturbance by about 9% and 18% with respect to the 6-h and 12-h spin-up times, respectively. The accuracy gain results from the closer frequency consistency, which was shown in Sec. 3.3.

Table 1

Standard deviation in mm of the interferograms Φi, the residuals after APS correction derived by 6 h Φ^i6 and 12 h Φ^i12 predictions and the related spatial baselines of the interferograms in m.