Mass conservation in atmospheric continuity equation numerical solution

My phd project is heavily related to numerical modeling of planetary atmospheres. In particular now I am dealing with a particular expression of the continuity equation, involving a thermodynamic flux. Equations read

$$\frac{\partial n}{\partial t} = -\frac{\partial \Phi}{\partial z}$$

$$\Phi = -nD\left(\frac{1}{n}\frac{\partial n}{\partial z} + \frac{1+\alpha}{T}\frac{\partial T}{\partial z} + \frac{1}{h}\right) -nK\left(\frac{1}{n}\frac{\partial n}{\partial z} + \frac{1}{T}\frac{\partial T}{\partial z} + \frac{1}{H}\right)$$

Where of course $n = n\left(z,t\right)$ refers to the particular chemical species density the equation solves. $h$ is the atmospheric scale height for the particular molecule while H is the average atmospheric scale height corresponding to the average atmospheric molecule. $D$ and $K$ are the molecular and eddy diffusion coefficients respectively. $T$ is the temperature. $\alpha$ is a parameter related to heat diffusion.

I use a standard implicit method to find a system of linear equations so that its solution yields the density at a new time step. I solve a tridagonal matrix in Fortran basically. For 1º order derivatives I take a 1º order central finite differences and 2º order central finite differences for second order derivatives. At the boundaries I use forward and backward approaches when necessary.

As I stated the solution seems correct when studying the density profiles. But there is something very annoying when at the boundaries fixed flux conditions are applied. Let's see by integrating the equation over all the domain

$$\int_{z_{bot}}^{z_{top}}dz\left(z,t\right) = \rho \left(t\right) = - \int_{z_{bot}}^{z_{top}}dz\frac{\partial \Phi}{\partial z}$$

$$\implies \rho \left(t\right) = - \Phi_{top} + \Phi_{bottom}$$

Being $\rho\left(t\right)$ the magnitude known as column density. After this long introduction here my doubt comes. I fix $\Phi_{bottom} = \Phi_{top} = 0$ so that we should expect $\rho\left(t\right) \equiv constant$. Nevertheless this is not the case at all, mass is not conserved in my numerical system by very big numbers, integrating for each molecule. The processing of results is undertaken by Python scripts using numpy when necessary.

I have been struggling with this for a month and I draw the following few guesses;

Could be that some parameters are not well treated. For instance there is a point where I need to take the numerical derivative of a big set of magnitudes. When studying the outcome it is very very fluctuating for a smooth curve.

Since the system deals with an atmosphere, the background profile used is a hydrostatic kind one, with densities ranging from $n = 10^{20} cm^{-3}$ to $n = 10^{8} cm^{-3}$. Using double precision numbers in fortran could make that during computation numerical fluctuations at the bottom part of the atmosphere may "mask" or screen what happens at the top, yielding a fake production of particles out of the blue.

This method does not conserve mass.

Am I on the correct way to find an answer?

Thank you for your time and attention.

• My first guess is that you use a finite difference method for the numerical solution of your equation. Such method is not necessary mass conservative. I may write a finite volume method for your equation that will fulfill the mass conservation. To do it, one needs to know which of your data, additionally to $n$ and $T$, are depending on $z$, and, moreover, if $T$ is available only approximately, and if yes where are these approximative values available - at the same location as $n$? P.S. By the way, $n$ is missing in the integral to compute $\rho$. – Peter Frolkovič Nov 23 '16 at 15:25