1 Introduction

In this study, we apply the Bayesian data assimilation model on updatding the global GIA process by using GPS data. We use the soluton 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; then we remove the stationarity constraint by using the information of expert defined pseudo polygons where the GIA should be zero in theroy. 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 modelled as stationary process — Subset model

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

In the following, we demonstrate the global statinary approach step by step.

2 Data and prior

First we setup the data and prior. Note that the ICE6G solution is given on a 1 degree grid. For any given location, we define the GIA prior to be the value of the grid in which the location falls.

2.1 GIA forward model solution

Load the GIA forward solution and define it as a spatial polygon data set.

library(sp); library(INLA); library(GEOmap)
library(ggplot2); library(grid); library(gridExtra)
source("functions.R")
ice6g <- read.table("ice6g.txt", header = T)
polycoords <- ice6g[,c(6:13, 6,7)] 
plist <- lapply(ice6g$ID, 
                function(x) Polygons(list(Polygon(cbind(lon = as.numeric(polycoords[x, c(1,3,5,7,9)]), 
                                                        lat = as.numeric(polycoords[x, c(2,4,6,8,10)])))), ID = x))
Plist <- SpatialPolygons(plist, proj4string = CRS("+proj=longlat"))
## Note that in this data set the grid boundaries are defined on [-0.5, - 359.5] and the centre points on [0, 359].
## Corresponding transformation is needed in the following in geometery operations.

The GIA process will be approximated on a triangular mesh with piece-wise linear basis functions, we will only use the GIA values at the triangle vertices to represent the processes. The following generate a semi-regular mesh on a sphere with approximately 1 degree resolution and we define the GIA prior to be GIA values at the mesh nodes. The semi-regular mesh is generated by using Fibonacci points as intial 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\).

fibo_points <- fiboSphere(N = 3e4, LL = FALSE)
mesh <- inla.mesh.2d(loc = fibo_points, cutoff = 0.01, max.edge = 0.5)
summary(mesh)
## 
## Manifold:    S2
## CRS: N/A
## Vertices:    30000
## Triangles:   59996
## Boundary segm.:  0
## Interior segm.:  0
## xlim:    -0.999934 0.999934
## ylim:    -0.9999457 0.9998971
## zlim:    -0.9999667 0.9999667

Now find the GIA value at the mesh nodes.

meshLL <- Lxyz2ll(list(x=mesh$loc[,1], y = mesh$loc[,2], z = mesh$loc[,3]))
meshLL$lon <- ifelse(meshLL$lon >= -0.5, meshLL$lon,meshLL$lon + 360)
mesh_sp <- SpatialPoints(data.frame(lon = meshLL$lon, lat = meshLL$lat), proj4string = CRS("+proj=longlat")) 
mesh_idx <- over(mesh_sp, Plist)
GIA_prior <- ice6g$trend[mesh_idx]

2.2 GPS data

Load the GPS data.

GPSV4b <- read.table("GPS_Versionx.txt", header = T)
## The original GPS data use lon in [0, 360]

Detrend the GPS data by the GIA prior.

GPS_data <- GPSV4b
GPS_loc <- do.call(cbind, Lll2xyz(lat = GPS_data$lat, lon = GPS_data$lon))
GPS_sp <- SpatialPoints(data.frame(lon = ifelse(GPS_data$lon>359.5, GPS_data$lon - 360, GPS_data$lon), 
                                   lat = GPS_data$lat), proj4string = CRS("+proj=longlat"))
GPS_idx <- over(GPS_sp, Plist)
GPS_mu <- ice6g$trend[GPS_idx]
GPS_data$trend0 <- GPS_data$trend - GPS_mu

2.3 Prior setup for the parameters

Setup the priors for the parameters for the Gaussian process. Assume the prior distributions of both \(\rho\) and \(\sigma^2\) are log normal. The prior mean of the correlation length is set to be 500km based on expert prior opinion on the residual process and the mean for the variance is 20 which is about the range of the GIA values. The variances of both distributions are set to be large for a vague enough prior.

## Priors mean and variance for the parameters: rho and sigma
mu_r <- 500/6371
v_r <- (1000/6371)^2
mu_s <- 20
v_s <- 40^2

## Transform the parameters for the SPDE_GMRF approximation
trho <- Tlognorm(mu_r, v_r)
tsigma <- Tlognorm(mu_s, v_s)

3 Inference by INLA

We use INLA to do the estimation and prediction of the global stationary model. First set up the INLA model.

Mesh_GIA <- mesh
## The SPDE model
lsigma0 <- tsigma[1]
theta1_s <- tsigma[2]
lrho0 <- trho[1]
theta2_s <- trho[2]
lkappa0 <- log(8)/2 - lrho0
ltau0 <- 0.5*log(1/(4*pi)) - lsigma0 - lkappa0
GIA_spde <- inla.spde2.matern(Mesh_GIA, B.tau = matrix(c(ltau0, -1, 1),1,3), B.kappa = matrix(c(lkappa0, 0, -1), 1,3),
                              theta.prior.mean = c(0,0), theta.prior.prec = c(sqrt(1/theta1_s), sqrt(1/theta2_s)))

## Link the process to observations and predictions
A_data <- inla.spde.make.A(mesh = Mesh_GIA, loc = GPS_loc)
A_pred <- inla.spde.make.A(mesh = Mesh_GIA, loc = rbind(GPS_loc, Mesh_GIA$loc))

## Create the estimation and prediction stack
st.est <- inla.stack(data = list(y=GPS_data$trend0), A = list(A_data),
                     effects = list(GIA = 1:GIA_spde$n.spde), tag = "est")
st.pred <- inla.stack(data = list(y=NA), A = list(A_pred),
                      effects = list(GIA=1:GIA_spde$n.spde), tag = "pred")
stGIA <- inla.stack(st.est, st.pred)

## Fix the GPS errors
hyper <- list(prec = list(fixed = TRUE, initial = 0))
formula = y ~ -1 +  f(GIA, model = GIA_spde)
prec_scale <- c(1/GPS_data$std^2, rep(1, nrow(A_pred)))

Then we run the INLA model. Note that this will take more than 10min and require memory larger than 32GB. We ran this on a server with 56 cores and 256GB memory.

## Run INLA
res_inla <- inla(formula, data = inla.stack.data(stGIA, spde = GIA_spde), family = "gaussian",
                 scale =prec_scale, control.family = list(hyper = hyper),
                 control.predictor=list(A=inla.stack.A(stGIA), compute =TRUE))

INLA_pred <- res_inla$summary.linear.predictor

4 Analyse results

Now assemble the inla inference and prediction results.

INLA_pred <- res_inla$summary.linear.predictor

## Extract and project predictions
pred_idx <- inla.stack.index(stGIA, tag = "pred")$data
GPS_idx <- pred_idx[1:nrow(GPS_data)]
GIA_idx <- pred_idx[-(1:nrow(GPS_data))]

## GPS 
GPS_u <- INLA_pred$sd[GPS_idx]
GPS_pred <- data.frame(lon = GPS_data$lon, lat = GPS_data$lat, u = GPS_u)

## GIA
GIA_diff <- INLA_pred$mean[GIA_idx] 
GIA_m <- GIA_diff + GIA_prior
GIA_u <- INLA_pred$sd[GIA_idx]
proj <- inla.mesh.projector(Mesh_GIA, projection = "longlat", dims = c(360,180), xlim = c(0,360), ylim = c(-90, 90))
GIA_grid <- expand.grid(proj$x, proj$y)
GIA_pred <- data.frame(lon = GIA_grid[,1], lat = GIA_grid[,2],
                       diff = as.vector(inla.mesh.project(proj, as.vector(GIA_diff))),
                       mean = as.vector(inla.mesh.project(proj, as.vector(GIA_m))),
                       u = as.vector(inla.mesh.project(proj, as.vector(GIA_u))))

ress <- list(res_inla = res_inla, spde = GIA_spde, st = stGIA, 
            mesh = Mesh_GIA, GPS_pred = GPS_pred, GIA_pred = GIA_pred)

4.1 Plot the posteriors of the hyper parameters.

res_inla <- ress$res_inla
GIA_spde <- ress$spde
pars_GIA <- inla.spde2.result(res_inla, "GIA", GIA_spde, do.transf=TRUE)
theta_mean <- pars_GIA$summary.theta$mean
theta_sd <- pars_GIA$summary.theta$sd

## Find the mode of rho and sigma^2
lrho_mode <- pars_GIA$summary.log.range.nominal$mode
lrho_mean <- pars_GIA$summary.log.range.nominal$mean
lrho_sd <- pars_GIA$summary.log.range.nominal$sd
rho_mode <- exp(lrho_mean - lrho_sd^2)

lsigma_mode <- pars_GIA$summary.log.variance.nominal$mode
lsigma_mean <- pars_GIA$summary.log.variance.nominal$mean
lsigma_sd <- pars_GIA$summary.log.variance.nominal$sd
sigma_mode <- exp(lsigma_mean - lsigma_sd^2)

plot(pars_GIA$marginals.range.nominal[[1]], type = "l",
     main = bquote(bold(rho("mode") == .(round(rho_mode, 4))))) # The posterior from inla output

plot(pars_GIA$marginals.variance.nominal[[1]], type = "l", 
     main = bquote(bold({sigma^2}("mode") == .(round(sigma_mode, 4))))) # The posterior from inla output

## The estimated correlation length is about 568km
rho_mode*6371
## [1] 644.5193

4.2 Plot the predictions.

GPS_pred <- ress$GPS_pred
GIA_pred <- ress$GIA_pred

map_prior <- map_res(data = ice6g, xname = "x_center", yname = "y_center", fillvar = "trend", 
                    limits = c(-7, 22), title = "Prior GIA mean field")

## Plot the GIA predicted mean
map_GIA <- map_res(data = GIA_pred, xname = "lon", yname = "lat", fillvar = "mean", 
                  limits = c(-7, 22), title = "Predicted GIA")

## Plot the GIA difference map
map_diff <- map_res(data = GIA_pred, xname = "lon", yname = "lat", fillvar = "diff", 
                  limits = c(-8, 8), title = "GIA difference: Updated - Prior")

## Plot the GIA difference map
map_sd <- map_res(data = GIA_pred, xname = "lon", yname = "lat", fillvar = "u", 
                    colpal = colorRamps::matlab.like(12), title = "Predicted uncertainties")

 
## Display
print(map_prior)

print(map_GIA)

print(map_diff)

print(map_sd)