Reduction of apparent temporal variations of tidal parameters by a proper local response model

We describe a new harmonic tidal analysis method, which constrains the solution to be near a reference model. This regularization stabilizes the linear regression, allowing us to infer model parameters for each tidal harmonic. This overcomes the need to create a priori groupings of harmonics. The inversion is done iteratively by adjusting the reference model to reduce the data misfit. The frequency dependence of the solution is thus data-driven. We find models for the different spherical degrees independently. Our procedure allows narrow-band variations of the tidal admittance. We test the hypothesis that some of the temporal variations of tidal parameters found in previous studies were caused by inappropriate body tide models in combination with a priori wave grouping. We determine a local response model from 11.5 years of data recorded by the superconducting gravimeter SG056 at Black Forest Observatory (BFO, Schiltach). Using this as an a priori model in a non-regularized moving window analysis of wave groups composed from summed harmonics, we find that periodic variations of groups M1, K1, μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}2, N2, L2, and S2 are reduced by up to a factor of 7 compared to earlier studies. Some variations previously seen in the M2 group are captured as well.

(tidal harmonics; Doodson 1921), using a linear regression to determine the complex amplitude factors (tidal parameters).Cartwright (1999) gives a historical overview of the development of the harmonic method, and Wenzel (1997a) provides the theoretical basis of recent implementations.Dierks and Neumeyer (2002) present a review and comparison of three implementations.In the gravity community, two widely used programs are BAYTAP-G (Tamura and Agnew 2008;Tamura et al. 1991) and Eterna (Wenzel 1996(Wenzel , 2022)); the latter has become a de facto standard in the European tidal gravity community.Schüller (2015Schüller ( , 2020) ) presented an extended and improved version called ETERNA-x.All of these approaches seek model parameters for a priori defined groups (sums) of harmonics which minimize the residual between the model predicted tidal signal and the recording without additional constraints.Meurers (2004, his Fig. 5) presents a moving window analysis (MWA) of the gravity data from the superconducting gravimeter (SG) at Vienna for O 1 and M 2 , showing an annual variation of gravimetric factor for M 2 .The variation of both factors is in the order of magnitude of 0.2 per mille.The admittance of the solid Earth is not expected to vary on these time scales.Meurers (2004) discusses model errors in the air pressure correction as a possible cause.Similar observations are reported by Jahr (2015).Meurers et al. (2016) extended this investigation for M 2 to nine European SG stations.The relative amplitude of variation of gravimetric factor there is in the range from 0.14 per mille to 0.33 per mille, where amplitude increases with increasing amplitude of ocean loading for M 2 (Meurers et al. 2016, their Table 3 and Fig. 3).They confirm the annual variation of parameters and suggest an additional variation with 8.85 years period.They discuss ocean loading as a possible cause, along with numerical artifacts in the analysis, and artifacts caused by unstable calibration and improper preprocessing.According to these authors, the 8.85-year period variation might point to ocean loading being different for forcing at spherical degree 2 and degree 3, a problem recently investigated by Sulzbach et al. (2022).The variation they observe could be explained by tidal parameters depending on frequency and spherical degree.The available length of the time series does not allow for the required frequency resolution in the tidal analysis.Merriam (1995) already reported satellite harmonics of M 2 in tidal analysis residuals, which are equivalent to an annual cycle of M 2 admittance with a relative amplitude of 0.54 per mille.They are found in an analysis of data from the SG at Cantley using the response method.He suspected large-amplitude ocean loading in the Bay of Fundy to cause this signal.Systematic variations of tidal parameters for the SG056 at BFO Schiltach with different periodicities are reported by Schroth (2013) for different tidal wave groups.Schroth et al. (2018) present a detailed catalog of temporal variations seen for 19 European and global SG stations in moving window analyses of 12 wave groups.For European stations, they report annual variations for M 2 of up to 0.8 per mille and for K 1 in the order of one per mille, as well as semiannual variations for S 2 of a few per mille.Most remarkable are the variations with a period of 8.8 years which they find for N 2 in the order of 1.5 per cent and for L 2 in the order of 10 per cent.Data from some stations close to the sea are found to produce even larger and irregular variations (e.g., Onsala, Ny-Ålesund, and Syowa).
In their discussion, Schroth et al. (2018) rule out effects of instrumental origin.As possible causes, they consider shortcomings of the body tide model used in the analysis, a time-dependent response of the entire Earth including ocean loading, and cross talk between analysis parameters.Variations seen with a period of 8.8 years strongly suggest an inappropriate ratio of tidal admittance to forcing at spherical degree 2 and degree 3, since the affected groups contain satellite harmonics of significant amplitude and different degrees.Eterna 3.40, as used by Schroth et al. (2018), assumes the ratio between the admittance of the degree 2 tides and the degree-3 tides given by the Wahr-Dehant-Zschau model (Wahr 1981;Zschau and Wang 1985;Dehant 1987;Dehant et al. 1999).This is a theoretical body tide model and by definition does not account for ocean loading or atmospheric signals.A ratio not matching the total response of the Earth would create an incorrect modulation in the synthetic tidal gravity time series.Schroth (2020) investigates loading contributions from the oceans and the atmosphere as well as non-tidal radiation driven signals as a source of temporal variations of the apparent admittance.
Most tidal analysis methods do not support the separate handling of tidal harmonics of different degrees.Investigations of the admittance for degree 3 tides, as presented by Melchior et al. (1996) and Ducarme (2012), are rare.ETERNA-x (Schüller 2015(Schüller , 2020) ) provides means to separate the degree 3 admittance and is used by Sulzbach et al. (2022) for this reason.
The duration of the time series limits the obtainable frequency resolution.With conventional implementations, it is impossible to handle all tidal harmonics listed in catalogs like those by Tamura (1987) and Hartmann and Wenzel (1995) individually.A priori wave grouping is required, usually based on the Rayleigh criterion, though not taking into account the actual resolution limited by the signal-to-noise ratio of the data (Munk and Hasselmann 1964;Ducarme and Schüller 2018).Tidal analysis then seeks only an admittance for the sum of harmonics in a group, whose relative sizes are as predicted by the a priori tide model.An incorrect frequency dependence of the model will cause beats (time variations) that are also incorrect (Schroth 2020).Schüller (2015) argues that wave grouping should always occur at a maximum possible resolution.
We introduce an approach to tidal analysis that uses a constraint with respect to a reference model.This allows us to abandon a priori wave grouping.We thus can search for harmonics which show a different dependence on frequency and spherical degree than assumed in the a priori body tide model.We use this method to create a local response tidal model, without time variations in its parameters, for SG056 at BFO Schiltach using 11.5 years of data.In this model, the admittances to tidal forcing at different harmonic degrees are kept separate.We also separate tidal parameters directly affected by solar radiation at S 1 and allow for a detailed frequency dependence in the diurnal (K 1 ) and semi-diurnal (M 2 ) band.We use these results as our a priori model in a MWA, using 12 wave groups without any prior constraint.This analysis shows much smaller time variations of tidal parameters than found by Schroth et al. (2018), demonstrating that the local response model removes causes of false temporal variations in tidal parameters.

Theory and method
Based on a tidal catalog giving amplitude and phase for a harmonic development, we compute the gravity signal for a rigid oceanless Earth.The approach we take can be applied to any tidal observation, be it gravity, strain, tilt, or surface displacement.The theory introduced in the following uses symbols appropriate for observation of tidal gravity.For measurements of gravity, the complex admittance of the measurements referred to g rig (t), is expressed by a pair of tidal parameters; the 'gravimetric factor' δ( f ) and the 'phase lead' ( f ), both functions of frequency f .We use δ l and l for these tidal parameters for tidal harmonic g rig l (t) at frequency f l .
In a linear regression approach, we fit the synthetic tidal gravity signal g syn (t) to both 'observed gravity signal' g(t) and local barometric pressure p(t), the latter to account for gravity fluctuations caused by changes in the local atmosphere.Signal g syn (t) is computed by a development of tidal harmonics g rig l (t) scaled by δ l and shifted by l .We compute the synthetic data for an a priori body tide model and seek for adjustment factors to these parameters, by minimizing the data misfit in a least-squares sense.
Unlike conventional approaches such as those described by Wenzel (1997b), we add a model constraint to the objective function, which is commonly known as 'Tikhonov Regularization' or 'Ridge Regression.'The parameters are constrained to be close to reference values, deviating from the reference only if justified by a significant decrease of data misfit.This removes the instability created by trying to determine parameters for tidal harmonics closely spaced in frequency or of very small amplitude.In conventional approaches, this instability is removed by summing harmonics closely spaced in frequency, a so-called wave group.Wave grouping is a model parameterization where tidal harmonics are organized in frequency bands (Venedikov 1961;Chojnicki 1973).Harmonics within the same frequency band are scaled by a common factor.Individual harmonics cannot be distinguished in the signal.
We next define the mathematical symbols and elementary equations for the new approach, and some are summarized in the glossary (Appendix A.1).

Tidal signals
With a tidal catalog of L harmonics, the rigid Earth tide at a given location is with harmonic frequency f l .The amplitude A l and the phase φ l are computed from the information given in the tidal catalog based on the station's location on the Earth's surface.
For simplicity in the description, we neglect the slight time dependence (e.g., Bartels 1957;Cartwright and Tayler 1971;Tamura 1987;Hartmann and Wenzel 1995) of amplitude A l (t) and the exact astronomical arguments (e.g., Simon et al. 1994;Tamura 1987), which results in a nonlinear time scale in the cos-function.The essential point is that g rig l (t) can be computed precisely for each given point in time.
The synthetic tidal gravity signal for a given model accounts for the admittance of the solid Earth by the tidal parameters δ l (gravimetric factor) and l (phase lead).To express the tidal analysis problem with linear parameters, we rewrite the synthetic signal as with the rigid Earth tide signal and its quadrature We use symbols similar to those used by Wenzel (1997a).The 'tidal model' is specified by the linear parameters The conventional model parameters then are gravimetric factor

Objective function
To set up the inverse problem, we express the signals by the time series samples where k = 1 . . .K and p(t) is a recording of local air pressure and g(t) is recorded gravity.All of these signals are filtered and tapered (application of a time-domain window function) prior to the analysis in a consistent way.We compute the synthetic tides for an a priori (initial) model with parameters X ini l and Y ini l .Quantities X and Y are real and imaginary parts of the complex admittance to which the adjustment factors x l and y l are applied.The final model parameters then are and Factors x l and y l initially are set x l = 1 and y l = 0. Hence, the factors X l and Y l mostly vary smoothly with frequency f l , but the initial model also captures other frequency dependence, notably the resonance due to the free core nutation (FCN) in the diurnal band.The structure of Eq. ( 9) is such that the adjustment factors can introduce an additional phase shift even in the case where Y ini l might be zero.In that sense, the initial parameters define an initial synthetic tidal signal and its quadrature The synthetic tidal gravity signal in Eq. ( 2) thus is Beside minimizing the data misfit only, as in conventional tidal analysis (Wenzel 1997a, b), we also add a model constraint to the objective function Here, R ini is an a priori admittance to air pressure and r is the adjustment factor to this.
The second term (in curly braces) of Eq. ( 13) expresses a model constraint normalized by the number of parameters 2L + 1.When minimizing with respect to x l , y l , and r , model parameters that are close to their reference values x ref l , y ref l , and r ref , respectively, will be preferred.Model parameters will only deviate from the reference if the data misfit expressed in the first term of Eq. ( 13) is significantly reduced.The solution of the optimization problem seeks a compromise between minimizing the data misfit on the one hand and keeping model parameters within range of the reference values on the other hand.The trade-off parameter α adjusts the emphasis on one criterion over the other.Wave grouping is not needed.
The σ should specify the expected rms amplitude of the noise level in the data in units of data which makes the expression dimensionless.A normalized data misfit near unity indicates that the fit is successful with respect to the noise level.A smaller value would indicate overfitting of the data.

Finding the minimum
The parameters x l , y l , and r which minimize the objective function in Eq. ( 13) are the solution of a system of linear equations, which is found by searching the stationary point with as is common for least-squares problems.To set up the system of linear equations, we collect the regressors in the so-called forward operator matrix G G G with Model adjustment factors are collected in parameter vector m with We let R ini have units of 1 nm s −2 hPa −1 so that all columns in G G G have units of 1 nm s −2 and all m l (l = 1 . . .2L + 1) are dimensionless.The signal vector of synthetic gravity then is and where the elements g k of vector g are the samples of observed gravity and δ m = m − m ref is the update of the adjustment factors (model parameters) with respect to their reference values.We abbreviate the scaling factors By substituting we find the system of linear equations Its solution a satisfies the least-squares condition in Eq. ( 15).We find the solution by computing the singular value decomposition of the matrix where U U U and V V V are the orthonormal matrices of eigenvectors in data space and model space, respectively, and where is the diagonal matrix of singular values, of which N ≤ (2L + 1) are nonzero with λ j > 0. The factorization has to be computed only once to find the solution for any value α, which controls the filter factors in the diagonal matrix diag(η j ).From this solution, the vector of model parameter adjustment factors is obtained by

Finding the optimal parameters
The solution in Eq. ( 29) to the optimization problem in Eq. ( 14) not only depends on the recorded gravity g and air pressure p.It is deliberately controlled by the choice of the reference model m ref and the trade-off parameter α, in particular.We commonly choose the value of α by trial and error from a set of solutions computed for different values.
For large α, the improvement of the model will be small.At the opposite end, for small α the data will be overfitted by using unreasonable model parameters.This can be expressed by the trade-off between distance to the reference model and data misfit Because of the minimization condition in Eq. ( 14), d (α) necessarily decreases monotonically with decreasing α, while m (α) increases.While d (α) by definition decreases with decreasing α, this is not the case for the misfit of the untapered synthetics G G G m(α) with respect to the untapered data g .Where d (α) increases with decreasing α while d (α) further decreases, further update of the model is not appropriate to better fit the actual tidal signal.Instead the model starts to fit a significant fraction of the noise in the data.Quantity σ in Eq. ( 32) is the expected rms noise level in the (untapered) data, while σ in Eq. ( 31) and used in Eq. ( 13) accounts for the taper loss.
To search for the optimum, we chart m (α) against d (α) as shown in Fig. 1.Generally, at large α the model parameters are tied to the reference values m ref , m (α) approaches zero and d (α) saturates (Fig. 2).With decreasing α, an increase in m (α) indicates that model parameters are adjusted in order to reduce data misfit.There is typically a range of α where d (α) significantly decreases with moderate increase in m (α).Eventually d (α) approaches 1, and the rms level of residuals equals that of the expected estimated noise.With further decreasing α, the value of m (α) monotonically increases, while there appears a threshold above which d (α) starts to increase again, which indicates that the procedure increasingly fits noise in the tapered data (overfitting).Values of α below this threshold (trade-off point at m ≈ 0.01 in Fig. 1) must be avoided.Usually values of α will be chosen larger than this threshold.In the cases displayed in Fig. 1, the optimal value is chosen at the first trade-off point, ( m ≈ 10 −4 , d ≈ 0.88) because a strong scatter of model parameters is observed for smaller values of α as discussed below for Fig. 3. Figure 2 compares the curve for d (α) with that for d (α), which clearly indicates overfitting for α < 0.25.
In addition to the misfit curves, we investigate the amplitude spectrum of g − G G G m(α), which should indicate that a decrease of d (α) is due to a reduced amplitude of remainders of tidal harmonics in the residual time series.

Unconstrained analysis and wave grouping
In cases where we aim to find the optimal parameters based on the recorded gravity only (which is the case in moving window analysis), we set α = 0 and resort to a wave grouping approach.The solution nevertheless is computed with the above given expressions.However, we then compose wave groups by summing over columns of matrices C C C ini and S S S ini (their matrix elements are C ini kl and S ini kl , respectively) in Eqs. ( 11a) and (11b), thus reducing the number of columns.The column index l from that point on (in particular for the derived factors x l and y l ) does no longer specify the tidal harmonic, but the wave group of tidal harmonics.
Columns contributing to the same wave group are not necessarily adjacent to each other.In particular, we may choose all columns for tidal harmonics of degree 2 to go into one group and those of degree 3 to go into a different group within a given frequency interval.In the data analysis discussed below, we use an unconstrained approach only in the framework of a moving window analysis and use a rather conventional definition of wave groups there.

The rigid Earth tide signal
Although in principle any tidal catalog in the form published with Eterna 3.40 (Wenzel 1996(Wenzel , 2022) could be used for this method, we decided on the tidal catalog compiled by Tamura (1987).This catalog lists 1200 harmonics and is sufficiently accurate to capture all tidal signals resolvable in the recorded data.Using this catalog speeds up the computation of the singular value decomposition in comparison with using the 12,935 harmonics listed in the catalog compiled by Hartmann and Wenzel (1995).
Based on the tidal catalog, we compute the 'rigid Earth tides' time series with a modified version of predict from the Eterna 3.40-package (Wenzel 1996(Wenzel , 2022)), predict rigid.The program is modified such that it outputs all C rig kl and its quadrature S rig kl as defined in Eq. ( 8) and not just their sum.

The initial tidal parameters
As the initial model X ini l and Y ini l , defined in Eqs. ( 3) and ( 11), we frequently use the theoretical body tide model as implemented in Eterna 3.40 (Wenzel 1996(Wenzel , 2022)).It is based on the contributions by Wahr (1981), Zschau and Wang (1985), Dehant (1987), andDehant et al. (1999).In Eterna, the amplitude factors are taken for a rotating, elliptical Earth with parameters of the Preliminary Reference Earth Model (PREM Dziewonski and Anderson 1981).We call this tidal model 'WDZe.'Though the theory behind it is viscoelastic, Eterna 3.40 takes account of the relaxation of elastic moduli, but neglects the phase; thus, Y WDZe l = 0. We use the updated FCN (free core nutation) period derived by Krásná et al. (2013) from VLBI, which is consistent with the theory for a non-hydrostatic inelastic Earth as presented by Dehant et al. (1999).This updates the original value used in Eterna 3.40, which was taken from Wahr (1981, his Table 1).
When body tide parameters (gravimetric factor and phase) for a station specific model are available, we use specific values for X ini l and Y ini l , as derived by a previous tidal analysis for the respective station (e.g., BF21 ).As will be demonstrated below, this can have significant consequences in a moving window analysis.

The solution of the inverse problem
The program used to solve the inverse problem in the sense of Eqs. ( 25)-( 29) is newly implemented in C++.We call it RATA (Regularization Approach to Tidal Analysis).Essential tasks of linear algebra are handled by the Eigen-library provided by Guennebaud et al. (2010).
A filter is applied to all time series consistently.The synthetic time series sequence as well as observations is convolved with the same zero-phase finite-impulse-response (FIR) band-pass or high-pass filter.While in Eterna 3.40 the synthetic tidal signals C C C ini , and S S S ini are only scaled by the filter gain corresponding to the frequency, we apply the full FIR filter.Thus, we avoid residuals, which otherwise would result from irregularities in the time scale, which are mainly due to leap seconds.
RATA offers means to run a moving window analysis.For a longer time series, the regressors p, C C C ini , and S S S ini as defined in Eq. ( 8) need to be computed only once.The program then subdivides all long time series into shorter, overlapping segments and computes the inversion result for each segment separately.
Prior to the actual tidal analysis for the respective data segment, the program applies a taper to all time series consistently.Tapering in tidal analysis may be discussed controversially.In a regression problem where data can be completely represented by a set of harmonic functions with known frequency and phase, there is no risk of spectral leakage and tapering is dispensable.Tapering, however, helps to reduce spectral leakage in the recorded data from signals that are not present in the regressors (Schüller 1976).As in Eterna 3.40 we apply a cosine window taper to the time series, which is the square root of a Hann window taper (Harris 1978).We consider this a compromise between loss of signal energy and mitigation of spectral leakage.

Data
We take level 2 data (60 s sampling interval, correction code 32) for SG056 sensor 1 at Black Forest Observatory [BFO] (1971) and locally recorded air pressure as provided in the IGETS data base (International Geodynamics and Earth Tide Service 2017; Voigt et al. 2016).The instrument is described by Forbriger and Heck (2018).Zürn (2014) gives a comprehensive overview of BFO.The data are preprocessed by EOST (Boy 2019) and cover the time period from 2009-10-01 to 2021-03-31 (i.e., 11.5 years).The preprocessing comprises the scaling to units of acceleration and the removal of steps.Gaps in the original recording as well as transient signals (earthquakes, glitches) exceeding a threshold in the time derivative of the signal are flagged and replaced by synthetic tides in case of gravity and synthetic data from the MERRA-2 model in the case of air pressure.Boy et al. (2020) discuss consequences of this approach.We use the gap-filled version of the data.We carefully checked the flagged time windows (where gaps where filled) in the residuals of our analyses and could not find signatures of synthetic data above the usual noise level.For tidal analysis, we decimated the data to a sampling interval of 15 min after application of the finite impulse response (FIR) low-pass filter n20s5m02, which is provided together with Eterna 3.40 (Wenzel 1996(Wenzel , 2022)).The gain of the filter is less than 10 −7 in the stop band and deviates less than 10 −7 from one at frequencies of 4 cpd and smaller.After decimation each time series then comprises 403,170 samples.Before tidal analysis we apply the n10m10m2 FIR high-pass filter (provided with Eterna 3.40) consistently to observations and regressors.The corner frequency of the high-pass is 0.48 cpd.We hence need not account for drift in the observations.The analyzed frequency band thus is 0.48 cpd to 7.2 cpd and well captures tidal constituents from the diurnal, semi-diurnal, and ter-diurnal band as well as higher frequencies.
We use local air pressure to at least partly correct for the impact of mass fluctuations in the atmosphere as suggested by Warburton and Goodkind (1977), Merriam (1992) and Crossley et al. (1995).groups only, because of the window length of three months 1 .Schroth et al. (2018) point out that if the ratio of harmonics within these groups should be inappropriate, the beating (amplitude modulation) of the synthetic tidal signal would differ from the observations and in consequence cause a periodic variation of tidal parameters as demonstrated by Schroth (2020, 4.1.2).Such variations are reported for some of the groups.Possible inadequacies of the model, as discussed by these authors, in particular are: (1) The fixed ratio between the admittance of Earth to tidal forcing of degree 2 with respect to admittance at degree 3 might be the cause for long period (8.8 years) variations found for L 2 and N 2 .(2) A modulation of ocean loading with annual period, as proposed by Müller et al. (2014), might cause the observed variation in M 2 .(3) The diurnal gravity signal contains a component of masses moved by thermal radiation rather than tidal forces (the socalled radiation tides) and commonly causing an anomaly in the parameters for S 1 , which is subsumed with the K 1 group, which in turn shows a strong variation of tidal parameters.(4) As well an inappropriate model of the free core nutation (FCN) resonance is mentioned by Schroth et al. (2018) as a possible inadequacy in K 1 .We rule out the latter in this case, because the FCN parameters obtained from VLBI (Krásná et al. 2013) and used by Schroth et al. (2018) are much more accurate than what could be inferred from gravity recordings.
In the discussed cases, time-invariant properties in the tidal signal, which are not well represented by the uniform admittance assumed for a wave group, can cause a temporal variation in the MWA.If the model of the tidal admittance within a group would be adjusted appropriately, these variations would be captured by a local model of tidal admittance.In a first step, we seek such a model for the gravity data spanning 11.5 year, by adjusting parameters for each harmonic listed in the tidal catalog by Tamura (1987).

Approach with model constraint
Even with 11.5 years of data, many of the harmonics in the catalog are too close to each other in frequency to be discriminable or are too small in amplitude to be well constrained by the observations.The regression problem would be ill-conditioned, and many parameters would be strongly affected by noise.To overcome the instability, we regularize the regression as introduced by the objective function in Eq. ( 13).With a large α, the regression would just reproduce the reference model.By gradually reducing α, we explore for which of the parameters a deviation from the reference is rewarded by a decrease in data misfit.
The first parameters to depart from their reference are those for harmonics of large amplitude.The parameters of their small-amplitude (by orders of magnitude) neigh-1 In fact, 83 days or less due to filter length and gaps.
bors stay at the reference values or account for noise in the data.We would not accept a tidal model if the parameters of minor harmonics differ from their immediate large-amplitude neighbors.Due to their small amplitude, a difference to the admittance with respect to their neighbors is not likely to be resolved since they are prone to noise.Moreover, there is no physical reason to allow these harmonics to let their parameter estimates depart from those in proximity, except at known resonances.This is in agreement with the 'credo of smoothness' (Munk and Cartwright 1966).We therefore adjust the values for the reference model in a finite frequency band to the new parameters assumed by the large-amplitude harmonics.We do this by defining 'reference groups,' i.e., groups for the reference model.Such a group is then a definition which harmonics share the same values of reference model parameters x ref l and y ref l .No grouping in the sense of summation takes place, and the number of free parameters remains constant throughout the whole procedure.Such parameterization still allows parameters (adjustment factors) x l and y l of any harmonic to deviate if further distance to the reference model should be justified by the data in the following analysis, and hence, the more significant harmonics could be inferred and the resolution is not lost.
In the subsequent iteration, we solve the regression problem with a constraint to this updated reference model.Where the parameters of the reference model are appropriate to express the observed signal, the parameters of the harmonics will align at the reference values.Where the large-amplitude harmonics keep further departing from the new reference level, the reference model gets updated once again.If parameters of two significant harmonics within one reference group depart in opposite direction, we split the reference group into two or more subgroups, to account for the conflicting constraints imposed by the observed signal.This way we iteratively approach a detailed model of tidal admittance, where adjusting reference groups based on constraints imposed by the data is part of the procedure.In this procedure, groups are not defined solely along frequency.We separate harmonics of degree 2, three, and four into different reference groups.The frequency bandwidths of the reference groups for degree 3 and four are much larger than for degree 2. The data do not call for a stronger frequency dependence for admittance of degree 3 and four.
In brief, we manually design a tidal model based on the solution of the model constraint regression.In each iteration, we use this model as a reference and in that test, whether the optimized tidal parameters stay with the value of the model.If not, we update the model accordingly and enter the next iteration.

Tidal catalog and parameters
We use the full Tamura (1987) catalog with L = 1200 harmonics.Throughout the entire procedure, we set σ = 1 nm s −1 , which corresponds to σ ≈ 1.4 nm s −1 in the case of a cosine taper function.We use the WDZe-model for X ini l and Y ini l (l = 1 . . .L) and R ini = 1 nm s −2 hPa −1 as the a priori model to setup G G G in Eq. ( 16).All parameter adjustments will be expressed with respect to this model.For the initial iteration, we set the reference parameters to We explore the misfit curve with 51 logarithmically spaced samples for values of α = 10 −5 . . . 10 5 .The misfit curve for this initial reference is displayed in Fig. 1 by green dots.For large values of α (bottom-right side of the plot), the distance to the reference model vanishes and the data misfit approaches the value d ≈ 13 (off-scale in Fig. 1) as obtained for the WDZe-model.By gradually reducing α, the model parameters start to depart from the reference, and distance to the reference model m (α) increases but d (α) decreases.We meet a first trade-off point at m (25) ≈ 10 −4 .Figure 3 (left panel) displays the parameters for O 1 and τ 1 at this trade-off point.The parameters for harmonics of small amplitudes still stick at the reference values (black lines), while the parameters of harmonics of large amplitude start to deviate from the reference.Further decreasing α certainly decreases data misfit to a second trade-off point, at m (0.25) ≈ 0.01 in Fig. 1.Beyond this value, the misfit computed for the untapered data d (α) increases again, although the misfit of the tapered data d (α) in the objective function in Eq. ( 19) necessarily further decreases, as shown in Fig. 2. The model hence is only more appropriate for tapered than for untapered data, which clearly indicates overfitting.
Despite finding the smallest value of data misfit at the second trade-off point, we choose the solution for the first.The reason is demonstrated in the right diagram of Fig. 3.At the second trade-off point the small-amplitude harmonics scatter strongly, they are mainly driven by noise in the signal, and with α ≈ 0.25 the regularization is too weak already.

Reference grouping
The wave groups used by Melchior et al. (1996) serve us as an initial set.They used a time series of similar duration (12 years).In the following, we will specify the groups by the Darwin names of the harmonics with the largest amplitude within the reference group.A full list of the groups in the final model from the current study is given in Tables 1 and 2. To distinguish tidal parameters for different degrees, we append * and + to names of degree 3 and four, respectively.There is no standardized nomenclature of tides.Other notations exist as well in the literature.Tables S1-S3 in the supplement compare some of them.Where we had to introduce additional groups, we took a Darwin name from Bartels (1985).In very few cases, we had to introduce new names.Other than Melchior et al. (1996) we completely separate degree 3 from degree 2 tides within groups.Narrow groups, defined by Melchior et al. (1996) only to separate the degree 3 tides from their degree 2 neighbors, were merged into larger degree 3 groups.
Figure 3 (left) shows that the parameters for the largeamplitude harmonics O 1 , τ 1 , and O 1 * depart from the reference, while their small-amplitude neighbors stay with the reference due to the regularization.We update the degree three reference for the entire displayed frequency band to O 1 * uniformly.O 1 and τ 1 depart in opposite direction.For this reason, we split the frequency range into two reference groups, O 1 at the lower frequencies and τ 1 at the higher frequencies, each of them being adjusted to a new reference value, which is uniform within the reference group.Opposite departure of significant amplitude parameters from a uniform reference level always indicates that the data call for the reference group to be split in two or more subgroups.In a few cases, we subsequently join (merge) reference groups, if the new parameter updates indicate that all harmonics in the larger frequency range could well go with the same parameter value.This was applied to degree 2 L 2 and KNO 2 to constitute L 2 (see Table S2 in the supplement for KNO 2 ).All degree 4 harmonics form one single reference group from the very beginning.Most degree 3 tides were joined into larger reference groups since due to their small amplitudes the grouping scheme was too fine and most small tides were not significant.Tables 1 and 2 in Appendix B specify the groups in the final reference model.The group symbols used in the following text refer to these tables.Where a group spans a larger frequency band, we use parentheses like (γ 2 , M 2 , δ 2 ), which refers to the entire M 2 group in Table 5 presented by Melchior et al. (1996).The symbol M 2 in bold font indicates the largest amplitude harmonic in the group, the harmonic which usually drives the update of the reference level.

Groups treated specifically
From the very beginning, we split tidal group (γ 2 , M 2 , δ 2 ) into three: (γ 2 , α 2 ), (ω 2 , M 2 ) and (β 2 , δ 2 ) in order to capture the annual modulation discussed by Schroth et al. (2018).As the iteration of the model update proceeds, it turns out that the data call for different levels of all groups from γ 2 to δ 2 .Therefore, we decided to split two M 2 satellites from their previous groups: γ 2 (from α 2 ) and δ 2 (from β 2 ).Moreover, we split (K 1 , κ 1 ) from the beginning, where κ 1 is the major amplitude nodal satellite to K 1 with Doodson number 165.565.S 1 contains the frequency of 1 cpd at which the recording contains gravity signals driven by solar radiation rather than tidal forces (the so-called radiation tides).For this reason, we put 164.556 and 164.554 into the separate group S 1 .For the reference parameters of their small-amplitude neighbors 164.544 and 164.566 (3% of already small S 1 ), we adjust to the average of the parameters of P 1 and K 1 .
In cases where the parameters of the major amplitude harmonic in a group remain at the reference value, we consider this signal not resolved by the data.In such cases, we merge smaller groups to larger groups.This primarily took place for parameters of degree 3 harmonics, and for degree 4, we end up with one single group M 4 + covering the entire tidal frequency band (all species).
For the very same reason, we initially merged (R 2 , K 2 ), because the parameters for R 2 remained at the level of the WDZe-model.The subsequent inversion, however, clearly indicated that the parameter value previously set for K 2 is not appropriate for R 2 .The parameter for the latter clearly departed in direction of the WDZe-model.We then kept R 2 and K 2 separate in subsequent iterations.
In the last iteration, we identified four additional groups that might be split: ζ 1 (from SGQ 1 , see Table S1 in the supplement and Table 1), o 1 (from Q 1 ), ι 1 (from J 1 ) and ω 2 (from M 2 ).The new groups do not have established names.They are treated in a special way in the sense, that after splitting these groups, we did not reset the factors to the initial values, and thus, this action should rather be considered as the 'final adjustment.'The reason for this is the more conservative assumption due to their medium amplitude (2-10 times smaller than harmonics in very close frequency vicinity, differing by fourth digit in the Doodson number, except ζ 1 ).These harmonics could potentially have their estimates Fig. 4 Gravimetric factor (middle) and phase (bottom) of the BF21-model.The dots are: the initial reference model (WDZe, small black), the degree 2 parameter estimate (blue), the degree 3 parameter estimate (red), the degree 4 parameter estimate (green).The phase of S 1 is almost 10 • (Table 3) and therefore is off-scale.Estimates of degree 4 parameter in the quarter-diurnal band are not displayed.The values of all other groups are within the axes-range of the diagrams.More detailed plots, showing all species separately, are in Fig. 9.The top panel shows the logarithm of catalog amplitude after multiplication with 1 × 10 10 s 2 nm −2 affected by noise and showed 'unclear deviations' from the reference in the previous iterations.

Termination of the iteration
At each iteration, we investigate the misfit curve, the departure of model parameters from the reference as α decreases, and the Fourier amplitude spectrum of the data residual.As we proceed, the necessary updates of the reference model become smaller and the first trade-off point of the misfit curve becomes weaker and finally almost vanishes as shown in Fig. 1.In the first iterations, we clearly see signals at tidal frequencies vanishing from the residual.As the iterative procedure progresses, the reduction in data misfit becomes smaller as shown in Fig. 1, while some narrow-band signals remain, which apparently are not of tidal origin.We terminate the iteration when no significant reduction in data misfit is achieved and the model parameters have converged to a local reference model.We needed 5 iterations (i.e., inversion runs) including the initial iteration to obtain (and confirm) the final model

The final model BF21
The final tidal model, which we call BF21, is specified by the numerical values of its parameters in Appendix B (Tables 3  and 4).A graphical display is shown in Figs. 4, 5, and 6.Estimates for confidence intervals are discussed in Appendix C 'Precision and accuracy' and are displayed in Fig. 9.
The dots in Figs. 4, 5, and 6 are displayed for harmonics in diurnal, semi-diurnal, and ter-diurnal bands in Tamura (1987) catalog.The factors x l and y l are adjusted to be uniform within each of the reference groups specified in Tables 1 and  2. Gravimetric factor and phase are computed by Eq. ( 7) from the final model parameters Figure 4 gives an overview of the entire model in terms of gravimetric factor (left) and phase (right).The parameters of the WDZe-model are displayed by black dots for comparison.Due to the high-pass filter, no signal is available at the frequencies below 0.5 cpd.Hence, the parameters there stick to the WDZe-model and for this reason are not shown.In the diurnal band, the very narrow FCN resonance dominates the degree 2 response (blue dots) in gravimetric factor and correspondingly in phase.In the semi-diurnal band, a kind of resonant character is seen as well in the gravimetric factor and the corresponding phase lead of degree 2. The resonance is less sharp and is presumably due to the signal contributed by ocean loading.The frequency dependence of the response at degree 3 is much weaker, with no clear resonance being present.For the degree 4 response, the model is uniform at all frequencies by definition of the reference groups.Neither the degree 3 nor the degree 4 parameters follow the frequency dependence of the degree two parameters.Further, the degree 3 and the degree 4 parameters depart largely from the reference level in different directions.A common, uniform ratio between parameters of different degrees obviously is not appropriate.
As an example for the update of the ratio between degree two and degree 3 admittance, we take the gravimetric factors for N 2 , N 2 *, L 2 , L 2 *, their update with respect to the WDZemodel and their 2 σ confidence intervals from Tables 3 and 4. Gravimetric factors are determined to better than 0.2% (2 σ level).While the gravimetric factors for N 2 * and L 2 * are slightly reduced (by less than 0.5%) with respect to WDZe, the factors for N 2 and L 2 are significantly increased (almost 3% for N 2 and about 1.4% for L 2 ).The ratio between degree 2 and degree 3 admittance hence changed by more than 10 σ in both cases.Both models, BF21 and WDZe, will thus produce different beating patterns in the tidal signal, where only one of them can match the data.This, we expect, is the cause of the periodic variation of parameters in the MWA, as discussed by Schroth et al. (2018) and Schroth (2020).Schroth (2020) discusses the temporal variation of the parameters for the M 2 group and the hypothesis that they are caused by ocean loading in particular.Müller et al. (2014) investigate apparent annual modulation of the oceanic M 2 -response by adjusting the satellites α 2 and β 2 in their time-invariant oceanic tidal model.Likewise we seek for a finer resolution of the model in the M 2 group and split the band from 1.923 cpd to 1.943 cpd into six groups γ 2 , α 2 , ω 2 , M 2 , β 2 , and δ 2 (Table 3).Figure 5 provides diagrams for the final model parameters focused on this frequency range.During the iterative procedure, it turned out that the param-eters of all four satellite groups γ 2 , α 2 , β 2 , and δ 2 should be updated with different values.This pattern produces annual as well as semi-annual modulation of the signal.
The variation of the gravimetric factor with frequency (Fig. 4) in the diurnal band is dominated by the FCN resonance.This resonance is specified in the WDZe-model based on precise results from VLBI observations.The tidal analysis just applies an adjustment to this resonance as shown in Fig. 6.The strongest deviation from the WDZe-parameters appears in the gravimetric factor and phase (off-scale in Fig. 4 right) of S 1 .This is due to non-tidal (thermal) forcing of masses in the atmosphere, and the parameters do not represent the tidal response of Earth's body.For the other major-amplitude groups in the diurnal band there is little frequency dependence seen in the adjustment, which indicates that the FCN resonance in the WDZe-parameters is appropriate.A rather uniform factor is applied, which appears different for the degree 2 and the degree 3 response.
Amplitude factors for degree 3 tides were investigated by Ducarme (2012).A collection of M 1 *, N 2 *, L 2 * and M 3 * (there: M 1 , 3MK 2 , 3MO 2 and M 3 ) estimates for 16 SG stations is presented for comparison with different Earth models.The estimates of these harmonics in BF21 align with his results from Strasbourg, Bad Homburg, and Wettzell.The detailed model of degree 3 tides for Moxa presented there is very similar to the BF21 model in terms of the observed response monotonicity.Even though most of the tidal parameters from BF21 model differ more than confidence intervals, in these Central European stations the ratio between tides of degree 2 and 3 changes in the same way along frequency, which is visible as a systematic trend.All these estimates, including BF21 confirm that the a priori ratio in all the studied observatories differ from the theoretical Earth models, which should cause apparent modulation in the MWA procedure.The cause for the ratio of degree 2 to degree 3 admittance being different from the prediction in the WDZemodel most probably is the contribution of ocean loading, which by definition is not considered in the WDZe-model.Munk and Cartwright (1966) already pointed to the need to separate degree 2 admittance from degree 3 admittance for ocean tides.Sulzbach et al. (2022) present a recent study of degree 3 ocean tide models.
Accuracy and precision of the BF21-model parameters are discussed in Appendix C.

Implementation
We expect that the BF21-model partly captures the timeinvariant signal modulation (beat) that is reported by Schroth et al. (2018).To investigate the remaining components in the apparent admittance of the Earth (including oceans and atmosphere), we run a moving window analysis (MWA) like Schroth et al. (2018) do.We subdivide the total time series into segments of 90 days length, which overlap such that the moving window progresses in intervals of 1 day.The full length time series are prepared in advance by decimation and filtering the recordings as well as C rig kl and S rig kl as defined in Eq. ( 8).With the parameters of the BF21-model, we compute the samples (index k) of the initial synthetic tidal signals of all harmonics (index l) and its quadrature as defined in Eq. ( 11).C BF21 kl and S BF21 kl are the elements of matrices C C C BF21 and S S S BF21 , respectively.
In order to compute results which are directly comparable to those presented by Schroth et al. (2018), we use a priori grouping and dismiss regularization by setting α = 0. We take the sum over the synthetic signals for all harmonics within each wave group in matrices C C C BF21 and S S S BF21 .The matrices C C C ini and S S S ini , which result from this summation, contain only 12 columns each.The wave groups are defined in Table 5.
Within each time window of the MWA, we compute factors x l and y l with respect to the elements C ini kl and S ini kl of matrices C C C ini and S S S ini , respectively, by minimizing the objective function in Eq. ( 13) without model constraint, i.e., α = 0. Thus, we obtain time depending tidal parameters with a value of gravimetric factor where and The index l specifies the wave group, not the time window.We run the very same analysis a second time, but with the WDZe-model as a reference.In Eqs. ( 34)-( 36), BF21 is replaced by WDZe then.The comparison of the obtained MWA time series for δ l and l shows in how far the BF21model appropriately adjusts the beating in the tidal signal.

The MWA results
Figures 7 and 8 show the MWA results for selected wave groups.We overlay the values for the analysis which uses the WDZe-model (in green) for comparison by the results shown by Schroth et al. (2018) (in black) to demonstrate the consistency with their results.Small differences between the green and the black curves are due to the different handling of gaps in the raw data.
For all parameters shown in Fig. 7, the amplitude of temporal variation is significantly reduced when running the MWA with respect to the BF21-model rather than the WDZemodel.While the results with respect to the WDZe-model have clearly periodic character, the results with respect to the BF21-model almost do not show this.For K 1 group, the annual variation was captured by the more appropriate ratio between K 1 , P 1 , and S 1 in the BF21-model.A detailed analysis showed that the K 1 /S 1 is the more important in this case.Some periodic variations remain, but with reduced amplitude, and they are less systematic.We suppose they are caused by thermally driven signal components at the S 1 frequency.Schroth et al. (2018, their Table 4) discuss possible causes of periodic variations in the MWA results, based on the frequency distance of major tidal contributions.For K 1 , they discuss the consequences of the group parameter being not appropriate for the signal at S 1 .For μ 2 (their 2N 2 ), they discuss the satellites at a distance of 1/0.56 per year (μ 2 , the larger amplitude variational tide from M 2 , Bartels 1985).The parameters for these contributions are adjusted separately in the BF21-model, and the semi-annual variation consequently is reduced in Fig. 7c and d.Similarly variations at a period of 0.56 years are expected for N 2 (ν 2 , the larger amplitude evectional tide from M 2 ) and for L 2 (λ 2 , the amplitude smaller evectional tide from M 2 ).Both are adjusted independently from the rest of the groups in the BF21-model, which reduces the semi-annual variations in Fig. 7e-h.The groups N 2 and L 2 contain significant contributions at degree 3, which appear as satellites to the main degree 2 harmonic at a distance of 1/8.8 per year (N 2 * and L 2 *, respectively, due to the revolution of the Moon's perigee).A modulation with the corresponding period is obvious in Fig. 7e-h for the parameters computed with respect to the WDZe-model.The BF21-model uses independent parameters for the degree 3 tides, adjusted to the observations.Consequently, the variations with 8.8 years period are gone, in the MWA with respect to the BF21-model.
The BF21-model in the frequency bands of the MWA groups as listed in Appendix D (Table 5) is in fact more detailed than what is discussed at the above paragraphs (cf.Tables 3 and 4).In how far further details contribute to the Fig. 7 The results of moving window analysis for tidal groups K 1 , μ 2 , L 2 , and N 2 (top to bottom) for gravimetric factor (left) and phase (right).Tidal parameters correspond to the main harmonic in the group, which gives the group its name.Colors: green: MWA with respect to the WDZe-model, orange: MWA with respect to the BF21-model, black: results by Schroth et al. (2018) for comparison Fig. 8 The results of moving window analysis for tidal group M 2 .Top: gravimetric factor (left) and phase (right).Tidal parameters correspond to the main harmonic in the group, which gives the group its name.Colors are: green: MWA with respect to the WDZe-model, orange: MWA with respect to the BF21-model, black: results by Schroth et al. (2018) for comparison.Bottom: difference between results with respect to BF21-model and with respect to WDZe-model changes between the MWA results with respect to the WDZemodel on the one hand and the BF21-model on the other hand is not investigated here.
For the M 2 group (Fig. 8a and b), the differences between the two results are not as clear.However, the variation amplitude is clearly reduced in the MWA with respect to the BF21-model.The variations of annual periodicity are captured by the fine-grained structure presented by α 2 , β 2 , M 2 , γ 2 , and δ 2 .This becomes obvious, when computing the difference between both results, as shown in Fig. 8c and d.Significant variations of rather random nature remain in the case of M 2 .The cause for the frequency dependence of tidal parameters within the M 2 group should not be searched in the Earth body.It might be due to a resonance in the ocean's admittance to tidal forcing within the M 2 band, which is part of the loading signal.

Conclusions
We implement and demonstrate a robust harmonic tidal analysis method, which uses a constraint to a reference model in order to regularize the regression problem.It is robust in the sense that parameters for each harmonic in the tidal catalog can be inferred without encountering a singularity.No a priori wave grouping is needed.The solution is found by updating the reference model, where parameters for harmonics deviate from the reference for specific harmonics.This update is done uniformly for all harmonics in a finite frequency band in order to account for the credo of smoothness.Different adjustments are only applied, where results of the linear regression problem indicate that harmonics significantly require different parameters to reduce the data misfit.In that sense grouping is applied, but to the reference model only, not to the parameters in the linear regression.This grouping is not set a priori or by an automatic algorithm.It is the investigators choice based on data constraints seen in the inversion and in that sense is data driven.The obtained spectral resolution is not limited by the time series length and signal-to-noise ratio and can provide the super-resolution of tides as suggested by Munk and Hasselmann (1964).
We apply this approach to 11.5 years of gravity data recorded with the superconducting gravimeter SG056 at Black Forest Observatory (BFO, Schiltach).With this analysis, we demonstrate that the robust approach can exploit the super-resolution of tides.According to the Rayleigh criterion, the fundamental frequency would be 1/11.5 cpy, and the analysis indicates the mean noise level at ≈ 0.88 nm s −1 .This tidal model captures the local response of the Earth as a whole (including oceans and atmosphere).We separate the response parameters for tides of degrees two, three, and four and demonstrate that the optimal model differs from the ratio between degree 2 and three body tides as hard-coded in programs like Eterna 3.40, in terms of WDZe.For L 2 and N 2 , this change is larger than 10σ and has strong consequences in the MWA, where the periodic variation of parameters is strongly reduced.Further, we demonstrate that the constrained approach allows us to identify variations of tidal parameters in rather narrow frequency bands, which would be lost in a priori grouping schemes.Then, again we find that groups (like KNO 2 ) being kept separate in other studies (e.g., Melchior et al. 1996;Calvo 2015) do not receive specific constraints from the data and should be kept together with their major neighbors (L 2 in the case of KNO 2 ).
The final model for sensor 1 of SG056, which we call BF21, has a very fine structure.In total, we identify 46 degree 2 groups, 14 (purely) degree 3 groups, and 1 degree 4 group, as given in Tables 1 and 2. This goes beyond what is possible with traditional implementations, which only allow for wave grouping being done along frequency.All parameters estimated by analysis of the 11.5 years of gravity for both spheres of SG056 are listed in Tables 3 and 4. The estimate for the air pressure (AIR) factor R BF21 is also displayed there.
The too coarse frequency resolution in traditional, unconstrained tidal analysis as reported, for example, by Schroth et al. (2018) was suspected to be the cause of some of the temporal variations found for tidal parameters.The results presented in the current contribution corroborate this hypothesis.We use the local response model BF21 as the a priori model in an unconstrained moving window analysis (MWA).In that, we replace the WDZe-model, which by definition is purely theoretical and only accounts for the body response, by the BF21-model, which additionally contains the timeinvariant ocean-loading response and atmospheric signals.
The MWA uses the same wave grouping as applied by Schroth et al. (2018).We demonstrate that the quasi-periodic temporal variations for some wave groups are reduced by up to a factor of seven in amplitude.For groups M 1 , K 1 , μ 2 , N 2 , L 2 , and S 2 , periodic temporal variations are clearly reduced or even vanish.For K 1 , it is not only essential to keep the radiation tides (S 1 ) separate, but also to allow for a small adjustment of the frequency dependence within the NDFW resonance.The latter most likely is not due to the FCN model used in the a priori body tide model being wrong, but due to loading effects taking place in this frequency band.Only part of the variations previously seen for M 2 are captured by the local response model as well.This leaves room for the remaining variations of the parameters for M 2 being caused by time varying ocean loading, as discussed by Schroth (2020).Some of the constituents (like α 2 , β 2 , or S 1 ), which are responsible for the temporal variation of larger groups, turn out to be not constant themselves.In the BF21-model, they do not exactly capture the time-invariant tidal response of the Earth, but rather the specific response (in case of α 2 and β 2 ) or a specific non-tidal signal contribution (in case of S 1 ) in the analyzed time window.
This way our analysis clearly shows that in a search for time-variant components of Earth's admittance, a detailed local reference model must be used as a basis to which adjustments are applied.To find this detailed model of the local response, regularization and data-driven grouping are needed.Approaches based on a priori grouping schemes, like based on the Rayleigh criterion, likely miss essential properties of the admittance function.Response properties previously described as time dependent are now described as frequency dependent, which is more appropriate.This does not imply any conclusion as to their physical causes.However, this is an essential step, for example, in the direction to identify ocean-loading signals in on-shore gravity records.Tidal ocean models become increasingly sophisticated.Sulzbach et al. (2022), for example, investigate the difference in tidal admittance for degree 3 and degree 2 in ocean tides.The results of the current study indeed point in the direction of ocean loading being different at degree 2 and degree 3. The results shown by Schroth (2020), however, indicate that current ocean models might not yet be able to correct for ocean-loading signals beforehand.
The proposed robust approach to tidal analysis allows us to exploit the super-resolution of tides.It supports a truly data-driven wave grouping and could in principle be further developed to infer a continuous and smooth frequencydependent model.Moreover, with this approach it is possible to deliberately test models of Earth's admittance by setting a specific reference model and testing whether the data constraint drives the parameters of specific harmonics away from the reference.The trade-off curve (Fig. 1), which uses the misfit for untapered data, effectively helps to avoid overfitting.The model constraint causes a general tendency toward the reference model, which is unfavorable.For this reason, deliberate testing of parameters is an essential part of the iterative procedure.Ciesielski (2023) applies RATA to five European SG stations and another eight SG stations distributed all around the globe.The length of these time series vary from 4 years to 22.5 years and they come with different signal-to-noise ratios.For all of them, Ciesielski (2023) demonstrates the exploitability of super-resolution and the reduction of periodic variations in the MWA results, if the improved local response model is used.
Adjustment factors x and y (and r ), regression factors applied to the a priori model that are used to adjust the implemented complex admittance factors in a priori model to the data.kl and S ini kl of matrices C C C ini and S S S ini , theoretical tidal forcing assumed for Earth at particular location and time, rigid Earth tide Group merge, the reference groups that used to have different common pairs of factors has one common pair of factors applied to account for parameters that have no significance in model constraint imposed by the observed signal Group split, the reference group that used to have one common pairs of factors is divided into two or more subgroups, to account for the conflicting constraints imposed by the observed signal Initial model, the reference model with which the iterative procedure starts, here it is local WDZe model, derived from PREM.In principle it can be any reasonable model.In phase and quadrature, harmonic signals that are offset in phase by one-quarter cycle (π/2), two components to which a harmonic signal with angle modulation can be decomposed into; its purpose is linearization of the problem, mapping amplitude, and phase on complex plane Local response model (time-invariant), The final solution of the iterative procedure, to which no further distance to the reference model is justified by data misfit decrease; it is expected to capture the significant features of the locally recorded signal Misfit curve, curve that displays distance to the reference model versus (untapered) data misfit for solutions under various trade-off parameter α Objective function E, here: also a loss function or cost function; expresses the terms that have to be minimized in the least-squares problem Regularization, imposing additional constraints that bias the solution to stabilize the inversion; is essential to estimate a usable solution to an otherwise intractable ill-posed or ill-conditioned inverse problem.Reference groups, the sets of harmonics to which the common reference model factors are applied Reference model, the reasonable model which is close to the expected solution, its difference to estimates is used as an additional penalty term in LS objective function Rigid Earth, theoretical Earth model, rigid oceanless Earth with no atmosphere that does not undergo deformation when tidal field is applied Taper, a window function by which signal (and regressors) are multiplied to avoid spectral leakage in time series analysis Time-invariant properties, signal or model properties which do not change in time; e.g., FCN resonance, expected solid Earth response to tidal forcing; Time-variant properties, that may periodically repeat or are systematic, mainly refers to beat or particular phenomena in ocean tides or radiation tides; periodic time-variant model properties often can be expressed by a frequency dependence of time-invariant properties Trade-off parameter α, a number that defines quantitatively whether LS fit seeks more for the data misfit minimum (variance) or for an additional constraint (distance to the reference model) minimum (bias) Trade-off point, point on the misfit curve where its behavior significantly changes.Unconstrained inversion, least squares regression fit without regularization term, tidal harmonics are a priori grouped (A priori) wave groups, sums of tidal signals, usually in defined frequency bands.The a priori body response is applied and hence the common factors are estimated for tidal parameters of the group; the purpose of that parameterization is to stabilize the otherwise ill-posed problem error estimates are given in Tables 3 and 4 together with the model parameters.

C.1 Exploring constraints
Significant information on the constraints put by the gravity data on the model parameters is already obtained during the iterative procedure discussed above.Changes in model parameters are repeatedly tested against the regularization constraint.A deviation from the initial WDZe reference model will only be accepted, if the recorded data call for it.This deviation is represented by a deviation of factors x WDZe l and y WDZe l from 1 and 0, respectively.This procedure and its consequences for the wave grouping chosen in the final model are discussed above.

C.2 Results for the second gravity sensor
The difference between the analysis results for two colocated sensors may represent a lower limit for the confidence intervals.The SG056 is a dual-sphere instrument which operates two different gravity probes within the same sensor unit.Forbriger and Heck (2018) describe the essential properties of this instrument.The probe mass of the lower sensor (G1) is 17.7 g, and the upper sensor (G2) in the unit uses a sphere of 4.34 g.Data recorded by the former typically appear less noisy, although it more likely saturates during stronger ground motion.The parameters for the BF21-model (Tables 3  and 4) are computed for data from sensor G1.
We additionally analyze data from sensor G2 in the very same (model-regularized, iterative) way and for the same time window as described above for the analysis of G1 data.In Tables 3 and 4, the parameters for this BF21G2-model are referred to as δ 2 and 2 .Because both sensors operate in the same unit they practically experience the same environmental noise source (including gravity noise from the atmosphere and hydrosphere as well as some instrumental noise sources).The inferred tidal model would be expected to be identical, though it is not.The difference δ 2 − δ 1 and 2 − 1 between the model parameters of both models is a measure for uncertainty and bias being produced by trade-off in the inversion as well as different instrumental noise and gap handling.The accuracy of the model parameters must not be assumed to be better than these differences, which are given in Tables 3 and  4 as well.

C.3 Unconstrained inversion
Using the groups defined in Tables 1 and 2 for an a priori grouping, we run an unconstrained inversion for the data of sensor G1.The values for gravimetric factor δ u and phase u for this BF21u-model together with confidence intervals

Fig. 1 Fig. 2
Fig. 1 Misfit curves obtained during the inversion for the stationary model BF21.Data misfit d (α) is computed for the untapered data by Eq. (32).Distance m (α) to the reference model is computed by Eq. (30).Both are a function of α.The curves for three stages in the iterative procedure are displayed: initial reference model (green), first iteration (orange), final model (black)

Fig. 3
Fig. 3 Gravimetric factors for the harmonics in the groups O 1 and τ 1 .The circles refer to: the degree 2 parameter estimate (blue), the degree 3 parameter estimate (red), the actual new reference model (small orange; top), and the hypothetical new reference model (small orange; bottom) for the next iteration.The diameter of circles corresponds to the logarithm of the harmonic forcing amplitude.Because some of the minor-amplitude harmonics are obscured by larger dots, each of the harmonics is indicated by a white dot overlaid to the larger blue ones.The

Fig. 5 Fig. 6
Fig.5Gravimetric factor (middle) and phase (bottom) of the BF21-model in M 2 tidal group.The dots are: the initial reference WDZe-model (small black), the degree 2 parameter estimate (blue), the degree 3 parameter estimate (red).The top panel shows the logarithm of catalog amplitude after multiplication with 1 × 10 10 s 2 nm −2 Complex admittance factors X and Y , the a priori model factors, the Earth body response that given theoretical model accounts for Constrained inversion, least squares regression fit with regularization term, harmonics are treated as separate regressors Data error σ , the expected value of the mean noise level (RMS) Data misfit d , a normalized residual, a term in the objective function that expresses how the estimated synthetic time series deviates from the recorded data; it refers to the tapered forward operator and data Data misfit d , a data misfit that results from applying an estimated model to the untapered time series and untapered forcing operator Distance m to the reference model, a term in the objective function that expresses how the estimated model deviates from the reference model.Forcing signals C rig l (t k ) and S rig l (t k ), that are elements C rig kl and S rig kl of matrices C C C rig and S S S rig , theoretical rigid Earth tidal forcing assumed for Earth at particular location and time Forcing signals C ini l (t k ) and S ini l (t k ), that are elements C ini

Fig. 9
Fig. 9 Parameters for the model BF21u as obtained from unconstrained linear regression after a priori grouping (see text).Their confidence intervals at the 1σ -level are displayed as error bars.The horizontal bars specify the frequency interval of the respective group.They are

Table 1
RATA Regularized approach to tidal analysis WDZe purely elastic version of body tide Earth model by Wahr-Dehant-Zschau in which phase shift is neglected but relaxation of elastic parameters is included Definition of degree 2 wave groups in the final BF21-model.See Table 2 for a definition of column heads.In fact, only two harmonics form S 1 group, 164.556 and 164.554.Number of harmonics in tidal catalog constituting a reference group; Lead acc: RMS of acceleration of the most significant harmonic; Gr acc: RMS of acceleration of the total signal of the group; Lead Amp: Catalog tidal potential V amplitude of the most significant harmonic; Gr Amp: Catalog tidal potential V amplitude of the total signal of the group; Min/Max freq: The smallest and the highest frequency in a group, respectively; and Doodson: Doodson number of the most significant harmonic.Logarithms are taken to the base of 10.Horizontal lines indicate boundaries of MWA groups as specified in Table 5, where commonly M 4 + constitutes one group together with 3cpd tides

Table 3
Degree 2 parameters for the BF21-model