Scavenging processes in multicomponent medium with first-order reaction kinetics: Lagrangian and Eulerian modeling

Abstract

A process of the removal of dissolved elements in the ocean by adsorption onto settling particulate matters (scavenging) is studied analytically and using the Lagrangian and Eulerian numerical methods. The generalized model of scavenging in a multicomponent reactive medium with first-order kinetics consisting of water and multifraction suspended particulate matter was developed. Two novel numerical schemes were used to solve the transport–diffusion–reaction equations for transport dominated flows. The particle tracking algorithm based on the method of moments was developed. The modified flux-corrected transport method for the Eulerian transport–diffusion–reaction equations is a flux-limiter method based on a convex combination of low-order and high-order schemes. The flux limiters in the developed approach are obtained as an approximate solution of a corresponding optimization problem with a linear objective function. This approach allows the construction of the flux limiters with desired properties. The similarity solutions of the model equations for an idealized case of instantaneous release of reactive radionuclide on the ocean surface were obtained. It was found analytically that the dispersion of reactive contamination caused by reversible phase transition with increase of settling velocity, concentration of particulate matter and distribution coefficient can be much greater than that caused by diffusion, whereas an increase in the desorption rate results in a decrease of the dispersion caused by the phase transfer. The solutions using both numerical schemes are consistent with the analytical similarity solution even at zero diffusivity. The scavenging of the $^{239,240}$Pu that was introduced to the ocean surface due to the fallout from past nuclear weapon testing was simulated. The simulation results were in agreement with the observations in the northern Pacific. It was shown that even if the concentration of the $^{239,240}$Pu on the particulate matter does not exceed 2% of the total concentration, settling of particulate matter plays a crucial role in the vertical transport and dispersion of the reactive radionuclide. The importance of the scavenging by both the large fast-settling particles and small particles slowly settling and dissolving with depth due to the biochemical processes was demonstrated. For large particles, the “pseudodiffusivity” caused by phase transfer was 60 times greater than the diffusivity.

Appendix: Procedure for solving the system of equations (48)

We provide below a brief description of the shown in Fig. 3a procedure for solving the system of equations (48). Assume that $\varDelta t$ satisfies the inequalities

$$\begin{aligned}&\varDelta t\mathop {\max }\limits _i \left[ {(1 - \theta _2 ){r _i} + (1 - \theta _1 ){a_{ii}}} \right] \le 1 \end{aligned}$$

(73)

$$\begin{aligned}&- (1 - \theta _1 )\;\varDelta t\;\mathop {\min }\limits _i \left[ {a_{ii-1} + a_{ii + 1}} \right] \le 1 \end{aligned}$$

(74)

where $a_{ij}$ is calculated using relations (49)–(50). Then, to solve the system of linear equation (48) in one time step, we use the following iterative process:

Step 1.:

Initialize positive numbers $\delta$, $\epsilon _1$, $\epsilon _2$. Set $p=0$, $\varvec{y^{n+1,0}} = {\varvec{y}}^n$, $\varvec{\psi ^{n,0}=\psi ^{n+1,0}} = 0$ .

Step 2.:

Compute the $\varvec{\psi }^{n,p+1}$ and $\varvec{\psi }^{n+1,p+1}$ by the formulas

$$\begin{aligned} \psi _{i +1/2} = {{\left\{ \begin{array}{ll}{\min (W_i^ + ,W_{i + 1}^ - ), \quad b_{ii} > 0}\\ {\min (W_i^- ,W_{i+1}^+ ), \quad b_{ii} < 0} \end{array}\right. }} \end{aligned}$$

(75)

where $b_{ij}$ is calculated using relations (51),

$$\begin{aligned} W_i^ \pm&= {} \min \left( 1, \frac{Q_i^ \pm }{P_i^ \pm } \right) \end{aligned}$$

(76)

$$\begin{aligned} Q_i^ +&= {} \frac{1}{\varDelta t} \left( {\mathop {\max }\limits _{j \in (i - 1,i,i + 1)} y_j^n - y_i^n} \right) + (1 - \theta _1 ) \left[ {a_{ii - 1}^n \left( {y_{i - 1}^n - y_i^n} \right) } \right. \nonumber \\&\quad+ \left. { a_{ii + 1}^n \left( {y_{i + 1}^n - y_i^n} \right) } \right] \end{aligned}$$

(77)

$$\begin{aligned} Q_i^ -&= {} \frac{1}{\varDelta t} \left( {\mathop {\min }\limits _{j \in (i - 1,i,i + 1)} y_j^n - y_i^n} \right) + (1 - \theta _1 ) \left[ {a_{ii - 1}^n\left( {y_{i - 1}^n - y_i^n} \right) } \right. \nonumber \\&\quad+ \left. { a_{ii + 1}^n\left( {y_{i + 1}^n - y_i^n} \right) }\right] \end{aligned}$$

(78)

$$\begin{aligned} P_i^ +&= {} \left[ {(1 - \theta _1 )\max \left( {0,b_{ii - 1}^n} \right) + \theta _1 \max \left( {0,b_{ii - 1}^{n + 1,p}} \right) } \right] \nonumber \\&\quad+ \left[ {(1 - \theta _1 )\max \left( {0,-b_{ii}^n} \right) + \theta _1 \max \left( {0,-b_{ii}^{n + 1,p}} \right) } \right] \end{aligned}$$

(79)

$$\begin{aligned} P_i^ -&= {} \left[ {(1 - \theta _1 )\min \left( {0,b_{ii - 1}^n} \right) + \theta _1 \min \left( {0,b_{ii - 1}^{n + 1,p}} \right) } \right] \nonumber \\&\quad+ \left[ {(1 - \theta _1 )\min \left( {0,-b_{ii}^n} \right) + \theta _1 \min \left( {0,-b_{ii}^{n + 1,p}} \right) } \right] \end{aligned}$$

(80)

Step 3.:

For the ${\varvec{\psi } ^{n,p+1}}$ and ${\varvec{\psi } ^{n+1,p+1}}$, we find ${{\varvec{y}}^{n+1,p+1}}$ from the system of linear equations

$$\begin{aligned} \begin{aligned}&\left[ {E + \varDelta t \theta _2 {R ^{n + 1}} + \varDelta t \theta _1 {A^{n + 1}}} \right] {\varvec{y}}^{n+1,p+1} \\&\quad =\left[ {E - \varDelta t(1 - \theta _2 ){R ^n} - \varDelta t\,(1 - \theta _1 ){A^n}} \right] {\varvec{y}}^n \\&\qquad - \varDelta t \left( (1 - \theta _1 ){B^n}\varvec{\psi } ^{n,p + 1} + \theta _1 {B^{n + 1,p}}\varvec{\psi } ^{n+1,p+1} \right) \\&\qquad + \varDelta t \left( {\varvec{g}}^{(\theta _1 ),n} + \varvec{\varphi }^{(\theta _2 ),n} \right) \end{aligned} \end{aligned}$$

(81)

Step 4.:

Algorithm stop criterion

$$\begin{aligned} \begin{aligned}&\mathop {\max } \limits _i \frac{{\left| {y_i^{n+1,p+1} - y_i^{n+1,p}} \right| }}{{\max \left( {\delta ,\left| {y_i^{n+1,p+1}} \right| } \right) }}< {\varepsilon _1}, \quad \mathop {\max } \limits _i \left| {\psi _{i +1/2}^{n+1,p+1} - \psi _{i +1/2}^{n+1,p}} \right|< {\varepsilon _2} \quad \mathrm{and} \\&\mathop {\max } \limits _i \left| {\psi _{i +1/2}^{n,p + 1} - \psi _{i +1/2}^{n,p}} \right| < {\varepsilon _2} \end{aligned} \end{aligned}$$

(82)

If conditions (82) hold, then ${\varvec{y}}^{n+1} = {{\varvec{y}}^{n+1,p+1}}$. Otherwise, $p = p + 1$ and go to Step 2.

Remark. Note that under restriction (73), the difference scheme (46) is monotone and it satisfies the local maximum principle for $\varDelta t$ from (74) [19].

In the calculations, we can use the following second- or third-order advective and diffusive numerical fluxes

$$\begin{aligned} F_{i +1/2}^{H2}&= {} \frac{1}{2}{u_{p,i +1/2}} \; (y_i + y_{i+1}) \end{aligned}$$

(83)

$$\begin{aligned} F_{i +1/2}^{H3}&= {} u_{p,i +1/2}^ + {y_i} + u_{p,i +1/2}^ - {y_{i+1}} + \left[ {\frac{3}{8} \left| {u_{p,i +1/2}} \right| (y_{i+1} -y_i) } \right. \nonumber \\&\quad+ \left. { \frac{1}{8}u_{p,i +1/2}^ + (y_i - y_{i-1}) - \frac{1}{8}u_{p,i +1/2}^ - (y_{i+2} - y_{i+1})} \right] \end{aligned}$$

(84)

$$\begin{aligned} Q_{i +1/2}^{H3}&= {} \frac{D_{i +1/2}}{\varDelta _{i +1/2}x} \left[ {y_{i + 1} - {y_i} + \frac{1}{24}(y_{i-1} - y_{i+2}) + \frac{1}{8}(y_{i+1} - y_i)} \right] \end{aligned}$$

(85)

The numerical flux $F_{i + {1/2}}^{H2}$ has second-order accuracy in space at ${x_{i + {1/2}}}$. The numerical flux $F_{i + {1/2}}^{H3}$ is the QUICK numerical flux [17] and it has third-order accuracy in space on a uniform mesh at ${x_{i + {1/2}}}$. The numerical flux $Q_{i + {1/2}}^{H3}$ is also third-order accuracy.

The corresponding antidiffusive fluxes second- and third-order accuracies are of the form

$$\begin{aligned}&\begin{aligned} h_{i +1/2}^{ad2}&= F_{i +1/2}^{H2} + Q_{i +1/2}^L - F_{i +1/2}^L - G_{i + 1/2}^L \\&\quad= - \min \left( {0,\frac{D_{i +1/2}}{\varDelta _{i +1/2}x} - \frac{\left| {u_{p,i +1/2}} \right| }{2}} \right) (y_{i+1} - y_i) \end{aligned} \end{aligned}$$

(86)

$$\begin{aligned}&\begin{aligned} h_{i +1/2}^{ad3}&= F_{i +1/2}^{H3} + Q_{i +1/2}^{H3} - F_{i +1/2}^L - G_{i + 1/2}^L \\&\quad= - \min \left( {0,\frac{D_{i +1/2}}{\varDelta _{i +1/2}x} - \frac{\left| {u_{p,i +1/2}} \right| }{2}} \right) (y_{i + 1} - y_i) \\&\quad\quad + \frac{1}{8} \left[ {-\left| {u_{p,i +1/2}} \right| (y_{i + 1} - y_i) + u_{p,i +1/2}^ + (y_i -y_{i-1}) - u_{p,i +1/2}^ - (y_{i+2} - y_{i+1})} \right] \\&\quad\quad - \frac{D_{i +1/2}}{\varDelta _{i +1/2}x} \left[ {\frac{1}{24}(y_{i-1} - y_{i+2}) + \frac{1}{8}(y_{i+1} -y_i)} \right] \end{aligned} \end{aligned}$$

(87)

In the numerical solution of the system of Eqs. (1)–(2) $\varphi$ in (31) is a function of $C_d$ and $C_{p,i}$. In this case, Eqs. (1)–(2) are solved sequentially by a simple iterative method.