<jupyter_start><jupyter_text>title = "(Advection-)diffusion-sorption-decay problem"date = "2022-03-09"author = "<NAME>, <NAME>, <NAME>, <NAME>, <NAME>"notebook = "DiffusionSorptionDecay.ipynb"[menu]  [menu.benchmarks]    parent = "Hydro-Component"<jupyter_code>import os
import ogs6py
import vtuIO
import numpy as np
from scipy import special
import matplotlib.pyplot as plt
from matplotlib.pyplot import cm<jupyter_output><empty_output><jupyter_text>1. Diffusion-sorption-decay problem  1.1 Problem description This benchmark is meant to model radionuclide migration through a semi-infinite column of undisturbed Opalinus Clay (OPA), which could potentially serve as a geological repository for radioactive waste. In this benchmark, Cesium ($^{135}$Cs) is chosen for investigation because it may be present in the repository if the disposal canister is ruptured. In view of the low hydraulic conductivity of OPA, Cesium migration is considered to be dominated by diffusion. Two additional processes are also considered, including (i) retardation resulting from solute sorption on the clay, and (ii) concentration attenuation caused by radioactive decay over time.In the one-dimensional model setup, the mass balance equation for the diffusion-sorption-decay process is given as follows:$$\frac{\partial \left( \phi R c\right)}{\partial t} = \frac{\partial}{\partial x} \left(\phi \mathrm{D_p} \frac{\partial c}{\partial x} \right) - \phi \alpha R c,$$where $c$ [mol/m$^3$] represents the $^{135}$Cs concentration in this example, $\mathrm{D_p}$ [m$^2$/s] is the pore diffusion coefficient for $^{135}$Cs, $\phi$ [-] is the porosity of OPA, $R$ [-] is the retardation factor defined as$$R = 1 + \rho k_{d} / \phi,$$with the dry density of OPA $\rho$ [kg/m$^3$] and the distribution coefficient $k_d$ [m$^3$/kg], and $\alpha$ [1/s] is the first-order decay constant, i.e.:$$\alpha = \mathrm{ln} 2 \, / \, t_{1/2},$$where $t_{1/2}$ [s] is the half life of the diffusing radionuclide.  1.2 Model setups Our simulation domain is set to be 20 meters long, which is sufficient to ensure that the concentration profile does not reach the right-hand boundary during simulation. The simulated time is one million years. Initially, the entire domain is assumed to be solute free, i.e., $c_{\mathrm{ini}}(x, 0) = 0$. The inlet concentration is held at 1 mol/L throughout the simulation, i.e., $c(0, t) = 1$ mol/L. The spatial domain is discretized by linear line elements with a length of 0.01 meter each. The time step size of one thousand years is used in the simulation. The linearized governing equation is iteratively solved using the Newton-Raphson method.The table below summarizes the parameters used in the simulation.| Parameter | Value | Unit || :-: | :-: | :-: || Porosity $\phi$ | 0.12 | - || Pore diffusion coefficient $\mathrm{D_p}$ | 8.33e-11 | m$^2$/s || OPA bulk density $\rho$ | 2394 | kg/m$^3$ || Distribution coefficient $k_{d}$ | 0.5 | m$^3$/kg || $^{135}$Cs half life $t_{1/2}$ | 2.3e6 | year || Time step size $\Delta t$ | 1e3 | year || Grid size $\Delta x$ | 0.01 | m|Notes: The parameter values are sourced from Nagra (2002).  1.3 Results **Analytical solution**The time-dependent concentration profile can be analytically calculated as a function of $t$ and $x$. The analytical solution for this diffusion-soprtion-decay system is given in Carslaw and Jaeger (1959):$$c(x, t) = \frac{c_0}{2} \left\{ \mathrm{exp} \left[-x \left(\frac{\alpha R}{\mathrm{D_p}}\right)^{\frac{1}{2}} \right] \mathrm{erfc} \left[ \frac{x}{2} \left(\frac{R}{\mathrm{D_p}t}\right)^{\frac{1}{2}} - \left(\alpha t\right)^{\frac{1}{2}}\right] + \mathrm{exp} \left[x \left(\frac{\alpha R}{\mathrm{D_p}}\right)^{\frac{1}{2}} \right] \mathrm{erfc} \left[ \frac{x}{2} \left(\frac{R}{\mathrm{D_p}t}\right)^{\frac{1}{2}} + \left(\alpha t\right)^{\frac{1}{2}}\right]\right\}$$Below a python script is presented, through which the concentration profiles can be analytically calculated at any time. Here the concentration profiles are illustrated at $t$ = 10$^3$, 10$^4$, 10$^5$, and 10$^6$ years.<jupyter_code>###Model parameters###
#Effective diffusion coefficient [m2/s]
De = 1e-11
#Porosity [-]
phi = 0.12
#Pore diffusion coefficient [m2/s]
Dp = De / phi
#Porous medium bulk density [kg/m3]
rho = 2.394e3
#Distribution coefficient [m3/kg]
Kd = 0.5
#135-Cs Half-life [year]
half_life = 2.3e6
#Decay constant [1/s]
alpha = np.log(2)/half_life/3.1536e7 # unit conversion from year to second

###Spatial and temporal discretization###
#Distance [m]
x = np.linspace(0, 2, num=201)
#Time [year]
time = np.array([1e3, 1e4, 1e5, 1e6])

###Initial condition and boundary conditions###
#Initial condition [mol/L]
c_ini = 0
#Inlet concentration [mol/L]
c_inlet = 1

###Intermediate parameters###
#Retardation factor [-]
R = 1 + rho*Kd/phi

###Analytical solution###
c = np.empty((0, x.size))
for t in time*3.1536e7: #unit conversion from year to second
    c_t = c_inlet/2*(np.exp(-x*(alpha*R/Dp)**0.5)*special.erfc(x/2*(R/Dp/t)**0.5-(alpha*t)**0.5) \
               + np.exp(x*(alpha*R/Dp)**0.5)*special.erfc(x/2*(R/Dp/t)**0.5+(alpha*t)**0.5))
    c = np.vstack([c, c_t])<jupyter_output><empty_output><jupyter_text>The analytically computed $^{135}$Cs concentration profiles at $t$ = 10$^3$, 10$^4$, 10$^5$, and 10$^6$ years are plotted as shown in the figure below.<jupyter_code>#plot analytical solution
def plot_analytical_solutions():
    fig, ax = plt.subplots()
    
    ax.set_xlim((0,2))
    ax.set_ylim((0,1))
    
    plt.xlabel('Distance [m]')
    plt.ylabel('$^{135}$ Cs concentration [mol/L]')

    color_map=iter(cm.rainbow(np.linspace(0,1,len(time))))
    
    for c_t, t, color in zip(c, time, color_map):
        ax.plot(x, c_t, linestyle='-', lw=1.5,
                label=str(np.format_float_scientific(t))+' years',
                c=color, zorder=10, clip_on=False)
    
    ax.legend(frameon=False, loc='upper right', numpoints=1, 
              fontsize=12, ncol=1)
    
    ax.xaxis.grid(color='gray', linestyle='dashed')
    ax.yaxis.grid(color='gray', linestyle='dashed')
    
plot_analytical_solutions()<jupyter_output><empty_output><jupyter_text>**Numerical solution** Correspondingly, the OGS input files of this 1D mass transport model can be found here.Then, the numerical solution of OpenGeoSys is plotted against the analytical solution for comparison.<jupyter_code>#Run OGS simulation
prj_name = "1D_DiffusionSorptionDecay"
data_dir = os.environ.get('OGS_DATA_DIR', '../../../Data')
prj_file = f"{data_dir}/Parabolic/ComponentTransport/DiffusionSorptionDecay/{prj_name}.prj"
out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')

if not os.path.exists(out_dir):
    os.makedirs(out_dir)
os.chdir(out_dir)

print(f"ogs {prj_file} > out.txt")
! ogs {prj_file} > out.txt

#Read simulation results
pvdfile = vtuIO.PVDIO(f"{prj_name}.pvd", dim=1)

def plot_simulation_results():
    fig, ax = plt.subplots()
    
    ax.set_xlim((0,2))
    ax.set_ylim((0,1))
        
    plt.xlabel('Distance [m]')
    plt.ylabel('$^{135}$ Cs concentration [mol/L]')
    
    color_map = iter(cm.rainbow(np.linspace(0,1,len(time))))
    #Plot analytical solutions 
    for c_t, color in zip(c, color_map):
        ax.plot(x, c_t, linestyle='-', lw=1.5,
                c=color, zorder=10, clip_on=False)
        
    #Add simulation results
    color_map=iter(cm.rainbow(np.linspace(0,1,len(time))))
    for t, color in zip(time, color_map):
        c_t = pvdfile.read_set_data(t*3.1536e7, 'Cs', data_type="point", pointsetarray=[(i,0,0) for i in x])
        plt.plot(x, c_t, label="Sim. "+str(np.format_float_scientific(t))+' years', 
                 color=color, marker='o', markevery=5, linestyle="", zorder=10, clip_on=False)
    
    ax.legend(frameon=False, loc='upper right', numpoints=1, 
              fontsize=12, ncol=1)
    
    ax.xaxis.grid(color='gray', linestyle='dashed')
    ax.yaxis.grid(color='gray', linestyle='dashed')
    
plot_simulation_results()<jupyter_output>ogs ../../../Data/Parabolic/ComponentTransport/DiffusionSorptionDecay/1D_DiffusionSorptionDecay.prj > out.txt
/bin/bash: ogs: command not found<jupyter_text>In the figure above, it can be seen that the numerical profiles and analytical solutions agree well at different times. **Error estimation**To evaluate the accuracy of numerical solution, the Euclidean norm of the approximate solution error $|| \mathbf{c} - \mathbf{c}^{\mathrm{exact}}||_{2} = \left(\sum\mathop{}_{\mkern-5mu i} |c_i - c_{i}^{\mathrm{exact}}|^2 \right)^{1/2}$ over the entire domain is calculated based on the exact analytical solution $\mathbf{c}^{\mathrm{exact}}$.<jupyter_code>l2_norm_error = np.empty((0, 1))

for c_ext, t in zip(c, time):
    c_sim = pvdfile.read_set_data(t*3.1536e7, 'Cs', data_type="point", pointsetarray=[(i,0,0) for i in x])
    
    l2_norm_error_t = np.log10(np.sum((c_sim - c_ext)**2)**0.5)
    l2_norm_error = np.vstack([l2_norm_error, l2_norm_error_t])
def plot_l2_norm_error():
    fig, ax = plt.subplots()
    
    ax.set_xlim((0,1e6))
    ax.set_ylim((-4,0))
    
    plt.xlabel('Time [year]')
    plt.ylabel('Log $||\mathbf{c}-\mathbf{c^{exact}}||_{2}$')
    
    ax.plot(time, l2_norm_error, linestyle='-', lw=1.5, 
            marker='o', zorder=10, clip_on=False)
    
    ax.xaxis.grid(color='gray', linestyle='dashed')
    ax.yaxis.grid(color='gray', linestyle='dashed')
plot_l2_norm_error()<jupyter_output><empty_output><jupyter_text>As the simulation time increases, the L$_2$ error norm decreases rapidly and evenually drops down to 10$^{-3}$.  2. Advection-diffusion-sorption-decay problem  2.1 Problem description From the preceding example, we further incorporate advective transport within OPA. The mass balance equation is thus extended into$$\frac{\partial \left( \phi R c\right)}{\partial t} = \frac{\partial}{\partial x} \left(\phi \mathrm{D_p} \frac{\partial c}{\partial x} - qc \right) - \phi \alpha R c,$$with the Darcy velocity $q$ [m/s]. Note that the effect of velocity-dependent dispersion is neglected.  2.2 Model setups The parameters used in this example are kept the same as in the preceding one (see the table below), except we addtionally consider a constant Darcy velocity of 2e-11 m/s through OPA.| Parameter | Value | Unit || :-: | :-: | :-: || Porosity $\phi$ | 0.12 | - || Pore diffusion coefficient $\mathrm{D_p}$ | 8.33e-11 | m$^2$/s || OPA bulk density $\rho$ | 2394 | kg/m$^3$ || Distribution coefficient $k_{\mathrm{d}}$ | 0.5 | m$^3$/kg || $^{135}$Cs half life $t_{1/2}$ | 2.3e6 | year || Darcy velocity $q$ | 2e-11 | m/s | | Time step size $\Delta t$ | 1e3 | year || Grid size $\Delta x$ | 0.01 | m|  2.3 Results **Analytical solution**Likewise, we find the analytical solution for this advection-diffusion-soprtion-decay system in Van Genuchten (1981), which is given as:$$c(x,t) = c_0 H(x,t) + M(x,t),$$where$$H(x,t) = \frac{1}{2}\mathrm{exp}\left[\frac{\left(v-u\right)x} {2D_{\mathrm{p}}}\right]\mathrm{erfc}\left[\frac{Rx-ut}{2\left(D_{\mathrm{p}}Rt\right)^{1/2}}\right] + \frac{1}{2}\mathrm{exp}\left[\frac{\left(v+u\right)x} {2D_{\mathrm{p}}}\right]\mathrm{erfc}\left[\frac{Rx+ut}{2\left(D_{\mathrm{p}}Rt\right)^{1/2}}\right],$$$$M(x,t) = -c_{\mathrm{ini}} \mathrm{exp}\left(-\frac{\mu t}{R}\right)\left\{\frac{1}{2}\mathrm{erfc}\left[\frac{Rx-vt}{2\left(D_{\mathrm{p}}Rt\right)^{1/2}}\right] + \frac{1}{2}\mathrm{exp}\left(\frac{vx}{D_{\mathrm{p}}}\right)\mathrm{erfc}\left[\frac{Rx+vt}{2\left(D_{\mathrm{p}}Rt\right)^{1/2}}\right]\right\} + c_i \mathrm{exp}\left(-\frac{\mu t}{R}\right),$$and$$u = v \left(1 + \frac{4\mu D_{\mathrm{p}}}{v^2} \right)^{1/2},$$with the pore water velocity $v$ [m/s], i.e., $v = q \,/ \phi$, and $\mu = \alpha R$. With the python script provided below, the $^{135}$Cs concentration profiles at $t$ = 10$^3$, 10$^4$, 10$^5$, and 10$^6$ years are analytically computed.<jupyter_code>#Auxilary functions
def H(x, t):
    return 0.5*np.exp((v-u)*x/2/D)*special.erfc((R*x-u*t)/2/(D*R*t)**0.5) \
    + 0.5*np.exp((v+u)*x/2/D)*special.erfc((R*x+u*t)/2/(D*R*t)**0.5)

def M(x, t):
    return - c_ini*np.exp(-mu*t/R)*(0.5*special.erfc((R*x-v*t)/2/(D*R*t)**0.5)\
    + 0.5*np.exp(v*x/D)*special.erfc((R*x+v*t)/2/(D*R*t)**0.5)) \
    + c_ini*np.exp(-mu*t/R)

###Input parameters###
#Effective diffusion coefficient [m2/s]
De = 1e-11
#Porosity [-]
phi = 0.12
#Pore diffusion coefficient [m2/s]
D = De / phi
#Porous medium bulk density [kg/m3]
rho = 2.394e3
#Distribution coefficient [m3/kg]
Kd = 0.5
#Retardation factor [-]
R = 1 + rho*Kd/phi
#135-Cs Half-life [year]
half_life = 2.3e6
#Decay constant [1/s]
k = np.log(2)/half_life/3.1536e7 # unit conversion from year to second
#Include advective mechansim
#Darcy velocity [m/s]
q = 2e-11
#Pore water velocity [m/s]
v = q / phi

###Spatial and temporal discretization###
#Distance [m]
x = np.linspace(0, 2, num=201)
#Time [year]
time = np.array([1e3, 1e4, 1e5, 1e6])

###Initial condition and boundary conditions###
#Initial condition [mol/L]
c_ini = 0
#Inflow concentration [mol/L]
c0 = 1

###Intermediate parameters###
mu = k * R
u = v*(1+4*mu*D/v**2)**0.5

###Analytical solution###
c = np.empty((0, x.size))
for t in time*3.1536e7: #unit conversion from year to second
    c_t = c0*H(x, t) + M(x, t)
    c = np.vstack([c, c_t])

plot_analytical_solutions()<jupyter_output><empty_output><jupyter_text>**Numerical solution** Also, a 1D mass transport model is established. The corresponding OGS input files can be found here.Then, the numerical solution is plotted against the analytical solution for comparison.<jupyter_code>#Run OGS simulation
prj_name = "1D_AdvectionDiffusionSorptionDecay"
data_dir = os.environ.get('OGS_DATA_DIR', '../../../Data')
prj_file = f"{data_dir}/Parabolic/ComponentTransport/AdvectionDiffusionSorptionDecay/{prj_name}.prj"
out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')

if not os.path.exists(out_dir):
    os.makedirs(out_dir)
os.chdir(out_dir)

print(f"ogs {prj_file} > out.txt")
! ogs {prj_file} > out.txt

#Read simulation results
pvdfile = vtuIO.PVDIO(f"{prj_name}.pvd", dim=1)

#Plot simulation results
plot_simulation_results()<jupyter_output>ogs ../../../Data/Parabolic/ComponentTransport/AdvectionDiffusionSorptionDecay/1D_AdvectionDiffusionSorptionDecay.prj > out.txt
/bin/bash: ogs: command not found
{'fraction_non_alphanumeric': 0.17267246214245655, 'fraction_numerical': 0.03287997756590017, 'mean_word_length': 4.44881588999236, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '322886', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5482, 'n_tokens_neox': 5126, 'n_words': 1341, 'path': 'Tests/Data/Parabolic/ComponentTransport/DiffusionSorptionDecay/DiffusionSorptionDecay.ipynb', 'repo_name': 'garibay-j/ogs', 'size': '14306'}