Ocean tide loading displacements in western Europe: 2. GPS-observed anelastic dispersion in the asthenosphere

GPS-observed vertical ocean tide loading displacements show in Cornwall, southwest England, and in Brittany, northwest France, discrepancies of 2–3mm with predicted values based on the isotropic Preliminary Reference EarthModel for themain tidal harmonicM2, yet in central Europe the agreement is better than 0.5mm. By comparison of ocean tide models and validation with tide gauge observations, we demonstrate that the uncertainties in the former are too small to cause this disagreement. Furthermore, we find that different local models of the crust and different global elastic reference models derived from seismological observations can only reduce the observed discrepancies to 1–2mm, which still exceeds the GPS observational uncertainty of 0.2–0.4mm. It is customary to use the elastic properties of the Earth as given by seismic models. Previously, there has been insufficient evidence to determine how tomodify these properties during the transformation from seismic to tidal frequencies to account for possible anelastic dispersion in the asthenosphere, and so this effect has been ignored. If we include this effect, then our discrepancies reduce further to 0.2–0.4mm. This value is of the same order as the sum of the remaining errors due to uncertainties in the ocean tide models and in the GPS observations themselves. This research provides evidence in western Europe of a reduction of around 8–10% of the seismic shearmodulus in the asthenosphere at tidal frequencies. In addition, we find that the asthenosphere absorption band frequencies can be represented by a constant quality factor Q.


Introduction
The solid Earth is deformed due to the gravitational attraction of the Moon and Sun and due to the varying weight of the ocean tides [Baker, 1984;Agnew, 2007]. These two effects are called the solid Earth tide (or Earth body tide) and ocean tide loading (OTL), respectively, and can be observed by all precise space geodetic techniques such as very long baseline interferometry [e.g., Petrov and Ma, 2003] and GPS [e.g., Khan and Tscherning, 2001;Allinson et al., 2004]. For most geophysical research, for example, measurement of tectonic velocities, glacial isostatic adjustment, and vertical land motion at tide gauges, these tidal displacement signals are removed. Such corrections applied in space geodetic analysis need to be of the highest accuracy since Penna et al. [2007] and Fu et al. [2012] have shown that inaccurate values can create spurious fortnightly, semiannual, and sometimes annual effects. In this research we focus on the OTL displacement, which shows a greater degree of spatial variability largely on account of the more complicated distribution of surface loading compared with the Moon and Sun's direct gravitational forces.
The elastic properties of the Earth used in the OTL computations are normally based on global seismic models. This is not necessarily correct since at the lower tidal frequencies the elastic properties could be slightly different due to dissipation effects, mainly in the asthenosphere. This effect is called anelastic dispersion. Due to dissipation, there is also a delay between the applied stress and the resulting strain. Zschau [1978] estimated that this could cause a few degrees of phase delay in vertical OTL displacement values. Nevertheless, anelasticity has usually been ignored because this phenomenon has always been smaller than the uncertainty in the OTL values caused by errors in the ocean tide models. However, since the launch of the TOPEX/Poseidon altimeter satellite in 1992, and its successors Jason 1 and 2, the accuracy of global ocean tide models has improved substantially [Stammer et al., 2014]. At the same time the number of continuous GPS stations at which OTL can be observed has increased to a few thousands through national reference networks and international bodies such as the International GNSS Service (IGS) and the European Reference Frame (EUREF). Ito [Dziewonski and Anderson, 1981] for computing OTL values for the western United States. Ito and Simons [2011] used the discrepancies between GPSobserved and predicted OTL displacements to study the elastic properties of the lithosphere and asthenosphere. In this research we follow a similar approach for western Europe.
Parts of western Europe (Cornwall, southwest England, and Brittany, northwest France) experience exceptionally large vertical OTL displacements, up to 45 mm amplitude for the main tidal harmonic M 2 , based on both GPS observations [e.g., Vergnolle et al., 2008] and modeled values [e.g., Penna et al., 2008]. These arise mainly due to nearby large-amplitude ocean tides in the Celtic Sea and the Atlantic Ocean which, as a very rough first approximation, can be represented by a disc of seawater with a radius of 375 km and a thickness of 1.2 m. When we represent this area by a half-space with a compressional velocity V P = 8 km/s, shear velocity V S = 4 km/s, and density ρ = 3 g/cm 3 , then we compute near the edge of the disc (representing the coast of Cornwall and Brittany), a vertical displacement of around 40 mm, close to what is observed. Furthermore, if we compute the corresponding vertical and horizontal stress underneath the center of the disc load as a function of depth [Sigmudsson et al., 2010], then we find that at a depth of 400 km the horizontal stress is almost zero while the vertical stress already has reduced by half. This highlights the fact that our GPS-based OTL observations presented in this paper are mostly influenced by the properties of the lithosphere and asthenosphere.
In this paper we first present discrepancies between GPS observations of vertical M 2 OTL displacement and predictions using Green's functions computed using the commonly adopted isotropic version of PREM and thereafter investigate their cause. We assess in detail the quality of the ocean tide models and their subsequent contribution to the OTL error budget. We then evaluate the use of different Earth models with different elastic and anelastic properties to determine optimal OTL predictions for western Europe.

Observed OTL Displacements
To obtain GPS observations of OTL displacement, data from 259 GPS stations across western Europe (see Figure 1) were analyzed using the (GNSS-Inferred Positioning System) GIPSY v6.1.2 software in kinematic precise point positioning mode. The data set comprised all Natural Environment Research Council British Isles continuous GNSS Facility (NERC BIGF) (http://www.bigf.ac.uk), EUREF (http://www.epncb.oma.be), IGS (http://www.igs.org), and coastal Institut Géographique National (IGN) stations (http://rgp.ign.fr) that had at least 1100 days of data available during the window 2007.0 to 2013.0. The only exception was for the 11 stations Journal of Geophysical Research: Solid Earth 10.1002/2015JB011884 subject to large OTL effects but which had slightly fewer data days (at least 2.1 years of data and typically greater than 2.5 years, shown by Penna et al. [2015] to result in only very small degradations compared with using 4 years of data). Seven additional stations met the 1100 days criterion, but results were not available due to GPS data processing problems. The stations used, their data center, approximate longitude and latitude, and the number of days of GPS data are provided in Table S1 in the supporting information, and the processing strategy follows that of Penna et al. [2015]. Jet Propulsion Laboratory (JPL) "repro1" fiducial orbits, clocks, and Earth orientation parameters were held fixed, and data were processed in 30 h sessions centered on noon on the UTC day to minimize daybreak effects. Ambiguities were fixed to integers, and station coordinates were estimated every 5 min, extracting the values from the central 24 h of each session. Also estimated every 5 min were receiver clock offsets, zenith wet delays (mapped from slant to zenith using the VMF1 gridded mapping function), and north-south and east-west tropospheric delay gradients. A priori zenith hydrostatic and wet delays were applied using European Centre for Medium-Range Weather Forecasts-based values. The solid Earth tides were modeled according to International Earth Rotation Service (IERS) 2010 conventions [Petit and Luzum, 2010], and first-order OTL corrections were applied using the Gutenberg-Bullen Earth model, the FES2004 ocean tide model [Lyard et al., 2006], and the hardisp program available in Petit and Luzum [2010]. Thus, our GPS analyses provided residual deviations from these predicted OTL values, which were then added back to obtain the total GPS-observed OTL displacement, per tidal harmonic and coordinate component.
In order to estimate the tidal parameters, the 5 min coordinate estimates were averaged to one value every 30 min and then time series formed by concatenation. Four to six years of data were deemed to provide a large sample from which the principal tidal harmonics could be estimated. The analyses, with optimal coordinate and zenith wet delay process noise values of 3.2 × 10 À6 km/√s and 1.0 × 10 À7 km/√s respectively (with a carrier phase observational standard deviation of 10 mm), are described in more detail in our companion paper [Penna et al., 2015], and we demonstrate that the maximum errors of the obtained tidal parameters are about 0.2-0.4 mm, by performing rigorous testing with synthetic tidal signals.
The residual vertical position time series were analyzed using the ETERNA software package [Wenzel, 1996] to estimate the amplitude and phase lag at the M 2 harmonic. The standard deviations of the estimated tidal amplitudes lie between 0.1 and 0.2 mm for the vertical component of harmonic M 2 . Comparison with the modified Lomb-Scargle periodogram [Scargle, 1982] shows amplitude agreements of around 0.1 mm.
The GPS results are given in the IGS08 realization of ITRF2008: a reference frame that has its origin at the center of mass of the deformed solid Earth plus the applied load. Following Blewitt [2003], we will call this the CM frame. However, OTL predictions presented in this paper have been computed in the frame that has its origin at the center of mass of the deformed solid Earth only, called the CE frame. The transformation between the two frames depends on the global distribution of the ocean tides. Using the ocean tide models described in section 3, we estimate that the transformation of the GPS observations from the CM to the CE frame has been undertaken with an error of around 0.15 mm for harmonic M 2 . The total uncertainty is the sum of the uncertainty in the CM-CE transformation and the uncertainty of the tidal analysis output from ETERNA, using simple propagation of errors.
The GPS-observed vertical M 2 OTL phasors, in the CE frame, are listed in Table S2 and shown in Figure 1. The very large OTL values over Cornwall, southwest England, and Brittany, France, are clearly visible as well as the low values over central Europe. A general anticlockwise rotation of the OTL phasors can be seen along the coast from south to north. This is caused by the tidal Kelvin wave for M 2 in the Atlantic Ocean that propagates from south to north along the coast of Europe. As a result, there is an increasing phase lag from south to north, although the OTL generated by the ocean tides in the North Sea disrupts this pattern. Also note the large OTL values over Scotland, southern Ireland, and along the west coast of Iberia.

Predicted OTL Displacements
In the previous section we discussed the observed OTL, but we have also computed the OTL displacements using the recent global ocean tide model FES2012 [Carrère et al., 2012] on a spherical Earth with a commonly used Green's function, based on the isotropic version of PREM [Francis and Mazzega, 1990;Guo et al., 2004] which we will call PREM-iso. More details of the computations will be presented in section 4. The differences between our GPS-observed OTL values and the predicted OTL values are shown in Figure 2 for the vertical component for Journal of Geophysical Research: Solid Earth 10.1002/2015JB011884 the M 2 harmonic, in which the blue circles represent the 95% observational confidence interval using the aforementioned propagation of errors. The radius of the error circle varies between 0.4 and 0.6 mm, which is smaller than the 1-2 mm listed by Vergnolle et al. [2008], who also studied OTL in Brittany using a shorter GPS data set.
In Figure 2 one can observe that the residuals in central Europe are less than the observational error circles (0.4 mm), which shows that the uncertainties in the OTL and the solid Earth tide models in that area are very small. Furthermore, it can be seen from Figure 2 that in Cornwall and in Brittany the discrepancies can reach 3 mm, corresponding to 7% of the total vertical M 2 OTL signal, which we will investigate in the following sections. It is probably due to the large uncertainty in their OTL estimates that Vergnolle et al. [2008] failed to observe the discrepancies between observations and model predictions that are shown here.

Ocean Tides
In this section we assess the accuracy of the ocean tide models around western Europe, which in section 4 will help to assess the accuracy of the predicted OTL values.
We have used six recent global tide models: DTU10 [Cheng and Andersen, 2010], EOT11a [Savcenko et al., 2011], FES2012 [Carrère et al., 2012], GOT4.7 [Ray, 1999], HAMTIDE [Zahel, 1991[Zahel, , 1995, and TPXO8 [Egbert and Erofeeva, 2002]. FES2012, HAMTIDE, and TPXO8 are barotropic tide models into which tide gauge and satellite altimetry observations have been assimilated. DTU10 and EOT11a took the FES2004 tide model [Lyard et al., 2006] as their first estimate and then used satellite observations to improve it. For GOT4.7 a similar approach was adopted, but instead of FES2004, a combination of several global, regional, and local tide models was used, blended across their mutual boundaries [Ray, 2013]. The six global tide models were used to compute a mean ocean tide model. The standard deviations of the M 2 vector differences of each of these models from the mean model are shown in Figure 3. For this plot, all tide models were interpolated onto a grid with a spatial resolution of 1/16°. It can be seen that in the open ocean all tide models differ by a few centimeters at worst. However, near the coast their differences increase to a few tens of centimeters with the notable exception of the Bristol Channel where the differences reach the meter level. However, the resulting error in predicted vertical M 2 OTL displacement at distances of more than 120 km from this area will be less than 0.3 mm. Another problematic area is the Wash estuary, along the east coast of England. It is mainly the GOT4.7 model that differs from the other models in this region.
For 76 shelf water pelagic tide stations, Stammer et al. [2014] show that the standard deviation of their difference with the tide models is for M 2 around 3-6 cm on the European shelf. We have performed a similar comparison using the coastal gauges and bottom pressure recorders shown in Figure 3. The harmonic constants observed at these locations were taken from the tidal bank formerly maintained by the International Hydrographic Organization and from the values given by Davies et al. [1997]. Pelagic tidal constants were taken from the International Association for Physical Sciences of the Oceans (IAPSO) [Smithson, 1992] and Global Undersea Pressure (GLOUP) data banks. The comparison was performed for each of the areas shown in Figure 3 and was also extended to tide gauges along the west coast of Iberia [Quaresma and Pichon, 2013]. The values for the harmonic constants used from the different data sources are listed in Table S3. Dividing the region into separate areas helps to indicate where the accuracy of the ocean tide models is high and where some problems still exist. As was mentioned in section 2, OTL is large in Cornwall and Brittany, with vertical M 2 observed versus model discrepancies of around 2-3 mm. For that reason, we have also defined a small area around Cornwall that will illustrate the OTL contribution from the very near ocean tides.
The results of the comparison are shown in Table 1, which lists the mean value, and the standard deviation, of the length of the vector difference between the tide model and the observations. There are only two tide gauges in the Cornwall sea area, and while insufficient to obtain reliable statistics, the mean difference is nevertheless only a few centimeters for most tidal models. The value is larger for the GOT4.7 model because of its coarser grid size. For the interpolation of this model to the position of the gauge, grid cells north of Cornwall were used which are not equal to the tides south of Cornwall. . The standard deviations of the vector differences between six global ocean tide models (DTU10, EOT11a, FES2012, GOT4.7, HAMTIDE, and TPXO8) and their mean for harmonic M 2 . The white labeled polygons define the areas for which the OTL contributions have been determined, and the white dots represent the locations of the tide gauges and bottom pressure recorders. The Wash is also labeled as the ocean tide models exhibit large differences there.

10.1002/2015JB011884
For the Celtic Sea area, we divided the observations into pelagic and tide gauge observations. The pelagic observations compare very well (subcentimeter level) with all tide models. The tide gauge observations are all located near the entrance of the Irish Sea and show differences of a few centimeters. Thus, the accuracy of the tide models is quite good over most of this area but deteriorates near the northern boundary.
Agreements of several centimeters between tide gauge observations and the tide models are also obtained in the English Channel, North Sea, and Irish Sea. The DTU10 and EOT11a models perform rather poorly in the Irish Sea, which can be explained by the older FES2004 model (on which they are both based, with the addition of limited satellite altimetry data in this region) having several large incorrect tidal values near the eastern boundary of the Irish Sea.
In the Bristol Channel the global tide models differ not only substantially among themselves but also from the tide gauge observations. The relatively good performance of GOT4.7 is illusionary because it only has tidal information near the entrance, where the tides are smaller. In this area the fine grid resolution of FES2012 and TPXO8 helps to obtain a much better fit with the observations than DTU10, EOT11a, and HAMTIDE.
The comparison between tide gauges and models is rather poor along the west coast of France. This is caused by gauges near the Pertuis Breton basin in the northeastern part of the Bay of Biscay where the tides change rapidly [Nicolle and Karpytchev, 2007]. Finally, all models perform well in the waters west of Scotland and around the Atlantic coast of Iberia.
In this section we have demonstrated that the ocean tide models exhibit no systematic biases and that they have the same order of magnitude of differences among themselves (see Figure 3), as with the tide gauge observations (see Table 1). By considering actual ocean tide models, the spatial correlation of the tides is correctly taken into account, which as Scherneck [1993] noted is necessary to obtain realistic propagation of errors in the estimation of the corresponding OTL uncertainty. Therefore, we may select any of these recent models as a representative candidate and the intermodel standard deviation as representative of the error.
In the literature a variation of 1% for the mean seawater density value can be found [Bos and Baker, 2005]. We use a mean value of 1030 kg/m 3 and estimate that the error in the vertical M 2 OTL displacements caused by its uncertainty is less than ± 0.2 mm.
Another source of error in the OTL values for stations near the coast is the accuracy of the coastline. For our calculations we used the coastline of Wessel and Smith [1996], digitized to a spatial resolution of around 150 m. When we used a mask file with a resolution of 300 m to compute the OTL values, we found that the differences were less than 0.1 mm.
We have computed our own Green's function for the isotropic version of PREM, termed PREM-iso (used in the computation of the residuals displayed in Figure 2), by solving the differential equations of Alterman et al. [1959] using the Chebyshev collocation method [Guo et al., 2001]. In PREM the upper 3 km consists of water. We have replaced this water layer with the density and elastic properties from the underlying rock layer since our GPS stations are on land. Using this Green's function and the six ocean tide models mentioned in the previous section, the vertical M 2 OTL displacement values over western Europe, with a spatial resolution of 0.2°, were computed. Their standard deviations at each point are shown in Figure 4 and provide an indication of how errors in the ocean tide models propagate into the OTL values. Figure 4 shows large standard deviations in the Bristol Channel but also some notable differences in the Wash estuary. As was noted in section 3, the latter is mainly due to differences in GOT4.7 compared to the other models. Figure 4 also demonstrates that the standard deviation is not larger than 0.3 mm in the areas around Cornwall and Brittany, illustrating that the close agreements between the ocean tide models (few centimeters in the dominant Cornwall, Celtic Sea, and English Channel defined areas), which are in keeping with the close agreements between all ocean tide models and the tide gauge and bottom pressure data shown in section 3, also result in very close OTL agreements. Meanwhile, as suggested in section 3, the large ocean tide model errors in the Bristol Channel area have negligible effect on the OTL outside of that immediate area. Therefore, the 2-3 mm discrepancies in OTL displacements in Cornwall and Brittany presented in section 2 cannot be caused by ocean tide model errors, and given that the GPS estimation error is only 0.2-0.4 mm, the discrepancies must arise from errors in the Green's function. To investigate this, in the rest of this paper we will mainly use and show results from one ocean tide model only (FES2012).

Elastic Green's Functions
In the construction of PREM various types of observation were used such as free oscillations of the Earth, surface waves, traveltime data of body wave phases, and the mass and radius of the Earth and gravity [Dziewonski and Anderson, 1981]. In the early 1990s a new traveltime table for seismic phases was generated that led to a new seismic model called IASP91 [Kennett and Engdahl, 1991] and its successor AK135 [Kennett et al., 1995]. To be precise, both models only provide seismic velocities as a function of depth and they need to be augmented with a density model. In this case the density is, except for the crust, very similar to PREM [Kennett et al., 1995]. The AK135 and IASP91 models have a Moho depth at 35 km, while PREM's is at 24.4 km. According to Chadwick and Pharaoh [1998], the Moho depth is around 27-30 km in the southwest of England and in the north of Scotland. To investigate the influence of the crust on model OTL predictions in more detail, we have replaced the top of PREM-iso with the crust model for the southwest of England of Holder and Bott [1971]. This local seismic model of the crust only provides a depth profile of the compressional velocity V P . We have used the empirical relations of Brocher [2005] to determine the corresponding values of the density ρ and shear velocity V S that are necessary to compute the Green's function. We call this Green's function PREM-HB.
So far, we explicitly mentioned that we use the fully isotropic version of PREM, since there also exists a transversely isotropic version that for depths between 24.4 and 220 km provides different shear wave velocity values for the horizontal and vertical directions [Dziewonski and Anderson, 1981]. Following the approach of Pagiatakis [1990], we have modified our differential equations to include this effect. We call the resulting Green's function PREM-trans. Next, we consider the more recent global reference model STW105 [Kustowski et al., 2008] that is also transversely isotropic. In addition, Kustowski et al. [2008] developed a transversely isotropic seismic tomographic model for the upper mantle, called S362ANI. Using the region between longitudes 8°W and 2°W and latitudes 43°N and 51°N, covering southwest England and northwest France, we computed the mean density, compressional velocity, and shear velocity of S362ANI and used these values to compute a Green's function.
The Green's functions based on PREM-iso, AK135, PREM-HB, PREM-trans, and STW105 are shown in Figure 5. Due to its thicker crust, the AK135 Green's function differs with respect to PREM-iso at small distances up to 100 km from the point load, which is in agreement with the conclusions of Wang et al. [2012]. The same behavior is observed for PREM-HB because the crust model of Holder and Bott [1971] has lower seismic velocities than PREM-iso. The differences between PREM-trans and PREM-iso are smaller: around 2.5% at a distance of 36 km. STW105 deviates from PREM-iso at distances between 10 and 300 km, with a maximum difference of 4.5%. All the Green's functions are similar to PREM-iso for distances greater than 300-500 km.

Loading Contributions From Different Water Areas
Using the PREM-iso, AK135, PREM-HB, and STW105 Green's functions, and the six ocean tide models, we have computed the contributions to the vertical M 2 OTL displacement at the GPS station at Newlyn (NEWL) in Cornwall, southwest England, for each of the water areas defined and discussed in section 3. The results for PREM-iso are given in Table 2 and Figure 6 (top row), which also contains the GPS-observed value with its 95% confidence error bound. Figure 6 confirms that the largest contributions come from the Celtic Sea and Cornwall sea areas and that although the tides are large in the Bristol Channel (and exhibit the greatest intermodel differences of any water area), their contribution to OTL at NEWL is rather small. This further confirms that the different ocean tide models result in very similar OTL displacement values.
Figure 6 (middle row) shows similar results for the GPS station at Brest (BRST), Brittany, although there is a larger spread between the different ocean tide models especially as a result of the English Channel contribution. Figure 6 (bottom row) shows predicted and observed OTL for the station at Zimmerwald (ZIMM), Switzerland, which, due to its large distance from the ocean, has much lower OTL values, again with a small spread.

Differences Between Elastic Green's Functions Across Western Europe
The spatial distribution of the differences between the GPS-observed and predicted vertical M 2 OTL displacement for western Europe, using FES2012 and different elastic Green's functions, is shown in Figure 7. The predicted OTL values have been computed using Green's functions based on PREM-trans, AK135, and STW105. The large misfits of 2-3 mm obtained with PREM-iso in Cornwall and Brittany ( Figure 2) are still clearly visible for these Green's functions, with slight reductions in amplitude to around 1-2.5 mm. In addition, Figure 7 shows some substantial discrepancies in Scotland, southern Ireland, and along the west coast of Iberia. The plots for the other ocean tide models look similar.
The mean moduli of the vector differences shown in Figure 7, together with their minimum and maximum values and the values at NEWL and BRST, are listed in Table 3. All GPS stations have been used for these calculations except for the stations SWAS and CARI in the Bristol Channel and the station SKEE, along the east coast of England just north of the Wash estuary, which have been omitted because of their proximity to localized regions of disagreement between the ocean tide models.
Since the agreement between observed and predicted OTL values is better than 0.5 mm in most of western Europe, the mean value of the size of the differences in Table 3 is rather small (0.4-0.7 mm for the elastic Green's functions considered in this section). Nevertheless, including all stations in the statistics ensures that our selection of Green's function results in an overall improvement.
From Table 3 we can conclude that using the PREM-iso Green's function and FES2012 produces the largest residuals, whereas including the transverse isotropic effects (PREM-trans) reduces slightly the misfits at NEWL and BRST, to 2.4 and 1.7 mm, respectively. The results obtained using the Green's function based on the older Gutenberg-Bullen Earth model [Farrell, 1972] are also listed in Table 3 and are very similar to PREM-iso. Furthermore, Table 3 also demonstrates that the use of the Green's function based on PREM-HB only reduces the misfit at NEWL and BRST by about 0.4 mm compared to PREM-iso. Therefore, the influence of the crust on the predicted OTL values is small. Better results are obtained with IASP91, its successor AK135, STW105, and S362ANI. For the latter Green's function the misfits at NEWL and BRST are 1.7 and 1.1 mm, respectively.

Optimal Semiempirical Elastic Green's Function
In the previous section we demonstrated that modifying the elastic Green's function for distances of less than a few hundred kilometers from the point load improves the agreement between GPS-observed and predicted OTL values but does not explain all of the discrepancy. This corresponds roughly to modifying the elastic properties of the seismic Earth models to a depth of a few hundred kilometers, in agreement with the conclusions drawn from our simple half-space model loaded by a disc.
On the other hand, we also saw that the influence of the elastic properties of the crust was minimal; i.e., the discrepancies were not substantially reduced. This suggests that the elastic properties of the upper mantle are more critical for obtaining the most accurate predicted OTL values.
To examine this further, we have computed an optimal semiempirical elastic Green's function that minimizes the sum of the squared misfits between the modeled and observed OTL values at each GPS station. The parameters varied were the shear modulus within the asthenosphere (by reducing the shear modulus in PREM), the depth to the top of the asthenosphere (D1), and its bottom depth (D2), as shown in Figure 8. This figure also shows the depth profiles of the density, shear modulus, and bulk modulus for various Earth models. The isotropic version of PREM was taken as the reference model that has its asthenosphere between the depths of 80 and 220 km. For each set of values for these three parameters, new Love numbers and Green's functions were computed and used to predict the vertical M 2 OTL values using the FES2012 and GOT4.7 ocean tide models. Using the Nelder-Mead (also known as the downhill simplex) numerical optimization scheme, we found that both for FES2012 and GOT4.7, a reduction of 10-11% of the shear modulus in the asthenosphere gave the best agreement with the observations (mean misfit 0.3 mm, maximum 1.5 mm). The statistics obtained using these two optimal Green's functions are also included in Table 3, labeled as PREM-modFES and PREM-modGOT, respectively. The estimated values of D1 and D2 were around 50-80 and 340-350 km, respectively, which implies a thicker asthenosphere than that in PREM-iso. PREM-modFES is shown in Figure 8, and the reduction of the shear modulus in the asthenosphere, between 50 and 340 km depth, is clearly visible. We have shown that by empirically altering the properties of the asthenosphere, the discrepancies can be reduced further to 0.4-0.5 mm at NEWL and BRST, clearly demonstrating the large influence of this layer on our OTL results.

Anelasticity Effects
In section 4 we have shown that changing the elastic properties of the crust and using more recent reference Earth models, or including transversely isotropic effects, reduces the discrepancy between observed and predicted vertical M 2 OTL values, but there remain misfits of around 1-2 mm in Cornwall and Brittany. Furthermore, we have demonstrated in section 5 that if the shear modulus of the asthenosphere is reduced by 10-11%, a better agreement between predicted and observed OTL values is obtained, with the NEWL and BRST discrepancies reduced to 0.4-0.5 mm.
As was mentioned by Dziewonski and Anderson [1981], the listed seismic velocities of PREM are only valid at 1 Hz. For other frequencies the effect of dissipation, which causes a reduction of the shear modulus, needs to be taken into account. The amount of dissipation is mostly represented by the quality factor Q that is inversely related to the average energy E dissipated per cycle: Bulk dissipation Q À1 κ is much smaller than dissipation in the shear modulus Q À1 μ and will be neglected. We will assume that the dissipation takes place in an absorption band from seismic down to tidal frequencies with a  constant Q μ quality factor. Furthermore, the reference frequency is 1 Hz (ω 0 = 1), which gives us the following perturbation in the shear modulus μ at a tidal frequency ω [Dahlen and Tromp, 1998;Lambeck, 1988]: The real part of equation (3) is dependent on frequency, and this effect is called anelastic dispersion. Note that equation (3) is only valid for frequencies ω and ω 0 inside the absorption band with constant Q. The relation between stress and strain is still linear although the elastic constants are now complex numbers. Material which exhibits this behavior is called anelastic. However, it is the anelastic dispersion which causes the largest changes in the Green's functions.
In the asthenosphere, Q μ reaches a minimum of around 80. At tidal frequencies equation (3) shows that this results in a reduction of the shear modulus of about 8.5%. This value is in agreement with the value of 10% we found in section 5, where we searched for the optimal reduction of this parameter. For the other layers above the asthenosphere, Q μ has a value of around 600 and, as a result, the effect of dissipation will be less. For the layers underneath the asthenosphere Q μ is around 150, but here the depth is becoming too large to have a notable effect on the OTL values in western Europe.
Equation (3) also shows that we now have a complex-valued shear modulus [Zschau, 1978] that will produce a slight OTL phase lag. Nevertheless, the effect on anelasticity is small and at NEWL and BRST including the imaginary part of the Green's functions changes the phase by only 0.2°.
Using these modifications of the shear modulus, we have computed new complex-valued Green's functions for AK135, PREM, STW105, and S362ANI, with those for STW105 and S362ANI shown in Figure 9, again normalized  using the original PREM-iso Green's function. To indicate that the new Green's functions now refer to the period of harmonic M 2 , the label "(M 2 )" has been added to their names. As was explained in section 4.2, only the model S362ANI regional mean values (corrected using equation (3)) around western Europe were used to compute the Green's function, which is tabulated in the appendix. Figure 9 shows that for the real part (G Re ) the deviations for AK135 and STW105 with respect to PREM-iso are larger than at a period of 1 s (shown in Figure 5), and the Green's functions are about 10% larger than PREM-iso at distances of around 100 km from the point load. Figure 9 also shows that the imaginary parts are around 100 times smaller than the real part and as such, their influence on the results is small. Next, Figure 9 shows that the elastic (realvalued) PREM-modFES is in good agreement with the real part of S362ANI(M 2 ).
In considering Figure 8 again, we can see a good agreement between the shear modulus of S362ANI(M 2 ) and PREM-modFES in the asthenosphere. PREM-modFES has a deeper asthenosphere, but as we explained in section 1, our OTL values are not very sensitive to the elastic properties at those large depths.
The Green's functions PREM-iso(M 2 ), PREM-modFES, STW105(M 2 ), and S362ANI(M 2 ) have been used to compute new vertical M 2 OTL predictions that are shown in the phasor plots for NEWL and BRST in Figure 10. The phasors have slightly larger amplitudes than those for PREM-iso ( Figure 6), and this helps to close the misfit with the observations. The closest agreement with the GPS observations at NEWL and BRST is obtained using the PREM-modFES and S362ANI(M 2 ) Green's functions.
The spatial distribution of the misfit between the observed and predicted OTL values, computed using the Green's functions PREM-modFES and S362ANI(M 2 ) and the tide model FES2012, is plotted in Figure 11. Note that now the agreement over Cornwall and Brittany is excellent, and there have been substantial reductions in the misfits over Scotland and the west coast of Iberia, while the agreement has also improved slightly for the rest of western Europe.
The statistics of the results shown in Figure 11 (and those for the other (M 2 ) Green's functions not plotted) have been added to Table 3, which shows that for S362ANI(M 2 ) the misfit at NEWL and BRST has been reduced to 0.2-0.4 mm. These values are similar, although slightly better, than those obtained using the PREM-modFES Green's function.
The modifications of the elastic properties of the asthenosphere and lithosphere presented in this and previous sections affect not only the predicted OTL values but also the solid Earth tide. The amplitude of the vertical M 2 solid Earth tide in this region is around 60 mm, and our different set of elastic properties of the Earth causes a variation in this amplitude of at most 0.6 mm, with the largest variation caused by the dissipation effect. However, the solid Earth tide corrections applied in our GIPSY

Conclusions
We have used GPS to estimate the vertical OTL displacements at harmonic M 2 for 259 stations across western Europe with confidence limits of around 0.2-0.4 mm. In central Europe, far from the ocean, the agreement between predicted and observed OTL values is better than 0.5 mm when using the isotropic, elastic version of PREM. However, in Cornwall, southwest England, and in Brittany, France, where the vertical M 2 OTL displacement amplitudes can reach 45 mm, discrepancies of 2-3 mm were found. Using an intercomparison of six recent global ocean tide models (DTU10, EOT11a, FES2012, GOT4.7, HAMTIDE, and TPXO8) and a validation with tide gauge data, we show that this cannot be caused by errors in our knowledge of the ocean tides in the area. Therefore, the elastic properties of the Earth used in the OTL computations need to be modified. However, the use of various global seismic isotropic elastic Earth models such as Gutenberg-Bullen, PREM, IASP91, and AK135 does not help to explain the observed misfit significantly although the latter two models produce a better agreement. Next, replacing the crust with a local model of Holder and Bott [1971] only has a small beneficial effect of no more than 0.4 mm on the predicted OTL values.
We found that using the transversely isotropic version of PREM to compute the OTL values reduced the discrepancy by around 0.3-0.4 mm compared to its isotropic version in the areas with large OTL values. The transversely isotropic seismic reference model STW105 and the seismic tomography model S362ANI of Kustowski et al. [2008] performed better than the aforementioned Earth models, but there still remained around 1-2 mm discrepancy over Cornwall and Brittany.
Since our observations are most sensitive to the elastic properties of the asthenosphere, we estimated two optimal semiempirical Green's functions by varying the depth and thickness of the asthenosphere of the isotropic version of PREM and scaling its shear modulus. We found that a reduction of 10-11% of the shear modulus produces the best result and that the thickness of the asthenosphere had to be increased by around 100 km. This reduced the discrepancies to around 0.4-0.5 mm over Cornwall and Brittany using the FES2012 tide model, clearly demonstrating the importance of the elastic properties of the asthenosphere to our OTL values.
Next, we included the changes in elastic properties of the seismic model at tidal frequencies due to dissipation effects that cause anelastic dispersion. Using an absorption band with a constant quality factor Q from seismic to tidal frequencies, the shear modulus is reduced by 8.5% in the asthenosphere at the M 2 frequency. This is in close agreement with the results we obtained using our optimal semiempirical Green's function.
The dissipation effect was applied to the aforementioned Earth models, and new Green's functions were computed at the M 2 tidal frequency. Using the corrected S362ANI(M 2 ) Green's function, the discrepancies reduce to around 0.2-0.4 mm over Cornwall and Brittany using the FES2012 tide model, which is of the same order as the sum of the remaining errors due to errors in the ocean tide models, numerical errors in the computation of the loading values, and the uncertainty in the observations themselves. a The first column gives the angular distance in degrees. The second and third columns list the real part of the radial and tangential Green's function, multiplied by 10 12 (aθ), where a = 6371 km and θ is in radians. The fourth and fifth columns provide the corresponding imaginary parts.
Journal of Geophysical Research: Solid Earth 10.1002/2015JB011884 from the NERC BIGF (http://www.bigf.ac. uk), IGN (http://rgp.ign.fr), EUREF (http:// www.epncb.oma.be), and IGS (www.igs. org) data providers; tide gauge harmonic constants were available from the tidal data bank formerly maintained by the International Hydrographic Organization and Davies et al. [1997], and pelagic tidal constants are available from the IAPSO and GLOUP data banks. These data providers are gratefully acknowledged. We thank NASA JPL for the provision of the GIPSY software and precise orbit and clock products used and thank the developers of the ocean tide models used. Careful reviews by Duncan Agnew and Abbas Khan are acknowledged. M.S.B. was financially supported by national funds through FCT in the scope of the project IDL-FCT-UID/GEO/50019/2013 and grant number SFRH/BPD/89923/2012. Part of this work was funded by NERC grant NE/C003438/1 to P.J.C. and N.T.P., with early contributions from Maxim Keshin (now at Trimble Terrasat Gmbh) acknowledged. Some of the figures were generated using the GMT software [Wessel et al., 2013].
Journal of Geophysical Research: Solid Earth 10.1002/2015JB011884