Numerically solving the Advection-diffusion equation with no-flux boundary condition leads to violation of mass conservation

Question

I am trying to solve numerically the advection-diffusion equation of the following form

$$\frac{\partial C}{\partial t}=\alpha\frac{\partial^2 C}{\partial x^2}+\beta \frac{\partial C}{\partial x}$$

with no flux BCs namely

$$\alpha\frac{\partial C}{\partial x}+\beta C=0$$

at $x=0$ and $x=L$.

I solve this as one would normally solve the diffusion equation, namely by using finite differences and the forward Euler algorithm.

$$C_N(t+dt)=C_N(t)+dt\left(\alpha\frac{C_{N+1}(t)-2C_N(t)+C_{N-1}(t)}{dx^2}+\beta\frac{C_{N+1}(t)-C_{N-1}(t)}{2dx}\right)$$

The BC's I implement via the ghost point method i.e. $N=n+1$ and $N=0$ are excluded and eliminated via the boundary conditions. For $N=1$ we get

$$\alpha\frac{C_{2}-C_{0}}{2dx}+\beta C_1=0$$

or

$$C_{0}=C_{2}+\frac{2dx\beta}{\alpha}C_1$$

and for $N=n$

we get similarly

$$C_{n+1}=C_{n-1}-\frac{2dx\beta}{\alpha}C_n$$

which gives us the time-evolution equations at the boundary as

$$C_1(t+dt)=C_1(t)+dt\left(\alpha\frac{2C_{2}(t)-2C_1(t)+\frac{2dx\beta}{\alpha}C_1(t)}{dx^2}-\frac{\beta^2}{\alpha}C_1\right)$$

and

$$C_n(t+dt)=C_n(t)+dt\left(\alpha\frac{2C_{n-1}(t)-2C_n(t)-\frac{2dx\beta}{\alpha}C_n(t)}{dx^2}-\frac{\beta^2}{\alpha}C_n\right)$$.

Now when I try to solve this numerically mass is not conserved, but the mass non-conservation can also easily be seen by considering the system initially being of homogenous concentration $C_0$. During the first iteration, all the interior points remain fixed while at the boundary we get

$$C_n(t+dt)=C_0+dt\left(\alpha\frac{-\frac{2dx\beta}{\alpha}C_0}{dx^2}-\frac{\beta^2}{\alpha}C_0\right)$$

and

$$C_1(t+dt)=C_0+dt\left(\alpha\frac{\frac{2dx\beta}{\alpha}C_0}{dx^2}-\frac{\beta^2}{\alpha}C_0\right)$$.

Now this clearly does not conserve the mass of the system as it should per the no-flux boundary conditions. What is causing this?

Edit: By mass here I mean the integral of the concentration

$$\int_0^L C(x) dx$$

which is calculated using the trapezoidal rule.

Edit 2: Here is how the system currently evolves over time

The code can also be found here

# 1. Import libraries
import numpy as np
import matplotlib.pyplot as plt

# 2. Set up parameters
nx = 101
dx = 0.01
nt = 20000
alpha = 0.1
beta=-1
dt = 0.0001

# 3. Initial conditions

q = np.ones(nx)
qn= q.copy()

# 4. Solving the PDE
for n in range(nt):

    qn=q.copy()

    q[1] = qn[1] + alpha * dt / dx ** 2 * (qn[2]-2*qn[1]+qn[2]+2*dx*(beta/alpha)*qn[1]) +\
          ( beta*dt/(2*dx))*(qn[2] - (qn[2]+2*dx*(beta/alpha)*qn[1]))
    for i in range(2, nx-2):
        
        
        q[i] = qn[i] + alpha * dt / dx ** 2 * (qn[i+1]-2*qn[i]+qn[i-1]) +\
          ( beta*dt/(2*dx))*(qn[i+1] - qn[i-1])
        
    q[nx-2] = qn[nx-2] + alpha * dt / dx ** 2 * (qn[nx-3]-2*qn[nx-2]+(qn[nx-3]-(beta/alpha)*2*dx*qn[nx-2])) +\
           beta*dt/(2*dx)*(-qn[nx-3] +(qn[nx-3]-(beta/alpha)*2*dx*qn[nx-2]) )


# 5. Plot results
plt.plot(np.linspace(0, 1, 101), q)
plt.xlim(dx, 1-dx)
plt.show()

Answer 1

For $\alpha, \beta = \mathrm{const}$ you can rewrite your equation

$$ \frac{\partial C}{\partial t}=\alpha\frac{\partial^2 C}{\partial x^2}+\beta \frac{\partial C}{\partial x}, \qquad (1), $$

to

$$ \frac{\partial C}{\partial t}=\frac{\partial\left(\alpha \frac{\partial C}{\partial x} +\beta C\right)}{\partial x},\qquad (2). $$

With $(2)$ you can see that a conservative numerical method requires a consistent approximation of flux and boundary condition

$$ f=-\left(\alpha\frac{\partial C}{\partial x}+\beta C\right), $$

in order to preserve

$$ \frac{\partial C}{\partial t}+\frac{\partial f}{\partial x} = 0. $$