Anelasticity and Lateral Heterogeneities in Earth's Upper Mantle: Impact on Surface Displacements, Self‐Attraction and Loading, and Ocean Tide Dynamics

Surface displacement and self‐attraction and loading (SAL) elevation induced by ocean tides are known to be affected by material properties of the solid Earth. Recent studies have shown that, in addition to elasticity, anelasticity considerably impacts surface displacements due to ocean tide loading (OTL). We employ consistent 3D seismic elastic and attenuation tomography models to construct 3D elastic and anelastic earth models, and derive corresponding averaged 1D elastic/anelastic models. We apply these models to systematically study the impact of anelasticity and lateral heterogeneity on M2 OTL displacements and SAL elevation. We find that neglecting lateral heterogeneities highly underestimates displacements and SAL elevation in mid‐ocean‐ridge regions and in some coastal areas of North and Central America. In comparison to PREM, 3D anelastic models can increase the predicted amplitudes of the vertical displacement and SAL elevation by up to 1.5 mm. The increased amplitudes reduce the discrepancy between GPS‐observed OTL displacements and their predictions based on PREM in places like Cornwall (England), Brittany (France), and the Ryukyu Islands (Japan). Applying our results to ocean tides, we discover that the impact on ocean tide dynamics exceeds the predicted SAL elevation correction with an RMS of about 1 mm, reaching an RMS of more than 5 mm in areas like North Atlantic or East Pacific. Due to the fact that such a value reaches the accuracy of modern data‐constrained tidal models, we regard the impact of anelastic shear relaxation as significant in tidal modeling.

systematically explored the sensitivity of load tides to perturbations in the density and elastic structures of the solid Earth. They found that surface displacements only weakly depend on perturbations in the density, but strongly on perturbations in the elastic moduli down to depth 500 km. Furthermore, they stated that perturbations in the density and elastic moduli of the crust mainly affect surface displacements in the vicinity of the load. Based on these findings, this study will focus on perturbations in elasticity of the upper mantle (down to depth 650 km) by considering lateral heterogeneities associated with 3D elastic Earth models. We will employ advanced 3D seismic elastic tomography models published recently by Karaoğlu and Romanowicz (2018) and Debayle et al. (2020) to derive the shear and bulk moduli of the mantle.
It is well known that attenuation, or the inverse seismic quality faction Q, causes the dispersion of seismic velocities and, correspondingly the reduction of elastic moduli from high frequencies to low frequencies.
As pointed out by Karaoğlu and Romanowicz (2018), early studies on body wave, surface wave, and normal model measurements (e.g., Anderson et al., 1965;Resovsky et al., 2005;Sailor & Dziewonski, 1978;Widmer et al., 1991) have identified a number of robust attenuation features of the Earth. Bulk attenuation is much weaker than shear-wave attenuation in the mantle. Shear-wave attenuation is low in the lithosphere and lower mantle, but high in a zone approximately coinciding with the low-velocity channel of upper mantle, and decreases sharply across the transition zone. These features lead us to consider anelasticity in a similar depth range to that of 3D elasticity, that is, the depth range of the upper mantle from 80 to 650 km.
The classical formulation for Q based on seismic and laboratory measurements is the so-called α-law (Equation 10 later in text), which says that Q has a dependence on frequency with power α (≠ 0) or is a constant (α = 0). Anelasticity associated with this α-law (Equation 7 later in text) can be constructed by a superposition of the standard linear solid models in a specific absorption band (Anderson & Minster, 1979;Zener, 1948). However, based on mineral physics and laboratory experiments of rocks under mantle temperature and pressure conditions, Sundberg and Cooper (2010) argued that Q and anelasticity could also be considered through the adoption of the Andrade model, while Lau and Faul (2019) as well as Ivins et al. (2020) favored the Extended Burgers Model (EBM). Lau and Holtzman (2019) argued that regardless of the viscoelastic model adopted, the frequency dependence of Q can be a significant effect for loading process acting on timescales longer than seismic timescales.
Applying the classical α-law over the seismic to tide frequency band, Benjamin et al. (2006), Kang et al. (2015), and Ray and Egbert (2012) showed that anelasticity has an impact on the geoid (body tide), and a value of α in the range of 0.15-0.30 can fit geodetic observations reasonably well. Different from the traditional normal mode (TNM) theory used to study body tides, Benjamin et al. (2006), Lau et al. (2015), and Wahr and Bergen (1986) developed an extended normal mode (ENM) approach, which predicts a more accurate perturbation of body tides due to anelasticity (Lau et al., 2016). Bos et al. (2015) and Wang et al. (2020) found that anelasticity based on the constant Q model can help bridging the gap between observations and predictions of surface displacements induced by the M2 tide loading. We will follow studies applying the classical α-law and adopt both constant and frequency-dependent Q models to consider anelasticity. How- The study of load tides began in 1960s (e.g., Longman, 1962Longman, , 1963, following the approach used in the study of Earth's free oscillation (Pekeris & Jarosch, 1958). Similar to the response to tidal forcing, the response of the Earth to surface loading is expressed by load Love numbers (LLNs). In a first approach, the LLNs as characteristics of a specific Earth model are used to derive the corresponding Green's function, allowing the convolution with OTL forcing in the spatial domain. In a second approach, the convolution of OTL forcing and response function is performed in the spectral domain of spherical harmonics, meaning a direct multiplication of the OTL with the LLN spectrum. In both approaches, a 1D Earth structure is assumed, as the Love numbers, by definition, are only dependent on the Legendre degree, meaning the Green's functions depend exclusively on the distance between the observation and the loading point. Recently, the spectral-finite element method (Martinec, 2000) has been modified by Tanaka et al. (2019) to calculate the responses of a 3D elastic Earth. Martinec (2000) formulated the surface loading problem of a viscoelastic gravitating sphere with an initial-value approach allowing for lateral variations in viscosity in the spectral domain. This formulation has been extended to include the effects of compressibility (Tanaka et al., 2011), and further to include the effects of lateral variability in elastic parameters (Tanaka et al., 2019). In this study, we extend the code to fully account for 3D variations and apply it to a global loading process like ocean tides.
In Section 2, we describe the numerical methods, the Earth models, and the ocean tide models in detail. In Section 3, we discuss the deviations in surface displacements and SAL elevation. In Section 4, we estimate the impact of the described effect on measurable quantities by comparing the GPS-data to modeled solid Earth motion (Section 4.1) and explore the impact on ocean tidal dynamics observable via altimetry (Section 4.2). The conclusions are presented in Section 5.

Methodology
We adopt the spectral-finite element method developed by Martinec (2000). The validity of this method for calculating the response of an elastic Earth with 3D structure has been demonstrated in Tanaka et al. (2019) who consider the effects of sphericity, self-gravitation, compressibility, and advection of pre-stress and lateral heterogeneity in elasticity. The authors applied the method in a sensitivity study considering rotational symmetry of the model and of the surface load. Here, the numerical code is extended for a fully 3D case and for anelasticity as a further perturbation of elastic parameters. The solution domain extends from the core-mantle boundary (CMB) to the surface of the Earth. The cut-off degree for spherical harmonics is 680. The corresponding resolution of the spatial longitude-latitude Gauss-Legendre grids is 2,048 E  1,022. The vertical number of layers (finite elements) from the CMB to the surface amounts to 146 with a layering thickness of 10 km in the upper mantle. The cut-off degrees for the seismic elastic and attenuation tomography models are up to 96 and 16, respectively, and the layering thickness used to invert for these models in the mantle is larger than 10 km. Hence, the settings of the spectral-finite element method are sufficient to resolve the features of the Earth's structure, as resolved by our chosen seismic elastic and attenuation models. In addition, as it will be shown in Figure 7a, the LLNs for the constructed Earth models mainly differ between degree 4 and 400; therefore, the chosen truncation degree 680 allows us to compute sufficiently accurate solutions presented in Section 4. Last but not least, to accelerate calculations, parallel programming with OpenMP has been implemented.

The Elastic and Anelastic Models
To consider the lateral heterogeneity and anelasticity in the upper mantle, we constructed four types of models in terms of elastic moduli: 1D elastic model (1D-e), 3D elastic model (3D-e), 3D anelastic model (3D-ae), and 1D anelastic model (1D-ae). We start with constructing 3D elastic/anelastic models and then derive the associated 1D models by taking the spatial average of 3D models on a sphere. HUANG ET AL. For the 3D elastic model, the shear and bulk moduli (μ and Κ, respectively) are derived from 3D S and P wave velocities (Vs and Vp, respectively). We adopt two 3D isotropic S wave velocity models, the companion Vs model of SEMUCB-UMQ (Karaoğlu & Romanowicz, 2018) and the DR2020s (Debayle et al., 2020), and derive perturbations in P wave velocity (δV p ) using the scaling factor from Montagner and Anderson (1989), that is, (1) which is consistent with French and Romanowicz (2014), Karaoğlu and Romanowicz (2018), and Panning and Romanowicz (2006).
For the 1D elastic model, the shear and bulk modulus, denoted by μ 0 and Κ 0 , respectively, at the seismic reference frequency w 0 = 2π rad/s, are derived from the shear and bulk modulus for the 3D elastic models (μ 3D-e and Κ 3D-e , respectively) by where 〈 〉 denotes spatial averaging within a layer. μ 3D-e and Κ 3D-e can be expressed further as where δµ 0 and δΚ 0 denote 3D perturbations with respect to the 1D elastic model.
For the 3D-ae models, we consider isotropic anelasticity, which can be accounted for by frequency-dependent perturbations of the elastic parameters (e.g., Dahlen & Tromp, 1998;Lau et al., 2016;Wahr & Bergen, 1986). As stated in the introduction, bulk attenuation is neglected. Then, the elastic parameters for the M2 tide frequency (w) are derived from those for the seismic frequency (w 0 ) as where δµ(w) is the anelastic perturbation of shear modulus with respect to the 3D-e model. For the 1D-ae models, the shear and bulk moduli are gained by spatial averaging where Equations 6 and 7 are applied. Equations 8 and 9 imply that the bulk modulus of the 1D-ae model is equal to that of the 1D-e model while the shear modulus differs by an amount equal to the average of δµ(w).
10.1029/2021JB022332 5 of 18 Taking the limit α → 0 in Equation 11, the constant-Q model reads as (e.g., Kanamori & Anderson, 1977) To be consistent with the seismic evidence that Q weakly depends on frequency over the seismic band, Benjamin et al. (2006) assumed Q is a constant across seismic frequencies but varies in the form of Equation 10 below the seismic frequency band. Correspondingly, they obtained Here, w m is the frequency above, which Q is a constant. Benjamin et al. (2006) found that when w m takes the value of 2 π × 3.09 × 10 −4 rad/s (the frequency of the longest-period free oscillation of the Earth) and α ranges between 0.2 and 0.3, body-tide and Chandler Wobble measurements can be fitted best.
We follow the parameterization of Benjamin et al. (2006) and take a value of α = 0.25 in Equation 13 to determine δµ. To study the impact of frequency dependence on anelasticity, we compare the value of δµ calculated by Equation 13 with that calculated by Equation 12. Bos et al. (2015) showed that the imaginary part of δµ produces a displacement Green's function, which is much smaller than that produced by the real part; thus we neglect the imaginary part of δµ. Also, in contrast to above mentioned studies where Q is assumed to be constant throughout the mantle or only depends on radial distance, we consider lateral variations of Q and adopt 3D Q models of the upper mantle.
In fact, there are very few global 3D attenuation models of the upper mantle published (Debayle et al., 2020). The spectral resolution of such Q models does not exceed 16 in terms of spherical harmonic degree, which is low compared to 680 used in our Earth models but is able to produce robust Q values with strong lateral variations (Adenis et al., 2017;Karaoğlu & Romanowicz, 2018). Among these models, we adopt two 3D Q models, SEMUCB-UMQ (Karaoğlu & Romanowicz,, 2018) developed at the University of California, Berkeley and QsADR17 (Adenis et al., 2017;Debayle et al., 2020) developed at the ENS de Lyon. Both Q models have a maximum degree of 16 and are jointly inverted with the corresponding 3D Vs models. The Berkeley model has been developed with a hybrid waveform inversion approach and with datasets consisting of fundamental modes as Love and Rayleigh waves in the period range from 60 to 400 s. The Lyon model has been built with an automated waveform inversion approach and with Rayleigh wave data in periods spanning from 50 to 250 s. Both Q models consider elastic focusing effects, source, and station terms. Although these two Q models are based on independent approaches and datasets, they are both well correlated with the Vs models in the upper part of the upper mantle, and both exhibit strong correlation with surface tectonic features down to depths from 200 to 250 km with low attenuation beneath continents and high attenuation beneath oceans. Figure 1a shows the depth variations of the shear modulus of the 1D-e and 1D-ae models, obtained from aforementioned seismic elastic and attenuation models. The top panel shows that, in general, μ of the 1D-e Berkeley model is larger than that of the Lyon model (except for depths greater than ∼500 km), but the difference between them is less than ∼3% with respect to PREM (Dziewonski & Anderson, 1981). The middle panel shows the reduction of shear modulus of the 1D-ae model with respect to the corresponding 1D-e model. We can observe that the two different α models (α = 0.0 and 0.25) predict similar reductions of the shear modulus, with differences less than 1%. In addition, the shear modulus reduction predicted by the Berkeley and by the Lyon models are similar (both staying at ∼ 4%) at depths from ∼280 to 650 km, but differ above 280 km, that is, in the high-attenuation zone, where the Berkeley model predicts a reduction of μ reaching up to ∼11% and the Lyon model reaching up to ∼5%. The bottom panel shows the variation of μ 6 of 18 of the 1D-ae models with respect to PREM. We can see that the shear modulus of the Lyon model is smaller than that of the Berkeley model at depths from ∼300 to ∼500 km and larger in other depths. In Figures 1b  and 1c, we exemplarily show lateral variations in the shear modulus of the 3D-ae models with respect to PREM at depth 135 km for α = 0.25. Clearly, we see that for both the Berkeley and Lyon models, the shear modulus is larger than the mean in the continental areas, which is associated with low attenuation, and smaller than the mean in the ocean areas, which is associated with high attenuation. Over the oceans, the Berkeley model is characterized by pronounced shear modulus variations along the mid-ocean ridge systems in the Indian, Pacific, and the Atlantic oceans while the Lyon model shows broad and strong variations beneath hot spots in Hawaii, the South Pacific, Indian, and Atlantic oceans. The similarities and differences of these two sets of models allow us to produce both comparable and differentiated results on the effects of anelasticity and heterogeneity. Last but not least, since α has a minor impact on the reduction of shear modulus, we only discuss the case of α = 0.25 in the following sections.

The Ocean Tide Model
The loading of ocean water mass anomalies induced by the M2 partial-tide is taken from the global ocean tide model TiME (Sulzbach et al., 2021;Weis et al., 2008) that has been run with an updated bathymetry scheme deduced by conservative interpolation from the RTOPO2-dataset (Schaffer et al., 2016). This bathymetry includes Antarctic sub-ice-shelf cavities that are important to realistically capture tidal resonance systems. The barotropic model employs temporal finite differences with a time step of approximately 180 s (1/240 M2-period) and has a high horizontal resolution of 1/12°. Energy dissipation is carried out by quadratic bottom friction, parametrized eddy-viscosity, and parametrized generation of internal waves.
The effect of SAL acceleration induced in response to the redistribution of water masses is included by a spectral approach employing spherical harmonical functions Y lm (Schindelegger et al., 2018). The ocean surface elevation Z = Σ lm Z lm Y lm undergoes a spectral decomposition, where individual spectral coefficients Z lm are multiplied with a combination of load Love numbers 1+k l -h l (Wang et al., 2012) before being reassembled by evaluating Here, ρ sw = 1,025 kg/m 3 and ρ E = 5,510 kg/m 3 are the mean density of seawater and the solid Earth, respectively. In contrast to the decomposition of the astronomical lunisolar tidal potential into spherical harmonics, this sum converges slowly and cannot be truncated at a low degree without introducing significant errors. Thus, the sum is evaluated up to a maximum degree/order of l max = 1,024, which is sufficient to precisely include the effects of SAL potential on ocean tide dynamics even in the most-shallow water areas (Schindelegger et al., 2018). Hereby, the barotropic acceleration g ∇ [SAL[static]), which is self-consistently included by repeated evaluation of Equation 14 at each time-step, directly affects the flow of water masses. Additionally, it is possible to apply a barotropic acceleration exerted by a static forcing term SAL [static], which is not influenced by the ad-hoc sea-surface elevation. Recently obtained M2-solutions have an open ocean (depth > 1,000 m) RMS of 3.39 cm when compared to the data-constrained FES14-atlas (Lyard et al., 2006(Lyard et al., , 2021, which corresponds to approximately 10%-15% of the expected M2-signal. After an initially transient model ocean state, the sea-surface modulations converge to a temporally harmonic oscillation. Analytically, the harmonic M2 partial-tide elevation Z can be written as (e.g., Schwiderski, 1980) where Z is the tidal elevation, ξ, δ, and σ are the amplitude, phase, and frequency of the tide, respectively, θ and φ are colatitude and longitude, respectively. For the purpose of OTL modeling, it is natural to calculate the responses of the Earth at two specific epochs t 1 and t 2 : 1.
To conserve ocean masses, a uniform layer of water mass has been subtracted from Z 1 (−0.38 mm) and Z 2 (−0.34 mm) of Sulzbach et al. (2021). The obtained results are shown in Figure 2. It is notable that large tides are mainly located along the coasts and in the central Pacific Ocean. The amplitude of the response field (R) at these two epochs can be expressed as

The Influence of Lateral Heterogeneity on Load Tide: Experimental Studies
Before discussing load tides due to realistic OTL, we perform three experiments to explore how the locations of lateral heterogeneity affect the elastic response and how large the impact will be. The locations of the water load and structural change are indicated by red and yellow circles in the left panels of Figure 3, respectively. Results of surface vertical displacement induced by the tidal loading are shown in the right panels of Figure 3. The water load disc has a positive height of 1.5 m, which is comparable to the tidal amplitude.
The structural change has values taken from the 1D anelastic Berkeley model. From the figure we see that, when the lateral heterogeneity is situated right below or in the vicinity of the loading (Experiments 1 and 2), there is a decrease of vertical displacement below the loading and an increase in the vicinity of the load. The change is in general up to 0.7 mm and around 1% of the vertical displacement of the 1D elastic model. When the lateral heterogeneity is far from the loading, for example, 60° in Experiment 3, we see almost no change of vertical displacement below the loading and its surrounding but some tiny change (less than 0.005 mm) in the region with lateral heterogeneity.

The Impact of Lateral Heterogeneity and Anelasticity on Load Tide and SAL Elevation
We now focus on the perturbations in the load tide, that is, the vertical displacement (U), and SAL elevation, caused by the perturbations in elastic moduli. We calculate the differences in amplitude of U and SAL elevation between the 1D-ae Berkeley/Lyon model and PREM and the difference between the 3D-ae and 1D-ae Berkeley/Lyon model. We found that the results of SAL elevation are rather similar to those of U, an indication that the geoid undulation in SAL elevation is small compared to the vertical displacement. Therefore, we only show U results in Figure 4 while SAL elevation results are referenced in Figure S1.
Above all, we concentrate on the deviations of 1D-ae models shown in Figures 4a (Berkeley model) and 4c (Lyon model). We can observe that the amplitude of U mostly exhibits an increase over oceans, whereas a decrease over continents adjacent to tide loads, which is consistent with findings from Experiments 1 and 2. Moreover, Figures 4a and 4c show that the large differences in amplitude of U mainly occur in coastal regions, for example, the western coastal area of France, Colombia, and New Zealand, with the magnitude of above 0.5 mm. In open ocean areas, the differences are in general below 0.4 mm. The fact that significant deviations occur in shelf regions, that is, the near field of large tides is due to large local tidal amplitudes ( Figure 2) and a weak shallow structure of the upper mantle caused by high attenuation (Figure 1a, bottom panel). Furthermore, the Lyon model produces larger deviations in open ocean areas than the Berkeley model, which is attributed to a larger shear modulus reduction at larger depths from ∼300 to ∼500 km (bottom panel of Figure 1a) where perturbations in elastic moduli have a big impact on far-field (open ocean area) solutions. Setting of experiments (left panels) to explore the impact of the relative location of lateral heterogeneity with respect to ocean water load on load tides, and profile of vertical displacements with distance from the load (right panels). The load is indicated as a red circle in the map and the structure change is as a yellow circle. The profile of vertical displacement is marked as a white dashed line. Figures 4b and 4d show the differences in amplitude of U between 3D-ae and 1D-ae models. It is clear that if the lateral heterogeneity in elasticity and anelasticity is not considered, both the Berkeley and Lyon 1D-ae model will underestimate the amplitude of U in mid-ocean ridge areas or beneath hot spots of the Indian, Pacific, and Atlantic oceans by 0.1-0.2 mm. This amount is quite significant considering that the amplitude of U for the 1D-ae model in these places is less than 0.4 mm. In addition, evident underestimations also occur in the coastal areas of New Zealand, North and Central America with a magnitude of ∼0.2-0.3 mm. Moreover, the 3D-ae Berkeley model predicts a notable underestimation (0.3-0.4 mm) by the 1D-ae model in the coastal area of Iceland. On the other hand, significant overestimations of around 0.4-0.6 mm by both 1D-ae models occur in the northwestern coastal area of Australia, the shelf-region of North Europe, the Labrador Sea, and the mouth of Hudson Bay. In contrast to the Lyon model, the 3D-ae Berkeley model also predicts notable overestimations of 0.3-0.4 mm by the 1D-ae model in the Antarctic shelf region. Table 1 lists the peak (maximum and minimum) amplitude differences of U and SAL elevation between the 1D-ae models and PREM, between the 3D-ae models and 1D-ae models, and between the 3D-ae models and PREM. There, the Berkeley and Lyon models both predict that the peak amplitude differences in U and SAL  elevation between the 3D-ae and 1D-ae model are of the order of ∼1 mm and comparable to those between 1D-ae and PREM. In addition, the maximum amplitude differences of U and SAL elevation between the 3D-ae model and PREM are ∼1.5 and ∼1.1 mm for the Berkeley and Lyon model, respectively. However, it should be noted that the Berkeley model predicts the maximum difference occurring between Iceland and Greenland, while the Lyon model points to the Bay of Biscay.

The Root Mean Square of Difference in U and SAL Elevation Between Different Model Pairs
In the previous section, we have shown that significant deviations in the amplitude of U and SAL elevation can occur due to the consideration of anelasticity and lateral heterogeneity. To quantify the impact on amplitude, as well as phase of U and SAL elevation simultaneously, we calculate the RMS of differences in U and SAL elevation for the same model pairs as shown in Table 1. Again, the results of SAL elevation are similar to those of U in the ocean areas. This time we place the results of U in Figure S2 but show the results of SAL elevation in Figure 5 because the latter can be used as a reference in Section 4.2, where we discuss the feedback to ocean tide dynamics.
We first look at the RMS between the 1D anelastic model and PREM shown in Figures 5a and 5c. Features similar to the amplitude differences shown in Figures 4a and 4c can be observed: (1) large RMS values above 0.3 mm occur in the coastal regions where the tidal amplitude is large, and the SAL elevation is sensitive to the weak shallow structure of the upper mantle; (2) the RMS in the open ocean areas is relatively small (less than 0.25 mm) due to small local tidal amplitudes and a relatively strong deeper structure of the upper mantle; (3) the Lyon model predicts larger RMS values in the open ocean areas due to a weaker deeper structure of the upper mantle (from ∼300 to 500 km). For the results between the 3D-ae and 1D-ae models, that is, Figures 5b and 5d, we see that the RMS in the mid-ocean ridge systems or below the hot  Table 2 shows the maximum RMS values. We see that the Berkeley and Lyon models both predict that the maximum RMS values of U and SAL elevation between the 3D-ae and 1D-ae models are comparable to those for the 1D-ae model and PREM, both at a range of ∼0.7-0.9 mm. The Berkeley model predicts the maximum RMS for U and SAL elevation between the 3D-ae model and PREM around 1 mm, which is slightly higher than ∼0.7-0.8 mm predicted by the Lyon model. Bos et al. (2015) demonstrated that GPS-observed vertical displacements due to ocean tide loading show a maximum discrepancy of 2-3 mm in comparison with predicted values based on PREM at Newlyn in Cornwall, Southwest of England and Brest in Brittany, Northwest France. Wang et al. (2020) also noticed a maximum discrepancy of over 1.5 mm on the Ryukyu Islands even with the best regional ocean tide model NAO99Jb. In Figure 6, we show the difference in amplitude of vertical displacement between the 3D-ae model and PREM for both regions. It can be seen that consideration of anelasticity and lateral heterogeneity increases the amplitude of vertical displacement by ∼0.7-0.8 mm in Newlyn and Brest, 0.12-0.36 mm on Ryukyu Islands, predicted by both the Berkeley and Lyon model. The sum of uncertainties in predicted vertical displacements caused by ocean tide models and in GPS observations themselves is ∼0.2-0.4 mm in Newlyn and Brest (Bos et al., 2015) and 0.0-∼0.8 mm on Ryukyu Islands (Wang et al., 2020). Obviously, the effects of anelasticity and lateral heterogeneity in the upper mantle based on global seismic models are not able to completely reduce the maximum discrepancy in these regions.

Feedback to Ocean Tide Dynamics
The ocean tide model TiME introduced above is now used to study the feedback of laterally variable elasticity and anelasticity of the Earth's upper mantle structure on the tidal dynamics in the world's ocean.
To estimate the tidal ocean-surface height change induced by SAL potential, individual sets of load Love numbers for the PREM and the 1D Lyon/Berkley anelastic Earth models have been calculated. The differences in LLNs (see Figure 7a) are most prominent between degrees l = 4-400 and express different Earth responses to ocean tide loading. Using different sets as input for the SAL-module, see Equation 14, we repeat the initial simulation that provided the loading pattern introduced in Section 2 and compare the resulting tidal elevations by calculating their respective RMS. As the SAL-forcing for PREM and 1D-ae models is self-consistently included to full extent by utilizing the respective sets of LLNs, the corresponding ocean response to anelasticity can be directly computed. On the other hand, the effects of lateral heterogeneity cannot be encoded in LLNs. Therefore, the ocean response patterns arising for the 3D-ae earth models are  Table 2 The Even though the SAL-fields of the PREM and Lyon/Berkeley 1D-ae solutions differ only on a scale of 1 mm RMS (see Table 2), the impact on tidal dynamics is much stronger pronounced and reaches RMS values of between 5 and 10 mm in some areas (e.g., North Atlantic, East Pacific), also showing relevant magnitudes in some open ocean regions (see Figures 7b and 7d).
While the spatial patterns of the RMS-deviations are similar for both 1D-ae models, the induced deviations by the Lyon 1D-ae model are larger. This can be traced back to larger deviations of the respective LLNs to the PREM-LLNs for low degrees (see Figure 7a), since those cause the most prominent feedback on ocean tide dynamics.

of 18
A consideration of residual 3D-ae effects (see Figures 7c and 7e) causes additional minor deviations with an RMS that reaches amplitudes of up to 2 mm in some areas but is below 1 mm in most regions. Interestingly, the spatial deviation-patterns differ from the ones discussed above, showing augmented response in the Pacific region while patterns in the Indian and Atlantic Ocean are strongly diminished. Intercomparing between the two model families, the response patterns strongly resemble each other while 3D effects in the Berkeley model are stronger pronounced. This can only be partially attributed to apparently larger 3D-SAL deviations for the Berkeley model (compare Figure 5) as the ocean does not respond locally but collectively to periodic forcing patterns. For example, the resulting ocean response to two individual point-excitations might interfere constructively or destructively depending on several factors including the relative phase of the two disturbances. Due to this, it is not possible to naively relate the expected ocean response to the 15 of 18 apparent magnitude of a certain SAL-deviation. To understand the ocean response mechanisms to a certain differential SAL-forcing one would need to intermediate this effect by decomposing the excitation pattern into oceanic normal modes (e.g., Müller, 2008).
Concerning the observed response patterns for 1D-ae and 3D-ae solutions the intermediation of the ocean tide response by only a few near-resonant normal modes explains both the similarities between the patterns and the amplification of their magnitude (∼5 mm) with respect to the original deviations in the SAL elevation (∼0.5 mm).

Conclusion
Applying the jointly inverted 3D seismic tomography model and attenuation models from two research groups, we constructed four types of Earth models in terms of shear and bulk moduli: the 1D elastic model, 1D anelastic model, 3D elastic model, and 3D anelastic model. In calculating the dispersion of shear modulus from the seismic to M2 tide frequency, we considered two values of α, that is, 0.00 and 0.25 that are typically used in literature (e.g., Benjamin et al., 2006;Bos et al., 2015). It turns out that α reduces the shear modulus by less than 1% when increasing from 0 to 0.25. We adopt 0.25 for after-mentioned calculations, which could be seen as an estimation of the maximum effect of anelasticity. We found that, compared to PREM, the 1D anelastic models experience a high reduction of shear modulus in the high-attenuation zone. For the Berkeley model, this reduction is as large as 12%. Due to strong lateral variations in the seismic tomography model and attenuation model, the 3D anelastic models compared to the 1D anelastic models experience a significant decrease of shear modulus (up to ∼30% for the Lyon model, 19% for the Berkeley model) in the mid-ocean ridges or beneath hot spots where the attenuation is high. In addition, we found that the 1D anelastic models in general underestimate the shear modulus on the continents but overestimates it over the oceans.
For numerical modeling, we applied the spectral-finite element method capable of handling lateral heterogeneities in viscosity of a spherical gravitating and compressible Earth model (Tanaka et al., 2011). This approach has been modified recently for lateral heterogeneities in elastic moduli (Tanaka et al., 2019) and further extended here for 3D variations and anelasticity.
To estimate the impact of lateral heterogeneity, we performed a number of synthetic tests for which we showed a deviation in the response at the level of 1% at maximum mainly associated with regions where loading and structural changes coincide. In the next step, we considered the 1D and 3D anelastic models for the ocean tide loading of the M2 astronomical tide. Applying the M2 solution of the global ocean tide model TiME, we assessed the impact of anelasticity and lateral heterogeneity on the ocean tide loading.
We found that, compared to PREM, the 1D anelastic models experience a significant increase in the amplitude of U and SAL elevation in coastal regions due to a large local tide amplitude and the weak shallow structure of the upper mantle. If lateral heterogeneity is neglected, a considerable underestimation of U and SAL elevation occurs in places like the mid-ocean ridges or beneath hot spots while a comparable overestimation occurs in coastal areas of Australia, northern Europe, eastern Canada, etc. These features are confirmed by both the Berkeley model and Lyon model. In addition, the 3D anelastic Berkeley model predicts a maximum increase of ∼1.5 mm between Iceland and Greenland compared to PREM while the Lyon model predicts ∼1.1 mm in the Bay of Biscay.
The 1D anelastic Lyon model predicts higher RMS values in the open ocean areas due to a weaker structure at depths from ∼300 to 500 km than the Berkeley model. Moreover, the 1D anelastic models predict more widespread high RMS values than the 3D anelastic models, an indication that 1D anelasticity effect is more significant at global scale than the effect of lateral heterogeneity. However, both Lyon and Berkeley models confirm that the maximum RMS values of U and SAL elevation between the 3D and 1D anelastic model are comparable to those between the 1D anelastic model and PREM, both at the range of ∼0.7-0.9 mm. In addition, the Berkeley model predicts the maximum RMS for U and SAL elevation between 3D anelastic model and PREM around 1 mm, which is slightly higher than ∼0.8 mm for U and 0.7 for SAL elevation of the Lyon model.

of 18
Our results based on global seismic models can be applied to reduce the discrepancy of GPS-observed vertical displacements with their predictions based on PREM. For example, the 3D anelastic Berkeley and Lyon models predict that anelasticity and lateral heterogeneity in the upper mantle increase the amplitude of vertical displacement by ∼0.7-0.8 mm in Newlyn and Brest, 0.12-0.36 mm on Ryukyu Islands. The obtained amplitude increase here is not able to reduce the maximum localized discrepancy in the corresponding region to the range of the sum of uncertainties in GPS observations and in vertical displacement predictions caused by errors in ocean tide models. This implies that regional elastic and attenuation models with higher resolutions might be worthy of further considerations. On the other hand, we neglect perturbations in density and elastic moduli of the crust and perturbations in density of the mantle associated with 3D structure, which has an impact on displacements as well.
The interpretation of the tidal signal showed some consequences in ocean tide dynamics when modifying the SAL elevation according to the deviating solid-earth response. Inserting the SAL elevation back into the tidal equation, we calculated the RMS in the tidal elevation between PREM and 1D anelastic models and the RMS between 3D and 1D anelastic models. The obtained results have a number of implications: As a secondary tide-generating force, the SAL potential induced by anelasticity is found to modify the ocean tide dynamics. While the RMS in tidal elevation between PREM and 1D anelastic models is below 5 mm in most regions, pointing to a small perturbation compared to present M2 amplitudes, it exhibits the same order of magnitude as the accuracy of modern data-constrained tidal models (e.g., Stammer et al., 2014). Especially when considering applications like satellite gravimetry that requires even more precise tidal solutions (Flechtner et al., 2016), this perturbation can be regarded as relevant for data-unconstrained tidal modeling. The same holds for satellite-data constrained tidal modeling where only the sea-level E = U + Z is observable by satellite measurements. Here, a correction of E by the expected Load-tide height U is necessary to obtain ocean tide solutions Z. Therefore, the deviations of around 1 mm in U between elastic and anelastic models directly translate into biased M2-elevations. This can be regarded as important because modern data-constrained tidal models are able to reach an accuracy of only 3 mm in the open ocean. Based on these findings, the mean relaxation of the shear modulus should be considered for tidal modeling, whereas the effect of additional 3D variabilities in mantle elasticity appears to be of second order for global studies.

Data Availability Statement
The modelling data discussed in this paper are published separately at GFZ as Huang, P., Sulzbach, R., Tanaka, Y., Klemann, V., Dobslaw, H., Martinec, Z., & Thomas, M. (2021) "Anelasticity and lateral heterogeneities in Earth's upper mantle: impact on surface displacements, self-attraction and loading and ocean tide dynamics," https://doi.org/10.5880/GFZ.1.3.2021.003, under the Creative Commons Attribution 4.0 International Licence (CC BY 4.0). The data publication is containing the amplitude and root-mean-square fields of surface vertical displacement and self-attraction and loading due to ocean tide loading − the M2 tide derived from model TiME (Sulzbach et al., 2021), and the root-mean-square fields of M2 tide. The fields have been calculated for the 1D elastic solid Earth model PREM and 3D and 1D anelastic models. Figures 4-7, and S1, and Tables 1 and 2 can be easily reproduced from these data fields applying the calculus discussed in the paper. The anelastic Earth models can be determined using the calculus discussed in Section 2.2 by making use of the elastic and attenuation tomography models from the University of California, Berkeley (Karaoğlu, H. & Romanowicz, 2018) and the École Normale Supérieure (ENS) de Lyon (Debayle et al., 2020), respectively. All response fields (U and SAL) are calculated with the spectral-finite element method (Martinec, 2000;Tanaka et al., 2019).