The Continuity Equation: From Continuous to Discrete

1. The Integral (Global) Form

Consider a conserved extensive quantity with volume density $\rho(\mathbf{x}, t)$ and flux density $\mathbf{J}(\mathbf{x}, t)$. Conservation demands that the rate of change of the total amount inside an arbitrary fixed control volume $\Omega$ equals the net influx through its boundary $\partial\Omega$:

$$\frac{d}{dt} \int_\Omega \rho \, dV = -\oint_{\partial\Omega} \mathbf{J} \cdot \hat{\mathbf{n}} \, dA + \int_\Omega S \, dV$$

where $\hat{\mathbf{n}}$ is the outward unit normal and $S(\mathbf{x}, t)$ is a source/sink term. The minus sign reflects the convention that $\hat{\mathbf{n}}$ points outward: a positive $\mathbf{J} \cdot \hat{\mathbf{n}}$ corresponds to an outflow, hence a decrease in the quantity inside $\Omega$.

2. The Differential (Local) Form

Since $\Omega$ is fixed in space, the time derivative passes under the integral. Applying the divergence theorem to the surface integral:

$$\oint_{\partial\Omega} \mathbf{J} \cdot \hat{\mathbf{n}} \, dA = \int_\Omega \nabla \cdot \mathbf{J} \, dV$$

we obtain:

$$\int_\Omega \left( \frac{\partial \rho}{\partial t} + \nabla \cdot \mathbf{J} - S \right) dV = 0$$

Since this must hold for every arbitrary volume $\Omega$, and the integrand is assumed continuous, the integrand itself must vanish pointwise:

$$\boxed{\frac{\partial \rho}{\partial t} + \nabla \cdot \mathbf{J} = S}$$

In the source-free case ($S = 0$), this reduces to the familiar form:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \mathbf{J} = 0$$

2.1 Expanded in Cartesian Coordinates

Writing $\mathbf{J} = (J_x, J_y, J_z)$:

$$\frac{\partial \rho}{\partial t} + \frac{\partial J_x}{\partial x} + \frac{\partial J_y}{\partial y} + \frac{\partial J_z}{\partial z} = 0$$

3. Discretisation on a Uniform Cartesian Grid

We now derive the discrete form in detail. Consider a uniform grid in $d$ dimensions with spacing $\Delta x$, $\Delta y$, $\Delta z$ and time step $\Delta t$. We label grid points by integer indices $(i, j, k)$ and time levels by $n$:

$$x_i = i \, \Delta x, \quad y_j = j \, \Delta y, \quad z_k = k \, \Delta z, \quad t^n = n \, \Delta t$$

The discrete density at grid point $(i,j,k)$ and time $t^n$ is denoted $\rho^n_{i,j,k}$.

3.1 Finite-Volume Approach (Conservative Discretisation)

The most natural discretisation proceeds from the integral form. Consider the cell $\Omega_{i,j,k}$ centred at $(x_i, y_j, z_k)$ with volume $\Delta V = \Delta x \, \Delta y \, \Delta z$. Integrating the continuity equation over this cell:

$$\frac{d}{dt} \int_{\Omega_{i,j,k}} \rho \, dV = -\oint_{\partial \Omega_{i,j,k}} \mathbf{J} \cdot \hat{\mathbf{n}} \, dA$$

Step 1: Approximate the volume integral. Assuming $\rho$ is approximately uniform over the cell:

$$\int_{\Omega_{i,j,k}} \rho \, dV \approx \rho_{i,j,k} \, \Delta V$$

Step 2: Approximate the surface integral. The cell has six faces. The flux through each face is approximated by the flux at the face centre times the face area. For the two faces perpendicular to $x$:

• Right face at $x_{i+1/2}$: area $\Delta y \, \Delta z$, outward normal $+\hat{\mathbf{x}}$, flux contribution $+J_x \big|_{i+1/2,j,k} \, \Delta y \, \Delta z$
• Left face at $x_{i-1/2}$: area $\Delta y \, \Delta z$, outward normal $-\hat{\mathbf{x}}$, flux contribution $-J_x \big|_{i-1/2,j,k} \, \Delta y \, \Delta z$

Similarly for $y$ and $z$. The total surface integral becomes:

$$\oint_{\partial \Omega} \mathbf{J} \cdot \hat{\mathbf{n}} \, dA \approx \left( J_{x,i+1/2} - J_{x,i-1/2} \right) \Delta y \, \Delta z + \left( J_{y,j+1/2} - J_{y,j-1/2} \right) \Delta x \, \Delta z + \left( J_{z,k+1/2} - J_{z,k-1/2} \right) \Delta x \, \Delta y$$

where we have suppressed the constant indices for brevity (e.g. $J_{x,i+1/2} \equiv J_x\big|_{i+1/2,j,k}$).

Step 3: Approximate the time derivative. Using a forward Euler (explicit) discretisation:

$$\frac{d}{dt}\left(\rho_{i,j,k} \, \Delta V\right) \approx \frac{\rho^{n+1}_{i,j,k} - \rho^n_{i,j,k}}{\Delta t} \, \Delta V$$

Step 4: Assemble. Dividing through by $\Delta V$:

$$\boxed{\frac{\rho^{n+1}_{i,j,k} - \rho^n_{i,j,k}}{\Delta t} + \frac{J^n_{x,i+1/2} - J^n_{x,i-1/2}}{\Delta x} + \frac{J^n_{y,j+1/2} - J^n_{y,j-1/2}}{\Delta y} + \frac{J^n_{z,k+1/2} - J^n_{z,k-1/2}}{\Delta z} = 0}$$

This is the discrete continuity equation in conservative (finite-volume) form.

3.2 Explicit Update Formula

Solving for $\rho^{n+1}_{i,j,k}$:

$$\rho^{n+1}_{i,j,k} = \rho^n_{i,j,k} - \Delta t \left[ \frac{J^n_{x,i+1/2} - J^n_{x,i-1/2}}{\Delta x} + \frac{J^n_{y,j+1/2} - J^n_{y,j-1/2}}{\Delta y} + \frac{J^n_{z,k+1/2} - J^n_{z,k-1/2}}{\Delta z} \right]$$

Defining the Courant numbers:

$$\nu_x = \frac{\Delta t}{\Delta x}, \quad \nu_y = \frac{\Delta t}{\Delta y}, \quad \nu_z = \frac{\Delta t}{\Delta z}$$

the update becomes:

$$\rho^{n+1}_{i,j,k} = \rho^n_{i,j,k} - \nu_x \left( J^n_{x,i+1/2} - J^n_{x,i-1/2} \right) - \nu_y \left( J^n_{y,j+1/2} - J^n_{y,j-1/2} \right) - \nu_z \left( J^n_{z,k+1/2} - J^n_{z,k-1/2} \right)$$

3.3 Specifying the Face Fluxes

The discretisation is completed by a prescription for the face fluxes. Several standard choices exist.

3.3.1 Central Differences (Second-Order, Non-Dissipative)

If $\mathbf{J} = \rho \, \mathbf{v}$ with a known velocity field $\mathbf{v}$, a simple choice is:

$$J_{x,i+1/2} = \frac{1}{2}\left(\rho_i v_{x,i} + \rho_{i+1} v_{x,i+1}\right)$$

This is second-order accurate in space but introduces no numerical dissipation, which can lead to spurious oscillations.

3.3.2 Upwind Scheme (First-Order, Dissipative)

For $v_x > 0$:

$$J_{x,i+1/2} = \rho_i \, v_{x,i+1/2}$$

For $v_x < 0$:

$$J_{x,i+1/2} = \rho_{i+1} \, v_{x,i+1/2}$$

Compactly:

$$J_{x,i+1/2} = v^+_{x,i+1/2} \, \rho_i + v^-_{x,i+1/2} \, \rho_{i+1}$$

where $v^+ = \max(v, 0)$ and $v^- = \min(v, 0)$.

3.3.3 Lax–Friedrichs Flux

A simple alternative that adds artificial diffusion:

$$J_{x,i+1/2}^{\text{LF}} = \frac{1}{2}\left(J_{x,i} + J_{x,i+1}\right) - \frac{\Delta x}{2\Delta t}\left(\rho_{i+1} - \rho_i\right)$$

3.4 One-Dimensional Summary

In 1D with $\mathbf{J} = \rho \, v$ and uniform velocity $v$, the upwind scheme with forward Euler gives the clean form:

$$\rho^{n+1}_i = \rho^n_i - \frac{v \, \Delta t}{\Delta x} \left( \rho^n_i - \rho^n_{i-1} \right) \quad \text{for } v > 0$$

which is stable under the CFL condition:

$$\left| \frac{v \, \Delta t}{\Delta x} \right| \leq 1$$

4. Conservation Property of the Discrete Scheme

A key feature of the finite-volume discretisation is that it is telescoping: summing the update equation over all cells, all internal face fluxes cancel in pairs, leaving only boundary contributions. Explicitly, for a 1D domain with cells $i = 1, \dots, N$:

$$\sum_{i=1}^{N} \rho^{n+1}_i \, \Delta x = \sum_{i=1}^{N} \rho^n_i \, \Delta x - \Delta t \left( J_{N+1/2} - J_{1/2} \right)$$

The total discrete "mass" changes only through the boundary fluxes $J_{N+1/2}$ and $J_{1/2}$, exactly mirroring the integral conservation law. This is the discrete analogue of conservation and holds regardless of the flux approximation chosen.

5. Higher-Order Time Integration

The forward Euler discretisation is only first-order in time. Higher-order methods can be used while preserving the conservative spatial structure.

Crank–Nicolson (implicit, second-order):

$$\frac{\rho^{n+1}_i - \rho^n_i}{\Delta t} + \frac{1}{2}\left[\left(\frac{\delta J_x}{\delta x}\right)^n + \left(\frac{\delta J_x}{\delta x}\right)^{n+1}\right] = 0$$

SSP-RK2 (explicit, second-order, strong-stability preserving):

$$\rho^{(1)}_i = \rho^n_i - \Delta t \, L(\rho^n)_i$$

$$\rho^{n+1}_i = \frac{1}{2} \rho^n_i + \frac{1}{2}\left[\rho^{(1)}_i - \Delta t \, L(\rho^{(1)})_i\right]$$

where $L(\rho)_i$ denotes the discrete spatial divergence operator evaluated at $\rho$.

6. Summary Table

| Aspect | Continuous | Discrete (Finite Volume) | |---|---|---| | Domain | $\Omega \subset \mathbb{R}^d$ | Grid cells $\Omega_{i,j,k}$ | | Density | $\rho(\mathbf{x}, t)$ | $\rho^n_{i,j,k}$ | | Flux | $\mathbf{J}(\mathbf{x}, t)$ | $J_{x,i\pm 1/2}$, etc. | | Divergence | $\nabla \cdot \mathbf{J}$ | $\frac{J_{i+1/2} - J_{i-1/2}}{\Delta x} + \cdots$ | | Time derivative | $\partial_t \rho$ | $\frac{\rho^{n+1} - \rho^n}{\Delta t}$ | | Conservation | $\frac{d}{dt}\int \rho\,dV = -\oint \mathbf{J}\cdot\hat{\mathbf{n}}\,dA$ | Telescoping sum: internal fluxes cancel |