You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
381 lines
14 KiB
381 lines
14 KiB
#============================================================================= |
|
# Remineralization with reactivity distribution |
|
# |
|
# Script corresponds to the second application in the manuscript. |
|
# It produces figure 2 in the manuscript and some additional plots. |
|
# |
|
# Coded by J. Rooze |
|
#============================================================================= |
|
|
|
library(marelac) |
|
library(ReacTran) |
|
library(rootSolve) |
|
|
|
#============================================================================= |
|
# Units used in the program: |
|
#============================================================================= |
|
|
|
#normalize |
|
time_scale <- 1 / 1e-3 #time scale = 1/k [yr] |
|
length_scale <- 15 * 1e-2 #[m] |
|
|
|
|
|
#============================================================================= |
|
# Model domain and grid definition |
|
#============================================================================= |
|
L <- 50 * 1e-2 / length_scale # depth of sediment domain [-] |
|
N <- 50 # number of grid layers |
|
grid <- setup.grid.1D(x.up = 0, L = L, N = N) |
|
|
|
# A handy function for plotting profiles |
|
plot_profile <- function(...) { |
|
plot(y = grid$x.mid * length_scale, ylim = c(L,0) * length_scale, type = "l", |
|
ylab = "depth (m)", ...) |
|
} |
|
|
|
#============================================================================= |
|
# Model parameters: |
|
#============================================================================= |
|
|
|
PerSecToPerYr <- 365*24*60^2 |
|
dens_sed <- 2.6 * 1e6 #g/m3 |
|
|
|
SAR <- 0.001 * (time_scale / length_scale) #[m/yr] [yr/m] = [-] |
|
|
|
parms <- c(Db_max = 1e-10 * PerSecToPerYr * time_scale / length_scale^2, #[m2/yr] [yr/m^2] = [-] |
|
Db_efold = 0.02 / length_scale, #[m]/ [m] = [-] |
|
TOC.swi = 3/100) #[g C / g sed] |
|
|
|
Db <- setup.prop.1D(func = p.exp, grid = grid, y.0 = parms["Db_max"], |
|
y.inf = 0, x.att = parms["Db_efold"]) |
|
|
|
#============================================================================= |
|
# Distribution functions |
|
#============================================================================= |
|
initial_uniform_dist <- function(k, k_min, k_max, val = 1) { |
|
val * (k >= k_min) * (k <= k_max) |
|
} |
|
|
|
#The distribution function |
|
#alternative expressions also work well, for instance: exp(b0 + b1 * k + b2*k^2) |
|
polynom_dist <- function(k, b0, b1, b2) { |
|
C0 * exp(b0 * sqrt(k) + b1 * k) + b2 * k * exp(k) |
|
} |
|
|
|
#mom = integer for moment |
|
polynom_highermom <- function(k, k_mean, mom, b0, b1, b2) { |
|
if (mom == 1) { |
|
return(k * polynom_dist(k, b0, b1, b2)) |
|
} else { |
|
return((k-k_mean)^mom * polynom_dist(k, b0, b1, b2)) |
|
} |
|
} |
|
|
|
multiroot_func <- function(x, k_min, k_max, conc, k_mean, k_var) { |
|
b0 <- x[1] |
|
b1 <- x[2] |
|
b2 <- x[3] |
|
|
|
F1 <- conc - integrate(polynom_dist, b0 = b0, b1 = b1, b2 = b2, |
|
lower = k_min, upper = k_max)$value |
|
F2 <- k_mean - integrate(polynom_highermom, b0 = b0, b1 = b1, b2 = b2, |
|
mom = 1, k_mean = k_mean, |
|
lower = k_min, upper = k_max)$value / conc |
|
F3 <- k_var - integrate(polynom_highermom, b0 = b0, b1 = b1, b2 = b2, |
|
mom = 2, k_mean = k_mean, |
|
lower = k_min, upper = k_max)$value / conc |
|
|
|
return(c(F1, F2, F3)) |
|
} |
|
|
|
polynom_dist_coefs <- function(conc, k_mean, k_var, k_min, k_max, guess, ...) { |
|
sol <- multiroot(multiroot_func, k_min = k_min, k_max = k_max, |
|
conc = conc, k_mean = k_mean, k_var = k_var, |
|
start = guess, ...) |
|
} |
|
|
|
#function to numerically integrate distribution function |
|
int_kpowC <- function(k_min = 0, k_max = Inf, b0, b1, b2, pow_ord) { |
|
integrate(f = function(k, ...) { |
|
k^pow_ord * polynom_dist(k, ...) |
|
} , lower = k_min, upper = k_max, b0 = b0, b1 = b1, b2 = b2)$value |
|
} |
|
|
|
|
|
#============================================================================= |
|
# Choose distribution (with values for working examples) |
|
#============================================================================= |
|
|
|
k_min <- 0 * time_scale |
|
k_max <- 5e-2 * time_scale |
|
conc.swi <- integrate(initial_uniform_dist, |
|
k_min = k_min, k_max = k_max, |
|
lower = k_min, upper = k_max)$value |
|
#correct for concentration in parms |
|
uniform_conc <- parms[["TOC.swi"]]/conc.swi |
|
conc.swi <- integrate(initial_uniform_dist, |
|
k_min = k_min, k_max = k_max, |
|
lower = k_min, upper = k_max, |
|
val = uniform_conc)$value |
|
|
|
#check: k_min + (k_max - k_min) * 0.5 |
|
kmean.swi <- integrate(function(k, k_min, k_max) |
|
k * initial_uniform_dist(k, k_min, k_max, val = uniform_conc), |
|
k_min = k_min, k_max = k_max, |
|
lower = k_min, upper = k_max)$value / conc.swi |
|
|
|
#check: sd(seq(k_min, k_max, length = 1e6))^2 |
|
kvar.swi <- integrate(function(k, k_min, k_max, k_mean) { |
|
(k-k_mean)^2 * initial_uniform_dist(k, k_min, k_max, val = uniform_conc) |
|
}, k_min = k_min, k_max = k_max, k_mean = kmean.swi, |
|
lower = k_min, upper = k_max)$value / conc.swi |
|
|
|
guess.swi <- c(0,0,0) #guess for an uniform distribution |
|
C0 <- uniform_conc #needed for polynom_dist function |
|
|
|
#============================================================================= |
|
# Check initial distribution |
|
#============================================================================= |
|
|
|
#calculate coefficients of polynomial distriburion |
|
coefs_init <- polynom_dist_coefs(conc = conc.swi, k_mean = kmean.swi, k_var = kvar.swi, |
|
k_min = k_min, k_max = k_max, guess = guess.swi, |
|
rtol = 1e-8 )$root |
|
|
|
kaxs <- seq(k_min, k_max, length = 100) |
|
recon_dist <- polynom_dist(k = kaxs, b0 = coefs_init[1], b1 = coefs_init[2], b2 = coefs_init[3]) |
|
|
|
#plot initial distribution |
|
plot(y = recon_dist, x = kaxs/time_scale, |
|
xlab = "k [1/y]", ylab = "dens", type = "l", ylim = c(0, max(recon_dist))) |
|
abline(v = kmean.swi/time_scale, col = 3, lty = 2) #mean |
|
|
|
#============================================================================= |
|
# Model formulation |
|
#============================================================================= |
|
|
|
model <- function (t, state, parms) { |
|
with(as.list(parms), { |
|
# Initialization of state variables |
|
TOC <- state[1:N] # #g C / g sed |
|
Ck <- state[(N+1):(2*N)] |
|
Cvar <- state[(2*N+1):(3*N)] |
|
Cage <- state[(3*N+1):(4*N)] |
|
|
|
k_mean <- Ck / TOC |
|
k_var <- Cvar / TOC |
|
|
|
# Transport and rean terms |
|
tranTOC <- tran.1D(C = TOC, C.up = TOC.swi, D = Db$int, |
|
v = SAR, VF = 1, dx = grid) |
|
tranCk <- tran.1D(C = Ck, C.up = TOC.swi * kmean.swi, D = Db$int, |
|
v = SAR, VF = 1, dx = grid) |
|
tranCvar <- tran.1D(C = Cvar, C.up = TOC.swi * kvar.swi, D = Db$int, |
|
v = SAR, VF = 1, dx = grid) |
|
tranCage <- tran.1D(C = Cage, C.up = 0, D = Db$int, |
|
v = SAR, VF = 1, dx = grid) |
|
|
|
#extra transport term for the variance |
|
kmean_bnd <- c(kmean.swi, |
|
0.5 * (k_mean[1:(N-1)] + k_mean[2:N]), |
|
k_mean[N]) #N+1 vector length |
|
dkmean_dx <- (kmean_bnd[2:(N+1)] - kmean_bnd[1:N]) / grid$dx #N vector length |
|
Cvar_mutran <- 2 * Db$mid * TOC * dkmean_dx^2 |
|
|
|
# reactions [mol/m3/yr] |
|
R <- - k_mean * TOC |
|
|
|
if(any(range(k_mean) < 0 )) |
|
browser() |
|
|
|
#Following line goes over all cells and integrates k^2 C(k) for the chosen distribution |
|
guess <- guess.swi |
|
int_kpow2C <- int_kpow3C <- rep(0, length.out = N) |
|
for (i in 1:N) { |
|
polynom_coefs <- polynom_dist_coefs(conc = TOC[i], k_mean = k_mean[i], k_var = k_var[i], |
|
k_min = k_min, k_max = k_max, guess = guess)$root |
|
|
|
b0 <- polynom_coefs[1] |
|
b1 <- polynom_coefs[2] |
|
b2 <- polynom_coefs[3] |
|
|
|
int_kpow2C[i] <- int_kpowC(k_min, k_max, b0, b1, b2, 2) |
|
int_kpow3C[i] <- int_kpowC(k_min, k_max, b0, b1, b2, 3) |
|
|
|
#improve estimate coefficient for next cell |
|
guess <- c(b0, b1, b2) |
|
} |
|
|
|
R.TOC <- R |
|
R.Ck <- - int_kpow2C |
|
R.Cvar <- - int_kpow3C + 2 * k_mean * int_kpow2C - k_mean^3 * TOC |
|
R.Cage <- TOC + Cage / TOC * R #aging + reaction |
|
|
|
dTOCdt <- tranTOC$dC + R.TOC |
|
dCkdt <- tranCk$dC + R.Ck |
|
dCvardt <- tranCvar$dC + R.Cvar + Cvar_mutran |
|
dCagedt <- tranCage$dC + R.Cage |
|
|
|
# Assemble the total rate of change, and return fluxes and reaction rates |
|
return(list(c(dTOCdt = dTOCdt, dCkdt = dCkdt, dCvardt = dCvardt, dCagedt = dCagedt))) |
|
}) |
|
} |
|
|
|
discrete_validation_model <- function(t, state, parms) { |
|
with(as.list(parms), |
|
with(as.list(transport_parms), { |
|
|
|
dCdt <- rep(0, length.out = length(state)) |
|
for (i_bin in 1:n_bins) { |
|
conc <- state[((i_bin-1)*N+1):(i_bin*N)] |
|
|
|
trans <- tran.1D(C = conc, C.up = Cup_bins[i_bin], D = Db$int, |
|
v = SAR, VF = 1, dx = grid) |
|
|
|
R <- - k_bins[i_bin] * conc |
|
|
|
dCdt[((i_bin-1)*N+1):(i_bin*N)] <- trans$dC + R |
|
} |
|
return(list(dCdt)) |
|
}) |
|
) |
|
} |
|
|
|
#============================================================================= |
|
# Model solution |
|
#============================================================================= |
|
|
|
# Initialize state variables for solver |
|
TOC.in <- rep( unname(parms["TOC.swi"]), length.out = N) #wght C / wght sed |
|
Ck.in <- kmean.swi * TOC.in |
|
Cvar.in <- kvar.swi * TOC.in |
|
Cage.in <- 0 * TOC.in |
|
state <- c(TOC.in, Ck.in, Cvar.in, Cage.in) |
|
|
|
#solve the model |
|
#std <- steady.1D(y = state, func = model, parms = parms, |
|
# names = c("TOC", "Ck", "Cvar"), nspec = 3) |
|
|
|
#times <- seq(0, 0.32, length.out = 10) |
|
times <- seq(0, 1000, by = 1) / time_scale #[yr] [1/yr] = [-] |
|
|
|
|
|
print(system.time( |
|
trans <- ode.1D(y = state, times = times, func = model, parms = parms, |
|
names = c("TOC", "Ck", "Cvar", "Cage"), nspec = 4, method = "vode", |
|
verbose = TRUE) |
|
)) |
|
|
|
if (attributes(trans)$istate[1] != 2) |
|
stop("numerical solution did not converge") |
|
|
|
state_trans <- trans[length(times), 2:(length(state)+1)] |
|
TOC <- state_trans[1:N] |
|
Ck <- state_trans[(N+1):(2*N)] |
|
Cvar <- state_trans[(2*N+1):(3*N)] |
|
Cage <- state_trans[(3*N+1):(4*N)] |
|
|
|
#============================================================================= |
|
# Validation model with discrete bins of OM with different reactivities |
|
#============================================================================= |
|
n_bins <- 30 #number of bins |
|
bin_width <- (k_max-k_min)/n_bins #delta k for each bin |
|
k_bins <- bin_width * ((1:n_bins) - 0.5) #center of bins |
|
Cup_bins <- vapply(1:n_bins, function(i_bin) { #upper boundary for each bin |
|
integrate(f = polynom_dist, b0 = coefs_init[1], b1 = coefs_init[2], b2 = coefs_init[3], |
|
lower = k_min + (i_bin-1) * bin_width, |
|
upper = k_min + i_bin * bin_width)$value |
|
}, FUN.VALUE = 1 ) |
|
|
|
validation_parms <- list(n_bins = n_bins, transport_parms = parms, |
|
k_bins = k_bins, Cup_bins = Cup_bins) |
|
|
|
C_init_bins <- as.vector(sapply(1:n_bins, function(i) rep(Cup_bins[i], N) )) |
|
|
|
val_state <- C_init_bins |
|
|
|
valtrans <- ode(y = val_state, times = times, func = discrete_validation_model, |
|
parms = validation_parms, method = "vode") |
|
|
|
state_val <- valtrans[length(times), 2:(length(val_state)+1)] |
|
stateval_mat <- matrix(ncol = n_bins, state_val) |
|
|
|
TOC_val <- rowSums(sapply(1:n_bins, function(i) |
|
return(state_val[ ((i-1)*N+1):(i*N)]))) |
|
Ck_val <- rowSums(sapply(1:n_bins, function(i) |
|
return(state_val[ ((i-1)*N+1):(i*N)] * k_bins[i] ))) |
|
Cvar_val <- vapply(1:N, function(i) { |
|
conc_row <- sum(stateval_mat[i,]) |
|
k_mean_row <- sum(k_bins * stateval_mat[i,]) / conc_row |
|
sum( stateval_mat[i,] * (k_bins - k_mean_row)^2 ) |
|
}, FUN.VALUE = 1) |
|
|
|
|
|
|
|
|
|
#============================================================================= |
|
# Plot results (Figure 2 in manuscript) |
|
#============================================================================= |
|
par(mfrow = c(1,4)) |
|
plot_profile(x = TOC * 100, xlim = c(0,parms["TOC.swi"])*1e2, |
|
xlab = paste("TOC", "(%)")) |
|
lines(TOC_val * 100, y = grid$x.mid * length_scale, lwd=2, lty = 3, col = 2) |
|
legend("topleft", c("continuous", "discrete"), lty = c(1,3), |
|
lwd = c(1,2), col = c(1,2), bty = "n" ) |
|
|
|
plot_profile(x = Ck/TOC / time_scale, xlim = c(k_min, k_max)/time_scale, |
|
xlab = paste("k mean", "(y-1)")) |
|
lines(Ck_val/TOC_val/time_scale, y = grid$x.mid * length_scale, col = 2, lty = 3, lwd = 2 ) |
|
|
|
plot_profile(x = sqrt(Cvar/TOC / time_scale^2), #xlim = c(0, max(sqrt(Cvar/TOC))) / time_scale, |
|
xlab = paste("sd k (y-1)")) |
|
lines(sqrt(Cvar_val/TOC_val/time_scale^2), y = grid$x.mid * length_scale, col = 2, lty = 3, lwd = 2 ) |
|
|
|
|
|
coef_top <- polynom_dist_coefs(conc = TOC[1], k_mean = Ck[1]/TOC[1], |
|
k_var = Cvar[1]/TOC[1], |
|
k_min = k_min, k_max = k_max, guess = guess.swi, |
|
rtol = 1e-8 )$root |
|
top_dist <- polynom_dist(kaxs, coef_top[1],coef_top[2], coef_top[3]) |
|
coef_bot <- polynom_dist_coefs(conc = TOC[N], k_mean = Ck[N]/TOC[N], |
|
k_var = Cvar[N]/TOC[N], |
|
k_min = k_min, k_max = k_max, guess = guess.swi, |
|
rtol = 1e-8 )$root |
|
bot_dist <- polynom_dist(kaxs, coef_bot[1], coef_bot[2], coef_bot[3]) |
|
|
|
|
|
plot(y = top_dist, x = kaxs/time_scale, |
|
xlab = "k (y-1)", ylab = "distribution", type = "l", col = "green", |
|
ylim = c(0, max(top_dist)), |
|
xlim = c(k_min, k_max) / time_scale ) |
|
|
|
lines(y = bot_dist, x = kaxs/time_scale, col = "brown" ) |
|
|
|
legend("right", legend = c("top", "bottom"), lty = c(1,1), |
|
col = c("green", "brown"), bty = "n") |
|
|
|
|
|
#============================================================================= |
|
# Additional plots |
|
#============================================================================= |
|
|
|
par(mfrow = c(1,4)) |
|
#Plot the age of organic matter over depth |
|
plot_profile(x = Cage/TOC * time_scale, |
|
xlab = paste("age", "(y)")) |
|
|
|
#Plot the Peclet number |
|
if (parms["Db_max"] > 0 ) { |
|
plot_profile(x = SAR * parms["Db_efold"] / Db$mid, xlab = "Pe" ) |
|
abline(v = 1, lty = 2) |
|
} |
|
|
|
#Plot mean reactivity vs age |
|
plot(x = Ck/TOC / time_scale, y = Cage/TOC * time_scale, type = "l", |
|
ylab = "age (y)", xlab = "mean k (y-1)", |
|
ylim = rev(range( Cage/TOC * time_scale))) |
|
|
|
#Plot variance reactivity vs age |
|
plot(x = Cvar/TOC / time_scale^2, y = Cage/TOC * time_scale, type = "l", |
|
ylab = "age (y)", xlab = "var age (y^2)", |
|
ylim = rev(range( Cage/TOC * time_scale)))
|
|
|