Modeling
Numerical models
Numerical models are widely used to solve physical or chemical problems by describing processes using equations and numbers.
We can generally describe a numerical model as y=F(x), where x are the causes of the process, F are the algorithms that describe the process, and y are the consequences of the process.
In Earth Sciences, we typically know the consequences (y) of the process we are studying (e.g. a concentration of something as a result of time, that will be the main cause in a geochronology problem), and we want to know the causes x. Therefore, finding or approximating the inverse model x=F’(y) would be very useful for our purposes.
The way we solve a problem involving numerical models will depend on whether or not we can find or approximate the inverse model.
Forward problem
Forward problem
Forward models as y=F(x) are used to make informed predictions. However, when we can find the inverse of our model $x=F’(y)$, we can solve our problem directly.
In geochronology, this x=F’(y) typically means expressing the true age of a sample as a function of concentrations and other known parameters.
Mathematically, a forward problem is a wellposed problem, where a unique solution exists and the solutions change continuously with the known parameters.
Radiocarbon calibration
A good example of a forward problem is the calculation of a calibrated C14 age from an apparent C14 age, as the C14 ages calculated in the previous session.
Apparent radiocarbon ages are calibrated using calibration lines. You
can download the data corresponding to a calibration line from here
You should get a file called 16947259732SP.14c
. Save it in your
working folder and open it as text: you will find out that is a comma
separated file (csv).
If we have a radiocarbon age (e.g. 2000±20), we should be able to
create a script to calibrate the entire probability distribution of the
age (see previous session) along the calibration curve using interp1
and produce an output like this:
%% import calibration curve
fid = fopen('16947259732SP.14c');
imported = textscan(fid, '%f %f %f %f %f',...
'HeaderLines', 12,'Delimiter',',');
fclose(fid);
% select the data from the calibnration curve
calBP=imported{1};
C14=imported{2};
%% interpolate calibration curve to arrays of 1 position per calibrated year
refcalBP=min(calBP):1:max(calBP);
refC14=interp1(calBP,C14,refcalBP);
%% input our data
C14age=2000;
C14ageerr=20; % one sigma error
% calculate probabilities of my data
C14probs=normalprobs(refC14,C14age,C14ageerr);
%% Plot the calibration curve and our data
% select only the "most" probable data to plot
sel=(C14probs>max(C14probs)/1000);
figure % start figure
% plot the part of the calibration curve that is relevent for us
subplot(3,3,[2 6])
plot(refcalBP(sel),refC14(sel))
title('Calibration curve')
% plot our data
subplot(3,3,[1 4])
plot(C14probs(sel),refC14(sel))
ylabel('C14 age')
xlabel('P')
% plot the data calibrated
subplot(3,3,[8 9])
plot(refcalBP(sel),C14probs(sel))
xlabel('calibrated C14 age')
ylabel('P')
Note that:

We are ignoring the uncertainty of the calibration curve (
imported{3}
). To know how to incorporate all the errors, check the script “MatCal” by Lougheed & Obrochta (2016) 
To represent several subplots in the same window we are using
subplot(r,c,[a b])
, wherer
andc
are the number of rows and columns, anda
andb
are corners of the area where we want to plot. E.g.subplot(3,4,[7 12])
would start plotting in the blue area:
Inverse problem
Inverse problem
Sometimes it is not possible to express our problem as x=F’(y), because sometimes it is impossible to get the inverse of F(x).
Most geochemical models used in geochronology allow the calculation of a theoretical concentration or measurable signal C as a function of time t and other parameters: C=f(t,C_0,…). However, some of these problem cannot be solved for t=f(C,C_0,…). In this cases, we will need to guess the age $t$ corresponding to our known concentrations C,C_0, etc.
Illposed problem
Mathematically, this kind of problems are often illposed problems. Therefore, we cannot assume that they have a unique solution and we should check the sensitivity of our results to a change in our known parameters.
In geochronology, this means that apart of answering the main question:
Which ages are compatible with my data?
but we should also answer the question:
Which ages are not compatible with my data?
Cosmogenic depthprofile dating
Cosmogenic depthprofile dating
The accumulation of Be10 under a sedimentary surface depends on the inherited Be10 concentration (C0), the different Be10 production rates (Psp, Pfμ and Psμ) and attenuation lengths (Λsp, Λfμ andΛsμ), the Be10 decay constant (λ), the density of the sediment (ρ), the depth (z), the erosion rate of the surface (ε) and the age of the landform (t):
C=C0+
Psp./(λ+ε.*ρ./Λsp).*exp(z.*ρ./Λsp).*(1exp((λ+ε.*ρ./Λsp).*t))+
Pfμ./(λ+ε.*ρ./Λfμ).*exp(z.*ρ./Λfμ).*(1exp((λ+ε.*ρ./Λfμ).*t))+
Psμ./(λ+ε.*ρ./Λsμ).*exp(z.*ρ./Λsμ).*(1exp((λ+ε.*ρ./Λsμ).*t))
This equation cannot be solved for t. Also, when we have a dataset of Be10 concentrations under a surface (a Be10 depthprofile), we want to solve the problem for C0, ε andt. How can we do this?
Be10 accumulation model
The following function calculates theoretical Be10 concentrations:
function [ C ] = exposure_model(P,L,l,density,z,C0,erosion,t)
C=C0+...
P(1)./(l+erosion.*density./L(1)).*exp(z.*density./L(1)).*...
(1exp((l+erosion.*density./L(1)).*t))+...
P(2)./(l+erosion.*density./L(2)).*exp(z.*density./L(2)).*...
(1exp((l+erosion.*density./L(2)).*t))+...
P(3)./(l+erosion.*density./L(3)).*exp(z.*density./L(3)).*...
(1exp((l+erosion.*density./L(3)).*t));
end
Be10 data
The following code defines all the known parameters and the Be10 concentrations from the sampled depthprofile for an alluvial fan in Almería (Spain):
%% Production rates
P=[4.35,0.0985,0.0855]; % production rates in at/g/a
L=[160,1137,1842]; % attenauation lengths in g/cm^2
l=4.9975E7; % decay contant in a^(1)
%% Field data
density=1.8; % g/cm^3
z=[267,195,141,95,46,3]; % depth of the sameples in cm
Be10=[91000,184000,265000,430000,732000,1070000]; % Be10 concentrations in atoms/g
Be10error=[9100,16000,18000,29000,61000,81000]; % Be10 uncertainties in atoms/g
Try exposure_model(P,L,l,density,10,0,0.0001,10000)
to calculate the
Be10 concentration accumulated in a sample 10 cm below a 10 ka old
surface being eroded at a rate of 1 mm/ka (0.0001 cm/a).
We can reproduce the theoretical depth profile for these conditions along the first 3 m under the surface:
zref=0:300; % depth reference in cm
concentrations=exposure_model(P,L,l,density,zref,0,0.0001,10000);
plot(concentrations,zref,'b')
Now we just need to find which theoretical values of inheritance, age
and erosion rates (the last 3 parameters in the exposure_model
function) match our Be10
concentrations within Be10error
uncertainties!
Monte Carlo methods
Monte Carlo methods
The simplest way of guessing the values for C0,erosion,t
that fit our
data Be10
at our depths z
could be just trying a lot of random
values of C0,erosion,t
and check which theoretical concentrations are
closer to our data. This is called a Monte Carlo experiment.
To perform this Monte Carlo experiment, we should define a way of measuring how close is our model to our data. A Χ² function would do the job:
The following code runs a MonteCarlo experiment of 100 000 models
assuming that $\epsilon$ is between 0 and 50 m/Ma (0.005 cm/a), the
landfom age is between less than 3 Ma (3E6
a), and $C_0$ is smaller
than the lowest concentration.
%% Monte carlo experiment
nummodels=100000; % define how many models
C0i=rand(1,nummodels)*min(Be10); % random inheritences
ti=rand(1,nummodels)*3e6; % random ages
erosioni=rand(1,nummodels)*0.005; % random erosion rates
chisquarevalues=rand(1,nummodels)*NaN; % allocate memory for the chi square array
% calculate the chi squared vales for each model
for n=1:nummodels
% calculate the model concetratios for the depths z
Cmodel=exposure_model(P,L,l,density,z,C0i(n),erosioni(n),ti(n));
% calculate the chi squared for this model
chisquarevalues(n)=sum(((CmodelBe10)./Be10error).^2);
end
Which models should we consider to represent the uncertainty of the results?
When fitting a model to data, we have to report how many parameters are we trying to fit and how many data we have. The number of parameters should be lower that the data points and the difference between them are the Degrees of Freedom of our model. In our model we have
DOF = 6  3 = 3 degrees of freedom.
When performing this kind of inverse modeling, the models that fit the data with a Χ² value below the minimum Χ² value plus the degrees of freedom are often considered to fit the data within one sigma confidence level.
Therefore, we can calculate which of the models fit our data within onesigma, assuming that this is defined by the models with Χ² values between the minimum Χ² and Χ²+DOF:
DOF=3; % degrees of freedom (# of samples  # of parameters)
minchi=min(chisquarevalues); % minimum chisquare value (best model)
best=find(chisquarevalues==minchi); % location of the best model
% location of the models fitting the data within onesigma
onesigma=find(chisquarevalues<minchi+DOF);
% display previous infrormation
disp(['Min chisquared value = ' num2str(minchi)])
disp([num2str(length(onesigma)) ' modles fitting one sigma'])
disp(['Age: ' num2str(min(ti(onesigma))/1e3) '  '...
num2str(max(ti(onesigma))/1e3) ' ka'])
disp(['Erosion: ' num2str(min(erosioni(onesigma))*1e4) '  '...
num2str(max(erosioni(onesigma))*1e4) ' m/Ma'])
disp(['Inheritance: ' num2str(min(C0i(onesigma))) '  '...
num2str(max(C0i(onesigma))) ' atoms/g'])

What is the best result?

Does it fit the data well? (Χ²~0)

Plot randomized values against their corresponding Χ² values to get a idea of the distribution of the results. Tip: plot only the Χ² values below the best value+10.

Plot the sample Be10 concentrations and the theoretical Be10 of the best fit.

Select the models fitting the data within onesigma confidence level and plot these models in grey (
’Color’,[0.7 0.7 0.7]
).
We should get a high number of fitting models to get an idea of which the distribution of the parameters values that fit our data. Try increasing the number of random models to get at least 300 fitting models.
Convergence methods
Convergence
Another way of getting more models fitting our data is changing the limits of the randomized parameters while generating more models. For example, until now we have been considering ages that are between 0 and 3 Ma, but after running about 100000 models, it is pretty clear that the ages fitting the data are lower than 1 Ma, and that the erosion rates should be lower than 5 m/Ma (0.0005 cm/a). Rerun your solver applying better limits.
We can even program our solver to start converging after a learning process of a certain number of models, automatizing what we have just done manually. However, we should be cautious making our random models to converge very fast because we can miss solutions that fit our data. To avoid that, we could make them converge to Χ² values Χ²min+10·DOF. Actually, we should also allow our random models to diverge out of the initial limits when we find good solutions close to our limits.
Goalseeking algorithms
Goalseeking algorithms
Until now, we have been solving the question “Which ages are compatible with my data?”, but are not explicitly answering the question “Which ages are not compatible with my data?”.
The cosmogenic depthprofile models often show that the fitting models are scattered towards old ages. This is because the fitting area in the $\epsilont$ space is a narrow valley that we can easily miss when randomizing the $\epsilon$ and $t$ values.
To avoid this, we can randomize only the $C_0$ and $t$ parameters and make our program to seek actively which is the best $\epsilon$ that fit the data for each of the models. This operation is sometimes called “Χ² minimization”. Χ² minimization slows down our program but will guarantee that the models outside the fitting age range do not fit our data.
To minimize the value of Χ² we can use some builtin functions as
fminunc
or fminsearch
(type help fminsearch
for more information).
However, the way these algorithms work change in the different versions
of MATLAB and Octave, so we will never be sure that our program is going
to work the same way in someone else computer. Therefore, it is highly
recommended to build our own minimization algorithm.
An easy solution could be to use interp1
as a goal seeker of the
deviations. The piece of code in the next slide includes this goal
seeker in the modeling loop to force getting always the best erosion
rate.
How many models do you need to run now to get 300 fitting resutls?
% start testing models
for n=1:nummodels
erosionref=[0,logspace(5,2,100),10^10]'; % define an array with erosion rates
% calculate the deviations corresponding to each erosion rate
% deviations are defined as the sum of (CmodelCsample)/Uncertainty
% for all the samples
deviations=...
sum(...
(exposure_model(P,L,l,density,z,C0i(n),erosionref,ti(n))Be10)./...
Be10error...
,2);
% Interpolate the erosion rate values to find the one the model that
% fit the data better (for the age and Co corresponding to this random
% model). THen store the result at erosioni(n), ooverwriting the
% previously defined value.
erosioni(n)=interp1(deviations,erosionref,0);
% calculate the model concetratios for the new erosion rate at the depths z
Cmodel=exposure_model(P,L,l,density,z,C0i(n),erosioni(n),ti(n));
% calculate the chi squared for this model
chisquarevalues(n)=sum(((CmodelBe10)./Be10error).^2);
end
Summary
Summary

Using convergence (and divergence) algorithms help our inverse modeling program to run faster and allow us to set wide starting limits.

Using goalseeking algorithms slows down the calculation of each individual simulation. However, it usually allows us to run less models, and also guarantee the reproducibility and accuracy of our results. Compare the two plots representing the solutions in the $\epsilont$ space with and without using Χ² minimization:
vs.
Exercise
Use previous models to solve the age of a landform with the following Be10 depth profile:
Sample depth Be10
cm atoms/g
 
250 25±2
163 45±3
113 60±5
73 100±7
43 140±10
11 200±15
 
Start testing ages between 0 and 10 Ma and try introducing some convergence code. This will allow you to create a code that will work on any Be10 database.