%\usepackage{graphicx} % Required for inserting images
%\usepackage{hyperref} 
%\usepackage[normalem]{ulem}
%\usepackage{amsmath,amssymb,bm}

%% Maths operators
\newcommand{\de}      {\mathrm{d}}
\newcommand{\De}      {\mathrm{D}}

\newcommand{\grad}    {\bm{\nabla}}
\newcommand{\Div}     {\bm{\nabla}\cdot}

%% Vectors
\newcommand{\uv}      {\bm{u}}

\section*{}
The current implementation of correction to the energy equation in the presence of mass diffusion follows \href{https://doi.org/10.1080/03091929.2019.1634705}{ Gent et al. 2020, GAFD, 114, 77}. The basic approach followed in the article is to treat $$ \frac{\De }{\De t}(\rho e)-e\frac{\De\rho}{\De t}=\text{RHS} $$ (see below) as invariant under the correction for mass diffusion. However, it is also possible to follow a different line of reasoning to arrive at the required correction for the energy equation in the presence of mass diffusion. Here, we present such an alternate approach that leads to modification of the correction terms in the temperature equation, while the entropy equation remains unchaged from that with no correction.

We know that the first law of thermodynamics can be expressed mathematically for a fluid medium as follows:
\begin{displaymath} 
   \frac{\De q}{\De t}
   = \frac{\De e}{\De t}+p\frac{\De} {\De t} \left(\frac{1} {\rho}\right) 
\end{displaymath}
or 
\begin{displaymath} 
  \rho \frac{\De q}{\De t}
   = \rho \frac{\De e}{\De t}-\frac{p}{\rho}\frac{\De\rho}{\De t} 
\end{displaymath}
or 
\begin{equation} \label{eq:first-law}
  \rho T \frac{\De s}{\De t}
   = \rho \frac{\De e}{\De t}-\frac{p}{\rho}\frac{\De\rho}{\De t} ,
\end{equation} 
where $q$ denotes the heat gained per unit mass, $t$ denotes time, $e$ is the internal energy per unit mass, $p$ is the pressure, $\rho$ is the density, $T$ is the temperature, and $s$ is the entropy per unit mass. As usual, $\De/\De t \equiv \partial/\partial t + \uv\cdot\grad$ is the advective derivative, where $\uv$ is the fluid velocity.

Based on the above, we shall use the following two equivalent forms of the first law (equation \ref{eq:first-law}) for our subsequent discussion:
\begin{equation} \label{eq:entropy}
  \rho T\frac{\De s}{\De t}
   = \mathcal{L}
\end{equation}
or 
\begin{equation} \label{eq:internal-energy}
  \rho \frac{\De e}{\De t}-\frac{p}{\rho}\frac{\De\rho}{\De t}
   = \mathcal{L}.
\end{equation}
Here, $\mathcal{L}$ represents the rate of heating per unit volume, which is the net effect of all the sources and sinks of heat for the system. 

Furthermore, the continuity equation in the absence of any mass diffusion looks like:

\begin{equation}\label{eq:continuity}
  \frac{\De\rho}{\De t}
  = - \rho \Div\uv.
\end{equation}
However, in the presence of mass diffusion, this gets modified to
\begin{equation}\label{eq:continuity-diffusion}
  \frac{\De\rho}{\De t}
  = - \rho \Div\uv + \zeta_D\nabla^2\rho\; ,
\end{equation}
where $\zeta_D$ is the mass diffusion coefficient.

The first law of thermodynamics in the form expressed as Eq.~\ref{eq:entropy} or \ref{eq:internal-energy} is independent of whether mass diffusion is present (equation \ref{eq:continuity-diffusion}) or not (equation \ref{eq:continuity}). 

There are basically two forms of the energy equation implemented in the Pencil Code, one in terms of $s$ and the other in terms of $\ln T$. Let's look at them one by one.


\begin{enumerate}
\item 
\textbf{Entropy equation:}
The Pencil Code solves the equation for specific entropy $s$ in the form:
\begin{equation}\label{eq:entropy-code}
  \frac{\partial s}{\partial t} = -\uv\cdot\grad s
  +\frac{\mathcal{L}}{\rho T },
\end{equation}
which can be seen to be directly equivalent to equation (\ref{eq:entropy}) after rearrangement. This is true whether the continuity equation has the form expressed in equation (\ref{eq:continuity}) or that in equation (\ref{eq:continuity-diffusion}).


\item 
\textbf{Temperature equation:}
The other form of the energy equation, expressed in terms of the temperature, is:
\begin{equation}\label{eq:temperature}
\frac{\partial\ln T}{\partial t} = -\uv\cdot\grad \ln T
  +\frac{\mathcal{L}}{\rho c_v T} -\left(\gamma-1\right) \Div \uv.
\end{equation}
Here, $c_v$ is the specific heat at constant volume, and $\gamma=c_p/c_v$ is the ratio of specific heats. Let us now examine whether equation (\ref{eq:temperature}) is equivalent to either of the  equations \eqref{eq:entropy} or \eqref{eq:internal-energy}. After slight rearrangement, equation \eqref{eq:temperature} becomes
\begin{displaymath}
\rho c_v T \frac{\De \ln T}{\De t} +\left(\gamma-1\right) \rho c_v T \Div \uv = 
  \mathcal{L} 
\end{displaymath}
or 
\begin{displaymath}
\rho c_v  \frac{\De T}{\De t} +\left(\gamma-1\right) \rho c_v T \Div \uv = 
  \mathcal{L} 
\end{displaymath}
or 
\begin{equation} \label{eq:energy-nodiffusion}
\rho  \frac{\De e}{\De t} +p \Div \uv = 
  \mathcal{L} ,
\end{equation}
where we have used $\De e=c_v \De T$ and $p=\left(\gamma-1\right) \rho c_v T$ for an ideal gas.

When mass diffusion is absent, we can substitute for $\Div \uv$ from equation \eqref{eq:continuity} into equation \eqref{eq:energy-nodiffusion} to recover equation \eqref{eq:internal-energy}. However, this is not possible in the presence of mass diffusion (equation \ref{eq:continuity-diffusion}).
Therefore, we need to include a correction term in equation \eqref{eq:temperature} in that case. If we add the term $[ \left(\gamma-1\right) / \rho]\zeta_D\nabla^2\rho\ $ to the RHS of equation \eqref{eq:temperature} when mass diffusion is present, equation \eqref{eq:energy-nodiffusion} gets modified to 
\begin{equation} \label{eq:energy-diffusion}
\rho  \frac{\De e}{\De t} +p \Div \uv - \frac{p}{\rho} \zeta_D \nabla^2\rho= 
  \mathcal{L} ,
\end{equation}
by virtue of $p=\left(\gamma-1\right) \rho c_v T$. 
Now, when we combine equation \eqref{eq:energy-diffusion} with equation \eqref{eq:continuity-diffusion}, we can recover equation \eqref{eq:internal-energy} representing the conservation of energy. 

Therefore, we require a correction term $[ \left(\gamma-1\right) / \rho]\zeta_D\nabla^2\rho\ $ to be added to the RHS of equation \eqref{eq:temperature} in the presence of mass diffusion, as under:    %MR: under?
\begin{align}
\frac{\partial\ln T}{\partial t}=&-\uv\cdot\grad \ln T
  +\frac{\mathcal{L}}{\rho c_v T} -\left(\gamma-1\right) \Div \uv
\nonumber \\
&+ \frac{\left(\gamma-1\right)} {\rho} \zeta_D\nabla^2\rho.
\end{align}

\end{enumerate}

\noindent We note that the basic difference between the above approach and the one used in Gent et al. (2020) is that instead of equation \eqref{eq:internal-energy}, the equation 
\begin{equation} \label{eq:gent20}
 \frac{\De }{\De t}(\rho e)-e\frac{\De\rho}{\De t} =\text{RHS},
\end{equation}
where $\text{RHS}=\mathcal{L}+\frac{p}{\rho}\frac{\De\rho}{\De t}$, is used as an invariant under the correction for mass diffusion, and the correction is meant to be suitably applied to both forms of the energy equation \eqref{eq:entropy-code} and \eqref{eq:temperature}. Note that in Gent et al (2020), the term $\frac{p}{\rho}\frac{\De\rho}{\De t}$ on the RHS is not taken into account when treating equation \eqref{eq:gent20} as invariant. The correction in their paper turns out to be subtracting the term $( 1 / \rho)\zeta_D\nabla^2\rho $ from the RHS of equation \eqref{eq:temperature} and the term $( c_v / \rho)\zeta_D\nabla^2\rho $ from the RHS of equation \eqref{eq:entropy-code}, as follows:
\begin{equation}
\frac{\partial\ln T}{\partial t} = -\uv\cdot\grad \ln T
  +\frac{\mathcal{L}}{\rho c_v T} -\left(\gamma-1\right) \Div \uv- \frac{1} {\rho} \zeta_D\nabla^2\rho\,; 
\end{equation}
\begin{equation}
  \frac{\partial s}{\partial t} = -\uv\cdot\grad s
  +\frac{\mathcal{L}}{\rho T }-\frac{c_v} {\rho}\zeta_D\nabla^2\rho.
\end{equation}

In future, it will be useful to consider the effect of implementing the mass diffusion correction in different ways in the code.

