23 Analytical solution to the general radionuclide transport equation
General form of the transport equation
The general form of the radionuclide transport equation is:
[latex]\frac{\partial}{\partial t}(\epsilon K N)+ \nabla (\epsilon v N) - \nabla^2 (\epsilon D N)+\epsilon K \lambda N =0.[/latex]
This governing equation is used not only for radionuclides, but for any dissolved contaminant in a saturated, porous geological media. This is the general version of the advection-diffusion-reaction transport equation. Depending on the discipline, different symbols may be used in the equations. For example, the radionuclide concentration in the groundwater (N) may be shown as (c). For nuclear engineering problems, the reaction constant (λ) is simply the decay constant of the radionuclide of interest.
Simplifying assumptions
In reality, parameters such as fluid porosity or advection may be both spatially and temporally dependent. The transport equation would then be solved numerically. It is preferred, however, that an analytical solution be obtained if possible, if reasonable, simplified assumptions can be made. We can assume a one dimensional spatial geometry, for example, to represent a parallel, planar fracture. If we further assume that the length of the fracture is sufficiently short such that there is not much variation in the parameters, then they can also assumed to be constant.
The transport equation reduces to:
[latex]\epsilon K \frac{\partial N}{\partial t}+\epsilon v \frac{\partial N}{\partial x}-\epsilon D \frac{\partial^2 N}{\partial x^2}+\epsilon K \lambda N=0.[/latex]
For a semi-infinite medium, the following side conditions are applied:
[latex]\displaylines{ N(x,0)=0,~~ 0 \lt x \lt \infty, \\ N(0,t)=N^*,~~ t \gt \infty, \\ N(\infty,t)=0,~~ t \gt \infty. }[/latex]
Integral transformation
For this governing equation, with these conditions, an analytical solution can be obtained. It is admittedly somewhat of an abstraction of reality, but results do give an indication of the physical behavior of radionuclide migration through a porous medium. This is very important to repository performance assessment, as we are concerned with release rates into the far field.
Apply the following integral transformation:
[latex]N(x,t)\equiv\lambda\int_0^te^{-\lambda\tau}u(x,\tau)d\tau+e^{-\lambda t}u(x,t).[/latex]
Then, compute the derivatives of the above expression:
[latex]\displaylines{ \frac{\partial N}{\partial t}=e^{-\lambda t}\frac{\partial u}{\partial t} \\ \frac{\partial N}{\partial x}=\lambda\int_0^te^{-\lambda\tau}\frac{\partial u}{\partial x}d\tau+e^{-\lambda t}\frac{\partial u}{\partial x} \\ \frac{\partial^2 N}{\partial x^2}=\lambda\int_0^te^{-\lambda\tau}\frac{\partial^2 u}{\partial x^2}d\tau+e^{-\lambda t}\frac{\partial^2 u}{\partial x^2} }[/latex]
The subsequent procedure is to substitute the derivatives into the governing equation and then apply some basic algebra and rearrange the terms.
[latex]e^{-\lambda t}K\frac{\partial u}{\partial t}+ \lambda\int_0^te^{-\lambda\tau}v\frac{\partial u}{\partial x}d\tau+ e^{-\lambda t}v\frac{\partial u}{\partial x}- \lambda\int_0^te^{-\lambda\tau}D\frac{\partial^2 u}{\partial x^2}d\tau- e^{-\lambda t}D\frac{\partial^2 u}{\partial x^2} +K\lambda^2\int_0^te^{-\lambda\tau}u d\tau+ e^{-\lambda t} K\lambda u=0.[/latex]
Then:
[latex]\displaylines{ e^{-\lambda t}(K\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial x}-D\frac{\partial^2 u}{\partial x^2}+K\lambda u) +\lambda\int_0^t[e^{-\lambda\tau}(v\frac{\partial u}{\partial x}-D\frac{\partial^2 u}{\partial x^2})]d\tau \\ +K\lambda^2\int_0^te^{-\lambda\tau}ud\tau=0. }[/latex]
Evaluate the integral for the last term:
[latex]K\lambda^2\int_0^te^{-\lambda\tau}ud\tau=- e^{-\lambda t}K\lambda u+\lambda\int_0^te^{-\lambda\tau} K\frac{\partial u}{\partial t} d\tau.[/latex]
Substitute back into the governing equation and apply some more algebra:
[latex]\displaylines{ e^{-\lambda t}(K\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial x}-D\frac{\partial^2 u}{\partial x^2}+K\lambda u-K\lambda u) \\ +\lambda\int_0^t[e^{-\lambda\tau}(K\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial x}-D\frac{\partial^2 u}{\partial x^2})]d\tau=0. }[/latex]
[latex]e^{-\lambda t}(K\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial x}-D\frac{\partial^2 u}{\partial x^2}) +\lambda\int_0^t[e^{-\lambda\tau}(K\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial x}-D\frac{\partial^2 u}{\partial x^2})]d\tau=0.[/latex]
The solution to:
[latex]K\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial x}-D\frac{\partial^2 u}{\partial x^2}=0,[/latex]
will satisfy the above mathematical expression with transformed side conditions:
[latex]\displaylines{ u(x,0)=0,~~ 0 \lt x \lt \infty, \\ u(0,t)=u^*,~~ t \gt \infty, \\ u(\infty,t)=0,~~ t \gt \infty. }[/latex]
Integral transformations are common in higher level applied math. This procedure for the radionuclide transport equation was applied because the transformed expression for the transformed variable (u) is known to have an analytical solution that can be readily obtained with the Laplace transform. It is a common mathematical model that is valid for a wide range of physical problems.
Using the Laplace transform
By definition, the Laplace transform is:
[latex]\tilde{f}(x,s)\equiv\int_0^\infty f(x,t)e^{-st}dt.[/latex]
The Laplace transform is applied to the time derivative in the transformed radionuclide transport equation to obtain:
[latex]sK\tilde u+v\frac{\partial \tilde u}{\partial x}-D\frac{\partial^2 \tilde u}{\partial x^2}=0.[/latex]
It is typical to rearrange the differential equation as:
[latex]D\frac{\partial^2 \tilde u}{\partial x^2}-v\frac{\partial \tilde u}{\partial x}-sK\tilde u=0.[/latex]
It is typical to rearrange the differential equation as:
[latex]\frac{\partial^2 \tilde u}{\partial x^2}-\frac{v}{D}\frac{\partial \tilde u}{\partial x}-\frac{sK}{D}\tilde u=0.[/latex]
The transformed side conditions are:
[latex]\displaylines{ \tilde{u}(0,s)=\frac{u^*}{s}, \\ \tilde{u}(\infty,s)=0. }[/latex]
What are the limits on s?
For this mathematical model, the solution is known to be:
[latex]\tilde u=Ae^{x(\frac{v}{2D}+\sqrt{\frac{v^2}{4D^2}+\frac{sK}{D}})}+Be^{x(\frac{v}{2D}-\sqrt{\frac{v^2}{4D^2}+\frac{sK}{D}})},[/latex]
where A and B are the constants of integration.
Because the solution must be finite it is intuitive that A = 0. Students should be able to prove this mathematically. Subsequently, with A = 0, at x = 0:
[latex]B=\frac{u^*}{s},[/latex]
the solution to the Laplace transformed model is then:
[latex]\tilde u=\frac{u^*}{s}e^{x(\frac{v}{2D}-\sqrt{\frac{v^2}{4D^2}+\frac{sK}{D}})}.[/latex]
Inverting the Laplace solution into the real time domain
Obtaining the solution to the Laplace transformed differential equation is fairly straightforward because these forms of differential equations have well known solutions. Inversion to the real time domain is not exactly as straightforward. This is why math is considered an art. Obtaining the solution to the transport equation or any similarly structured mathematical problem requires some knowledge of the actual solution in order to arrange the Laplace solution into a form that is readily invertible.
The inverse of the Laplace transform is defined as:
[latex]f(t)\equiv\frac{1}{2\pi i}\lim_{\tau\rightarrow\infty}\int_{\gamma -i\tau}^{\gamma +i\tau}e^{st}\tilde{f}(t)ds.[/latex]
This is a line integral in the complex plane (sort of), where the line integral must be computed to avoid singularities. Fortunately, for engineers, there are tables of inverse Laplace transforms, so the integral does not have to be computed all the time. However, some skill is needed in order to manipulate any given Laplace solution into a recognizable form such that the tables can be used. Which, again, requires some knowledge of the real solution.
Manipulating the Laplace solution into an invertible form
Rearrange the Laplace solution as:
[latex]\tilde u=u^*e^{(\frac{xv}{2D})} \cdot \frac{1}{s}e^{(-x\sqrt{\frac{v^2}{4D^2}+\frac{sK}{D}})}.[/latex]
Note that the first factor is only x dependent. The second factor contains the Laplace variable (s). This is a common first step in setting up the solution into an invertible form.
The next step is to eliminate the square root that contains the Laplace variable (s). Yet another integral transformation is needed.
[latex]\frac{\sqrt\pi}{2}e^{-2\sigma}=\int_0^\infty e^{(-\xi^2-\frac{\sigma^2}{\xi^2})}d\xi.[/latex]
Note that the squared terms will help to this end.
The argument in the exponential function that contains the Laplace variable (s) then can be conveniently defined as:
[latex]\sigma\equiv \frac{x}{2} \sqrt{\frac{v^2}{4D^2}+\frac{sK}{D}}.[/latex]
Substitute this all into the Laplace solution.
[latex]\tilde u=u^*e^{(\frac{xv}{2D})} \cdot \frac{2}{\sqrt{\pi}}\int_0^\infty \frac{1}{s}e^{(-\xi^2-\frac{x^2}{4\xi^2}(\frac{v^2}{4D^2}+\frac{sK}{D}))}d\xi.[/latex]
Similarly, group the terms containing the Laplace variable (s) together.
[latex]\tilde u=u^*e^{(\frac{xv}{2D})} \cdot \frac{2}{\sqrt{\pi}} \int_0^\infty e^{(-\xi^2-\frac{x^2}{4\xi^2}\frac{v^2}{4D^2})}\cdot \frac{1}{s}e^{(-\frac{Kx^2}{4D\xi^2}s)}d\xi.[/latex]
Now, the model is in an invertible form.
Inverting the Laplace solution
The last factor in the Laplace solution is a known form of the Heaviside step function.
[latex]\displaylines{ \tilde{f}(s)=\frac{1}{s}e^{-\mu s} \\ f(t)=h(t-\mu) }[/latex]
Then, the second factor in the integral can be inverted.
[latex]u=u^*e^{(\frac{xv}{2D})} \cdot \frac{2}{\sqrt{\pi}} \int_0^\infty e^{(-\xi^2-\frac{x^2}{4\xi^2}\frac{v^2}{4D^2})}h(t-\frac{Kx^2}{4D\xi^2})d\xi.[/latex]
[latex]u=u^*e^{(\frac{xv}{2D})} \cdot \frac{2}{\sqrt{\pi}} \int_{\sqrt{\frac{Kx^2}{4Dt}}}^\infty e^{(-\xi^2-\frac{x^2}{4\xi^2}\frac{v^2}{4D^2})}d\xi.[/latex]
[latex]u=u^*e^{(\frac{xv}{2D})} \cdot \frac{2}{\sqrt{\pi}} \int_{\sqrt{\frac{Kx^2}{4Dt}}}^\infty e^{(-\xi^2-\frac{x^2v^2}{16D^2\xi^2})}d\xi.[/latex]
The solution has been transformed back into the real space (x,t), but the integral needs to be evaluated now in order to obtain the formal solution for the transformed radionuclide transport solution (u).
Obtaining the transformed radionuclide transport solution
The integral in the transformed solution is also a known form:
[latex]\int e^{(-a^2\xi^2-\frac{b^2}{\xi^2})}d\xi=\frac{\sqrt{\pi}}{4a}(e^{2ab}~erf(a\xi+\frac{b}{\xi})+e^{-2ab}~erf(a\xi-\frac{b}{\xi})).[/latex]
For this problem:
[latex]a=1,~b=\frac{xv}{4D}.[/latex]
The general form of the transformed solution is then:
[latex]u=[u^*e^{(\frac{xv}{2D})}~\frac{2}{\sqrt{\pi}}][\frac{\sqrt{\pi}}{4}(e^{\frac{xv}{2D}}~erf(\xi+\frac{vx}{4D\xi})+e^{-\frac{xv}{2D}}~erf(\xi-\frac{xv}{4D\xi}))].[/latex]
Evaluate the expression by applying the integration limits. This only involves straightforward algebra, and left to the student. The transformed solution is then:
[latex]u=\frac{u^*}{2}e^{(\frac{xv}{2D})}~[e^{\frac{xv}{2D}}~erfc(\frac{Kx+vt}{\sqrt{4KDt}})+e^{-\frac{xv}{2D}}~erfc(\frac{Kx-vt}{\sqrt{4KDt}})],[/latex]
Making use of:
[latex]erfc(\mu)\equiv 1-erf(\mu).[/latex]
Solution to the radionuclide transport equation
The first step in obtaining the radionuclide transport solution makes use of the boundary condition at x = 0.
[latex]u=\frac{N^*}{2}e^{(\frac{xv}{2D})}~[e^{\frac{xv}{2D}}~erfc(\frac{Kx+vt}{\sqrt{4KDt}})+e^{-\frac{xv}{2D}}~erfc(\frac{Kx-vt}{\sqrt{4KDt}})][/latex]
Recall the initial integral transformation:
[latex]N(x,t)\equiv\lambda\int_0^te^{-\lambda\tau}u(x,\tau)d\tau+e^{-\lambda t}u(x,t).[/latex]
Therefore, the solution for (N) can be obtained.
[latex]e^{-\lambda t}u(x,t)=\frac{N^*}{2}e^{(\frac{xv}{2D}-\lambda t)}~[e^{\frac{xv}{2D}}~erfc(\frac{Kx+vt}{\sqrt{4KDt}})+e^{-\frac{xv}{2D}}~erfc(\frac{Kx-vt}{\sqrt{4KDt}})],[/latex]
[latex]e^{-\lambda t}u(x,t)=\frac{N^*}{2}[e^{(\frac{xv}{D}-\lambda t)}~erfc(\frac{Kx+vt}{\sqrt{4KDt}})+e^{-\lambda t}~erfc(\frac{Kx-vt}{\sqrt{4KDt}})].[/latex]
The solution is therefore:
[latex]N(x,t)=\frac{\lambda N^*}{2}\int_0^t [e^{(\frac{xv}{D}-\lambda\tau)}~erfc(\frac{Kx+v\tau}{\sqrt{4KD\tau}})+e^{-\lambda\tau}~erfc(\frac{Kx-v\tau}{\sqrt{4KD\tau}})]d\tau+\frac{N^*}{2}[e^{(\frac{xv}{D}-\lambda t)}~erfc(\frac{Kx+vt}{\sqrt{4KDt}})+e^{-\lambda t}~erfc(\frac{Kx-vt}{\sqrt{4KDt}})].[/latex]
