Recall RHS of Eq(\ref{mass})

\begin{equation}\label{mass_R}
V_R\frac{\partial}{\partial t}\left[\phi\left(X_c\xi_oS_o+Y_c\xi_gS_g+W_c\xi_wS_w\right)\right]
\end{equation}

Let's define $A_c=\phi\left(X_c\xi_oS_o+Y_c\xi_gS_g+W_c\xi_wS_w\right)$ and discretize Eq(\ref{mass_R}) in time space we get

\begin{equation}\label{RHS}
	\frac{V_R}{\Delta t}\left[A_c^{n+1}-A_c^{n}\right]
\end{equation}

Using Taylor's series to expand $A_c^{n+1}$ as follows

\begin{equation}\label{R_taylor}
A_c^{n+1}=A_c^l+\left(\frac{\partial A_c}{\partial P_o}\right)^l\delta P_o+\left(\frac{\partial A_c}{\partial S_o}\right)^l\delta S_o+\left(\frac{\partial A_c}{\partial S_g}\right)^l\delta S_g+\sum_{i=1}^{2}\left[\left(\frac{\partial A_c}{\partial X_i}\right)^l\delta X_i+\left(\frac{\partial A_c}{\partial Y_i}\right)^l\delta Y_i\right]
\end{equation}

Let's expand Eq(\ref{R_taylor}) term by term

\begin{equation}\label{R_p}
	\begin{array}{cl}
		\left(\frac{\partial A_c}{\partial P_o}\right)^l\delta P_o & = \left(\frac{\partial}{\partial P_o}\left(\phi\left(X_c\xi_oS_o+Y_c\xi_gS_g+W_c\xi_wS_w\right)\right)\right)^l\delta P_o \\
		& = \phi\left(X_cS_o\left(c_\phi\xi_o+\frac{\partial \xi_o}{\partial P_o}\right)+Y_cS_g\left(c_\phi\xi_g+\frac{\partial \xi_g}{\partial P_o}\right)+W_cS_w\left(c_\phi\xi_w+\frac{\partial \xi_w}{\partial P_o}\right)\right)^l\delta P_o
	\end{array}
\end{equation}

\begin{equation}\label{R_so}
	\begin{array}{cl}
		\left(\frac{\partial A_c}{\partial S_o}\right)^l\delta S_o & = \left(\frac{\partial}{\partial S_o}\left(\phi\left(X_c\xi_oS_o+Y_c\xi_gS_g+W_c\xi_w(1-S_o-S_g)\right)\right)\right)^l\delta S_o \\
		& = \phi\left(X_c\xi_o-W_c\xi_w\right)^l\delta S_o
	\end{array}
\end{equation}

\begin{equation}\label{R_sg}
	\begin{array}{cl}
		\left(\frac{\partial A_c}{\partial S_g}\right)^l\delta S_g & = \left(\frac{\partial}{\partial S_g}\left(\phi\left(X_c\xi_oS_o+Y_c\xi_gS_g+W_c\xi_w(1-S_o-S_g)\right)\right)\right)^l\delta S_g \\
		& = \phi\left(Y_c\xi_g-W_c\xi_w\right)^l\delta S_g
	\end{array}
\end{equation}

For $c=1$

\begin{equation}\label{R_x1}
	\begin{array}{cl}
		\sum_{i=1}^{2}\left(\frac{\partial A_1}{\partial X_i}\right)^l\delta X_i &= \left(\frac{\partial A_1}{\partial X_1}\right)^l\delta X_1+\left(\frac{\partial A_1}{\partial X_2}\right)^l\delta X_2\\
		& = \phi\left(\frac{\partial X_1\xi_oS_o}{\partial X_1}\right)^l\delta X_1+\phi\left(\frac{\partial X_1\xi_oS_o}{\partial X_2}\right)^l\delta X_2\\
		& = \phi S_o^l\left(\left(\xi_o+X_1\frac{\partial \xi_o}{\partial X_1}\right)^l\delta X_1+\left(X_1\frac{\partial \xi_o}{\partial X_2}\right)^l\delta X_2\right)
	\end{array}
\end{equation}

\begin{equation}\label{R_y1}
	\begin{array}{cl}
		\sum_{i=1}^{2}\left(\frac{\partial A_1}{\partial Y_i}\right)^l\delta Y_i &= \left(\frac{\partial A_1}{\partial Y_1}\right)^l\delta Y_1+\left(\frac{\partial A_1}{\partial Y_2}\right)^l\delta Y_2\\
		& = \phi\left(\frac{\partial Y_1\xi_gS_g}{\partial Y_1}\right)^l\delta Y_1+\phi\left(\frac{\partial Y_1\xi_gS_g}{\partial Y_2}\right)^l\delta Y_2\\
		& = \phi S_g^l\left(\left(\xi_g+Y_1\frac{\partial \xi_g}{\partial Y_1}\right)^l\delta Y_1+\left(Y_1\frac{\partial \xi_g}{\partial Y_2}\right)^l\delta Y_2\right)
	\end{array}
\end{equation}

For $c=2$

\begin{equation}\label{R_x2}
	\begin{array}{cl}
		\sum_{i=1}^{2}\left(\frac{\partial A_2}{\partial X_i}\right)^l\delta X_i &= \left(\frac{\partial A_2}{\partial X_1}\right)^l\delta X_1+\left(\frac{\partial A_2}{\partial X_2}\right)^l\delta X_2\\
		& = \phi\left(\frac{\partial X_2\xi_oS_o}{\partial X_1}\right)^l\delta X_1+\phi\left(\frac{\partial X_2\xi_oS_o}{\partial X_2}\right)^l\delta X_2\\
		& = \phi S_o^l\left(\left(X_2\frac{\partial \xi_o}{\partial X_1}\right)^l\delta X_1+\left(\xi_o+X_2\frac{\partial \xi_o}{\partial X_2}\right)^l\delta X_2\right)
	\end{array}
\end{equation}

\begin{equation}\label{R_y2}
	\begin{array}{cl}
		\sum_{i=1}^{2}\left(\frac{\partial A_2}{\partial Y_i}\right)^l\delta Y_i &= \left(\frac{\partial A_2}{\partial Y_1}\right)^l\delta Y_1+\left(\frac{\partial A_2}{\partial Y_2}\right)^l\delta Y_2\\
		& = \phi\left(\frac{\partial Y_2\xi_gS_g}{\partial Y_1}\right)^l\delta Y_1+\phi\left(\frac{\partial Y_2\xi_gS_g}{\partial Y_2}\right)^l\delta Y_2\\
		& = \phi S_g^l\left(\left(Y_2\frac{\partial \xi_g}{\partial Y_1}\right)^l\delta Y_1+\left(\xi_g+Y_2\frac{\partial \xi_g}{\partial Y_2}\right)^l\delta Y_2\right)
	\end{array}
\end{equation}

For $c=3$

\begin{equation}\label{R_x3}
	\begin{array}{cl}
		\sum_{i=1}^{2}\left(\frac{\partial A_3}{\partial X_i}\right)^l\delta X_i &= \left(\frac{\partial A_3}{\partial X_1}\right)^l\delta X_1+\left(\frac{\partial A_3}{\partial X_2}\right)^l\delta X_2\\
		& = \phi\left(\frac{\partial (1-X_1-X_2)\xi_oS_o}{\partial X_1}\right)^l\delta X_1+\phi\left(\frac{\partial (1-X_1-X_2)\xi_oS_o}{\partial X_2}\right)^l\delta X_2\\
		& = \phi\left[S_o\left(\frac{\partial \xi_o}{\partial X_1}-\xi_o-X_1\frac{\partial \xi_o}{\partial X_1}-X_2\frac{\partial \xi_o}{\partial X_1}\right)+S_o\left(\frac{\partial \xi_o}{\partial X_2}-X_1\frac{\partial \xi_o}{\partial X_2}-\xi_o-X_2\frac{\partial \xi_o}{\partial X_2}\right)\right]\\
		& = \phi S_o^l\left[\left(\left(1-X_1-X_2\right)\frac{\partial \xi_o}{\partial X_1}-\xi_o\right)^l\delta X_1+\left(\left(1-X_1-X_2\right)\frac{\partial \xi_o}{\partial X_2}-\xi_o\right)^l\delta X_2\right]	
		\end{array}
\end{equation}

\begin{equation}\label{R_y3}
	\begin{array}{cl}
		\sum_{i=1}^{2}\left(\frac{\partial A_3}{\partial Y_i}\right)^l\delta Y_i &= \left(\frac{\partial A_3}{\partial Y_1}\right)^l\delta Y_1+\left(\frac{\partial A_3}{\partial Y_2}\right)^l\delta Y_2\\
		& = \phi\left(\frac{\partial (1-Y_1-Y_2)\xi_gS_g}{\partial Y_1}\right)^l\delta Y_1+\phi\left(\frac{\partial (1-Y_1-Y_2)\xi_gS_g}{\partial Y_2}\right)^l\delta Y_2\\
		& = \phi\left[S_g\left(\frac{\partial \xi_g}{\partial Y_1}-\xi_g-Y_1\frac{\partial \xi_g}{\partial Y_1}-Y_2\frac{\partial \xi_g}{\partial Y_1}\right)+S_g\left(\frac{\partial \xi_g}{\partial Y_2}-Y_1\frac{\partial \xi_g}{\partial Y_2}-\xi_g-Y_2\frac{\partial \xi_g}{\partial Y_2}\right)\right]\\
		& = \phi S_g^l\left[\left(\left(1-Y_1-Y_2\right)\frac{\partial \xi_g}{\partial Y_1}-\xi_g\right)^l\delta Y_1+\left(\left(1-Y_1-Y_2\right)\frac{\partial \xi_g}{\partial Y_2}-\xi_g\right)^l\delta Y_2\right]	
		\end{array}
\end{equation}

Substitute Eq (\ref{R_taylor}) - (\ref{R_y3}) into Eq (\ref{RHS})

for $c=1$

\begin{equation}\label{mass_RHS1}
	\phi\frac{V_R}{\Delta t}\left(A_1^l-A_1^n\right)+\phi\frac{V_R}{\Delta t}\left[
	\begin{array}{cc}
		\left(X_1S_o\left(c_\phi\xi_o+\frac{\partial \xi_o}{\partial P_o}\right)+Y_1S_g\left(c_\phi\xi_g+\frac{\partial \xi_g}{\partial P_o}\right)\right)^l\delta P_o\\
		+\left(X_1\xi_o\right)^l\delta S_o+\left(Y_1\xi_g\right)^l\delta S_g \\
		+S_o^l\left(\left(\xi_o+X_1\frac{\partial \xi_o}{\partial X_1}\right)^l\delta X_1+\left(X_1\frac{\partial \xi_o}{\partial X_2}\right)^l\delta X_2\right)\\
		+S_g^l\left(\left(\xi_g+Y_1\frac{\partial \xi_g}{\partial Y_1}\right)^l\delta Y_1+\left(Y_1\frac{\partial \xi_g}{\partial Y_2}\right)^l\delta Y_2\right) 
	\end{array}
	\right]
\end{equation} 

for $c=2$

\begin{equation}\label{mass_RHS2}
	\phi\frac{V_R}{\Delta t}\left(A_2^l-A_2^n\right)+\phi\frac{V_R}{\Delta t}\left[
	\begin{array}{cc}
		\left(X_2S_o\left(c_\phi\xi_o+\frac{\partial \xi_o}{\partial P_o}\right)+Y_2S_g\left(c_\phi\xi_g+\frac{\partial \xi_g}{\partial P_o}\right)\right)^l\delta P_o\\
		+\left(X_2\xi_o\right)^l\delta S_o+\left(Y_2\xi_g\right)^l\delta S_g \\
		+S_o^l\left(\left(X_2\frac{\partial \xi_o}{\partial X_1}\right)^l\delta X_1+\left(\xi_o+X_2\frac{\partial \xi_o}{\partial X_2}\right)^l\delta X_2\right)\\
		+S_g^l\left(\left(Y_2\frac{\partial \xi_g}{\partial Y_1}\right)^l\delta Y_1+\left(\xi_g+Y_2\frac{\partial \xi_g}{\partial Y_2}\right)^l\delta Y_2\right) 
	\end{array}
	\right]
\end{equation}
 
for $c=3$

\begin{equation}\label{mass_RHS3}
	\phi\frac{V_R}{\Delta t}\left(A_3^l-A_3^n\right)+\phi\frac{V_R}{\Delta t}\left[
	\begin{array}{cc}
		\left(X_3S_o\left(c_\phi\xi_o+\frac{\partial \xi_o}{\partial P_o}\right)+Y_3S_g\left(c_\phi\xi_g+\frac{\partial \xi_g}{\partial P_o}\right)\right)^l\delta P_o\\
		+\left(X_3\xi_o\right)^l\delta S_o+\left(Y_3\xi_g\right)^l\delta S_g \\
		+S_o^l\left[\left(\left(1-X_1-X_2\right)\frac{\partial \xi_o}{\partial X_1}-\xi_o\right)^l\delta X_1+\left(\left(1-X_1-X_2\right)\frac{\partial \xi_o}{\partial X_2}-\xi_o\right)^l\delta X_2\right]	\\
		+S_g^l\left[\left(\left(1-Y_1-Y_2\right)\frac{\partial \xi_g}{\partial Y_1}-\xi_g\right)^l\delta Y_1+\left(\left(1-Y_1-Y_2\right)\frac{\partial \xi_g}{\partial Y_2}-\xi_g\right)^l\delta Y_2\right]
	\end{array}
	\right]
\end{equation} 

for Water Equation

\begin{equation}\label{mass_RHSW}
	\phi\frac{V_R}{\Delta t}\left(A_c^l-A_c^n\right)+\phi\frac{V_R}{\Delta t}\left[
	\begin{array}{cc}
		\left(S_w\left(c_\phi\xi_w+\frac{\partial \xi_w}{\partial P_o}\right)\right)^l\delta P_o\\
		+\left(\xi_w\right)^l\delta S_o+\left(\xi_w\right)^l\delta S_g \\
	\end{array}
	\right]
\end{equation} 


