TCNtools_eulerian_solver.Rmd 6.56 KB
 GODARD Vincent committed Nov 24, 2020 1 ---  GODARD Vincent committed Dec 03, 2020 2 title: 'Introducing the Eulerian solver for concentration calculations'  GODARD Vincent committed Nov 27, 2020 3 4 5 output: html_document: df_print: paged  GODARD Vincent committed Nov 24, 2020 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