noise is a good statistical representation of the day-to-day variability in the IWV time series around the mean seasonal cycle. The unknown parameters of this model are: $\mu$, the mean IWV, $\omega$ the slope of the linear trend, the $a_i$ and $b_i$ coefficients of the Fourier Series, $\phi$ the autocorrelation of the noise, and $\sigma_e^2$ the variance of the noise. We use a Generalized Least Squares (GLS) algorithm to estimate all the parameters and their formal errors.

## 3. Results

### 3.1. Segmentation Results

This section discusses the segmentation results and how they are impacted by four factors: (1) the GNSS data processing (IGS versus CODE), (2) the length of time series (short period, 1994–2010, and long period, 1994–2018), (3) the reference data set (ERAI versus ERA5), and (4) the auxiliary data used in the ZTD to IWV conversion (ERAI versus ERA5).

Table 2 summarizes the segmentation results for all four comparisons. Statistics are given for 81 common stations in both GNSS data sets. They include average data properties, such as the mean of the estimated monthly variances and the standard deviation of the estimated functional (stdf). The segmentation results are compared by means of the total number of detected change-points, both before and after the screening, the number of outliers, the number of metadata validations, and the number of similar detections.

<table><caption>**Table 2.** Summary of pairwise comparisons of segmentation results from various data sets used in this work. The validation with respect to GNSS metadata and the similar detection statistics used a closeness window of ±62 days. (a) Segmentation is run over the full time series (1994–2018), but change-points are compared for the time-limited period (1994–2010). (b) This CODE data set uses auxiliary data from ERA5. (c) This CODE data set uses auxiliary data from ERAI.</caption><thead><tr><th rowspan="2">Data Set</th><th colspan="2">(1) Impact of Processing</th><th colspan="2">(2) Impact of Time Length</th><th colspan="2">(3) Impact of Reference</th><th colspan="2">(4) Impact of Auxiliary</th></tr><tr><th>IGS—ERAI Time-Matched</th><th>CODE—ERAI Time-Matched</th><th>CODE—ERAI Time-Limited</th><th>CODE—ERAI</th><th>CODE (b) —ERAI</th><th>CODE (b) —ERA5</th><th>CODE (b) —ERA5</th><th>CODE (c) —ERA5</th></tr></thead><tbody><tr><td>Time span</td><td>1995–2010</td><td>1995–2010</td><td>1994–2010</td><td>1994–2018 (a)</td><td>1994–2018</td><td>1994–2018</td><td>1994–2018</td><td>1994–2018</td></tr><tr><td>Mean of the monthly variances (kg m<sup>-2</sup>)</td><td>0.68</td><td>0.62</td><td>0.62</td><td>0.63</td><td>0.61</td><td>0.46</td><td>0.46</td><td>0.46</td></tr><tr><td>Standard deviation of the functional (kg m<sup>-2</sup>)</td><td>0.26</td><td>0.24</td><td>0.24</td><td>0.23</td><td>0.23</td><td>0.17</td><td>0.17</td><td>0.17</td></tr><tr><td>No. detections</td><td>231</td><td>257</td><td>296</td><td>249</td><td>364</td><td>398</td><td>398</td><td>392</td></tr><tr><td>No. outliers</td><td>36</td><td>38</td><td>73</td><td>40</td><td>60</td><td>71</td><td>71</td><td>87</td></tr><tr><td>No. detections after screening</td><td>211</td><td>235</td><td>252</td><td>227</td><td>333</td><td>359</td><td>359</td><td>343</td></tr><tr><td>Validations after screening</td><td>63</td><td>68</td><td>77</td><td>78</td><td>114</td><td>131</td><td>131</td><td>125</td></tr><tr><td>Validations after screening (%)</td><td>29.9</td><td>28.9</td><td>30.6</td><td>34.4</td><td>34.2</td><td>36.5</td><td>36.5</td><td>36.4</td></tr><tr><td>Similar detections</td><td colspan="2">103~48.8%</td><td colspan="2">185~81.5%</td><td colspan="2">151~45.3%</td><td colspan="2">243~70.9%</td></tr></tbody></table>

#### 3.1.1. Impact of GNSS Data Processing

In comparison to CODE, IGS is noisier with respect to the reanalysis with about 10% higher monthly variances in the time-matched data on average (Table 2). Inspection of the monthly variances by station shows that the noise is systematically larger in the IGS data set for all stations except at USUD (Figure 3c). All the stations located in Antarctica (DAV1, MAW1, SYOG, MCM4, MAC1, CAS1) show a significantly larger relative monthly variance in the IGS data set, with a maximum of 43% at DAV1. We believe that the use of the VMF1 mapping function in the CODE solution makes a major difference as it accounts for the high-resolution (6-hourly) temporal variations in the atmospheric (dry and wet) layering, whereas the GMF used in the IGS solution only models a smooth seasonal variation. On the other hand, the stdf, which quantifies the magnitude of the periodic bias, does not show such a systematic difference, although CODE shows smaller values again on average. A majority of stations in the IGS data set have slightly larger stdf (54 stations, i.e., 67%,
larger versus 27 stations, i.e., 33%, smaller). For a few stations, the difference in stdf is relatively big, exceeding ±50% (Figure 3d,e). The difference in stdf can be explained by the difference in mapping functions, a priori ZHD, and elevation cutoff angle. Specific cases are discussed below.

<figure><img src="image_18.png" /><figcaption>**Figure 3.** Comparison of segmentation results from two different GNSS data sets, IGS repro1 (IGS), and CODE REPRO2015 (CODE), at 81 common stations: (a) number of detected change-points, (b) mean of the monthly variance, (c) relative difference in the mean of monthly variance (IGS – CODE), (d) standard deviation of the functional modeling the periodic bias (stdf), (e) relative difference in the standard deviation of the functional (IGS – CODE). All the parameters are estimated from time-matched series for the period 1995–2010. The relative differences are computed relative to the mean (IGS + CODE)/2.</figcaption></figure>

The segmentation found more change-points in the CODE data set, both before and after the screening, compared to IGS, and slightly more outliers (38 versus 36). After screening, the validation percentages of the two data sets are very similar (28.9 and 29.9%).
However, the number of similar change-points is relatively small (48.8%), indicating a strong dependence of the segmentation results on the GNSS processing method. There are several reasons that can explain the difference in the number, as well as the position, of change-points. Firstly, the segmentation is sensitive to the magnitude and stationarity of the noise. Since the noise is higher for the IGS data set, a number of small offsets are not detected there. Secondly, IGS includes mapping function changes in 2008 and 2009, which are not included in CODE. They lead to 15 extra-detections in IGS. Thirdly, the improved tropospheric model and updated antenna/radome calibrations used in the CODE processing are likely to reduce some mean biases and seasonal variations in addition to the reduced noise, which may lead to a difference in the number and positions of change-points (e.g., higher periodic biases in IGS may lead to extra change-points). To summarize, some of the data features may increase the number of change-points in IGS, while others may increase the number in CODE. Below, we examine a few special cases.

Figure 4 shows a pathological example (station MCMC4, McMurdo, Antarctica). For this station and all other stations in Antarctica, the difference in a priori ZHD leads to a reduced bias (mean and seasonal) in the CODE IWV estimates with respect to ERAI. The mean difference between CODE and IGS for MCM4 is about −0.5 kg m⁻². Large oscillations are seen in both solutions during the period 2002–2006. These oscillations are believed to arise in the GNSS measurements due to snow intrusion between the antenna and radome as discussed by Koulali and Clarke [49]. Most of the stations in Antarctica show this feature for some periods. In the CODE solution, the oscillations are slightly smaller, which may be explained by the use of a lower cutoff angle which would mitigate the anomalous path delay variations partly. In the IGS series, these oscillations lead to several extra change-points, although they are not due to equipment changes and, thus, are not validated by the metadata. The mapping function change in 2008 (from GMF to NMF) and 2009 (from NMF to GMF) in IGS also leads to 2 additional change-points. As a consequence, the IGS series at MCM4 has 17 change-points, including 6 outliers due to high noise, whereas the CODE solution has only 5 change-points. The number of validations with respect to metadata are 5 and 2, for IGS and CODE, respectively, with one validated change-point in 2006 (due to a receiver change) similar in both cases.

Next, we examine another example (station GOPE, Czech Republic), where the equip-
ment changes have a strong and different impact on the CODE and IGS solutions. Over the
study period (1995–2010), many changes occurred at this station. Figure 5 gives a visual
description of these changes based on the IGS metadata, where each color represents a
specific equipment type (producer and product), and the vertical dashed lines indicate a
sub-type change (serial number or firmware for receivers). Five different antenna/radome
types were used at this site: TRM14532.00/NONE (from 13 May 1995 to 4 November 1999
and from 24 July 2000 to 4 October 2000), ASH701073.3/SNOW (from 4 November 1999 to
24 July 2000), ASH701946.3/NONE (from 4 October 2000 to 14 July 2006), and
TPSCR3_GGD/CONE (from 14 July 2006 to 14 December 2009), TPSCR.G3/TPSH (since
14 December 2009). Three different receiver types were used: TRIMBLE_4000SSE (from
13 May 1995 to 4 November 1999, and 24 July 2000 to 4 October 2000), ASHTECH_Z18
(from 4 November 1999 to 24 July 2000, and 4 October 2000 to 14 December 2009), and
TPS_NETG3 (from 14 December 2009 to 8 August 2018). The figure shows that the same
receivers were removed and reinstalled for some periods, as well as that they got several
firmware updates (see the dashed vertical lines). The cutoff angle was also changed twice,
going from 15 degrees to 5 degrees on 4 November 1999, and from 5 degrees down to
0 degree on 14 December 2009. The impact of these equipment and setting changes on the
GPS-ERAI IWV differences is obvious from Figure 6. For example, there is a strong change
in the mean IWV difference for both solutions before and after 4 October 2000. This change
is detected (and validated) by the segmentation in both data sets. Before that date, CODE
has 2 other change-points, while IGS has only one, although IGS shows similar changes in
the mean to CODE but they are not detected. Neither of these change-points is explained by
the metadata. After 4 October 2000, CODE has 3 more change-points, while IGS has only 2.
All 3 are validated in the case of CODE, and none is validated in the case of IGS. However, the IGS detected change-points seem correct as they well represent the changes in the mean of that series. The difference in the number and positions of change-points between the IGS and CODE segmentation results is actually due to the difference in the two GNSS IWV data, as highlighted in the lower plot of Figure 6. This difference can be explained by the different antenna/radome calibration models used in the two GNSS processing. The most striking bias between the two series is for the period from 14 July 2006 to 14 December 2009 where the TPSCR3_GGD/CONE antenna/radome pair is used. In the igs05.atx calibration file used in the IGS solution, this calibration was of “Field” type (i.e., an older relative calibration dating back to April 2005), while the CODE solution used a more recent “Robot” type calibration (from April 2013). For all other antenna/radome types used at this station, both GNSS solutions used similar calibration types (“Robot” or “Copied”), although the CODE calibrations were more recent in all cases, which may explain the small residual biases between the two solutions. Another abrupt shift in the CODE–IGS series is observed on 1 January 2002 which coincides with the introduction of GLONASS observations in the CODE solution. A closer look reveals that, for some reason, the shift is merely in the IGS–ERAI series and that the CODE–ERAI series is, rather, drifting.

The GOPE example shows that antenna/radome changes at a site are likely to produce an abrupt shift in the mean IWV, even though calibration models are used during the data processing. Different calibration types and versions have different biases, leading to different inhomogeneities in the two reprocessed GNSS data sets analyzed in this work and to different segmentation results. Many such cases could be found by inspecting the time series and segmentation results. In some cases, receiver changes and observation cutoff angle changes also induce abrupt shifts of different amplitudes for the two data sets. This feature might be due to the use of different processing cutoff angles (3 degrees in CODE and 7 degrees in IGS). In other cases, we can see a small drift (below 0.5 kg m⁻²) between the two solutions, which might be due to the change in the number of satellites used in the processing, especially with the inclusion of GLONASS observations in the CODE solution starting from 2002, although not many GNSS stations were recording GLONASS data in the period from 2002 to 2010. Since the amplitude of such small drifts is generally small compared to the noise, it would not induce extra change-points in the segmentation; although this cannot be guaranteed for sure, it may impact the trend estimates.

### 3.1.2. Impact of the Length of Time Series

The impact of the length of time series is inspected from the comparison of CODE-ERAI segmentation results where the segmentation was run both on the full (1994–2018) series and on the series limited to 1994–2010. The change-points statistics (number of detections, outliers, and validations) reported in the second section of Table 2 are given for the common period (1994–2010).

It can be deduced from Table 2 that the estimated monthly mean variance and the estimated stdf are nearly similar on average over all stations. However, inspecting individual stations shows that small differences exist (Figure 7b–e). Slightly more than half of the stations have smaller variance in the full series (57 versus 43%) and smaller stdf (52 versus 48%). This difference can be explained by the fact that the GNSS IWV estimates become more accurate with time; thus, the agreement with the reanalysis improves in the more recent years included in the full series. On the other hand, the cases where the noise is larger in the full time series may either be due to situations where the GNSS measurements get noisier, or the ERA-Interim reanalysis becomes less accurate with time. This is, for example, the case at station COCO where the monthly variance is 25% higher in the full series. The same situation is observed at many tropical stations, although the quality of the GNSS data has generally improved at these sites. We suspect that the problem is in the quality of the ERAI reanalysis, which may indeed have degraded in recent years due to the end of several satellite missions [18]. One obvious case is with station COCO discussed in the next sub-section.
<figure><img src="image_19.png" /><figcaption>**Figure 4.** Time series of IWV differences at station MCM4 (McMurdo, Ross Island, Antarctica). The two top plots show the results for IGS-ERAI and CODE-ERAI, with segmentation results superposed (the red line shows the weighted mean as estimated by the segmentation method, where the black dashed lines indicate the change-point positions, the cyan line is the estimated monthly mean variance, and the purple line is the estimated functional modeling the periodic bias). The validated change-points are indicated by a marker at the lower end of the black dashed line. The 3rd plot from the top shows the two IWV difference series and detected change-points for clarity of the comparison. The lower plot shows the difference, IGS-CODE, with change-points superposed. In all plots, the markers at the top indicate the equipment changes. The color code is given in the lower plot. The processing changes apply only to IGS.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 5.** Schematic view of the equipment changes that occurred at station GOPE (Ondrejov, Czech Republic) during the period 1995–2010. Each color symbolizes a different equipment type (producer and product). Vertical dashed lines show additional changes with the same type (serial number or firmware upgrade in the case of receivers). Information extracted from the IGS site logs (ftp://igs.ign.fr/pub/igs/igscb/station/log/, accessed on 30 July 2021).</figcaption></figure>

The number of change-points detected from the long time series in the common period is smaller than for the short time series, both before and after screening (Table 2). The same is found for the number of outliers, while the percentage of validations is consistently increased. These results indicate that, on average, the segmentation yields more accurate results with the long time series. The difference in the number of change-points can be understood from the fact that the penalty criterion is conservative and tends to choose a small number of change-points. For instance, the total number of change-points found from the long time series is 321 (not shown), i.e., not much larger than the total number for the short time series (296, as reported in Table 2). Because of this conservative property, the number and positions of the change-points in the two series cannot be expected to be the same at most stations. Figure 7 shows that 39 stations (48%) have the same number of change-points (although the positions may not be the same), 16 stations (20%) have a larger number in the long series, and 26 stations (32%) have a smaller number in the long series. Nevertheless, 80.8% of the change-points are similar (i.e., within ± 62 days).

Figure 8 shows the example of station VILL (Villafranca, Spain), where the mean noise and stdf are very similar for the two-time series, but the segmentation results are quite different: 6 change-points are detected in the short time series and 2 in the long time series. Two change-points are similar in both series and are validated. The additional change-points found in the short time series capture short but significant variations in the mean. These change-points are not retained in the long time series because other such variations are seen all along in the long time series. The penalty criterion avoids selecting all those change-points, and the final optimal solution eventually has only one change-point. This solution seems more reasonable.

### 3.1.3. Impact of the Reference Data Set

The third section of Table 2 summarizes the segmentation results when either ERAI or ERA5 is used as a reference, i.e., the segmentation is run for CODE-ERAI or CODE-ERA5 IWV differences, when the same auxiliary data is used in the GNSS ZTD to IWV conversion (here from ERA5). Globally, both the mean noise and the stdf with ERA5 are reduced by 25% compared to ERAI. Figure 9c,e show that the reduction in noise and stdf is observed at 95% and 75% of all stations, respectively. This difference can be explained by the lower representativeness error in ERA5 due to higher spatial resolution, as well as higher quality of the IWV temporal variations in this reanalysis, probably due to the assimilation of more satellite observations. Figure 10 shows the most striking case of increase in noise with time for ERAI. The impact on the segmentation results is quite significant. Only one change-point is detected in the more noisy series, while seven change-points are detected in the less noisy one. At many stations, there is also an excess noise in ERAI during the moist period of each year; see the example of station KIRU (Kiruna, Sweden) in Figure 11. At this station, the CODE-ERAI series has much larger seasonal variations in the noise and in the functional than CODE-ERA5. This leads to more change-points in more noisy series (6 versus 2) because of the sharp increase in the noise during some years, which is not well represented by the periodic bias and monthly variance. As a result, four outliers are detected in the CODE-ERAI segmentation results. In the CODE-ERA5, the two change-
points are validated by the metadata, which is, again, better than in the CODE-ERAI results. Other examples with large changes in stdf are stations KIT3, POL2, and SANT (Figure 9d). These stations were also pointed as extreme cases of representativeness errors in ERAI by Bock and Parracho [17] due to steep orography.

<figure><img src="image_21.png" /><figcaption>**Figure 6.** Similar to Figure 4, but for station GOPE (Ondrejov, Czech Republic).</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 7.** Similar to Figure 3, but comparing segmentation results from two different lengths of the CODE data set (long time series, from 1994 to 2018, short time series, from 1994 to 2010). The mean noise and standard deviation of the functional (stdf) are representative of the full times series, but the number of breakpoints are counted in the common period only (1994–2010).</figcaption></figure>

From Table 2, we see that the total number of change-points is larger when ERA5 is used as a reference. This can be understood as the consequence of the general decrease of the noise and periodic bias with this reanalysis, as discussed in the case of station COCO above. When the noise is decreased, it is easier to detect a small offset in the time series. With the increase in the number of change-points, the number of outliers is increased, as well. However, after screening, the percentage of validations is higher with ERA5 as a reference. So, there is a clear benefit of using the more recent reanalysis as a reference.
<figure><img src="image_23.png" /><figcaption>**Figure 8.** Similar to Figure 4, but for station VILL (Villafranca, Spain).</figcaption></figure>

Figure 9a shows that 37 stations (46%) in CODE-ERA5 have a higher number of change-points than CODE-ERAI, 29 stations (36%) have a smaller number, and 15 stations have a similar number. From Figure 9c,e, we see that MKEA (Mauna Kea volcano, Mauna Kea, USA) and USUD (Usuda, Japan) are two cases where the mean noise or stdf increased with ERA5 as a reference by 40% and 107%, respectively. Both stations are located in regions of steep topography where both reanalyses have significant representativeness errors compared to the GNSS observations. In the case of MKEA, the station is located at an altitude of 3729 m, whereas the altitudes of the surrounding grid points from both reanalyses are much lower. In the case of USUD, the situation is opposite, with the station being closer to the sea level than the surrounding grid points from the reanalyses.

### 3.1.4. Impact of the Auxiliary Data Set

The auxiliary data used in the conversion of GNSS ZTD to IWV impacts the quality of the GNSS IWV data and may lead to different segmentation results in a similar way as the processing and reference data sets. Table 2 shows that on average the mean noise and stdf are the same, but the segmentation statistics are slightly different (number of change-points, outliers, and validations). Figure 12 shows that the noise and stdf results actually change for many stations. In general, the absolute values of the noise are very close, but the relative differences are not that small. At 60% of the stations, ERAI induces larger noise than ERA5, with values up to 10–20%, while, at 40% of the stations, ERA5 yields similar or higher mean noise in ERAI, but the relative increase there is small (2.5% at maximum). These results are consistent with the representativeness differences between the reanalyses discussed above, although the pressure and temperature data are much less subject to small-scale variations than IWV.
<figure><img src="image_24.png" /><figcaption>**Figure 9.** Similar to Figure 3 but comparing the segmentation results using two different reference data sets, ERA-Interim (ERAI) and ERA5.</figcaption></figure>

The results are similar for stdf (60% of the stations have a larger periodic bias with ERAI), but the relative difference can be much higher (up to ±80%). This is because many stations are located in complex regions, such as the mountains and near the oceans. In some cases, ERA5 induces a larger periodic bias compared to the ERAI, for instance, at CHUR (Churchill, MB, Canada), KERG (Port aux Français, French Southern Territories), and TABL (Wrightwood, CA, USA).
<figure><img src="image_25.png" /><figcaption>**Figure 10.** Similar to Figure 4, but for station COCO (Cocos (Keeling) Island, Australia).</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure 11.** Similar to Figure 4, but for station KIRU (Kiruna, Sweden).</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Figure 12.** Similar to Figure 3, but comparing segmentation results from GNSS data sets that used two different auxiliary data, from ERA-Interim (ERAI) and ERA5.</figcaption></figure>

Although the total number of change-points in the two data sets are very similar before and after screening (see Table 2), the number of change-points can be quite different from station to station (Figure 12a). The results with auxiliary data from ERAI show 18% more outliers than with ERA5. From this perspective, it is better to use ERA5. However, the percentage of similar change-points is still quite high (around 71%), which points to a moderate impact of the auxiliary data on the segmentation results in the end.
## 3.2. IWV Trend Estimates

### 3.2.1. Impact of GNSS and Reanalysis Data Set Properties on Trend Estimates

Table 3 summarizes the trend results obtained with the different GNSS data sets and the two reanalyses discussed in Section 3.1. The numbers report the mean and standard deviation of the trend estimates (in kg m⁻² year⁻¹) over the 81 stations, as well as the number of significant trends at the 0.05 level (using a Student's t-test), and the standard error in the trend estimate (1 − σ). Following from Section 3.1, three-time periods, with lengths 16, 17, and 25 years, are presented, respectively.

<table><caption>**Table 3.** Summary of IWV trends from various data sets used in this work. The number of stations with significant trends at level $\alpha = 0.05$ is given in brackets. (a) GNSS data converted with auxiliary data from ERAI and segmentation applied the CODE—ERA5 IWV difference. (b) GNSS data converted with auxiliary data from ERA5 and segmentation applied the CODE—ERAI IWV difference. (c) GNSS data converted with auxiliary data from ERA5 and segmentation applied the CODE—ERA5 IWV difference.</caption><thead><tr><th>Time Span</th><th>1995–2010</th><th>1994–2010</th><th>1994–2018</th></tr></thead><tbody><tr><td>Std error<br>(kg m<sup>-2</sup> year<sup>-1</sup>)</td><td>0.035</td><td>0.033</td><td>0.018</td></tr><tr><td>ERAI<br>(kg m<sup>-2</sup> year<sup>-1</sup>)</td><td>0.018 ± 0.055<br>(9)</td><td>0.013 ± 0.049<br>(10)</td><td>0.027 ± 0.034<br>(37)</td></tr><tr><td>ERA5<br>(kg m<sup>-2</sup> year<sup>-1</sup>)</td><td>0.011 ± 0.052<br>(8)</td><td>0.008 ± 0.047<br>(8)</td><td>0.027 ± 0.031<br>(35)</td></tr><tr><th>GPS</th><th>IGS<br>time-matched</th><th>CODE<br>time-matched</th><th>CODE<br>time-limited</th><th>CODE (a)</th><th>CODE (b)</th><th>CODE (c)</th></tr><tr><td rowspan="2">Raw data</td><td>IWV trend<br>(kg m<sup>-2</sup> year<sup>-1</sup>)</td><td>0.024 ± 0.059<br>(20)</td><td>0.018 ± 0.060<br>(18)</td><td>0.016 ± 0.060<br>(23)</td><td>0.033 ± 0.032<br>(46)</td><td>0.030 ± 0.031<br>(41)</td><td>0.030 ± 0.031<br>(41)</td></tr><tr><td>RMSE wrt ERA5<br>(kg m<sup>-2</sup> year<sup>-1</sup>)</td><td>0.044</td><td>0.046</td><td>0.046</td><td>0.033</td><td>0.033</td><td>0.033</td></tr><tr><td rowspan="2">corrected IWV<br>by validations</td><td>IWV trend<br>(kg m<sup>-2</sup> year<sup>-1</sup>)</td><td>0.015 ± 0.052<br>(12)</td><td>0.014 ± 0.052<br>(11)</td><td>0.011 ± 0.052<br>(15)</td><td>0.027 ± 0.027<br>(34)</td><td>0.025 ± 0.030<br>(34)</td><td>0.027 ± 0.026<br>(34)</td></tr><tr><td>RMSE wrt ERA5<br>(kg m<sup>-2</sup> year<sup>-1</sup>)</td><td>0.038</td><td>0.039</td><td>0.040</td><td>0.019</td><td>0.022</td><td>0.019</td></tr><tr><td rowspan="2">corrected IWV by<br>all breakpoints</td><td>IWV trend<br>(kg m<sup>-2</sup> year<sup>-1</sup>)</td><td>0.017 ± 0.053<br>(9)</td><td>0.016 ± 0.054<br>(9)</td><td>0.012 ± 0.048<br>(13)</td><td>0.027 ± 0.030<br>(33)</td><td>0.027 ± 0.032<br>(35)</td><td>0.027 ± 0.030<br>(34)</td></tr><tr><td>RMSE wrt ERA5<br>(kg m<sup>-2</sup> year<sup>-1</sup>)</td><td>0.021</td><td>0.022</td><td>0.022</td><td>0.006</td><td>0.012</td><td>0.006</td></tr></tbody></table>

From the two reanalyses, we see that the mean trends are positive, indicating a net moistening, globally, with slightly different values between the three periods. This reminds us that the mean linear trends from different periods may not generally agree because they are strongly influenced by inter-annual to inter-decadal variability. However, the decrease in the standard deviation is noticeable from the shorter to the longer period (e.g., from 0.052 kg m⁻² year⁻¹ to 0.031 kg m⁻² year⁻¹ for the ERA5 data set), which indicates a decreasing influence of the inter-annual variability with time, as well as more consistent trend estimates from the global network with long time series. This decrease is also seen in the GNSS data sets, raw and corrected. It is also consistent with a decrease in the standard error with the longer time series, from 0.035 to 0.018 kg m⁻² year⁻¹, and the subsequent increase in the number of significant trends, e.g., from 8 to 35 with ERA5. ERAI and ERA5 show different means and standard deviations in the short periods, which is not surprising according to the differences in the IWV time series extracted from the two reanalyses (see Section 3.1). The mean difference is negligible on the longer period, but the variability is slightly smaller in ERA5. The RMS difference amounts to 0.013 kg m⁻² year⁻¹ (not reported in Table 3), which indicates that there are substantial local differences in the trends from the two reanalyses. Note that the mean positive (moistening) IWV trend of 0.027 kg m⁻² year⁻¹ from the reanalyses (and also from the GNSS data) is fairly consistent with the prediction from Clausius–Clapeyron law of 7% IWV increase per 1 °C induced by the global increase in temperature of ~1 °C over the past four decades, given a global mean IWV of 18 kg m⁻² [50].
Next, we examine the results for the two GNSS data sets, IGS and CODE, before and after homogenization, and their differences with respect to ERA5. In this section, we consider two different corrected (homogenized) data sets. In the first one, we use only the change-points validated from the metadata, while, in the second one, we use all the detected change-points. We will refer to these data sets as “partially corrected” and “fully corrected” data sets, respectively. Ideally, we could also consider a third version where only the change-points attributed to GNSS are included, but this is not possible here because no nearby stations are available in many cases (see the discussion in Section 2). The raw GNSS trends show quite a large difference in the mean (0.024 versus 0.018 kg m⁻² year⁻¹) and a RMS difference of 0.016 kg m⁻² year⁻¹ (not reported Table 3). This difference is not unexpected given the differences in the data processing (Section 2) and the inhomogeneities they induce, as discussed from the segmentation results in Section 3.1. Especially, the inhomogeneities in the IGS data set due to the older antenna/radome calibration models may be a significant cause of uncertainty in the trends. The mean difference is reduced in both corrected data sets, although the RMS difference is not reduced in the partially corrected data set (0.019 kg m⁻² year⁻¹), contrary to what one would expect after both data sets are homogenized. This result can be understood from the fact that the segmentation results of the IGS and CODE data sets are sometimes very different and the validated change-points may not coincide in both solutions (see Figures 4 and 6). On the other hand, the IGS and CODE GNSS fully corrected data sets are much more consistent (RMS difference of 0.006 kg m⁻² year⁻¹). The latter result gives good confidence that the segmentation method is able to detect all the significant change-points in either data set. However, we know that not all these change-points may come from the GNSS time series, but some of them may be due to inhomogeneities in the reference reanalysis (in this case, ERAI). As a consequence, the fully corrected GNSS trends will be very close to the trends from the segmentation reference data set in the end. In Table 3, we give the RMS difference between the GNSS data sets and ERA5 (which is taken as another reference, although not independent from ERAI). In general, the RMS differences between the fully corrected GNSS trends and ERA5 are significantly smaller than between the raw or the partially corrected trends and the ERA5 trends. We note also that the number of significant trends is drastically reduced for both corrected GNSS data sets (from 20 to 9 for IGS and from 18 to 9 for CODE). This indicates that a large portion of uncorrected stations had significant trends because of inhomogeneities in their time series.

The results from the CODE time-limited data set are quite similar, although the change in the mean trend is reduced between the uncorrected and the corrected series, and the agreement with ERA5, in the end, is improved. The number of significant trends is also higher than with the time-matched data sets. Recall that the main difference between this data set and the CODE time-matched is only one more year and fewer gaps, but these differences can have a sensible impact on the segmentation results and trend estimates.

The results from the long time series using either ERAI or ERA5 as auxiliary data or reference data for the segmentation are presented in the rightmost part of Table 3. The mean trends from the uncorrected GNSS data are slightly larger than the trends from the reanalyses (0.030–0.033 kg m⁻² year⁻¹ for GNSS compared to 0.027 kg m⁻² year⁻¹ for ERAI and ERA5). The corrected GNSS series achieve closer mean trends to the reanalyses, as well as reduced RMS differences. The smallest RMS difference with respect to ERA5 is found with the fully corrected GNSS data using ERA5 as a reference, which is to be expected (ERA5 is not independent in this case). In terms of variability, the partially corrected trends show slightly smaller standard deviation, i.e., smaller spatial variability, which suggests more homogeneous and consistent trends in the global network. The slightly higher variability in the fully corrected GNSS trends might be due to some inhomogeneities in the reference series (ERA5 or ERAI). This point is further discussed in the next subsection. The impact of the segmentation reference (ERAI versus ERA5) on the corrected GNSS trends is small in terms of mean, but the RMS difference between the GNSS trends at individual stations amounts to 0.015 kg m⁻² year⁻¹ for the partially corrected series
and 0.012 kg m⁻² year⁻¹ for the fully corrected series (not shown). The impact of the auxiliary data is significantly smaller, with a RMS difference between the trends using ERAI and ERA5 as auxiliary of 0.008 kg m⁻² year⁻¹ for the partially corrected series and 0.002 kg m⁻² year⁻¹ for the fully corrected series (not shown).

### 3.2.2. Impact of Homogenization on GNSS Trend Estimates

In this sub-section, we analyze in more details our 'best' data set at hand, i.e., the long CODE GNSS series (1994–2018) using ERA5 as auxiliary data and reference for the segmentation. Figure 13a shows the IWV trend estimates for the GNSS data, uncorrected and corrected, and ERA5. A majority of the GNSS stations have positive trends (89% with the partially corrected data, 86% with the fully corrected data, versus 80% for the ERA5 data), consistent with the overall positive mean trend of 0.027 kg m⁻² year⁻¹ reported in Table 3. Among these, only ~41% of the corrected GNSS trends are significant at *p* = 0.05 (or *t* = 1.99, see Figure 13b), versus 49% with the uncorrected GNSS data, and 41% with ERA5. The largest positive trend is found for KOUR (Kourou, French Guiana), reaching a value of ~0.110 kg m⁻² year⁻¹ (*t* = 5) for the uncorrected and partially corrected trends, and 0.150–0.153 kg m⁻² year⁻¹ (*t* ~ 7) for the fully corrected GNSS data and ERA5. This station, as well as a few others in the Tropics (KOKB, GUAM, BRMU, CRO1, etc.), shows consistent and significant moistening over the past 2.5 decades. Strong moistening is also found at several Mediterranean stations (MEDI, CAGL, MATE), confirming the strong warming in this area [51]. A few stations show consistent drying from GNSS data (uncorrected and corrected) and ERA5, mostly in arid regions, such as JPLM (Pasadena, CA, USA), HRAO (Krugersdorp, South Africa), ALIC (Alice Springs, Australia), and SANT (Santiago, Chile).

A striking feature in Figure 13 is that many uncorrected GNSS trends are quite large, and significant, while the corrected trends are much smaller, and are, in some cases, insignificant. The trends of the partially corrected series are significantly different from the uncorrected trend (with a RMS difference of 0.022 kg m⁻² year⁻¹), although they agree with each other better than with ERA5 (0.033 kg m⁻² year⁻¹). On the other hand, the trends of the fully corrected GNSS data are much closer to ERA5 (with a RMS difference of 0.006 kg m⁻² year⁻¹) than to the uncorrected trends (RMS difference of 0.031 kg m⁻² year⁻¹). The RMS difference between the partially corrected and the fully corrected trends is 0.016 kg m⁻² year⁻¹. It is actually expected that the corrected trends will tend to align with the reference data (in this case, ERA5). However, the partially corrected trends are relatively independent from the reference since only about 36.5% of the detected change-points are validated. It is also noticeable that the partial correction has a strong impact on the trends at many sites, such as a change in the sign (e.g., GOPE, KARR, VILL, etc.), a change from significant to insignificant (e.g., HERS, SHAO, SYOG, WUHN, etc.), or vice versa (GOPE, KOKB, TIDB, etc.). The trends of the fully corrected data get closer to ERA5.

Four cases will be discussed in more detail in the following, where the IWV trends change from significant to insignificant, or vice versa. The corresponding time series and change-points are shown in Figure 14.

Firstly, we examine the case of station WUHN (Wuhan, China), where the uncorrected GNSS trend is positive and significant (0.115 kg m⁻² year⁻¹, *t* = 2.44), while the fully corrected trend and the ERA5 trend are of opposite sign (−0.031 and −0.037 kg m⁻² year⁻¹, respectively) but not significant (*t* = −0.65 and −0.79). At this site, the change in sign of the trend is suspected to arrive from a strong downward shift in the ERA5 time series, which was already commented by Parracho et al. [14] and attributed to a change in the radiosonde type in Wuhan in 2006 impacting the ERA-Interim reanalysis. Parracho et al. [14] also noticed an extended region of a negative trend in eastern China, in contradiction with the GNSS observations and the MERRA-2 reanalysis. A similar shift is actually detected in the nearby station SHAO (Sheshan, China), where the uncorrected GNSS trend of 0.086 kg m⁻² year⁻¹ is decreased to 0.043 kg m⁻² year⁻¹ in the fully corrected data. The
ERA5 trend is even smaller (0.037 kg m⁻² year⁻¹), although not negative. At both WUHN and SHAO, the GNSS trend changes from significant to insignificant after correction.

<figure><img src="image_28.png" /><figcaption>**Figure 13.** Comparison of: (a) trend estimates kg m⁻² year⁻¹, and (b) t-values from different uncorrected (raw) and corrected GNSS data (CODE with auxiliary ERA5 and reference ERA5) and ERA5 at 81 common stations.</figcaption></figure>

Secondly, let us examine the case of station HERS (Hailsham, UK). The segmentation detected six change-points, among which only two are validated with the metadata. Overall, the mean shifts are going upwards, meaning that the inhomogeneities induce a spurious positive trend in the GNSS series. The correction of the GNSS series for two validated change-points has a strong impact on the trend, decreasing it from 0.081 to 0.024 kg m⁻² year⁻¹ and from significant to insignificant (t-value from 5.1 to 1.5). Including the 4 additional change-points has a further, although small, effect, leading to a final GNSS trend of 0.030 kg m⁻² year⁻¹, close to the ERA5 trend (0.031 kg m⁻² year⁻¹). Two nearby stations (HERT and HRM1) could be used in the attribution step to confirm the two validated change-points but not the other ones. The other change-points could not be tested. The impact of the correction is substantial and seems justified at this station, with a final trend reduced by -0.051 kg m⁻² year⁻¹, i.e., a factor of 2.7.
<figure><img src="image_29.png" /><figcaption>**Figure 14.** Similar to the upper plot in Figure 4, but for different stations: HERS (Hailsham, United Kingdom), GOPE (Ondrejov, Czech Republic), KOKB (Waimea, United States), and GUAM (Dededo, Guam). The IWV differences are computed as GNSS - ERA5, where GNSS is converted using auxiliary data from ERA5, and the segmentation is run with ERA5 as a reference.</figcaption></figure>

Next, we examine station GOPE (Ondrejov, Czech Republic), which has a strong significant trend after correction but insignificant before. GOPE is a special case, which has a negative trend in the raw data in contradiction was many surrounding stations in Europe, such as ZIMM (Switzerland), WTZR (Germany), and BOR1 (Poland), which have positive trends. This feature was already noticed by Parracho et al. [14] from the uncorrected IGS data set over the shorter period (1995–2010). Figure 14 shows that the mean shifts are going
downwards, so inducing a negative trend in the GNSS series compared to ERA5. Two change-points are validated with the metadata. After correction of these change-points, the trend goes from a insignificant drying of −0.020 kg m⁻² year⁻¹ to a significant moistening of 0.046 kg m⁻² year⁻¹. Three other change-points have a minor impact (the fully corrected trend is 0.044 kg m⁻² year⁻¹) because the most important break in 2000 is validated. For this station, we could also test the attribution with several nearby stations collocated with station WTZR (distant by 162 km). The two validated change-points, as well as the one in 2001, could be attributed to GOPE.

The last example is station GUAM (Dededo, Guam), in the western tropical Pacific, which has a similar large trend in ERA5 to KOKB, another station in the Pacific Ocean. The trends are very different between the partially and fully corrected GNSS series at GUAM because only one change-point is validated, and it is located near the beginning of the series. The last three change-points have a strong impact on the GNSS correction, although their origin is questionable. Indeed, they are located quite far away from any known equipment change reported in the metadata. The last change-point (on 26 September 2017) could be checked in the attribution step with the nearby station GUUG (Mangilao, USA), located at a distance of 18 km from GUAM. Comparing the GNSS series at GUUG to the ERA5 series at GUAM revealed a significant change in mean on this date. From this result, we should attribute this change-point to the ERA5 series and not the GNSS series. At this site, thus, we also suspect the other unvalidated change-points to be due to ERA5. This assumption may be further checked by inspecting observation statistics from the assimilation system, but this is left for future work.

## 4. Conclusions

This study investigated the sensitivity of the segmentation method and the IWV trend estimates to different GNSS and reference data properties. It was shown that the GNSS processing methods and the reference data (ERAI versus ERA5) have the strongest impact on the segmentation results (i.e., number and positions of change-points), while the impact is weaker when the length of the time series (17 versus 25 years) or the auxiliary data (ERAI versus ERA5) are changed. Changing the latter two was shown to achieve 81% and 71% similar change-points, but only 45–49% when the GNSS data set or the reference was changed. These discrepancies in the results indicate that the segmentation is sensitive to small changes in the mean signal, in the noise, and in the periodic bias. These features are determined by several aspects of the GNSS processing procedure, especially: the a priori ZHD correction, the antenna/radome calibration model, the mapping functions, and the elevation cutoff angle. The more recent and more accurate ZHD correction, antenna/radome calibration model, and mapping functions used in the CODE reprocessing achieve smaller noise and periodic bias but, at the same time, lead to an increase in the number of detected change-points because the segmentation is able to detect smaller changes in the mean. The same is observed when the reference is changed from ERAI to ERA5. The reduction in the noise and periodic bias in ERA5 is partly due to the better spatial resolution and partly to the more recent atmospheric model and assimilation of more observations. In terms of the percentage of validation of the detected change-points with respect to the GNSS metadata, the length of the time series has the strongest impact. The average percentage is ~30% with the short time series (16 or 17 years), while it rises to 34–36% with the long time series (25 years), for a window of ±62 days. This difference might be due to the particular penalty function that we used here [46,47]. The penalty function might be calibrated differently to achieve a better balance between the probability of detection and the probability of false detection to our particular data [6].

This study also points to the strong impact of the inhomogeneity correction on the estimated trends. When only the change-points validated with the metadata are used, the impact of the correction is substantial on the global mean trend, on the dispersion, and on the RMS difference with respect to ERA5. When all the detected change-points are used, the mean and dispersion are again modified, but to a lesser extent, and the RMS
difference with ERA5 is further reduced. The latter feature is expected when ERA5 or ERAI are used as a reference in the segmentation, especially since a few change-points may be due to the reanalyses (e.g., suspected at stations WUHN, GUAM). The uncertainty in the estimated trends from the use of the different reprocessed GNSS data sets or reference reanalysis in the segmentation is about 0.015–0.019 kg m⁻² year⁻¹ (RMS difference between the tested pairs of data sets) when only the validated change-points are used and 0.002 to 0.012 kg m⁻² year⁻¹ when all the change-points are used. The auxiliary data has a marginal impact on the trends. The longer time series (25 years) provides higher accuracy in the trend estimates, with a mean standard error of 0.018 kg m⁻² year⁻¹ and a dispersion of ~0.03 kg m⁻² year⁻¹ throughout the global network. The homogenized GNSS trends and both reanalyses agree in the global mean IWV trend estimate of 0.027 kg m⁻² year⁻¹, which indicates a global moistening over the past 2.5 decades at a rate close to the 7% per degree of warming predicted by Clausius Clapeyron equation [50].

The main limitation of our current homogenization method is the attribution step (Figure 2), which could not be used to check all the change-points because of the lack of nearby stations. In future work, we plan to include reprocessed GNSS data from additional stations in regions where dense networks are available (e.g., Europe, USA, Japan, etc.). Some limitations of the current segmentation method were highlighted, as well, especially the dependence on the length of the time series, which may be mitigated by careful calibration of the penalty function. The segmentation method is also sensitive to the long-term variations in the noise and bias, which are not specifically modeled here (we assume both are periodic with a fundamental period of 1 year). At some sites, strong inter-annual or decadal variations in the noise and/or bias were captured by the segmentation. These spurious change-points need to be detected and removed, unless the long-term variations are reduced, e.g., from the use of a better reference data set, such as a nearby GNSS station, rather than a reanalysis.

Finally, this study also helped to better understand the performance and limitations of the GNSS data processing procedures. This knowledge can be helpful in the future to reprocess the GNSS data and produce better “homogenized” daily IWV time series at the processing level.

**Author Contributions:** E.L., O.B. and K.N.N. developed the statistical conceptualization methods of the work. Data and climatological aspects were dealt with by O.B., and computational issues were handled by O.B. and K.N.N., all shared in the writing of the manuscript. A.Q. developed the segmentation method and R package used in this study. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the CNRS program LEFE/INSU through the project VEGAN.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The daily IWV data set based on the IGS repro1 ZTD data (120 stations, period 1995–2010) is available from https://doi.org/10.14768/06337394-73a9-407c-9997-0e380dac5591, accessed on 30 July 2021 [43]. The daily IWV data set based on the CODE REPRO2015 ZTD data (436 stations, period 1994–2018) is available from https://dx.doi.org/10.25326/18, accessed on 30 July 2021 [44].

**Acknowledgments:** The authors are grateful to AERIS, the French data and service center for atmosphere, for providing the ERA-Interim and ERA5 reanalysis data, and hosting the GNSS IWV data sets. The contribution of the fourth author has been conducted as part of the project Labex MME-DII (ANR11-LBX-0023-01) and within the FP2M federation (CNRS FR 2036).

**Conflicts of Interest:** The authors declare no conflict of interest.
## Abbreviations

The following abbreviations are used in this manuscript:

<table><tr><td>CODE</td><td>Center for Orbit Determination in Europe</td></tr><tr><td>ECMWF</td><td>European Center for Medium-Range Weather Forecasts</td></tr><tr><td>ERP</td><td>Earth Rotation Parameter</td></tr><tr><td>ERAI</td><td>ECMWF reanalysis Interim</td></tr><tr><td>ERA5</td><td>ECMWF reanalysis v5</td></tr><tr><td>GNSS</td><td>Global Navigation Satellite System</td></tr><tr><td>GPS</td><td>Global Positioning System</td></tr><tr><td>GMF</td><td>Global Mapping Function</td></tr><tr><td>IGS</td><td>International GNSS Service</td></tr><tr><td>IWV</td><td>Integrated water vapor</td></tr><tr><td>JPL/NASA</td><td>NASA's Jet Propulsion Laboratory</td></tr><tr><td>NMF</td><td>Neill Mapping Function</td></tr><tr><td>PC0</td><td>Phase Center Offset</td></tr><tr><td>PCV</td><td>Phase Center Variation</td></tr><tr><td>PWV</td><td>Precipitable Water Vapor</td></tr><tr><td>VMF</td><td>Vienna Mapping Function</td></tr><tr><td>ZHD</td><td>Zenith Hydrostatic Delay</td></tr><tr><td>ZTD</td><td>Zenith Tropospheric Delay</td></tr><tr><td>ZWD</td><td>Zenith Wet Delay</td></tr></table>

## References

1. Caussinus, H.; Mestre, O. Detection and correction of artificial shifts in climate series. *J. R. Stat. Soc. Ser. C Appl. Stat.* **2004**, *53*, 405–425. [CrossRef]

2. Reeves, J.; Chen, J.; Wang, X.L.; Lund, R.; Lu, Q.Q. A Review and Comparison of Changepoint Detection Techniques for Climate Data. *J. Appl. Meteorol. Climatol.* **2007**, *46*, 900–915. [CrossRef]

3. Jones, P.D.; Raper, S.C.B.; Bradley, R.S.; Diaz, H.F.; Kellyo, P.M.; Wigley, T.M.L. Northern Hemisphere Surface Air Temperature Variations: 1851–1984. *J. Clim. Appl. Meteorol.* **1986**, *25*, 161–179. [CrossRef]

4. Easterling, D.; Peterson, T. A new method for detecting undocumented discontinuities in climatological time series. *Int. J. Climatol.* **1995**, *15*, 369–377. [CrossRef]

5. Peterson, T.C.; Easterling, D.R.; Karl, T.R.; Groisman, P.; Nicholls, N.; Plummer, N.; Torok, S.; Auer, I.; Boehm, R.; Gullett, D.; et al. Homogeneity adjustments of in situ atmospheric climate data: A review. *Int. J. Climatol. A J. R. Meteorol. Soc.* **1998**, *18*, 1493–1517. [CrossRef]

6. Menne, M.J.; Williams, C.N. Detection of Undocumented Changepoints Using Multiple Test Statistics and Composite Reference Series. *J. Clim.* **2005**, *18*, 4271–4286. [CrossRef]

7. Szentimrey, T. Development of MASH homogenization procedure for daily data. In Proceedings of the Fifth Seminar for Homogenization and Quality Control in Climatological Databases, WCDMP-No. 71, Budapest, Hungary, 29 May–2 June 2008; pp. 123–130.

8. Costa, A.C.; Soares, A. Homogenization of Climate Data: Review and New Perspectives Using Geostatistics. *Math. Geosci.* **2009**, *41*, 291–305. [CrossRef]

9. Venema, V.K.C.; Mestre, O.; Aguilar, E.; Auer, I.; Guijarro, J.A.; Domonkos, P.; Vertacnik, G.; Szentimrey, T.; Stepanek, P.; Zahradnicek, P.; et al. Benchmarking homogenization algorithms for monthly data. *Clim. Past* **2012**, *8*, 89–115. [CrossRef]

10. Bevis, M.; Bussinger, S.; Herring, T.A.; Rocken, C.; Anthes, R.A.; Ware, R.H. GPS Meteorology: Remote Sensing of Atmospheric Water Vapor Using the Global Positioning System. *J. Geophys. Res.* **1992**, *97*, 15787–15801. [CrossRef]

11. Bock, O.; Bosser, P.; Bourcy, T.; David, L.; Goutail, F.; Hoareau, C.; Keckhut, P.; Legain, D.; Pazmino, A.; Pelon, J.; et al. Accuracy assessment of water vapour measurements from in situ and remote sensing techniques during the DEMEVAP 2011 campaign at OHP. *Atmos. Meas. Tech.* **2013**, *6*, 2777–2802. [CrossRef]

12. Van Malderen, R.; Brenot, H.; Pottiaux, E.; Beirle, S.; Hermans, C.; De Maziere, M.; Wagner, T.; De Backer, H.; Bruyninx, C. A multi-site intercomparison of integrated water vapour observations for climate change analysis. *Atmos. Meas. Tech.* **2014**, *7*, 2487–2512. [CrossRef]

13. Ning, T.; Wickert, J.; Deng, Z.; Heise, S.; Dick, G.; Vey, S.; Schöne, T. Homogenized Time Series of the Atmospheric Water Vapor Content Obtained from the GNSS Reprocessed Data. *J. Clim.* **2016**, *29*, 2443–2456. [CrossRef]

14. Parracho, A.C.; Bock, O.; Bastin, S. Global IWV trends and variability in atmospheric reanalyses and GPS observations. *Atmos. Chem. Phys.* **2018**, *18*, 16213–16237. [CrossRef]

15. Dee, D.P.; Uppala, S.M.; Simmons, A.J.; Berrisford, P.; Poli, P.; Kobayashi, S.; Andrae, U.; Balmaseda, M.A.; Balsamo, G.; Bauer, P.; et al. The ERA-Interim reanalysis: Configuration and performance of the data assimilation system. *Q. J. R. Meteorol. Soc.* **2011**, *137*, 553–597. 10.1002/qj.828. [CrossRef]
16. Gelaro, R.; McCarty, W.; Suárez, M.J.; Todling, R.; Molod, A.; Takacs, L.; Randles, C.A.; Darmenov, A.; Bosilovich, M.G.; Reichle, R.; et al. The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2). *J. Clim.* **2017**, *30*, 5419–5454. [CrossRef] [PubMed]

17. Bock, O.; Parracho, A. Consistency and representativeness of integrated water vapour from ground-based GPS observations and ERA-Interim reanalysis. *Atmos. Chem. Phys.* **2019**, *19*, 9453–9468. [CrossRef]

18. Hersbach, H.; Bell, B.; Berrisford, P.; Hirahara, S.; Horányi, A.; Muñoz-Sabater, J.; Nicolas, J.; Peubey, C.; Radu, R.; Schepers, D.; et al. The ERA5 Global Reanalysis. *Q. J. R. Meteorol. Soc.* **2020**, *146*, 1999–2046. [CrossRef]

19. Vey, S.; Dietrich, R.; Fritsche, M.; Rülke, A.; Steigenberger, P.; Rothacher, M. On the homogeneity and interpretation of precipitable water time series derived from global GPS observations. *J. Geophys. Res. Atmos.* **2009**, *114*. [CrossRef]

20. Domonkos, P. Adapted Caussinus-Mestre Algorithm for Networks of Temperature series (ACMANT). *Int. J. Geosci.* **2011**, *2*, 293–309. [CrossRef]

21. Quarello, A. Development of New Homogenisation Methods for GNSS Atmospheric Data. Application to the Analysis of Climate Trends and Variability. Ph.D. Thesis, Sorbonne Universite, Paris, France, 2020.

22. Van Malderen, R.; Pottiaux, E.; Klos, A.; Domonkos, P.; Elias, M.; Ning, T.; Bock, O.; Guijarro, J.; Alshawaf, F.; Hoseini, M.; et al. Homogenizing GPS Integrated Water Vapor Time Series: Benchmarking Break Detection Methods on Synthetic Data Sets. *Earth Space Sci.* **2020**, *7*, e2020EA001121. [CrossRef]

23. Teunissen, P.J.; Montenbruck, O. (Eds.) *Springer Handbook of Global Navigation Satellite Systems*; Springer International Publishing: Berlin/Heidelberg, Germany, 2017. [CrossRef]

24. Petit, G.; Luzum, B. *IERS 2010 Conventions*; Technical Report; IERS: Frankfurt-am-Main, Germany, 2010.

25. Guerova, G.; Jones, J.; Douša, J.; Dick, G.; de Haan, S.; Pottiaux, E.; Bock, O.; Pacione, R.; Elgered, G.; Vedel, H.; et al. Review of the state of the art and future prospects of the ground-based GNSS meteorology in Europe. *Atmos. Meas. Tech.* **2016**, *9*, 5385–5406. [CrossRef]

26. Davis, J.L.; Herring, T.A.; Shapiro, I.I.; Rogers, A.E.E.; Elgered, G. Geodesy by radio interferometry: Effects of atmospheric modeling errors on estimates of baseline length. *Radio Sci.* **1985**, *20*, 1593–1607. [CrossRef]

27. Bock, O. Standardization of ZTD screening and IWV conversion. In *Advanced GNSS Tropospheric Products for Monitoring Severe Weather Events and Climate: COST Action ES1206 Final Action Dissemination Report*; Jones, J., Guerova, G., Douša, J., Dick, G., de Haan, S., Pottiaux, E., Bock, O., Pacione, R., van Malderen, R., Eds.; Springer International Publishing AG: Cham, Switzerland. 2020; Chapter 5, pp. 314–324. [CrossRef]

28. Bock, O.; Bosser, P.; Flamant, C.; Doerflinger, E.; Jansen, F.; Fages, R.; Bony, S.; Schnitt, S. Integrated water vapour observations in the Caribbean arc from a network of ground-based GNSS receivers during EUREC⁴A. *Earth Syst. Sci. Data* **2021**, *13*, 2407–2436. [CrossRef]

29. Tregoning, P.; Herring, T.A. Impact of a priori zenith hydrostatic delay errors on GPS estimates of station heights and zenith total delays. *Geophys. Res. Lett.* **2006**, *33*, 10.1029/2006GL027706. [CrossRef]

30. Schmid, R.; Steigenberger, P.; Gendt, G.; Ge, M.; Rothacher, M. Generation of a consistent absolute phase-center correction model for GPS receiver and satellite antennas. *J. Geod.* **2007**, *81*, 781–798. [CrossRef]

31. Schmid, R.; Dach, R.; Collilieux, X.; Jäggi, A.; Schmitz, M.; Dilssner, F. Absolute IGS antenna phase center model igs08.atx: Status and potential improvements. *J. Geod.* **2015**, *90*, 343–364. [CrossRef]

32. Ning, T.; Wang, J.; Elgered, G.; Dick, G.; Wickert, J.; Bradke, M.; Sommer, M.; Querel, R.; Smale, D. The uncertainty of the atmospheric integrated water vapour estimated from GNSS observations. *Atmos. Meas. Tech.* **2016**, *9*, 79–92. [CrossRef]

33. Zumberge, J.F.; Heflin, M.B.; Jefferson, D.C.; Watkins, M.M. Precise point positioning for the efficient and robust analysis of GPS data from large networks. *J. Geophys. Res.* **1997**, *102*, 5005–5017. doi:10.1029/96JB03860. [CrossRef]

34. Byun, S.H.; Bar-Sever, Y.E. A new type of troposphere zenith path delay product of the international GNSS service. *J. Geod.* **2009**, *83*, 1–7. [CrossRef]

35. Boehm, J.; Niell, A.E.; Tregoning, P.; Schuh, H. The Global Mapping Function (GMF): A new empirical mapping function based on numerical weather model data. *Geophys. Res. Lett.* **2006**, *33*, L07304. [CrossRef]

36. Niell, A.E. Global mapping functions for the atmosphere delay at radio wavelengths. *J. Geophys. Res. Solid Earth* **1996**, *101*, 3227–3246. [CrossRef]

37. Susnik, A.; Dach, R.; Villiger, A.; Maier, A.; Arnold, D.; Schaer, S.; Jäggi, A. CODE Reprocessing Product Series. 2016. Available online: https://boris.unibe.ch/80011/ (accessed on 25 August 2021). [CrossRef]

38. Dach, R.; Lutz, S.; Walser, P.; Fridez, P. Bernese GNSS Software Version 5.2; User Manual. 2015. Available online: https://boris.unibe.ch/72297/ (accessed on 8 May 2012). [CrossRef]

39. Dach, R.; Schaer, S.; Arnold, D.; Orliac, E.; Prange, L.; Susnik, A.; Villiger, A.; Jäggi, A. CODE Final Product Series for the IGS; 2018. Available online: https://boris.unibe.ch/119490/ (accessed on 25 August 2021). [CrossRef]

40. Boehm, J.; Werl, B.; Schuh, H. Troposphere mapping functions for GPS and very long baseline interferometry from European Centre for Medium-Range Weather Forecasts operational analysis data. *J. Geophys. Res.* **2006**, *11*, 2406. [CrossRef]

41. Bock, O.; Pacione, R.; Ahmed, F.; Araszkiewicz, A.; Bałdysz, Z.; Balidakis, K.; Barroso, C.; Bastin, S.; Beirle, S.; Berckmans, J.; et al. Use of GNSS Tropospheric Products for Climate Monitoring (Working Group 3). In *Advanced GNSS Tropospheric Products for Monitoring Severe Weather Events and Climate*; Jones, J., Guerova, G., Douša, J., Dick, G., de Haan, S., Pottiaux, E., Bock, O., Pacione, R., van Malderen, R., Eds.; Springer International Publishing: Cham, Switzerland, 2020; pp. 267–402
42. Stepniak, K.; Bock, O.; Wielgosz, P. Reduction of ZTD outliers through improved GNSS data processing and screening strategies. *Atmos. Meas. Tech.* **2018**, *11*, 1347–1361. [CrossRef]

43. Bock, O. GPS Data: Daily and Monthly Reprocessed IWV Data from 120 Global GPS Stations, Version 1.2. 2016. Available online: https://observations.ipsl.fr/espri/metadata/global_gps_iwv_v1.2.html/ (accessed on 25 August 2021). [CrossRef]

44. Bock, O. Global GNSS IWV Data at 436 Stations over the 1994–2018 Period. 2019. Available online: https://www.aeris-data.fr/metadata/metadata/?5829bd3d-4593-4e66-bce2-c8c6311360af (accessed on 25 August 2012). [CrossRef]

45. Quarello, A.; Bock, O.; Lebarbier, E. A new segmentation method for the homogenisation of GNSS-derived IWV time-series. *arXiv* **2020**, arXiv:2005.04683.

46. Birgé, L.; Massart, P. Gaussian model selection. *J. Eur. Math. Soc.* **2001**, *3*, 203–268. [CrossRef]

47. Lebarbier, E. Detecting Multiple Change-Points in the Mean of Gaussian Process by Model Selection. *Signal Process.* **2005**, *85*, 717–736. [CrossRef]

48. Weatherhead, E.C.; Reinsel, G.C.; Tiao, G.C.; Meng, X.; Choi, D.; Cheang, W.; Keller, T.; DeLuisi, J.; Wuebbles, D.J.; Kerr, J.B.; et al. Factors affecting the detection of trends: Statistical considerations and applications to environmental data. *J. Geophys. Res. Atmos.* **1998**, *103*, 17149–17161. [CrossRef]

49. Koulali, A.; Clarke, P.J. Effect of antenna snow intrusion on vertical GPS position time series in Antarctica. *J. Geod.* **2020**, *94*. [CrossRef]

50. Dunn, R.J.H.; Stanitski, D.M.; Gobron, N.; Willett, K.M. Global Climate [in "State of the Climate in 2019"). *Bull. Am. Meteorol. Soc.* **2020**, *101*, S9–S127. [CrossRef]

51. Ulbrich, U.; Lionello, P.; Belušić, D.; Jacobeit, J.; Knippertz, P.; Kuglitsch, F.G.; Leckebusch, G.C.; Luterbacher, J.; Maugeri, M.; Maheras, P.; et al. 5—Climate of the Mediterranean: Synoptic Patterns, Temperature, Precipitation, Winds, and Their Extremes. In *The Climate of the Mediterranean Region*; Lionello, P., Ed.; Elsevier: Oxford, UK, 2012; pp. 301–346.
# Chapter 4

# Development of the attribution method

## 4.1 Paper No. 2: A statistical method for the attribution of change-points in segmented IWV difference time series

This chapter comprises two main sections. The first one corresponds to a paper, submitted to the International Journal of Climatology (presently under review) which describes the attribution method developed in this thesis and the results on a real data set of 81 main stations. The second section presents additional research conducted alongside the paper's topic. It includes namely a comprehensive assessment of the noise model identification and parameter estimation tools used in the paper, based on numerical simulations. It also explores some alternative methodologies that may help to further increase the efficiency of several parts (regression, prediction) of the attribution method.

### 4.1.1 Abstract

Many segmentation or change-point detection methods used for the homogenization of climate time series from station data use a reference series against which the station data is compared. The main advantage of this approach is to remove the common climate signal and thus improve the detection power. One drawback is that it is difficult to decide whether the detected change-point is due to the main series or to the reference. A so-called attribution procedure is typically applied in a post-processing step for each detected change-point and each main station. This paper describes a new statistical method for the attribution of detected change-points in GNSS minus reanalysis series of Integrated Water Vapour (IWV). It works by combining the GNSS and reanalysis data from the main station with similar data from one or several nearby stations. The paired data from one main station and one nearby station form a set of four base series. The six series of differences formed from these four base series are tested for a significant jump at the time of the detected change-point in the main station. A statistical predictive rule is used to attribute the change-points in the four base series from the six test results. Original aspects of our method are: 1) the significance test, which is based on a generalized linear regression approach, taking both heteroscedasticity and autocorrelation into account; 2) the predictive rule, which uses a machine learning method and is constructed from the test results obtained with the real data by using a resampling strategy. Four popular machine learning methods have been compared using cross-validation and the best one was applied to a real data set (49 main stations with 114 change-points). The results depend on the choice of the test significance level and the aggregation method combining the prediction results when
several nearby stations are available. We find that 62% of the change-points are attributed to GNSS, 19% to the reanalysis, and 10% are due to coincident detections.

### 4.1.2 Introduction

Long records of climate observations are crucial for monitoring regional and global climate change and understanding the underlying climate processes (Trenberth et al., 2013; Dunn et al., 2021). However, many observational climate data are affected by inhomogeneities due to changes in instrumentation, in station location, in observation and processing methods, and/or in the measurement conditions around the station (Peterson et al., 1998a; Mitchell and Jones, 2005; Menne et al., 2009). Inhomogeneities often take the form of abrupt changes, which are detrimental to estimating trends and multi-scale climate variability (Jones et al., 1986; Easterling and Peterson, 1995). Various homogenization methods have been developed for the detection and the correction of such change-points in the context of climate data analysis (Peterson et al., 1998b; Reeves et al., 2007; Costa and Soares, 2009; Venema et al., 2012). The change-point detection step, also called segmentation, can be performed in two classical ways, using either a statistical test (e.g. (Alexandersson, 1986; Easterling and Peterson, 1995; Menne and Williams, 2005; Menne and Williams, 2009; Szentimrey, 2008; Wang et al., 2010)) or a penalized likelihood approach (e.g. (Caussinus and Mestre, 2004; Lu et al., 2010; Domonkos, 2011; Mestre et al., 2013)). The former proceeds sequentially and detects one change-point at a time, which leads inevitably to a sub-optimal solution. On the other hand, the second approach estimates all the change-points at once, and is thus optimal or sub-optimal, depending on the search algorithm. When the whole parameter space is explored, such as with the dynamic programming algorithm, the method is optimal. Many climate segmentation methods are used on differenced data, where the target series is differenced with respect to a reference series. Using differenced series helps to remove the common climate signal and improves thus the detection power of the segmentation method. However, one drawback of this approach is that any detected change-point can be either due to the target series or to the reference series, if the latter is not homogeneous. In this so-called relative homogenization approach, the reference series has been traditionally constructed by compositing the series from several nearby stations (Alexandersson, 1986; Menne and Williams, 2005; Guijarro, 2011). Compositing relaxes the need for homogeneous reference series thanks to the averaging from many nearby stations, such that the detected change-points can be attributed with good confidence to the target series. Unfortunately, in practice, composited reference series often contain non-negligible inhomogeneities. Another approach based on the pairwise comparison of individual series has been shown to be an interesting alternative (Caussinus and Mestre, 2004; Menne and Williams, 2009; Mestre et al., 2013; Domonkos et al., 2021). In this approach, the change-points from the target and reference series are disentangled in a post-segmentation step, referred to as "attribution". In Caussinus and Mestre (2004), the attribution step is done manually, by using both statistical inference and historical information (station metadata) in an iterative way. In Menne and Williams (2009), an automatic procedure is proposed that attributes a change-point to the station with the highest overall count of detections. This method also uses station metadata when available. It assigns the detected change-points to the nearest known event from the station history within some confidence limit. In Mestre et al. (2013), both a semi-automatic method similar to Caussinus and Mestre (2004) and a fully automatic method based on the joint detection of all series at once are implemented, but the latter is not a relative homogenization method and is thus not recommended (Domonkos, 2021).
The above-mentioned attribution methods generally require many nearby stations in order to find out which station is the cause of the detected change-point. They also operate in an iterative way, alternating the segmentation and attribution steps, and perform better when history information is included. In this work, we propose a new attribution method which significantly relaxes these constraints. Firstly, it works even if only one nearby station is available, which makes it usable in data sparse networks. Secondly, it operates in a post-processing mode, meaning that it uses as input the results from the segmentation step and does not need to iterate, although iterations may possibly help to make it more robust. Thirdly, it uses a predictive rule based on machine learning to attribute the cause of change-point among the tested series. The latter is trained in a preliminary step based on real data and is thus optimized for the specific data of interest.

The data of interest in this work is Integrated Water Vapor (IWV) derived from Global Navigation Satellite System (GNSS) measurements (Bock, 2019) and from the fifth ECMWF reanalysis (ERA5) (Hersbach et al., 2020). Because the global GNSS data set is relatively sparse, the reanalysis is used as a reference to form the target minus reference difference series necessary to the segmentation step (Ning et al., 2016; Bock et al., 2019; Van Malderen et al., 2020; Nguyen et al., 2021; Quarello et al., 2022). The primary goal of the attribution is thus to determine whether any change-point detected by the segmentation is due to the GNSS series or to the reanalysis series. Although the long-term stability of reanalyses is sometimes questioned (Thorne and Vose, 2010), they are rather homogeneous, especially in recent years (roughly after 2000) (Sterl, 2004; Kozubek et al., 2020). Inhomogeneities in reanalyses are mainly suspected when changes occur in the global observing system, e.g. the start or end of satellite missions which data are assimilated (Rienecker et al., 2011; Schroeder et al., 2016).

In this study, the GNSS minus reanalysis data are segmented with the "GNSSseg" method developed by Quarello (2020) which is based on a penalized likelihood approach. It detects abrupt changes in the mean in the presence of a periodic (seasonal) bias and a periodic variance (on a monthly basis) and is available on the CRAN (https://cran.r-project.org/web/packages/GNSSseg/index.html). It has been used in a benchmark exercise where it was ranked one of the best among 8 segmentation tools (Van Malderen et al., 2020). The attribution step was not necessary in that simulation study because the reference series was homogeneous.

Figure 4.1 helps to explain the idea of the attribution method proposed in this paper. Let us denote by G and E the GNSS and ERA5 reanalysis series of the main station, respectively, and G' and E' those from a nearby station. We denote by $t_1$, $t_2$, and $t_3$, the change-points detected by the segmentation method in the G-E series, and by $t'_1$ and $t'_2$, the change-points detected in the G'-E' series. These change-points have jumps in the mean of +1, -1, and -0.5 signal unit for the G-E series, and +1 and -0.5 signal unit for the G'-E' series. Note that in this sketch, the time period of the G'-E' series covers all the change-points of the main station, but in practice, several nearby stations may be necessary. The positions of the change-points illustrate different typical situations encountered in practice with our data. The first change-point in the nearby station, $t'_1$, is quite far from all the change-points detected in the main station. This illustrates the fact that the causes of inhomogeneities in GNSS data are primarily station-specific, i.e. coincident change-points in G and G' are expected to be rare. On the other hand, $t'_2$ is close in time to $t_3$ which illustrates an inhomogeneity in the reanalysis data with a large spatial extension, i.e. impacting both E and E'. Real data often contain data gaps which are due to instrumental failures
leading eventually to an equipment change and possibly to an inhomogeneity. This situation is illustrated with a gap after $t_2$ in the G-E series. The likelihood of these different situations is summarized in the following "empirical" rules which will help interpreting the features seen in the difference series:

(R1) it is unlikely that change-points in two different GNSS series (here G and G') occur at the same time because they are station-specific in nature (e.g. hardware failure, equipment change, local environmental change).

(R2) on the other hand, it is likely that change-points in the reanalysis occur simultaneously at the main and nearby sites (impacting E and E') because they are expected to have a large spatial extent (e.g. due to a change in assimilation of satellite measurements).

<figure><img src="image_30.png" /></figure>

Figure 4.1 – Schematic view of three paired series of differences, G-E, G'-E', and G-G', where G and E are the series from the main station, and G' and E' the series from the nearby station. Change-points detected by the segmentation method in the main (nearby) station are noted $t_k (t'_k)$ and are indicated by the vertical solid (dotted) lines. By convention, $t_0 (t'_0)$ and $t_K (t'_K)$ refer to the time of the first and last observation, respectively, in the main (nearby) station. The colored horizontal lines with arrows indicate the segments on the left and the right of the change-points that are used to estimate the deterministic and stochastic parameters of the regression model. This figure is discussed the Introduction.

Inspection of the first two series of differences in Figure 4.1 in the light of these rules suggests that $t_1$ is likely due to a +1 jump in G, $t'_1$ is likely due to a +1 jump in G', $t_2$ is likely due to a -1 jump in G, and $t'_2$ and $t_3$ are likely due to a -0.5 jump in both E and in E'. However, to confirm these guesses, we need to inspect additional series of differences combining more of the four base series (G, E, G', and E'). The lower plot in Figure 4.1 shows the G-G' series. It is straightforward, by the same reasoning, to confirm the guessed interpretation of the former two series. In a more general procedure, we would use all six combinations of the four base series and by deduce which of the four base series is/are the cause of the jumps observed in the multiple differenced series.
<table><caption>Table 4.1 – BS Table: Theoretical configurations of the jumps in the four base series (G, E, G', E'), coded as: 0=no jump, -1=downward jump, +1=upward jump, and the associated conditional and joint probabilities. SD Table: resulting jumps in the series of differences, coded on five levels (-2, -1, 0, 1, 2). RSD Table: similar to SD Table, but values are restricted to three levels only (-1, 0, +1). In theory the six test results could be searched in the RSD Table and the corresponding configuration of the jumps in the four based series attributed. Duplicated results in the SD and RSD Tables are highlighted with colored background; the configurations of lower joint probabilities are removed, leaving 46 cases in the SD Table and 38 cases in the RSD Table. See Appendix B for further details.</caption><thead><tr><th colspan="5">BS Table</th><th colspan="5">SD Table</th><th colspan="5">RSD table</th></tr><tr><th colspan="3">Jump in the base series</th><th>Conditional probability</th><th>Joint probability</th><th colspan="5">Jump in the series of difference</th><th colspan="5">Restricted jump in the series of difference</th></tr><tr><th>G</th><th>E</th><th>G'</th><th>E'</th><th>P(G', E', G, E)</th><th>P(G, E, G', E')</th><th>G-E</th><th>G-G'</th><th>G-E'</th><th>E-E'</th><th>G'-E'</th><th>G'-E'</th><th>G-E</th><th>G-G'</th><th>G-E'</th><th>E-E'</th><th>G'-E'</th><th>G'-E'</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0.8</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>2</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0.045</td><td>2</td><td>1</td><td>0</td><td>-1</td><td>-1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>3</td><td>1</td><td>0</td><td>0</td><td>-1</td><td>0.045</td><td>3</td><td>1</td><td>2</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>4</td><td>1</td><td>0</td><td>1</td><td>0</td><td>0.045</td><td>4</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>5</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0.0025</td><td>4</td><td>1</td><td>0</td><td>0</td><td>-1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>6</td><td>1</td><td>0</td><td>1</td><td>-1</td><td>0.0025</td><td>5</td><td>1</td><td>2</td><td>1</td><td>0</td><td>-1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>7</td><td>1</td><td>0</td><td>-1</td><td>0</td><td>0.045</td><td>5</td><td>1</td><td>2</td><td>1</td><td>0</td><td>-1</td><td>-1</td><td>-1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><td>8</td><td>1</td><td>0</td><td>-1</td><td>1</td><td>0.0025</td><td>6</td><td>1</td><td>2</td><td>0</td><td>-1</td><td>-2</td><td>-1</td><td>-1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><td>9</td><td>1</td><td>0</td><td>-1</td><td>-1</td><td>0.0025</td><td>6</td><td>1</td><td>2</td><td>1</td><td>0</td><td>-1</td><td>-1</td><td>-1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><td>10</td><td>0</td><td>-1</td><td>0</td><td>0</td><td>0.045</td><td>8</td><td>1</td><td>0</td><td>0</td><td>-1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>11</td><td>0</td><td>-1</td><td>0</td><td>1</td><td>0.045</td><td>9</td><td>1</td><td>0</td><td>-1</td><td>-2</td><td>-1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>12</td><td>0</td><td>-1</td><td>0</td><td>-1</td><td>0.81</td><td>10</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>13</td><td>0</td><td>-1</td><td>1</td><td>0</td><td>0.0025</td><td>11</td><td>1</td><td>-1</td><td>0</td><td>-1</td><td>1</td><td>2</td><td>2</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><td>14</td><td>0</td><td>-1</td><td>1</td><td>1</td><td>0.0025</td><td>12</td><td>1</td><td>-1</td><td>-1</td><td>-2</td><td>0</td><td>2</td><td>2</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><td>15</td><td>0</td><td>-1</td><td>1</td><td>-1</td><td>0.045</td><td>13</td><td>1</td><td>-1</td><td>1</td><td>0</td><td>2</td><td>2</td><td>2</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><td>16</td><td>0</td><td>-1</td><td>-1</td><td>0</td><td>0.0025</td><td>14</td><td>1</td><td>1</td><td>0</td><td>-1</td><td>-1</td><td>-2</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>17</td><td>0</td><td>-1</td><td>-1</td><td>1</td><td>0.045</td><td>14</td><td>1</td><td>1</td><td>-1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>19</td><td>-1</td><td>0</td><td>0</td><td>0</td><td>0.81</td><td>15</td><td>-1</td><td>-1</td><td>-1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>-1</td><td>-1</td><td>0</td><td>0</td></tr><tr><td>20</td><td>-1</td><td>0</td><td>0</td><td>-1</td><td>0.045</td><td>16</td><td>-1</td><td>-1</td><td>-2</td><td>-1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>-1</td><td>0</td><td>0</td></tr><tr><td>21</td><td>-1</td><td>0</td><td>0</td><td>1</td><td>0.045</td><td>17</td><td>-1</td><td>-1</td><td>0</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>22</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0.045</td><td>18</td><td>-1</td><td>-2</td><td>0</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>23</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0.0025</td><td>19</td><td>-1</td><td>-2</td><td>2</td><td>-1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>24</td><td>-1</td><td>0</td><td>1</td><td>-1</td><td>0.0025</td><td>20</td><td>-1</td><td>-2</td><td>0</td><td>1</td><td>2</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>25</td><td>-1</td><td>0</td><td>-1</td><td>1</td><td>0.045</td><td>21</td><td>-1</td><td>0</td><td>-2</td><td>-1</td><td>-2</td><td>-1</td><td>-1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>26</td><td>-1</td><td>0</td><td>-1</td><td>1</td><td>0.0025</td><td>21</td><td>-1</td><td>0</td><td>-2</td><td>-1</td><td>-2</td><td>-1</td><td>-1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>27</td><td>-1</td><td>0</td><td>-1</td><td>-1</td><td>0.0025</td><td>22</td><td>-1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>28</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0.045</td><td>22</td><td>-1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>-1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>29</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0.81</td><td>23</td><td>-1</td><td>0</td><td>-1</td><td>0</td><td>-1</td><td>-1</td><td>-1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>30</td><td>0</td><td>1</td><td>0</td><td>-1</td><td>0.045</td><td>24</td><td>-1</td><td>0</td><td>1</td><td>2</td><td>1</td><td>-1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><td>31</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0.0025</td><td>24</td><td>-1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><td>32</td><td>0</td><td>1</td><td>1</td><td>-1</td><td>0.045</td><td>25</td><td>-1</td><td>-1</td><td>-1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>33</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0.0025</td><td>26</td><td>-1</td><td>-1</td><td>1</td><td>2</td><td>2</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>34</td><td>0</td><td>1</td><td>-1</td><td>0</td><td>0.0025</td><td>27</td><td>-1</td><td>1</td><td>0</td><td>1</td><td>-1</td><td>-2</td><td>-2</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>35</td><td>0</td><td>1</td><td>-1</td><td>1</td><td>0.045</td><td>28</td><td>-1</td><td>1</td><td>1</td><td>2</td><td>0</td><td>-2</td><td>-2</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>36</td><td>0</td><td>1</td><td>-1</td><td>-1</td><td>0.0025</td><td>29</td><td>1</td><td>1</td><td>-1</td><td>0</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>37</td><td>1</td><td>-1</td><td>0</td><td>0</td><td>0.045</td><td>30</td><td>2</td><td>1</td><td>0</td><td>-2</td><td>-1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><td>38</td><td>1</td><td>-1</td><td>0</td><td>1</td><td>0.045</td><td>31</td><td>2</td><td>1</td><td>2</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><td>39</td><td>1</td><td>-1</td><td>0</td><td>-1</td><td>0.81</td><td>32</td><td>2</td><td>0</td><td>0</td><td>-2</td><td>0</td><td>2</td><td>2</td><td>2</td><td>1</td><td>1</td><td>0</td></tr><tr><td>40</td><td>1</td><td>-1</td><td>1</td><td>0</td><td>0.0025</td><td>33</td><td>2</td><td>0</td><td>0</td><td>-2</td><td>0</td><td>2</td><td>2</td><td>2</td><td>1</td><td>1</td><td>0</td></tr><tr><td>41</td><td>1</td><td>-1</td><td>1</td><td>1</td><td>0.0025</td><td>34</td><td>2</td><td>0</td><td>0</td><td>2</td><td>0</td><td>2</td><td>2</td><td>2</td><td>1</td><td>1</td><td>0</td></tr><tr><td>42</td><td>1</td><td>-1</td><td>-1</td><td>1</td><td>0.045</td><td>35</td><td>2</td><td>2</td><td>1</td><td>-1</td><td>-1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>43</td><td>1</td><td>-1</td><td>-1</td><td>0</td><td>0.0025</td><td>36</td><td>2</td><td>2</td><td>0</td><td>-2</td><td>-2</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>44</td><td>1</td><td>-1</td><td>-1</td><td>1</td><td>0.0025</td><td>37</td><td>2</td><td>2</td><td>2</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>45</td><td>1</td><td>-1</td><td>-1</td><td>-1</td><td>0.045</td><td>38</td><td>-2</td><td>-1</td><td>-1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>46</td><td>-1</td><td>1</td><td>0</td><td>0</td><td>0.045</td><td>39</td><td>-2</td><td>-1</td><td>-2</td><td>0</td><td>-1</td><td>-1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>47</td><td>-1</td><td>1</td><td>0</td><td>1</td><td>0.81</td><td>40</td><td>-2</td><td>-1</td><td>0</td><td>2</td><td>1</td><td>-1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>48</td><td>-1</td><td>1</td><td>0</td><td>-1</td><td>0.045</td><td>41</td><td>-2</td><td>-2</td><td>-1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>49</td><td>-1</td><td>1</td><td>1</td><td>0</td><td>0.0025</td><td>42</td><td>-2</td><td>-2</td><td>0</td><td>2</td><td>2</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>50</td><td>-1</td><td>1</td><td>1</td><td>1</td><td>0.045</td><td>43</td><td>-2</td><td>-2</td><td>0</td><td>1</td><td>-1</td><td>-1</td><td>-1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>51</td><td>-1</td><td>1</td><td>1</td><td>-1</td><td>0.0025</td><td>44</td><td>-2</td><td>-2</td><td>0</td><td>1</td><td>-1</td><td>-1</td><td>-1</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>52</td><td>-1</td><td>1</td><td>-1</td><td>1</td><td>0.0025</td><td>45</td><td>-2</td><td>0</td><td>0</td><td>2</td><td>0</td><td>-2</td><td>-2</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>53</td><td>-1</td><td>1</td><td>-1</td><td>1</td><td>0.045</td><td>46</td><td>-2</td><td>0</td><td>0</td><td>2</td><td>0</td><td>-2</td><td>-2</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>54</td><td>-1</td><td>1</td><td>-1</td><td>-1</td><td>0.0025</td><td>47</td><td>-2</td><td>0</td><td>0</td><td>2</td><td>0</td><td>-2</td><td>-2</td><td>1</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

The BS Table in Table 4.1 displays all the relevant combinations of jumps/no jumps in the four base series. The corresponding "theoretical" test results for the series of difference are given in the SD Table. The RSD Table shows the corresponding "practical" test results where the values are restricted to ± 1 (meaning an upward/downward jump) and 0 (meaning no jump, see Appendix A for further details). The latter two tables contain some duplicate configurations highlighted by the colored background which can be distinguished based on the joint probabilities (see Appendix B for the computation of the probabilities). In the end, we can distinguish 38 configurations in the RSD Table. In theory, we could thus use this table to attribute the jumps in the four base series based on the test results of the six series of differences. However, in practice, some of the test results may be wrong due to false negatives and/or false positives, and the combination of the six test results
would not be in the table. To overcome this difficulty, our attribution procedure builds on two main bricks. First, it uses an efficient test based on the Generalized Least Squares (GLS) method. This method is known to have higher detection power than most traditional regression methods in the presence of heteroscedastic and autocorrelated noise as is the case with our data. Second, a predictive rule is constructed based on the tests results from real data, using a machine learning algorithm. This is an efficient way to predict the most likely solution when the combination of the six test results is not in the table.

In subsection 4.1.3, we describe the stochastic properties of our data set composed of IWV time series from ground-based GNSS data and the ERA5 data. We highlight the embedded heteroscedasticity and autocorrelation of these differenced IWV series. In subsection 4.1.4, we evaluate the power of several test methods for testing a fixed change in mean with simulated data mimicking the heteroscedasticity and autocorrelation properties of our real data. We show that the GLS approach is superior to the other methods. In subsection 4.1.5, we describe the method for the construction of the predictive rule, compare the performance of four popular machine learning methods, and present the results on a real data set. Subsection 4.1.6 discussed the results and concludes.

### 4.1.3 Data characterization

#### 4.1.3.1 Data sets

In this study, we use daily IWV data from GNSS observations and from the ERA5 reanalysis which have been segmented beforehand using the GNSSseg segmentation method (Quarello et al., 2022). Our main goal here is to attribute each of the change-points detected with GNSSseg either to the GNSS series or to the ERA5 series. Several past studies using similar data highlighted the presence of abrupt changes in the mean of GNSS series (Vey et al., 2009; Bock et al., 2014; Ning et al., 2016; Parracho et al., 2018) or in the reanalysis data (Schroeder et al., 2016; Ning et al., 2016; Parracho et al., 2018; Nguyen et al., 2021). Inhomogeneities in the GNSS data are mainly due to equipment changes and changes in the station's environment, but the magnitude of the jumps may also depend on the data processing procedure (Nguyen et al., 2021). In this work, we use reprocessed GNSS data from Center for Orbit Determination in Europe (CODE), covering the period from 1994 to 2014 (REPRO2015), extended until the end of 2018 by a consistent operational processing. Details of the processing are described in Nguyen et al. (2021) and references therein, and the data set is available from Bock (2019). These data are from a global network of 436 stations. Data from the ERA5 reanalysis have been extracted at the location of each station and the difference series, G-E, have been segmented using the GNSSseg package. For the purpose of the present study we selected 81 stations with the longest time series. These will be our "main stations". Nearby stations were searched with a distance limit of 200 km in horizontal and 500 m in vertical, but very few were found in the CODE data set. So we used instead the GNSS data reprocessed by the Nevada Geodetic Laboratory (NGL) which comprises nearly 20 thousand stations (Blewitt et al., 2018). The NGL data were converted to IWV, differenced with respect to ERA5, and passed through the GNSSseg segmentation as well. We ended up with 114 detected change-points in 49 main stations that can be tested with respect to 312 nearby stations, resulting in a total number of 494 main/nearby pairs.
### 4.1.3.2 Pre-processing

Before we form the six difference series and estimate the change in the mean, the IWV data are adjusted for the station height difference and screened for outliers. The IWV adjustment is done following the method described in Bock et al. (2022) with correction model coefficients estimated from ERA5 on a global grid. This step is important when the main station and nearby station are not at the same altitude. It impacts both the GNSS data and the ERA5 data as the latter are extracted at the height of each GNSS station, either main or nearby. Once the G' and E' data are made consistent with the G and E data, the six difference series are formed.

The screening consists in two steps. The first one is a classical outlier detection procedure in which the data points exceeding three standard deviations from the median are removed. In this procedure, the median and a robust standard deviation estimator (Yohai and Zamar, 1988) are computed in a sliding window of length +/- 60 days around the current point. To guarantee the representativeness and accuracy of the estimates in the presence of data gaps, a minimum number of 20 data points are required.

The second step consists in removing data in short segments (less than 80 days) pertaining to a cluster of change-points detected in the main station (see an example in Figure 3 in Quarello (2020)). This problem occurs occasionally in regions where the GNSS data and reanalysis data have a significant representativeness difference (Bock and Parracho, 2019). In such situation, we keep the first change-point but remove data between the first and the last change-points of the cluster. Similarly, when a change-point in the nearby station is very close (less than 10 days) to a change-point in the main station, we consider that they are both due to the same cause (most likely a change in the reanalysis) and we remove the data points in the nearby series between the two change-points. This case is illustrated in Figure 4.1 where the data between $t_3$ and $t'_2$ have been removed. In Figure 4.1, we also illustrate the case of a gap in the G-E series just after $t_2$ which may be due to a screened cluster. Note that, as a result of the screening, the number of data points in a difference series combining the main and nearby sites (e.g. G-G') is always less or equal to that of a collocated difference series (e.g. G-E or G'-E'). In the proposed test procedure, a minimum number of 200 consecutive points is required on each side of the change-points.

### 4.1.3.3 Characterization of the data and model building

The GNSS minus reanalysis difference series show usually strong heteroscedasticity and periodic (seasonal) biases, along with weak autocorrelation (Quarello et al., 2022). In the following, a series is modelled using the following regression model:

$$z_t = \mu_L + \delta x_t + s_t + e_t, \quad (4.1)$$

where $t$ refers to the time, $\mu_L$ is the mean of the signal on the left of the change-point, $\delta$ is the amplitude of the jump, $x_t$ is a step function ($x_t = 0$ if $t \le t_k$ and $1$ if $t > t_k$, where $t_k$ is the time of the change-point detected by the segmentation method), $s_t$ is the Fourier series, and $e_t$ is the noise term. For ease of notation, we use $t$ as the time index, with $t = 1, ..., n$, but in reality the data may contain gaps and the time values are not consecutive. To account for this, $t$ can be replaced by $t(i)$, with $i = 1, ..., n$. To account for both heteroscedasticity and autocorrelation, we follow José C. Pinheiro (2000) and represent $e_t$ as the product of two factors:
$$e_t = e_t^* \sigma_t, \qquad (4.2)$$

where $e_t^*$ represents a stationary autocorrelated process of unit variance and $\sigma_t^2$ is the time-varying variance of $e_t$, i.e. $var[e_t] = \sigma_t^2$. Preliminary investigation of our data showed that most of the time the noise model is well approximated by an AR(1). Other noise models such as MA(1), ARMA(1,1), and pure white noise occur sometimes. We tested also for higher order ARMA(p,q) models and they are very rare. We limit thus ourselves to the four possible ARMA(p,q) models, with $p, q \in \{0, 1\}$. Recall that an ARMA(1,1) model writes (Shumway and Stoffer, 2017):

$$e_t^* = \phi e_{t-1}^* + \theta w_{t-1} + w_t, \qquad (4.3)$$

where $w_t$ is a Gaussian white noise. The noise model identification and parameter estimation methods are described in the next subsection. Note that other stochastic models including periodic variations in the mean, heteroscedasticity and autocorrelation have been proposed by Lund et al. (1995).

<figure><img src="image_31.png" /><figcaption>Figure 4.2 – Top: GNSS minus ERA5 time series at station ALBH (Victoria, Canada), in gray, and estimated Fourier series, in black, for a long, homogeneous, segment (no change-point detected by the segmentation method). Middle: FGLS regression residuals (jagged curve) and moving median (smooth curve). Bottom: moving standard deviation illustrating the strong heteroscedasticity in the data.</figcaption></figure>
Figure 4.2 shows an example of a time series (jagged gray curve), with the estimated Fourier series (smooth black curve), the estimated standard deviation (SD), $\hat{\sigma}_t$ (black curve at bottom), and the regression residuals (jagged orange curve). The strong heteroscedasticity is obvious, and because it is not stationary, we used a moving window approach (similar to the outlier screening procedure described above) to estimate it.

Table 4.2 and 4.3 summarize the characteristics of our data set in terms of heteroscedasticity and noise structure, respectively, for all six series of difference (G-E, G-G'), for the main and all nearby stations. The results are sorted according to the distance between the main and the nearby stations (smaller or larger than 50 km). Regarding the heteroscedasticity, three groups can be identified when the distance between sites is small. The first group (G1) includes G-E and G'-E', i.e. the series with collocated data, which have moderate mean SD of 0.7 kg m⁻². The second group (G2) includes E-E' and G-G', i.e. the series comparing the same technique, which have the smallest mean SD (0.5 kg m⁻²). The last group (G3) involves data from non-collocated data and mixed techniques, and gets the largest mean SD. As the distance increases, the mean SD of series involving different sites increases, as expected from increased representativeness differences. Another striking feature is that the half-range of the variation in SD is around 70% for all six series, indicating that heteroscedasticity is a strong feature in our data.

<table><caption>Table 4.2 – Characterization of the heteroscedasticity in the real data from 494 main/nearby series. The table reports the mean and the half-range of the standard deviation for each of the six paired difference series. The mean values are sorted by distance.</caption><thead><tr><th rowspan="2">Distance</th><th colspan="2">Mean of SD</th><th rowspan="2">Half-range of SD (%)</th></tr><tr><th>&lt; 50 km</th><th>≥ 50 km</th></tr></thead><tbody><tr><td>G-E</td><td>0.7 ± 0.26</td><td></td><td>72 ± 20</td></tr><tr><td>G'-E'</td><td>0.66 ± 0.24</td><td></td><td>67 ± 19</td></tr><tr><td>G-G'</td><td>0.52 ± 0.17</td><td>1.31 ± 0.47</td><td>63 ± 21</td></tr><tr><td>E-E'</td><td>0.41 ± 0.17</td><td>1.26 ± 0.47</td><td>73 ± 26</td></tr><tr><td>G-E'</td><td>0.82 ± 0.21</td><td>1.38 ± 0.46</td><td>67 ± 21</td></tr><tr><td>G'-E</td><td>0.83 ± 0.26</td><td>1.39 ± 0.46</td><td>66 ± 20</td></tr></tbody></table>

<table><caption>Table 4.3 – Characterization of the autoregressive noise structure of the real data. The table reports the mean estimated coefficients of the noise model for each of the six paired difference series, sorted by distance.</caption><thead><tr><th rowspan="3">Distance</th><th colspan="6">&lt; 50km</th><th colspan="6">≥ 50km</th></tr><tr><th colspan="2">series</th><th colspan="2">ARMA(1,1)</th><th rowspan="2">phi</th><th rowspan="2">theta</th><th rowspan="2">AR(1)</th><th rowspan="2">MA(1)</th><th rowspan="2">ARMA(1,1)</th><th rowspan="2">phi</th><th rowspan="2">theta</th></tr><tr><th>Coefficients</th><th>phi</th><th>theta</th><th>phi</th><th>theta</th></tr></thead><tbody><tr><td>G-E</td><td>0.30</td><td>0.00</td><td>0.57</td><td>-0.32</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>G'-E'</td><td>0.31</td><td>0.22</td><td>0.59</td><td>-0.34</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>G-G'</td><td>0.33</td><td>0.19</td><td>0.65</td><td>-0.31</td><td>0.30</td><td>0.22</td><td>0.11</td><td>0.12</td><td></td><td></td><td></td><td></td></tr><tr><td>E-E'</td><td>0.31</td><td>0.21</td><td>0.34</td><td>0.23</td><td>0.29</td><td>0.20</td><td>0.25</td><td>0.20</td><td></td><td></td><td></td><td></td></tr><tr><td>G-E'</td><td>0.33</td><td>0.24</td><td>0.59</td><td>-0.24</td><td>0.31</td><td>0.21</td><td>0.29</td><td>0.08</td><td></td><td></td><td></td><td></td></tr><tr><td>G'-E</td><td>0.32</td><td>0.21</td><td>0.57</td><td>-0.28</td><td>0.30</td><td>0.22</td><td>0.18</td><td>0.21</td><td></td><td></td><td></td><td></td></tr></tbody></table>
Figure 4.3 shows the distributions of the noise models and of the estimated coefficients for the six differences, again sorted according to the distance. The AR(1) model is the dominant model, with a proportion between 50% and 80%, independently of the distance, while the white noise model is extremely rare. The proportion of MA(1) and ARMA(1,1) depends on the distance and the series: ARMA(1,1) is dominant for the collocated series (similar to the noise group G1), as well as for the series comparing the same technique (group G2), when the distance is small. On the opposite, when the distance is large, MA(1) becomes more frequent, like for the series mixing techniques and sites (group G3). The increase of the distance does thus not only increase the variance of the noise but changes also its nature. Another interesting aspect is the values of the coefficients. For the AR(1), they are very similar (around 0.3) for all series, regardless of the distance. Similarly, for the MA(1), they are very similar (around 0.2). More surprisingly, the estimated coefficients of the ARMA(1,1) model for the non-collocated series depend somewhat on the distance, with the exception of E-E'. Values of $\hat{\phi}$ and $\hat{\theta}$ are around 0.6 and -0.3, respectively, for the collocated series and when the distance is small, and around 0.2 for both coefficients, when the distance is large. For E-E', the values are always around 0.2. The ARMA(1,1) models with coefficients of opposite sign found at short distance suggest that in these cases the noise is a mixture of AR(1) and white noise (Shumway and Stoffer, 2017). When the distance increases, the moving average part becomes more important, which may be interpreted as a spatial/temporal averaging of the variability in the difference series. The mean values of the estimated coefficients are reported in Table 4.3.

<figure><img src="image_32.png" /></figure>

Figure 4.3 – Results of noise model identification in the real data. (a, b) Histogram of model types (white noise, AR(1), MA(1), ARMA(1,1)) selected with `auto.arima` function for each of the six series of differences (G-E, G-E', etc.); the bar heights show the percentage (y-axis) of cases for each series (x-axis), the number of cases is indicated on the top of each bar. (c, d) noise model coefficients, $\hat{\phi}$ and $\hat{\theta}$, estimated with `arima` function, for each model. Results are sorted according to the distance between the main and the nearby stations, (a, c) smaller than 50 km, (b, d) larger than 50 km.
### 4.1.4 Proposed tests for a fixed change-point

#### 4.1.4.1 Regression model and different tests

In this subsection, we compare different procedures to test the significance of the jumps associated to change-points detected by the segmentation (i.e. considered here at known position). The series modelled by Eq. (4.1), with specifications (4.2) and (4.3), is now rewritten in matrix form:

$$ z = X\beta + e, \qquad (4.4) $$

where $\beta$ includes the coefficients of the deterministic part of the model, $\beta = (\mu_L, \delta, a_1, ..., a_4, b_1, ..., b_4)'$, and $X$ includes the corresponding regressors. Here, the $a_l$ and $b_l$, $l = 1, ..., 4$, are the coefficients of a Fourier series of order 4, and the corresponding regressors are $\cos(2\pi lt(i)/T)$ and $\sin(2\pi lt(i)/T)$, with $T = 365$ days, and $t(i)$ is the time of the $i^{th}$ observation, $z_i$, $i = 1, ..., n$. The noise vector, $e$ is assumed to be distributed as $N(0, \Sigma_0)$, where $\Sigma_0$ is the variance-covariance matrix describing the noise model. A correct specification of this matrix is crucial for the subsequent tests of the estimated coefficients of the deterministic model.

When $\Sigma_0$ is known, either Ordinary Least Squares (OLS) or Generalized Least Squares (GLS) methods can be used, for which the solutions for $\hat{\beta}$ and the corresponding variance-covariance matrices write:

$$ \hat{\beta}_{OLS} = (X'X)^{-1}X'z \quad \text{and} \quad \text{var}[\hat{\beta}_{OLS}] = (X'X)^{-1}(X'\Sigma_0X)(X'X)^{-1}, \qquad (4.5) $$

$$ \hat{\beta}_{GLS} = (X'\Sigma_0^{-1}X)^{-1}X'\Sigma_0^{-1}z \quad \text{and} \quad \text{var}[\hat{\beta}_{GLS}] = (X'\Sigma_0^{-1}X)^{-1}. \qquad (4.6) $$

Recall that, in the presence of heteroscedasticity and/or autocorrelation, GLS is the Best Linear Unbiased Estimator (BLUE) while OLS solution is not, although OLS remains unbiased.

In practice, $\Sigma_0$ is typically unknown and needs to be estimated. In the classical linear model (CLM) framework it is generally assumed that the data are independent and homoscedastic, e.g. $\Sigma_0 = \sigma_0^2 I_n$, where $I_n$ is the identity matrix. An estimator of $\text{var}[\hat{\beta}_{OLS}]$ is then simply

$$ \overline{\text{var}}[\hat{\beta}_{OLS}]_{CLM} = \hat{\sigma}_0^2 (X'X)^{-1}. \qquad (4.7) $$

where $\hat{\sigma}_0^2 = \hat{e}'\hat{e}/(n-k)$ is an unbiased estimator of the noise variance $\sigma_0^2$.

The CLM assumptions are not satisfied with our data and, despite the OLS solution (4.5) remains unbiased, the variance estimator is strongly biased and leads to significant inference errors. To solve this problem, some methods have been proposed in the literature. The two main ones, which we considered in this work, are:

- the so-called OLS-HAC that consists in still using the OLS solution but to consider a consistent estimator for the variance that is the Heteroscedasticity and Autocorrelation Consistent estimator (HAC) (White, 1980; Newey and West, 1986). This estimator is robust to the presence of heteroscedasticity and serial correlations of unknown form and has good asymptotic properties. The key idea is to estimate $M = X'\Sigma_0X$ instead of $\Sigma_0$, which is difficult to estimate due to its large size (it contains nominally $n(n+1)/2$ parameters).
The variance estimator writes:

$$ \widehat{var}[\hat{\beta}_{OLS}]_{HAC} = (X'X)^{-1}\hat{M}(X'X)^{-1}. \qquad (4.8) $$

A large class of HAC estimators is the non-parametric kernel estimators. A drawback of this type of estimator is that its performance varies with the choices of the kernel function and its bandwidth, but the advantage is that it does not require to specify the covariance structure and is computationally very fast. Here we consider the "Quadratic Spectral" kernel (Andrews, 1991) with its proposed optimal bandwidth value. This method is available the R package sandwich (Zeileis, 2006).

— the Feasible GLS (FGLS) which consists in estimating $\Sigma_0$, assuming it has a specific structure with a reduced number of parameters to be estimated, in the GLS solution Eq. (6). Denoting $\hat{\Sigma}_n$ the estimator of $\Sigma_0$, Eq. (6) becomes:

$$ \hat{\beta}_{FGLS} = (X'\hat{\Sigma}_n^{-1}X)^{-1}X'\hat{\Sigma}_n^{-1}z \quad \text{and} \quad \widehat{var}[\hat{\beta}_{FGLS}] = (X'\hat{\Sigma}_n^{-1}X)^{-1}. \qquad (4.9) $$

Following José C. Pinheiro (2000), we decompose $\Sigma_0$ as $\Sigma_0 = VCV$, where $V = diag(\sigma_t)$ and $C$ is the correlation matrix associated to the noise model (see Shumway and Stoffer (2017) for the formulation of matrix $C$ as a function of $\phi$ and $\theta$ for an ARMA(1,1)). With this parameterization, $\Sigma_0$ is described by a maximum of $n + 2$ coefficients.

Since both $\hat{\beta}_{FGLS}$ and $\widehat{var}[\hat{\beta}_{FGLS}]$ depend on $\hat{\Sigma}_n$, an iterative procedure is implemented. In order to stabilize the convergence and also to avoid over-fitting, we first select the best noise model among the four possible models (ARMA($p,q$) with $p, q \in \{0, 1\}$) using the auto.arima function in the R package forecast (Hyndman and Khandakar, 2008). The BIC criterion is used and is complemented by a test of the significance of the final model coefficients. Next, we use the following iterative procedure:

1. fit the OLS solution (4.5);

2. compute the noise variance $\hat{V}_n = diag(\hat{\sigma}_t^2)$ from the OLS residuals in a moving window (see subsection 4.1.3.3);

3. fit a preliminary FGLS solution (4.6) where $\Sigma_0$ is replaced by $\hat{V}_n$;

4. fit the ARMA model coefficients, $\hat{\phi}$ and $\hat{\theta}$, from the FGLS residuals;

5. compute $\hat{\Sigma}_n = \hat{V}_n \hat{C}_n \hat{V}_n$, where $\hat{C}_n$ is the correlation matrix of the ARMA model with the fitted coefficients;

6. fit the final FGLS solution (4.9);

7. repeat steps 3 to 6 until convergence.

In step 4, the parameters $\hat{\phi}$ and $\hat{\theta}$ of the ARMA noise structure are estimated by Maximum Likelihood with the function arima in R. In step 7, the convergence is tested from the difference of $\hat{\beta}$, $\hat{\phi}$, and $\hat{\theta}$, of two successive iterations. The maximum number of iterations is set to 10 and is never reached. For
the preliminary selection of the noise model, a simplified scheme including steps 1 to 4 is used, and auto.arima is used in step 4 instead of arima. Numerical simulations showed that auto.arima is able to select the correct model in 99 % of the cases when it is white noise, 93 % when it is MA(1), 90 % when it is AR(1), and 63 % when it is ARMA(1,1), for a sample size n ≥ 1000 with typical values for the coefficients (Figure 4.3).

Note that, compared to the OLS-HAC approach, the FGLS procedure, due to its iterative scheme, is computa-
tional time-demanding. However, when the noise model is correctly specified, the FGLS estimate of $\overline{var}[\hat{\beta}]$ is
more accurate than its OLS-HAC counterpart, and the power of the test is improved.

In the following we are interested in testing the null hypothesis $H_0 : \delta = 0$. The associated test statistic is:

$$
\tau_{\delta,\text{OLS-HAC}} = \frac{\hat{\delta}_{\text{OLS}}}{\hat{\sigma}_{\hat{\delta}_{\text{HAC}}}} \quad \text{for the OLS-HAC, and} \quad \tau_{\delta,\text{FGLS}} = \frac{\hat{\delta}_{\text{FGLS}}}{\hat{\sigma}_{\hat{\delta}_{\text{FGLS}}}} \quad \text{for the FGLS,} \quad (4.10)
$$

where $\hat{\sigma}_{\delta}$ is the estimated standard error of $\hat{\delta}$ which is extracted from $\overline{var}[\hat{\beta}_{\bullet}]$, and $\bullet = \text{OLS-HAC}$ or FGLS.
For the OLS-CLM and GLS estimators considered in the simulation study, similar statistics are computed using
their respective estimators for $\hat{\delta}$ and $\hat{\sigma}_{\delta}$ given by Eqs. (5) and (6).

In contrast to the HAC estimator, the asymptotic properties (unbiasedness and efficiency) of the FGLS estimator
are typically not known. Numerical simulations with the FGLS procedure described above show that $\hat{\beta}_{FGLS}$ is
not biased and its variance is very close to GLS, although slightly larger, but still smaller than HAC. Finally,
the distributions of both statistics are very close to $\mathcal{N}(0, 1)$ under $H_0$. The critical value, $\tau_{\alpha/2}$, associated to a
given significance level, $\alpha$, will thus be computed using the standard normal distribution.

4.1.4.2 Evaluation based on simulations

We conducted a large number of simulations, for different types of noise characteristics and sample sizes, to
assess the test methods introduced above. In these simulations, we modelled the noise heteroscedasticity by a
"raised cosine" function, $\sigma_t^2 = \sigma_m^2 - \sigma_v^2 \cos 2\pi t/T$, where $\sigma_m^2$ and $\sigma_v^2 \le \sigma_m^2$
modulation of the variance over one period, $T$, respectively.

The results are compared in Figure 4.4 for significance level, $\alpha = 0.05$. It is seen that the False Positive
Rate (FPR) stays generally fairly close to the nominal significance level for GLS, FGLS, and OLS-HAC, for
most autocorrelation and heteroscedasticity characteristics, except when the autocorrelation is very strong. In
contrast, the OLS-CLM method performs very badly when the data is autocorrelated or heteroscedastic, due
to the bias in its variance estimator (Eq. 4.7). The power of the test is measured by the True Positive Rate
(TPR). It depends on the jump amplitude, the sample size, and the noise characteristics, in addition to the fixed
significance level. When $\phi$ increases, it is generally observed that the TPR decreases, for all four methods. This
is due to an increase in $\hat{\sigma}_{\delta}$ which is sometimes interpreted as a reduction of the "equivalent sample size" (Zwiers
and Storch, 1995). The higher TPR of OLS-CLM is actually a consequence of its higher FPR and does not
indicate a good performance per se. In contrast, when the heteroscedasticity increases at constant $\phi$, GLS and
FGLS clearly outperform OLS-HAC. For OLS-HAC the TPR remains actually constant while for GLS and, to a
lesser extent FGLS, the TPR increases with stronger heteroscedasticity. This is explained by the fact that GLS
and FLGS take the weight of every observation into account; when heteroscedasticity is strong, observations with small errors have high weight which leads to a decrease in $\hat{\sigma}_{\delta}$. When both heteroscedasticity and auto-correlation are strong, the power of the test is low for all methods, but FLGS still performs better than OLS-HAC.

<figure><img src="image_33.png" /><figcaption>Figure 4.4 – False Positive Rate (FPR) and True Positive Rate (TPR) of jump detections with the regression model (1) and four different estimation methods, for three scenarios: (a, b) AR(1) noise of unit variance with $\phi = 0, ..., 0.90$; (c, d) heteroskedastic and AR(1) noise, with $\phi = 0.3$ and half-range of variance from 0 to 80 %; (e, f) heteroskedastic and AR(1) noise, with $\phi = 0, ..., 0.90$, and half-range of variance of 80 %. For TPR, the amplitude of the jump is fixed to 0.356, which corresponds to TPR=0.95 when $\phi = 0$. The sample size is $n=400$ and the number of simulated series is $m=1000$.</figcaption></figure>

Figure 4.5 provides additional insight into the power of the FGLS test as a function of the jump amplitude and sample size, for a typical noise configuration. A 75 % probability of detection is expected for jumps of 0.25 or larger with a sample size of $n=600$. Most of the G-E, G-E', and G-G' series in the real data actually fulfill this requirement (see Figure 4.6, further discussed in the next subsection). For stronger noise (usually due to larger distance) or smaller jumps, a larger sample size would be required to maintain a high detection probability.
<figure><img src="image_34.png" /><figcaption>Figure 4.5 – True Positive Rate (TPR) of jump detections with the FGLS method, as a function of jump amplitude and sample size, in case of heteroskedastic and AR(1) noise with $\phi = 0.3$, mean variance of 1, and half-range of variance of 80%.</figcaption></figure>

### 4.1.4.3 Application to real data

The FGLS procedure was applied to the series of differences from the 494 main/nearby pairs.

Figure 4.6 shows the distribution of estimated jump amplitudes, their standard errors, and the associated absolute t-values computed from Eq. 4.10, where the non-collocated series (G-G', G-E', E-E', and G'-E) are sorted by distance. Notably, the three series involving G have significantly larger median jump amplitudes (around 0.3 kg m⁻²) than the other three series, regardless of the distance. This result suggests that large jumps are occurring more often in the G series than in the E, E', or G' series. In G'-E' and G'-E, the median jump is small, as expected and expressed in our first rule stating that it is unlikely to have a coincident change-point in a nearby GNSS station when there is one detected in the main station. Additionally, a notable observation is that the median jump in E-E' is much larger at larger distance, which may be due to errors in the estimated jumps induced by a increased noise at larger distance. The variation of the SD of the noise with distance directly impacts also the jump standard error. The standard error of estimated jumps is notably smaller in collocated series, such as G-E and G'-E', as well as in non-collocated series at short distance. Furthermore, the standard errors in G-E' and G'-E (non-collocated series from different techniques) are slightly larger than in G-G' and E-E' (non-collocated series from the same technique) even at short distance, as also noticed in the three noise groups discussed in subsection 4.1.3.3. Finally, the t-values can be interpreted by considering the jump magnitudes and their standard errors. It is evident that the three series involving G yield larger t-values due to higher jump magnitudes. In contrast, the other three series have much smaller t-values, mainly because some of the large jumps at larger distance are damped by the larger standard errors. A common feature to all non-collocated series is that the t-values decrease with distance.

Figure 4.7 shows the corresponding test results with a significance level $\alpha=0.05$. At short distance, the three series involving G are almost all significant. Especially, all the G-E jumps are significant, which demonstrates a high consistency between our FGLS tests and the segmentation results. Almost all G-E' jumps are significant as well, while almost all E-E' are not significant. This latter result confirms our second rule (E and E' are expected to be consistent, i.e. either no jump or a jump in both, simultaneously). Most G-G' jumps are also
<figure><img src="image_35.png" /><figcaption>Figure 4.6 – Distribution of absolute jump amplitudes and their standard errors, and the associated t-values computed from the FGLS estimates of the real data (494 main/nearby pairs). The results are sorted based on the main/nearby distance (< 50 km and ≥ 50 km).</figcaption></figure>

significant, which confirms the idea that most jumps are in G. Finally, the G'-E' and G'-E jumps are most of the time insignificant, which again supports of our first rule (G and G' are unlikely to change simultaneously). As the distance increases, the proportion of insignificant jumps also increases due to higher standard errors.

<figure><img src="image_36.png" /><figcaption>Figure 4.7 – Distribution of test results associated to the estimated amplitudes of jumps shown in Figure 4.6, sorted by distance: (a) < 50 km, (b) ≥ 50 km. Test results are color coded as: green for insignificant and red/blue for significant downward/upward jump.</figcaption></figure>

### 4.1.5 Predictive rule

The objective is to build a classifier $\hat{\psi}(x)$ representing the prediction of the configuration y, i.e. the quadruplet composed of G, E, G', and E', given x, the vector composed of the test statistics from the series of differences. In
the development of this classifier, we are confronted with two principal challenges. First, the true configurations are unknown, resulting in the unavailability of y for training, evaluation and prediction. Second, the presence of all configurations in the data is nearly improbable. This is due to the rarity of occurrence for certain configurations according to the probabilities indicated in Table 4.1. To address these challenges, we propose to generate a synthetic dataset based on the $N = 494$ test results of the real data using a bootstrapping technique as introduced by (Efron and Tibshirani, 1993). This dataset would ensure each configuration is represented through a set of $(x, y)$ pairs. We then compare the performance of four popular classifiers for the prediction objective. The reader can refer to subsection 2.2.4 of Chapter 2 for more details about the supervised classification, two resampling methods and a presentation of the four algorithms we considered here.

### 4.1.5.1 Preliminary considerations

**Considered test results.** In this task, we employed the test statistics from five series: G-G', G-E', E-E', G'-E', and G'-E', given that the considered cases pertain exclusively to situations where the G-E test holds significance. Denote by $z_\ell = (z_{\ell1}, \dots, z_{\ell 5})$ the vector of the five test statistics for the $\ell$th test in the five series of difference and by $Z = (z_\ell)_{\ell=1,\dots,N}$ the formed data set of size $N \times 5$, with $N = 494$, which will be called the original data set in the sequel. Note that the results of all the nearby stations for a given change-point in a given main station can be viewed as replicates reducing the real information to 114, i.e. the number of couples of (main station, change-point); the available information is thus quite small.

**Considered configurations.** In Table 4.1, among the 38 configurations of the RSD Table, there are two doubles of the five coded test results with same prior probabilities: configurations (7,28) and (12,19). We decide to keep the configurations 7 and 19 which contains a change-point in G. This reduces the total number of configurations to $C = 36$.

**The four considered learning algorithms are:** the linear discriminant analysis (LDA) (Fisher, 1936), the classification and regression trees (CART) (Breiman et al., 1984), the Random Forest (RF) (Ho, 1995) and the k Nearest Neighbors (k-NN) (Fix and Hodges, 1989; Cover and Hart, 1967). The latter three involve parameters that need to be tuned. They are here automatically optimized by K-fold cross-validation with $K = 10$ using the generic function 'train' of the R package caret.

**Building of the complete synthetic data set.** As mentioned earlier, the bootstrapping technique has been employed to construct the synthetic dataset, operating on the principle of random sampling with replacement from the original data. More precisely, for each configuration $y$ and each series of difference $j$, we create $N_y$ vectors of the five test statistics or t-values (the sample $x$) by resampling among the test statistic values $(z_{\ell j})_\ell$ that lead to the test conclusion of $y$. The correspondence is made with respect to the test outcome ($-1, 0$ or $1$) for a given significance level $\alpha$. As an example, let's consider the configuration 1 in Table 4.1. Each test statistic (t-value) is randomly selected from the respective series of difference, ensuring that the significance levels of these five t-values are (1, 1, 0, 0, 0). The constructed data set is noted $D = \{(y_\ell, x_\ell)_\ell\}_{\ell=1,\dots,n}$ of size $n = \sum_y N_y$.

A notable consideration is the potential severe imbalance in the configurations within the data, as discussed in subsection 2.2.4.3 of Chapter 2. This feature is well-known to produce biased classifiers for the minor
configurations. Two options can be considered: if we consider the imbalanced aspect as a problem, we can use the same number of replicates $N_y = R$ for each configuration y. This case is called the 'balanced sample case'. On the contrary, if we want to take into account this reality, we can take a number of replicates proportional to the prior probability of each configuration given in Table 4.1, i.e. $N_y = np_y$ if n is the total number of data. This case is called the 'imbalanced sample case'.

### 4.1.5.2 The proposed Cross-Validation Bootstrap (CVB) procedure

Cross-validation is a popular statistical technique to test a classifier. It involves splitting the data into two subsets: the learning set, on which the classifier is constructed, and a test set, on which the classifier is tested. Since observations of the complete data set D are replicated from the original data set Z, which is small and repeated, there is a risk of overlap between the learning and test data set, inducing inevitably a bias and leading to an underestimation of misclassification error. This is why, we propose here a so-called cross-validation bootstrap (CVB) strategy which consists in first splitting the original data set Z into the learning and test subsets before constructing the complete data set D. The proposed CVB procedure is described in Algorithm 1.

**Data:** the original data Z

for b = 1 to B do

1. sample a learning data set $Z^{b,L}$ from Z with probability 0.8, and form the test data set $Z^{b,T}$ with the remaining 20% of data. The random sampling is performed on the rows of Z, i.e. on each test;

2. form the two associated complete data sets $D^{b,L}$ and $D^{b,T}$ from $Z^{b,L}$ and $Z^{b,T}$ by preserving the learn/test proportion of 80%/20%, i.e. for each configuration y, $D^{b,L}$ contains 0.8$N_y$ samples, and $D^{b,T}$ 0.2$N_y$. In the 'balanced sample case', we chose $N_y = R = 100$ and in the 'imbalanced sample case', the smallest value $N_y$ is chosen to 5, leading to a learn sample containing 4 data and a test sample containing only one data;

3. construct the four classifiers on the learning data set $D^{b,L}: \hat{\psi}^{b,k}, k \in \text{LDA, CART, RF, k-NN}$;

4. compute the misclassification error of the classifiers on the test data set $D^{b,T}$ with $n_T$ rows:

$$ \text{err}^{b,k} = \sum_{\ell=1}^{n_T} \mathbb{1}_{\{\hat{\psi}^{b,k}(x_{\ell}^{b,T}) \neq y_{\ell}^{b,T}\}} \quad \text{for } k \in \text{LDA, CART, RF, k-NN} $$

end

**Evaluation:** compute the mean and the standard deviation of misclassification error for each classifier, e.g. the mean is given by:

$$ \overline{\text{err}}^k = \frac{1}{B} \sum_{b=1}^{B} \text{err}^{b,k} \quad \text{for } k \in \text{LDA, CART, RF, k-NN} $$

**Algorithm 1:** The CVB procedure.

Table 4.4 gives the mean and the standard deviation of misclassification error for the four considered classifiers with $B = 20$. The table presents results for three scenarios: the first and second one involve constructing the
complete dataset using different sampling (balanced vs. imbalanced), both with $\alpha = 5\%$, while the last one employs a balanced sampling with $\alpha = 1\%$. Compared to the balanced sample, the misclassification error is lower for the imbalanced sample in the case of LDA and k-NN, but slightly higher for CART and RF. Similar behavior is observed when comparing learning with $\alpha = 5\%$ and $\alpha = 1\%$ for the balanced sample. Overall, the Random Forest algorithm outperforms the other classifiers in all three scenarios, with the best performance achieved when trained with a balanced sample with $\alpha = 5\%$. We thus choose the Random Forest algorithm and select as the final predictive rule, $\hat{\psi}$, the one with smallest error among the B classifiers. The predictive power of the five series of difference based on the accuracy criterion (the percentage of correct predictions) are in the decreasing order: E-E', G-G', G-E', G'-E and G'-E'.

<table><caption>Table 4.4 – Mean misclassification error ± one standard deviation, for the four classifiers in three scenarios: 'balanced sample' with $N_y = R = 100$ and $\alpha = 0.05$, 'balanced sample' and $\alpha = 0.01$, 'imbalanced sample' with $N_y = np_y$ and $\alpha = 0.05$.</caption><thead><tr><th>c</th><th>test level</th><th>sample case</th><th>LDA</th><th>CART</th><th>KNN</th><th>RF</th></tr></thead><tbody><tr><td rowspan="3">err<sup>c</sup></td><td>0.05</td><td>balanced</td><td>0.1463 ± 0.021</td><td>0.0142 ± 0.011</td><td>0.1412 ± 0.018</td><td>0.0049 ± 0.003</td></tr><tr><td>0.05</td><td>imbalanced</td><td>0.1108 ± 0.004</td><td>0.0165 ± 0.010</td><td>0.0351 ± 0.004</td><td>0.0054 ± 0.004</td></tr><tr><td>0.01</td><td>balanced</td><td>0.1424 ± 0.029</td><td>0.0210 ± 0.0417</td><td>0.1301 ± 0.022</td><td>0.0106 ± 0.033</td></tr></tbody></table>

### 4.1.5.3 Application to the real data set

The objective is to predict the configuration for each change-point of each main station. When several nearby stations are available for a given change-point the results are aggregated using a weighted prediction score. For a configuration c, the prediction score writes:

$$ \hat{P}(y_{(\text{main,change-point})}) = c|\text{nearby station}(ns)) = \frac{\sum_{ns} w_{ns} \mathbb{1}_{\{\hat{\psi}(x_{ns})=c\}}}{\sum_{ns} w_{ns}} $$

where $w_{ns}$ denotes the weight of the nearby station ns, and the final configuration is the one with the highest score

$$ \hat{y}_{(\text{main,change-point})} = \arg \max_{c} \hat{P}(y(\text{main, change-point}) = c | \text{nearby station}) $$

We compared two different weightings:

- inverse distance weighting: $w_{ns} = 1/d_{ns}$, where $d_{ns}$ is the distance between the nearby ns and the main station,

- weighting proportional to the joint probability given in Table 4.1: $w_{ns} = p_c$.

In the probability-based weighting, when the highest score is reached by two different configurations (e.g. $c = 1$ and 10), the one with the shortest distance is selected.

Figure 4.8 presents the distribution of predicted configurations, after aggregation, for four variants: (a) balanced sampling with $\alpha = 5\%$, aggregated with distance; (b) balanced sampling with $\alpha = 1\%$, aggregated with distance; (c) imbalanced sampling with $\alpha = 5\%$, aggregated with distance; and (d) balanced sampling with $\alpha = 5\%$, aggregated according to prior probability. Across all figures, four predominant groups emerge consistently: G
(c=1 and 15), (G, E, E') (c=31 and 35), (E, E') (c = 10 and 23), and E (c = 8 and 22). Remarkably, these configurations correspond to the highest joint probabilities, *p*, as indicated in Table 4.1: G and (E, E') with *p* = 0.18, (G, E, E') with *p* = 0.04, and E with *p* = 0.01. This demonstrates that the classifier actually predicts the configurations which we believe are the most likely in the real data, even when these probabilities are not directly used in the procedure such as in variants (a) and (b).

<figure><img src="image_37.png" /><figcaption>Figure 4.8 – Distribution of the final predicted configurations, after aggregation, from the real data with the Random Forest method. The numbers in color bars refer to the configuration number *c* among the 38 cases displayed in RSD Table (Table 4.1). Results are plotted for four cases: a) balanced sample learning with *α* = 0.05 and aggregated by distance, b) balanced sample with *α* = 0.01 and aggregated by distance, c) imbalanced sample with *α* = 0.05 and aggregated by distance, and d) balanced sample with *α* = 0.05 and aggregated by prior probability.</figcaption></figure>

In variant (a), 47 of the change-points (i.e. 41%) are attributed to group G and 29 (i.e. 25%) to group (G, E, E'), after the aggregation. Analysis of the six test results before and after the prediction helps to understand the relatively high frequency of these two groups. In general, the test results can be of two sorts: either the six results correspond to a configuration in the Table 4.1, and in this case the predictive rule predicts the same result (as expected), or the result is initially not in the Table, and the predictive rule will select a configuration that is "close" to the initial configuration. Among all the test results going to group G, i.e. (1,1,1,0,0,0) for *c* = 1 and (-1,-1,-1,0,0,0) for *c*=15, about 75% are initially in the Table. This high percentage is consistent with the observation that many jumps are significant in the first three tests and insignificant in the last three, as seen
in Figures 4.6 and 4.7. The 25% of cases which are not initially in the Table differ from these configurations by one or two elements, e.g. case (1,0,1,0,0,0) differs from $c=1$ by only the 2nd element (the G-G' test). This case is then attributed to $c=1$ by the predictive rule when the absolute value of the t-value of the estimated jump in the G-G' series is close to the critical value, $\tau_{\alpha/2}=1.96$, in combination with smaller t-values in E-E', G'-E', and G'-E. For group (G, E, E'), the percentage of cases that are not in the Table is slightly more than 50%. Almost all these cases are either (1,1,1,0,1,0) or (-1,-1,-1,0,-1,0), which differ only by the 6th element (the G'-E test) from the final configurations $c=31$ (1,1,1,0,1,1) or $c=35$ (-1,-1,-1,0,-1,-1), respectively. Contrary to the G-G' series, the G'-E series has smaller t-values on average, hence the frequent 0 in the initial test results. The fact that almost all these cases are finally predicted as $c=31$ or 35 can be explained by the simultaneous occurrence of: high t-values in G-G' and G-E', a small t-values in E-E', a high t-values in G'-E', and a value close to the critical value for G'-E. An example is provided in Figure 4.9. In this example, one may also suspect that the t-value of G'-E' is excessively large, given the small value of the corresponding jump (-0.12) compared to the jumps in G-E, G-G', and G-E'. One way to reduce the occurrence of excessively high t-values in G'-E' is to increase the critical value of the test. For example if we set $\tau_{\alpha/2}=2.58$ ($\alpha = 0.01$), the test result here becomes 0 and the initial configuration becomes $c = 15$, which has a higher probability in Table 4.1 and is thus preferred.

<figure><img src="image_38.png" /></figure>

Figure 4.9 – Example of test result for station FAIR (Fairbanks, Alaska) with nearby station CLGO at a distance of 21 km. The series of IWV differences are shown in gray. The black vertical solid line shows the change-point detected in G-E by the segmentation (04 October 2017). The blue triangles indicate known equipment changes in the main station from the GNSS metadata. The horizontal red lines show the means estimated by the FGLS regression on the left and the right of the change-point in each series.

The impact of using $\alpha = 0.01$ is further illustrated on all tests with Figure 4.8b. Only 9 change-points are now assigned to group (G, E, E'), which is considerably smaller than with $\alpha = 0.05$. Actually, 10 change-points
moved to group G and 9 to group (E, E'). This difference can be understood by inspecting the distribution of t-values with respect to the corresponding critical values (2.58 vs. 1.96). Figure 4.6 shows that many t-values for E-E', G'-E, and G'-E' are smaller than 2.58 in absolute value. When these tests become insignificant, while the other three stay significant, the predicted configuration becomes $c=1$ or 15, and when G'-E' remains significant or close to 2.58, the predicted configuration becomes $c=31$ or 35. Additionally, many configurations with low probabilities ($c=3, 7, 17, 18, 32, 33$) have also disappeared.

Figure 4.8c shows the impact of the (probability-based) imbalanced sampling in the learning procedure. Overall, the results for the main groups are not much different compared to the balanced sampling. Two noticeable differences emerge, however. Firstly, group (G, E, E') reduces only slightly in size, from 29 to 20. The smallness of the impact is explained by the fact almost half of test results are initially in the Table 4.1 and are not changed by the prediction. Secondly, almost all the configurations with the lowest probabilities, such as $c = 7, 21, 32, 33, 38$, with $p \le 5.6 \times 10^{-4}$, are removed. Other configurations, with slightly higher probabilities, but still with $p \le 0.01$, such as $c = 3, 16, 17$ (G, E') and $c = 29$ and 34 (G, E), emerge or are reinforced, which is not wanted.

Figure 4.8d shows the variant (d) where the aggregation is based on the prior probabilities. The distribution is quite different from that based on distance (Figure 4.8a): more change-points are attributed to the preferred groups, 62% in G and 19% in (E, E'), fewer to other groups such as (G, E, E') and E, and many configurations of low probability disappear. The distribution is actually quite similar to that of variant (b), but in contrast to the latter, this variant keeps a high power in the test (thanks to $\alpha=0.05$). As a result, with variant (d), group E is much smaller than with variant (b). In this respect, variant (d) is preferred among all four variants. Note, however, that there is a limitation in the usage of the aggregation procedure which holds for all variants: when there is only one nearby station (30 % of the cases), the aggregation has no impact and the final result is the one selected by the prediction rule. In variant (d), this explains why there are still configurations with a low probability ($c=8/22, 3, 31/35, 21, 38$).

### 4.1.6 Conclusions and perspectives

We proposed a post-processing method for the attribution of change-points detected by a segmentation scheme involving multiple series of differences (target minus reference). In our application, each of the stations provides a GNSS series (G) and a reanalysis series (E). The segmentation is run on the G-E series and the goal of the attribution method is to predict if the inhomogeneity (jump in the mean) is in G or E. The method proceeds along the following steps:

1. **Data selection and pre-processing.** For each detected change-point in a main station (hereafter, the "main change-point"): (a) select nearby stations with a horizontal distance smaller than 200 km and height difference smaller than 500 m; (b) run the segmentation method on the G'-E' for each nearby data and select only homogeneous segments from the nearby to compare with the main; (c) correct the nearby series, G' and E', for the height difference with respect to the main station, so that all four series (G, E, G', and E') are representative of the same height; (d) form the six series of differences (G-E, G-G', etc.) and remove the outliers.

2. **Test the significance of the jumps.** For each main change-point and each of the six series: (a) identify the noise model; (b) fit a regression model including a jump at the position of the main change-point
when at least $n = 200$ consecutive points are available on the left and right of the change-point, using an iterative FGLS procedure; (c) test the significance of the jump at the significance level $\alpha = 0.05$.

3. **Use a predictive rule to predict the configuration.** For each nearby, the learned classifier will predict the configuration (i.e. which of the G, E, G', and E' series have a significant jump) corresponding to the six test results. When several nearby series are used, a weighted prediction score is computed to select the final configuration.

The method has been applied to a real data set of 494 cases (114 change-points from 49 main stations compared with 312 nearby stations). The data characterization showed that the data have a strong heteroscedasticity, with mean annual seasonal variation in the standard deviation around 70% (half-range), and a moderate autocorrelation, with a typical lag-1 correlation coefficient of 0.3. A FGLS test procedure was implemented to ensure an accurate inference. The predictive rule has been trained on the real data. Several classifiers have been compared for the predictive rule and the Random Forest was selected.

To our knowledge, both the FGLS regression test approach and the Random Forest classifier have never been used in the context of climate series homogenization.

The FGLS tests and the classification results of the studied data set have been assessed using i) our expertise of the data set (formulated out in two probabilistic rules) and ii) metadata informing about known equipment changes at the GNSS sites. Very consistent and plausible results were found from both the FGLS tests and the classification. With a significance level of 5% and employing a balanced sample for the learning step in the predictive rule, as well as aggregating results from nearby sources based on prior probability, the findings clearly indicate predominance of significant jumps in the series involving G (62%), (E, E') (19%), and (G, E, E') (10%) as expected. The remaining 9% of unexpected results are thought be linked with low detection power of the FGLS test when the noise is large (e.g. due to large distance between the main station and the nearby) and possibly random errors in the classification due to the smallness of the learning sample.

Some possibilities to further improve the method are: i) to use a bigger data set to improve the predictive method, ii) to refine the nearby selection rules to improve the robustness and the power of the test procedure (e.g., select nearby series with smaller percentage of gaps, shorten the distance between the main station and nearby), iii) compute the critical value used in the FGLS test from a more realistic empirical distribution. These options will be tested in a future work.

### 4.1.7 Appendix

#### 4.1.7.1 Test table

Table 4.1 shows the theoretical test results one would obtain for the six paired difference series (G-E, G'-E', G-G', G-E', G'-E, E-E') with a perfect test method, for 54 different combinations of jumps in the four base series (G, E, G', E'). The jumps in the base series are coded on three values: 0 (no jump), +1 (upward jump), and -1 (downward jump). Here, only cases where either G or E, or both, have a jump are considered, because we always start we a change-point detected by the segmentation method in the G-E series (so either in G, in E, or in both). The corresponding results in the SD Table are thus theoretically coded on five different values (0,
1, 2, -1, and -2). However, in practice, a test result will be either reject (-1 or +1, where the sign indicates if the jump is upward or downward) or fail to reject (0) the null hypothesis of no jump. The RSD Table shows the corresponding practical results where the results are coded on three values only. The SD Table contains actually 46 unique combinations of the six test results and 8 duplicates highlighted by a colored background, while the RSD Table contains only 38 unique combinations. The duplicates are sorted out depending on probabilities (see Appendix B), i.e. those with lower probabilities are not counted.

### 4.1.7.2 Computation of prior probabilities

This section explains how we compute the conditional and joint probabilities reported in the Table 4.1. The notation hereafter is *A* (italics) for events and *A* (straight) for time series.

Let *G*, *E*, *G'*, and *E'* represent the event "there is a jump" in each of the four base series, with the possible outcomes: 0=no jump, -1=the jump is downward, +1=the jump is upward. Let $P(A|B)$ represent the conditional probability of *A* given *B*. Let us assume that the events *G* and *E* are independent, and that *G'* and *E'* are independent. Now we reformulate the two rules that we stated in the Introduction and associate them with probabilities:

$$ (R1) \quad P(G' \neq 0 | G \neq 0 \cup E \neq 0) = 0.1. $$

$$ (R2) \quad P(E = E') = 0.9. $$

The values of 0.1 and 0.9 are chosen in a way to reflect the contrast between the corresponding events. Note that the first rule is stated more generally than in the Introduction as here it is now also conditional on a jump in series *E*. Indeed, in practice, it is unlikely that a jump in series *G'* occurs at the same time as in any other series. The first probability accounts thus for two incompatible events: ($G' \neq 0 | G \neq 0, E = 0$) and ($G' \neq 0 | G = 0, E \neq 0$), which we can suppose to be of similar probability $p_1 = 0.05$. It follows from the first rule that $P(G' = 0 | G \neq 0 \cup E \neq 0) = 0.9$. It follows from the second rule that $P(E' \neq E) = 0.1$, which also corresponds to two incompatible events ($E' = 0 | E = -1$) and ($E' = 0 | E = +1$), the probability of which we will also assume to be equal to $p_2 = 0.05$.

The conditional probabilities $P(G', E' | G, E)$ are obtained from the products of the individual probabilities $P(G', E' | G, E) = P(G' | G, E) \times P(E' | G, E)$. These probabilities are sufficient to distinguish the duplicates in the first 36 rows of the SD Table and the RSD Table (e.g. rows 1 and 18, in the numbering going from 1 to 54, have the same 6 values but are associated with different conditional probabilities). However, more duplicates appear in the RSD Table in the last 18 rows (row numbers 37 to 54). In order to distinguish them, we introduce prior probabilities associated to the events *G* and *E*. There are six different combinations: ($G, E) \in \{(+1, 0), (0, -1), (-1, 0), (0, -1), (+1, -1), (-1, +1)\}$. In the first 4 combinations, a single change occurs in either *G* or *E*, while in the latter 2, two changes occur simultaneously. The former is more likely, so we attribute it a probability of $p_3 = 0.225$, from which it results that the latter has a probability of $p_4 = 0.05$ which is deduced from the equation $4 \times p_3 + 2 \times p_4 = 1$.

Finally, the joint probabilities reported in Table 4.1 are obtained from the product of the conditional and the prior probabilities: $P(G, E, G', E') = P(G', E' | G, E) \times P(G, E)$.
## 4.2 Additional studies

### 4.2.1 Assessment of the stochastic model identification and parameter estimation

In the work presented in the previous section, we used the `auto.arima` function of the R package `forecast` (Hyndman and Khandakar, 2008) to identify the noise model in the time series of IWV differences and the `arima` function to estimate the model parameters. Results are presented in Figure 4.3 of the paper, for all the change-points tested from the main-nearby pairs, for each of the six series of differences.

The utilization of `auto.arima` in time series forecasting is widely acknowledged in the literature (Shaub, 2020; Yan et al., 2022; Radke et al., 2020; Lai and Dzombak, 2020). However, its skill for model identification has not been extensively discussed. In this subsection, we specifically evaluate its effectiveness within our application framework, i.e. the identification of four particular stochastic models: white noise (WN), AR(1), MA(1), and ARMA(1,1). We recall that, in a preliminary investigation, we tested ARMA($p, q$) models up to order $p, q = 2$, but we found that only a negligible proportion of the time series were identified with order 2 and few of them had actually significant coefficients ($\alpha = 0.05$), so we limited ourselves to order $p = 1$ and $q = 1$.

Our examination below investigates the impact of the length of the series, $N$, and the coefficient values, $\phi$ and $\theta$, on the variance of the coefficient estimates. We have seen from equations (2.20 to 2.23), introduced in subsection 2.2.3.8, that the variance of the coefficient estimates increases when the coefficients values decrease. To which extent the coefficient values impact the model identification performance is, however, not well known. The subsections below first assess the model identification performance of `auto.arima` in terms of the probability of identifying the true model. Then, the uncertainty of the coefficients estimated with the `arima` function is examined assuming the model identified is correct. We also investigate the influence of gaps in the time series on both model identification and parameter estimation quality. This is a situation frequently encountered with real data. Finally, we discuss the possible impact of the uncertainty in the model identification and parameter estimation in the case of the real data presented in the paper.

#### 4.2.1.1 Simulation set-up

The simulations involve varying the length of time series, ranging from $N = 200$ to 2000 data points. Coefficients for the AR(1) and MA(1) models are chosen from the range of 0.1 to 0.6, covering typical values encountered in real data. For the ARMA(1,1) model, the study examines six cases where the sum $s = \phi + \theta$ goes from 0.1 to 0.6, where $\phi$ is varied from 0.8 to -0.8 and $\theta$ is adjusted accordingly. Importantly, cases where $\phi = 0$ or $\theta = 0$, representing a MA(1) or an AR(1) model, respectively, are excluded from the ARMA(1,1) discussion. The analysis comprises 1000 simulation runs and the variance of the white noise (or innovation) is set to 1.

#### 4.2.1.2 Noise model identification results

The relationship between series length and the efficiency of model identification is presented in Figure 4.10 and Figure 4.11. All the plots demonstrate a consistent trend: the True Positive Rate (TPR), which quantifies the fraction of time series for which the true model is correctly identified, improves as the series length increases. However, the slope of the TPR curves varies depending on the model type and coefficient values. Three interesting features can be mentioned. Firstly, the dependency of TPR on series length for the AR(1)
and MA(1) models appear to be strongly resembling for similar coefficient values. The slope of the TPR curves is gentler for larger coefficient values and becomes steeper as the coefficient magnitude reduces. This behavior is certainly in connection with the impact of coefficient values on estimation uncertainty (smaller coefficients leading to larger uncertainty, discussed in Chapter 2.2.3). Notably, for coefficient values of 0.3 or larger, the TPR exceeds 70% when $N = 200$ and 95% when $N = 1400$, for both models. However, with smaller coefficients, a TPR value of 95% is never reached, even for a series length of up to 2000. Secondly, for the ARMA(1,1) model, the shape of the TPR curve is clearly dependent on the sum of the coefficients, $s = \phi + \theta$, where inverting the values of $\phi$ and $\theta$ yields similar TPR results. It is also seen that the slope of the TPR curves if steeper when either the absolute values of coefficients are larger or when the sum is larger. For example, for a scenario wherein the sum of the two coefficients is $s = 0.3$, the TPR goes from roughly 40% for $N = 200$ to 95% for $N = 800$ when $\phi = 0.7$ and $\theta = -0.4$ (and similarly when $\phi = -0.4$ and $\theta = 0.7$). However, with the combination $\phi = 0.6$ and $\theta = -0.3$, the TPR of 95% is only achieved for a longer series $N = 1800$. Lastly, several cases exist where the TPR remains notably low even with extended series lengths, suggesting inherent challenges in model identification. For example, when the sum of the coefficient values is 0.3, coefficient pairs such as (0.1, 0.2), (0.2, 0.1), (-0.1, 0.4), and (0.4, -0.1) result in TPRs below 10%, even when $N$ reaches 2000.

In Figures 4.12 and 4.13, the same results are plotted in an alternative way to enhance the influence of model coefficients on the efficiency of the model identification.

For the AR(1) and MA(1) models shown in Figure 4.12, the TPR ascends as coefficients increase, reaching a value near 1 for $\phi$ or $\theta$ values of approximately 0.6 for all series lengths. For smaller coefficients, the TPR displays greater dispersion across series lengths. It is interesting to examine which alternative model is actually identified when TPR is smaller than 100%. Figure 4.14 shows that a negligible fraction of cases are identified as ARMA(1,1) when the true series is AR(1) or MA(1). This can be easily understood by the parsimony effect embedded in the BIC model selection. However, when the coefficient of the AR(1) model becomes small, a non negligible fraction of cases are identified as MA(1) and vice-versa. In the extreme case when the coefficients are close to zero, they may be identified as WN (e.g., about 30% when $\phi$ or $\theta=0.1$).

For the ARMA(1,1) model (Figure 4.13), the TPR behavior is very different depending on whether the sum of $\phi$ and $\theta$ is smaller or larger than 0.3, showing either single or dual dips in the TPR curve respectively. It is noteworthy that these dips consistently appear when one of the coefficients nears zero. For example, with $s = 0.3$, when either $\phi = 0.1$ and $\theta = 0.2$ or its converse, the TPR goes to zero. Conversely, for $s > 0.3$, e.g. when $s = 0.5$, dual dips emerge for combinations $\phi = 0.4$ and $\theta = 0.1$, and $\phi = 0.1$ and $\theta = 0.4$. The reason why the TPR is so low in this case can be understood by examining the details of the model identification outcomes shown in Figure 4.15. Indeed, in situations where one parameter is small and one is large, it is logical that a simpler model is preferred, i.e. a AR(1) or a MA(1). This is clearly observed on this Figure: an ARMA(1,1) is correctly identified only when both $\phi$ and $\theta$ are not close to 0, whereas it is predominantly identified as AR(1) when $\theta$ is small and as MA(1) when $\phi$ is small. Moreover, the fraction of correct identification not only depends on the values of $\phi$ and $\theta$ but also on their sum, $s = \phi + \theta$. Notably, Figure 4.15e shows that for the latter example with $s = 0.5$, over 85% of the series are classified under the AR(1) model (note that this simulation is specific to the case $N = 1000$). However, a special situation is observed in Figure 4.15a for $s = 0.1$, where a
significant fraction of cases are identified as WN even when $\phi$ and $\theta$ get large. Indeed, when $s$ is close to zero, the values of $\phi$ and $\theta$ are nearly opposite and the stochastic process is actually close to a white noise process. Identifying an ARMA model when $\phi = -\theta$ is an ill-posed problem which leads to infinite variance in both estimated parameters as reflected from Eqs. (2.22 and 2.23). This issue is referred to as parameter redundancy or over-parameterization (Shumway and Stoffer, 2011).

<figure><img src="image_39.png" /><figcaption>Figure 4.10 – True Positive Rate (TPR) against series length for AR(1) model (left) and MA(1) model (right).</figcaption></figure>
<figure><img src="image_40.png" /><figcaption>(a) s = 0.1</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>(b) s = 0.2</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>(c) s = 0.3</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>(d) s = 0.4</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>(e) s = 0.5</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>(f) s = 0.6</figcaption></figure>

Figure 4.11 – True Positive Rate (TPR) against series length for ARMA(1,1) model. Plots (a) to (g) depict the results for six variations, with coefficient sums $s = \phi + \theta$ ranging from 0.1 to 0.6.
<figure><img src="image_46.png" /><figcaption>Figure 4.12 – True Positive Rate (TPR) as a function of the coefficient value for the AR(1) model (left) and the MA(1) model (right).</figcaption></figure>
<figure><img src="image_47.png" /><figcaption>Figure 4.13 – True Positive Rate (TPR) as a function of the coefficient values for ARMA(1,1) model. Plots (a) to (f) depict the results for six variations, with coefficient sums $s = \phi + \theta$ ranging from 0.1 to 0.6.</figcaption></figure>
<figure><img src="image_48.png" /><figcaption>Figure 4.14 – Outcomes from the model identification for AR(1) (left) and MA(1) (right) as a function of coefficient values, specifically for series lengths of 1000.</figcaption></figure>
<figure><img src="image_49.png" /><figcaption>(a) s = 0.1</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>(b) s = 0.2</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>(c) s = 0.3</figcaption></figure>

<figure><img src="image_52.png" /></figure>
<figure><img src="image_53.png" /><figcaption>(e) s = 0.5</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>(f) s = 0.6</figcaption></figure>

Figure 4.15 – Outcomes from the model identification of ARMA(1,1) as a function of coefficient values, specifically for series lengths of 1000. Plots (a) to (f) depict the results for six variations, with coefficient sums $s = \phi + \theta$ ranging from 0.1 to 0.6.

### 4.2.1.3 Parameter estimation results

Figures 4.16 and 4.17 display the mean standard error of the estimated coefficients for AR(1), MA(1), and ARMA(1,1) models obtained from the simulations. These values are in good agreement with the estimates computed from equations (2.20) to (2.23) presented in subsection 2.2.3. From these plots, it is convenient to predict the uncertainty in the estimated parameters when designing an experiment or analysis with real data. For example, for an assumed ARMA(1,1) process with $\phi = 0.6$ and $\theta = -0.3$, as reported in the paper, a standard error close to 0.21 is expected on both $\hat{\phi}$ and $\hat{\theta}$ when $N = 200$, which can be decreased to 0.08 for both estimates when $N = 1000$.
<figure><img src="image_55.png" /></figure>

Figure 4.16 – Mean coefficient standard error as a function of time series length for models: AR(1) with $\phi = 0.3$ on the left, MA(1) with $\theta = 0.3$ on the right.

<figure><img src="image_56.png" /></figure>

Figure 4.17 – Mean coefficient standard error as a function of time series length for the ARMA(1,1) model with $s = \phi + \theta = 0.3$.

<figure><img src="image_57.png" /></figure>

Figure 4.18 – Mean coefficient standard error as a function of time series length for the ARMA(1,1) model with $s = \phi + \theta = 0.5$.

### 4.2.1.4 Impact of gaps in the time series

In practice, observational data series frequently exhibit measurement gaps. While some studies have investigated the impact of gaps on stochastic model parameter estimations, and have proposed modified
estimation procedures assuming a known model, the impact of gaps on model identification has not been much studied. Intuitively, one can expect that missing data will weaken the correlation structure, making it more challenging to correctly identify the model. The simulations below investigate the performance of the `auto.arima` function in the presence of gaps for the special case of an AR(1) process with a value of $\phi = 0.3$. In this simulation setup, the length of the time series remained constant, while a fraction of data points were replaced with NA values.

Figure 4.19 shows the percentage of identified models when the true model is an AR(1) (left) and the corresponding coefficient estimates (right) when the fraction of missing data varies between 0% and 50%.

The model identification results show that the AR(1) process is generally well identified, with a small fraction being MA(1). The TPR of the AR(1) drops from 91% to 76% when the fraction of gaps rises from 0 to 50%. Interestingly, even with half the data missing, the ability to detect a serial correlation in the series (either as AR(1), MA(1) or ARMA(1,1)) remains high at around 99%, while the rate of misclassifying AR(1) as white noise is under 1% (when $\phi = 0.3$).

Figure 4.19b presents the coefficient estimates obtained with the `arima` function, which takes missing data (in the form of NA values) into account. The plot displays the distribution of $\hat{\phi}$ and $\hat{\theta}$ when either an AR(1) or a MA(1) model is identified while the data are simulated from an AR(1) with $\phi=0.3$. Firstly, for the correctly identified AR(1) cases, the estimates average consistently around 0.3; however, the scatter increases with increasing data gaps. Secondly, in cases where MA(1) is identified, the median of coefficients is about 0.33. The slightly larger coefficient value of the MA(1) process compared to the AR(1) process is easily explained by the fact that these processes have a similar value of the autocorrelation function (ACF) at lag 1 for these specific coefficient values.

<figure><img src="image_58.png" /><figcaption>Figure 4.19 – Percentage of identified models and the associated coefficient estimate against the percentage of missing data when simulating the AR(1) model with $\phi = 0.3$.</figcaption></figure>

### 4.2.1.5 Interpretation of results from real data

The simulation results presented in the previous subsection can help to assess the reliability of the data characterization and to interpret the results presented in the paper.
First, we examine the uncertainty in the identification of the three main stochastic models with their typical (average) coefficients found in the paper. Recall that the minimum series length considered in the paper is $N=400$, composed of 200 points on each side of the change-point, and in practice, it often exceeds $N=1000$. Beginning with the AR(1) model with $\phi = 0.3$, we expect a correct identification of this model with a TPR=80% for a sample size of $N = 400$, and a TPR exceeding 90% when $N > 1000$ (green line in Figure 4.10a). For the MA(1) model with $\theta = 0.2$, TPR is about 60% for $N=400$ and rises to 83% for $N = 2000$ (olive line in Figure 4.10b). For the ARMA(1,1) model, we consider two coefficient combinations: (a) $\phi = 0.6$ and $\theta = -0.3$, and (b) $\phi = 0.3$ and $\theta = 0.2$. For the former, we expect a TPR of 40% at $N=400$ which rises to 95% by $N=1800$ (purple line in Figure 4.11c). For the latter, TPR initially lingers 10% when $N = 400$ but approaches 95% only at $N = 2000$ (cyan-green line in Figure 4.11e). The standard errors for these two ARMA(1,1) cases are $std[\hat{\phi}] = 0.069$ and $std[\hat{\theta}] = 0.082$ in case (a), and $std[\hat{\phi}] \sim std[\hat{\theta}] \sim 0.065$ in case (b), according to Equations (2.22) and (2.23).

<figure><img src="image_59.png" /><figcaption>Figure 4.20 – Coefficient estimates, $\hat{\theta}$ vs. $\hat{\phi}$, for all six series of IWV differences and all (494) main-nearby pairs in the real dataset, with series length $N \ge 1000$. The colors correspond to different models (WN in blue, AR(1) in orange, MA(1) in yellow, and ARMA(1,1) in purple).</figcaption></figure>

Next, we examine the dispersion of the estimated coefficients in the real data. Figure 4.20 shows the coefficient estimates, $\hat{\phi}$ and $\hat{\theta}$, from the data analysed in the paper, restricted to series with $N>1000$ data. This plot complements Figure 4.3 of the paper, but without separating the results for the data combinations. The
representation of the results in the $(\hat{\theta}, \hat{\phi})$ plane gives a special insight into the characteristics of the identified ARMA(1,1) models.

The distribution of the $\hat{\phi}$ and $\hat{\theta}$ coefficients for MA(1) and AR(1) models show values extending over a quite large range, from about 0.05 (for both coefficients) to 0.45 for $\hat{\theta}$ and 0.9 for $\hat{\phi}$. The TPR for the coefficient values below 0.1 is smaller than 40% (green curves in Figure 4.12a and b), i.e. more than 60% of these time series could be mis-identified as white noise (30%) or MA(1) instead of AR(1) and vice-versa (30%) as revealed by Figure 4.14a and b. For coefficient values above 0.3, the TPR exceeds 90%, so these case are correctly identified and estimated. The uncertainty in the coefficient estimates (~ 0.03, Figure 4.16) certainly adds to the scatter in the observed range of values.

As for the identified ARMA(1,1) models, the estimated coefficients exhibit a very large scatter, covering almost all possible values, from -1 to +1, except around zero. The $(\hat{\phi}, \hat{\theta})$ points appear actually to fall into two major and one minor groups, depending primarily on $\hat{s} = \hat{\phi} + \hat{\theta}$. the first major group aligns more or less on the $\hat{s} = 0.2$ line and includes $\hat{\phi}$ values in the range 0.4 to 0.99. Referring to Figure 4.15b for the case $s = 0.2$, the values of $\hat{\phi}$ at the lower end have very low probability to be identified as ARMA(1,1), which means that this class may be under-represented in our analysis. At the higher end, the values align more on the $\hat{s} < 0.1$ line, and these series, although identified as ARMA(1,1), are in reality close to a white noise. The standard error of the coefficients for the latter cases becomes very large, with makes their interpretation difficult. The few cases in the minor group where $\hat{\phi}$ takes negative values and $\hat{\theta}$ positive values are close to $\hat{s} \sim 0.1$, and may just be some examples of biased estimates from the first major group.

The second major group of the ARMA(1,1) results is centred on $(\hat{\phi}, \hat{\theta}) = (0.3, 0.2)$. As was shown in Figure 4.3 of the paper, these cases correspond to all the E-E' series plus some other series when the distance between main and nearby is larger than 50 km. The scatter in this group is relatively small, which can be explained by the higher accuracy of the $\hat{\phi}$ and $\hat{\theta}$ estimates in this group compared to the previous (see 4.18). However, the probability of correct identification for this group remains only around 60% (Figure 4.15e), with almost 40% going to the AR(1) group. It is thus probable that the ARMA(1,1) cases from both groups are under-represented in our results, and that the observed AR(1) and MA(1) groups include the latter.

### 4.2.2 Understanding of test and prediction results

In this subsection, we complement the discussion on the test and prediction results given in Section 4.1.4.3 of the paper. First, our attention is on the results of the balanced sample learning (variant (a) of Figure 4.26. We will focus on cases where combinations of the six tests initially not in the Test Table (Table 4.1 presented in the paper). Second, we discuss the impact of the distance on the test results. Next, we analyse the distribution of the predicted configurations before aggregation and, last, we discuss the uncertainty in the final classifier.

**Further discussion of the results of variant (a).** Among the 494 tested cases, 225 correspond to configurations that are not in the Test Table. These cases can be assembled in different groups, of which the 3 main are presented in Table 4.5. Group 1 contains 52 cases, of which 32 correspond to test results of (-1,-1,-1,0,-1,0) and 20 of (1,1,1,0,1,0). Group 2 contains 37 cases, among which 22 correspond to test results of (-1,0,0,0,0,0)
and 15 of (1,0,0,0,0). Group 3 contains 14, where 10 cases correspond to test results of (-1,0,-1,0,0) and 4 of (1,0,1,0,0).

For the cases in group 1, the initial three tests have usually large t-values, while the subsequent three tests present lower, frequently insignificant t-values (which is a characteristic of configuration $c=1$). However, the cases in group 1 have a significant G'-E' test (this configuration is not in Table 4.5). This is due to the fact that the G' and E' series are co-located, resulting in a lower noise level compared to the other two series (E-E' and G'-E), which remain influenced by the distance (larger noise at larger distance). The reason why these cases are predicted into different configuration is discussed below.

The groups 2 and 3 are most frequently observed at greater distances, averaging around 114 km and 126 km, respectively. Group 2 becomes more prominent when the t-values in G-E' are insignificant, while Group 3 corresponds to cases when the t-values in G-E' are only slightly larger than the critical value (the average values of these t-values is 2.24 for the upward shift and -2.33 for the downward shift).

Examining the prediction outcomes for these three groups reveals two points. First, a general observation is that predictions typically adjust only one or two elements among the six. Second, the way the six tests relate to each other and to the critical value, determines which elements are adjusted. These underlying rationales is evident across all groups.

In the first group, with test results such as (-1,-1,-1,0,-1,0), configurations are predicted such as $c=15$ (where G'-E' goes from -1 to 0), $c=16$ (E-E' goes from 0 to -1), and $c=35$ (G'-E goes from 0 to -1). The clear predominance of prediction configuration $c = 35$ is further discussed in the subsection 4.1.5.3 of the paper. Analyzing further, the configuration predicted as $c=15$ has the following t-values for the 6 series: (-7.5, -4.9, -7.5, -1.4, -2, -0.6). Here, the G'-E' is close to the critical value (1.96) but considerably lower than the first two tests (G-G' and G-E'), while G'-E is far to the critical value. Consequently, it is intuitive to shift G'-E' to being insignificant and predicting $c = 15$ instead of 35.

In group 2, five test results are insignificant. The classifier adjusts either the G-E' and G-G' result or the E-E' and G'-E' result. This leads to the predicted configurations 1, 15, 8, or 22, mainly because they are the only configurations with three zeros in the Test Table (Table 4.1).

Finally, the classifier predicts the cases of the group 3 to configurations 1 and 15, which differ from the test results only by the G-G' result. The average absolute t-values for G-G' in this group stands at 1.6, nearing the critical value (1.96) more than the G-E', E-E', and G'-E' results which have average absolute t-values of 0.4, 0.7, and 0.8, respectively.
<table><caption>Table 4.5 – (left part) Three main groups of test results which are not in the Table 4.1; (right part) configurations predicted by the classifier and number of cases (rightmost column).</caption><thead><tr><th rowspan="2">Group</th><th colspan="7">Test results</th><th colspan="7">Predicted</th><th rowspan="2">Count</th></tr><tr><th>G-E</th><th>G-G'</th><th>G-E'</th><th>E-E'</th><th>G'-E'</th><th>G'-E'</th><th>Config</th><th>G-E</th><th>G-G'</th><th>G-E'</th><th>E-E'</th><th>G'-E'</th><th>G'-E'</th></tr></thead><tbody><tr><td rowspan="3">1</td><td>-1</td><td>-1</td><td>-1</td><td>0</td><td>-1</td><td>0</td><td>15</td><td>-1</td><td>-1</td><td>-1</td><td>0</td><td>0</td><td>0</td><td>1</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>16</td><td>-1</td><td>-1</td><td>-1</td><td>-1</td><td>-1</td><td>0</td><td>1</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>35</td><td>-1</td><td>-1</td><td>-1</td><td>0</td><td>-1</td><td>-1</td><td>30</td></tr><tr><td></td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>2</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>31</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>18</td></tr><tr><td rowspan="3">2</td><td>-1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>15</td><td>-1</td><td>-1</td><td>-1</td><td>0</td><td>0</td><td>0</td><td>17</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>22</td><td>-1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>-1</td><td>5</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>2</td></tr><tr><td></td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>8</td><td>1</td><td>0</td><td>0</td><td>-1</td><td>0</td><td>1</td><td>13</td></tr><tr><td>3</td><td>-1</td><td>0</td><td>-1</td><td>0</td><td>0</td><td>0</td><td>15</td><td>-1</td><td>-1</td><td>-1</td><td>0</td><td>0</td><td>0</td><td>10</td></tr><tr><td></td><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>4</td></tr></tbody></table>

<figure><img src="image_60.png" /><figcaption>Figure 4.21 – Example of test results for station POTS (Postdam, Germany) with nearby station LDB at a distance of 74 km. The series of IWV differences are shown in gray. The black vertical solid line shows the change-point detected in G-E by the segmentation (19 June 2009).</figcaption></figure>

**Impact of distance on test results.** In the subsection 4.1.4.3 of the paper, we have observed that the distance between the main and nearby stations has an impact on the standard error of the jump, and
therefore on the t-value. In particular, the standard error of estimated jumps is clearly smaller for short distance (lower than 50 km). Here, we discuss in more details this point according to the noise level of the series.

First, recall that the variance of the GLS jump estimator, $\hat{\delta}_{GLS}$, is given by Equation (4.6): $var[\hat{\beta}_{GLS}] = (X'\Sigma_0^{-1}X)^{-1}$, where $\Sigma_0$ is the variance-covariance matrix of the noise. In the simplest scenario of independent and identically distributed (IID) noise with variance $\sigma_0^2$, we have $\Sigma_0 = \sigma_0^2 I_n$, and it can be shown that, for a simple model with no coefficients for the Fourier series, when the change-point is located in the middle of a series of total length $n$, the variance of the estimator of the jump then writes:

$$var[\hat{\delta}_{GLS}] = \frac{4}{n^2}\sigma_0^2 \qquad (4.11)$$

In the case when the noise is autocorrelated but remains homoscedastic, then $\Sigma_0 = VCV = C\sigma_0^2$, and $var[\hat{\delta}_{GLS}]$ is still proportional to $\sigma_0^2$. Finally, when the noise is autocorrelated and heteroscedastic, it can be shown that $var[\hat{\delta}_{GLS}]$ is proportional to the time-mean variance $<\sigma_t^2>$. In the real data, this variance increases with the distance.

Figure 4.22 presents the graph of the estimated standard error of the jump estimator as a function of the time mean of the moving standard deviation for each series of difference. The relationship between the two variables is almost linear, as expected. The (positive) correlation between these two variables ranges from $r = 0.76$ to 0.84. For non-collocated series, the range of values is larger, where the larger values correspond to larger distances (distinguished by red and blue points). The collocated series, G-E and G'-E', have a maximum standard deviation of about 1.6, whereas for non-collocated series, it lies between 3 and 3.5.

Figure 4.23 plots the noise level of the four non-collocated series as a function of the distance. As expected, the further away the nearby station is, more noise is introduced into the difference series. Consequently, a strong dependence of the uncertainty in the jump estimates on the distance becomes evident as depicted in Figure 4.24.

Figure 4.25 presents the t-values of the four non-collocated series plotted against distance. Among these, G-G' and G-E' display a declining trend in t-values as distance increases. In contrast, G'-E and E-E' show a near uniform t-value distribution. This pattern is consistent with the findings observed in the t-values box plot presented in the paper. The decreasing trend in G-G' and G-E' is attributed to the rise in standard error with increased distance, while the uniform distribution in G'-E and E-E' is due to their relatively smaller jump amplitudes compared to the others.
<figure><img src="image_61.png" /><figcaption>Figure 4.22 – Estimated standard error of the jump estimator as a function of the moving standard deviation for the four non-collocated series. The value *r* is the correlation coefficient between the two variables and the points are colored according to whether the distance between the series is less (red) or more (blue) than 50 km.</figcaption></figure>
<figure><img src="image_62.png" /><figcaption>Figure 4.23 – Mean of the moving variance as a function of the distance in four non-collocated series.</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>Figure 4.24 – Estimated standard error of the jump estimator as a function of the distance in four non-collocated series.</figcaption></figure>
<figure><img src="image_64.png" /><figcaption>Figure 4.25 – T-value as a function of the distance in four non-collocated series.</figcaption></figure>

**Predictive rule.** Further insight into the attribution results presented in the subsection 4.1.5.3 of the paper can be gained by looking into the distribution of predicted configurations of the whole 494 cases before aggregation. Here we analyse three of the variants presented in the paper: (a) balanced and (c) imbalanced sampling with $\alpha = 0.05$, and (b) balanced sampling with $\alpha = 0.01$. The results are presented in Figure 4.26. For each configuration, we distinguish whether the test results from 6 series of differences corresponded to a configuration in the table or not.

At a significance level of $\alpha = 0.05$, 269 cases correspond to configurations in Table 4.1 (yellow bars in Figure 4.26a), predominantly in configurations $c=15, 1, 31, 35, 23$, and 10, which hold the highest conditional probabilities in the Test Table (Table 4.1). The green bars in the figure highlight cases that are not in the table and are predicted, for variant (a), as $c = 35$ (41 cases), $c = 15$ (37 cases), $c = 23$ (27 cases), $c = 8$ (21 cases), $c = 10$ (19 cases), $c = 31$ (18 cases), etc. The dominance of those configurations was previously explained through the three main groups in Table 4.5.

Modifying the significance level from 0.05 to 0.01 has a pronounced effect on the prediction results. We observe an increase in the cases predicted under configurations $c = 1, 15, 8$ and 22, while predictions for configurations $c = 31$ and $c = 35$ show a decline. This is consistent to changes observed in the final results in paragraph 4.1.5.3.

The final variant (c) is obtained by training with an imbalanced sample. To discern the influence of this imbalance consideration, we contrasted them with results from balanced sample. The results of both are given
in Figure 4.26c and 4.26a respectively. While the count of the test results corresponding to a configuration, represented by orange bars, remains unchanged (since the significance level is consistent), the effects of sample balance manifest in three distinct points showed by green bars. Firstly, the G group ($c = 1, 15$) and E+E' group ($c = 10, 23$) maintain their values with only slight change. Secondly, configurations with simultaneous change-points at G and E possessing small marginal probabilities (e.g., configurations 31, 33 and 37) are absent in predictions, as expected by using imbalanced sample (see the discussion in the subsection 2.2.4.3). Thirdly, specific configurations, such as $c = 3$ and $c = 16$ (with test configurations (1,1,1,0,1,1) for 3 and (-1,-1,-1,-1,-1,0) for 16), show unexpectedly high prediction rates. These predictions arise from test combinations (-1,-1,-1,0,-1,0) and (1,1,1,0,1,0), which were predicted to configurations 35 and 31 in balanced sample training respectively (with test configurations (-1,-1,-1,0,-1,-1) for 35 and (1,1,1,0,1,1) for 31). The predictive rule of this variant modifies the E-E' from 0 to $\pm 1$ diverging from the variant (a) that adjusted the G'-E from 0 to $\pm 1$. This phenomenon arises due to the higher probability of the $c = 31, 35$ configurations.
<figure><img src="image_65.png" /><figcaption>Figure 4.26 – Distribution of the predicted configuration for each of the 496 cases. Colors distinguish between test results aligning with table configurations (Yes in yellow) and those that do not (No in green).</figcaption></figure>

**Discussion on the choice of the final classifier.** The training process is applied $B = 20$ times (see subsection 4.1.5 of the paper) and the final classifier is chosen as the one with the best predictive power among the 20 classifiers. It is noteworthy that our dataset's limited size, with a lot of sparsity in terms of configurations, can lead to strong replication during the formation of the completed dataset, potentially causing overfitting issues
in the construction of the classifier. Consequently, 2 different classifiers among the 20 can both show similar good performance in terms of mean misclassification error, but they may predict exactly the same results. The question here is to assess to what extent choosing one classifier instead of another impacts the prediction results.

Figure 4.27 illustrates the outcomes from two distinct Random Forest classifiers (for the repetitions $b = 3$ and $b = 4$) on the real dataset, both achieving an identical mean error of 0. These classifiers were trained with an equal sample size of $R = 100$ for each configuration, at a significance level of $\alpha = 0.01$. We observed some difference: configurations 1, 8, 15, and 23 exhibit varying numbers of cases, and some other configurations, such as 26, 33, and 37, appear or disappear after the classification obtained in $b = 4$. This highlights how sensitive the results are to the sampling process, especially in our case where the dataset is of small size.

<figure><img src="image_66.png" /><figcaption>Figure 4.27 – Distribution of the predicted configuration for each of the 496 cases for two predictors constructed during two repetitions $b$ of the algorithm in the case of an equal training sample size of $R = 100$ and with a significance level of $\alpha = 0.01$ for the test results. Colors distinguish between test results aligning with the Test Table (Table 4.1) (Yes in yellow) and those that do not (No in green).</figcaption></figure>
### 4.2.3 Strategies for further improving the prediction results

#### 4.2.3.1 Correction of distance bias

As highlighted in the previous subsection, the distance strongly influences the noise level in difference series and consequently the jump uncertainty. We propose a method to correct biases in the t-values of the non-collocated E-E' series. According to the rule n°2 presented in the paper, the series E and E' should be the same due to the large spatial correlation in the reanalysis. Any jump observed in E-E' is a manifestation of climatic differences between the locations of the two stations which may induce false detections in the non-collocated series. The idea is to estimate this bias from the E-E' series and correct the other series for it. This correction can potentially change the outcome of the significance test. If effective, we would expect an increase in significant jumps in G-E' and a decrease in G'-E, assuming that the change-point actually comes from G. Figure 4.28 shows the evolution of significance levels in three tests: after correction, G-E' sees an increase in significant tests, while on the contrary G'-E has more insignificant tests. These are slight changes but coincide with our expectations.

<figure><img src="image_67.png" /><figcaption>Figure 4.28 – Test count in three groups: 0 indicates no jump, 1 signifies a significant upward jump, and -1 represents a significant downward jump, from original results and after bias correction.</figcaption></figure>

Using this strategy reduces drastically the number of initial possible configurations from 36 to 10 because the E-E' test is now removed. The remaining configurations are: 1, 5, 10, 13, 15, 18, 23, 27, 31 and 35 (configurations containing a 0 for E-E'). This modification has a small positive effect on the test results. Figure 4.29 presents the configurations from the test results before and after this correction. The number of cases where the test results correspond to a configuration in the table increases despite the substantial reduction in rows. Number of cases initially match to configuration 1 and 15 increase. This shift can be attributed to the previously mentioned increase in significant tests in G-E' and a decrease in significant tests in G'-E.
<figure><img src="image_68.png" /><figcaption>Figure 4.29 – Configuration count before prediction: when considering the series E-E' and when ignoring it. A value of 0 signifies test results not present in the table.</figcaption></figure>

### 4.2.3.2 Multiple testing correction

The attribution method is based on the on simultaneous testing results from six series, which poses a multiple testing problem. This well-known issue posits that an increase in the number of tests made leads to a higher probability of erroneous conclusions. Simply, under the null hypothesis, the probability of a false rejection is denoted as $\alpha$ and thus, the probability of a true detection is $1 - \alpha$. For $k$ independent tests, the rate of correctly retaining the null hypothesis is $(1 - \alpha)^k$. Consequently, the probability of rejecting at least one test when all null hypotheses are true is expressed as $1 - (1 - \alpha)^k$, which called family wise error rate (FWER). The commonly used Bonferroni correction controls the FWER. If we test each hypothesis at a significance level of $\alpha$, we guarantee that the probability of having one or more false positives is less than $\alpha$. However, this control is often too strict. An other less restrictive measure is the False Discovery Rate (FDR) proposed by Benjamini and Hochberg (1995). It is defined as the proportion of false discoveries:

$$FDR = E[V/R], \qquad (4.12)$$

where $R$ represents the total number of rejected null hypotheses or positive and $V$ represents the number of erroneously rejected null hypotheses or false positive. When $R$ is equal to $V$, FDR becomes FWER. Control for FDR at level $\alpha$ consists in: (1) calculating the p-values obtained from the $k$ tests ordered from smallest to largest ($p_{(1)} \le p_{(2)} \le \dots \le p_{(k)}$), (2) finding the test with the highest rank $j$ for which $p_{(j)} \le \frac{j}{k}\alpha$ and (3) declaring the tests of rank 1, 2, ..., $j$ as significant.

Figure 4.30 compares the significance of the tests for each series of difference before and after the FDR correction. As expected, we observe an increase, but a slight one, of the insignificant test, except obviously for G-E. The effect is minor, but it can have a positive impact on the configuration obtained from the test results. These configuration numbers are given in Figure 4.31. Firstly, after the correction, there is a slight increase in the number of cases falling into configurations 1 and 15, while configurations 8, 22, and 31 show a slight decrease. An advantageous feature of FDR is its consideration of G'-E', which exhibits more significant jumps due to lower noise in comparison to other non-collocated series (depicted in Figure 4.22). Secondly, subsequent to the correction, there is a rise in the count of cases where the six tests align with a configuration in the table.
<figure><img src="image_69.png" /></figure>

Figure 4.30 – Test count in three groups: 0 indicates no jump, 1 signifies a significant upward jump, and -1 represents a significant downward jump, from original results and after FDR correction.

<figure><img src="image_70.png" /></figure>

Figure 4.31 – Configuration count before prediction from original results and after FDR correction. A value of 0 signifies test results not present in the table.

### 4.2.4 Conclusions

In this section, we have carried out additional analyzes in order to detail and discuss in more depth some of the results of the article (in particular on the identification of models and on the different stages of attribution, i.e. on the tests and the predictive rule).

Firstly, we detailed the results for studying the performance of the method used to characterize the data through simulations. The results revealed a relationship between the accurate identification of the noise model and the length of the series or the parameter values. Specifically, the ARMA(1,1) noise model requires a longer series for accurate identification compared to AR(1) and MA(1) models. Additionally, when the parameter is small (<0.3), there is a higher likelihood of misidentifying the AR(1) model as MA(1) or vice versa. Furthermore, the presence of data gaps can intensify model identification confusion and parameter estimate dispersion, but it has been shown that the impact remains small for reasonable fraction of gaps (e.g. < 20%, see 4.19).

Secondly, we analyzed in more details the prediction results from the significance tests, categorizing them based on their presence in the Test Table (Table 4.1 of the paper). Notably, the standard error of the jump estimates, and thus the test result, strongly depends on the data noise, which is influenced by the distance between the two series involved in the test. This relationship is linear for the G-E, G-G', and G-E' but not for
the E-E', G'-E, and G'-E' series. This is simply because jumps are more likely to be due to G rather than E, G', and E', as previously mentioned in the paper result.

Finally, we tested potential improvements to the prediction step by making bias correction and applying multiple testing correction. Both proposed modification yield positive, although, marginal improvements, such as a slight increase in the number of high probability configurations and a slight decrease in some low probability configurations. Because the improvements are only marginal, they are not further considered.
# Chapter 5

# Conclusions and perspectives

## 5.1 Conclusions

Water vapor plays a pivotal role in climate change due to its powerful feedback mechanism. Understanding and quantifying the change in water vapor within the context of climate change represent a prominent scientific challenge, as emphasized in recent research (Colman and Soden, 2021; Douville et al., 2022). While there is a global consensus regarding an upward trend in water vapor, several studies have noted discrepancies in trend estimates, particularly at the regional level, when comparing data from climate models, reanalyses, and observational sources (Parracho et al., 2018; Santer et al., 2021; Allan et al., 2022).

In this context, GNSS provides reliable long-term IWV data that are available under all weather conditions and independent from the reanalyses, which are useful to assess and validate the accuracy of reanalyses and climate models, thereby facilitating enhancements in these models. However, the presence of abrupt offsets in the GNSS data is problematic since trend estimates are sensitive to such inhomogeneities.

Within the overarching goal of providing a homogenized long term record of IWV data from GNSS, this study utilized a relative homogenization approach. The comprehensive homogenization process comprises four fundamental steps:

* firstly, generating the differences in IWV between GNSS and reanalysis data (used as a reference);

* secondly, segmenting this series;

* thirdly, attributing the identified change-points to either GNSS or reanalysis series;

* finally, correcting the IWV GNSS time series.

For the second step of detecting change-point, we applied the GNSSseg segmentation method developed by Quarello (2020) specifically for these data.

This thesis made contributions in two key areas. First, it provided enhanced insight on the sensitivity of the GNSSseg segmentation method and trend estimates to the GNSS and reference (reanalysis) data properties. Second, an innovative automatic attribution method that was developed, which is a crucial step, with the segmentation, in the global homogenization process.
**Sensitivity of the segmentation method and trend estimates to data properties.** This study was presented in the first paper (Chapter 3). The comparison of segmentation results of 81 common stations in two GNSS datasets (IGS repro1 and CODE REPRO2015) and two reanalyses (ERA Interim and ERA5) revealed the sensitivity of the segmentation method to different data properties. In particular, the segmentation method was found to be sensitive to variations in noise magnitude (heteroscedasticity) and the presence of periodic bias in the data. When we exchanged either the GNSS IWV dataset or the reference dataset, it was observed that only about half of the detected change-points remained similar. These sensitivities are mainly explained by the differences in the heteroscedasticity and biases in the GNSS and reanalysis data. Especially, it was shown that more recent GNSS processing or reanalysis data exhibit reduced noise and biases, and thus smaller representativeness differences between them. This consequently enabled the segmentation method to identify smaller jumps in the IWV difference series. Modification in the series length or auxiliary data, which includes surface pressure and the weighted mean temperature used in the conversion from ZTD to IWV, resulted in a much more substantial similarity of 70–80% in change-point detection. These factors subsequently influenced trend estimates, as evidenced by the mean trend and its dispersion.

We also observed a substantial impact of jump corrections on estimated trends. Notably, the uncertainty in the estimated trends resulting from the use of different reprocessed GNSS datasets or reference reanalyses in the segmentation ranged from 0.015-0.019 kg m⁻² year⁻¹ (RMS difference between the tested pairs of data sets) when only the validated change-points were used and 0.002 to 0.012 kg m⁻² year⁻¹ when all the change-points were used. However, on a global scale, homogenized GNSS data (corrected using either validated or all change-points) and both reanalyses agreed on a global mean IWV trend estimate of 0.027-0.030 kg m⁻² year⁻¹, indicating a global moistening trend close to 7% K⁻¹ over the past 2.5 decades. The longer time series (25 years) provided higher accuracy in trend estimates, with a mean standard error of 0.018 kg m⁻² year⁻¹ and a dispersion of approximately 0.03 kg m⁻² year⁻¹ across the global network.

This study also addresses two main limitations in the homogenization process. First, regarding the segmentation method, we observed that some stations exhibit interannual and longer-term variations in the noise and bias characteristics which are not currently modeled in the segmentation method (both the noise variance and the bias are assumed periodic with a fundamental period of 1 year). At those sites, interannual and/or decadal variations in bias especially were captured by the segmentation. The long-term variations could be reduced, for example, through the use of a better reference reanalysis or a nearby GNSS station, rather than a reanalysis. Second, only about 35% of change-points are validated by metadata, meaning that 65% of the detections can be undocumented changes in GNSS or changes in reanalysis data. This highlights the critical role of the attribution step, which aims at determining the origins of identified change-points.

**Development of an automatic attribution method.** The method developed in this thesis relies on a the combination of GNSS and reanalysis data from a main station with similar data from nearby stations. For each detected jump in the main series, and each associated nearby station, six series of differences are formed from the four base series and a test of significance is performed on the jump. The development of the proposed method involves several meticulous steps. Firstly, we characterize the heteroscedasticity and autocorrelation structure in the six series of differences. The considered noise models are white noise, AR(1), MA(1), and ARMA(1,1). Notably, we show that only a tiny number of series are identified as white noise. Secondly,