1 Introduction

In this study, we apply the Bayesian data assimilation model on updating the global GIA process by using GPS data. We use the solution from the ICE6G-VM5 model as our prior mean for the GIA process. The GPS data are processed into yearly vertical bedrock movement from a selected global network.

Here we first assume the updating process is stationary on the sphere (Stationary model); then we remove the stationarity constraint by using the information of expert defined pseudo polygons where the GIA should be zero in theory. Three different approaches are used in modelling the non-stationarity:

  1. Use pseudo-polygon observations — Polygon observation

  2. Use pseudo-polygons to define the subset where GIA can be modeled as stationary process — Subset model

  3. Use pseudo-polygons to define a mixture Gaussian process where the updating process has different spatial properties inside and outside of the pseudo-polygons — Mixture model

Detailed description of each approach and codes of how they are implemented in R-INLA can be through the above links. In the following, we summarize major assumptions and compare the results.

2 BHM and assumptions

In this framework, we define the GIA process as a continuous spatial process on the sphere and define it by \(X\). Meanwhile, we have physical GIA model simulations \(m\) as our prior mean of the process. The actual value of the process can be measured by the GPS observations \(Y\) with measurement errors \(\varepsilon\) that are mutually independent. We use the GPS measurement to update the deviance between the prior and unobserved true process by through the linking the process to the observations \[Y = AX + \varepsilon\].

By removing the prior mean from both side we have \[\tilde{Y} = A\tilde{X} + \epsilon\] where \(\tilde{Y} = Y- Am\) and \(\tilde{X} = X - m\). Then it is reasonable to assume that the discrepancy process \(\tilde{X}\) is a stationary spatial Gaussian process on the sphere with hyper parameter \(\theta\). If the forward physical model fit the true process well, we should expect the discrepancy process to capture only small scale variability and the spatial correlations range should not be long. Hence, the BHM can be written as \[ \tilde{Y} | \tilde{X}, \theta \sim \mathcal{N} (A\tilde{X}, Q_Y^{-1})\\ \tilde{X} |\theta \sim \mathcal{GP}(0, \theta)\\ \theta \sim \pi(\theta)\]

The purpose is to produce a map of the updated GIA mean and the corresponding uncertainty at one degree resolution.

2.1 GMRF approximation and generating mesh

The BHM as defined above can be computationally expensive (even not feasible) when the dimension of \(\tilde{X}\) is large in generating \(\tilde{Y}\). We follows the SPDE approach to approximate \(\tilde{X}\) by a GMRF and use INLA for approximate Bayesian inference.

For the SPDE approach, we need to represent the process by a GMRF defined on the vertices of a triangulation (mesh) of the study domain. General guide lines about choosing mesh can be found here. In this experiment, we choose the mesh to be approximately regular and evenly distributed on the sphere with about 1 degree resolution.

The semi-regular mesh is generated by using Fibonacci points as initial locations. The Fibonacci points are approximately evenly distributed on the sphere. The goal is to produce a map with approximately one degree resolution, hence we also need a similar resolution for the triangulation. Note that the triangle edge be no larger than the process correlation length. A one degree grid is about \(110 \times 110 km\) which corresponds to a correlation length about \(110/6371 = 0.017\). The number of triangles should be \(f = 180 * 360 = 64800\) corresponding to number of vertices of \(v \approx (f+2)*2/5= 25921\).

## Generate mesh
fibo_points <- fiboSphere(N = 12960, L0 = TRUE)
fibo_points_xyz <- do.call(cbind, Lll2xyz(lat = fibo_points[,2], lon = fibo_points[,1]))
mesh1 <- inla.mesh.2d(loc = fibo_points_xyz, cutoff = 0.01, max.edge = 0.5)
summary(mesh1)

2.2 Pseudo Polygons

In fact the discrepancies process is non-stationry on the sphere. For certain regions we are more sure about the GIA value that they should be zero with high chance. Therefore we define the “zero regions” where the process should constrained to be zero. Below we show the final polygons representing the zero regions and details can be found here. The guide lines are:

  1. Ensemble mean of 8 GIA solutions smaller than 0.3

  2. Ensemble standard error of the 8 GIA solutions smaller than 0.4

  3. Area of the region larger than 200

We can using a single polygon observation to update the process while still assuming the process to be stationary on the sphere. We call this approach the polygon observation model. Implementation can be found here.

2.3 Non-stationarity

We take two approach for modelling non-stationarity.

2.3.1 Modelling subset

Here we only model the GIA process on the subset of sphere that is outside the zero region. We call this the subset model and details can be found here

2.3.2 mixture models

Here we assume the process have different spatial correlation length within and outside the zero region. The mixture model is similar to the barrier model. Details of the implementation can be found here.

One major difficulty is to separate the mesh in the subset of interest from the zero region. This done by finding the triangles in the mesh inside the interested region. Details can also be found in the above link and below we show the mesh that differentiate the triangles in the two regions.

2.3.3 pseduo observations

For the mixture model, we also need to define some pseudo-observations to constraint zero region. We simple choose these observations to be at the triangle vertices and all equal to zero with measurement error equal to 0.1 which is smaller than the smallest of the real observation.

3 Comparing results

We summarize the results from the four approaches in the following.

3.1 Hyper parameters

Below shows the posterior densities of the hyper parameters \(\sigma^2\) and \(\rho\) from the 4 models. The marginal variances reflect the marginal variability of a single vertices and the correlation lengths \(\rho\) describes the distance at which two points are almost independent to each other. The refined model shows smaller variability and short correlation length estimate of the discrepancy process.

Below shows the predicted GIA mean from the four models. For the subset model, the zero region are masked in grey which means there is no prediction available.

Below shows the predicted GIA uncertainty from the four models.

Below shows the difference between the predicted GIA and the prior from the four models.