Markov Chain Monte Carlo Algorithm for Taylor Glacier 14C data

This readme file presents basic descriptions for the Markov Chain Monte Carlo (MCMC) algorithm written in MATLAB. The MCMC method is used to constrain two parameters in cosmogenic 14C production model by muons. For simplicity, the two model parameters that are optimized with the MCMC method are "fneg" and "ffast" – corresponding scaling factors for negative muon capture and fast muon reaction that are constant with respect to depths. The MCMC method aims to optimize these 2 parameters given the observations – which in this case represent 14 total 14C measurements from 7 unique depth levels (6.85m, 15m, 19.5m, 40m, 51m, 61.5m, 72m) obtained from ice cores drilled in Taylor Glacier, Antarctica.

File descriptions

External data
a. flowpath_MC.mat – contains the pool of 1000 flowpaths (see section 1)
b. flowpath_trim.mat – contains the flowpaths that didn't crash into bedrock (see section 2)
c. P_neutron.mat – production rate from neutron as function of depth (see section 3)
d. P_muon.mat – production rate from muons as function of depth (see section 3)
e. all_data.mat – contains measurement data, sample depths and uncertainties
Functions
a. trim_trajectory.m – this function removes part of the flow trajectory that goes out of bounds (goes beyond 72km), see section 2
b. calc14C_v3.m – the main integration function for 14C concentration in the ice parcel (see section 3)
c. interp1qr_m – function for faster 1D interpolation
Run scripts
a. trim_bedrock.m – trim flowpath_MC.mat into flowpath_trim.mat (see section 2)
b. run_MCMC.m – run the MCMC algorithm (see section 4)

1. Ice flow model

Buizert et al. (2012) developed a 2D ice flow model for Taylor Glacier. Using estimates of ablation rates along the glacier and bedrock profile from Kavanaugh et al. (2009a) the model generates an ice parcel back-trajectory for each given sample depth. For the MCMC algorithm I perturbed the ablation rates within their uncertainties and generated 1000 flow trajectories for 8 unique depth levels (6.85m, 10m, 15m, 19.5m, 40m, 51m, 61.5m, 72m). I added the 10m sample depth because we are only missing 14CH4 data on that depth level, and 14CH4 account for less than 1% of total 14C signal. At the moment, the 10m data are not included in the MCMC fit algorithm, but later they can be if necessary. The 8x1000 flow trajectories are saved in " flowpath_MC.mat" file.

Important variables in the " flowpath_MC.mat" file include

age (1x861) – this is the age of the ice parcel trajectory, going from 0 at the drill location to 6000 years before present.
h_age_save (8x861x1000) – this is the depth of the ice parcel trajectory (y-location in the glacier). The first dimension (8) represents the 8 unique depth levels. The 2nd dimension (861) corresponds to the age and the 2nd dimension of x_age_save (x-location in the glacier). The 3rd dimension (1000) represents the pool of 1000 possible ice parcel trajectories given the perturbation in ablation rates.
x_age_save (8x861x1000) – this is the x-location of the ice parcel trajectory. The first dimension (8) represents the 8 unique depth levels. The 2nd dimension (861) corresponds to the age and the 2nd dimension of h_age_save (y-location in the glacier). The 3rd dimension (1000) represents the pool of 1000 possible ice parcel trajectories given the perturbation in ablation rates.
flow_rand (1x1000) – this is the random number (η) that corresponds to the likelihood of each 1000 possible ice parcel trajectories. The probability for a given flow_rand is
P(η|I) = exp(-0.5*((η- ηbest)/σ)^2)
where ηbest = mean of η and σ = stdev of η. because the flow_rand variable is generated with MATLAB's randn function, ηbest = mean = 0 and σ = stdev = 1, thus P(η|I) simplifies into
P(η|I) = exp(-0.5*(flow_rand)^2)

2. Trimming the trajectories.

Under extreme ablation rate scenarios, the ice parcel trajectories crash into bedrock and become unphysical afterwards. I use the "trim_bedrock.m" run script to simply remove the trajectories that crash into bedrock and resave the remaining trajectories in "flowpath_trim.mat" file.

Additionally, the bedrock data only extend to 72km upstream of the glacier terminus. Given different ablation rates scenario, some ice parcel trajectories go "out of bounds" beyond 72km faster than others and the flow trajectories also become unphysical afterwards since there's no bedrock.

To account for the limitation in the extent of the bedrock data, we assume that the amount of 14C at the beginning of the trajectory (corresponding to 72km in the x-direction) is at steady state (production equals decay). The age of our 72m deep sample is approximately 70 kyr BP. Survey by Morse et al. (1998) showed that the depth of ~70 kyr BP ice at northern flank of Taylor Dome is ~575m, thus we assumed that the depth of the long term transport for the 72m sample is 575m. At 72km, the 72m sample under the "best" ablation scenario is at 699m deep, thus we use 699m as the baseline depth "z_baseline." For other ice parcel trajectories, the steady state depth of long-term transport (z_deep) is calculated following

z_deep = 575 – (z_baseline – zx=72km)

where zx=72km is the depth of the corresponding sample at 72km x-direction.

We want to trim the unphysical parts out so that the integration of 14C content in the model can be consistent. MATLAB doesn't like having a matrix with different number of elements, thus the function "trim_trajectory.m" takes an individual ice parcel trajectory [h_age (1x861), x_age(1x861), age(1x861),z_baseline(1x1)] and transform them into cell {} format where it is possible to save individual vectors with different elements in them.

3. 14C production model

The baseline 14C production rate by neutrons, in unit of 14C atoms g ice-1 yr-1 as function of depth (m) is obtained from LSD model by Lifton et al. (2014) and saved in "P_neutron.mat" file. The baseline 14C production rate by muons, in unit of 14C atoms g ice-1 yr-1 as function of depth (m) is obtained from LSD model by Balco et al. (2008) and saved in "P_muon.mat" file. Both baseline 14C production rate files are scaled to the altitude and latitude of Taylor Glacier drill site.

From the baseline production rates, the amount of 14C in the sample for a given flow trajectory is calculated by the differential equation dC/dt = fneg*Pneg(z) + ffast*Pfast(z) – λC fneg = constant scaling factor for negative muon capture (unitless) Pneg = 14C production from negative muon capture (14C atoms g ice-1 yr-1) fneg = constant scaling factor for fast muon reaction (unitless) Pneg = 14C production from fast muon reaction (14C atoms g ice-1 yr-1) z = sample depth (m) – note that z is also a function of time z(t) λ = radiocarbon decay constant (1.216e-4 yr-1) C = 14C concentration in the ice parcel (14C atoms g ice-1)

As an initial condition (C0) before the integration, we assumed that at the depth of long-term transport (z_deep), the 14C concentration in the ice parcel is at equilibrium dC/dt (at z_deep) = 0 = fneg*P_neg*z_deep + ffast* P_fast*z_deep - C0*λ currently we are ignoring production from neutron because we're only fitting to data below 6.85m.

In the MATLAB code, the expected total 14C in the ice parcel at the end of flow trajectory (Cexp) is calculated with the function "calc14C_v3.m". The "calc14C_v3.m" function takes P_n (production rate from neutron), P_neg (production rate from negative muon), P_muf (production rate from fast muon), z_P (depth corresponding to the production rate), h_age (depth vector of the ice trajectory), age (age vector of the age trajectory), C0 (initial condition), and f (3x1 matrix that linearly scale the production rates, [fn fneg ffast]).

4. Markov Chain Monte Carlo algorithm**

The MCMC is run with the "run_MCMC.m" script. The general recipe for the MCMC is