Geodetic evidence for a buoyant mantle plume beneath the Eifel volcanic area, NW Europe

SUMMARY

The volcanism of the Eifel volcanic field (EVF), in west-central Germany, is often considered an example of hotspot volcanism given its geochemical signature and the putative mantle plume imaged underneath. EVF's setting in a stable continental area provides a rare natural laboratory to image surface deformation and test the hypothesis of there being a thermally buoyant plume. Here we use Global Positioning System (GPS) data to robustly image vertical land motion (VLM) and horizontal strain rates over most of intraplate Europe. We find a spatially coherent positive VLM anomaly over an area much larger than the EVF and with a maximum uplift of ∼1 mm yr⁻¹ at the EVF (when corrected for glacial isostatic adjustment). This rate is considerably higher than averaged over the Late-Quaternary. Over the same area that uplifts, we find significant horizontal extension surrounded by a radial pattern of shortening, a superposition that strongly suggests a common dynamic cause. Besides the Eifel, no other area in NW Europe shows significant positive VLM coupled with extensional strain rates, except for the much broader region of glacial isostatic adjustment. We refer to this 3-D deformation anomaly as the Eifel Anomaly. We also find an extensional strain rate anomaly near the Massif Central volcanic field surrounded by radial shortening, but we do not detect a significant positive VLM signal there. The fact that the Eifel Anomaly is located above the Eifel plume suggests that the plume causes the anomaly. Indeed, we show that buoyancy forces induced by the plume at the bottom of the lithosphere can explain this remarkable surface deformation. Plume-induced deformation can also explain the relatively high rate of regional seismicity, particularly along the Lower Rhine Embayment.

Satellite geodesy, Europe, Continental tectonics: extensional, Dynamics: convection currents, and mantle plumes, Intra-plate processes

1 INTRODUCTION

Intraplate volcanic activity in west-central Europe has long been associated with mantle upwellings (Granet et al. 1995; Hoernle et al. 1995). Most effort has focused on the Eifel Volcanic Field (EVF) where a period of late Quaternary volcanism (Fig. 1) continued until ∼11 ka and included the explosive eruption of Laacher See in the eastern EVF at 12.9 ka (Nowell et al. 2006; Schmincke 2007). Although the volcanism does not exhibit a clear space-time progression indicative of a hotspot track, geochemical analyses (Griesshaber et al. 1992; Hoernle et al. 1995; Aeschbach-Hertig et al. 1996; Wedepohl & Baumann 1999; Buikin et al. 2005; Bräuer et al. 2013; Caracausi et al. 2016) have shown that EVF (and some other central European) volcanic rocks and gases have the characteristics of a mantle source, while seismological studies have shown evidence for a mantle plume underneath the Eifel (Ritter 2007). Specifically, there exists a low seismic velocity anomaly down to ∼410 km depth (Ritter et al. 2001; Keyser et al. 2002; Pilidou et al. 2005; Budweg et al. 2006; Montelli et al. 2006; Zhu et al. 2012). Some studies have imaged a low seismic velocity anomaly underneath this location in the lower mantle (Goes et al. 1999; Grunewald et al. 2001; Zhao 2007), but the 660 km discontinuity seems unaffected (Budweg et al. 2006), which suggests that there is no physical connection between the lower and upper mantle anomalies. There is evidence for a broader low velocity zone at ∼50 km depth that could be interpreted as the plume head (Budweg et al. 2006; Mathar et al. 2006), which is consistent with the lithosphere-asthenosphere boundary (LAB) underneath the Eifel being relatively shallow at ∼40–50 km depth (Plomerová & Babuška 2010; Seiberlich et al. 2013). For the Massif Central (France), the other major Quaternary volcanic field in western Europe (with the most recent eruption at ∼7 ka of Lac Pavin (Juvigné & Gilot 1986; Nowell et al. 2006; Chapron et al. 2010)), some earlier tomographic studies (Granet & Trampert 1989; Granet et al. 1995; Sobolev et al. 1997) suggested an underlying plume, but only few studies have convincingly reproduced those findings (Spakman & Wortel 2004; Barth et al. 2007; Koulakov et al. 2009). Limited geochemical analyses have revealed a mantle signature there (Matthews et al. 1987; Aeschbach-Hertig et al. 1999; Zangana et al. 1999).

Figure 1.

Black dots are epicentres of earthquakes between 1000 and 2006 in the SHEEC-SHARE database (Grünthal et al. 2013; Stucchi et al. 2013). Colours are epicentre density for circular areas with radius (R) of 30 km (with each event in a circle weighted by (1–D/R), where D is distance from event to centre of circle). Blue dots are centres of Quaternary EVF activity, and blue triangle is location of latest activity in Massif Central. Light blue lines are seismogenic faults (Basili et al. 2013), and the Rhenish Massif is outlined in green. B, Belgium; NL, The Netherlands; LRE, Lower Rhine Embayment; URG, Upper Rhine Grabren. Inset: red solid/dashed polygon is extent of data and model, respectively, blue polygon is extent of results presented here, blue dashed polygon is area of Fig. 7.

Despite the above indications suggesting EVF volcanism being the result of decompression melting of a buoyant mantle upwelling, the fact that volcanism seems to flare up at the end of glacial periods (also at Massif Central) has been interpreted as volcanism instead being caused (or modulated) by glacial unloading (Nowell et al. 2006). Alternatively, the location of the EVF near the major European Cenozoic Rift System has been used to argue that volcanism is the product of fluid/gas pathways caused by Alpine-collision induced dilatancy along shear bands in the upper mantle and lower crust (Regenauer-Lieb 1998) and/or passive partial melting of the asthenospheric mantle induced by lithospheric stretching (Wilson & Downes 1992; Lustrino & Carminati 2007).

If a thermally buoyant plume is present, it is predicted to cause significant domal uplift at the surface, although the amount of uplift and spatial extent thereof depend on the stage of upwelling, temperature contrast, viscosity, width of the plume (head) (Griffiths & Campbell 1991; Hill 1991), the plume's inherent composition/density (Dannberg & Sobolev 2015), and lateral strength variations in the lithosphere induced by thermal weakening (Garcia-Castellanos et al. 2000). Vertical strength stratification within the lithosphere may even yield topographic undulations rather than a singular dome (Burov & Guillou-Frottier 2005; Burov & Gerya 2014), and one study argued this to be the case for the Eifel (Guillou-Frottier et al. 2007). By measuring the rate and spatial extent of present-day surface motion, the specifics and mere existence of a buoyant Eifel or Massif Central plume can be assessed. More generally, northwestern Europe provides a unique natural laboratory to contrast current surface deformation with plume model predictions, because other active hotspots are either underneath oceans or in areas that are tectonically and/or volcanically too active to assess the secular deformation (i.e. Yellowstone and Afar).

Hints of significant present-day uplift of the central Rhenish Massif (i.e. the area surrounding the EVF, Fig. 1) have existed for decades (Mälzer et al. 1983), but inconsistencies between different geodetic studies in the surrounding area have put into question whether measured VLM even reflected tectonic movement (Demoulin & Collignon 2000, 2002; Camelbeeck et al. 2002; Campbell et al. 2002; Francis et al. 2004). Long-running absolute gravity measurements revealed significant regional uplift at a number of stations in the Belgium part of the Rhenish Massif (i.e. Ardennes) (Van Camp et al. 2011), but that study lacked the density and broad spatial extent of measurements, nor the horizontal sensitivity, that our study provides.

The results of the 3-D surface motion in NW Europe presented here are part of a larger study encompassing most of intraplate Europe (Fig. 2). By investigating both vertical and horizontal deformation, by using robust estimations techniques (of the model and their uncertainties), and by doing this over most of intraplate Europe, our study is better suited to detect significant regional anomalies than previous studies, which had limited spatial reach and/or focused on only vertical (Kontny & Bogusz 2012; Serpelloni et al. 2013; Husson et al. 2018; Bogusz et al. 2019) or horizontal deformation (Ward 1998; Marotta et al. 2004; Nocquet et al. 2005; Tesauro et al. 2006; Bogusz et al. 2014; Keiding et al. 2015; Craig et al. 2016; Neres et al. 2018; Masson et al. 2019). This paper focuses mostly on the only significant 3-D deformation anomaly observed in our model (i.e. the ‘Eifel Anomaly’), but results for the Massif Central will be presented as well.

Figure 2.

Locations of GPS stations used, coloured by the length of the time-series. Black triangles are locations of velocities added from the literature (Kierulf et al. 2013, 2014; Lahtinen et al. 2019).

2 GPS DATA ANALYSIS

2.1 Processing details

We processed all the available geodetic-quality GPS data from continuously operating stations in our study area (red polygon inset Fig. 1). Our study benefited tremendously from the data made available to us from many regional state and commercial networks and that made this type of study possible for the first time for all of NW Europe (see Acknowledgements). The GPS data were retrieved from archives in the form of daily station-specific RINEX files, which contain raw, dual-frequency carrier phase and pseudo-range data, typically for every 15- or 30-s epoch. We then reduced the data for each station to a time-series of daily precise point positions (Zumberge et al. 1997) using the GipsyX version 1.0 software (released January 2019) licensed by the Jet Propulsion Laboratory (JPL), together with JPL's high-precision GPS orbit and clock products. Data were processed for each station independently, and thus the results are insulated from potential data problems at other individual stations. Metadata required for correct processing options were provided by the RINEX header information, corrected by custom-software and alias tables for commonly known errors such as field misplacement, typos, and non-standard receiver-type and antenna-type fields. One exception is that metadata from the International GNSS Service (IGS) was used for antenna types if available, which is the case for approximately 4 per cent of the data. Another exception is that we used our own metadata for approximate coordinates, which should be accurate to < 10 m to ensure the validity of linearized observation equations. After editing the data using the TurboEdit algorithm (Blewitt 1990) available with GipsyX, the GPS data were reduced to 300-s epochs using decimation for carrier phase, and using carrier-smoothing for pseudo-range.

In addition to the RINEX files, we used essential input data produced by JPL, including precise GPS orbits, Earth orientation, eclipse shadow events, clock offsets, and WLPB biases (Bertiger et al. 2010) for carrier phase ambiguity resolution. The satellite orbits were minimally constrained, which requires that we transform the estimated station coordinates every day by a global, seven-parameter transformation into the IGS14 reference frame, using daily values provided by JPL. For modelling ocean tidal loading, we used coefficients of the FES2004 model (Lyard et al. 2006) computed by Chalmers University, Sweden, produced by the email interface to the server described at http://holt.oso.chalmers.se/loading/. For modelling neutral atmospheric delay (commonly known as ‘tropospheric delay’, but includes stratospheric delay), we used VMF1GRID gridded map products from University of Vienna, Austria (Boehm et al. 2006), which is based on the ECMWF numeric weather model based on global meteorological data. VMF1GRID allowed us to calibrate for nominal so-called dry and wet delays, and epoch-dependent mapping function parameters for each station. For first-order ionospheric calibration, GipsyX implements the so-called ionosphere-free linear combination of observables. For higher-order ionospheric calibration, we used JPL's IONEX gridded map product available in daily files, together with the NOAA's 12^th generational magnetic field model IGRF12. The inclusion of these higher-order calibrations has been shown to substantially improve geodetic position estimates (Kedar et al. 2003; Fritsche et al. 2005; Hernández-Pajares et al. 2007). For antennas, we used calibrations made available by the IGS (Schmid et al. 2007). Following IGS standards, if a calibration was not available for a particular antenna, we did not use the data. If a calibration was not available for particular radome, we used the no-radome calibration for an antenna. Internally, the GipsyX software implements JPL's planetary ephemeris to compute tidal effects.

The estimation strategy was set by default according to the distributed version of GipsyX Version 1.0, which we now describe. Satellite positions and clocks are held to their nominal values, whereas station clocks are estimated freely as white noise every 300-s epoch. Carrier phase biases are initially estimated freely as constants in between detected cycle-slips (integer-wavelength discontinuities), then double-difference carrier phase biases to IGS stations (using the WLPB file produced daily by JPL) are constrained to integer-wavelength values on both frequencies (which are linearly combined for the ionospheric-free bias). Parameters of the neutral atmosphere are estimated as random walk processes at every 300-s epoch. These parameters include wet zenith delay (with constraints at 5 cm per square-root hour), and two gradients parameters (with constraints at 5 mm per square-root hr). The three station Cartesian coordinates were estimated as a constant over each 24 hr period. Data are processed with several iteration cycles to clean up the data using a series of post-fit residual outlier tests set by default. Our official data analysis strategy summary can be found in the Supporting Information.

2.2 Time-series analysis

The output of each station-day process includes a 24-hr estimate of a constant station position in the IGS14 frame which is obtained by using daily 7-parameter transformations produced by JPL. The IGS14 reference frame is computed by the International GNSS Service (IGS) as a GPS-compatible realization of the International Terrestrial Reference Frame ITRF2014 (Altamimi et al. 2016). For this study we consider 3-D position time-series between 2000 January 1 and 2019 October 5 , but use only stations for which the time-series span at least 2.5 yr and for which there are at least 304 data points (which translates into a theoretical minimum 33 per cent and 4.2 per cent completeness for the shortest and longest time-series considered). We estimate velocities for 2383 stations (Fig. 2), of which 1414 are in the area presented here (blue polygon inset Fig. 1).

The GPS velocities are estimated in a multistep process. First, we fit to the position time-series a station motion model that includes an intercept, trend, annual and semi-annual periodic signal, and offsets due to known equipment changes or other unknown causes. We then iteratively identify position outliers in the residual time-series and remove them from the position time-series. Outliers are being defined as a position that exceeds 3σ from the median, and σ is defined as |${1.4826 \times MAD}$|⁠, where MAD (‘median absolute deviation’) is the absolute deviation around the median, that is, median(|res(t)-median(res(t))|) (Gauss 1816; Hampel 1974). With this definition we follow the general suggestion by Leys et al. (2013) and a GPS time-series specific practice advocated by Klos et al. (2016). The 1.4826 factor is there so that σ equals the traditional standard deviation describing the variance in res(t) in case res(t) is normally distributed (Huber 1981). Our algorithm removes a median 0.76 per cent, 0.72 per cent, 0.85 per cent of the positions in the original East, North and Up time-series, respectively.

Next, the cleaned residual time-series (i.e. the time-series minus the station motion model and minus outliers) are used to construct time-series of regional common-mode components (CMC), which are removed from the original time-series. CMC removal is important because otherwise strain rates or VLM anomalies can arise from combining velocities inferred from different time-spans for which the time-series can be biased by long-period CMC signals (Márquez-Azúa & DeMets 2003, Serpelloni et al. 2013) (Supporting Information Figs S3 and S8). The CMC estimation procedure used here is described in Kreemer & Blewitt (in preparation); it uses cross-correlation-based weights, is explicitly robust against outlier data, and takes advantage of all data (including stations with short time-series). The key features are: 1) the CMC of stations with longer time-series are used to correct the shorter time-series of other stations in a hierarchical scheme, 2) as a result, all stations can potentially be ‘filter stations’ (although we identify outliers) and the precision and the spatio-temporal appropriateness of the CMC is dictated by the spatio-temporal density of the data, and 3) each step in the procedure is based on robust median statistics.

After we fit the station model once more to the ‘filtered’ time-series, and we correct the offsets, we estimate the final trend and uncertainty using the robust MIDAS algorithm (Blewitt et al. 2016). MIDAS-derived velocity uncertainties are typically similar to, and tend to be more conservative than those derived from stochastic models (Santamaría-Gómez et al. 2017; Simon et al. 2018; Murray et al. 2019).

To help constrain the far-field glacial isostatic adjustment (GIA) signal, we add 80 published velocities in southern Scandinavia (Kierulf et al. 2013, 2014; Lahtinen et al. 2019). We also exclude three stations located above abandoned coal mines in the Kempen (Belgium) and Limburg (The Netherlands) areas and which experience anomalously fast uplift due to underground groundwater flooding (Caro Cuenca et al. 2013; Vervoort 2016).

3 METHODOLOGY

3.1 VLM imaging

We interpolate the GPS station vertical velocities to obtain a robust VLM map using a method that is based on the GPS Imaging algorithm of Hammond et al. (2016), but differs in some aspects. The major difference is that we do not use a spatial structure function to determine how the weight the rates of stations as function of their distance from a station/point-of-interest. Such function implies a similar spatial coherence in vertical rates across the entire study area. This is an inappropriate assumption here, because our study area includes both potential long-wavelength (e.g. GIA (e.g. Nocquet et al. 2005; Husson et al. 2018)) and short wavelength (e.g. subsidence above the Groningen gas field (Ketelaar 2009)) signals. Instead, we aim to consider a set of stations that are as local as possible. We do this by applying a Delaunay triangulation of the station locations (Renka 1997). To reduce the influence of stations that are relatively far away, we try to find more local stations. We do this by calculating the median distance from the station/point-of-interest to all Delaunay-connected stations. We then include the rates of additional stations in the median VLM estimation when the distance to those stations is less than the median distance calculated above. In case the rate of a relatively far away station is anomalous, the addition of these extra stations would effectively down-weight the outlier (even more so than a typical median approach of the neighbouring stations would already do). Moreover, the addition of other nearby stations would make the median as local as possible and constrained by the density of stations and spatial variation therein.

In practice, when station density is rather uniform, zero to few extra stations will be identified. In our case this approach matters only for places where the station density changes laterally, such as along coast lines. It is also important to note that in the algorithm of Hammond et al. (2016) one needs to subjectively choose the weight of the station itself, while here the station itself has objectively the same weight as the other nearby stations (i.e. in case the standard deviations in the observations are the same, see below). Another important difference between our algorithm and that of Hammond et al. (2016) is how we estimate uncertainties and how we estimate the weighted median (see below)

In short, our alternative GPS Imaging approach (which we call ‘Robust Network Imaging’ (RNI) since it could be applied to any type of geospatial network data) includes the following steps. First, the vertical rate of each station i (⁠|${{x_i}}$|⁠) is replaced by a weighted median (or ‘despeckled’) value (⁠|${{{\hat x}_i}}$|⁠) derived from the rates of the N most local stations. These most local stations include the station itself, those connected to it in a Delaunay scheme, and any other station that is closer to the station itself than the median distance between the station itself and those connected to it. The weights are given by |${{\rm{1}}/{\sigma _i}}$|⁠, where |${{\sigma _i}}$| is the standard deviation in the observed rate. The weighted median algorithm is based on that most recently presented by Bowden et al. (2016), and also implemented by Kreemer & Blewitt (in preparation). For the standard deviation in |${\hat x_i}$| we use the MAD, similarly as was defined for the outlier detection.

In a second step we use these despeckled rates and standard deviation therein (⁠|${{{\hat x}_i}}$|⁠, |${\hat \sigma _i}$|⁠) to estimate the weighted median vertical rate at a gridded set of M evaluation points j within the convex hull of the station locations. A new Delaunay triangulation is calculated on a set of points that includes the evaluation point together with all station locations. The weighted median |${{{\hat x}_j}}$| is then estimated from |${{{\hat x}_i}}$| and |${\hat \sigma _i}$| of all local stations defined by the same procedure as in the first step. The standard deviation for the median at the evaluation point |${\hat \sigma _j}$| is estimated similarly as above. For this we use the absolute deviation from typically the original observed rates (i.e. x_i) which would yield conservative uncertainties. However, in a few instances (i.e. Supporting Information Fig. S2(b) and the dark shaded blue in Fig. 6) we use the absolute deviation from the despeckled values (⁠|${{{\hat x}_i}}$|⁠), which would typically yield a smaller uncertainty that reflects the spatial variation in the imaged VLM rather than in the underlying data. Either way, these robust estimates of standard deviation differ from those in the method of Hammond et al. (2016), which were either defined as the root-mean-square residual scatter in x_i, or as the weighted mean of x_i.

In practice, the Delaunay triangulation does not work if two stations are exactly collocated. We therefore allow for station coordinates and rates to be averaged when stations are within 1 m distance from each other. This is done before RNI starts. Also, when we add stations that are within a median distance, we add 100 m to that distance so that, when the median distance is in fact the distance to a station that is nominally collocated with others, the rates of those stations are considered as well.

3.2 Strain rate imaging

The horizontal strain rates are estimated with a modification of the ‘MELD’ (Median Estimation of Local Deformation) imaging algorithm of Kreemer et al. (2018). As was shown by Kreemer et al. (2018), MELD is very robust against outlier velocities, and gives realistic uncertainties. MELD-derived strain rates are derived at an evaluation point from the multivariate median of a set of strain rates from a number of station-based local triangles (see Kreemer et al. (2018)) and requires two parameters. The first one is the minimum number of triangles to use for estimating the strain rate (N_min). We use here N_min = 56, which can theoretically be reached with as few as 6 stations. The second MELD input parameter is based on the theoretical standard deviation in strain rate in a triangle of stations based on the triangle's geometry and size. The maximum allowed value (σ_max) is used to exclude triangles that are too small or skinny, which would make the strain rate estimate uncertain. We set σ_max = 11.1364, which was chosen such that noise of the amplitude of the median 1-D standard deviation in the GPS velocities (∼0.11 mm yr⁻¹) does not result in strain rates > 1 × 10⁻⁹ yr⁻¹. That is, no strain rates from a triangle of stations are considered if the standard deviation in strain rate in any one component exceeds σ_max. In general, σ_max ≈ 1.225/σ_GPS. An equilateral triangle that yields a standard deviation in strain rate of 11.1364 × 10⁻⁹ yr⁻¹ has an area with equivalent circle radius of 47.2 km, which is thus our minimum theoretical spatial resolution. The general relationship between σ_max and this minimum radius (R_min) is R_min ≈ 525.2/σ_max.

We make one change to the MELD algorithm as presented by Kreemer et al. (2018). This change is applied after the strain rate at each evaluation point (typically part of a grid of points) is estimated. Because the strain rates are based on triangles with areas equivalent to circles with radii typically larger than the distance between evaluation points, we replace the results at the evaluation points with a weighted spatial average based on all evaluation points within a distance R_median from the evaluation point considered (including the point itself). Here R_median is the corresponding radius for the median area of all triangles considered at a point. R_median is by definition ≥ R_min. The spatial smoothing not only better aligns the spatial resolution of the model results to a corresponding spatial smoothness, the smoothed field also minimizes significant model differences between two neighbouring evaluation points which could exist if each of those points are inside a different station triangle which could have yielded a significantly different set of stations to be considered. The Appendix details the estimation of standard deviation in each spatial average given that model parameters between neighbouring grid points are derived from many common velocity data yielding a large degree of correlation.

4 RESULTS

4.1 Vertical land motion

Observed VLM (Fig. 3a) for individual GPS stations shows some regional patterns, but also large variability. In a first step, we pre-process the VLM by applying our RNI algorithm. Based on this despeckled VLM field (Fig. 3b), we then use RGI again to estimate the VLM at a grid of evaluation points (Fig. 3c). Offshore areas are clipped in the figures, because chequerboard tests (Fig 4. and Supporting Information Fig. S1) reveal an expected lack of resolution there. The same tests show that at the Eifel area 200 × 200 km and 100 × 100 km VLM patches (i.e. Fig 4. and Supporting Information Fig. S1, respectively) can be significantly recovered at the 2σ, but that at the Massif Central area only the larger patches can be resolved only at the 1σ level, because of there being fewer stations there compared to the Eifel area. Because the results presented here are part of a larger study that covers most of intraplate western Europe (Fig. 2), they do not suffer from boundary effects. Numerical results are given in the Supplementary Material.

Figure 3.

(a) Observed VLM, (b) despeckled VLM, (c) imaged VLM (outline of Rhenish Massif (green polygon), centre of Eifel plume (Ritter et al. 2001) (black triangle) and Massif Central (purple triangle) are shown for reference), (d) imaged VLM corrected for GIA and shown only where corrected VLM > 2σ.

Figure 4.

(a) Input chequerboard with alternating +1 and −1 mm yr⁻¹ VLM in 200 × 200 km cells. Open/closed triangle is location of Eifel and Massif Central, (b) Imaged value at GPS locations, (c) Imaged VLM at 0.1° grid, (d) Imaged VLM where absolute value is ≥1σ of imaged VLM, (e) Imaged VLM where absolute value is ≥2σ of imaged VLM.

We observe an area of positive VLM over most of the Rhenish Massif, including the EVF, and it is centred on the central Rhenish Massif (Fig. 3c). This VLM signal is anomalous given that most of intraplate Europe south of Scandinavia is subsiding, which likely reflects forebulge collapse related to GIA. Given that GIA models can differ significantly, we choose to subtract from our results the VLM predicted by the GIA model that is most consistent with our regional observations of the rate of forebulge collapse (Husson et al. 2018). When we present our result relative to this GIA model (Fig. 3d), we find that the EVF uplift anomaly is the only coherent significant signal (at the 2σ level) in NW Europe. The area of anomalous uplift covers a roughly circular/oval area and includes most of the west-central Rhenish Massif as well as southeastern Netherlands (i.e. Limburg). The highest uplift is slightly above 1 mm yr⁻¹ relative to the regional GIA-associated subsidence, and is located right at the EVF. While VLM at the Massif Central is higher than in surrounding areas in France, the VLM signal is insignificant (i.e. at the 2σ level, but even so at the 1σ level, which suggest there either does not exist a significant VLM anomaly or the anomaly is much less than 200 × 200 km in scale (Supporting Information Fig. S1)).

Note that our study also contains a couple areas of anomalous subsidence (i.e. > 1 mm yr⁻¹). Examples are the western Paris Basin, the western part of The Netherlands, and the northernmost part of The Netherlands (Groningen). Subsidence in the western part of The Netherlands has previously been observed with InSAR and is explained by peat decomposition (Caro Cuenca & Hanssen 2008). Subsidence in Groningen has previously also been detected by InSAR (Ketelaar 2009) and is due to gas extraction (van Thienen-Visser & Breunese 2015; Jagt et al. 2017). This anthropogenically induced subsidence is the only significant negative VLM signal (at the 2σ level) in our model, with imaged subsidence up to ∼6 mm yr⁻¹ (not corrected for GIA).

4.2 Horizontal deformation

We use the velocities derived from the horizontal time-series of the same set of GPS stations used in the VLM analysis to infer the horizontal strain rate field (a subset of the velocities is shown in Fig. 5a). For all models the results are estimated on a 0.1° grid of evaluation points. For almost all areas in our preferred model the strain rates represent deformation for an area with ∼50–60 km radius (Supporting Information Fig. S4a) and this is inherently controlled by the station spacing and, more importantly, our choice to exclude strain rate estimates based on stations that are close together. The uncertainty in the dilatational strain rates is shown in Supporting Information Fig. S4(b).

Figure 5.

(a) Observed horizontal velocities relative to the extension anomaly that encompasses the Eifel (panel b). For clarity, velocities which differ > 0.25 mm yr⁻¹ from the local velocity gradient are not shown. (b) Colours and vectors show style and (normalized) principal axes of strain rate tensor (averaged over non-overlapping equal areas), respectively. Outline of Rhenish Massif and centre of Eifel and Massif Central are shown (black and purple triangles, respectively). Also shown for reference are green dashed lines, which are ellipses around Scandinavia and to which contractional strain rates in the northern part of the model are oriented orthogonally. Orange dashed circles indicates deviation from that pattern and highlights that contractional axes around the Eifel and Massif Central are oriented radially to their respective extension anomalies, (c) contours are rate of dilatation (red is extension, blue is contraction) and model velocities (based on inferred strain and rotation rates) are in same reference frame and with same scale as in (a). (d) Same contours as in (c), but only where dilatation rate is > 2σ.

Numerical results are given in the Supplementary Material, and we visualize the results in terms of style (Fig. 5b) and the dilatational amplitude of the strain rate field (Figs 5c and d). We show the dilatational strain rates, because in most places the strain rate tensor in the study area is dominated by either extension or contraction. A large part of the area is dominated by contractional strain rates, which can be explained by contraction of the GIA forebulge, as evidenced by the systematic rotation of the contraction direction in the regions directly surrounding the former ice-sheet. In the southeastern part of our model area, i.e. southern Germany, we see enhanced contraction that could be associated with Alpine shortening protruding into intraplate Europe. The pattern of the wide-scale contraction appears to be interrupted by a dilatational strain rate anomaly centred NNW of the EVF and defined by significant ∼N-S oriented extension and a maximum dilatation rate of ∼3.6 ± 0.9 × 10⁻⁹ yr⁻¹. Furthermore, the orientations of the contractional strain rates in most of the areas directly surrounding the extension anomaly exhibit a radial pattern oriented towards that anomaly. For the Massif Central we find a significant (at the 2σ level) extensional strain rate anomaly just west of the area of the most recent volcanic activity, with extension being bi-axial and a magnitude of ∼1‒2 × 10⁻⁹ yr⁻¹. The pattern of contractional strain rate around this extensional anomaly also shows a radial pattern, as seen for the Eifel, but it is less convincing.

To illustrate the robustness of our results, we also show the results for models with minimum triangle sizes at 50 per cent (i.e. minimum ∼23 km radius, Supporting Information Fig. S5) and 150 per cent (i.e. ∼71 km, Supporting Information Fig. S5) compared to our preferred model. When the spatial scale is reduced, other strain rate features appear but the only significant feature is a high-dilatation feature centred at the same place as the extensional anomaly NNW of the Eifel in our preferred model. When the spatial scale is increased, the same anomaly persists but at a lower rate than in our preferred model. We do observe extensional strain rates > 2 × 10⁻⁹ yr⁻¹ just west of recent volcanism in the Massif Central when we reduce the spatial scale (Supporting Information Fig. S5b), but for that model (as well as for the one with increased spatial scale) the dilatational strain rate anomaly is not significant at the 2σ level (Supporting Information Figs S5c and S6c).

To illustrate that our preferred MELD-derived model is robust against outlier velocities, we also present a model based on a set of velocities minus outliers and that looks very similar (Supporting Information Fig. S7) to the one presented in Fig. 5. Supporting Information Fig. S8 shows a model that is based on velocities derived from time-series that do not have CMC removed. That model shows a significant contractional area of 2–3 × 10⁻⁹ yr⁻¹ (with contraction in NS orientation) in northern France, which is also present in the model with the reduced spatial resolution (Supporting Information Fig. S5b), and this anomaly was also recently observed by Masson et al. (2019). We now show this anomaly to be mostly insignificant. Masson et al’s observation may either be an artefact of CMC in the time-series possibly not entirely removed by their first-order filtering and/or the fact that they assumed a spatial resolution that is too small (given the data noise) and an underestimation of their uncertainties (given that spatial scale). Finally, in Supporting Information Fig. S9 we show, for comparison, a result based on the original MELD algorithm that does not include the added spatial averaging and which inherently more scattered, but still reveals the major anomalies.

5 THE EIFEL ANOMALY

To further illustrate the significance of the VLM result across the Eifel, we plot the results across the profile defined in Fig. 3d (Fig. 6). The profile (corrected for GIA) highlights the broad nature and significance of positive VLM across the broader Eifel/Rhenish Massif area. We define the ‘Eifel Anomaly’ by the area where we see superimposed significant uplift and significant extensional strain rate, which is the same area where there is enhanced seismicity, the EVF, and which sits above the imaged mantle plume (Fig. 7a). The area undergoing the most significant extension is a bit offset NNW from the area of largest uplift. The centre of the uplift anomaly is, in turn, slightly offset to the north from the projected mantle plume and coincides with the EVF. Across the anomaly, the maximum horizontal separation rate is ∼0.33 mm yr⁻¹ (in a roughly NS direction), compared to the maximum GIA-corrected uplift of ∼1 mm/yr.

Figure 6.

GIA-corrected VLM for profile shown in Fig 3(d). Open circles are original VLM (inside profile box), orange squares are despeckled VLM (errors bars are 1σ), red line is imaged VLM along profile, dark and light blue outline is 1σ in imaged VLM defined, respectively, by deviation from original VLM and despeckled VLM.

Figure 7.

(a) Colours are imaged VLM corrected for GIA, vectors are principal axes of horizontal strain rate tensor (averaged over non-overlapping equal areas), purple dashed line outlines area of significant dilatation rate (at 2σ level). Focal mechanisms of regional studies (Hinzen 2003; Camelbeeck et al. 2007), colour coded for the implied extensional (red), contractional (blue), or strike-slip (green) deformation. Blue dots are centres of Quaternary EVF activity, and black lines are seismogenic faults (Basili et al. 2013). Green outline is Rhenish Massif and triangle is the projected centre of Eifel plume (Ritter et al. 2001). (b) VLM and strain rate predicted from our best-fitting plume model. Green star is centre of gain function applied to vertical force.

We do not correct the horizontal strain rate model for the effect of GIA. Available GIA models do either not present/predict horizontal motions (e.g. Husson et al. 2018; Simon et al. 2018) or their horizontal prediction is questionable, as Kreemer et al. (2018) showed to be the case for the ICE6G_C(VM5a) model (Argus et al. 2014; Peltier et al. 2015) in North America. From looking elsewhere in our model we assess that there could be ∼1 × 10⁻⁹ yr⁻¹ NNE-oriented contraction related to GIA at the latitude of the Eifel Anomaly. So, if we would have attempted a GIA correction, it would increase the extensional strain rate at the Eifel anomaly by that amount. We furthermore note that in the plume modelling, discussed below, we do account for gradients in strain rates, which could be there partly due to GIA.

6 PLUME MODELLING

We use a simple model to test if a buoyant plume can explain the long-wavelength 3-D surface deformation. A rising plume will exert dynamic and buoyancy forces on the elastic lithosphere and we focus here on the buoyancy forces. In order to mimic those forces coming from a mantle plume head, we model them as a bi-modal Gaussian areal distribution of half-space vertical forces (Mindlin 1936; Anderson 1937) exerted on a plane at depth (corresponding to the elastic lithosphere). The relationship between surface deformation is schematically given as: [u,v,w] = A(x,y:σ₁,σ₂,θ,c) × F(u,v,w:x,y,c). The 3-D surface deformation is represented here as displacement rates (u,v,w, being east, north and vertical velocities, respectively), but we also invert for strain rates. The gain function, A(x,y), is a normalized bi-modal Gaussian amplitude distribution with standard deviations (σ₁,σ₂), centred on x₀,y₀, with a rotation relative the East axis (θ), and a fixed depth of c = 50 km, and F is the vertical force. This 6-parameter model is fit to the measured displacement and strain rates, with regional trends removed. Removing the regional trends involves linear fits that added another 9 parameters, for a total of 21. The gain function has a maximum value of 1 at its origin and is evaluated on a grid of 40 km spacing that spans the area (515 × 422 km). Given it is a coarse model, smoothing of the horizontal displacements is necessary, and the resultant strains. A smoothing boxcar of 110 km was used on both horizontal velocities and strain data and corresponding model values, but not on the vertical data. The kinematic indicators are evaluated on an ∼10 km grid. The associated matlab program is made in available in the Supplementary Material.

For our best-fitting model, we find the half-widths of the bi-modal Gaussian distribution to be 174 and 98 km in the roughly EW and NS directions, respectively, and the centre is found at 50.5° N 6.4° E. Table 1 lists all the model parameters, and Fig. 8 graphically shows function A. Our simple model adequately fits the long-wavelength measured displacement and strain rates (Fig. 7b and Supporting Information Fig. S10). We constrain the plume head to be at 50 km depth (consistent with imaged depth of LAB and plume head (Budweg et al. 2006; Mathar et al. 2006)), because of trade-offs between the Gaussian widths and the depth of the model. Supporting Information Fig. S11 shows the predicted surface deformation for the best-fitting model that forces the depth to be at 27 and 100 km.

Figure 8.

Graphical representation of gain function A(x,y). See the text for explanation.

Table 1.

Plume model parameters (see the text for explanations).

Parameter	Value	st. dev.
Vertical Force F (1e11 N)	15.6438	0.1194
σ1 (km)	186.5662	1.6380
σ2 (km)	105.4722	1.2793
Lon0 (° E)	6.5000	0.0101
Lat0 (° N)	50.5000	0.0050
θ (radians CCW from East)	0.2542	0.0082
u0 (regional offset, mm)	0.3502	0.0064
v0 (regional offset, mm)	−0.0113	0.0001
w0 (regional offset, mm)	−0.0016	0.0002
du/dE (regional slope, mm km−1)	0.2421	0.0069
du/dN (regional slope, mm km−1)	−0.7635	0.0124
dv/dE (regional slope, mm km−1)	−0.0001	0.0002
dv/dN (regional slope, mm km−1)	0.0048	0.0002
dw/dE (regional slope, mm km−1)	−0.0130	0.0002
dw/dN (regional slope, mm km−1)	−0.0103	0.0002
dexx/dE (regional slope, mm km−1)	0.0136	0.0005
dexx/dN (regional slope, mm km−1)	−0.0230	0.0005
deyy/dE (regional slope, mm km−1)	−0.0158	0.0004
deyy/dN (regional slope, mm km−1)	0.0297	0.0006
dexy/dE (regional slope, mm km−1)	−0.0083	0.0004
dexy/dN (regional slope, mm km−1)	0.0026	0.0005

7 DISCUSSION AND CONCLUSIONS

The remarkable superimposition of significant uplift, horizontal extension, and volcanism in the Eifel area strongly suggests a causal relationship with the putative underlying mantle plume. The circular VLM pattern is consistent with the Quaternary uplift but at odds with studies, some specific to the Eifel area, that predict an undulating pattern (Burov & Guillou-Frottier 2005; Guillou-Frottier et al. 2007; Burov & Gerya 2014). Those predicted undulations are at a shorter wavelength than when just having a central dome, but if those undulations were present we would have expected to see those in the VLM field as we can resolve 100 km wide anomalies around the Eifel (Supporting Information Fig. S1).

To first order, most model studies would predict to find uplift (Campbell 2005; Dannberg & Sobolev 2015) and extension (Burov et al. 2007; Cloetingh et al. 2013) above a buoyant mantle plume. Indeed, we obtain a good regional fit to the long-wavelength aspects of the surface deformation by applying buoyancy forces related to the plume head at the bottom of the lithosphere. This is the simplest model consistent with seismic evidence, but it should be noted that an actual inversion for the depth is very poorly constrained. For example, placing the force distribution at crustal level could fit our observations too (Supplementary Supporting Information Fig. S11). There is, however, no evidence of a regional magmatic sill, although seismic velocities in the lowermost crust underneath most of the Rhenish Massif are found to be significantly reduced (Prodehl et al. 2006). We also note that the surface deformation contains some details (such as some sharp edges in the uplift anomaly) which our simple model does not fit. The fit may improve when considering lateral variations in the strength of the lithosphere as caused by plume-induced thermal weakening, which would also require a shallower source (Garcia-Castellanos et al. 2000). Furthermore, the asymmetry of the deformation pattern and its offset from the imaged plume could be a consequence of the plume being tilted and/or the interaction between a rising plume and a moving plate (Wüllner et al. 2006).

Based on river terrace elevation data, the Rhenish Massif is known to have uplifted since ∼700–800 ka (i.e. the same time as Quaternary Eifel volcanism commenced), with maximum uplift between ∼140 m (Westaway 2001; Demoulin & Hallot 2009) and 250 m (Van Balen et al. 2000; Meyer & Stets 2002) centred on the Eifel area, where we see the highest uplift. These data imply ∼0.1‒0.3 mm yr⁻¹ of average uplift since ∼700‒800 ka. Such rate is considerably lower than we find here and could either be evidence that uplift has increased since the onset of the volcanism in the late Quaternary or that the loading/unloading effect related to glacial periods causes the net VLM to vary considerably over time.

From our spatially integrated force and the first-order assumption that the plume has effectively been buoyant since between 250 ka (to explain 250 m Quaternary uplift) and 800 ka (at today's rate) ago, we estimate that a 360 km (i.e. 410 minus 50 km) high plume requires a ∼57‒184 kg m⁻³ density reduction (i.e. ∼0.7‒5.6 per cent of a 3300 kg m⁻³ dense reference mantle), which is consistent with observed seismic velocity reductions (Ritter 2007).

Although the EVF is found in what is typically referred to as intraplate Europe, the area near the EVF is seismically unusually active (Fig. 1). Much of that seismicity is attributed to faulting within the Lower Rhine Embayment (LRE) (Hinzen & Reamer 2007), and those faults appear to have increased their slip rate since the same time as the onset of late Quaternary uplift and volcanism (∼700‒800 ka) (Gold et al. 2017). In this area of highest seismic activity we also find the highest extension rates. We particularly note that the roughly N‒S oriented extension we find above the LRE is favourably orientated to generate extensional earthquakes along the normal faults within the LRE, as evidenced by the orientation of the mostly extensional earthquakes there (Plenefisch & Bonjer 1997; Hinzen 2003; Camelbeeck et al. 2007) (Fig. 7a). Our findings suggest that the surface deformation imposed by the Eifel plume explains why the LRE is so much more seismically active than many of the faults associated with other failed rifts in Europe (such as the Upper Rhine Graben (URG) between the Eifel and Alps). In fact, we see no deformation anomaly (horizontal or vertical) along the URG that could be interpreted as localized extension; a finding that is consistent with some prior studies (Rózsa et al. 2005; Tesauro et al. 2005), but inconsistent with others (Campbell et al. 2002; Fuhrmann et al. 2013, 2015).

Recently, it was also found that low-frequency seismic swarms occur in the lower crust underneath the Laacher See (Hensch et al. 2019). This activity was interpreted as the vertical migration of magma or magmatic fluids. Those findings, when combined with ours (as well as observations of continual gas emissions (Griesshaber et al. 1992; Aeschbach-Hertig et al. 1996; Buikin et al. 2005; Bräuer et al. 2013; Caracausi et al. 2016) and degassing events being correlated with seismicity (Berberich et al. 2019)), strongly suggest that the EVF is an active dynamic system.

SUPPLEMENTARY INFORMATION

supplementary-material.pdf

Figure S1. (a) Input chequerboard with alternating +1 and –1 mm yr⁻¹ VLM in 100 × 100 km cells. Open/closed triangle is location of Eifel and Massif Central; (b) imaged value at GPS locations; (c) imaged VLM at 0.1° grid; (d) imaged VLM where absolute value is ≥1σ of imaged VLM; (e) Imaged VLM where absolute value is ≥ 2σ of imaged VLM.

Figure S2. (a) Standard deviation defined relative to observed station VLM. (b) Standard deviation defined relative to despeckled VLM.

Figure S3. (a–d) Same as Fig. 3 and (e) same as Fig. 6 but all results based on VLM derived from time-series not filtered for common-mode errors

Figure S4. (a) Spatial scale of strain rate estimate expressed as the radius of a circle whose area is the median area of all triangles considered in the strain rate estimation. (b) Standard deviation in the dilatation rate.

Figure S5. Same as Fig. 5(b)–(d) but with a spatial scale that is 50 per cent of preferred model.

Figure S6. Same as Fig. 5(b)–(d) but with a spatial scale that is 150 per cent of preferred model.

Figure S7. Same as Fig. 5(b)–(d) but with outlier velocities (>1.5 mm yr⁻¹) removed before modelling.

Figure S8. Same as Fig. 5(b)–(d) but based on velocities derived from time-series that had no common-mode errors removed.

Figure S9. This result is based on the MELD algorithm as originally presented (Kreemer et al. 2018).

Figure S10. Left: observed (stars) and modelled vertical and horizontal velocities middle of the model in E–W direction: blue is vertical, red is NS and black is EW velocities, respectively. Right: Same in N–S direction.

Figure S11. Same as Fig. 7(b) but for different depths: (a) 27 km and (b) 100 km

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

We are extremely grateful to the many agencies, companies and networks that have made GPS RINEX data available. In many cases, the data were provided specifically to make our study possible. We specifically thank the following networks and institutions for raw RINEX data: EUREF (Bruyninx et al. 2019), IGS and the many agencies working under its umbrella (Johnston et al. 2017), CDDIS (Noll 2010), Technical University Delft (The Netherlands), NETPOS (The Netherlands), 06-GPS (The Netherlands), LNR Globalcom (The Netherlands), NAM/Shell (The Netherlands), ESTEC (The Netherlands), Royal Observatory of Belgium (Bruyninx & Defraigne 2018), FLEPOS (Belgium), WALCORS (Belgium), Institut Géographique National (France), Réseau GNSS Permanent (France), ORPHEON (France), RENAG (France) (RESIF 2017), SONEL (France) (Gravelle et al. 2019), SPSLux (Luxembourg), Natural Environment Research Council's British Isles continuous GNSS Facility (United Kingdom), Leica Smartnet (Poland), ASG-EUPOS (Poland), Pecny Observatory (Czech Republic), GEONAS (Czech Republic) (Schenk et al. 2010), VESOG (Czech Republic) (Kostelecký & Kostelecký 2005), CZEPOS (Czech Republic) (Kostelecký & Kostelecký 2005), FreDNet (Italy), Regione Autonoma Friuli Venezia Giulia (Italy), STPOS Bolzano (Italy), Politecnico di Torino (Italy), ARPA Piemonte (Italy), Rete GPS Veneto (Italy), Rete Dinamica Nazionale (Italy), Instituto Nazionale di Geofisica e Vulcanologia (Italy), ITALPOS (Italy), Geodetski Institute (Slovenia), NIEP (Romania), ROMPOS (Romania), ESTPOS (Metsar et al. 2018) (Estonia), LATPOS (Latvia), LITPOS (Lithuania), Danish GPS Center Aalborg (Denmark), Danish Geodata Agency (Denmark), Instituto Geográfico Nacional (Spain), Topcon (Spain), CATNET Catalonia (Spain), GNSS Activa del Principado de Asturias (Spain), Red Activa de Estaciones GNSS de Cantabria (Spain), Red de estaciones GNSS de Castilla y León (Spain), Red Extremeña de Posicionamiento (Spain), Carto Galicia (Spain), REGAM Murcia (Spain), Meristenum Murcia (Spain), Government of La Rioja (Spain), Universidad de Oviedo (Spain), Ayuntamiento de Leganes (Spain), ETSI Topografia Geodesia y Cartografia (Spain), Leica Smartnet (Poland), ASG-EUPOS (Poland), ReNEP—Instituto Geografico Portugues (Portugal), SERVIR (Portugal), Norwegian Mapping Authority (Norway), Norwegian Metrology Service (Norway), SWEPOS (Sweden), Finnish Geodetic Institute (Finland), LIPOS (Austria), SIGNAL (Slovenia), Bundesamt für Kartographie und Geodäsie (Germany), GFZ German Research Centre for Geosciences (Germany) (Uhlemann et al. 2016), Bundesanstalt für Gewässerkunde (Germany), Deutsches Geodätisches Forschungsinstitut (Germany) (Seitz et al. 2014), ASCOS (Germany) and SAPOS networks operated by various German States (i.e. State Office for Spatial Information and Land Development Baden-Württemberg, Hessian State Office for Land Management and Geoinformation, Landesamt für innere Verwaltung Mecklenburg-Vorpommern, Landesvermessung und Geobasisinformation Brandenburg, Staatsbetrieb Geobasisinformation und Vermessung Sachsen, Landesamt für Vermessung und Geoinformation Thüringen, Landesamt für Vermessung und Geobasis Information Rheinland-Pfalz, Bezirksregierung Köln, Geoinformation und Landentwicklung Saarland, and Landesamt für Digitalisierung, Breitband und Vermessung Bayern). ORPHEON GNSS data were provided to the authors for scientific use in the framework of the GEODATA-INSU-CNRS convention. The services of the UK Natural Environment Research Council (NERC) British Isles continuous GNSS Facility (BIGF), www.bigf.ac.uk, in providing archived GNSS data for this station to this project, are gratefully acknowledged. A special thanks to H. van der Marel for helping to make the Kadaster/06-GPS/NAM data in the Netherlands available, A. Araszkiewicz for helping to make the Polish Smartnet data available and T. van Dam for helping to make the SPSLux data in Luxembourg available. The final GPS velocities and GPS position time-series used in the analysis (filtered and unfiltered, and corrected for offsets) can be found at https://doi.org/10.7910/DVN/ONATFP, and the original time-series at http://geodesy.unr.edu (Blewitt et al. 2018). Model results are given in the Supporting Information. We thank J. Harvey for the locations of the volcanic centres and I. Zaliapin for his help with the material presented in the Appendix. All figures were made with the Generic Mapping Tools (Wessel et al. 2013). This work got started through financial support to CK from the Royal Dutch Academy of Sciences. Additional support for CK was provided by United States Geological Survey (USGS) National Earthquake Hazard Reduction Program (NEHRP) award G18AP00019 and for CK and GB by National Aeronautics and Space Agency (NASA) grants NNX16AK89G and 80NSSC19K1044. We thank J. Freymueller and and anonymous reviewer for constructive reviews.

APPENDIX: WEIGHTED AVERAGE AND ITS VARIANCE

For each evaluation point we wish to obtain the weighted average (W) of the N strain rate estimates m_i, i = 1,…, N, consisting of the estimate at the evaluation point itself and those at other evaluation points within a certain distance. The weights (w_i) are defined as |${\sigma _{ii}^{ - {\rm{1}}}}$|⁠, normalized by |${\sum\limits_{i = {\rm{1}}}^N {{w_i}} }$|⁠, and |${{\sigma _{ii}} = \sigma _i^{\rm{2}}}$| with σ_i being the standard deviation in m_i. The expected value of W is given by

$$\begin{eqnarray*} {E\left[ W \right] = \sum\limits_{i = {\rm{1}}}^N {{w_i}} {m_i}} \end{eqnarray*}$$

(A1)

and the variance in W by

$$\begin{eqnarray*} {\mathrm{ Var}\left[ W \right] = \sum\limits_{i = {\rm{1}}}^N {w_i^{\rm{2}}} {\sigma _{ii}} + {\rm{2}}\sum\limits_{i \ne j} {{w_i}} {w_j}{\sigma _{ij}}}. \end{eqnarray*}$$

(A2)

The question is what is σ_ij, that is, the covariance between m_i and m_j? In our case, each estimate of m_i and m_j can be thought of as an average of strain rates (although it is really part of a multivariate median) for T_i and T_j number of triangles, respectively, and those are based on GPS velocities at P_i and P_j points, respectively. Unlike the triangles, the sets of used GPS points provide a straightforward measure of interdependency between m_i and m_j. We therefore propose the following definition of the covariance between two estimates (E₁ and E₂) based on P₁ and P₂ observations, respectively,

$$\begin{eqnarray*} {\mathrm{ Cov}\left[ {{E_{\rm{1}}},{E_{\rm{2}}}} \right] \approx \frac{{{p_{\rm{0}}}}}{{\rm{2}}}\left( {\frac{{{P_{\rm{1}}}\mathrm{ Var}\left( {{E_{\rm{1}}}} \right) + {P_{\rm{2}}}\mathrm{ Var}\left( {{E_{\rm{2}}}} \right)}}{{{P_{\rm{1}}}{P_{\rm{2}}}}}} \right)} \end{eqnarray*}$$

(A3)

where p₀ is the number of common GPS stations that were used in both estimations.

The proposed covariance (A3) has the following natural properties:

1. If the estimations do not use common stations, they are uncorrelated: i.e. p₀ = 0 and thus Corr(E₁,E₂) = 0.

2. If the estimations use the same stations, they are fully correlated: i.e. p₀ = P₁ = P₂. and thus Corr(E₁,E₂) = 1.

3. If the two estimations have very different variances (and the station number is approximately the same), the examined covariance is mostly affected by the largest variance.

The approximation in (A3) is based on the following general and conceptual case. Consider a set of independent random variables |${{X_{\rm{1}}},....,{X_{{p_{\rm{0}}} + {p_{\rm{1}}} + {p_{\rm{2}}}}}}$|⁠. Assume that the estimation E₁ at the first location uses P₁ = p₀ + p₁ observations, and the estimation E₂ at the first location uses P₂ = p₀ + p₂ observations in such a way that p₀ observations (and only those ones) are used at both locations. Specifically, we may assume that both locations use the common observations

$$\begin{equation*} {\left\{ {{X_{\rm{1}}},....,{X_{{p_{\rm{0}}}}}} \right\}}. \end{equation*}$$

In addition, the first location used the unique observations

$$\begin{equation*} {\left\{ {{X_{{p_{\rm{0}}} + {\rm{1}}}},....,{X_{{p_{\rm{0}}} + {p_{\rm{1}}}}}} \right\}} \end{equation*}$$

and the second location uses the unique observations

$$\begin{equation*} {\left\{ {{X_{{p_{\rm{0}}} + {p_{\rm{1}}} + {\rm{1}}}},....,{X_{{p_{\rm{0}}} + {p_{\rm{1}}} + {p_{\rm{2}}}}}} \right\}}. \end{equation*}$$

For convenience, we denote the observations used at the first locations Y_i, i = 1,…, P₁, and the observations used at the second location Z_i, i = 1,…, P₂, in such a way that Y_i = Z_i = X_i for all i = 1,…, p₀. Let us assume that the examined estimations are the sample averages:

$$\begin{equation*} {{E_{\rm{1}}} \equiv {{\bar X}_{\rm{1}}} = \frac{{\rm{1}}}{{{P_{\rm{1}}}}}\sum\limits_{i = {\rm{1}}}^{{P_{\rm{1}}}} {{Y_i}} }\quad {{E_{\rm{2}}} \equiv {{\bar X}_{\rm{2}}} = \frac{{\rm{1}}}{{{P_{\rm{2}}}}}\sum\limits_{i = {\rm{1}}}^{{P_{\rm{2}}}} {{Z_i}} } \end{equation*}$$

The covariance between |${{{\bar X}_{\rm{1}}}}$| and |${{{\bar X}_{\rm{2}}}}$| is then given by:

$$\begin{equation*} \begin{array}{@{}*{2}{l}@{}} {Cov\left( {{{\bar X}_{\rm{1}}},{{\bar X}_{\rm{2}}}} \right)}&{{ = \frac{{\rm{1}}}{{{P_{\rm{1}}}{P_{\rm{2}}}}}Cov\left( {\sum\limits_{i = {\rm{1}}}^{{P_{\rm{1}}}} {{Y_i}} ,\sum\limits_{i = {\rm{1}}}^{{P_{\rm{2}}}} {{Z_i}} } \right)}}\\ {}&{ = \frac{{\rm{1}}}{{{P_{\rm{1}}}{P_{\rm{2}}}}}Cov\left( {\sum\limits_{i = {\rm{1}}}^{{n_{\rm{0}}}} {{Y_i}} + \sum\limits_{i = {n_{\rm{0}}} + {\rm{1}}}^{{P_{\rm{1}}}} {{Y_i}} ,\sum\limits_{i = {\rm{1}}}^{{n_{\rm{0}}}} {{Z_i}} + \sum\limits_{i = {n_{\rm{0}}} + {\rm{1}}}^{{P_{\rm{2}}}} {{Z_i}} } \right)}\\ {}&{ = \frac{{\rm{1}}}{{{P_{\rm{1}}}{P_{\rm{2}}}}}Cov\left( {\sum\limits_{i = {\rm{1}}}^{{n_{\rm{0}}}} {{X_i}} + \sum\limits_{i = {n_{\rm{0}}} + {\rm{1}}}^{{P_{\rm{1}}}} {{Y_i}} ,\sum\limits_{i = {\rm{1}}}^{{n_{\rm{0}}}} {{X_i}} + \sum\limits_{i = {n_{\rm{0}}} + {\rm{1}}}^{{P_{\rm{2}}}} {{Z_i}} } \right)}\\ {}&{ = \frac{{\rm{1}}}{{{P_{\rm{1}}}{P_{\rm{2}}}}}Cov\left( {\sum\limits_{i = {\rm{1}}}^{{n_{\rm{0}}}} {{X_i}} ,\sum\limits_{i = {\rm{1}}}^{{n_{\rm{0}}}} {{X_i}} } \right)}\\ {}&{ = \frac{{\rm{1}}}{{{P_{\rm{1}}}{P_{\rm{2}}}}}Var\left( {\sum\limits_{i = {\rm{1}}}^{{n_{\rm{0}}}} {{X_i}} } \right) = {\sigma ^{\rm{2}}}\frac{{{p_{\rm{0}}}}}{{{P_{\rm{1}}}{P_{\rm{2}}}}}} \end{array}{} \end{equation*}$$

Observe that

$$\begin{equation*} {{\sigma ^{\rm{2}}} = Var\left( {{{\bar X}_{\rm{1}}}} \right){P_{\rm{1}}} = Var\left( {{{\bar X}_{\rm{2}}}} \right){P_{\rm{2}}}} \end{equation*}$$

and accordingly

$$\begin{equation*} {\mathrm{ Cov}\left( {{{\bar X}_{\rm{1}}},{{\bar X}_{\rm{2}}}} \right) = Var\left( {{{\bar X}_{\rm{1}}}} \right)\frac{{{p_{\rm{0}}}}}{{{P_{\rm{2}}}}} = \mathrm{ Var}\left( {{{\bar X}_{\rm{2}}}} \right)\frac{{{p_{\rm{0}}}}}{{{P_{\rm{1}}}}}}. \end{equation*}$$

However, the assumption of equal variance σ² might be violated. We hence suggest, as a simple heuristic device, to estimate the covariance as the following average:

$$\begin{equation*} {\mathrm{ Cov}\left[ {{{\bar X}_1},{{\bar X}_2}} \right] \approx \frac{1}{2}\left( {\mathrm{ Var}\left( {{{\bar X}_1}} \right)\frac{{{p_0}}}{{{P_2}}} + \mathrm{ Var}\left( {{{\bar X}_2}} \right)\frac{{{p_0}}}{{{P_1}}}} \right)} \end{equation*}$$

which is equivalent to (A3).