The sensitivity of ocean tide loading displacements to the structure of the upper mantle and crust of Taiwan Island

Ocean tide loading (OTL) displacements are sensitive to the shallow structure of the solid Earth; hence, the high-resolution spatial pattern of OTL displacement can provide knowledge to constrain the shallow Earth structure, especially in coastal areas. In this study, we investigate the sensitivity of the modeled M2 OTL displacement over Taiwan Island to perturbations of three physical quantities, namely, the density, bulk modulus, and shear modulus in the upper mantle and crust. Then, we compare the sensitivity of the modeled M2 OTL displacement to Earth models with the sensitivity to ocean tide models using root mean square (RMS) differences. We compute the displacement Green’s function and OTL displacement relative to the center of mass of the solid Earth (CE) reference frame, analyze the sensitivity to the three physical quantities in the CRUST1.0 model and the Preliminary Reference Earth Model (PREM), and present their spatial patterns. We find that displacement Green’s functions and OTL displacements are more sensitive to the two elastic moduli than the density in the upper mantle and crust. Moreover, their distinctive sensitivity patterns suggest that the three physical quantities might be constrained independently. The specific relationships between the perturbed structural depths and the distance ranges of peak sensitivities from the observation points to the coastline revealed by the shear modulus can mitigate the nonuniqueness problem in inversion. In particular, the horizontal tidal components observed by the Global Positioning System (GPS) can yield better results in inversions than the vertical component owing to the smaller OTL model errors and the higher structural sensitivity (except for the shear modulus in the asthenosphere).


Introduction
Ocean tide loading (OTL), which includes the influences of displacement, gravity, strain and tilting, is the response of the solid Earth to the periodic redistribution of ocean water due to the gravitational effects from the Moon and Sun. The material properties of the upper mantle and crust inside the Earth control the OTL response at Earth's surface that can be observed by geodetic techniques, such as very long baseline interferometry (VLBI), the Global Positioning System (GPS), gravimeters, strainmeters, and tiltmeters. Measurements of the OTL response can therefore be expected to provide independent constraints on the material properties in the upper mantle and crust of the solid Earth.
Over the past several decades, such studies relying on gravity, strain, and tilt tidal data to infer Earth's structure have been explicitly attempted, yet there are some deficiencies in these methods, such as a sparse station coverage, an insufficient quality or accuracy of the data, or local structural sensitivity to strain and tilt (Zschau 1976;Baker 1980Baker , 1984Baker et al. 1991). For instance, OTL displacements were observed by VLBI with millimeter precision (Thomas et al. 2007), where the deficiency was a sparse station coverage. Owing to the development of You and Yuan Earth, Planets and Space (2021) 73:193 GPS techniques, GPS estimations of OTL displacements can reach sub-millimeter precision with a dense network (Yuan et al. 2009(Yuan et al. , 2010(Yuan et al. , 2013Yuan and Chao 2012;Bos et al. 2015;Penna et al. 2015;Martens et al. 2016b). Since the spatial distribution of OTL displacements can vary finely, especially in coastal areas, due to the complexity of regional ocean tides, GPS-observed tidal displacements with their present-day accuracy and spatial resolution show great advantages in inferring the structure of the upper mantle and crust inside the Earth (Ito and Simons 2011;Yuan and Chao 2012;Yuan et al. 2013;Bos et al. 2015;Martens et al. 2016b). Ito and Simons (2011) first constrained the density, shear modulus, and bulk modulus in the lithosphere and asthenosphere beneath the western United States using GPS-derived tidal displacement data. However, they did not consider the geocenter motion correction and used the residual OTL displacements with the assumption of an errorless model for the body tides (or solid Earth tides) (Yuan and Chao 2012). Bos et al. (2015) found that considering anelastic dispersion at tidal frequencies in transversely isotropic Earth models, with the shear modulus reduced by 8.5% in the asthenosphere, decreased the vertical discrepancies of OTL displacements over Cornwall and Brittany to the same order of magnitude as the sum of other errors, including the errors in ocean tide models and numerical computations, and GPS-related errors. More recently, a similar study focused on the region around the East China Sea (Wang et al. 2020).
Although GPS-estimated OTL displacements might be useful for investigating the elastic moduli and density structures in the shallow Earth beneath coastal regions, this approach needs to be explored further. To study the possibility of inverting GPS observations of OTL displacements for the two elastic moduli and density, the sensitivity of OTL displacements to perturbations of the Earth's structure in the upper mantle and crust needs to be analyzed exhaustively (Ito and Simons 2011;Martens et al. 2016a).
Recent studies compared displacement Green's functions and OTL displacements based on spherically symmetric, non-rotating, elastic, and isotropic (SNREI) Earth models at regional or global scales (Penna et al. 2008;Wang et al. 2012;Yuan and Chao 2012;Yuan et al. 2013;Bos et al. 2015;Dill et al. 2015;Martens et al. 2016b). Yuan and Chao (2012) and Yuan et al. (2013) attempted to discriminate shallow structures from three Earth models using GPS-observed OTL displacements, which seems to be difficult. Martens et al. (2016a) also compared displacement Green's functions and OTL displacements based on various SNREI Earth models and analyzed the OTL displacement sensitivities for the M 2 tidal harmonic, but the analysis was performed only at 64. 7°N and 341°E in Iceland. Martens et al. (2016b) and Martens and Simons (2020) computed the vector discrepancies and root mean square (RMS) differences of OTL displacements between different Earth models and ocean tide models and then compared and examined the OTL displacement sensitivities, revealing that OTL displacements are more sensitive to ocean tide models than to Earth models. Wang et al. (2012) and Dill et al. (2015) considered the factors of the structural depth and angular distance from the station to the mass load in computing displacement Green's functions and mass loading displacements and especially showed that discrepancies in the horizontal components of mass loading displacements are much larger than those in the vertical component. However, these studies can provide only a general knowledge of sensitivity to the Earth's structure and cannot reflect the sensitivity of a single physical quantity, including the shear modulus, bulk modulus, or density. Ito and Simons (2011) and Martens et al. (2016a) formulated and analyzed the displacement Green's function sensitivities to perturbations of the density, bulk modulus, and shear modulus as functions of the structural depth and angular distance. Nevertheless, Ito and Simons (2011) did not consider the effect of layer thickness on the displacement Green's functions, and the density sensitivities were contaminated by extraneous fluctuations to the two elastic moduli from density perturbations (Martens et al. 2016a). Both studies found that the sensitivities to these three physical quantities are distinctive, and the density sensitivity is weak in thin layers of the Preliminary Reference Earth Model (PREM). More importantly, Martens et al. (2016a) suggested that the nonuniqueness of inversion problems might be mitigated with a strategically distributed and dense GPS network. Furthermore, the possibility of inversion at present depends on the sensitivity to the structure of interest and the accuracy of the OTL model, which is related to the regional ocean tide model and computation method.
We aim to explore the possibility of probing the density, bulk modulus, and shear modulus in the upper mantle and crust beneath the Taiwan region using OTL displacement observations. This analysis focuses on the general sensitivity to Earth models and the individual sensitivity to the density, bulk modulus, and shear modulus. Since constraining the Earth's structure using OTL displacement observations requires that the sensitivity level of OTL displacements to the Earth's structure exceeds the OTL model errors, the sensitivity level of OTL displacements to ocean tide models also needs to be considered.
In this paper, we first investigate the sensitivity of displacement Green's functions as a function of angular distance to linear perturbations of the elastic moduli and density in the different regions of Earth models. Then, we examine the sensitivity of the modeled M 2 OTL displacements to linear perturbations of the three physical quantities in different layers of Earth models and finally compare the sensitivity of the modeled M 2 OTL displacements to Earth models with their sensitivity to ocean tide models.
We choose Taiwan Island to discuss the sensitivity of OTL displacements to linear perturbations of the density and two elastic moduli for the following reasons: (1) Taiwan Island is surrounded by the sea, where the OTL response is strong and the spatial spectrum of OTL displacements for the M 2 tidal harmonic can be fully shown.
(2) Our study shows that the vector discrepancies and RMS differences of the modeled M 2 OTL displacements for the horizontal components using different Earth models are generally larger than those using different ocean tide models. Taiwan Island may be an appropriate area to probe the bulk modulus, shear modulus, and density using GPS-observed OTL displacement data.

Methods
Load Love numbers describe the elastic response of the solid Earth to surface mass loads (Longman 1962), such as oceans, glaciers, lakes, and the atmosphere. The load Love numbers can be derived by solving the six following first-order linear differential equations (Longman 1962(Longman , 1963Farrell 1972): dy 5 dr = −4π Gρy 1 − (n + 1) r y 5 + y 6 dy 6 dr = − 4πGρ(n + 1) r y 1 + 4π Gρn 1 r y 2 + (n − 1) r y 6 where λ and μ are the Lamé parameters, r is the geocentric distance, ρ is the density, g is the gravity acceleration, G represents the universal gravitational constant, and y 1 , y 2 , y 3 , y 4 , y 5 , and y 6 represent the radial and tangential displacements, the radial and tangential stresses, the perturbed potential, and a variation as defined by the fourth formula in Eq.
(2), respectively. We apply the MATLAB package ELLN (Chen et al. 2015(Chen et al. , 2018 to solve Eq. (1) from the Earth's center to the surface with the corresponding boundary conditions, such as Eq. (3), for the surface boundary conditions: where R is the Earth's radius and g R is the gravity acceleration at the surface.
Since OTL displacements are modeled and measured at the surface, the radial dependence of the load Love number is not considered here; we compute the load Love numbers only at the surface, which are expressed as h ′ n , l ′ n and k ′ n . The load Love numbers can be defined by the following relations (Farrell 1972): The load Green's functions describe the response of the solid Earth to a unit point mass load and are derived by the infinite sums of the load Love numbers and Legendre functions (Longman 1963;Farrell 1972): (2) where G U and G V denote the vertical and horizontal displacement Green's functions, respectively; θ denotes the angular distance from an observation station to the unit point mass load; m e is the Earth's mass; and P n denotes the nth degree Legendre polynomial. Usually, asymptotic expressions and Kummer's transformation are applied to speed-up the convergence of load Green's functions (Farrell 1972;Guo et al. 2004). If the degree n is large enough, the load Love numbers h ′ n , nl ′ n and nk ′ n become constant (Farrell 1972), and the asymptotic values of the load Love numbers are given as The formulas of the displacement Green's functions can be expressed as follows: We compute the vertical and horizontal displacement Green's functions for angular distances θ from 0.001° to 180°: an interval of 0.0002° is used from 0.001° to 0.02°, an interval of 0.001° is used from 0.02° to 0.05°, an interval of 0.01° is used from 0.05° to 1°, an interval of 0.1° is used from 1° to 10°, an interval of 0.5° is used from 10° to 90°, and an interval of 1° is used from 90° to 180°. In this way, the effects of ocean mass loads in the near field can be included during the modeling of OTL displacements.
According to Farrell (1972), the OTL displacements can be modeled as follows: where U denotes the OTL displacement; λ and φ are the longitude and co-latitude, respectively, of the observation point; ρ w is the mean seawater density; H denotes the height of the ocean tide point; G denotes the displacement Green's function; A denotes the azimuth of the ocean tide point relative to the observation point; and S G denotes the integral area of the ocean.
The OTL displacement values are predicted by the convolution of the displacement Green's functions with an ocean tide model. We evaluate Eq. (10) using the SPOTL package (Agnew 1997(Agnew , 2012, which can conveniently combine a global ocean tide model with regional ocean tide models.

Displacement Green's function sensitivities
We introduce the sensitivity kernel K (Ito and Simons 2011;Martens et al. 2016a) to quantify the sensitivities of displacement Green's functions to linear perturbations of the density, bulk modulus, and shear modulus: where p denotes the perturbed parameter, i.e., the density, shear modulus or bulk modulus; l denotes the perturbed layer; ∆m l p denotes the linear variation of the perturbed parameter p in layer l of the original Earth model m; and ∆G denotes the variations between the perturbed Green's functions G(m + ∆m l p ) and original Green's functions G(m). According to Tarantola and Albert (2005) and Martens et al. (2016a), we denote the three physical quantities in common log space, so a + 1% linear perturbation ∆m l p , i.e., log10(1.01), is equal to 0.0043.
When we compute the sensitivity kernels of displacement Green's functions to one certain parameter (one of the two elastic moduli or the density), we keep the other two parameters fixed. We perturb the density, bulk modulus, and shear modulus separately in an original Earth model to derive three perturbed Earth models and then compute the displacement Green's functions with the perturbed Earth models and original Earth model. The CRUST1.0 (Laske et al. 2013) model is combined with PREM (Dziewonski and Anderson 1981) as the original Earth model in our study. Since the sediment layers in the CRUST1.0 model are thin and their contributions to the OTL displacements are negligible, we average the sediment layers with the upper crust to become the outermost layer in each model (Guo et al. 2004). Taiwan Island can be divided into six regions ( Fig. 1) based on the CRUST1.0 model, which correspond to the six sets of structures shown in Fig. 2.
As shown in Figs. 3, 4, the sensitivity patterns can be distinguished between the two elastic moduli. The sensitivity patterns of the shear modulus exhibit transitions between positive values and negative values, and for the bulk modulus, they are mostly positive. We also find that the displacement Green's functions are more sensitive to perturbations of the elastic moduli in the near-surface layers, even though these layers are much thinner than the deeper layers in the upper mantle. As the depth of the perturbed layers increases, the sensitivities are still significant in the low-velocity zone (LVZ) and gradually decrease from the transition zone to the lower mantle at depths ranging from 670 to 771 km. However, the sensitivities are restored in the lower mantle at depths from 771 to 2741 km since the thickness of the lower mantle layer is massive, which means a high number of perturbations. Moreover, the deeper the layers in which we perturb the elastic moduli, the greater the angular distances the peak sensitivities show. In particular, the perturbation depths of the shear modulus are approximately equal to the angular distances of the sensitivity transitions converted into kilometers, as shown in Fig. 3. The sensitivities of the horizontal displacement Green's functions to density perturbations are always negative in both the near and the far fields, and the sensitivities of the vertical displacement Green's functions exhibit positive sensitivities in the far field. The density sensitivity can therefore be discerned from the sensitivities of the elastic moduli in both the near and the far fields, and the density sensitivity can be discerned in the vertical and horizontal components. In addition, the density sensitivity kernels are weak in the lithosphere above 80 km depth and in the lower mantle at depths from 2741 to 2891 km. However, the sensitivity kernels become large in the lower mantle at depths from 771 to 2741 km. Obviously, the density kernels are sensitive to the layer thickness, not the perturbation depth.
In this section, we compute the displacement Green's functions relative to the center of mass of the solid Earth (CE) reference frame, as well as the center of mass of the entire Earth system (CM), which includes the solid Earth and surface loads, reference frame (Blewitt 2003). As the sensitivities of the displacement Green's functions in the CE frame to linear perturbations of the two elastic moduli are consistent with the corresponding sensitivities in the CM frame, their sensitivities are independent of the reference frame. However, the displacement Green's functions are more sensitive to density perturbations at certain far fields in the CM frame, as shown in Additional file 1: Figure S1, than to those in the CE frame, especially to perturbations in much thicker layers.
For completeness, we investigate the sensitivities of the displacement Green's functions to perturbations of the two elastic moduli and density in the other models (TW2, TW3, TW4, and TW5) in the CE frame. The sensitivity characteristics of the displacement Green's functions for one certain parameter (one of the two elastic moduli or the density) are consistent among these models, although there are some differences in their sensitivity patterns. Not surprisingly, these differences between TW1 and the other models are mainly found in the thicknesses of the perturbed layers, and the displacement Green's functions are generally more sensitive to thicker layers with the same name in different models. As we mentioned above, perturbations of the elastic moduli at deeper depths result in the peak sensitivity corresponding to a larger angular distance. In addition, the influences caused by the magnitudes of the three physical quantities on the sensitivity are negligible.
Since the thickness of the perturbed layer can affect the sensitivity kernel, we further perturb the layers with a fixed thickness of 20 km in TW1 at depths from 0 to 700 km. Clearly, the peak sensitivities are much smaller when perturbations of the elastic moduli occur in deeper layers, and the corresponding angular distances are larger, which indicate that perturbations of the elastic moduli at certain depths give rise to variations in the displacement Green's functions at specific angular ranges. We note in Fig. 5 that the high sensitivity of the displacement Green's functions can arise from perturbations of the elastic moduli only within tens of kilometers near the surface. The sensitivity kernels corresponding to perturbations of the bulk modulus are great above 40 km depth, and for the shear modulus, they are approximately Although the sensitivity kernels for the shear modulus exhibit transitions at certain angular distances, beneath 60 km depth, the displacement Green's functions are generally more sensitive to perturbations of the shear modulus, especially for the vertical component. The density kernels are weak in thin perturbed layers and basically remain consistent in the near field.

OTL displacement sensitivities
Analyzing the OTL displacement sensitivity to the upper mantle and crustal structure can provide necessary knowledge of the independent constraints on the three physical quantities (density, shear modulus, and bulk modulus), even the structural depth. OTL displacements at the surface depend not only on the angular distances from observation stations to mass loads, as in the case of the displacement Green's functions assuming a point mass load, but also on specific positions since ocean tides exhibit complex distributions at both spatial and temporal scales. We further investigate the sensitivity characteristics of the OTL displacements on Taiwan Island (with a resolution of 0.25° × 0.25°) to linear perturbations of the three physical quantities, especially the corresponding relations between the depth of the structural perturbations and the angular distance from the station to the coastline.
We begin with the sensitivity characteristics of OTL displacements to linear perturbations of the bulk modulus, shear modulus, and density in different layers of the original Earth model. Since the sensitivity characteristics of OTL displacements to perturbations of one certain quantity in these models (TW1, TW2, TW3, TW4, and TW5) are consistent and the strongest responses of the OTL displacements shown in Figure 6 appear in the TW1 region, we only show the results of TW1. For the other models, the sensitivity patterns show slight graphic shifts of the peak sensitivities and small amplitude differences (the most sensitive model is TW4, as shown in Additional file 1: Figure S3). Because the vector discrepancies of the modeled M 2 OTL displacements to +1% linear perturbations in the lithosphere are smaller than 0.03 mm, which can be negligible with the present accuracy, we only show the sensitivity patterns for +10% linear perturbations of the three physical quantities. First, we compute the variations in the modeled M 2 OTL displacements to +10% linear perturbations of the density, bulk modulus, and  Planets and Space (2021) 73:193 shear modulus separately in each crust and mantle layer of TW1. As OTL displacements are more sensitive to the upper mantle and crustal structure (Farrell 1972;Ito and Simons 2011;Yuan et al. 2013), we only show the results for the perturbed layers above the transition zone (400-600 km depth) in Figure 7. Then, we calculate the sensitivities of the modeled M 2 OTL displacements based on Eq. (12) to perturbations of the three physical quantities. The sensitivity patterns shown in Additional file 1: Figure  S2 corresponding to +10% linear perturbations of the density, bulk modulus, and shear modulus are consistent with the patterns of the vector discrepancies shown in Figure 7: where K l U,p is the sensitivity of OTL displacements to parameter p perturbed linearly in layer l; U denotes the OTL displacement modeled as Eq. (10); r denotes the location of the station relative to the mass load; H denotes the ocean tide model; G denotes the displacement Green's functions; ∆U denotes the vector difference of OTL displacements considering both amplitudes and phases; and ∆m l p denotes the variation of the perturbed parameter p in perturbed layer l.
As shown in Fig. 7, the modeled M 2 OTL displacements are also sensitive to linear perturbations of the elastic moduli but not to those of the density. In particular, the horizontal OTL displacements are more sensitive to perturbations of the elastic moduli than the vertical displacements, except to those of the shear modulus in the upper mantle. We note that the horizontal displacement Green's functions vary at wider angular ranges than the vertical displacement Green's functions shown in Fig. 3, which may explain the higher sensitivity of the horizontal OTL displacements. The regions in which the horizontal Fig. 3 Differences between the original and perturbed displacement Green's functions. The differences of displacement Green's functions for the horizontal (left) and vertical (right) components to + 1% linear perturbations of the bulk modulus (red), shear modulus (green), and density (blue) in the upper crust at depths of 0-13.88 km (top row), the middle crust at depths of 13.88-24.52 km (second row), the lower crust at depths of 24.52-32.15 km (third row), and the region above the LVZ (LID) at depths of 32.15-80 km (bottom row) separately from TW1 as a function of angular distance. The displacement Green's functions are scaled by 10 12 Rθ in the CE frame, where R is the radius of the Earth and θ is the angular distance You and Yuan Earth, Planets and Space (2021) 73:193 components are sensitive to perturbations of the elastic moduli in the upper mantle can extend to most of Taiwan Island and even exhibit a uniform pattern. Since the LVZ is 80-220 km deep, considering the specific relation between the angular range of the peak sensitivity and the perturbation depth, which is near or exceeds the maximum distance across Taiwan Island, when we perturb deeper layers beneath the LVZ, their sensitivities still exhibit a uniform pattern. It is worth noting that since the sensitivities of the displacement Green's functions for the shear modulus exhibit transitions at certain angular distances, as shown in Figs. 3, 4, as we perturb the shear modulus from the upper crust to the region above the LVZ (LID), the regions of peak sensitivities gradually move inland from western Taiwan, which reveals that the shear modulus at certain depths mainly controls the OTL responses on the surface within the specific distance ranges from stations to the coastline. The most sensitive regions of the modeled M 2 OTL displacements are mainly located in northern and western Taiwan Island. In fact, Fig. 6 shows that these regions are exactly where the responses of OTL displacements for the M 2 tidal harmonic are strong. The variations of the OTL displacements in the horizontal components in response to + 10% perturbations of the bulk modulus can reach the sub-millimeter level except in the lower crust and for the north component in the region above the transition zone, consistent with the magnitudes of the estimation accuracy of GPS tidal displacements (Yuan and Chao 2012;Yuan et al. 2013), and for the vertical component, the variations can reach the sub-millimeter level in the upper mantle. The variations of the OTL displacements in the east component in response to perturbations of the shear modulus are also at the sub-millimeter level except in the lower crust and the region above the transition zone, and for the vertical component, the variations increase dramatically in the LVZ and the region above the transition zone. Because the thicknesses of the LVZ and the region above the transition zone are thicker than the crustal layers, the responses of OTL displacements to perturbations of the elastic moduli are stronger. Here, we reiterate that  Yuan Earth, Planets and Space (2021) 73:193 the variations of the OTL displacements controlled by perturbations of one certain physical quantity in a specific layer cannot be added numerically to the variations of the OTL displacements caused by perturbations of the other physical quantities in any layer.
Since the layer thickness can affect the sensitivity, we normalize the sensitivity of OTL displacements by the corresponding layer thickness to attenuate this effect. Moreover, by implementing this operation, we can more clearly show the sensitivity kernels to linear perturbations of a single physical quantity in different layers. Figure 8 shows that the unit sensitivities of the modeled M 2 OTL displacements at 24.5°N on Taiwan Island to perturbations of the elastic moduli from the upper crust to the LVZ are separated for the horizontal components, and the sensitivities gradually decrease with increasing perturbation depth, especially for the bulk modulus. The unit sensitivities to perturbations of the bulk modulus in the crustal layers dramatically change in the vertical component. It is worth noting that the peaks and troughs of the unit sensitivities to the shear modulus perturbations, which correspond to specific angular distances, significantly translate in the vertical component. The density sensitivities are weak. We also analyze the sensitivities of the modeled M 2 OTL displacements to linear perturbations of the three physical quantities in layers with a fixed thickness of 20 km, and the effects of the perturbed depth on the sensitivities are very clear above 80 km for perturbations of the elastic moduli, especially for the shear modulus within the specific ranges of angular distances.

Sensitivity comparisons between Earth models and ocean tide models
We further discuss the general sensitivity of OTL displacements to different Earth models. The Earth models used in our study include TW1, AK135 (Kennett et al. 1995), AK135-F (Kennett et al. 1995;Montagner and Kennett 1996), STW105 (Kustowski et al. 2008), Huang 1D (Huang et al. 2014, and 1066A (Gilbert et al. 1975). Figure 9 depicts the structural deviations of the different Earth models relative to TW1 in log space above 800 km depth. We note that the structural deviations are mainly above 220 km depth. Due to the magnitudes of the three physical quantities and different layer thicknesses, 1066A and Huang 1D are significant, especially for Huang 1D in the upper mantle. The deviations in the two elastic moduli are significantly larger than those in the density. In the crustal layers, the deviations of the elastic moduli can exceed 0.2 in log space, which corresponds to a 58% perturbation, and deviations form Huang 1D in the upper mantle reach 0.1 in log space, which corresponds to a 26% perturbation. Most deviations are less than 0.03 in log space, which corresponds to a 7% perturbation.
Then, we compute the displacement Green's functions in the CE frame using these Earth models, and their respective differences are still relative to TW1. The sediment layers in TW1 and Huang 1D, as well as the ocean layers in STW105, AK135-F, and PREM, are correspondingly averaged with their outermost crustal layers. Figure 10 shows that the differences in the displacement Green's functions between 1066A and TW1 in the very near field (angular distances smaller than 0.1°) can be attributed to structural deviations above approximately 13 km depth, and the differences among AK135-F, PREM, STW105, 1066A, and TW1 at angular distances from 0.1° to 0.4° probably result from structural deviations above 40 km depth. Deviations in the structure at depths from 80 to 350 km cause differences in the displacement Green's functions between Huang 1D and TW1 at angular distances from 0.4° to 3°.
We model the M 2 OTL displacements based on the six Earth models and calculate the vector and RMS differences of the OTL displacements among these different Earth models. For the comparison, we also calculate the vector and RMS differences of the OTL displacements among five ocean tide models, namely, DTU10 (Cheng and Andersen 2011), EOT20 (Hart-Davis et al. 2021), HAMTIDE11A (Zahel 1995;Taguchi et al. 2014), FES2014b (Lyard et al. 2006;Carrère et al. 2016), TPXO9-atlas (Egbert and Erofeeva 2002;Egbert et al. 2010), and FESJB, which combines FES2014b with the NAO99Jb regional model (Matsumoto et al. 2000), and we attempt to compare the sensitivity level of OTL displacements to the Earth models with the sensitivity level of OTL displacements to these ocean tide models. The RMS differences of the modeled M 2 OTL displacements between different ocean tide models and FESJB are shown in Fig. 11, as are the RMS differences between different Earth models and TW1. The RMS differences of the modeled M 2 OTL displacements with the results between different ocean tide models and FES2014b are shown in Fig. 12, as are the results between different Earth models and PREM. Figures 11 and 12 show that the RMS differences of the modeled M 2 OTL displacements for the horizontal components, especially the east component, among the different Earth models are generally larger than those among the ocean tide models, which indicates that the OTL displacements for the M 2 tidal harmonic in the horizontal components are generally more sensitive to the Earth's structure than to the errors in those ocean tide models. However, the sensitivity levels of OTL displacements to different Earth models in the vertical component are relatively low, except for the Huang 1D and 1066A models. We find that the sensitivities of displacement Green's functions to linear perturbations of the elastic moduli are independent of the reference frame. However, for linear perturbations of the density, the sensitivities at some far fields are significantly larger in the CM frame than those in the CE frame, especially in the lower mantle at depths from 771 to 2741 km. The sensitivities of displacement Green's functions and OTL displacements to density perturbations seem to be mainly governed by the thickness of the perturbed layer, not the depth, which is different from the sensitivities to the elastic moduli perturbations.
The mass variations of the solid Earth caused by density perturbations in the structure enhance the gravitational potential, and the increased gravitational force attracts mass loads, which explains the mostly negative sensitivity to density perturbations and the higher sensitivity from a much thicker layer (Martens et al. 2016a).
For the sensitivities of displacement Green's functions, the unique characteristics (weak for the density, positive for the bulk modulus, and transitions between positive and negative for the shear modulus) and the relationships between the perturbation depth and the distance range of the peak sensitivity are consistent with the prerequisites of inversion. However, OTL displacements are the superposed effects of the responses from spherical shells at different depths, so the combined positive and negative anomalies of the three physical quantities at different spherical shells can result in equivalent variations of OTL displacements. Hence, GPS-observed OTL displacement data cannot correspond to a unique model of the solid Earth, especially with the complicated spatial distribution of ocean tides, and inversions suffer from the nonuniqueness problem (Martens et al. 2016a). Fortunately, shear modulus perturbations at certain depths can induce trough and peak sensitivities of OTL displacements at specific angular distances from observation stations to the coastline, as shown in Fig. 7, which can help to mitigate the nonuniqueness problem in inversions. The sensitivity kernels of OTL displacements to linear perturbations of the elastic moduli can directly connect the perturbation depths with corresponding OTL displacements at specific GPS stations to refresh the a priori model during the inversion process. However, this phenomenon seems to be significant only for the shear modulus above approximately 80 km depth in all components and for the vertical displacement extending to the LVZ.
As mentioned previously, the inversions for the three physical quantities at distinct depths based on GPSderived OTL displacement data with the present accuracy still need to be explored, mainly to alleviate the contamination from errors in the OTL model (arising from ocean tide models, Earth models, coastline refinement, and the convolution scheme) (Bos and Baker 2005;Penna et al. 2008) and errors in the body tide model (Yuan and Chao 2012;Yuan et al. 2013).
The accuracy and resolution of ocean tide models are the main limitations to accurately modeling OTL displacements in the long term, especially in certain coastal areas (Penna et al. 2008;Yuan and Chao 2012;Yuan et al. 2013;Martens et al. 2016b). In our study, we evaluate the sensitivity level of the modeled M 2 OTL displacements to different ocean tide models and compare it with the sensitivity level of the OTL displacements to different SNREI Earth models. As a result, the modeled M 2 OTL displacements in the horizontal components, especially the east component, are more sensitive to the Earth models than to errors in the ocean tide models; however, for the vertical component, the regional ocean tide model around Taiwan Island still needs to be more accurate. Although this is beyond the scope of our study, for the horizontal components, the sensitivity level of the modeled M 2 OTL displacements to different Earth models generally exceeds the GPS tidal estimation accuracy reported by Yuan and Chao (2012) and Yuan et al. (2013), and the sensitivity level to different ocean tide models is close to that accuracy. Meanwhile, for the vertical component, the sensitivity level to ocean tide models is also close to the GPS tidal estimation accuracy, and the sensitivity level to different Earth models is lower than that accuracy.  Consequently, the GPS-observed OTL displacements for some tidal constituents with accuracies of approximately 0.1 mm and 0.3 mm in the horizontal and vertical components, respectively (Yuan and Chao 2012;Yuan et al. 2013), might be more suitable to constrain the two elastic moduli and density beneath the Taiwan region compared with constraints on the regional ocean tide models. We note that the modeled M 2 OTL displacements in the horizontal components are more sensitive to the elastic moduli perturbations in the lithosphere than those in the vertical component, even for the bulk modulus in the LVZ. Moreover, the sensitivity level of the OTL displacements for the horizontal components to different Earth models exceeds the sensitivity level to ocean tide models. Hence, for the horizontal components, GPS-observed OTL displacement data can yield better results than those for the vertical component in inversions with the present accuracy. The sensitivities of OTL displacements in the vertical component to the shear modulus become strong in the LVZ, which can help to better constrain the shear modulus in the LVZ Wang et al. 2020). In addition, the accuracy of the regional ocean tide model needs to be improved to predict OTL displacements more accurately, especially in the vertical component.
For the body tide residuals in GPS-estimated OTL displacements, if they exist, they should be removed from the coastal tidal residuals to obtain structural information in the shallow Earth (Yuan and Chao 2012). Since body tides mainly reflect structural information in the Earth's deep interior (Yuan and Chao 2012;Yuan et al. 2013) and are long-wavelength signals (Martens et al. 2016a,b), these would be common values for the spatial scale of Taiwan Island. Moreover, the sensitivity patterns of OTL displacements for the horizontal components to perturbations of the elastic moduli in the upper mantle are almost uniform over Taiwan Island, as shown in Fig. 7. Some of these structural signals and the body tide residuals may be mixed in the inland tidal residuals. Thus, removing inland tidal misfits between GPS-observed and modeled OTL displacements on Taiwan Island from coastal tidal residuals may reduce the structural signals of the upper mantle. This work should be treated with caution with regard to retaining the structural signals of the LVZ, and even the LID, at coastal stations as much as possible, and the constrained depth and resolution of the radial model in inversions should be further investigated.

Conclusions
We investigated the sensitivities of displacement Green's functions and modeled M 2 OTL displacements to linear perturbations of the density, bulk modulus, and shear modulus separately in the crust and upper mantle beneath the Taiwan region. The results show that displacement Green's functions and OTL displacements are more sensitive to the elastic moduli, and the sensitivity patterns of the three physical quantities are distinguishable, which suggest that the three physical quantities might be constrained independently. Moreover, the sensitivity characteristics of the shear modulus can mitigate the nonuniqueness problem in inversions. The GPSobserved OTL displacement data with the present accuracy, especially in the horizontal components, within a strategically distributed GPS network might be used to constrain the density, shear modulus, and bulk modulus in the lithosphere and asthenosphere beneath the Taiwan region. We recommend a more accurate regional ocean tide model around the region to reduce the related errors in the OTL model. In addition, the misfits between GPSobserved and modeled OTL displacements at inland stations of Taiwan Island need be analyzed carefully to distinguish the structural signals in the asthenosphere from possible errors in the body tide model.