Parameters of the Earth’s Free Core Nutation from Diurnal Strain Tides

Earth deformation at the diurnal tidal frequencies includes the resonant tidal-forcing response caused by the Free Core Nutation (FCN), a retrograde mode related to the slight misalignment of the rotation axes of the outer core and mantle. We analyse data from four underground high-sensitivity laser extensometers, whose signal-to-noise ratio in the diurnal tidal band is particularly high, and provide an alternative independent estimate of the FCN complex frequency with respect to more usual techniques (nutation and gravity). Firstly, we differentiate displacements due to diurnal solid tides to obtain extension along any azimuthal direction in terms of three complex parameters (A, S, C) which depend on latitude and frequency. Then, we demonstrate that we can invert the FCN complex frequency and the sensitivity of Im(A) and Re(S) to the resonance from our data. Lastly we obtain the probability distributions of those four parameters. Our results are in full agreement with those from nutation and gravity, as well as with reference IERS (International Earth Rotation and Reference Systems Service) values. Sensitivities of Im(A) and Re(S) to the resonance are estimated here for the first time and are in agreement with values computed using reference Love and Shida numbers from IERS.

The Free Core Nutation (FCN) is a retrograde mode related to the slight misalignment of the rotation axes of the Earh's fluid outer core and mantle, with a period T FCN ∼ 430 sidereal days in the celestial reference frame. It is influenced by the flattening of the core-mantle boundary and by any deformations on the surface of that boundary. These deformations can be induced by various mechanisms: among the most relevant, there are the fluid pressure on the outer core-mantle boundary (CMB) and the electromagnetic torque due to the differential rotation of the outer core and mantle in the presence of magnetic fields and visco-magnetic friction at the CMB, which leads to a partly dissipative coupling. The FCN causes resonance effects on the forced nutations of the Earth's figure axis 1 .
In the Earth-fixed reference frame, the FCN complex frequency is f FCN = (1 + 1/T FCN ) (1 + i/(2Q FCN )) cycles per sidereal day (cpsd), where Q FCN is the quality factor of the resonance. Since f FCN falls inside the diurnal tidal band, the FCN also causes a resonant response of some diurnal tides (e. g. P 1 , K 1 , Ψ 1 and Φ 1 ) to the tide-generating forces. Gravity, displacement, strain and tilt tides can be expressed in terms of the Love (h) and Shida (l) numbers, which characterize the Earth response to the tide-generating forces, so that the FCN can be described through the resonant behaviour of h and l in the diurnal band (for a general summary, see, e.g. 2 ). Diurnal Love and Shida numbers are usually given as a function of the tidal excitation frequency σ by the resonance formula (e.g. 3 )  . A fourth free rotational frequency of the Earth (the inner core wobble 5 ) is not in Eq. 1 because its resonant frequency (about 0.0004 cpsd) is very far from the diurnal tidal band and its effects on diurnal tides are negligible. Love and Shida numbers depend on latitude, because rotation and ellipticity couple the spheroidal deformation of a given degree to spherical harmonics of other degrees, and spheroidal to toroidal modes of deformation.
Investigation of the resonant modifications of diurnal tides and forced nutations sheds light on the FCN, its related resonance parameters, the acting physical mechanisms, the relevant Earth geophysical properties, and ultimately the Earth interior. Many different research and application fields (e.g., the Earth's orientation in the space and global reference systems) would benefit of an accurate determination of the FCN parameters. FCN University of Salerno, Department of Physics, Fisciano, 84084, Italy. ✉ e-mail: antonella.amoruso@sa.infn.it open effects on gravity tides (measured by gravimeters) and forced nutations (observed by very long baseline interferometry, VLBI) have been used since decades to estimate the FCN complex frequency (see e.g. 4,[6][7][8][9][10][11]. VLBI observations also provide celestial pole offset (CPO) series (i. e., the difference between the observed celestial pole position and the one obtained from conventional precession and nutation models 3 ) which include a FCN term. Empirical FCN models give temporal variations of the FCN term amplitude and phase, mainly related to geomagnetic jerks and atmospheric and oceanic excitations (e.g. [12][13][14][15][16]; temporal variations (few sidereal days) of T FCN have also been suggested (e.g. 17,18 ). Some recent estimates of T FCN and Q FCN are listed in Table 1.
Few observations of the FCN have been made using borehole water level, tilt and/or strain [19][20][21][22][23][24][25][26] . In this paper we use extension data (fractional change of instrumental length, i.e., linear strain) from two pairs of underground high-sensitivity laser strainmeters (interferometers, Fig. 1a, see Data and Methods for details); strainmeters and extensometers are often used as synonyms in the literature. Two strainmeters (BA and BC in what follows) were located inside the Gran Sasso underground laboratories (LNGS, Laboratori Nazionali del Gran Sasso, Italy) 27 . The other two strainmeters (LAB780 and GAL16 in what follows) are located close to the Canfranc underground laboratory (LSC, Laboratorio Subterràneo de Canfranc, Spain) 28 . Generally speaking, extension at the tidal frequencies includes the response of the solid Earth to the forcing tidal potential (solid tides) and to the water load oscillations caused by ocean tides (ocean loading), their sum giving tidal strain, as well as the response to environmental parameter changes (e. g., temperature and pressure) and other processes (e.g., hydrological processes like rain and snow melting). The interferometers measure local (tens-of-meter scale) strain directly, but local strain is biased by siting effects (e. g., cavity effects due to tunnel installation, surface topography, and rock inhomogeneities) 29,30 . Coupling between local extension ε (measured by the interferometers) and regional deformation ε ij can be effectively described by three coefficients α, β, and γ per interferometer so that where the x-axis is directed along the interferometer, the y-axis is perpendicular to that, and ε xy is shear strain 31 . Extension along any azimuthal direction η is: where ε θθ , ε φφ , ε θφ are the horizontal components of the strain tensor in spherical coordinates (r, θ, φ), unit vectors θˆ and φ being directed southward and eastward respectively. Complex amplitudes of diurnal strain solid tides (phase being given relative to the local tidal potential) are (see Data where H f is the amplitude of the tidal term of frequency f and A is the mean Earth radius. Complex parameters A, S, and C (A being related to areal strain, S and C being related to the two shear strain components) are expressed in terms of the latitude-dependent Love and Shida numbers and follow the same resonance formula (Eq. 1), thus: σ σ σ − ≈ − ) and focus on the FCN contribution to A, S and C, i.e. A 2 , S 2 and C 2 . The experimental data set consists of amplitudes and phases (i.e., complex amplitudes) of Q 1 , O 1 , P 1 , K 1 , 1 Ψ , 1 Φ , N 2 and M 2 at each interferometer, as retrieved from the tidal analysis of strain records. Although the semidiurnal tides N 2 and M 2 are not affected by the FCN resonance, they have been considered to help evaluating siting effects on deformation. Optimal model parameters are obtained through the comparison between the experimental data set and theoretical diurnal tides (Eqs. 3 and 4) corrected for ocean loading and siting effects (Eq. 2). The comparison is quantified both through the sum of the magnitudes of the differences between experimental and theoretical tides (L 1 -norm) and through the sum of their squared magnitudes (L 2 -norm). Optimization through the L 2 -norm is more usual (e.g. 11 ) but also much more affected by outliers than optimization through the L 1 -norm (e.g. 24 ). Because of that, as a drawback, the use of L 1 -norm generally produces wider probability density functions (PDFs) of the inverted model parameters. We estimate PDFs and possible correlations of inverted parameters following a frequentist approach based on the inversion of thousands of synthetic data sets which we obtain from the experimental one by adding a gaussian random noise whose standard deviation is given by the tidal analysis of the first-difference strain sequences (see Data and Methods).
In what follows we use the expression "reference value" to indicate the value of each parameter if computed by using reference IERS 3 frequency-dependent Love and Shida numbers at 42.6° latitude, i. e. between LNGS and LSC (Fig. 1a).

Results
In principle, we would like to retrieve the following model parameters: and C 2 at each station; • α, β, and γ for each interferometer.
Unfortunately, such a goal is practically unreachable, not only because of the large number of parameters. We assume that A j , S j , and C j are the same for both stations since the latitudes of LNGS and LSC are very close to each other (Fig. 1a). Moreover, it is intrinsically impossible to retrieve 1 , as well as α, β, and γ from strain tides simultaneously; thus, we fix A A /(1 ) to their reference values. A study of tidal strain sensitivity to A 2 , S 2 , C 2 indicates that the effects of Im(S 2 ), Re(C 2 ), and Im(C 2 ) are too small to be detectable and tidal strain is mainly sensitive to A Re( ) 2 (Fig. 2a). Firstly we carry out validation tests on a data set constituted by theoretical body tides to which we attribute the same uncertainties as the experimental tides, showing that the inversion technique produces reliable estimates of T FCN , Q FCN (actually, log 10 (Q FCN )), Re(A 2 ), A Im( ) 2 , S Re( ) 2 , and, for each interferometer, α, β, and γ (see Data and Methods and Fig. 2b,c). However, there is a strong trade-off (linear correlation) between Re(A 2 ) and T FCN , with a T FCN increase by 10 sidereal days when Re(A 2 ) decreases by about 5%. A similar behaviour is also observed when inverting the experimental data set; thus, it is difficult to estimate both Re(A 2 ) and T FCN without strong a priori constraints. Since PDF of Re(A 2 ) from the experimental data set agrees with its reference value, regardless of the ocean and Earth models used to compute ocean loading (see Data and Methods), we choose to fix Re(A 2 ) in the final inversions; such a choice does not affect PDFs and correlations related to the other parameters (Fig. 2d). No other trade-off related to T FCN is visible. The trade-off between Q FCN and A Im( ) 2 is limited to small values of Q FCN (Q FCN < 10 4 ), thus we leave A Im( ) 2 as a free parameter in the inversion. PDF of Q log ( ) is bimodal, exhibiting a much smaller secondary mode at large Q FCN values. This secondary mode is related to synthetic data sets where complex amplitudes of Ψ 1 and 1 Φ are very different from the starting (theoretical) ones; it lowers noticeably when uncertainties (and consequently widths of the random noise used to generate synthetic data sets) on Ψ 1 and 1 Φ are halved. Moreover, as discussed in 24 , the only appreciable effect of decreasing Q FCN from very large values down to 20,000 is a change of the Ψ 1 phase by less than 2°, i. e., much less than its uncertainty. As a consequence, median of Q FCN from the inversion of synthetic data sets (19,943) practically coincides with the value used to generate the theoretical data set (20,000), but the 95% confidence interval is 5,224 to 6 × 10 8 . When Re(A 2 ) is fixed to its reference value in the inversions, PDFs and correlations involving T FCN are narrower than when Re(A 2 ) is left free; all the other PDFs and correlations remain nearly unchanged (see Fig. 2d).
Secondly we invert the experimental data set, but, before showing related results, a couple of points have to be discussed.
(i) In principle, changes in environmental and hydrological parameters could affect our results. As for the LNGS interferometers (BA and BC), we already performed several attempts to correct strain for barometric pressure and air temperature, by comparing spectra of residuals from the joint tidal analysis of strain and air temperature or barometric pressure with those from the sole strain 24 . Noise level did not lower and we did not correct strain records for temperature and pressure fluctuations; indeed, the BA peak around S 1 (Fig. 3) is probably of anthropic origin. As for the LSC interferometers (LAB780 and GAL16), we already demonstrated that air temperature and barometric pressure effects on strain are actually negligible on time scales of one day or shorter 28,32 . (2020) 10:9756 | https://doi.org/10.1038/s41598-020-66426-7 www.nature.com/scientificreports www.nature.com/scientificreports/ However, strain records show clear signatures of late-spring snow melting, characterized by daily deformation cycles lasting few weeks, and heavy rain episodes, characterized by several-hour-long deformation pulses (manuscript in preparation, but see 33 for a preliminary study). Removing the affected time periods from both LAB780 and GAL16 data approximately halves the noise amplitude in the diurnal band (red lines in Fig. 3), but decreases the number of analysed data by almost 10% and introduces several additional breaks into the strain sequences; as a consequence, it is difficult to foresee whether results improve or do not. Those cut strain records are hereinafter referred to as "cleaned" data.
(ii) Ocean loading mainly affects semidiurnal tides at LSC (Fig. 1b), while giving a small even though not negligible contribution to the other tides. Although computed ocean loading depends on the ocean and Earth models that are used in the computation (see Data and Methods), results obtained with tested models are always within few percent in amplitude and few degrees in phase. These small differences affects inversion misfits (Fig. 1c) but have a small influence on PDFs of model parameters (see later in the text). In what follows, ocean loading is www.nature.com/scientificreports www.nature.com/scientificreports/ computed using the combination of ocean and Earth models that gives the smallest misfit (i. e., FES2014b 34 and continental PREM 35,36 ) unless otherwise specified. Figure 4a shows PDFs and correlations of the free model parameters when we invert the experimental data set obtained from the analysis of all strain data. All L 1 -norm PDFs and correlations are fully consistent with, although wider than, those obtained when inverting the data set based on theoretical body tides (Fig. 2d); PDF of Q log ( ) is somewhat shifted toward smaller values. Siting coefficients are well constrained (Fig. 4b); for the first time K 1 (the largest diurnal strain tide) is also used to estimate them: since K 1 is heavily influenced by the FCN, we could do that just because the siting coefficient estimate is included into the tidal parameter inversion procedure. The L 2 -norm PDFs of S Re( ) 2 and A Im( ) 2 are shifted toward positive and negative values respectively, probably as a consequence of the strong sensitivity of the L 2 -norm inversion to the presence of outliers. This strong difference between L 1 -and L 2 -norm results does not occur in the validation tests, where synthetic data sets used for the PDFs are obtained by adding a gaussian noise to "exact" tidal complex amplitudes.
To see how measured tides agree with modelled tides as computed using the median values of the L 1 -norm parameter PDFs (Fig. 4a), we compare the observed and theoretical complex transfer functions (H FCN ) of the FCN for each interferometer (Fig. 5a). Observed H FCN is the ratio of observed solid tides (i. e., observed tides corrected for ocean loading) on computed solid tides without resonance, i. e., putting A S C 0 2 2 2 = = = in Eq. 4. Theoretical H FCN is the ratio of modelled solid tides (computed using the median values of the FCN parameter PDFs from the L 1 -norm inversion of tides obtained from the analysis of all strain data and corrected for siting effects) on computed solid tides without resonance. For both observed and theoretical H FCN , solid tides without resonance are corrected for siting effects (Eq. 2) using the medians of the PDFs of α, β, γ at each interferometer (Fig. 4b). The fit is very good for all the tidal harmonics, except 1 Ψ and Φ 1 , at all the interferometers. However, taking into account the large uncertainties in the Ψ 1 and 1 Φ ocean loading effects, which are not included in the error bars and might bias theoretical tidal strain, the agreement is still reasonable for all the interferometers but Im(H ) FCN for LAB780. The FCN effect on K 1 strain is large (up to 25% in amplitude and 7° in phase, depending on the interferometer) thus robustness of retrieved model parameters is also grounded on the major K 1 tide.
Correlation plots between theoretical and experimental tides confirm the aforementioned results (Fig. 5b). In this comparison, theoretical tides are computed using the medians of the FCN parameter PDFs (L 1 -norm inversion of tidal complex amplitudes obtained from the analysis of all strain data, Fig. 4a), corrected for ocean loading and siting effects. Again, siting effects are computed using the medians of the PDFs of α, β, γ at each interferometer (Fig. 4b).
Using different ocean and Earth models or analysing "cleaned" data produces practically identical L 1 -norm PDFs (Fig. 4c); siting coefficients vary by no more than few 10 −2 . We obtain final T FCN and Q log ( ) FCN 10 PDFs by stacking all the histograms in Fig. 4c; medians and 95% confidence intervals are listed in Table 1.

Discussion and conclusions
We have inverted the complex amplitudes of six diurnal (Q 1 , O 1 , P 1 , K 1 , Ψ 1 , 1 Φ ) and two semidiurnal (N 2 , M 2 ) strain tides per interferometer by using both L 1 -norm and L 2 -norm optimizations. The L 1 -norm optimization gives almost indistinguishable results when we obtain the experimental data set from all strain data or after removing the time periods affected by hydrological processes from LSC records ("cleaned" data), as well as using different ocean (FES214 and TPXO9) and Earth (cPREM, GB, IASP91) models (Fig. 4c). Results obtained using the L 2 -norm optimization are less stable. This behaviour is not surprising because non-stochastic ill-modelled processes (e.g. ocean loading, hydrologically-induced deformation, changes in air temperature and pressure) can contribute to the mismatch between theoretical and measured tides, thus affecting L 2 -norm optimization to a larger extent than L 1 -norm optimization. Figure 4a,c indicate that using the L 1 -norm effectively reduces the impact of poorly known contributions to the tidal signal on inversion results. Because of that, we consider L 1 -norm results more reliable than L 2 -norm ones, although L 1 -norm PDFs are wider.
Amplitudes of Ψ 1 and Φ 1 are very small, thus their relative uncertainties on amplitudes and absolute uncertainties on phases are much larger than for the other tides involved in the inversion. If uncertainties on 1 Ψ and Φ 1 are halved when we invert the data set constituted by theoretical body tides (blue lines and dots in Fig. 2b), then correlation between A Re( ) 2 and T FCN decreases, PDFs of T FCN , A Re( ) 2 and A Im( ) 2 become almost symmetrical, and the secondary mode of Q log ( ) FCN 10 PDF lowers noticeably; PDF of S Re( ) 2 does not change appreciably. Medians of all PDFs also do not change appreciably and are quite similar to the reference values used to compute the theoretical body tides. It is noteworthy that halving Ψ 1 and Φ 1 uncertainties has no visible effect on the PDF of T FCN when A Re( ) 2 is fixed (Fig. 2d). Because of the monochromatic nature of tides, halving uncertainties on experimental Ψ 1 and 1 Φ would require analysing strain records four times as long. PDF of T FCN in this work is about half as wide as in 24 mainly because of the inclusion of LSC data in the analysis, while the use of more recent ocean and Earth models and improvements in the inversion procedure have produced no significant change in the T FCN median.
As previously mentioned, gravity data have been extensively used to investigate FCN through its resonance effects on tides. Strain tidal measurements have a lower signal-to-noise ratio than gravity, but relative perturbations on strain tides are about 10 times larger than on gravity tides and inversions of synthetic data sets suggest that the resolving power of gravity and strain tides are comparable 24 . There is a correlation between T FCN and the real part of the resonance strength in gravity tides Re( ) 2 δ 10,11 : since Re( ) 2 δ strongly depends on h 2 2 , this correlation is consistent with that between T FCN and A Re( ) 2 in strain tides (Fig. 2b) www.nature.com/scientificreports www.nature.com/scientificreports/ We are conscious that there are some (currently) unavoidable inconsistencies in our approach. As for theoretical solid tides, they are computed using Love and Shida numbers based on a Spherical Non Rotating Elastic Isotropic (SNREI) Earth model (PREM) and modified to account for ellipticity and rotation 3 . Computations of solid tides using a full 3D Earth model with mantle heterogeneities evidence differences <0.1% in gravity and vertical displacement with respect to a radially symmetric Earth model 37 : these differences are small and eventually counterbalanced by siting effect corrections, which are estimated during the inversion procedure. However, because of the trade-off between A Re( ) 2 and T FCN , A Re( ) 2 is fixed to its reference value from parameters in 3 , which were obtained from VLBI data using geophysical models available at the end of the 20th century, and any change in A Re( ) 2 would imply a change in T FCN . As for ocean loading, strain Green functions from SNREI Earth models (GB, cPREM, IASP91) are used. A more correct approach would involve computing deformation using a full 3D Figure 2. Resolving power of tidal strain data. (a) Amplitude ratio and phase difference between solid tides computed along LAB780 and GAL16 using reference A 2 , S 2 , C 2 values and after increasing one by one A Re( ) 2  www.nature.com/scientificreports www.nature.com/scientificreports/ Earth model, which should be particularly accurate at short distances from the stations. This is an important point, even though Fig. 4c suggests that PDFs of T FCN and Q log ( ) are robust, presumably because of the small contribution of ocean loading to diurnal tides at LSC and LNGS.
In conclusion, data analysed in this work (strain tides) are quite different from those (forced nutations and gravity tides) generally used to estimate FCN parameters (period, quality factor, resonance strengths) and have a different sensitivity to some of them. The inversion approach is also different, being based on robust L 1 -norm misfit minimization. All the T FCN PDFs, which we obtain using different ocean and Earth models to compute the ocean loading contribution to tides, agree with reference values 3 and estimates from previous analyses of nutations and gravity (Table 1). Median values from GB and cPREM Earth models are somewhat smaller than reference T FCN , but IASP91 results are in full agreement with most of recent estimates from nutation. Although the final PDF is wider than those from the analyses of nutations, this work provides an alternative independent estimate of T FCN . Accurate estimates of Q FCN come from nutations only, nevertheless medians of the distributions in Fig. 4c are close to (actually, somewhat smaller than) its reference value 3

. Parameters
A Im( ) 2 and S Re( ) 2 are estimated in this work for the first time; their PDFs are in agreement with reference values computed using Love and Shida numbers from 3 .

Data and Methods tidal displacements and strain.
Here we give the diurnal and semidiurnal degree-2 displacement tides for an elliptical rotating Earth with an imperfectly elastic mantle as in 3 , Equations 7.1b and 7.1c. Differently from 3 , we use spherical coordinates (r, θ, ϕ), where the unit vectors θˆ and φ are directed southward and eastward respectively. We get the three components of the surface strain tensor (ε θθ , ε ϕϕ and ε θϕ , Eq. 4) in spherical coordinates after taking the spatial derivatives of displacements; the Earth surface is approximated as spherical (radius a) in the derivatives.
Because of the Earth's rotation and flattening, spheroidal displacements of degree 2 order 1 (diurnal degree-2 tides) are coupled to spheroidal displacements of degree 4 and toroidal displacements of degree 1 and 3. If diurnal tidal displacements are given with respect to the potential evaluated on a spherical surface, spheroidal displacements are expressed in terms of latitude-dependent Love ( θ h( )) and Shida (l( ) θ ) numbers, and toroidal displacements are expressed in terms of the additional numbers l (1) and l′, this latter acting as a no-net-Earth-rotation correction 38 .
For a diurnal tide of frequency f : www.nature.com/scientificreports www.nature.com/scientificreports/ where: After taking the spatial derivatives of displacements, we get Eq. 4. Parameters A, S, and C are:  3cos  1)  3  5  2  cos  1  2  ,   5  4  (3cos  1)  1  2  (cos  1) , Parameter l′ does not appear in the expressions for the strain because it refers to purely rotational displacements.
Parameters h (0) , h (2) , l (0) , l (1) , l (2) , and l′ follow the resonance formula given in Eq. 1; numerical values based on the preliminary reference Earth model (PREM 35 ) corrected for mantle anelasticity 39 are in 3 . FCN effects on A, S and C are represented by A 2 , S 2 and C 2 : from Table 7.1 in 3 , we see that the main contribution to the (complex) A 2 term comes from h 2 (0) , the main contribution to the (complex) S 2 term comes from l 2 (0) , and the main contribution to the (complex) C 2 term comes from l 2 (1) . We also give expressions for semidiurnal tides because we use them to evaluate siting effects on deformation. Spheroidal displacements of degree 2 order 2 (semidiurnal degree-2 tides) are coupled to spheroidal displacements of degree 4 and toroidal displacements of degree 3.
For a semidiurnal tide of frequency f: 3 ( )sin cos ( 2 ) Strain data. We use strain records from two pairs of high-sensitivity laser strainmeters (interferometers, Strain record analysis. At first, we low-pass filter and decimate (1 sample per 1800 s) strain data from each interferometer, thus generating evenly spaced time series with gaps. Then, we pre-whiten each ungapped data segment by generating first-difference sequences, because deformation noise is red (i.e., its power spectral density increases as frequency decreases) while most tidal analysis techniques assume a gaussian residual distribution (i.e., a white noise spectrum) when fitting tidal harmonics and, e. g., polynomials to the data time series. Pre-whitening makes fitting more reliable as it lowers low-frequency noise, which is the most difficult to deal with having a statistically non-null mean value; first-difference strain sequences exhibit a noise spectrum which is almost white, e.g. 41 .
To estimate signal-to-noise (S/N) in the diurnal tidal band for each interferometer, we calculate the amplitude spectrum of the gapped first-difference time series through an iterative deconvolution of the spectral window www.nature.com/scientificreports www.nature.com/scientificreports/ in the frequency domain, by using CLEAN code 42 (Fig. 3). Overall, spectra exhibit very small diurnal contamination, mainly as for LSC data, and very good S/N in the diurnal tidal band. Unfortunately, analogous detailed spectra are seldom published and thus it is difficult to make a direct comparison with other strainmeters.
We apply the tidal analysis code VAV version 5 (VAV05) 43 to the first-difference time series to generate an experimental tidal parameter set, which consists of 10 diurnal and 6 semi-diurnal tidal groups, for each interferometer. The set includes amplitude and phase (with respect to the local potential) of Q 1 , O 1 , P 1 , K 1 , 1 Ψ , 1 Φ , N 2 , M 2 , which are used as experimental tidal data set in the inversions. Tidal parameters are corrected for the effects of pre-whitening as follows: where A PW (φ PW ) and A (φ) are the amplitude (phase) for the pre-whitened and real series respectively, f is frequency and ∆t is sampling time 24 . inversions, probability distributions and correlations of inverted parameters. We use 32 complex tidal amplitudes (Q 1 , O 1 , P 1 , K 1 , Ψ 1 , Φ 1 , N 2 , M 2 at each interferometer), thus the total number of data is 64; each amplitude and phase pair is converted into its cosine (real) and sine (imaginary) terms, so that their dependence on the local-effect coefficients is linear. The model parameter space is 16- 2 , and the siting coefficients α, β, γ at each interferometer). We compute the frequentist probabilities of free model parameters by inverting 100,000 synthetic sets, each consisting of 32 complex tidal amplitudes which we generate by adding a Gaussian random noise (whose standard deviation is given by the formal error of the tidal analysis of real strain data) to each experimental tide (i. e., those given by the tidal analysis of first-differenced strain records). We use a Gaussian random noise since the tidal analysis gives results which are approximately normally-distributed around their "real" values, as a consequence of the almost-white noise spectrum of the first-difference strain sequences 41 , the least squares fitting which tidal analysis is based on, and the large number of analysed data. We assign the same uncertainty to all complex tidal amplitudes in each synthetic data set (which acts as an "experimental" data set for the inversion procedure) and minimize misfit using both the L 1 norm (sum of the magnitudes of the differences between "experimental" and theoretical tides) and L 2 norm (sum of squared magnitudes). As theoretical tides depend on FCN parameters non-linearly, we cannot assume that the misfit function has a single minimum with a quadratic-like shape even when the L 2 norm is considered.
Because of the large number of free parameters, we follow a two-step approach, where the parameter space is divided into two subspaces. The former subspace consists of the four FCN parameters, while the latter consists of the twelve siting coefficients. We rely on forward modelling and minimize misfit by using the downhill simplex algorithm 44 , which requires only function evaluations making almost no special assumptions, in the 4-dimensional FCN subspace. At each sampled point, misfit is computed after estimating the siting coefficients through χ-squared minimization of the distances between "experimental" and theoretical tides; as theoretical tides depend linearly on the siting coefficients, this estimate turns into a matrix inversion. Several tests prove that this approach gives quite similar results to and is much more efficient than minimizing misfit in the 16-dimensional parameter space at once.
Validation of the inversion procedure. We validate the inversion procedure by using it to invert complex amplitudes of Q 1 , O 1 , P 1 , K 1 , Ψ 1 , 1 Φ , N 2 and M 2 at each interferometer, obtained (i) computing the theoretical tidal strain components by using the reference 3 frequency-dependent Love and Shida numbers, T FCN , and Q FCN , (ii) combining them by using reasonable local effect coefficients 24,32 , and (iii) assigning them the uncertainties we got from the tidal analysis of the real extension records. 2 , and S Re( ) 2 , as well as siting coefficients for each interferometer, are free in the inversions. The same procedure has been repeated after halving 1 Ψ and 1 Φ uncertainties and/or fixing A Re( ) 2 to its reference value. Results are shown in Fig. 2b-d. Appropriate tests, whose results are not shown for brevity, demonstrate that fixing α, β, and γ does not affect PDFs and correlations appreciably. ocean loading. We compute ocean loading effects on strain tides by means of SPOTL code 45 . SPOTL does not include the most recent ocean models and uses integrated load Green functions computed from tables by 46 for three old Earth models, namely the Gutenberg-Bullen Earth model A 47 (GB) and two variants obtained after replacing top 1000 km by the continental shield crust and mantle structure of 48 and the oceanic crust and mantle structure of 48 . As for ocean models, we add FES2014b 34 and TPXO9-atlas 49 , released in 2019, to SPOTL, because recent models are more accurate than older ones, mainly in coastal areas. As for Earth models, we add integrated load Green functions computed from tables by 36 for a PREM 35 variant obtained after replacing top 40 km by continental crust (cPREM), and tables by 50 for IASP91 51 . We retain GB for comparison.
We compute ocean loading strain for all the diurnal tidal components included in both FES2014 and TPXO9-atlas, i. e. Q 1 , O 1 , P 1 , and K 1 , along the direction of each interferometer (ε xx ), perpendicularly to it ( yy ε ), and as shear strain (ε xy ). We do not use S 1 in our analysis, because of possible contaminations from environmental effects.
Tides 1 Ψ and Φ 1 are not included in any ocean model and we estimate their loading effects as follows. It is well known that the relation between the tide generating potential and the complex amplitude of the height of the www.nature.com/scientificreports www.nature.com/scientificreports/ ocean tides is strongly dependent on frequency as well as on the local and regional conditions that affect ocean dynamics. However, the ratio between the amplitudes of P 1 and K 1 ocean loading strain is approximately equal (within 10%) to the ratio of the related potential amplitudes for each interferometer and for each strain component. The only exception is ε yy when x is directed along the GAL16 interferometer, but yy ε gives a small contribution to GAL16 strain because  β α (Fig. 4b). Phases of P 1 and K 1 ocean loading strain differ by few degrees. Moreover, ocean loading contribution to measured deformation in the diurnal tidal band is small for both LNGS and LSC, affecting tidal strain amplitude by less than 5% and phase by at most few degrees (Fig. 1c). We assume that 1 Ψ and Φ 1 ocean loading strain is related to the tidal potential as for K 1 , but consider that the FCN also affects ocean tides. Thus, we estimate the effects of 1 Ψ and Φ 1 ocean loading on strain by multiplying K 1 ocean loading strain by two factors: (i) the ratio of the related potential amplitudes and (ii) the FCN resonant factor, computed inside the inversion procedure using Equations 32 and 33 of 4 .
Although semidiurnal tidal strain is not affected by the FCN, we also compute ocean loading for M 2 and N 2 , which are used for estimating strain cross-coupling due to siting effects and retrieve the large-scale Earth strain.
Strain load Green functions approximately decline as r −2 with distance r from a point load; thus, ocean loading effects depend on local crust/mantle structure also. Extent of the continental shelf around the two stations (Fig. 1a) suggests that cPREM could be more appropriate than average Earth models. Moreover, PREM is used to compute reference Love and Shida numbers in 3 . Since model parameter PDFs obtained from different ocean and Earth models are very similar to each other (Fig. 4c), for brevity we show most detailed inversion results from FES2014 and cPREM only.

Data availability
Datasets are available upon request.