Spatial Second-Order Positive and Asymptotic Preserving Unified Gas Kinetic Schemes for Radiative Transfer Equations

Abstract

We propose two spatial second-order schemes for linear radiative transfer equations by using the idea of the unified gas kinetic scheme (UGKS) to construct the numerical boundary fluxes, and show that the proposed schemes are both positive and asymptotic preserving. The UGKS was proposed by Xu and Huang (J Comput Phys 229:7747–7764, 2010) for continuum and rarefied flows firstly, and was then applied to a linear radiative transfer equation by Mieussens in (J Comput Phys 253:138–156, 2013) where the asymptotic preserving property of UGKS is shown. Although it is asymptotic preserving, UGKS can not always keep the positivity of solutions. We first apply UGKS to discretize a linear radiative transfer equation to have a spatial second-order scheme. Then, by a detailed analysis of the numerical boundary fluxes, we are able to find the reasons why the positive preserving property of UGKS fails. Finally, we carefully employ a linear scaling limiter and a flux correction to make UGKS positive-preserving but still asymptotic-preserving. Consequently, we propose two spatial second-order positive and asymptotic preserving unified gas kinetic schemes for the linear radiative transfer equation, thus improving the earlier work (J Comput Phys 444:110546, 2021) where only a first-order positive and asymptotic scheme is developed. The proposed schemes can well capture the solution of the diffusion limit equation in optically thick regions without requiring the cell size being smaller than the photon’s mean free path, while the solution in optically thin regions can also be well resolved in a natural way. To our best knowledge, this is the first time that a spatial second-order positive and asymptotic preserving gas kinetic scheme for linear radiative transfer equations is constructed. Several numerical experiments are included to validate the spatial second-order accuracy, positive- and asymptotic-preserving properties of the proposed schemes.

Appendix

1.1 Derivation of the Time Step Constraint (4.5)

Denoting

$$\begin{aligned} f(x) =(1+e^{-x})-\frac{2}{x}(1-e^{-x}), \end{aligned}$$

we see that the coefficient function ${\tilde{d}}(\Delta t, \epsilon , \sigma _s,\sigma _a)$ can be rewritten as

$$\begin{aligned} {\tilde{d}}(\Delta t, \epsilon , \sigma _s,\sigma _a) = -\frac{\sigma _s}{2\pi \epsilon ^4 \nu ^2}f( \nu \Delta t ), \end{aligned}$$

where we recall $\nu = \frac{\sigma _s}{\epsilon ^2}+\sigma _a$. When $\nu$ is small (i.e., both $\sigma _s$ and $\sigma _a$ are small), we have by the Taylor expansion that

$$\begin{aligned} f(\nu \Delta t) \le \frac{1}{6} (\nu \Delta t)^2. \end{aligned}$$

Therefore,

$$\begin{aligned} -{\tilde{d}}(\Delta t, \epsilon , \sigma _s,\sigma _a) \le \frac{\sigma _s}{12\pi \epsilon ^4}(\Delta t)^2. \end{aligned}$$

Thus, for the optically thin regime, in order to make (4.4) hold, we can approximately require $\Delta t$ to satisfy

$$\begin{aligned} \max _{i,j} \left( \frac{\sigma _{s,i\pm 1/2,j\pm 1/2}(\Delta t)^2 }{12\pi \epsilon ^4 \sigma _{s,i,j}} \right) \le \frac{ (\Delta x \Delta y)^2 }{ 4\pi \epsilon ^2 \left( (\Delta x)^2+(\Delta y)^2 \right) }, \end{aligned}$$

that is,

$$\begin{aligned} \Delta t \le \sqrt{ \min _{i,j} \frac{\sigma _{s,i,j}}{\sigma _{s,i\pm 1/2,j\pm 1/2}} \frac{ 3\epsilon ^2 (\Delta x \Delta y)^2}{ (\Delta x)^2+(\Delta y)^2}} \le \epsilon \min _{i,j} \sqrt{ \frac{\sigma _{s,i,j}}{\sigma _{s,i\pm 1/2,j\pm 1/2}} } \cdot \frac{\sqrt{6}}{2} \min \{\Delta x, \Delta y \}. \end{aligned}$$

The time step constraint (4.5) is now verified.

1.2 Truncation Error Estimation of the UGKS

This subsection is devoted to estimating the truncation errors of the original spatial second-order UGKS. For simplicity, we make the following assumption:

Assumption 8.1

We assume that the discretization of the angular variable is exact, that is, the angular error is ignored. In addition, we assume that $\sigma _s$, $\sigma _a$ are constants, and that the solution of (2.4) is sufficiently smooth.

In the following, for convenience, we add subscript $'{_h}'$ to the previous numerical solution and fluxes so as to distinguish them from the exact ones. Therefore, the previous discretized micro and macro Eqs. (3.17)–(3.16) can be rewritten as follows:

$$\begin{aligned} I^{n+1}_{i,j,m,h}= & {} I^{n}_{i,j,m,h}+\frac{\Delta t}{\Delta x}\left( F_{i-1/2,j,m,h}-F_{i+1/2,j,m,h}\right) +\frac{\Delta t}{\Delta y}\left( H_{i,j-1/2,m,h}-H_{i,j+1/2,m,h}\right) \nonumber \\{} & {} + \Delta t \frac{\sigma _{s}}{\epsilon ^2}\left( \frac{1}{2\pi }\rho ^{n+1}_{i,j,h} - I^{n+1}_{i,j,m,h}\right) -\sigma _{a}\Delta t I^{n+1}_{i,j,m,h}, \end{aligned}$$

(8.1)

$$\begin{aligned} \rho ^{n+1}_{i,j,h}= & {} \rho ^{n}_{i,j,h}+\frac{\Delta t}{\Delta x} \left( \Phi ^{n+1}_{i-1/2,j,h}-\Phi ^{n+1}_{i+1/2,j,h}\right) +\frac{\Delta t}{\Delta y} \left( \Psi ^{n+1}_{i,j-1/2,h}-\Psi ^{n+1}_{i,j+1/2,h}\right) \nonumber \\{} & {} \quad -\sigma _{a} \Delta t\rho ^{n+1}_{i,j,h}. \end{aligned}$$

(8.2)

In order to estimate the truncation error, we first derive an alternative form of the equation (2.4), from which we can have a more intuitive understanding on the construction of the UGKS scheme. More importantly, the truncation error estimate will be based on the alternative form. For this, we integrate the Eq. (2.4) over a spatial cell (i, j) and time interval $[t_n, t_{n+1}]$, and also divide it by $\Delta x\Delta y$. For the term $\partial _x I$, it holds that

$$\begin{aligned}{} & {} \frac{1}{\Delta x\Delta y }\int _{t_n}^{t_{n+1}}\int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}\int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \partial _xI(x,y,{\varvec{\Omega }},t) dxdydt\\{} & {} \quad =\frac{1}{\Delta x\Delta y}\int _{t_n}^{t_{n+1}} \int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}} \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \partial _x I(x,y_j,{\varvec{\Omega }},t) +(y-y_j)\partial ^2_{xy}I(x,y_j,{\varvec{\Omega }},t) +O((\Delta y)^2) dxdydt\\{} & {} \quad = \frac{1}{\Delta x}\left( \int _{t_n}^{t_{n+1}} I(x_{i+1/2},y_j,{\varvec{\Omega }},t)dt- \int _{t_n}^{t_{n+1}}I(x_{i-1/2},y_j,{\varvec{\Omega }},t)dt \right) + O\left( (\Delta y)^2\Delta t\right) . \end{aligned}$$

And the term $\partial _yI$ can be handled similarly. Hence, we conclude that

$$\begin{aligned} I^{n+1}_{i,j}= & {} I^{n}_{i,j}+\frac{\xi }{\epsilon \Delta x} \int _{t_n}^{t_{n+1}} I(x_{i-1/2},y_j,{\varvec{\Omega }},t) - I(x_{i+1/2},y_j,{\varvec{\Omega }},t) dt \nonumber \\{} & {} + \frac{\eta }{\epsilon \Delta y} \int _{t_n}^{t_{n+1}} I(x_i,y_{j-1/2},{\varvec{\Omega }},t)dt- I(x_i,y_{j+1/2},{\varvec{\Omega }},t)dt\nonumber \\{} & {} + \Delta t \frac{\sigma _s}{2\pi \epsilon ^2} \breve{\rho }_{i,j} -\Delta t \nu \breve{I}_{i,j} + O\left( \frac{\eta }{\epsilon } (\Delta x)^2\Delta t + \frac{\xi }{\epsilon }(\Delta y)^2 \Delta t\right) , \end{aligned}$$

(8.3)

where $\breve{\rho }, \breve{I}$ are defined in the same way as in (3.6).

Now, in order to determine $I(x_{i\pm 1/2},y_j,{\varvec{\Omega }},t)$, we get back to (2.4). Take $y=y_j$ in (2.4), then the resulting equation can be written as

$$\begin{aligned} \frac{d}{dt} \left( e^{\nu (t-t_n)}I\left( x+\frac{\xi }{\epsilon }t, y_j, {\varvec{\Omega }}, t\right) \right) = e^{\nu (t-t_n)}\left( \frac{\sigma _s }{2\pi \epsilon ^2} \rho \left( x+\frac{\xi }{\epsilon } t, y_j, t\right) -\frac{\eta }{\epsilon }\partial _y I\left( x+\frac{\xi }{\epsilon } t, y_j, {\varvec{\Omega }}, t\right) \right) . \end{aligned}$$

Integrating the above equation from $t_n$ to t for some $t\in [t_n,t_{n+1}]$, and replacing $x+\frac{\xi }{\epsilon } t$ by x, we get

(8.4)

with

(8.5)

Letting $x=x_{i\pm 1/2}$ in (8.5) and utilizing the differential mean value theorem, we find that there is a $\zeta _x\in (x_{i-1/2}, x_{i+1/2})$, such that

$$\begin{aligned}{} & {} \int _{t_n}^{t_{n+1}}\frac{\eta }{\epsilon } \int _{t_n}^te^{-\nu (t-s)}\left( \partial _y I \left( x_{i+1/2}-\frac{\xi }{\epsilon } (t-s),y_j, {\varvec{\Omega }}, s\right) \right. \\{} & {} \qquad - \left. \partial _y I \left( x_{i-1/2}-\frac{\xi }{\epsilon } (t-s), y_j, {\varvec{\Omega }}, s\right) \right) dsdt \\{} & {} \quad = \frac{\eta }{\epsilon } \int _{t_n}^{t_{n+1}} \int _{t_n}^te^{-\nu (t-s)} \Delta x \partial ^2_{xy} I \left( \zeta _x-\frac{\xi }{\epsilon } (t-s), y_j, {\varvec{\Omega }}, s\right) dsdt\\{} & {} \quad = O(\eta ~{\tilde{g}}~ \Delta x), \end{aligned}$$

where

$$\begin{aligned} {\tilde{g}}=\frac{1}{\epsilon \nu }( \frac{1}{\nu }e^{-\nu \Delta t}+\Delta t -\frac{1}{\nu } ). \end{aligned}$$

Thus, one finds the following relation

(8.6)

In the same manner, the relation for $I(x_i,y_{j\pm 1/2},{\varvec{\Omega }},t)$ and can be obtained. Substituting (8.6) into (8.3), we obtain the desired form, which we summarize in the following lemma.

Lemma 8.1

Under Assumptions 8.1, we have the following alternative form of the Eq. (2.4) in a spatial cell (i, j) and time interval $[t_n, t_{n+1}]$:

$$\begin{aligned} I^{n+1}_{i,j}= & {} I^{n}_{i,j}+\frac{\Delta t}{\Delta x} (F_{i-1/2,j}-F_{i+1/2,j})+ \frac{\Delta t}{\Delta y}(H_{i,j-1/2}-H_{i,j+1/2})\nonumber \\{} & {} + \Delta t \frac{\sigma _s}{2\pi \epsilon ^2} \breve{\rho }_{i,j} -\Delta t \nu \breve{I}_{i,j} + O\left( \frac{\eta }{\epsilon } (\Delta x)^2\Delta t + \frac{\xi }{\epsilon }(\Delta y)^2 \Delta t + \frac{\xi \eta }{\epsilon } {\tilde{g}} \right) , \end{aligned}$$

(8.7)

where

(8.8)

Integrating (8.7) with respect to the angular variable, we obtain the macroscopic equation:

$$\begin{aligned} \rho ^{n+1}_{i,j}= & {} \rho ^{n}_{i,j}+\frac{\Delta t}{\Delta x} (\Phi _{i-1/2,j}^{n+1}-\Phi _{i+1/2,j}^{n+1})+ \frac{\Delta t}{\Delta y}(\Psi _{i,j-1/2}^{n+1}-\Psi _{i,j+1/2}^{n+1})\nonumber \\{} & {} -\Delta t \sigma _{a,i,j} \breve{\rho }_{i,j} + O\left( (\Delta x)^2\Delta t + (\Delta y)^2 \Delta t + {\tilde{g}} \right) , \end{aligned}$$

(8.9)

where

(8.10)

In view of (8.1), (8.2), and (8.7) and (8.9), we conclude that the original spatial second-order UGKS is obtained immediately by taking the discrete angular directions $\{{\varvec{\Omega }}_m = (\xi _m,\eta _m),\; m=1,2,\ldots ,M\}$, and constructing the approximation to $\breve{\rho }$ and $\breve{I}$, and also taking the piecewise linear reconstruction for both $I(\ldots ,t_n)$ and $\rho$ in (8.5).

With (8.7) and (8.9) in hand, we are able to derive the truncation error. Denoting by $\tau _I, \tau _{\rho }$ the truncation errors of the scheme (8.1), (8.2), i.e., the remainders obtained by inserting the exact solution I and $\rho$ into (8.1) and (8.2), we get

$$\begin{aligned}{} & {} \frac{ I^{n+1}_{i,j,m} - I^{n}_{i,j,m} }{\Delta t} + \frac{F_{i+1/2,j,m,h}(I,\rho )-F_{i-1/2,j,m,h}(I,\rho )}{\Delta x}\nonumber \\{} & {} \qquad +\frac{H_{i,j+1/2,m,h}(I,\rho )-H_{i,j-1/2,m,h}(I,\rho )}{\Delta y}\nonumber \\{} & {} \quad =\frac{\sigma _s}{2\pi \epsilon ^2} \rho _{i,j}^{n+1} - \nu I_{i,j,m}^{n+1} + \frac{\tau _I}{\epsilon ^2}, \end{aligned}$$

(8.11)

$$\begin{aligned}{} & {} \frac{ \rho ^{n+1}_{i,j} - \rho ^{n}_{i,j} }{\Delta t} + \frac{\Phi _{i+1/2,j,h}(I,\rho )^{n+1}-\Phi _{i-1/2,j,h}^{n+1}(I,\rho )}{\Delta x}\nonumber \\{} & {} \qquad +\frac{\Psi _{i,j+1/2,h}^{n+1} (I,\rho )-\Psi _{i,j-1/2,h}^{n+1}(I,\rho )}{\Delta y} \nonumber \\{} & {} \quad = -\sigma _a \rho _{i,j}^{n+1} + \tau _{\rho }, \end{aligned}$$

(8.12)

where we have used the notations $F_{i+1/2,j,m,h}(I,\rho )$, $H_{i,j+1/2,m,h}(I,\rho )$ and $\Phi _{i+1/2,j,h}^{n+1}(I,\rho )$, $\Psi _{i,j+1/2,h}^{n+1}(I,\rho )$ to highlight that the numerical solutions $I_h,\rho _h$ in the numerical fluxes are replaced by the exact solutions $I, \rho$. Then we can prove the following estimates of these truncation errors:

Theorem 8.1

Under Assumptions 8.1, then for the regime $\epsilon =O(1)$ and the regime $\epsilon \ll 1$, it holds that

$$\begin{aligned} \big | \tau _I \big | \le C \left( \Delta t+ (\Delta x)^2+(\Delta y)^2\right) ,\qquad \quad \big | \tau _{\rho } \big | \le C \left( \Delta t+ (\Delta x)^2+(\Delta y)^2\right) , \end{aligned}$$

(8.13)

where $C\ge 0$ is a constant independent of $\epsilon$, $\Delta x$, $\Delta y$ and $\Delta t$.

Proof

Taking the angular direction $\varvec{\Omega }_m$ in (8.7), and then subtracting (8.11) from (8.7), one has

$$\begin{aligned} \frac{\tau _I}{\epsilon ^2}= & {} \frac{ \left( F_{i+1/2,j,m} -F_{i+1/2,j,m,h}(I,\rho )\right) - \left( F_{i-1/2,j,m} -F_{i-1/2,j,m,h}(I,\rho )\right) }{\Delta x} \nonumber \\{} & {} + \frac{ \left( H_{i,j+1/2,m} -H_{i,j+1/2,m,h}(I,\rho )\right) - \left( H_{i,j-1/2,m} -H_{i,j-1/2,m,h}(I,\rho )\right) }{\Delta y} \nonumber \\{} & {} - \frac{\sigma _s}{2\pi \epsilon ^2}\left( \breve{\rho }_{i,j}-\rho _{i,j}^{n+1}\right) +\nu \left( \breve{I}_{i,j,m} - I_{i,j,m}^{n+1}\right) +O\left( \frac{\eta }{\epsilon } (\Delta x)^2 + \frac{\xi }{\epsilon }(\Delta y)^2 +\frac{\xi \eta }{\epsilon } \frac{{\tilde{g}}}{\Delta t} \right) .\nonumber \\ \end{aligned}$$

(8.14)

And also subtracting (8.12) from (8.9), we see that

$$\begin{aligned} \tau _{\rho }= & {} \frac{ \left( \Phi _{i+1/2,j}^{n+1}-\Phi _{i+1/2,j,h}^{n+1}(I,\rho )\right) - \left( \Phi _{i-1/2,j}^{n+1}-\Phi _{i-1/2,j,h}^{n+1}(I,\rho )\right) }{\Delta x} \nonumber \\{} & {} + \frac{ \left( \Psi _{i,j+1/2}^{n+1}-\Psi _{i,j+1/2,h}^{n+1}(I,\rho )\right) - \left( \Psi _{i,j-1/2}^{n+1} - \Psi _{i,j-1/2,h}^{n+1}(I,\rho )\right) }{\Delta y} \nonumber \\{} & {} +\sigma _a\left( \breve{\rho }_{i,j} - \rho _{i,j}^{n+1}\right) +O\left( (\Delta x)^2 +(\Delta y)^2 + \frac{{\tilde{g}}}{\Delta t} \right) .\nonumber \\ \end{aligned}$$

(8.15)

Recalling the definition of $\breve{\rho }$ and $\breve{I}$, and by employing the Taylor expansion, we infer that

$$\begin{aligned} \breve{\rho }_{i,j}- \rho _{i,j}^{n+1} = O(\Delta t + (\Delta x)^2 +(\Delta y)^2 ),\qquad \breve{I}_{i,j,m}- I_{i,j,m}^{n+1} = O(\Delta t + (\Delta x)^2 +(\Delta y)^2). \end{aligned}$$

(8.16)

Since under Assumptions 8.1, the macro fluxes are exact integrals of the corresponding micro fluxes with respect to the angular variable. Hence, we shall focus on the terms of the fluxes in (8.14). For the first term in (8.14), we consider the case $\xi _m<0$ only, while $\xi _m>0$ can be handled similarly. Keeping in mind that $F_{i-1/2,j,m}$ is defined in (8.8) with angular direction taken as $\varvec{\Omega }_m$. We take Taylor’s expansions to $I\left( x_{i-1/2}-\frac{\xi _m}{\epsilon }(t-t_n),y_j, {\varvec{\Omega }}, t_n \right)$ at $x=x_i$, and $\rho \left( x_{i-1/2}-\frac{\xi _m}{\epsilon }(t-s),y_j,s\right)$ at $x=x_{i-1/2}, s=t_{n+1}$ in (8.5), and then substitute (8.5) into (8.8). After a straightforward calculation, we obtain

(8.17)

where

$$\begin{aligned} \breve{F}_{i-1/2,j,m}= & {} {\tilde{\alpha }}~ \xi _m \left( I^n_{i,j,m}-\frac{\Delta x}{2}\left( \frac{\partial I}{\partial x}\right) _{i,j,m}^n\right) +{\tilde{b}}~\xi _m^2 \left( \frac{\partial I}{\partial x}\right) _{i,j,m}^n +{\tilde{c}}~\xi _m\rho _{i-1/2,j}^{n+1} + {\tilde{d}}~\xi _m^2 ~ \left( \frac{\partial \rho }{\partial x}\right) _{i,j}^{n+1}, \nonumber \\ R_{i-1/2,j,m}= & {} {\tilde{e}}\left( \frac{\partial ^2 I}{\partial x^2}\right) ^n_{i,j,m} + {\tilde{f}}\left( \frac{\partial \rho }{\partial t}\right) _{i-1/2,j}^{n+1} +{\tilde{h}}\left( \frac{\partial ^2 \rho }{\partial x^2}\right) _{i-1/2,j}^{n+1} + {\tilde{\beta }}\left( \frac{\partial ^2 \rho }{\partial x \partial t}\right) _{i-1/2,j}^{n+1} \nonumber \\{} & {} + {\tilde{\gamma }}\left( \frac{\partial ^2 \rho }{\partial t^2}\right) _{i-1/2,j}^{n+1}+ higher\text {-}order~terms, \end{aligned}$$

(8.18)

with

$$\begin{aligned} {\widetilde{e}}= & {} -\frac{\xi _m}{2\epsilon \nu \Delta t} \left\{ \left( \frac{\Delta x}{2}+\frac{\xi _m}{\epsilon }\Delta t\right) ^2 e^{-\nu \Delta t}- \frac{(\Delta x)^2}{4} + \frac{2\xi _m}{\epsilon \nu } \left( \frac{\Delta x}{2}+\frac{\xi _m}{\epsilon }\Delta t\right) e^{-\nu \Delta t} \right. \\{} & {} \left. - \frac{\xi _m\Delta x}{\epsilon \nu } +\frac{2\xi _m^2}{\epsilon ^2 \nu ^2} e^{-\nu \Delta t} - \frac{2\xi _m^2}{\epsilon ^2 \nu ^2} \right\} ,\\ {\tilde{f}}= & {} \frac{\xi _m\sigma _s }{2\pi \epsilon ^3\Delta t} \left( \frac{1}{\nu ^2} \left( \Delta t + \frac{1}{\nu }\right) (1- e^{-\nu \Delta t}) -\frac{1}{2\nu }(\Delta t)^2 - \frac{1}{\nu ^2}\Delta t \right) , \\ {\tilde{h}}= & {} \frac{\xi _m^3\sigma _s }{4\pi \epsilon ^5\Delta t} \left\{ \frac{1}{\nu ^2} (\Delta t)^2 e^{-\nu \Delta t}+ \frac{4}{\nu ^3}\Delta te^{-\nu \Delta t}+ \frac{6}{\nu ^4} e^{-\nu \Delta t} + \frac{2}{\nu ^3}\Delta t - \frac{6}{\nu ^4} \right\} ,\\ {\tilde{\beta }}= & {} \frac{\xi _m^2\sigma _s }{2\pi \epsilon ^4\Delta t} \left\{ \frac{1}{\nu ^2} (\Delta t)^2e^{-\nu \Delta t} + \frac{3}{\nu ^3}\Delta t e^{-\nu \Delta t} + \frac{3}{\nu ^4} e^{-\nu \Delta t} + \frac{(\Delta t)^2}{2\nu ^2} - \frac{3}{\nu ^4} \right\} ,\\ {\widetilde{\gamma }}= & {} \frac{\xi _m\sigma _s }{4\pi \epsilon ^3\Delta t} \left\{ \frac{(\Delta t)^3}{3\nu } + \frac{(\Delta t)^2}{\nu ^2} e^{-\nu \Delta t} +\frac{2}{\nu ^3} \Delta t e^{-\nu \Delta t}-\frac{2}{\nu ^4} + \frac{2}{\nu ^4}e^{-\nu \Delta t} \right\} . \end{aligned}$$

Based on (8.17), we make the following decomposition

$$\begin{aligned}{} & {} \frac{ \left( F_{i+1/2,j,m} -F_{i+1/2,j,m,h}(I,\rho )\right) - \left( F_{i-1/2,j,m} -F_{i-1/2,j,m,h}(I,\rho )\right) }{\Delta x}\nonumber \\{} & {} \quad = \frac{ \left( \breve{F}_{i+1/2,j,m} -F_{i+1/2,j,m,h}(I,\rho )\right) - \left( \breve{F}_{i-1/2,j,m} -F_{i-1/2,j,m,h}(I,\rho )\right) }{\Delta x} \nonumber \\{} & {} \qquad + \frac{ R_{i+1/2,j,m} - R_{i-1/2,j,m} }{\Delta x},\nonumber \\ \end{aligned}$$

(8.19)

where we further have

$$\begin{aligned} \breve{F}_{i-1/2,j,m}-F_{i-1/2,j,m,h}(I,\rho )= & {} \left( -{\tilde{\alpha }}~ \xi _m~\frac{\Delta x}{2} + {\tilde{b}}~\xi _m^2\right) ~ \left( \left( \frac{\partial I}{\partial x}\right) _{i,j,m}^n - \delta _x I_{i,j,m}^n \right) \\{} & {} +{\tilde{c}}~\xi _m \left( \rho _{i-1/2,j}^{n+1}-{\overline{\rho }}_{i-1/2,j}^{n+1} \right) + {\tilde{d}}~\xi _m^2 ~ \left( \left( \frac{\partial \rho }{\partial x}\right) _{i-1/2,j}^{n+1}-\delta _x \rho _{i-1/2,j}^{n+1} \right) , \end{aligned}$$

with ${\overline{\rho }}_{i-1/2,j}^{n+1} = \frac{\rho ^{n+1}_{i-1,j}+\rho ^{n+1}_{i,j} }{2}$.

Since the solutions are smooth enough, we find that

$$\begin{aligned} \left( \frac{\partial I}{\partial x}\right) _{i,j,m}^n - \delta _x I_{i,j,m}^n= & {} -\frac{1}{3!} \left( \frac{\partial ^3 I}{\partial x^3} \right) _{i,j,m}^n (\Delta x)^2 + O( (\Delta x)^4), \\ \rho _{i-1/2,j}^{n+1}-{\overline{\rho }}_{i-1/2,j}^{n+1}= & {} -\frac{1}{2}\left( \frac{\partial ^2 \rho }{\partial x^2}\right) _{i-1/2,j}^{n+1}\left( \frac{\Delta x}{2}\right) ^2 + O\left( (\Delta x)^4\right) , \\ \left( \frac{\partial \rho }{\partial x}\right) _{i-1/2,j}^{n+1}-\delta _x \rho _{i-1/2,j}^{n+1}= & {} -\frac{1}{3!}\left( \frac{\partial ^3 \rho }{\partial x^3}\right) _{i-1/2,j}^{n+1}\left( \frac{\Delta x}{2}\right) ^2 + O\left( (\Delta x)^4 \right) . \end{aligned}$$

Thus, the first term on the right-hand side of (8.19) can be written as follows.

$$\begin{aligned}{} & {} \frac{ \left( \breve{F}_{i+1/2,j,m} -F_{i+1/2,j,m,h}(I,\rho )\right) - \left( \breve{F}_{i-1/2,j,m} -F_{i-1/2,j,m,h}(I,\rho )\right) }{\Delta x} \nonumber \\{} & {} \quad = \left( -{\tilde{\alpha }}~ \xi _m~\frac{\Delta x}{2} + {\tilde{b}}~\xi _m^2\right) ~ \left( \frac{1}{3!} \left( \left( \frac{\partial ^3 I}{\partial x^3} \right) _{i,j,m}^n - \left( \frac{\partial ^3 I}{\partial x^3} \right) _{i+1,j,m}^n \right) \Delta x + O( (\Delta x)^3) \right) \nonumber \\{} & {} \qquad +{\tilde{c}}~\xi _m \left( \frac{1}{2} \left( \left( \frac{\partial ^2 \rho }{\partial x^2}\right) _{i-1/2,j}^{n+1}- \left( \frac{\partial ^2 \rho }{\partial x^2}\right) _{i+1/2,j}^{n+1}\right) \frac{\Delta x}{4} + O\left( (\Delta x)^3\right) \right) \nonumber \\{} & {} \qquad + {\tilde{d}}~\xi _m^2 ~ \left( \frac{1}{3!}\left( \left( \frac{\partial ^3 \rho }{\partial x^3}\right) _{i-1/2,j}^{n+1}- \left( \frac{\partial ^3 \rho }{\partial x^3}\right) _{i+1/2,j}^{n+1} \right) \frac{\Delta x}{4} + O\left( (\Delta x)^3 \right) \right) \nonumber \\{} & {} \quad = \left( -{\tilde{\alpha }}~ \xi _m~\frac{\Delta x}{2} + {\tilde{b}}~\xi _m^2\right) ~ O\left( (\Delta x)^2 \right) + {\tilde{c}}~\xi _m~ O\left( (\Delta x)^2 \right) + {\tilde{d}}~\xi _m^2 ~O\left( (\Delta x)^2\right) , \end{aligned}$$

(8.20)

where the differential mean value theorem is used. As for the second term on the right-hand side of (8.19), we use the differential mean value theorem again to obtain

$$\begin{aligned} \frac{ R_{i+1/2,j,m} - R_{i-1/2,j,m} }{\Delta x}= & {} {\tilde{e}}~O(1) +{\tilde{f}}~\xi _m~O(1) +{\tilde{h}}~\xi _m^3~O(1) +{\tilde{\beta }}~\xi _m^2~O(1)\nonumber \\{} & {} + {\tilde{\gamma }}~\xi _m~O(1) + higher\text {-}order~terms. \end{aligned}$$

(8.21)

We remark that here the “higher-order terms” can be handled in the same manner, and we shall not repeat to mention it in what follows.

For the regime $\epsilon = O(1)$, we take the Taylor expansion to $e^{-\nu \Delta t}$ in the coefficients to obtain

$$\begin{aligned}{} & {} {\tilde{\alpha }}=O(1), \qquad {\tilde{b}}=O\left( \Delta t\right) , \qquad {\tilde{c}} = O(\Delta t),\qquad {\tilde{d}} = O\left( (\Delta t)^2)\right) ,\quad {\tilde{e}} = O\left( (\Delta x)^2 + \Delta x\Delta t)\right) ,\qquad \\{} & {} \qquad {\tilde{f}} = O\left( (\Delta t)^2\right) ,\qquad {\tilde{h}} = O\left( (\Delta t)^3\right) , \qquad {\tilde{\gamma }} = O\left( (\Delta t)^3\right) , \qquad {\tilde{\beta }} = O\left( (\Delta t)^3\right) . \end{aligned}$$

Hence, combining (8.19) and (8.20) with (8.21), we conclude that

$$\begin{aligned} \Big |\frac{ \left( F_{i+1/2,j,m} -F_{i+1/2,j,m,h}(I,\rho )\right) - \left( F_{i-1/2,j,m} -F_{i-1/2,j,m,h}(I,\rho )\right) }{\Delta x} \Big | \le C\left( (\Delta t)^2 +(\Delta x)^2 \right) .\nonumber \\ \end{aligned}$$

(8.22)

Similarly, one can obtain

$$\begin{aligned} \Big |\frac{ \left( H_{i,j+1/2,m} -H_{i,j+1/2,m,h}(I,\rho )\right) - \left( H_{i,j-1/2,m} -H_{i,j-1/2,m,h}(I,\rho )\right) }{\Delta y} \Big | \le C\left( (\Delta t)^2 +(\Delta y)^2\right) .\nonumber \\ \end{aligned}$$

(8.23)

Substituting (8.16), (8.22) and (8.23) into (8.14), and keeping in mind that ${\tilde{g}} = O\left( (\Delta t)^2\right)$ in the regime $\epsilon =O(1)$, we deduce that $\big | \tau _I \big |\le C\left( \Delta t+ (\Delta x)^2+(\Delta y)^2\right)$. And the estimate $\big |\tau _{\rho }\big | \le C\left( \Delta t+ (\Delta x)^2+(\Delta y)^2\right)$ follows immediately.

On the other hand, for the regime $\epsilon \ll 1$, we can verify that

$$\begin{aligned}{} & {} {\tilde{\alpha }}=O(\epsilon ), \qquad {\tilde{b}}=O(\epsilon ^2), \qquad {\tilde{c}}=O(\frac{1}{\epsilon }),\qquad {\tilde{d}}=O(1),\qquad {\tilde{e}}=O(\epsilon ),\qquad \nonumber \\{} & {} \quad {\tilde{f}}=O(\frac{1}{\epsilon }),\qquad {\tilde{h}} = O\left( \epsilon \right) , \qquad {\tilde{\gamma }} = O\left( \frac{1}{\epsilon }\right) , \qquad {\tilde{\beta }} = O\left( \Delta t\right) . \end{aligned}$$

(8.24)

Thus,

$$\begin{aligned} \epsilon ^2 \Big |\frac{ \left( F_{i+1/2,j,m} -F_{i+1/2,j,m,h}(I,\rho )\right) - \left( F_{i-1/2,j,m} -F_{i-1/2,j,m,h}(I,\rho )\right) }{\Delta x} \Big | \le C\left( \Delta t +(\Delta x)^2\right) , \end{aligned}$$

where we have used the fact that $\epsilon < \Delta t$ in the regime $\epsilon \ll 1$. The estimate for the $y-$direction flux can be obtained in the same manner. Therefore, we have $\big | \tau _I \big | \le C\left( \Delta t+ (\Delta x)^2+(\Delta y)^2\right)$.

Now, it remains to prove the estimate for $\tau _{\rho }$. Using (3.1), (8.20) and (8.21), one gets

$$\begin{aligned}{} & {} \frac{ \left( \Phi _{i+1/2,j}^{n+1}-\Phi _{i+1/2,j,h}^{n+1}(I,\rho )\right) - \left( \Phi _{i-1/2,j}^{n+1}-\Phi _{i-1/2,j,h}^{n+1}(I,\rho )\right) }{\Delta x}\\{} & {} \quad = \left( -{\tilde{\alpha }}~\frac{\Delta x}{2} + {\tilde{b}} \right) ~ O\left( (\Delta x)^2 \right) + {\tilde{d}} ~O\left( (\Delta x)^2\right) + {\tilde{e}}~O(1) +{\tilde{\beta }}~O(1) + higher\text {-}order~terms. \end{aligned}$$

In view of (8.24) and the fact that ${\tilde{g}}=O(\epsilon \Delta t)$ in the regime $\epsilon \ll 1$, we conclude $\big |\tau _{\rho }\big | \le C\left( \Delta t+ (\Delta x)^2+(\Delta y)^2\right)$. This completes the proof. $\square$