Experiment 2a – Update Mass – use a different GRACE data set

Since there was mass conservation problem in the GRACE super mascon data set, we now use a different GRACE data set to update mass.

Before using the BHM, we first load the previous data and new data for some intial checks.

1 Check data

Check the integrated sum and mean of all the data. For the integrated mean, we define the following weights

1. global weight: area_i / Earth surface area

2. Ocean weight: area_i/ Ocean area

3. Land weight: area_i/Land area

1.1 Super Mascons data

The area of the super mascons are provided in the original data set, however they are different from the area calculated by R using the definition of the polygons. The given areas of the polor polygons are roughly double of the calculated areas. See the plot below.

The total sum of the two versions of areas are

## Give areas sum:  510064041 
##  Calculated areas sum:  509412228

The calculated area sum seems smaller. This might be caused by error in the definition of the polygons. In the following, we calculate the sum and mean using the areas given by MS.

First we calculate the sums.

## The global sum of the super Mascons is:          -470325.3 
##  The sum of the super Mascons over the ocean is:  327148283 
##  The sum of the super Mascons over the land is:   -327618609

Next the we calculate the mean over the entire Earth surface.

## Mean over the Earth 
##  Global:  -0.0009220908 
##  Ocean:   0.6413867 
##  Land:    -0.6423088

Next the we calculate the mean over the Ocean.

## Mean over the Ocean 
##  Global:  -0.001337383 
##  Ocean:   0.930255 
##  Land:    -0.9315924

Next the we calculate the mean over the Land.

## Mean over the Ocean 
##  Global:  -0.002969448 
##  Ocean:   2.065485 
##  Land:    -2.068454

1.2 GIA data

We do the same sum and mean checks for the ice6g ewh data.

First we calculate the sums.

## The global sum of the ice6g ewh is:  -173968040 
##  The ocean sum of the ice6g ewh is:   -469235184 
##  The land sum of the ice6g ewh is:    295267143

Next the we calculate the mean over the entire Earth surface.

## Mean over the Earth 
##  Global:  -0.3410699 
##  Ocean:   -0.9199506 
##  Land:    0.5788807

Next the we calculate the mean over the Ocean.

## Mean over the Ocean 
##  Global:  -0.4797513 
##  Ocean:   -1.294009 
##  Land:    0.8142575

Next the we calculate the mean over the Land.

## Mean over the Ocean 
##  Global:  -1.17989 
##  Ocean:   -3.182457 
##  Land:    2.002567

1.2.1 Compare MS GIA data and GSFC derived GIA

We compare the GIA in ewh to the one drived from GSFC by GRACE - (GRACE-GIA) data.

1.3 New GRACE data

Also check for the new GRACE data (GSFC provided by RB).

First we calculate the sums.

## Sum for the GRACE trend 
##  The global sum of the GRACE trend is:  -3353350 
##  The ocean sum of the GRACE trend is:   232709056 
##  The land sum of the GRACE trend is:    -236062406

Next the we calculate the mean over the entire Earth surface.

## GRACE trend mean over the Earth 
##  Global:  -0.00657435  (Rory's: -0.000907213) 
##  Ocean:   0.4562336 
##  Land:    -0.4628079

Next the we calculate the mean over the Ocean.

## GRACE trend mean over the ocean 
##  Global:  -0.009247526  (Rory's:  -0.001276672) 
##  Ocean:   0.6417412          ( 0.6895327) 
##  Land:    -0.6509888         (-0.6908094)

Next the we calculate the mean over the Land.

## GRACE trend mean over the land 
##  Global:  -0.02274316 
##  Ocean:   1.578285 
##  Land:    -1.601028 ( -1.696294)

1.4 GRACE - GIA (ice6g)

The GRACE - GIA(ice6g) trend is also provided by RB from the GSFC data. We will use this data to update the mass trend later. First We compare this result with GRACE(GSFC) - ice6g(MS).

We calculate the sums for both.

## Sum for the trend RB  
##  The global sum is:  -4544743 
##  The ocean sum is:    547285571 
##  The land sum is:    -551830314

## Sum for the trend RB-MS 
##  The global sum is:   170614690 
##  The ocean sum is:    701944240 
##  The land sum is:    -531329549

Next the we calculate the mean over the entire Earth surface.

## RB trend over the Earth 
##  Global:  -0.008910114  (Rory's: 0.001688293) 
##  Ocean:   1.072971 
##  Land:    -1.081881

## 
##  RB - MS trend over the Earth 
##  Global:  0.3344956 
##  Ocean:   1.376184 
##  Land:    -1.041689

Next the we calculate the mean over the Ocean.

## RB trend over the ocean 
##  Global:  -0.01253303  (Rory's:  0.002375843) 
##  Ocean:   1.509248          ( 1.118967) 
##  Land:    -1.521781         (-1.116591)

## 
##  RB - MS trend  over the ocean 
##  Global:  0.4705037 
##  Ocean:   1.93575 
##  Land:    -1.465246

Next the we calculate the mean over the Land.

## RB trend mean over the land 
##  Global:  -0.03082345 
##  Ocean:   3.711812 
##  Land:    -3.742636 ( -2.741808)

## 
##  RB - MS trend over the land 
##  Global:  1.157147 
##  Ocean:   4.760742 
##  Land:    -3.603595

1.5 Find error estimates for the new GRACE data

The is no error estimates in the new GRACE data, so we try to derive it from the GRACE data provided by MS. Below shows the standard error of MS’ data.

We sample the error from this data to the 1 degree grid.

2 BHM Analysis

2.1 Prior setup

Use the polygon data to have rough estimate of the correlation length and variance.

2.2 Generate mesh

We also need a mesh to represent the mass process. For the same reason as discussed in previous section for SSH, we need a mesh with 1 degree resolution. The GRACE data is all over the globe, so we can use the .

2.3 Link data to process

grace_loc <- do.call(cbind, Lll2xyz(lon = graceR@coords[,1], lat = graceR@coords[,2]))
A_GRACE_data <- inla.spde.make.A(mesh = mesh0, loc = grace_loc)
A_M_pred <- inla.spde.make.A(mesh = mesh0, loc = rbind(mesh0$loc))

We use the GIA corrected GRACE trend in the data.

## Create the estimation and prediction stack
st.est <- inla.stack(data = list(y=graceR$trendMgia), A = list(A_GRACE_data),
                     effects = list(M= 1:M_spde$n.spde), tag = "est")
st.pred <- inla.stack(data = list(y=NA), A = list(rbind(A_GRACE_data, A_M_pred)),
                      effects = list(M=1:M_spde$n.spde), tag = "pred")
stM <- inla.stack(st.est, st.pred)

2.4 INLA inference

Now we can run INLA for the Bayesian inference. Do not run on a desktop the process may use up to 64GB memory at peak. We ran this on a server with enough memory.

## Fix altimetry errors as they are known
hyper <- list(prec = list(fixed = TRUE, initial = 0))
prec_scale <- c(1/graceR$std^2, rep(1, nrow(A_GRACE_data) + nrow(A_M_pred)))

## The formular for modelling the SSH mean
formula = y ~ -1 +  f(M, model = M_spde)

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

2.5 Results

2.5.1 Assemble and save results

Aftern running the INLA procedure, we assemble the inference and prediction results and save them as “exp2a_M2.RData”.

2.5.2 Plot the posteriors of the hyper parameters

## [1] "The estimated correlation lengths are:  1050.47907714161"

## [1] "The estimated marginal variances are:  426.381999313185"

2.6 Plot the predictions

We first plot the predicted \(X_{mass}\) on a 1 degree grid and then the aggregated \(X_m\) on the same super mascons grid as the GRACE data.

3 Sanity checks on the updated mass

To check whether the updated mass is reasonable, we calculate the mass trend averaged over Earth sphere, ocean, and land.

## GRACE predicted trend mean over the Earth 
##  Global:  0.08466148 
##  Ocean:   1.106984 
##  Land:    -1.022323

## GRACE predicted trend mean over the ocean 
##  Global:  0.1190854 
##  Ocean:   1.557092 
##  Land:    -1.438006

## GRACE predicted trend mean over the land 
##  Global:  0.2928761 
##  Ocean:   3.829478 
##  Land:    -3.536602