## 1.

## Introduction

Synthetic aperture radar (SAR) is a popular remote sensing technique for observing the topography of the earth. The backscattering mechanism of a radar signal is independent of the actual weather conditions. However, the wave propagation velocity depends on water vapor, pressure, and temperature.^{1}

Differential interferometric synthetic aperture radar (DInSAR) images are obtained from the subtracted phase information of two synthetic aperture radar acquisitions and are therefore affected by weather changes. This effect is known as an atmospheric phase screen (APS). The absolute ranging technique applied in the second experiment is also affected.^{2} In recent years, numerical weather prediction (NWP) models became state of the art to mitigate this significant effect independently of the radar data. Different authors have successfully demonstrated the mitigation of the APS using NWP.^{3}4.5.6.7.8.9.^{–}^{10} The method is beneficial for monitoring topographic changes when using the DInSAR or the absolute ranging technique.^{11} Cong et al.^{12} correspondingly demonstrated a straightforward technique exclusively using the ECMWF ERA-interim data, which is a global assimilation data product. In contrast to that technique, the NWP product-related technique enables physical interpolation between time steps based on a physical model of the atmosphere. In the case of ECMWF ERA-interim data, the time sampling is 6 h and a prediction interpolates temporally as well as spatially in a physically correct way.

The ERA-interim NWP input data, which has to be physically interpolated in time and space, is of coarse resolution. The initial data are spatially interpolated at the beginning of the forecast in order to facilitate the higher prediction resolution. This interpolation is not physically balanced and does not inherit any turbulence. Correspondingly, the model needs time to stabilize and to build up structure. This time is called the spin-up time. The imbalance is measured by the Weather Research and Forecasting (WRF) model through the diagnostic variables *dspdt* and *dmudt*, which are the surface pressure– and dry air mass–tendency, respectively. Commonly, both variables decrease rapidly in time and converge against a balanced state. Naturally, the kinetic energy spectra follow a spatial invariant characteristic that is harmed at the beginning because of the missing convection structure.^{13} Gong et al.^{14} investigated the spin-up time for the operational high-resolution rapid refresh NWP model over the Alaska region. Contrary to our results, Gong et al.^{14} did not detect any positive impact of spin-up times larger than 6 h for the mitigation. This might be related to the different climatic region, the different NWP model, or the additional assimilated weather data. In contrast to the investigation by Gong et al.,^{14} this paper focuses on the WRF model and a finer granular spin-up time investigation. Furthermore, the test sites of Wettzell and Mexico City are considered, where water vapor is present at higher levels than in Alaska. It is a matter of fact that the spatial water vapor distribution is affected by the spin-up time, because the kinetic energy spectra are affected by it.^{13} Further, spin-up times that are too short lead to forecasts that are too rainy,^{15} i.e., the prediction inherits too much condensed water vapor and that results in a biased estimate of the corresponding delay. In contrast, if too-long spin-up times are used, the model drifts away from the real atmospheric state and compensates the improvement resulting from the spin-up time. The optimal spin-up time $I$ of both contrary working processes and a positive effect of the spin-up time with respect to the total delay have not yet been reported.

Similar to the interferometric phases from SAR, global navigation satellite system (GNSS) propagation paths are also affected by the atmosphere. The absolute GNSS zenith path delay (ZPD) measurements have a temporal sampling rate of 15 min, and are very accurate and well investigated.^{16}17.^{–}^{18} Due to the good agreement between GNSS and very long baseline interferometry (with mean bias of $-3.4\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ and a standard deviation of 5.1 mm),^{17} the GNSS ZPD provides the ground truth for this analysis. The root mean squared errors (RMSEs) are derived from the ZPD provided by the NWP models. The spin-up time–dependent RMSEs report a 28% accuracy increase for the 12-h prediction when compared to the initialization. The absolute ranging technique then estimates the positions of a corner reflector from radar scenes using the delay predictions from the NWP models. An accuracy improvement of 21% is reported for the 11-h predictions in comparison to the 5-h predictions. Further, spatial frequencies of APS predictions and interferograms are compared. The 6- and 12-h predictions provide different APS predictions. Spatial frequencies of the 12-h predictions show closer consistency with the observed interferogram frequencies than the 6-h predictions at test site Mexico City. In addition, interferograms are corrected by the APS predictions, which report that the 12-h APS prediction mitigation is twice as good as the 6-h APS prediction mitigation.

The first objective of this paper is to report on the 12-h intersection point $I$, observed at Wettzell, that produces the most accurate delay prediction. The second objective lies in reporting on the closer consistency of spatial APS frequencies and the better APS mitigation when considering 12-h APS predictions and interferograms in comparison to 6-h predictions at the Mexico City test site. The better performance of the 12-h predictions is coincident with other WRF investigations.^{13}^{,}^{15} The results of this work can therefore be generalized for a broad spectrum of regions.

## 2.

## Methods

In this section, information about each test site is given first. After this, the WRF setup and the WRF initialization data are described. Then, the electromagnetic wave delay computation is specified. Last but not least, the four experiments to derive the objectives are separately described in detail.

## 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).

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.1.2.

#### Test site Mexico City

The Mexico City test case demonstrates that realistic water vapor distribution prediction results in a realistic APS prediction. The east side of Mexico City was considered, as acquired by the Sentinel-1 SAR satellite at 12:25 GMT. The test site has a tropical climate, and some associated mountains are covered by the SAR acquisitions, as shown in Fig. 1(b). Hence, water vapor effects in the interferograms are highly developed. That is the motivation for using this test site. The master acquisition date of the interferograms that are considered is December 2, 2014, whereas the dates of the slave scenes are October 27, 2014; November 8, 2014; December 14, 2014; and December 26, 2014.

The large beam coverage of Sentinel-1, illustrated by the skewed white rectangle in Fig. 1(b), enables a spatial frequency comparison between the predicted APS and the observed APS of a $149\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{km}\times 143\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{km}$ area. This area is indicated by the red rectangle in Fig. 1(b). The methodology of the comparison is described in Sec. 2.6. In another experiment, the same interferograms are corrected by the predicted APS. A detailed description of this is given in Sec. 2.7.

## 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.

## 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

## (1)

$${z}_{\mathrm{a}}(t)={10}^{-6}{\int}_{{\overrightarrow{\mathbf{z}}}_{0}}^{{\overrightarrow{\mathbf{z}}}_{\mathrm{s}}}{N}_{\mathrm{a}}(\overrightarrow{\mathbf{z}})\mathrm{d}\overrightarrow{\mathbf{z}},$$^{23}

## 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^{16} 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.^{24}

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\text{\hspace{0.17em}\hspace{0.17em}}\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

## (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}},$$## 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

## (4)

$${E}_{s}=\sqrt{\frac{1}{\tilde{n}}\sum _{s\le t\le s+w}{[{z}_{\mathrm{g}}(t)-{z}_{\mathrm{a}}(t)]}^{2}},$$## 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.^{25} 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.^{26} 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}$

## (5)

$$\widehat{\mathrm{V}}({W}_{i})=\frac{1}{2n}\sum _{c\in {W}_{i}}[\mathrm{r}{(c)}^{2}+\mathrm{i}{(c)}^{2}],$$## 2.7.

### Interferogram Correction Application

Remember, APSs disturb interferograms. Correspondingly, predicted APSs from NWP are used for the APS mitigation of interferograms.^{3}4.5.^{–}^{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.

## Results

Four experiments report coincident results, while two of them are applications of the presented objectives. First, the GNSS experiment provides the best spin-up time configuration $I$ for absolute delay mitigation. Second, the absolute ranging technique is an application of this datum $I$ and confirms it. Third, the frequency comparison shows that the best spin-up time configuration $I$ results in realistic APS predictions. Further, this is also the explanation of the delay prediction advantage derived by the best spin-up time configuration $I$. Fourth, the APS mitigation experiment illustrates the second application of $I$ and confirms this result again.

## 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

## (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}}}],$$## 3.1.2.

#### Spin-up time–dependent root mean squared error

A 6-h time window is used to see the spin-up time effect in Fig. 6(b). Accordingly, each RMSE estimate in this figure was derived from at least 284 observations. During the first 12 h, the RMSE estimates are below 15 mm; thus, the error is better than 1.2 mm in this time slot. Of course, the daily trends are also included in Fig. 6(b). Figure 6(a) is generated such that it is comparable to the first 24 spin-up time hours of Fig. 6(b). This means that the used time window for the RMSE computations of Figs. 6(a) and 6(b) is the same for each corresponding point. In Wettzell, an accuracy increase of about 28% (RMSE goes down) is visible within the first 12 h. This is significant because the RMSE decrease is about 4 mm, whereas the error is below 1.2 mm. This effect is related to the spin-up time because the time-of-day–dependent RMSE goes up.

## 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\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ and $-24.0\pm 21.9\text{\hspace{0.17em}\hspace{0.17em}}\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.^{15} 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.

## 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\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{km}\times 143\text{\hspace{0.17em}\hspace{0.17em}}\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^{27} 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.^{13}

## 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.^{28} 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.

Φi | Φ^i6 | Φ^i12 | Spatial baseline |
---|---|---|---|

16.6 | 14.6 | 10.6 | 17.6 |

14.7 | 15.4 | 13.9 | 9.1 |

12.0 | 10.7 | 8.6 | 45.3 |

16.9 | 14.3 | 16.5 | 8.1 |

## 4.

## Conclusion

At Wettzell, a 12-h spin-up time, which reduces the RMSE by about 28% in the zenith direction, is recommendable when compared to the initialization. The accuracy of the absolute ranging technique improves by about 21% and demonstrates the practical benefit of the more precise delay predictions at Wettzell. After this lead-time, the water vapor is physically correctly distributed, because the kinetic energy spectra are well presented and because the outlier likelihood, which was caused by too-rainy predictions, is reduced. Due to the more realistic distribution of the water vapor, the 12-h APS prediction shows a closer spatial frequency consistency than the 6-h prediction compared to interferograms at the Mexico City test site. This results in a better APS mitigation, which is about twice as good for the 12-h APS prediction as for the 6-h APS prediction. Accordingly, it can be assumed that the 12-h spin-up time of WRF provides more realistic APS predictions. The better performance of the 12-h predictions is coincident with other WRF investigations and the applicable region of the found objectives is therefore extended.

## Acknowledgments

We thank E-GVAP, the E-GVAP processing centers, and the GNSS data-owners for access to the E-GVAP NRT GNSS delay data ( http://egvap.dmi.dk). Further, we thank Prof. Dr. Eineder, Prof. Dr. Trautmann, Christian Minet, and Marion Heublein for discussions and proofreading.

## References

## Biography

**Franz-Georg Ulmer** is a PhD student at the Deutsches Zentrum für Luft-und Raumfahrt (DLR). He received his diploma in computer science from the University of Passau in 2010. His current research interests include atmospheric phase screen mitigation for SAR interferometry.

**Ulrich Balss** received his degree in physics and the Dr. phil. nat. degree from the Johann Wolfgang Goethe University, Frankfurt am Main, Germany, in 1994 and 2004, respectively. Between 2003 and 2010, he worked at the Technische Universität München, Munich, Germany. In 2010, he joined DLR’s Remote Sensing Technology Institute. He is involved in the German TerraSAR-X and TanDEM-X missions. In particular, he participated in the development of the TerraSAR-X multimode SAR processor.