TCNtools_eulerian_solver.Rmd 6.56 KB
Newer Older
GODARD Vincent's avatar
GODARD Vincent committed
1
---
GODARD Vincent's avatar
dev app    
GODARD Vincent committed
2
title: 'Introducing the Eulerian solver for concentration calculations'
GODARD Vincent's avatar
GODARD Vincent committed
3
4
5
output:
  html_document:
    df_print: paged
GODARD Vincent's avatar
GODARD Vincent committed
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
---

# Objectives of this Notebook

## Background

We are going to consider simple computations of concentration under various conditions in terms of erosion, depth or age.
This will be done using an Eulerian point of view, which is the most straightforward and fastest way to perform such computation. 
In this case the quantity of interest (concentration) is computed at fixed depths below the surface, while the exhumed material is moving through this reference frame during its trajectory toward the surface.
More details on the differences between Eulerian and Lagrangian approaches, and their application to complex exposition/denudation histories, will be provided in another tutorial.

The relevant equation is the following,

$$
C=C_0e^{-\lambda t} + \sum_i \frac{P_i}{\frac{\rho \varepsilon}{\Lambda_i}+\lambda}e^{\frac{-\rho z}{\Lambda_i}}(1-e^{-(\frac{\rho \varepsilon}{\Lambda_i}+\lambda)t})
$$
with the following variables and parameters,

- $C$ the concentration (as a function of time $t$ and depth $z$)
- $C_0$ the inherited concentration
- $\lambda$ the decay constant for the considered nuclide
- $P_i$ the scaled surface production rate for the nuclide of interest and the $i$-th production pathway (spallation, stopped muons, fast muons)
- $\rho$ the density of the medium
- $\Lambda_i$ the attenuation length for the particules of the $i$-th production pathway
- $\varepsilon$ surface denudation

In order to stick with usual conventions in the following time will be measured in years (a), the unit of length will be cm and the depths will be expressed in g/cm$^2$ (i.e. actual depth $\times \rho$).

## Set up

The first thing we have to do is to load the **TCNtools** library (once it has been installed)

```{r ck_1}
library("TCNtools")
```

We should the define the basic parameters we are going to use for the computation, which are two vectors :

- a vector with the attenuation lengths for different particules (in g/cm$^2$)
    - neutrons for spallation reactions $\Lambda_{spal}$
    - stopping muons $\Lambda_{stop}$
    - fast muons $\Lambda_{fast}$    
- a vector (or matrix) with the SLHL production rates (in at/g/a), in this case for the *st* scaling scheme (@stone2000air), and decay constant $\lambda$ (in 1/a) for the nuclide(s) of interest.

```{r ck_2}
# Attenuation lengths
Lambda = c(160,1500,4320) # g/cm2
names(Lambda) <- c("Lspal","Lstop","Lfast") # we just give names to the element of the vector
# Production and decay parameters
prm = matrix(c( 4.01 , 0.012 , 0.039 , log(2)/1.36e6,
               27.93 , 0.84  , 0.081 , log(2)/0.717e6),
                nrow = 4,ncol=2 ) 
colnames(prm) <- c("Be10","Al26") # we just give names to the columns of the matrix               
rownames(prm) <- c("Pspal","Pstop","Pfast","lambda")  # we just give names to the rows of the matrix  
# material density
rho = 2.7
```


We also need to define the properties of our site of interest and compute the relevant scaling parameters.
```{r ck_3}
altitude = 1000 # elevation in m
latitude = 45 # latitude in degrees
P = atm_pressure(alt=altitude,model="stone2000") # compute atmospheric pressure at site
st = scaling_st(P,latitude) # compute the scaling parameters according to Stone (2000)
```

# Concentration along a profile

We first compute the changes in concentration with depth $z$ along a profile.
We are going to use the `solv_conc_eul` function. 
As always the notice of the function, including its various arguments, can be obtained  by typing `?solv_conc_eul` in the R console.
We consider no inheritance ($C_0$=0), so the evolution is starting from a profile with homogeneous zero concentration.

```{r ck_4}
z = seq(0,500,by=10) * rho # a vector containing depth from 0 to 500 cm by 10 cm increments,  and then converted into g/cm2 
C0 = 0 # inherited concentration
age = 10000 # the time in a
ero = 10 * 100/1e6*rho # denudation rate expressed in m/Ma and converted in g/cm2/a
C = solv_conc_eul(z,ero,age,C0,prm[,"Be10"],st,Lambda) # compute concentration
plot(C,z,type="l",ylim=rev(range(z)),lwd=3,xlab="Concentration (at/g)",ylab="Depth (g/cm2)")
```

Try to modify the `age` and `ero` (always keeping it in g/cm$^2$/a) parameters above, to see their influence on the profile. 

# Evolution of concentration with time

Now we are going to consider the evolution of concentration with time $t$.
The computation will be carried out at the surface ($z=0$), but this could be done at any arbitrary depth. 

```{r ck_5}
age =  seq(0,100e3,by=100) # a vector containing time from 0 to 100 ka by 100 a steps
z = 0 * rho # depth at which we are going to perform the calculation (cm converted to g/cm2)
C0 = 0 # inherited concentration
ero = 10 * 100/1e6*rho # denudation rate expressed in m/Ma and converted in g/cm2/a
C = solv_conc_eul(z,ero,age,C0,prm[,"Be10"],st,Lambda) # compute concentration
plot(age/1000,C,type="l",lwd=3,ylab="Concentration (at/g)",xlab="Time (ka)")
```
We see here the progressive build-up of concentration though time and the establishment of balance between gains (production) and losses (denudation and decay) leading to the concentration plateau at steady state.

Try to modify the  `ero` (always keeping it in g/cm$^2$/a) parameter above, to see its influence on time needed to reach steady state and the final concentration. 

# Evolution of concentration with denudation rate

Now we are going to consider the evolution of concentration with denudation rate $\varepsilon$.
The computation will be carried out at the surface ($z=0$), but this could be done at any arbitrary depth. 
We will consider that $t=+\infty$ and that we have reached the plateau concentration.

```{r}
ero = 10^seq(log10(0.1),log10(1000),length.out = 100) * 100/1e6*rho # a log-spaced vector for denudation rate expressed in m/Ma and converted in g/cm2/a
age =  Inf # infinite age
z = 0 * rho # depth at which we are going to perform the calculation (cm converted to g/cm2)
C0 = 0 # inherited concentration
C = solv_conc_eul(z,ero,age,C0,prm[,"Be10"],st,Lambda) # compute concentration
plot(ero/100*1e6/rho,C,log="xy",type="l",lwd=3,ylab="Concentration (at/g)",xlab="Denudation rate (m/Ma)")
```

This figure highlights the strong inverse relationship, at steady-state, between  denudation rate ($\varepsilon$) and concentration ($C$), which is the foundation of many geomorphological studies trying to establish landscape evolution rates.
Note the change in the relationship at very low denudation rates, which corresponds to the situation where the effects of radioactive decay become predominant.


# Dealing with inheritance



# References