Preprint title: Bayesian modelling of piecewise trends and discontinuities to improve the estimation of coastal vertical land motion

One of the major sources of uncertainty affecting vertical land motion (VLM) 10 estimations are discontinuities and trend changes. Trend changes are most commonly caused 11 by seismic deformation, but can also stem from long-term (decadal to multidecadal) surface 12 loading changes or from local origins. Although these issues have been extensively addressed 13 for Global Navigation Satellite System (GNSS) data, there is limited knowledge of how such 14 events can be directly detected and mitigated in VLM, derived from altimetry and tide-gauge 15 differences (SATTG). In this study, we present a novel Bayesian approach to automatically 16 and simultaneously detect such events, together with the statistics commonly estimated to 17 characterise motion signatures. Next to GNSS time series, for the first time, we directly 18 estimate discontinuities and trend changes in VLM data inferred from SATTG. We show that, 19 compared to estimating a single linear trend, accounting for such nonlinearities significantly 20 increases the agreement of SATTG with GNSS values (on average by 0.36 mm/year) at 339 21 globally distributed station pairs. 22 The Bayesian change point detection is applied to 606 SATTG and 381 GNSS time series. 23 Observed VLM, which is identified as linear (i.e. where no significant trend changes are 24 detected), has a substantially higher consistency with large scale VLM effects of Glacial 25 Isostatic Adjustment (GIA) and contemporary mass redistribution (CMR). The standard 26 deviation of SATTG (and GNSS) trend differences with respect to GIA+CMR trends is by 27 38% (and 48%) lower for VLM which is categorized as linear compared to nonlinear VLM. 28 Given that in more than a third of the SATTG time series nonlinearities are detected, the results 29 underpin the importance to account for such features, in particular to avoid extrapolation 30 biases of coastal VLM and its influence on relative sea level change determination. The 31 Bayesian approach uncovers the potential for a better characterization of SATTG VLM 32 changes on much longer periods and is widely applicable to other geophysical time series. 33 Julius Oelsmann Deutsches Geodätisches Forschungsinstitut der Technischen Universität München, Arcisstraße 21, 80333 Munich, Germany Tel.: +49 (89) 23031-1263 Fax: +49 (89) 23031-1240 E-mail: julius.oelsmann@tum ·Marcello Passaro · Laura Sánchez · Denise Dettmering · Christian Schwatke · Florian Seitz same address as first author 2 Julius Oelsmann et al.

To answer our research questions and to test our method, we apply the Bayesian model to 175 VLM time series from GNSS and SATTG, as well as to synthetically generated data. We use  In this research, we analyse VLM time series which are derived from SATTG differences 184 according to the recipe in Oelsmann et al. [2021]. In order to increase the quality and quantity 185 of altimetry data close to the coast, we use dedicated choices in terms of range and corrections 186 needed to estimate sea surface height (see Table 1). We use along-track altimetry data of combining the individual mission time series and monthly PSMSL data, the along-track 205 data is downsampled to monthly means to match the frequency of TG observations. The 206 correlations are computed independently for missions which share the same nominal track. 207 We spatially average the along-track data in the ZOI and compute the differences between 208 their monthly averages and the TG data. Furthermore, the following data selection criteria are 209 applied: We omit time series where the multi-mission, monthly SAT time series (averaged 210 in the ZOI) present a correlation with the TG data lower than 0.7 (i.e. ∼ 10th percentile of 211 all data) and a Root-Mean-Square (RMS) error higher than 5.5 cm (∼ 90th percentile of all 212 data). We only use SATTG time series with a minimum of 150 months of valid data, which 213 yields a number of 606 remaining SATTG estimates.   The modelled time series features are a trend, a harmonic annual cycle and a noise term.

237
[2021]). In contrast to the GNSS data, annual variations in SATTG data can however also 238 stem from discrepancies in the observations of the different techniques. As we show in the 239 following, these non-geophysical deviations can have much larger amplitudes than those 240 obtained from GNSS data and also influence the noise characteristics.

241
Several studies affirmed that a combination of white noise (WN) and power law noise . We use a spectral index of -0.9, which is close to flicker noise process, 246 and amplitudes of 2mm/year and 6mm/year −k/4 for white and coloured noise, respectively.  we take these different noise properties into account by testing the detectability of different 263 discontinuity-to-noise ratios, instead of absolute values of discontinuities. 264 We perform three experiments in which we vary (1) only the discontinuity-to-noise-ratio,

265
(2) the trend and (3) the number of change points, together with discontinuities and trends.

266
The full setup is described in Table 3.   Our overarching goal is to detect the most common time series features in GNSS and SATTG 299 data using a single comprehensive model. The major components considered in this study 300 are discontinuities o(t) (abrupt changes in height), trends g(t), a seasonal term seas and a 301 noise term η, which can also be identified in Fig. 2: Here, y(t) denotes either GNSS or SATTG observations at time t and is described with  of the trend component must be corrected for these discontinuities as follows: is: with x(t) ∈ 0, 1: Finally, the noise η in Eq.
(1) is approximated as a first order autoregressive process 335 AR(1). We emphasize that the presented model setup explicitly allows for trend changes, The positive autocorrelation coefficient φ and the white noise amplitude σ 2 w are both 362 drawn from halfnormal distributions with σ φ andσ w , respectively. Finally, we approximate 363 all the other parameters, the trend and discontinuities o, p, k, h and the monthly means m with 364 normal distributions. Hence we obtain the following set of unknown parameters of the model: parameters to be estimated.

369
In addition to the type of probability distribution P(Θ ), we also specify initial values of 370 the associated distribution parameters. Here, we make use of prior knowledge of common 371 GNSS and SATTG time series characteristics, to enhance accurate parameter estimation. As we define a so-called informative prior for q 0 , which expresses specific knowledge of the 376 expectation of a change point to occur. In this case a low probability of q 0 =10% is assigned. 377 We also define other initial settings, which are more thoroughly explained in the Appendix A.
378 Table 4 summarizes the complete model setup and initial assumptions. Note, that these initial 379 values are set for the normalized time series. 380 We use different MCMC samplers to generate inferences about the desired target distri-  We observe that the Bayesian ∆ PW estimates in the discontinuity and the trend ex-   Overall, we obtain relatively high TP detection rates (50% -100%), compared to previ-  The third mechanism which contributes to potential trend changes is nonlinear surface 532 deformation due to mass load changes. In the last row of Fig.4

574
Despite the abundance of time series, which are affected by both, discontinuities and trend 575 changes, in the majority of cases discontinuities are not necessarily associated with a trend 576 change (such as in Fig. 5(b) or in the GNSS time series in Fig.5(c)). In order to mitigate such

604
The percentage of improvements refers to the absolute deviations of trends as also listed in 605   Table 5.

606
There are only nine cases where more than two change points are detected. Here, the 607 accuracy of trends using the piecewise estimation decreases compared to the linear estimates.

608
This could be due to the increased fragmentation of the data and shortness of the time series 609 segments. Such small number of samples (9), however, hinders a robust assessment of the 610 significance of improvement/deterioration. In general, the lower consistency achieved in such The geographical distribution of the differences (mm/year) between ∆ LIN and ∆ PW is 622 illustrated in Fig.7. Improvement (deterioration) with respect to a linear trend estimation is

638
In some cases DiscoTimeS trend estimates yield a lower accuracy compared to the single 639 linear trend estimates. Some of these cases are located in Great Britain (Fig. 7(b)) and 640 Japan (Fig. 7(c)). There are various possible reasons which might explain such degradation.   for the case of no detected change points. We obtain the best agreement when also including 678 the CMR correction compared to using the GIA estimate only. high standard deviation in trend changes is detected (see Fig. 8b and Table 6 second section).

696
As for SATTG VLM estimates, the combined GIA+CMR effect improves the comparability 697 compared to the sole GIA VLM correction. σ w would need to be set individually. We normalize the data by the median of their 2-year 780 running-standard-deviation, hereinafter called σ obs . With this approach we avoid that extreme 781 discontinuities (in particular present in GNSS data), which present orders of magnitudes 782 larger than the 'true noise amplitude' influence the normalization. We also subtract the offset 783 of the first observation from the data.

784
Next to the initial probability of q 0 , which is explained in section 3.3 several other param- from the ensemble member with the best CV-score.

820
Using CV (or other criteria such as WAIC (widely applicable information criterion) or 821 DIC (deviance information criterion)) to select a single best-performing realization, can 822 however lead to overfitting of the data and introduce a significant selection bias [ using best loo to select the best GNSS realization. The fact that we obtain best results when we 847 choose the chain with the lowest number of effective parameters for SATTG ( (lowest p−loo )), 848 indicates that using best loo instead might lead to overfitting of the data, as also discussed by 849 Piironen and Vehtari [2017]. The necessity to apply different selection schemes is most likely In a similar way, we compute ∆ LIN to analyse the differences between single linear The example of the real data trend comparison can also be transferred to the sensitivity 881 experiments. Here, the piecewise SATTG fit can be thought of as the synthetic data fit and