\documentclass[11pt,letterpaper]{article}

\usepackage[pdftex]{graphicx}

\usepackage{natbib}

\usepackage{fullpage}

\usepackage{lineno}

\usepackage{multirow}

\usepackage{wrapfig}

\usepackage{amsmath}

\usepackage{amssymb}

\usepackage{sidecap}

\begin{document}

\setlength{\parindent}{0mm}

\setlength{\parskip}{0.4cm}

\bibliographystyle{apalike}

%\modulolinenumbers[5]

%\linenumbers

\title{\textbf{GEMINI}: \textbf{G}eospace \textbf{E}nvironment \textbf{M}odel of \textbf{I}on-\textbf{N}eutral \textbf{I}nteractions}

\author{Matthew D. Zettergren, PhD\\ Associate Professor of Engineering Physics\\ Center for Space and Atmospheric Physics\\ Physical Sciences Department \\Embry-Riddle Aeronautical University\\zettergm@erau.edu\\mattzett@gmail.com}

\maketitle

\tableofcontents

\pagebreak

%\section{TO DO LIST FOR THIS DOCUMENT:}

%Check section on artificial viscosity

%Check section with advection algorithm

%Check section with geomag. to geographic transformations (they differ from \citet{Huba:2000})

%Add advection boundary conditions section

\section{Preface}

\subsection{Purpose of this document}

The \textbf{G}eospace \textbf{E}nvironment \textbf{M}odel of \textbf{I}on-\textbf{N}eutral \textbf{I}nteractions (GEMINI) is a general-purpose, three-dimensional (3D) terrestrial ionospheric model capable of describing most processes relevant to the ionosphere at medium to small spatial scales (200 m to 10000 km). In the field of ionospheric physics, there are several well-established global models of the terrestrial ionosphere exist; however, there are comparatively fewer ``standard'' models for small-scale phenomena (as of 2018), which was a motivating factor behind developing GEMINI and making it available to our scientific community.

This document outlines the mathematical formulation of GEMINI and is intended to serve as a bridge between published descriptions of the model \citep{Zettergren:2012,Zettergren:2014,Zettergren:2015,Zettergren:2015b}, which are necessarily fairly terse, and the mathematical and (some) implementation details necessary to understand how the model functions. A separate README describing how to compile and run the model is included on the GitHub site for the GEMINI project, which currently houses some 8 branches of the code, all adapted to different specific problems of interest but each having the same core numerical components.

\subsection{Error reporting}

Please report any apparent errors to me at: mattzett@gmail.com

\subsection{Design philosophy}

During design of GEMINI top-level clarity and organization have generally been favored over speed and efficiency. This decision was largely driven by hope that GEMINI would prove useful to others by being easy to adapt to a wide variety of problems; hence, the model has been designed with source code modification in mind. Undoubtedly optimizations can and will be made in future branches/forks as needs arise.

GEMINI has also been designed to be easily extensible; e.g. it is feasible to add additional ion species into the model or to use a Cartesian or other type of orthogonal curvilinear grid. In principle, the GEMINI codebase could potentially be used to develop simulations of ionospheres (more generally, collisional plasmas) for many different types of of planets/environments.

As a final, general, note, GEMINI can run with either a 2D or 3D orthogonal, curvilinear grid. The 2D functionality was included to facilitate very high-resolution runs in two dimensions, which become overly burdensome in GEMINI's parent MATLAB code. Special cases in the mathematical formulation and numerical implementation arising when using a 2D grid are discussed in detail in Section \ref{sec:2D} of this document.

\subsection{GEMINI development history}

GEMINI is an evolution of the model first used in \citet{Zettergren:2012} to study terrestrial ionospheric plasma density structures generated by two-dimensional auroral current systems. It was later extended by \citet{Zettergren:2013} to study low-latitude forcing by infrasonic/acoustic wave and then in \citet{Zettergren:2014} to ingest highly processed data products to specify boundary conditions. These initial model prototypes were coded in MATLAB; \citet{Zettergren:2015b} translated the model into fortran with MPI parallelization and added in the third dimension in Cartesian coordinates. The present version extends the 3D fortran model further to include general orthogonal curvilinear coordinates (publication TBD).

\subsection{Executive summary}

The 3D ionospheric ``\underline{G}eospace \underline{E}nvironment \underline{M}odel for \underline{I}on-\underline{N}eutral \underline{I}nteractions'' (GEMINI), is based on the existing 2D model developed by \citet{Zettergren:2012} and later extended by \citet{Zettergren:2013,Zettergren:2014,Zettergren:2015b}. The present form of GEMINI functions in two or three dimensions, and uses general orthogonal curvilinear coordinates, usually either a tilted dipole \citep{Huba:2000} or Cartesian system. The model comprises a fluid system of equations \citep{Schunk:1977,Blelly:1993}, describing dynamics of the ionospheric plasma, self-consistently coupled to an electrostatic treatment of auroral and neutral dynamo currents. The fluid system is a set of three conservation equations (mass, momentum, and energy) for each ionospheric species $s$ relevant to the E-, F-, and topside regions ($s=\mathrm{O^+,NO^+,N_2^+,O_2^+,N^+,H^+}$).

\begin{linenomath*} \begin{equation}

\frac{\partial \rho_s}{\partial t} + \nabla \cdot \left( \rho_s \mathbf{v}_s \right) = m_s P_s - L_s \rho_s \label{continuity}

\end{equation} \end{linenomath*}

\begin{linenomath*} \begin{equation}

\left[ \frac{\partial }{\partial t} \left( \rho_s \mathbf{v}_s \right) + \nabla \cdot \left( \rho_s \mathbf{v}_s \mathbf{v}_s \right) \right] \cdot \hat{\mathbf{e}}_1 = \left[ -\nabla p_s + \rho_s \mathbf{g} + \frac{\rho_s} {m_s} q_s \mathbf{E} + \sum_t \rho_s \nu_{st} \left(\mathbf{v}_t - \mathbf{v}_s \right) \right] \cdot \hat{\mathbf{e}}_1 \label{momentum}

\end{equation} \end{linenomath*}

\begin{linenomath*} \begin{eqnarray}

\frac{\partial}{\partial t} \left( \rho_s \epsilon_s \right) + \nabla \cdot \left( \rho_s \epsilon_s \mathbf{v}_s \right) &=& - p_s (\nabla \cdot \mathbf{v}_s) - \nabla \cdot \mathbf{h}_s - \\

\nonumber&& \frac{1}{(\gamma_s - 1 )}\sum_t \frac{\rho_s k_B \nu_{st}}{m_s + m_t}\left[ 2 (T_s-T_t) - \frac{2}{3} \frac{m_t}{k_B}(\mathbf{v}_s-\mathbf{v}_t)^2 \right] \label{energy}

\end{eqnarray} \end{linenomath*}

The terms on the right-hand side of Equation \ref{continuity} encapsulate chemical production and impact ionization ($P_s$) and chemical loss ($L_s$). Source terms in the continuity equation for photoionization are calculated according the method presented in \citet{Solomon:2005} using solar fluxes from the EUVAC model \citep{Richards:1994}. Impact ionization is computed using the semi-empirical method of \citet[][and references therein]{Fang:2008}. Chemical reactions for the ionospheric model are taken from \citet[][and references therein]{Diloy:1996,StMaurice:1998}. In Equation \ref{momentum} $q_s$ is the charge of each species, $\mathbf{v}_t$ is the drift velocity of species $t$ (which can be either charged or neutral), and $\nu_{st}$ is the collision frequency of charged species $s$ with species $t$. Note that the momentum equation (\ref{momentum}) is solved in a time-dependent form only for the direction parallel to the geomagnetic field (denoted by the unit vector $\hat{\mathbf{e}}_1$). The partial pressure $p_s$ is related to the specific internal energy and temperature by an equation of state (i.e. the ideal gas law). Heat fluxes in Equation \ref{energy} ($\mathbf{h}_s$) are specified by a simple model of thermal conduction for the ions:

\begin{linenomath*} \begin{equation}

\mathbf{h}_s = - \lambda_s \nabla_{\parallel} T_s,

\end{equation} \end{linenomath*}

where $\lambda_s$ is the thermal conductivity, taken from \citet{Schunk:1974}. %Physical ion stresses are neglected in our model.

For the dimensions perpendicular to the geomagnetic field line, a steady-state momentum approximation is used:

\begin{linenomath*} \begin{equation}

\mathbf{v}_{s \perp} = \boldsymbol{\mu}_{s \perp} \cdot \left( \mathbf{E}_\perp + \frac{m_s \nu_s}{q_s} \mathbf{v}_{n\perp} \right). \label{drifts}

\end{equation} \end{linenomath*}

In this expression, $\boldsymbol{\mu}_{s\perp}$ is the ion mobility tensor and $\nu_s$ is the total ion-neutral collision frequency \citep{Zettergren:2012}. A static geomagnetic field is used to compute the ion mobilities in equation \ref{drifts}. %Note that this steady-state perpendicular drift assumption is consistent with the use of an electrostatic description for the fields, \citep[c.f.][ for a complete discussion]{Zettergren:2012}.

Mass and momentum density state variables for the electron species are solved by invoking quasi-neutrality and the definition of current density \citep[c.f.][equations 25 and 26]{Zettergren:2012}. A full transport equation is solved for the electron energy.

\begin{linenomath*} \begin{eqnarray}

\frac{\partial}{\partial t} \left( \rho_e \epsilon_e \right) + \nabla \cdot \left( \rho_e \epsilon_e \mathbf{v}_e \right) &=& - p_e (\nabla \cdot \mathbf{v}_e) - \nabla \cdot \mathbf{h}_e - \\

\nonumber&& \frac{1}{(\gamma_e - 1 )}\sum_t \frac{\rho_e k_B \nu_{et}}{m_e + m_t}\left[ 2 (T_e-T_t) - \frac{2}{3} \frac{m_t}{k_B}(\mathbf{v}_e-\mathbf{v}_t)^2 \right] + \frac{Q_e}{(\gamma_e - 1 )}\label{energye}

\end{eqnarray} \end{linenomath*}

For electron heat fluxes, both thermoelectric effects \citep{Schunk:1978} and thermal conduction are considered \citep{Banks:1973,Huba:2000}

\begin{linenomath*} \begin{equation}

\mathbf{h}_e = \left( - \lambda_e \nabla_{\parallel} T_e - \beta_e \mathbf{J}_{\parallel} \right).

\end{equation} \end{linenomath*}

The electron energy equation also differs in form from the ion equation above (equation \ref{energy}) in that it includes inelastic cooling terms and heating by photoelectrons, collectively denoted by $Q_e$. For the present work, cooling due to the excitation of rotational and vibrational modes of O$_2$ and N$_2$ are included, as is the excitation of fine structure of O. The photoelectron heating rate is calculated from the photoionization rate according to the method of \citet{Swartz:1972}.

Perpendicular components of the momentum equations are solved with a steady-state force balance approximation, as described in \citet{Zettergren:2012}. Electric fields are found by enforcing a divergence free current density, where the current is assumed to consist of both conduction currents (directly dependent on the electric field through the usual ionospheric Ohm's law) and polarization currents (from a spatio-temporally varying electric field) \citep[c.f.][ and references therein]{Mitchell:1985,Gondarenko:1999}. In essence, this approach assumes that the behavior of the ionospheric response currents is, to leading order, electrostatic. These currents, however, may be modified by a correction which accounts for polarization currents from slowly varying field structures. Hence, the ionospheric electric field is approximated by $\mathbf{E}=-\nabla \Phi$. Further assuming that the geomagnetic field lines are equipotentials and employing a field line integration yields an equation that can be solved for the electric potential, hence field:

\begin{equation}

\nabla_\perp \cdot \left( \boldsymbol{\Sigma}_\perp \cdot \nabla_\perp \Phi \right) + \nabla_\perp \cdot \left[ C_M \left( \frac{\partial}{\partial t} + \mathbf{v}_\perp \cdot \nabla_\perp \right) \left( \nabla_\perp \Phi \right) \right] = \nabla_\perp \cdot \left( \boldsymbol{\Sigma}_\perp \cdot \mathbf{E}_{0\perp} \right) \quad \label{eqn:potential}

\end{equation}

In this equation $\Phi$ is the electric potential generated by the ionospheric response to a magnetospherically imposed electric field $\mathbf{E}_{0\perp}$ and $\mathbf{v}_\perp = \left( \mathbf{E}_{0\perp} -\nabla \Phi \right) \times \mathbf{B}/B^2$. $\boldsymbol{\Sigma}_\perp$ is the field-line integrated conductance tensor, and $C_M$ is the inertial capacitance defined, e.g., in \citet{Mitchell:1985}. Note that Equation \ref{eqn:potential} implicitly enforces the condition that all current closes within the ionosphere through either conduction currents or polarization currents (viz. no field-aligned current flow through the top boundary of the model). Past work has generally found that polarization currents suppress GDI growth through nonlocal closure of currents, similar in effect to shorting out of GDI through a strongly conducting E-region. However, this term is also responsible for destabilizing the KHI mode \citep{Keskinen:1988}, which alternatively may be stabilized by conduction currents. The potential defined by Equation \ref{eqn:potential} does not include the ambipolar contribution; instead the ambipolar potential is computed during a separate calculation that resolves the parallel momentum equation (the only equation where the ambipolar term is needed). This ambipolar field evaluted, during this step, directly from the electron pressure.

\begin{linenomath*} \begin{equation}

\mathbf{E}_{a,\parallel} = \frac{1}{n_e q_e} \nabla_\parallel p_e \label{ambipolar}

\end{equation} \end{linenomath*}

This ambipolar field is added to the ``resistive'' part of the parallel field defined by $\mathbf{E}_{r,\parallel} = - \nabla_\parallel \Phi$ to form the total parallel electric field used in the ion momentum equations (c.f. Equation \ref{momentum}).

The most recent detailed descriptions of GEMINI are in recent publications by \citet[][Appendix A]{Zettergren:2015} and \citet{Zettergren:2015b}, though the model has been around for a while in various forms \citep[e.g.][]{Zettergren:2012}.

GEMINI includes capabilities to couple input from neutral dynamics models by way of files and automatically does interpolation, and rotation of neutral wind, density, and temperature perturbations (deviations from MSIS00) onto its curvilinear grid. These perturbations are then factored into all ionospheric calculations that use neutral parameters (dynamo, collisions, ionization, chemistry, etc.). Note that the coupling here is one-way (viz. neutrals affecting ions). Examples of using neutral dynamics information to drive GEMINI have been published by \citet{Zettergren:2013,Zettergren:2015,Zettergren:2017}.

Boundary conditions for electric potential, background electric fields, and field-aligned currents and electron precipitation can be included in GEMINI through input files. These inputs can be specified independent on a grid different from that used in GEMINI and are automatically interpolated to the simulation grid positions.

A branch of GEMINI exists that uses the suprathermal electron transport code GLOW to specify ionization and heating reates instead of \citet{Fang:2008,Swartz:1972}. This version has the ability to separate between thermal and suprathermal field-aligned currents and is able to more accurately specify both thermal electron heating by suprathermal electrons and ionization rates of major species for an arbitration input electron precipitation pattern and distribution.

\section{Mathematical formulation}

GEMINI comprises a fluid system of equations \citep{Schunk:1977,Blelly:1993,Huba:2000}, describing dynamics of the ionospheric plasma, self-consistently coupled to an electrostatic treatment of auroral and neutral dynamo currents that is based on a steady-state current continuity equation. These systems of equations are discussed in detail in this section.

\subsection{Fluid system of equations}

The fluid system is a set of three conservation equations (mass, momentum, and energy) for each ionospheric species $s$ relevant to the E-, F-, and topside regions ($s=\mathrm{O^+,NO^+,N_2^+,O_2^+,N^+,H^+}$):

\begin{linenomath*} \begin{equation}

\frac{\partial \rho_s}{\partial t} + \nabla \cdot \left( \rho_s \mathbf{v}_s \right) = m_s P_s - L_s \rho_s \label{continuity}

\end{equation} \end{linenomath*}

\begin{linenomath*} \begin{equation}

% \left[ \frac{\partial }{\partial t} \left( \rho_s \mathbf{v}_s \right) + \nabla \cdot \left( \rho_s \mathbf{v}_s \mathbf{v}_s \right) \right] \cdot \hat{\mathbf{e}}_1 = -\nabla p_s + \rho_s \mathbf{g} + \frac{\rho_s} {m_s} q_s \left( \mathbf{E} + \mathbf{v}_s \times \mathbf{B} \right) + \sum_t \rho_s \nu_{st} \left(\mathbf{v}_t - \mathbf{v}_s \right)\label{momentum}

\left[ \frac{\partial }{\partial t} \left( \rho_s \mathbf{v}_s \right) + \nabla \cdot \left( \rho_s \mathbf{v}_s \mathbf{v}_s \right) \right] \cdot \hat{\mathbf{e}}_1 = \left[ -\nabla p_s + \rho_s \mathbf{g} + \frac{\rho_s} {m_s} q_s \mathbf{E} + \sum_t \rho_s \nu_{st} \left(\mathbf{v}_t - \mathbf{v}_s \right) \right] \cdot \hat{\mathbf{e}}_1 \label{momentum}

\end{equation} \end{linenomath*}

\begin{linenomath*} \begin{eqnarray}

\frac{\partial}{\partial t} \left( \rho_s \epsilon_s \right) + \nabla \cdot \left( \rho_s \epsilon_s \mathbf{v}_s \right) &=& - p_s (\nabla \cdot \mathbf{v}_s) - \nabla \cdot \mathbf{h}_s - \\

\nonumber&& \frac{1}{(\gamma_s - 1 )}\sum_t \frac{\rho_s k_B \nu_{st}}{m_s + m_t}\left[ 2 (T_s-T_t) - \frac{2}{3} \frac{m_t}{k_B}(\mathbf{v}_s-\mathbf{v}_t)^2 \right] \label{energy}

\end{eqnarray} \end{linenomath*}

The terms on the right-hand side of Equation \ref{continuity} encapsulate chemical production and impact ionization ($P_s$) and chemical loss ($L_s$). Source terms in the continuity equation for photoionization are calculated according the method presented in \citet{Solomon:2005} using solar fluxes from the EUVAC model \citep{Richards:1994}. Impact ionization is computed by the semi-empircal method described in \citet{Fang:2008}. Chemical reactions for the ionospheric model are taken from \citet[][and references therein]{Diloy:1996,Stmaurice:1998}. In Equation \ref{momentum} $q_s$ is the charge of each species, $\mathbf{v}_t$ is the drift velocity of species $t$ (which can be either charged or neutral), and $\nu_{st}$ is the collision frequency of charged species $s$ with species $t$. Note that the momentum equation (\ref{momentum}) is solved in a time-dependent form only for the direction parallel to the geomagnetic field (denoted by the unit vector $\hat{\mathbf{e}}_1$). The partial pressure $p_s$ is related to the specific internal energy by an equation of state (for each species):

\begin{equation}

\rho_s \epsilon_s = \frac{p_s}{\gamma_s - 1} = \frac{n_s k_B T_s}{\gamma_s - 1}.

\end{equation}

$\gamma_s$ is the adiabatic index for species $s$ and $k_B$ is the Boltzmann constant. Heat fluxes in Equation \ref{energy} ($\mathbf{h}_s$) are specified by a simple model of thermal conduction for the ions:

\begin{linenomath*} \begin{equation}

\mathbf{h}_s = - \lambda_s \nabla_{\parallel} T_s,

\end{equation} \end{linenomath*}

where $\lambda_s$ is the thermal conducivity, taken from \citet{Schunk:1975}. Note that, as implied in the fluid equations above, physical ion stresses are neglected in our model. A steady-state momentum balance is used for the perpendicular ion drifts, as described in Section \ref{EM}.

%Ion momentum in the perpendicular direction uses a steady state approximation:

%\begin{linenomath*} \begin{equation}

%\mathbf{v}_{s \perp} = \boldsymbol{\mu}_{s \perp} \cdot \left( \mathbf{E}_\perp + \frac{m_s \nu_s}{q_s} \mathbf{v}_{n\perp} \right). \label{drifts}

%\end{equation} \end{linenomath*}

%In this expression, $\boldsymbol{\mu}_{s\perp}$ is the ion mobility tensor and $\nu_s$ is the total ion-neutral collision frequency \citep{Zettergren:2012}. A static geomagnetic field is used to compute the ion mobilities in equation \ref{drifts}. Note that this steady-state perpendicular drift assumption is consistent with the use of an electrostatic description for the fields. Note that time-dependence and ion inertia are ignored in Equation \ref{drifts}, although an exception to this is sometimes made in the electrodynamics solutions to accommodate small-scales and fast-varying systems- cf. Section \ref{sec:electrodynamics})

Mass and momentum density state variables for the electron species are solved by invoking quasi-neutrality and the definition of current density \citep[cf. also][equations 25 and 26]{Zettergren:2012}:

\begin{eqnarray}

n_e &=& \sum_{s \ne e} n_s \label{eqn:quasineu} \\

\mathbf{v}_e &=& \mathbf{J} - \frac{1}{n_e q_e} \sum_{s \ne e} n_s \mathbf{v}_s

\end{eqnarray}

A full transport equation is must be solved for the electron energy to include the effects of photoionization heating and heating by suprathermal electron precipitation.

\begin{linenomath*} \begin{eqnarray}

\frac{\partial}{\partial t} \left( \rho_e \epsilon_e \right) + \nabla \cdot \left( \rho_e \epsilon_e \mathbf{v}_e \right) &=& - p_e (\nabla \cdot \mathbf{v}_e) - \nabla \cdot \mathbf{h}_e - \\

\nonumber&& \frac{1}{(\gamma_e - 1 )}\sum_t \frac{\rho_e k_B \nu_{et}}{m_e + m_t}\left[ 2 (T_e-T_t) - \frac{2}{3} \frac{m_t}{k_B}(\mathbf{v}_e-\mathbf{v}_t)^2 \right] + \frac{Q_e}{(\gamma_e - 1 )}\label{energye}

\end{eqnarray} \end{linenomath*}

For electron heat fluxes, both thermoelectric effects \citep{Schunk:1978} and thermal conduction are considered \citep{Banks:1973,Huba:2000}

\begin{linenomath*} \begin{equation}

\mathbf{h}_e = \left( - \lambda_e \nabla_{\parallel} T_e - \beta_e \mathbf{J}_{\parallel} \right).

\end{equation} \end{linenomath*}

The electron energy equation also differs in form from the ion equation above (equation \ref{energy}) in the inclusion of inelastic cooling terms and heating by photoelectrons, collectively denoted by $Q_e$. For the present work, cooling due to the excitation of rotational and vibrational modes of O$_2$ and N$_2$ are included, as is the excitation of fine structure of O. The photoelectron heating rate is calculated according to the method presented in \citet{Swartz:1972} using photoionization rates computed in the model.

%Strictly speaking, these equations must be supplemented by Maxwell's equations for the electromagnetic fields. However, an electostatic assumption (or quasi-electrostatic assumptions) holds to a reasonable degree of accuracy, so a simpler formulation described in Section \ref{EM}, below, can be used in GEMINI.

\subsection{Electrostatic and quasi-electrostatic approximations} \label{EM}

For a coupled description of transport and electrodynamics in auroras, the fluid equations presented above must be supplemented by Maxwell's equations. However, the issue of finding the ionospheric electric field that results from a given magnetospheric boundary condition may be treated, to some degree of accuracy, as an electrostatic problem. The details of this approximation have been covered by several authors \citep[e.g.][and references therein]{Stmaurice:1996}, and the net effect is that the displacement current and electromagnetic induction may be neglected in situations where wave dynamics are not important or have negligible effects. This approximation precludes modeling of Alfv\'en waves and related inductive M-I coupling effects as discussed, for example in \citet{Lotko:2004}. Our focus in the following sections on slowly varying current systems supports this use of the electrostatic approximation. %These constraints reduce the Maxwell equations to the relations:

%\begin{linenomath*} \begin{eqnarray}

%\nabla^2 \Phi = - \frac{\rho_c}{\epsilon_0} &=& - \frac{\sum_s n_s q_s}{\epsilon_0} \\

%\nabla \cdot \mathbf{J} = \nabla \cdot \left( \sum_s n_s q_s \mathbf{v}_s \right) &=& 0 \label{divJSS}

%\end{eqnarray} \end{linenomath*}

Qausineutrality holds very strongly in the plasma (Equation \ref{eqn:quasineu}), and a steady-state momentum balance holds under certain conditions, which allows the velocities to be expressed in a simple way as a function of the electric field (e.g. as noted above in Equation \ref{drifts}). These facts facilitate the use of a steady state current continuity equation,

\begin{linenomath*} \begin{equation}

\nabla \cdot \mathbf{J} = \nabla \cdot \left( \sum_s n_s q_s \mathbf{v}_s \right) = 0, \label{divJSS}

\end{equation} \end{linenomath*}

to directly solve for electric potential as discussed in detail below.

%Several comments are in order regarding the connection of these equations to the fluid formulation of Equations \ref{continuity}, \ref{momentum}, and \ref{energy}. First, the current continuity equation is linearly dependent on the mass continuity equations (including all ions and electron). Because of this either all of the mass continuity and momentum equations must be solved, or all but one can be solved provided that the current continuity equation is added in to close the system. The former choice necessitates the use of Poisson's equation for specifying electric potential, hence electric field. The latter choice, however, affords a useful shortcut since the plasma is quasi-neutral ($n_e \approx \sum_{s \ne e} n_s$) and the drift velocities can be expressed simply as a function of electric field. In this case Equation \ref{divJSS} can directly be used to solve for electric potential, as discussed below.

\subsection{Momentum balance approximation and Ohm's law}

Ordinarily one would need to solve a full time dependent equation for ion momenrum perpendicular to the geomagnetic field, similar to Equation \ref{momentum}. However, in GEMINI, a steady state momentum balance assumption is used and admits a much simpler formulation of ionospheric electric fields. If the time-derivative and inertial terms of the perpendicular momentum equation are neglected, along with gravitational forces, then the drifts can be solved for explicitly in terms of the force per unit charge:

\begin{linenomath*} \begin{eqnarray}

\mathbf{v}_{s\perp} &=& \boldsymbol{\mu}_{s\perp} \cdot \left( \mathbf{E}_\perp - \frac{1}{n_s q_s} \nabla_\perp p_s + \frac{m_s \nu_s}{q_s} \mathbf{v}_{n\perp} \right) \label{mombalance} \\

\mathbf{v}_{e\parallel} &=& \mu_{e0} \left( \mathbf{E}_\parallel - \frac{1}{n_e q_e} \nabla_\parallel p_e \right) \label{mombalancepar}

\end{eqnarray} \end{linenomath*}

In this expression, $\boldsymbol{\mu}_{s\perp}$ is the ion mobility tensor and $\nu_s$ is the total ion-neutral collision frequency \citep{Zettergren:2012}, and $\mathbf{v}_n$ is the neutral wind velocity. A static geomagnetic field is used to compute the ion mobilities in Equation \ref{mombalance}. It is again noted that time-dependence and ion inertia are ignored in Equation \ref{mombalance}, although an exception to this is sometimes made in the electrodynamics solutions to accommodate small-scales and fast-varying systems- cf. Section \ref{sec:electrodynamics}). The ion and electron mobilities are defined using orthogonal coordinates with $\mathbf{B}$ in the $x_1-$direction by:

\begin{linenomath*} \begin{eqnarray}

\boldsymbol{\mu}_{s\perp}

\equiv

\frac{q_s}{m_s \nu_s} \left[ \begin{array}{cc}

\frac{\nu_{s}^2}{\nu_{s}^2 + \Omega_s^2} & \frac{\nu_{s} \Omega_s}{\nu_{s}^2+\Omega_s^2} \\

- \frac{\nu_{s} \Omega_s}{\nu_{s}^2+\Omega_s^2} & \frac{\nu_{s}^2}{\nu_{s}^2 + \Omega_s^2}

\end{array} \right]

&=&

\left[ \begin{array}{cc}

\mu_{sP} & -\mu_{sH} \\

\mu_{sH} & \mu_{sP}

\end{array} \right]

\label{mus} \\

\mu_{e0} & \equiv & \frac{q_e}{m_e \nu_e'}

\end{eqnarray} \end{linenomath*}

The cyclotron frequency is $\Omega_s \equiv q_s B_z/m_s$, and $\nu_s$ is the \emph{total} collision frequency. For the mobilities $\mu_{sP}$ and $\mu_{sH}$, the total collision frequency is calculated as a sum over ion-neutral collisions, which are dominant in regions where perpendicular transport is important.

\begin{linenomath*} \begin{equation}

\nu_s = \sum_n \nu_{sn}

\end{equation} \end{linenomath*}

Both electron-neutral and electron-ion collisions are used in computing the parallel mobility of the electrons, which dominates the parallel conductivity. Coulomb collisions are necessary since the field-aligned currents carried by electrons flow at high altitudes.

\begin{linenomath*} \begin{equation}

\nu_e' = \sum_n \nu_{en} + \sum_j \nu_{ej}

\end{equation} \end{linenomath*}

In this equation the quantity $\nu_{ej}$ refers to collisions of electrons with charged species $j$. As a final note, neutral and ion drifts are neglected in the computation of the parallel electron drift, which is assumed to be much larger than any other flows on account of the small electron mass (and corresponding very large parallel electron mobility).

The steady-state momentum balance (Equations \ref{mombalance} and \ref{mombalancepar}) is justified in the sense that the principle charge carriers, electrons in the field-aligned direction, ions and electrons perpendicular to $\mathbf{B}$ attain a steady-state quickly. The ions drift perpendicular to the field in a steady state for time scales $\tau \gg 2 \pi \Omega_s^{-1}$. Electrons are in a steady-state in field-parallel direction on time-scales much longer than the electron collision time (probably $10^{-2}$ s or less) \citep[e.g.][]{Stmaurice:1996}.

To obtain an Ohm's law expression for use with equation \ref{divJSS} the momentum balance is multiplied by $n_s q_s$ and summed over species:

\begin{linenomath*} \begin{eqnarray}

\mathbf{J}_\perp &=& \boldsymbol{\sigma}_\perp \cdot \mathbf{E}_\perp - \sum_s \boldsymbol{\mu}_{s\perp} \cdot \nabla_\perp p_s + \sum_s n_s m_s \nu_s \boldsymbol{\mu}_{s\perp} \cdot \mathbf{v}_{n\perp} \label{ohm} \\

\mathbf{J}_\parallel &=& \sigma_0 \mathbf{E}_\parallel - \mu_{e0} \nabla_\parallel p_e \label{ohmpar}

\end{eqnarray} \end{linenomath*}

The conductivities are defined by:

\begin{linenomath*} \begin{eqnarray}

\boldsymbol{\sigma}_\perp

\equiv

\sum_s n_s q_s \boldsymbol{\mu}_{s\perp}

&=&

\left[ \begin{array}{ccc}

\sigma_P & -\sigma_H \\

\sigma_H & \sigma_P

\end{array} \right] \\

\sigma_0 &\approx& n_e q_e \mu_{e0}

%=

%\left[ \begin{array}{ccc}

% \sum_s \frac{q_s^2 n_s}{m_s \nu_{s}} \frac{\nu_{s}^2}{\nu_{s}^2 + \Omega_s^2} & \sum_s \frac{q_s^2 n_s}{m_s \nu_{s}} \frac{\nu_{s} \Omega_s}{\nu_{s}^2+\Omega_s^2} & 0 \\

% - \sum_s \frac{q_s^2 n_s}{m_s \nu_{s}} \frac{\nu_{s} \Omega_s}{\nu_{s}^2+\Omega_s^2} & \sum_s \frac{q_s^2 n_s}{m_s \nu_{s}} \frac{\nu_{s}^2}{\nu_{s}^2 + \Omega_s^2} & 0 \\

%0 & 0 & \sum_s \frac{n_s q_s^2}{m_s \nu_s}

%\end{array} \right]

\end{eqnarray} \end{linenomath*}

It is worth emphasizing that GEMINI does not treat the parallel ion drift with a steady state formulation (cf. Equation \ref{momentum}), while Equation \ref{ohmpar} clearly enforces some sort of parallel momentum steady state. This is justified by the fact that typical time scales for electron momentum variations are much shorter than the time scale for ion momentum variations. The former is roughly given by the inverse electron-ion collision frequency, while the latter is given by the inverse ion-electron collision frequency. Because of conservation of momentum in collisions \citep[c.f.][]{Schunk:2009} the electron-ion collision frequency is larger than the ion-electron frequency by a factor of roughly the ion to electron mass ratio, and leads to a much smaller momentum relaxation time for electrons as compared to ions. Hence, one may treat the electrons as in a steady state, while solving a time-dependent equation for ions and not render our current continuity inconsistent, provided that the ions do not contribute much to the parallel current (which, based on simple mobility arguments outlined here, they should not).

Therefore one may use the steady state (for electrons) implied by Equation \ref{ohmpar}, while still treating ion parallel momentum in a time-dependent and nonlinear manner.

\subsection{Handling of ambipolar and large-scale, ``background'' electric fields}

In GEMINI the ambipolar and background parts of the electric field is removed prior to solving for potential, and then added back into the momentum equation when the forces computed. This is done as a matter of convenience to avoid extra source terms in the numerical solution of our potential equation (described in later sections of this document). Thus, the electric field is divided into two parts, the ambipolar field $\mathbf{E}_a$ which is present even in the absence of magnetospheric forcing and neutral winds, and the response field $\mathbf{E}_r$ which results from magnetospherically imposed currents, potential patterns, or neutral winds. The parallel ambipolar field can be found directly by zeroing out $\mathbf{J}$ and $\mathbf{v}_n$ in Equation \ref{ohmpar}.

\begin{linenomath*} \begin{equation}

\mathbf{E}_{a\parallel} \approx \sigma_0^{-1} \mu_{e0} \nabla_\parallel p_e = \frac{1}{n_e q_e} \nabla_\parallel p_e \label{Eambpar}

\end{equation} \end{linenomath*}

GEMINI presently neglects contributions of perpendicular pressure gradients to current, but for simulations with very high resolution these may need to be added. Ohm's law with $\mathbf{E}=\mathbf{E}_r+\mathbf{E}_a$ becomes an equation in the disturbance quantities only.

\begin{linenomath*} \begin{eqnarray}

\mathbf{J}_\perp &=& \boldsymbol{\sigma}_\perp \cdot \mathbf{E}_{r\perp} + \left( \sum_s n_s m_s \nu_s \boldsymbol{\mu}_{s\perp} \cdot \mathbf{v}_{n\perp} \right) + \boldsymbol{\sigma}_\perp \cdot \mathbf{E}_{0\perp} \label{ohm2} \\

\mathbf{J}_\parallel &=& \sigma_0 \mathbf{E}_{r\parallel}

\end{eqnarray} \end{linenomath*}

The third term on the right-hand side of Equation \ref{ohm2} represents currents generated by large-scale ``background'' electric fields and convection, which is treated as a separate term in GEMINI. Invoking the electrostatic assumption $\mathbf{E}_r = -\nabla \Phi$ gives an equation that can be solved directly for electric potential.

\begin{linenomath*} \begin{equation}

\nabla_\perp \cdot \left( \boldsymbol{\sigma}_\perp \cdot \nabla_\perp \Phi \right) + \nabla_\parallel \cdot \left( \sigma_0 \nabla_\parallel \Phi \right) = \nabla_\perp \cdot \left( \sum_s n_s m_s \nu_s \boldsymbol{\mu}_{s\perp} \cdot \mathbf{v}_{n\perp} \right) + \nabla_\perp \cdot \left( \boldsymbol{\sigma}_\perp \cdot \mathbf{E}_{0\perp} \right) \label{divJSS2}

\end{equation} \end{linenomath*}

Note that this problem becomes homogeneous if neutral winds and background electric fields are neglected, while inclusion of the inhomogeneous terms allows for modeling of neutral drag effects on current generation and inclusion of large-scale background ionospheric electric fields and drifts.

%\subsection{Background vs. response electric fields}

%Often we wish to specify a uniform background field, which drives an ionospheric response due to e.g. nonuniform conductivies or other processes. This can be accommodated easily in the above formulation and slightly altered potential equations to solve:

%\begin{linenomath*} \begin{equation}

%\nabla_\perp \cdot \left( \boldsymbol{\sigma}_\perp \cdot \nabla_\perp \Phi \right) + \nabla_\parallel \cdot \left( \sigma_0 \nabla_\parallel \Phi \right) = \nabla_\perp \cdot \left( \sum_s n_s m_s \nu_s \boldsymbol{\mu}_{s\perp} \cdot \mathbf{v}_{n\perp} \right)

%\end{equation} \end{linenomath*}

%\subsection{On the mixing of time-dependent and steady state formulations}

%The apparent inconsistency that the reviewer notes can be resolved by carefully considering the assumptions underlying Equations 2, 7, and 10. These were not adequately discussed in the original manuscript so additional text has been added to the revision to correct this oversight.

%The electric field is separated, using superposition, into two parts in the model: (1) the ambipolar electric field Ea and (2) the disturbance response field Er. When solving for the electric potential in equation 7 the ambipolar electric field is not included, so the potential represents the part of the electric field that is due to the ionospheric resistivity (viz. collisions). When solving momentum equations, the ambipolar electric field (computed from electron pressure) is added to the electric field which is a result of driving a field-aligned current through the resistive ionosphere. Hence, both the electron pressure and resistive parts of the electric field are resolved in our model. The separation of different contributions to the electric field (ambipolar vs. resistive) is done merely as a matter of convenience (i.e. the pressure term could be left in the current continuity equation), and tests indicate that this approach does not impact solutions for electric field.

%As indicated above, the conductivities in equation 7 (now equation 8 - listed below for reference) are calculated using all species for the perpendicular conductivities (Pedersen and Hall) and only the electrons for the parallel conductivity.

%Note that this equation encodes assumptions of a perpendicular momentum steady state and a parallel electron momentum steady-state (further-discussed below). The disturbance response portion of the electric field can be calculated from the potential and then added to the ambipolar field to yield an expression for total, parallel electric field. Once the parallel field is calculated it can be used in the parallel momentum equation (which does, as the reviewer notes, include the important inertial term):

%The electron parallel drift is obtained by computing the current density from Ohm?s law and then applying the ?definition? of current density as indicated by equation 10 (now equation 12).

%As indicated above, our current continuity equation assumes steady state ion momentum perpendicular to the geomagnetic field and assumes that electrons carry the parallel current. Under a steady state momentum assumption ion and electron drifts are determined by the mobility and forces on each population. Since the electron mobility is much higher than the ion mobility in the parallel direction (see Figure 1?, not to be confused with the manuscript Figure 1, included below), the parallel ion drift contributions to current are neglected and the electrons are assumed to carry all of the parallel current. Even though the ions are not in a steady state in our model (e.g. equation 2) the assumption that they contribute negligibly to the parallel current is still likely to apply. Hence, we are able to treat the parallel ion momentum in a time-dependent manner (which is necessary for any sensible description of time-dependent ion upflow) while still using a current continuity equation that encodes a steady-state momentum balance.

\subsection{Leading-order electrodynamics} \label{sec:electrodynamics}

Electrodyamic effects on ionospheric electric fields can be encapsulated, to leading order, by include the polarization drifts and current in the current continuity equation, but still retaining the electrostatic assumption, $\nabla \times \mathbf{E}=0$. The polarization current, due to slow variations in the electric field with time, is assumed to be only in the perpendicular direction and is given by:

\begin{equation}

\mathbf{J}_p = \frac{\sum_s \rho_s}{B^2} \left( \frac{\partial}{\partial t} + \mathbf{v}_\perp \cdot \nabla_\perp \right) \mathbf{E}_r

\end{equation}

When included in the steady-state current continuity equation, along with electrostatic assumptions we get a time-dependent potential equation:

\begin{eqnarray}

\nabla_\perp \cdot \left( \boldsymbol{\sigma}_\perp \cdot \nabla_\perp \Phi \right) + \nabla_\parallel \cdot \left( \sigma_0 \nabla_\parallel \Phi \right) + \nabla_\perp \cdot \left[ c_M \left( \frac{\partial}{\partial t} + \mathbf{v}_\perp \cdot \nabla_\perp \right) \left( \nabla_\perp \Phi \right) \right] &=& \nonumber \\ \nabla_\perp \cdot \left( \sum_s n_s m_s \nu_s \boldsymbol{\mu}_{s\perp} \cdot \mathbf{v}_{n\perp} \right) + \nabla_\perp \cdot \left( \boldsymbol{\sigma}_\perp \cdot \mathbf{E}_{0\perp} \right) \label{eqn:electrodynamic}

\end{eqnarray}

The parallel component of the polarization current can be neglected for our purposes due to the very small parallel electric field and its extremely slow variations in time. The inertial capacitance $c_M$ is defined by \citep{Mitchell:1985}:

\begin{equation}

c_M \equiv \frac{\sum_s \rho_s}{B^2}

\end{equation}

The polarization effects can be derived by considering the ion momentum equation, include both inertia and time dependent terms neglected in prior static analysis, which can be written in the form:

\begin{equation}

m_s \left( \frac{\partial }{\partial t} + \mathbf{v}_s \cdot \nabla \right) \mathbf{v_s} = - \frac{1}{n_s} \nabla p_s + m_s \mathbf{g} + q_s \mathbf{E} + q_s \mathbf{v}_s \times \mathbf{B} + \sum_t m_s \nu_{st} \left(\mathbf{v}_t - \mathbf{v}_s \right)

\end{equation}

Taking the cross product of both sides of this equation with the magnetic field gives (note that a constant magnetic field is also assumed in the development below):

\begin{equation}

\mathbf{B} \times \left( \mathbf{v}_s \times \mathbf{B} \right) = \left\{ - m_s \left( \frac{\partial }{\partial t} + \mathbf{v}_s \cdot \nabla \right) \mathbf{v}_s - \frac{1}{n_s} \nabla p_s + m_s \mathbf{g} + q_s \mathbf{E} + \sum_t m_s \nu_{st} \left(\mathbf{v}_t - \mathbf{v}_s \right) \right \} \times \frac{\mathbf{B}}{q_s}

\end{equation}

\begin{equation}

\mathbf{v}_s (\mathbf{B} \cdot \mathbf{B}) - \mathbf{B}( \mathbf{v}_s \cdot \mathbf{B} ) = \left\{ - m_s \left( \frac{\partial }{\partial t} + \mathbf{v}_s \cdot \nabla \right) \mathbf{v}_s - \frac{1}{n_s} \nabla p_s + m_s \mathbf{g} + q_s \mathbf{E} + \sum_t m_s \nu_{st} \left(\mathbf{v}_t - \mathbf{v}_s \right) \right \} \times \frac{\mathbf{B}}{q_s}

\end{equation}

Assuming that the primary drift is perpendicular to the geomagnetic field:

\begin{equation}

\mathbf{v}_s = \left\{ - m_s \left( \frac{\partial }{\partial t} + \mathbf{v}_s \cdot \nabla \right) \mathbf{v_s} - \frac{1}{n_s} \nabla p_s + m_s \mathbf{g} + q_s \mathbf{E} + \sum_t m_s \nu_{st} \left(\mathbf{v}_t - \mathbf{v}_s \right) \right \} \times \frac{\mathbf{B}}{B^2 q_s}

\end{equation}

Which can then be rearranged to emphasize the terms of interest:

\begin{equation}

\mathbf{v}_{s\perp} = \frac{\mathbf{E} \times \mathbf{B}}{B^2} - \left\{ m_s \left( \frac{\partial }{\partial t} + \mathbf{v}_{s\perp} \cdot \nabla \right) \mathbf{v}_{s\perp} \right\} \times \frac{\mathbf{B}}{B^2 q_s} + \left\{ - \frac{1}{n_s m_s} \nabla p_s + \mathbf{g} + \sum_t \nu_{st} \left(\mathbf{v}_t - \mathbf{v}_{s\perp} \right) \right \} \times \frac{\mathbf{B}}{B^2 q_s}

\end{equation}

The first term on the right hand side of this equation is by far the dominant term:

\begin{equation}

\mathbf{v}_{s\perp} \approx \mathbf{v}_{E} = \frac{\mathbf{E} \times \mathbf{B}}{B^2}

\end{equation}

A higher-order correction to the ion drift can then be derived by replacing all instances of $\mathbf{v}_s$ with $\mathbf{v}_E$ on right hand side of the full drift equation:

\begin{equation}

\mathbf{v}_{s\perp} = \frac{\mathbf{E} \times \mathbf{B}}{B^2} - \left\{ m_s \left( \frac{\partial }{\partial t} + \mathbf{v}_{E} \cdot \nabla \right) \mathbf{v}_{E} \right\} \times \frac{\mathbf{B}}{B^2 q_s} + \left\{ - \frac{1}{n_s m_s} \nabla p_s + \mathbf{g} + \sum_t \nu_{st} \left(\mathbf{v}_t - \mathbf{v}_{E} \right) \right \} \times \frac{\mathbf{B}}{B^2 q_s}

\end{equation}

The time dependent term can be expanded out:

\begin{eqnarray}

- m_s \left\{ \left( \frac{\partial }{\partial t} + \mathbf{v}_{E} \cdot \nabla \right) \mathbf{v}_{E} \right\} \times \frac{\mathbf{B}}{B^2 q_s} &=& \frac{m_s}{B^2 q_s} \left( \frac{\partial }{\partial t} + \mathbf{v}_{E} \cdot \nabla \right) \left( \frac{\mathbf{E} \times \mathbf{B}}{B^2} \right) \times \mathbf{B} \\ &=& -\frac{m_s}{B^2 q_s} \left( \frac{\partial }{\partial t} + \mathbf{v}_{E} \cdot \nabla \right) \left( \frac{-1}{B^2} ( \mathbf{E} (B^2) - \mathbf{B} (\mathbf{E} \cdot \mathbf{B}) ) \right) \\ &=& \frac{m_s}{B^2 q_s} \left( \frac{\partial }{\partial t} + \mathbf{v}_{E} \cdot \nabla \right) \mathbf{E}

\end{eqnarray}

Thus we have, ignoring other terms:

\begin{equation}

\mathbf{v}_{s\perp} \approx \frac{\mathbf{E} \times \mathbf{B}}{B^2} + \frac{m_s}{B^2 q_s} \left( \frac{\partial }{\partial t} + \mathbf{v}_{E} \cdot \nabla \right) \mathbf{E}

\end{equation}

At lower altitude, obviously collisional terms also matter, as well. The current density resulting from these drifts is the polarization current density used in the GEMINI electrodynamic equation:

\begin{equation}

\mathbf{J}_p = \sum_s n_s q_s \mathbf{v}_{s\perp} = \frac{\sum n_s m_s }{B^2} \left( \frac{\partial }{\partial t} + \mathbf{v}_{E} \cdot \nabla \right) \mathbf{E}

\end{equation}

\subsection{Response magnetic fields}

Magnetic fields generated by ionospheric currents from dynamo sources (described above) are calculated from the Biot-Savart Law:

\begin{equation}

\mathbf{B}(\mathbf{x}) = \frac{\mu_0}{4 \pi} \int \frac{\mathbf{J}(\mathbf{x}') \times (\mathbf{x} - \mathbf{x}')}{|\mathbf{x} - \mathbf{x}'|^3} d^3 x'

\end{equation}

where $\mathbf{x}$ and $\mathbf{x'}$ represent field and source coordinates, respectively.

\subsection{Organization of model}

The GEMINI model consists of a simultaneous solution to Equations \ref{divJSS2}, \ref{continuity}, \ref{momentum}, and \ref{energy} for the unknowns $\Phi$, $n_s$, $\mathbf{v}_s$, and $T_s$. As outlined above, momentum balance is assumed in the computation of $\Phi$ at any particular instant. The evolution of potential on longer times scales is captured by evaluating the steady state equation Equation \ref{divJSS2} each time step as the fluid moment variables (densities, etc. hence conductivities) are updated. The result is a mathematical formulation capable of describing current closure, electric fields, heating, chemistry, plasma structuring, and upwelling processes in the auroral ionosphere.

\section{Use of generalized curvilinear coordinates in GEMINI}

GEMINI makes use of generalized curvilinear coordinates, denoted $x_1,x_2,x_3$, having, respectively, associated differential lengths, $h_1 dx_1,h_2 dx_2,h_3 dx_3$. Specific forms of the derivative operations needed in the model are documented in detail in this section.

%The master version of the code, while it uses these naming conventions, only functions for Cartesian coordinate (i.e. it assumes all metric factor are equal to one). The curvilinear branch includes all of the differential geometry described in this document.

\subsection{Gradient and divergence}

The gradient and divergence operations needed for both the fluid and electrostatic equations are given, in curvilinear form, by:

\begin{equation}

\nabla \cdot \mathbf{A} = \frac{1}{h_1 h_2 h_3} \left[ \frac{\partial}{\partial x_1} \left( h_2 h_3 A_1 \right) + \frac{\partial}{\partial x_2} \left( h_1 h_3 A_2 \right) + \frac{\partial}{\partial x_3} \left( h_1 h_2 A_3 \right) \right]

\end{equation}

for a vector field $\mathbf{A}$, and:

\begin{equation}

\nabla \Phi = \frac{1}{h_1} \frac{\partial \Phi}{\partial x_1} \hat{\mathbf{e}}_1 + \frac{1}{h_2} \frac{\partial \Phi}{\partial x_2} \hat{\mathbf{e}}_2 + \frac{1}{h_3} \frac{\partial \Phi}{\partial x_3} \hat{\mathbf{e}}_3

\end{equation}

for a scalar field $\Phi$

\subsection{Tensor derivatives in curvilinear form}

The bulk momentum fluxes occuring in the the parallel ion momentum equations (Equation \ref{momentum}) require careful treatment as they include various geometric terms which encapsulate, e.g. centrifugal forces. To proceed we specify the 1-coordinate ($x_1$) to be along the geomagnetic field line. The following identity regarding divergence of a symmetric rank-2 tensor ($\mathbf{T}$) in orthogonal coordinates is useful for paring down the bulk momentum flux term:

\begin{equation}

\left[ \nabla \cdot \mathbf{T} \right]_1 = \frac{1}{h_1 h_2 h_3} \left[ \frac{\partial}{\partial x_1}(h_2 h_3 T_{11}) + \frac{1}{h_1} \frac{\partial}{\partial x_2}(h_1^2 h_3 T_{12}) + \frac{1}{h_1} \frac{\partial}{\partial x_3}(h_1^2 h_2 T_{13}) \right] - \frac{T_{22}}{h_1 h_2} \frac{\partial h_2}{\partial x_1} - \frac{T_{33}}{h_1 h_3} \frac{\partial h_3}{\partial x_1}. \label{divtens}

\end{equation}

For our specific case this gives:

\begin{equation}

[\nabla \cdot (\rho \mathbf{v} \mathbf{v})]_1 = \frac{1}{h_1 h_2 h_3} \left[ \frac{\partial}{\partial x_1}(h_2 h_3 \rho v_1^2) + \frac{1}{h_1} \frac{\partial}{\partial x_2}(h_1^2 h_3 \rho v_1 v_2) + \frac{1}{h_1} \frac{\partial}{\partial x_3}(h_1^2 h_2 \rho v_1 v_3) \right] - \frac{\rho v_2^2}{h_1 h_2} \frac{\partial h_2}{\partial x_1} - \frac{\rho v_3^2}{h_1 h_3} \frac{\partial h_3}{\partial x_1}.

\end{equation}

In GEMINI an artificial viscosity (yielding an artificial stress) is included in the momentum and energy equations to deal with situations where steep solutions may form. The stress divergence parallel component (see Section \ref{sec:artvisc} and Equations \ref{eqn:stressmom} and \ref{eqn:stressen}) may also be needed for the momentum equation. As outlined above, we retain only the $Q_{11}$ element of the stress tensor, termed the parallel stress. Accordingly we may again invoke the identity of equation \ref{divtens}.

\begin{equation}

(\nabla \cdot \mathbf{Q})_1 = \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1}(h_2 h_3 Q_{11})

\end{equation}

If the parallel stress is given by a Navier-Stokes formula then we have:

\begin{equation}

Q_{11} = \left[-\mu_s \left( \nabla \boldsymbol{\mathrm{v}} + (\nabla \boldsymbol{\mathrm{v}})^T - \frac{2}{3} (\nabla \cdot \boldsymbol{\mathrm{v}}) \mathbf{I} \right) \right]_{11} = -\mu_s \left( 2 [\nabla \boldsymbol{\mathrm{v}}]_{11} - \frac{2}{3} (\nabla \cdot \boldsymbol{\mathrm{v}}) \right)

\end{equation}

In order to pursue this further we need the $\hat{\mathbf{e}}_1 \hat{\mathbf{e}}_1$ component of the gradient of the velocity in generalized orthogonal coordinates. The required identity is:

\begin{equation}

[\nabla \boldsymbol{\mathrm{v}}]_{11} = \frac{h_1 \frac{\partial v_1}{\partial x_1} - v_1 \frac{\partial h_1}{\partial x_1}}{h_1^2} + \frac{1}{h_1} \left( \frac{v_1}{h_1} \frac{\partial h_1}{\partial x_1} + \frac{v_2}{h_2} \frac{\partial h_1}{\partial x_2} + \frac{v_3}{h_3} \frac{\partial h_1}{\partial x_3} \right)

\end{equation}

In principle the above equations allow us to calculate the effects of stress, though these are presently neglected in the model since past studies have shown that NS stress approximations are invalid at altitudes where stress is significant \citep{Schunk:1975}.

\subsection{Curvilinear form of fluid equations}

The conservation of mass, momentum, and energy equations for the ions can be represented in general curvilinear coordinates as with the magnetic field being aligned with the $x_1-$direction. The continuity equation reads:

\begin{equation}

\frac{\partial \rho_s}{\partial t} + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( h_2 h_3 \rho_s v_{s1} \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_2} \left( h_1 h_3 \rho_s v_{s2} \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_3} \left( h_1 h_2 \rho_s v_{s3} \right) = m_s P_s - L_s \rho_s \label{eqn:continuitycoord}

\end{equation}

This equation has the form of an advection equation with local sources and sinks corresponding to chemical and photo-production terms.

The velocity has 3 components, hence, 3 equations. However, the vertical and horizontal parts are treated differently, as discussed above, with the result that only the time-dependence parallel momentum equation needs to be solved. %The horizontal drifts are controlled by short time-scale electrodynamic processes and are well approximated by:

%\begin{equation}

%v_{s2} \hat{\mathbf{e}}_2 + v_{s3} \hat{\mathbf{e}}_3 = \boldsymbol{\mu}_{s\perp} \cdot \mathbf{E}_\perp

%\end{equation}

%The field-aligned component of velocity is found by solving the full parallel momentum equation.

\begin{eqnarray}

\frac{\partial}{\partial t} \left( \rho_s v_{s1} \right) &+& \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( h_2 h_3 \rho_s v_{s1}^2 \right) + \frac{1}{h_1^2 h_2 h_3} \frac{\partial}{\partial x_2} \left( h_1^2 h_3 \rho_s v_{s1} v_{s2} \right) + \frac{1}{h_1^2 h_2 h_3} \frac{\partial}{\partial x_3} \left( h_1^2 h_2 \rho_s v_{s1} v_{s3} \right) \nonumber \\

~ &~& - \frac{\rho_s v_{s2}^2}{h_1 h_2} \frac{\partial h_2}{\partial x_1} - \frac{\rho_s v_{s3}^2}{h_1 h_3} \frac{\partial h_3}{\partial x_1} = \rho_s g_1 - \frac{1}{h_1} \frac{\partial p_s}{\partial x_1} + \frac{\rho_s q_s}{m_s} E_1 + \sum_t \rho_s \nu_{st} \left( v_{t1} - v_{s1} \right) \label{eqn:momentumcoord}

\end{eqnarray}

This equation has the form of an advection equation with local sources and sinks. The negative, geometric terms on the left-hand side can are treated as source terms for purposes of numerical solution.

The ion energy equation in curvilinear forms is:

\begin{eqnarray}

\frac{\partial}{\partial t} \left( \rho_s \epsilon_s \right) &+& \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( h_2 h_3 \rho_s \epsilon_s v_{s1} \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_2} \left( h_1 h_3 \rho_s \epsilon_s v_{s2} \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_3} \left( h_1 h_2 \rho_s \epsilon_s v_{s3} \right) \nonumber \\

&=& -p_s \frac{1}{h_1 h_2 h_3} \left( \frac{\partial}{\partial x_1} \left( h_2 h_3 v_{s1} \right) + \frac{\partial}{\partial x_2} \left( h_1 h_3 v_{s2} \right) + \frac{\partial}{\partial x_3} \left( h_1 h_2 v_{s3} \right) \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( \frac{h_2 h_3}{h_1} \lambda_s \frac{\partial T_s}{\partial x_1} \right) \nonumber \\

&~& - \frac{1}{\gamma_s - 1} \sum_t \frac{\rho_s k_B \nu_{st}}{m_s + m_t} \left[ 2 \left( T_s - T_t \right) - \frac{2 m_t}{3 k_B} \left( \left( v_{s1} - v_{t1} \right)^2 + \left( v_{s2} - v_{t2} \right)^2 + \left( v_{s3} - v_{t3} \right)^2 \right) \right] \label{eqn:ionenergycoord}

\end{eqnarray}

The electron energy has an extra term corresponding to the divergence of the thermoelectric heat flux.

\begin{equation}

- \nabla \cdot \left( - \beta_e \mathbf{J} \right)

\end{equation}

This term may be expanded using the product rules and simplified by enforcing a divergence free current:

\begin{equation}

- \nabla \cdot \left( - \beta_e \mathbf{J} \right) = \nabla \beta_e \cdot \mathbf{J} + \beta_e \nabla \cdot \mathbf{J} = \nabla \beta_e \cdot \mathbf{J}

\end{equation}

Finally if we consider only the parallel components of the heat flux and use the definition $\beta_e = \frac{5}{2} \frac{k_B T_e}{|q_e|}$, this becomes:

\begin{equation}

\nabla \beta_e \cdot \mathbf{J} = \frac{5}{2} \frac{k_B J_1}{|q_e|} \frac{1}{h_1} \frac{\partial T_e}{\partial x_1}

\end{equation}

The electron energy equation can then be written, in curvilinear coordinates as:

\begin{eqnarray}

\frac{\partial}{\partial t} \left( \rho_e \epsilon_e \right) &+& \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( h_2 h_3 \rho_e \epsilon_e v_{e1} \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_2} \left( h_1 h_3 \rho_e \epsilon_e v_{e2} \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_3} \left( h_1 h_2 \rho_e \epsilon_e v_{e3} \right) \nonumber \\

&=& -p_e \frac{1}{h_1 h_2 h_3} \left( \frac{\partial}{\partial x_1} \left( h_2 h_3 v_{e1} \right) + \frac{\partial}{\partial x_2} \left( h_1 h_3 v_{e2} \right) + \frac{\partial}{\partial x_3} \left( h_1 h_2 v_{e3} \right) \right) \nonumber \\

&~& + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( \frac{h_2 h_3}{h_1} \lambda_e \frac{\partial T_e}{\partial x_1} \right) + \frac{5}{2} \frac{k_B J_1}{|q_e|} \frac{1}{h_1} \frac{\partial T_e}{\partial x_1} \nonumber \\

&~& - \frac{1}{\gamma_e - 1} \sum_t \frac{\rho_e k_B \nu_{et}}{m_e + m_t} \left[ 2 \left( T_e - T_t \right) - \frac{2 m_t}{3 k_B} \left( \left( v_{e1} - v_{t1} \right)^2 + \left( v_{e2} - v_{t2} \right)^2 + \left( v_{e3} - v_{t3} \right)^2 \right) \right] \nonumber \\

&~& + \frac{Q_e}{\gamma_e - 1} \label{eqn:electronenergycoord}

\end{eqnarray}

\subsection{Curvilinear form of electrostatic equation}

\subsubsection{Electrostatic equation}

In GEMINI, the $x_1$-direction is taken to correspond to the field-aligned dimensions of the model grid. This dimension ordering makes parallelization easier, since the grid is most easily split amongst worker processes along the 3rd dimension or the arrays (this way subdomain portions of full-grid arrays are contiguous in memory). The curvilinear form of Equation \ref{divJSS2} may be written as:

\begin{eqnarray}

\frac{1}{h_1 h_2 h_3} \left[ \frac{\partial}{\partial x_1} \left( \frac{h_2 h_3}{h_1} \sigma_0 \frac{\partial \Phi}{\partial x_1} \right) + \frac{\partial}{\partial x_2} \left( \frac{h_1 h_3}{h_2} \sigma_P \frac{\partial \Phi}{\partial x_2} \right) - \frac{\partial}{\partial x_2} \left( h_1 \sigma_H \frac{\partial \Phi}{\partial x_3} \right) \right. &+& \nonumber \\ \left. \frac{\partial}{\partial x_3} \left( h_1 \sigma_H \frac{\partial \Phi}{\partial x_2} \right) + \frac{\partial}{\partial x_3} \left( \frac{h_1 h_2}{h_3} \sigma_P \frac{\partial \Phi}{\partial x_3} \right)

\right] = \nabla_\perp \cdot \left( \sum_s n_s m_s \nu_s \boldsymbol{\mu}_{s\perp} \cdot \mathbf{v}_{n\perp} \right) + \nabla_\perp \cdot \left( \boldsymbol{\sigma}_\perp \cdot \mathbf{E}_{0\perp} \right)

\end{eqnarray}

The right-hand side (RHS) of this equation is treated as a source term that is computed with a direct numerical derivative. This equation can be simplified by noting the equivalence of cross partial derivatives.

\begin{eqnarray}

-\frac{\partial}{\partial x_2} \left( h_1 \sigma_H \frac{\partial \Phi}{\partial x_3} \right) + \frac{\partial}{\partial x_3} \left( h_1 \sigma_H \frac{\partial \Phi}{\partial x_2} \right) &=& -\frac{\partial}{\partial x_2} \left( h_1 \sigma_H \right) \frac{\partial \Phi}{\partial x_3} - h_1 \sigma_H \frac{\partial^2 \Phi}{\partial x_2 \partial x_3} \\ &+& \frac{\partial}{\partial x_3} \left( h_1 \sigma_H \right) \frac{\partial \Phi}{\partial x_2} + h_1 \sigma_H \frac{\partial^2 \Phi}{\partial x_3 \partial x_2} \nonumber \\ &=& -\frac{\partial}{\partial x_2} \left( h_1 \sigma_H \right) \frac{\partial \Phi}{\partial x_3} + \frac{\partial}{\partial x_3} \left( h_1 \sigma_H \right) \frac{\partial \Phi}{\partial x_2},

\end{eqnarray}

which yields:

\begin{eqnarray}

\frac{1}{h_1 h_2 h_3} \left[ \frac{\partial}{\partial x_1} \left( \frac{h_2 h_3}{h_1} \sigma_0 \frac{\partial \Phi}{\partial x_1} \right) + \frac{\partial}{\partial x_2} \left( \frac{h_1 h_3}{h_2} \sigma_P \frac{\partial \Phi}{\partial x_2} \right) - \right. \frac{\partial}{\partial x_2} \left( h_1 \sigma_H \right) \frac{\partial \Phi}{\partial x_3} &+& \nonumber \\ \left. \frac{\partial}{\partial x_3} \left( h_1 \sigma_H \right) \frac{\partial \Phi}{\partial x_2} + \frac{\partial}{\partial x_3} \left( \frac{h_1 h_2}{h_3} \sigma_P \frac{\partial \Phi}{\partial x_3} \right)

\right] = \nabla_\perp \cdot \left( \sum_s n_s m_s \nu_s \boldsymbol{\mu}_{s\perp} \cdot \mathbf{v}_{n\perp} \right) + \nabla_\perp \cdot \left( \boldsymbol{\sigma}_\perp \cdot \mathbf{E}_{0\perp} \right) \label{eqn:curvpot}

\end{eqnarray}

This particular form of the potential equation is useful, numerically, since it does not contain any cross-partial derivatives. It is possible to rewrite the right-hand side of this equation, which is numerically useful as it reduces the message-passing burden in the model by eliminating the need to exchange individual ion mobilities between worker processes.

\begin{equation}

\nabla_\perp \cdot \left( \sum_s n_s m_s \nu_s \boldsymbol{\mu}_{s\perp} \cdot \mathbf{v}_{n\perp} \right) = \nabla_\perp \cdot \left[ \boldsymbol{\sigma}_\perp \cdot \left( \mathbf{v}_{n\perp} \times \mathbf{B} \right) \right]

\end{equation}

\subsubsection{Equipotential field line (EFL) potential formulation}

In many cases it is possible (and useful, from the standpoint of computational efficiency) to solve a field-integrated form of the potential equation, corresponding to an assumption of equipotential field lines (the EFL approximation). In situations where the scale sizes perpendicular to the geomagnetic field are not too small, the field lines can effectively be treated as equipotentials \citep{Farley:1959,Huba:1988}. Rerranging and integrating Equation \ref{eqn:curvpot} along the x$_1$-direction) yields:

\begin{eqnarray}

&~& \frac{\partial}{\partial x_2} \left( \Sigma_P^{(2)} \frac{\partial \Phi}{\partial x_2} \right) + \frac{\partial}{\partial x_3} \left( \Sigma_P^{(3)} \frac{\partial \Phi}{\partial x_3} \right) - \frac{\partial \Sigma_H'}{\partial x_2} \frac{\partial \Phi}{\partial x_3} + \frac{\partial \Sigma_H'}{\partial x_3} \frac{\partial \Phi}{\partial x_2} \nonumber \\ &=& \int h_1 h_2 h_3 \left\{ \nabla_\perp \cdot \left[ \boldsymbol{\sigma}_\perp \cdot \left( \mathbf{v}_{n\perp} \times \mathbf{B} \right) \right] + + \nabla_\perp \cdot \left( \boldsymbol{\sigma}_\perp \cdot \mathbf{E}_{0\perp} \right) \right\} d x_1 + \left[ h_2 h_3 J_1 \right] \left|^{x_{1,max}}_{x_{1,min}} \right. \label{eqn:curvpotEFL}

\end{eqnarray}

where the following definitions have been used:

\begin{eqnarray}

\Sigma_P^{(2)} &\equiv& \int \frac{h_1 h_3}{h_2} \sigma_P dx_1 \\

\Sigma_P^{(3)} &\equiv& \int \frac{h_1 h_2}{h_3} \sigma_P dx_1 \\

\Sigma_H' &\equiv& \int h_1 \sigma_H dx_1

\end{eqnarray}

and the following simplifications were made:

\begin{equation}

\int \frac{\partial}{\partial x_1} \left( \frac{h_2 h_3}{h_1} \sigma_0 \frac{\partial \Phi}{\partial x_1} \right) d x_1= \left[ h_2 h_3 J_1 \right] \left|^{x_{1,max}}_{x_{1,min}} \right.

\end{equation}

Lateral boundaries (corresponding to the two field-perpendicular directions) are always assumed to be at a fixed potential, hence, Equation \ref{eqn:curvpotEFL} is always solved with Dirichlet conditions. Note also that if the user specifies a topside potential pattern that it is not necessary to solve Equation \ref{eqn:curvpotEFL} since the potential is already known and the fields and current densities may be computed directly. In contrast, with a specified field aligned current at the highest altitudes of the model, the potential is not known \emph{a priori} and Equation \ref{eqn:curvpotEFL} must to be solved. This the the mode of operation most typically used in GEMINI, particularly for 3D simulations where a field-resolved potential solution is not practical, especially in high-resolution simulations. %maybe move to numerical section?

\subsubsection{Electrodynamic equation in EFL form}

Electrodynamic contributions to electric fields and currents are approximated using the polarization current as outlined in Section \ref{sec:electrodynamics}. Including this current in the field-integrated form of the ionospheric potential equation reads:

\begin{equation}

\nabla_\perp \cdot \left( \boldsymbol{\Sigma}_\perp \cdot \nabla_\perp \Phi \right) + \nabla_\perp \cdot \left[ C_M \left( \frac{\partial}{\partial t} + \mathbf{v}_\perp \cdot \nabla_\perp \right) \left( \nabla_\perp \Phi \right) \right] = \nonumber \\ \nabla_\perp \cdot \left[ \boldsymbol{\sigma}_\perp \cdot \left( \mathbf{v}_{n\perp} \times \mathbf{B} \right) \right] + \nabla_\perp \cdot \left( \boldsymbol{\Sigma}_\perp \cdot \mathbf{E}_{0\perp} \right) \label{eqn:electro}

\end{equation}

Where $C_M$ is the field-line integrated inertial capacitance. In the current version of GEMINI, solution of this equation is only supported in Cartesian coordinates - attempts to use electrodynamics with other coordinate system will cause GEMINI to terminate the simulation without attempting a solution (an error will also be printed to the screen).

In Cartesian coordinates (the only system supported for electrodynamic simulations) the electrodynamic, EFL form of the polarization current is:

\begin{eqnarray}

- \nabla_\perp \cdot \left[ C_M \frac{\partial}{\partial t} \left( \nabla_\perp \Phi \right) \right] - \nabla_\perp \cdot \left[ C_M \left( \mathbf{v}_\perp \cdot \nabla_\perp \right) \left( \nabla_\perp \Phi \right) \right] &=& \nonumber \\

- \left\{ \frac{\partial}{\partial x_2} \left[ C_M \frac{\partial}{\partial t} \left( \frac{\partial \Phi}{\partial x_2} \right) \right] + \frac{\partial}{\partial x_3} \left[ C_M \frac{\partial}{\partial t} \left( \frac{\partial \Phi}{\partial x_3} \right) \right] \right\} &-& \nonumber \\

\frac{\partial}{\partial x_2} \left( C_M v_2 \frac{\partial^2 \Phi}{\partial x_2^2} + C_M v_3 \frac{\partial^2 \Phi}{\partial x_3 \partial x_2} \right) &-& \nonumber \\

\frac{\partial}{\partial x_3} \left( C_M v_2 \frac{\partial^2 \Phi}{\partial x_2 \partial x_3} + C_M v_3 \frac{\partial^2 \Phi}{\partial x_3^2} \right)

\end{eqnarray}

Inserting this into the full EFL potential PDE results in a rather more complex potential equation to be solved as compared to the purely electrostatic form:

\begin{eqnarray}

\frac{\partial}{\partial x_2} \left( \Sigma_P \frac{\partial \Phi}{\partial x_2} \right) + \frac{\partial}{\partial x_3} \left( \Sigma_P \frac{\partial \Phi}{\partial x_3} \right) - \frac{\partial \Sigma_H}{\partial x_2} \frac{\partial \Phi}{\partial x_3} + \frac{\partial \Sigma_H}{\partial x_3} \frac{\partial \Phi}{\partial x_2} &+& \nonumber \\

\left\{ \frac{\partial}{\partial x_2} \left[ C_M \frac{\partial}{\partial t} \left( \frac{\partial \Phi}{\partial x_2} \right) \right] + \frac{\partial}{\partial x_3} \left[ C_M \frac{\partial}{\partial t} \left( \frac{\partial \Phi}{\partial x_3} \right) \right] \right\} &+& \nonumber \\

\frac{\partial}{\partial x_2} \left( C_M v_2 \frac{\partial^2 \Phi}{\partial x_2^2} + C_M v_3 \frac{\partial^2 \Phi}{\partial x_3 \partial x_2} \right) &+& \nonumber \\

\frac{\partial}{\partial x_3} \left( C_M v_2 \frac{\partial^2 \Phi}{\partial x_2 \partial x_3} + C_M v_3 \frac{\partial^2 \Phi}{\partial x_3^2} \right) &=& \nonumber \\

= \int \left\{ \nabla_\perp \cdot \left[ \boldsymbol{\sigma}_\perp \cdot \left( \mathbf{v}_{n\perp} \times \mathbf{B} \right) \right] + + \nabla_\perp \cdot \left( \boldsymbol{\Sigma}_\perp \cdot \mathbf{E}_{0\perp} \right) \right\} d x_1 + \left[J_1 \right] \left|^{x_{1,max}}_{x_{1,min}} \right. \label{eqn:electrodynamic}

\end{eqnarray}

\subsubsection{Field-aligned current calculations}

When an EFL approximation is used, the field aligned current calculation differs from what it would normally be in the field-resolved formulation, i.e. $\mathbf{J}_\parallel=-\sigma_0 \nabla_\parallel \Phi$, since we explicity assume $\nabla_\parallel \Phi = 0$ in formulating the EFL equations. The field-aligned currents are instead computed from the divergence of the perpendicular current. Specifically the parallel current is computed to be whatever current is necessary to enforce the divergence free conditions given the perpendicular currents generated by the solution of Equation \ref{eqn:curvpotEFL}. Enforcing current continuity with a known perpendicular current density yields:

\begin{equation}

\left. \left[ h_2 h_3 J_1 \right] \right|_{x_1'}^{x_1''} = - \int_{x_1'}^{x_{1}''} h_1 h_2 h_3 \left( \nabla_\perp \cdot \mathbf{J}_\perp \right) d x_1 \label{eqn:FACcalc}

\end{equation}

where $x_1'$ and $x_1''$ represent two different locations on the computational grid.

Note that it is \emph{always} necessary to know the current density at one of the model $x_1$ locations (ostensibly at a boundary) to produce a solution for parallel current over the entired grid via Equation \ref{eqn:FACcalc}. Two different choices of ``top'' and "bottom" boundary conditions (i.e. those in the $x_1$-direction) can be used to specify the current somewhere on the grid: (1) a user-specified a potential pattern (``Dirichlet'' conditions) and (2) a user-specified field-aligned current pattern (``Neumann'' conditions). If ``Dirichlet'' conditions are the perpendicular current may be directly calculated and used with the assumption that the field-aligned current at the lowest altitude of the model is zero to produce an expression for the parallel current throughout the domain. In the case that the lowest altitude corresponds to the maximum $x_1$ value we have:

\begin{equation}

\left. \left[ h_2 h_3 J_1 \right] \right|_{x_1} = \left. \left[ h_2 h_3 J_1 \right] \right|_{x_{1,max}} + \int_{x_1}^{x_{1,max}} h_1 h_2 h_3 \left( \nabla_\perp \cdot \mathbf{J}_\perp \right) d x_1

\end{equation}

If the lowest altitude corresponds to the minimum $x_1$ value, the parallel current over the rest of the grid can be calculated from:

\begin{equation}

\left. \left[ h_2 h_3 J_1 \right] \right|_{x_1} = \left. \left[ h_2 h_3 J_1 \right] \right|_{x_{1,min}} - \int_{x_{1,min}}^{x_1} h_1 h_2 h_3 \left( \nabla_\perp \cdot \mathbf{J}_\perp \right) d x_1

\end{equation}

Note that in either case, we are explicitly using knowledge of the parallel current at one of the boundaries. On a closed dipole grid, either equation immediately above may be used since the current density is known to go to zero at both ends. Note that the above formulas may also be used when the user specifies a top boundary current (``Neumann'' conditions) as may be the case when an open grid (i.e. not spanning both hemispheres) is used.

\subsubsection{Magnetic fields}

The particular integral for the magnetic fields is always evaluated in a Cartesian form, irrespective of the coordinate system used in the simulations. To facilitate this, simulation results are interpolated and rotated onto first into a spherical ECEF coordinate system, then a local ($x_{1,2,3}$ = up, south, east) coordinate system. The component form of the magnetic field for a 3D simulation reads:

\begin{equation}

\mathbf{B}(\mathbf{x}) = \frac{\mu_0}{4 \pi} \left\{ \hat{\mathbf{e}}_1 \int \frac{J_2 R_3 - J_3 R_2}{R^3} d^3 x' - \hat{\mathbf{e}}_2 \int \frac{J_1 R_3 - J_3 R_1}{R^3} d^3 x' + \hat{\mathbf{e}}_3 \int \frac{J_1 R_2 - J_2 R_1}{R^3} d^3 x' \right\} \label{eqn:magint3D}

\end{equation}

where $\mathbf{R} \equiv \mathbf{x}-\mathbf{x'} = \sum_i R_i \hat{\mathbf{e}}_i$ is the displacement vector between source and field grid points.

\section{Dipole coordinates}

The most commonly used \emph{curvilinear} coordinate system for GEMINI is a dipole coordinate system \citep{Huba:2000}. For purposes of illustrating how one might set up an arbitrary orthogonal curvilinear system for GEMINI we present and end-to-end derivation of the dipole metric and associated basis vectors here. A similar procedure may be used to fully characterize any other curvilinear system that is orthogonal.

\subsection{Derivation of metric for dipole coordinates}

The definitions for dipole coordinates used in GEMINI are (switching to superscripts to denote contravariant vs. covariant):

\begin{equation}

x^1 = q; \quad x^2 = p; \quad x^3 = \phi

\end{equation}

with the $q,p,\phi$ coordinates defined as in \citet{Huba:2000}, and $\phi$ being the same coordinate as used in a spherical system. In terms of spherical coordinates (measured from the center of Earth's tilted dipole) and Cartesian coordinates:

\begin{eqnarray}

q &=& \frac{R_e^2}{r^2} \cos \theta = \frac{R_e^2 z}{\left( x^2+y^2+z^2\right)^{3/2}}\\

p &=& \frac{r}{R_e} \frac{1}{\sin^2 \theta} = \frac{\left( x^2+y^2+z^2\right)^{3/2}}{R_e \left( x^2 + y^2 \right) } \\

\phi &=& \tan^{-1} \left( \frac{y}{x} \right)

\end{eqnarray}

Inverting $q,p$ coordinates to obtain $r,\theta$ is accomplished by solving the relations:

\begin{eqnarray}

q^2 \left( \frac{r}{R_e}\right)^4 + \frac{1}{p} \left( \frac{r}{R_e}\right) &=& 1 \\

\cos \theta &=& q \frac{r^2}{R_e^2}

\end{eqnarray}

From these transformations, the basis vectors and metric for dipole coordinates can be derived. The contravariant basis vectors are defined in terms of contravariant differentials $d x^i$ as:

\begin{equation}

\boldsymbol{\varepsilon}^i = \frac{\partial x^i}{\partial x} \hat{\mathbf{e}}_x + \frac{\partial x^i}{\partial y} \hat{\mathbf{e}}_y + \frac{\partial x^i}{\partial z} \hat{\mathbf{e}}_z \label{eqn:conbasis}

\end{equation}

The covariant basis vectors (not used immediately) are defined by:

\begin{equation}

\boldsymbol{\varepsilon}_i = \frac{\partial x}{\partial x^i} \hat{\mathbf{e}}_x + \frac{\partial y}{\partial x^i} \hat{\mathbf{e}}_y + \frac{\partial z}{\partial x^i} \hat{\mathbf{e}}_z \label{eqn:conbasis}

\end{equation}

The necessary derivative for computing the contravariant basis be computed from the dipole coordinates as functions of the cartesian coordinates $x,y,z$ (listed above) and are given by:

\begin{eqnarray}

\frac{\partial x^1}{\partial x} &=& \frac{\partial q}{\partial x} = \frac{-3 R_e \sin \theta \cos \theta \cos \phi}{r^3} \\

\frac{\partial x^1}{\partial y} &=& \frac{\partial q}{\partial y} = \frac{-3 R_e \sin \theta \cos \theta \sin \phi}{r^3} \\

\frac{\partial x^1}{\partial z} &=& \frac{\partial q}{\partial z} = \frac{R_e^2 \left( 1 - 3 \cos^2 \theta \right)}{r^3} \\

\frac{\partial x^2}{\partial x} &=& \frac{\partial p}{\partial x} = \frac{\cos \phi}{ R_e \sin^3 \theta} \left( 1 - 3 \cos^2 \theta \right) \\

\frac{\partial x^2}{\partial y} &=& \frac{\partial p}{\partial y} = \frac{\sin \phi}{ R_e \sin^3 \theta} \left( 1 - 3 \cos^2 \theta \right) \\

\frac{\partial x^2}{\partial z} &=& \frac{\partial p}{\partial z} = \frac{3 \cos \theta}{R_e \sin^2 \theta}\\

\frac{\partial x^3}{\partial x} &=& \frac{\partial \phi}{\partial x} = \frac{- \sin \phi}{r \sin \theta} \\

\frac{\partial x^3}{\partial y} &=& \frac{\partial \phi}{\partial y} = \frac{\cos \phi}{r \sin \theta} \\

\frac{\partial x^3}{\partial z} &=& \frac{\partial \phi}{\partial z} = 0

\end{eqnarray}

The contravariant basis vectors are, therefore, from Equation \ref{eqn:conbasis}:

\begin{eqnarray}

\boldsymbol{\varepsilon}^1 &=& \frac{-3 R_e \sin \theta \cos \theta \cos \phi}{r^3} \hat{\mathbf{e}}_x + \frac{-3 R_e \sin \theta \cos \theta \sin \phi}{r^3} \hat{\mathbf{e}}_y + \frac{R_e^2 \left( 1 - 3 \cos^2 \theta \right)}{r^3} \hat{\mathbf{e}}_z \\

\boldsymbol{\varepsilon}^2 &=& \frac{\cos \phi}{ R_e \sin^3 \theta} \left( 1 - 3 \cos^2 \theta \right) \hat{\mathbf{e}}_x + \frac{\sin \phi}{ R_e \sin^3 \theta} \left( 1 - 3 \cos^2 \theta \right) \hat{\mathbf{e}}_y + \frac{3 \cos \theta}{R_e \sin^2 \theta} \hat{\mathbf{e}}_z \\

\boldsymbol{\varepsilon}^3 &=& \frac{- \sin \phi}{r \sin \theta} \hat{\mathbf{e}}_x + \frac{\cos \phi}{r \sin \theta} \hat{\mathbf{e}}_y

\end{eqnarray}

The metric tensors corresponding to our dipole curvilinear space are related to differential displacements in the new coordinate system. Based on the definition of the covariant unit vectors a differential displacement in a curvilinear system is:

\begin{equation}

d \mathbf{r} = d x^i \boldsymbol{\varepsilon}_i \label{eqn:dr}

\end{equation}

Differential distance squared, is then:

\begin{equation}

dr^2 \equiv d \mathbf{r} \cdot d \mathbf{r} = \left(d x^i \boldsymbol{\varepsilon}_i \right) \cdot \left( d x^j \boldsymbol{\varepsilon}_j \right) = \left( \boldsymbol{\varepsilon}_i \cdot \boldsymbol{\varepsilon}_j \right) dx^i dx^j = g_{ij} dx^i dx^j,

\end{equation}

where

\begin{equation}

g_{ij} \equiv \boldsymbol{\varepsilon}_i \cdot \boldsymbol{\varepsilon}_j

\end{equation}

is defined as the covariant metric. By analogy we can define a contravariant metric:

\begin{equation}

g^{ij} \equiv \boldsymbol{\varepsilon}^i \cdot \boldsymbol{\varepsilon}^j

\end{equation}

Note that, as a consequence of the chain rule applied to sequences of derivatives from the basis vectors, the following relations hold:

\begin{equation}

\boldsymbol{\varepsilon}_i \cdot \boldsymbol{\varepsilon}^j = \delta_i^j

\end{equation}

\begin{equation}

g_{ik} g^{kj} = \delta_i^j

\end{equation}

From the latter expression we see that the covariant and contravariant metrics can be understood as matrix inverses of each other. One may convert from covariant to contravariant bases by using these two metrics:

\begin{equation}

\boldsymbol{\varepsilon}_j = g_{ji} \boldsymbol{\varepsilon}^i

\end{equation}

\begin{equation}

\boldsymbol{\varepsilon}^j = g^{ji} \boldsymbol{\varepsilon}_i

\end{equation}

Components of vectors may similarly be transformed through the metric:

\begin{equation}

A_j = g_{ji} A^i

\end{equation}

\begin{equation}

A^j = g^{ji} A_i

\end{equation}

Because the dipole coordinate system is orthogonal, there are only three non-zero elements of the contravariant metric, $g^{ii}$ (no implied summation), and we have:

\begin{eqnarray}

g^{ij} \rightarrow \left[

\begin{array}{ccc}

\frac{R_e^4}{r^6} \left( 1 + 3 \cos \theta \right) & 0 & 0 \\

0 & \frac{1}{R_e^2} \frac{\left( 1+3 \cos \theta \right)}{\sin^6 \theta} & 0 \\

0 & 0 & \frac{1}{r^2 \sin^2 \theta}

\end{array}

\right]

\end{eqnarray}

This is related to the covariant metric $g_{ij}$ via the relation:

\begin{equation}

g_{ik} g^{kj} = \delta_i^j

\end{equation}

Thus the covariant metric is the inverse of the contravariant metric:

\begin{eqnarray}

g_{ij} \rightarrow \left[

\begin{array}{ccc}

\frac{r^6}{R_e^4 \left( 1 + 3 \cos \theta \right)} & 0 & 0 \\

0 & \frac{R_e ^2 \sin^6 \theta}{\left( 1+3 \cos \theta \right)} & 0 \\

0 & 0 & r^2 \sin^2 \theta

\end{array}

\right]

\end{eqnarray}

The metric factors, $h_i$, (needed for derivatives and intergrals) are defined by (no summation):

\begin{equation}

h_i \equiv \sqrt{g_{ii}}

\end{equation}

For the dipole system, applying this definition gives:

\begin{eqnarray}

h_1 = h_q &=& \frac{r^3}{R_e^2} \frac{1}{\sqrt{1+3 \cos^2 \theta}} \\

h_2 = h_p &=& \frac{R_e \sin^3 \theta}{\sqrt{1+3 \cos^2 \theta}} \\

h_3 = h_\phi &=& r \sin \theta

\end{eqnarray}

Finally, dipole unit vectors may be derived by noting equivalence of various expressions for $d \mathbf{r}$:

\begin{equation}

dr^2 = g_{ij} dx^i dx^j

\end{equation}

In orthogonal coordinates this simplifies to:

\begin{equation}

dr ^2 = g_{11} dx^1 dx^1 + g_{22} dx^2 dx^2 + g_{33} dx^3 dx^3 = \sum_i h_i^2 (dx_i)^2

\end{equation}

where in the last equality we have reverted to explicit summation and use of pure subscript indices. It is notable in this equation that the elements of the summation $h_i^2 dx_i^2$ can be understood as differential distance (squared) in the $i^{th}$ direction. Thus, we arrive at yet another expression for $d \mathbf{r}$:

\begin{equation}

d \mathbf{r} = \sum_i h_i dx^i \hat{\mathbf{e}}_i

\end{equation}

The $\hat{\mathbf{e}}_i$ vectors here are true unit vectors since the differentials accompanying them $h_i dx^i$ have units of true distance (units of meters). Comparing with Equation \ref{eqn:dr}, we see that (no summation):

\begin{equation}

\hat{\mathbf{e}}_i = \frac{1}{h_i} \boldsymbol{\varepsilon}_i

\end{equation}

Alternatively we may write down the units vectors in terms of the contravariant basis (which could be easier to compute) (no sum over $i$ - only $j$):

\begin{equation}

\hat{\mathbf{e}}_i = \frac{g_{ij}}{h_i} \boldsymbol{\varepsilon}^j

\end{equation}

In an orthogonal system only the diagonal elements of the metric are present, so that the unit vectors become (no summation):

\begin{equation}

\hat{\mathbf{e}}_i = \frac{g_{ii}}{h_i} \boldsymbol{\varepsilon}^i = h_i \boldsymbol{\varepsilon}^i

\end{equation}

Invoking this result with the above-derived contravariant basis vectors gives:

\begin{eqnarray}

\hat{\mathbf{e}}_1 &=& \frac{-3 \cos \theta \sin \theta \cos \phi}{\sqrt{1+3 \cos^2 \theta } } ~ \hat{\mathbf{e}}_x + \frac{-3 \cos \theta \sin \theta \sin \phi}{\sqrt{1+3 \cos^2 \theta } } ~ \hat{\mathbf{e}}_y + \frac{\left( 1 - 3 \cos^2 \theta \right) }{\sqrt{1+3 \cos^2 \theta } }~ \hat{\mathbf{e}}_z \\

\hat{\mathbf{e}}_2 &=& \frac{\cos \phi \left(1 - 3 \cos^2 \theta \right)}{\sqrt{1 + 3 \cos^2 \theta}} ~ \hat{\mathbf{e}}_x + \frac{\sin \phi \left(1 - 3 \cos^2 \theta \right)}{\sqrt{1 + 3 \cos^2 \theta}}~ \hat{\mathbf{e}}_y + ~ \frac{\cos \theta \sin \theta}{\sqrt{1+3 \cos^2 \theta }} \hat{\mathbf{e}}_z \\

\hat{\mathbf{e}}_3 &=& - \sin \phi ~ \hat{\mathbf{e}}_x + \cos \phi ~ \hat{\mathbf{e}}_y

\end{eqnarray}

\subsection{Tilted dipole approximations}

The $r,\theta,\phi$ coordinates of a point in space, are understood to be spherical coordinates measured from the center of the Earth's tilted magnetic dipole. Unit vectors used for grid-related calculations are given, in tilted dipole, Cartesian coordinates by:

\begin{eqnarray}

\hat{\mathbf{e}}_r &=& \sin \theta \cos \phi ~ \hat{\mathbf{e}}_x + \sin \theta \sin \phi ~ \hat{\mathbf{e}}_y + \cos \theta ~ \hat{\mathbf{e}}_z \\

\hat{\mathbf{e}}_\theta &=& \cos \theta \cos \phi ~ \hat{\mathbf{e}}_x + \cos \theta \sin \phi ~ \hat{\mathbf{e}}_y - \sin \theta ~ \hat{\mathbf{e}}_z \\

\hat{\mathbf{e}}_\phi &=& -\sin \phi ~ \hat{\mathbf{e}}_x + \cos \phi ~ \hat{\mathbf{e}}_y

%\hat{\mathbf{e}}_q &\equiv& \frac{\frac{d \mathbf{x}}{d q}}{\left| \frac{d \mathbf{x}}{d q} \right|} \nonumber \\

% &=& \frac{-3 \cos \theta \sin \theta }{\sqrt{1+3 \cos^2 \theta } } \cos \phi ~ \hat{\mathbf{e}}_x + \frac{-3 \cos \theta \sin \theta }{\sqrt{1+3 \cos^2 \theta } } \sin \phi ~ \hat{\mathbf{e}}_y + \frac{-2 \cos^2 \theta + \sin^2 \theta }{\sqrt{1+3 \cos^2 \theta } } ~ \hat{\mathbf{e}}_z \\

%\hat{\mathbf{e}}_p &=& \hat{\mathbf{e}}_\phi \times \hat{\mathbf{e}}_q

\end{eqnarray}

%In deriving the $\hat{\mathbf{e}}_q$ unit vector use was made of the fact that:

%\begin{equation}

%\frac{d \mathbf{x}}{d q} = \frac{d}{d q} \left( r \sin \theta \cos \phi \hat{\mathbf{e}}_x + r \sin \theta \sin \phi \hat{\mathbf{e}}_y + r \cos \theta \hat{\mathbf{e}}_z \right)

%\end{equation}

%Dervitives needed to evaluate this unit vector are give by \citet{Schunk:2009}:

%\begin{eqnarray}

%\frac{d r}{d q} &=& \frac{-2 r^3 \cos \theta}{R_e^2 \left( 1 + 3 \cos^2 \theta \right)}\\

%\frac{d \theta}{d q} &=& \frac{-r^2 \sin \theta}{R_e^2 \left( 1 + 3 \cos^2 \theta \right)} \\

%\frac{d \phi}{d q} &=& 0 \\

%\end{eqnarray}

The magnetic field is computed from dipole formulas in spherical coordinates, i.e.:

\begin{equation}

\mathbf{B}=\frac{\mu_0 m}{4 \pi r^3} \left(2 \cos \theta ~ \hat{\mathbf{e}}_r + \sin \theta ~ \hat{\mathbf{e}}_\theta \right)

\end{equation}

where $m=7.94 \times 10^{22}$ (A m$^2$) is the magnetic moment of the Earth. The magnetic field magnitude (and therefore $q$ and $x_1$ component is then given by:

\begin{equation}

B=\frac{\mu_0 m}{4 \pi r^3} \sqrt{1+3 \cos^2 \theta}

\end{equation}

THE FOLLOWING IS TAKEN DIRECTLY FROM THE SOURCE CODE AND IS NOT THE SAME AS WHAT IS PRESENTED IN HUBA'S PAPER... NEED TO LOOK INTO THIS...

To define a transformation from spherical dipole coordinates (i.e. geomagnetic coordinates) to geographic coordinates one must define the Earth's magnetic dipole tilt angles, taken in GEMINI to be:

\begin{eqnarray}

\theta_0 &=& 11.435^\circ \nonumber \\

\phi_0 &=& 290.24^\circ \nonumber

\end{eqnarray}

These are the zenith and azimuth angles, respectively, measured from the geographic (i.e. spin) axes to the magnetic axes. The polar angle in the tilted dipole coordinate system (i.e. $\phi$) is then given by:

\begin{equation}

\theta = \cos^{-1} \left( \cos \theta_0 \cos \theta_g + \sin \theta_0 \sin \theta_g \cos \left( \phi_g - \phi_0 \right) \right)

\end{equation}

where $\theta_g,\phi_g$ are the angle of the location of interest in geographic coordinates. The azimuth angle is then given by:

\begin{equation}

\alpha = \cos^{-1} \left( \frac{ \cos \theta_g - \cos \theta \cos \theta_0 }{\sin \theta \sin \theta_0} \right)

\end{equation}

\begin{equation}

\phi=

\begin{array}{c}

\pi - \alpha \quad (\phi_0 > \phi_g ~\land~ \phi_0-\phi_g > \pi ~~\lor~~ \phi_0 < \phi_g ~\land~ \phi_g - \phi_0 < \pi) \\

\alpha - \pi \quad (\mathrm{otherwise})

\end{array}

\end{equation}

The reverse transformations from geomagnetic to geographic are given by:

\begin{eqnarray}

\theta_g &=& \cos^{-1} \left( \cos \theta \cos \theta_0 - \sin \theta \sin \theta_0 \cos \phi \right) \\

\beta &=& \cos^{-1} \left( \frac{\cos \theta - \cos \theta_0 \cos \theta_g}{\sin \theta_0 \sin \theta_g} \right)

\end{eqnarray}

\begin{equation}

\phi_g=

\begin{array}{c}

\phi_0 - \beta \quad (\phi > \pi) \\

\phi_0 + \beta \quad (\phi \le \pi)

\end{array}

\end{equation}

\section{Numerical solution of fluid equations: general approach and ordering}

The multi-dimensional fluid equations describing the time-depend evolution of density, drift, and temperature (Equations \ref{eqn:continuitycoord}, \ref{eqn:momentumcoord}, \ref{eqn:ionenergycoord} and \ref{eqn:electronenergycoord}) are solved sequentially through the use of a time-step splitting technique. Generally speaking, each of these equations may be represented in the form:

\begin{equation}

\frac{\partial f}{\partial t} = \mathcal{L}_a \{f\} + \mathcal{L}_d \{f\} + \mathcal{L}_{ss} \{f\} + \mathcal{L}_{sn} \{f\},

\end{equation}

where $\mathcal{L}_{a,d,ss,sn}$ represent the advection, diffusion, stiff source, and non-stiff source operators, respectively. These operators are split apart and dealt with sequentially in the model (the time-step splitting, or operator splitting technique). The ordering used in GEMINI is that the advection operator is solved first.

\begin{equation}

\frac{\partial f}{\partial t} = \mathcal{L}_a \{f\}

\end{equation}

$\mathbf{v}_a$ is the advection velocity, and $\partial f / \partial t$ is an approximation of the time derivative. The partially updated value of $f$ (which has been advected) is then processed through a non-stiff source step:

\begin{equation}

\frac{\partial f}{\partial t} = \mathcal{L}_{ns} \{f\}

\end{equation}

Then a diffusion-source step is applied to the solution that has already been updated through advection and non-stiff source solutions:

\begin{equation}

\frac{\partial f}{\partial t} = \mathcal{L}_d \{f\}

\end{equation}

Typically the time scales associated with diffusion and stiff source processes require that a implicit scheme (based on backward differences) be used for efficiency. The advection terms are treated with explicit finite-volume methods since they only need to resolve time scales dictated by the ion sound speed. The last split step is a solution for the stiff source terms (which is applied to the solution which has been partially updated through the advection, non-stiff source, and diffusion operations):

\begin{equation}

\frac{\partial f}{\partial t} = \mathcal{L}_{ss} \{f\}

\end{equation}

This sequence of substeps approximates the solution to the full transport equation. Note that only the energy equations include all of these substeps; the continuity and momentum equations are solved in two steps: an advection substep and then a stiff sources substep.

Along with the time-step splitting one must choose the order in which to resolve the the different substeps. In principle any ordering will produce an approximate solution to the full equation, but we have found that solving the operators in order of increasing numerical stiffness (viz. the same order presented above) results in the most stable and accurate ordering. Specifically, first advection is solved, then nonstiff sources, diffusion, and, finally, stiff source terms. The solutions are organized such that the advection parts of all ion equations are solved together, followed by a solution for the non-stiff terms in each equation, followed by the diffusion terms for each equation. Finally the Energy stiff source terms are resolved, followed by the momentum equation stiff terms, and then the continuity equation stiff terms.

\subsection{Finite difference approximations used in GEMINI}

The numerical schemes for solving the potential equation and diffusion substep of the fluid equation solutions all rely on finite difference approximations to partial derivatives. This type of approximation is standard, but is modified slightly to deal with non-uniform grids required for efficient solution of the fluid equations.

A centered first derivative approximation at some point $x$ with neighboring grid points at $x-\Delta x_1$ and $x+\Delta x_2$ may be found from the expansions of some function $f$:

\begin{eqnarray}

f(x-\Delta x_1) &=& f(x) - f'(x)\Delta x_1 + f''(x) \frac{\Delta x_1^2}{2!} - f'''(x) \frac{\Delta x_1^3}{3!} + \dots \nonumber \\

f(x+\Delta x_2) &=& f(x) + f'(x)\Delta x_2 + f''(x) \frac{\Delta x_2^2}{2!} + f'''(x) \frac{\Delta x_2^3}{3!} + \dots \nonumber

\end{eqnarray}

Subtracting the first equation from the second gives an approximate first derivative.

\begin{equation}

f'(x) = \frac{f(x+\Delta x_2) - f(x-\Delta x_1)}{\Delta x_1+\Delta x_2} + \mathcal{O}(\Delta x_1-\Delta x_2) \label{1stderiv}

\end{equation}

Note that this is roughly second order accurate in space provided that $\Delta x_1 \approx \Delta x_2$, i.e. provided that the grid spacing does not change rapidly from one point to the next. In grid index form:

\begin{equation}

\left[ \frac{\partial f}{\partial x} \right]_i \approx \frac{f_{i+1} - f_{i-1}}{x_{i+1} - x_{i-1}}

\end{equation}

A second derivative approximation is obtained by iteratively applying Equation \ref{1stderiv}.

\begin{equation}

\left[ \frac{\partial^2 f}{\partial x^2} \right]_i \approx \frac{\left( \left[ \frac{\partial f}{\partial x} \right]_{i+1/2} - \left[ \frac{\partial f}{\partial x} \right]_{i-1/2} \right)}{x_{i+1/2} - x_{i-1/2}}

\end{equation}

which makes use of the grid midpoints $x\pm1/2$. We define, for use hereafter in this document, the forward, backward, and centered difference in $x$ by:

\begin{eqnarray}

\Delta x_{i,f} &=& x_{i+1} - x_i \nonumber \\

\Delta x_{i,b} &=& x_i - x_{i-1} \nonumber \\

\Delta x_{i,1/2} &=& x_{i+1/2} - x_{i-1/2} \nonumber

\end{eqnarray}

Using these definitions a finite difference form of the first derivative is:

\begin{equation}

\left[ \frac{\partial f}{\partial x} \right]_i \approx \frac{f_{i+1} - f_{i-1}}{\Delta x_{i,b}+\Delta x_{i,f}}

\end{equation}

The second derivative is:

\begin{equation}

\left[ \frac{\partial^2 f}{\partial x^2} \right]_i \approx \frac{f_{i+1} - f_i}{\Delta x_{i,f}\Delta x_{i,1/2}} - \frac{f_i - f_{i-1}}{\Delta x_{i,b}\Delta x_{i,1/2}}

\end{equation}

Cross partial derivative terms may be differenced as:

\begin{eqnarray}

\left[ \frac{\partial^2 f}{\partial x \partial y} \right]_{i,j} &\approx& \frac{1}{\Delta x_{i,b} + \Delta x_{i,f}}\left( \left[ \frac{\partial f}{\partial y} \right]_{i+1,j} - \left[ \frac{\partial f}{\partial y} \right]_{i-1,j} \right) \\

&\approx& \frac{1}{\Delta x_{i,b} + \Delta x_{i,f}} \frac{1}{\Delta y_{j,b} + \Delta y_{j,f}} \left( \left( f_{i+1,j+1} - f_{i+1,j-1} \right) - \left( f_{i-1,j+1} - f_{i-1,j-1} \right) \right)

\end{eqnarray}

where the indices $i$ and $j$ correspond to the $x_1-$ and $x_2-$ dimensions, respectively. Note that this differencing commutes, as the cross-partials should, e.g.:

\begin{equation}

\left[ \frac{\partial^2 f}{\partial x \partial y} \right]_{i,j} = \left[ \frac{\partial^2 f}{\partial y \partial x} \right]_{i,j}

\end{equation}

In several places in the model (e.g. in potential equation and heat conduction) it is necessary to difference a ``compound'' second order derivative term of the form:

\begin{equation}

\frac{\partial}{\partial x} \left( A \frac{\partial \Phi}{\partial x} \right) \label{eqn:compound}

\end{equation}

Two alternative finite difference approximations to this term may be derived by: (1) expanding the derivative using the product rule and then generated a finite-difference approximation using the above expressions, and (2) by differencing the equation directly. These approaches do not yield the same finite difference approximations. GEMINI uses the latter method for differencing the diffusion part of the energy equaiton as it has been found to produce slightly more stable results under conditions of very strong ion and electron heating. The FDE for this approach can be generated by sequentially applying first order differences to the expression, i.e:

\begin{equation}

\left[ \frac{\partial}{\partial x} \left( A \frac{\partial \Phi}{\partial x} \right) \right]_i \approx \frac{1}{\Delta x_{i,1/2}} \left( \left[ A \frac{\partial \Phi}{\partial x} \right]_{i+1/2} - \left[ A \frac{\partial \Phi}{\partial x} \right]_{i-1/2} \right)

\end{equation}

The derivative terms that appear in this expression can then be differenced as:

\begin{equation}

\left[ A \frac{\partial \Phi}{\partial x} \right]_{i+1/2} = \frac{A_{i+1/2}}{\Delta x_{i,f}} \left( \Phi_{i+1} - \Phi_{i} \right)

\end{equation}

where:

\begin{equation}

A_{i+1/2} \approx \frac{1}{2} \left( A_{i+1} + A_{i} \right)

\end{equation}

The final expression for the compound derivative of Equation \ref{eqn:compound} is, then, given by:

\begin{equation}

\left[ \frac{\partial}{\partial x} \left( A \frac{\partial \Phi}{\partial x} \right) \right]_i \approx \frac{1}{\Delta x_{i,1/2}} \left( \frac{A_{i+1/2}}{\Delta x_{i,f}} \left( \Phi_{i+1} - \Phi_{i} \right) - \frac{A_{i-1/2}}{\Delta x_{i,b}} \left( \Phi_{i} - \Phi_{i-1} \right) \right)

\label{eqn:compound2}

\end{equation}

Equation \ref{eqn:compound2} is used in differencing the potential and heat flow equations solved in GEMINI.

%An alternative to Equation \ref{eqn:compound2} this can be constructed by expanding the derivative prior to differencing:

%\begin{equation}

%\left[ \frac{\partial}{\partial x} \left( A \frac{\partial \Phi}{\partial x} \right) \right]_i = \left[ \frac{\partial A}{\partial x}\frac{\partial \Phi}{\partial x} + A \frac{\partial^2 \Phi}{\partial x^2} \right]_i \approx \left( \frac{A_{i+1}-A_{i-1}}{\Delta x_{i,b} + \Delta x_{i,f}} \right) \left( \frac{\Phi_{i+1}-\Phi_{i-1}}{\Delta x_{i,b} + \Delta x_{i,f}} \right) + A_i \left( \frac{\Phi_{i+1} - \Phi_i}{\Delta x_{i,f}\Delta x_{i,1/2}} - \frac{\Phi_i - \Phi_{i-1}}{\Delta x_{i,b}\Delta x_{i,1/2}} \right)

%\end{equation}

%%This form of differencing is used in the potential equation as it is most similar to the differencing scheme used for cross partial derivatives that appear in this equation.

%A third possibility would be:

%\begin{equation}

%\left[ \frac{\partial}{\partial x} \left( A \frac{\partial \Phi}{\partial x} \right) \right]_i \approx \frac{1}{\Delta x_{i,1/2}} \left( \left[ A \frac{\partial \Phi}{\partial x} \right]_{i+1} - \left[ A \frac{\partial \Phi}{\partial x} \right]_{i-1} \right) \approx ...

%\end{equation}

%but this form is not used at all in GEMINI.

%Compound derivatives with cross partial terms are slightly more complicated. For example a term of the form:

%\begin{equation}

%\frac{\partial}{\partial x} \left( A \frac{\partial \Phi}{\partial y} \right) \label{eqn:compoundcross}

%\end{equation}

%can be differenced as:

%\begin{eqnarray}

%\left[ \frac{\partial}{\partial x} \left( A \frac{\partial \Phi}{\partial y} \right) \right]_{i,j} &=& \left[ \frac{\partial A}{\partial x} \frac{\partial \Phi}{\partial y} + A \frac{\partial^2 \Phi}{\partial x \partial y} \right]_i \approx \left( \frac{A_{i+1,j}-A_{i-1,j}}{\Delta x_{i,b} + \Delta x_{i,f}} \right) \left( \frac{\Phi_{i,j+1}-\Phi_{i,j-1}}{\Delta y_{j,b} + \Delta y_{j,f}} \right) \nonumber \\ &+& \frac{A_{i,j}}{\left( \Delta x_{i,b} + \Delta x_{i,f} \right) \left( \Delta y_{j,b} + \Delta y_{j,f} \right)} \left( \left( f_{i+1,j+1} - f_{i+1,j-1} \right) - \left( f_{i-1,j+1} - f_{i-1,j-1} \right) \right) \label{eqn:compoundcross2}

%\end{eqnarray}

%Note that this operation has effectively a six point stencil. Other differencing possibilities exist for compound cross derivative, but GEMINI uses the form of Equation \ref{eqn:compoundcross2}. In the potential equation (where this term appears) we take the convention of expanding out all derivatives with the product rule (for consistency) \emph{before} generating the FDE.

\subsection{Mathematical structure of advection terms}

The advection substep for the solution of Equations \ref{eqn:continuitycoord}, \ref{eqn:momentumcoord}, and \ref{eqn:ionenergycoord} consists of three, 2D or 3D advection equations (3D is used here):

\begin{eqnarray}

\frac{\partial \rho_s}{\partial t} + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( h_2 h_3 \rho_s v_{s1} \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_2} \left( h_1 h_3 \rho_s v_{s2} \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_3} \left( h_1 h_2 \rho_s v_{s3} \right) &=& 0 \nonumber \\

\frac{\partial}{\partial t} \left( \rho_s v_{s1} \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( h_2 h_3 \rho_s v_{s1}^2 \right) + \frac{1}{h_1^2 h_2 h_3} \frac{\partial}{\partial x_2} \left( h_1^2 h_3 \rho_s v_{s1} v_{s2} \right) + \frac{1}{h_1^2 h_2 h_3} \frac{\partial}{\partial x_3} \left( h_1^2 h_2 \rho_s v_{s1} v_{s3} \right) &=& 0 \nonumber \\

\frac{\partial}{\partial t} \left( \rho_s \epsilon_s \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( h_2 h_3 \rho_s \epsilon_s v_{s1} \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_2} \left( h_1 h_3 \rho_s \epsilon_s v_{s2} \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_3} \left( h_1 h_2 \rho_s \epsilon_s v_{s3} \right) &=& 0 \nonumber

\end{eqnarray}

The geometric terms that appear in the momentum equation (which are from the tensor derivative) are lumped in with the stiff source terms and not considered here. The solution to these equations is advanced from time $t^n$ to time $t^{n+1}$ by adopting the previous time step values of drift as the set velocity with which the state variables $\rho_s,\rho_s v_{s1}$, and $\rho_s \epsilon_s$ are advected. This advection velocity is denoted $\mathbf{u}$ and defined by:

\begin{equation}

\mathbf{u} \equiv [\mathbf{v}_s]^n

\end{equation}

Hence, for a known advection velocity $\mathbf{u}$, our system of advection equations reads:

\begin{eqnarray}

\frac{\partial \rho_s}{\partial t} + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( h_2 h_3 \rho_s u_1 \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_2} \left( h_1 h_3 \rho_s u_2 \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_3} \left( h_1 h_2 \rho_s u_3 \right) &=& 0 \nonumber \\

\frac{\partial}{\partial t} \left( \rho_s v_{s1} \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( h_2 h_3 \rho_s v_{s1} u_1 \right) + \frac{1}{h_1^2 h_2 h_3} \frac{\partial}{\partial x_2} \left( h_1^2 h_3 \rho_s v_{s1} u_2 \right) + \frac{1}{h_1^2 h_2 h_3} \frac{\partial}{\partial x_3} \left( h_1^2 h_2 \rho_s v_{s1} u_3 \right) &=& 0 \nonumber \\

\frac{\partial}{\partial t} \left( \rho_s \epsilon_s \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( h_2 h_3 \rho_s \epsilon_s u_1 \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_2} \left( h_1 h_3 \rho_s \epsilon_s u_2 \right) + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_3} \left( h_1 h_2 \rho_s \epsilon_s u_3 \right) &=& 0 \nonumber

\end{eqnarray}

Which can be solved for partially updated values of the state variables $\rho_s,\rho_s v_{s1}$, and $\rho_s \epsilon_s$ (these still need to be further updated with solutions of the other operators).

\subsection{Numerical solution of the advection terms}

A finite-volume method is used to solve the advection equation listed above. This method is based on the high-resolution methods commonly used on modern computational fluid dynamics applications. The update formulas for an advected quantity $f$ in slope-limited form, for a one-dimensional advection problem on a \emph{uniform grid} are given by \citet[][ pg. 113]{Leveque:2002}:

\begin{equation}

f_i^{n+1} = f_i^{n} + \frac{\Delta t}{\Delta z} \left( \varphi_{i+1/2}^n - \varphi_{i-1/2}^n \right)

\end{equation}

This formula advances, by one time step, the solution of the equation:

\begin{equation}

\frac{\partial f}{\partial t} + \frac{\partial}{\partial z} \left( f u \right) = 0

\end{equation}

The fluxes, in slope-limited form, are given by (see Leveque's book for a derivation):

\begin{equation}

\varphi_{i-1/2}^{n+1/2} = \left\{ \begin{array}{c}

u_{i-1/2} f_{i}^n - \frac{1}{2} u_{i-1/2} \left( \Delta z + u_{i-1/2} \Delta t \right) \sigma_i^n \quad (u_{i-1/2} < 0) \\

u_{i-1/2} f_{i-1}^n + \frac{1}{2} u_{i-1/2} \left( \Delta z - u_{i-1/2} \Delta t \right) \sigma_{i-1}^n \quad (u_{i-1/2} \ge 0)

\end{array} \right.

\end{equation}

where $\sigma$ is an approximation to the slope of the function $f(z)$, and $u$ is the advection velocity. For a nonuniform grid the update formulas can be modified as:

\begin{equation}

f_i^{n+1} = f_i^{n} + \frac{\Delta t}{\Delta z_{i,1/2}} \left( \varphi_{i+1/2}^{n+1/2} - \varphi_{i-1/2}^{n+1/2} \right)

\end{equation}

The flux formulas become (THIS NEEDS TO BE CHECKED - IT'S NOT AT ALL CLEAR HOW I ARRIVE AT THESE FROM MY NOTES):

\begin{equation}

\varphi_{i-1/2} = \left\{ \begin{array}{c}

u_{i-1/2} f_{i}^n - \frac{1}{2} u_{i-1/2} \left( z_i - z_{i-1/2} + u_{i-1/2} \Delta t \right) \sigma_i^n \quad (u_{i-1/2} < 0) \\

u_{i-1/2} f_{i-1}^n + \frac{1}{2} u_{i-1/2} \left( z_{i-1/2} - z_{i-1} - u_{i-1/2} \Delta t \right) \sigma_{i-1}^n \quad (u_{i-1/2} \ge 0)

\end{array} \right.

\end{equation}

These can also be written in terms of backward differences:

\begin{equation}

\varphi_{i-1/2} = \left\{ \begin{array}{c}

u_{i-1/2} f_{i}^n - \frac{1}{2} u_{i-1/2} \left( \Delta z_{i,b} + u_{i-1/2} \Delta t \right) \sigma_i^n \quad (u_{i-1/2} < 0) \\

u_{i-1/2} f_{i-1}^n + \frac{1}{2} u_{i-1/2} \left( \Delta z_{i,b} - u_{i-1/2} \Delta t \right) \sigma_{i-1}^n \quad (u_{i-1/2} \ge 0)

\end{array} \right. \label{eqn:fluxes}

\end{equation}

In generalized curvlinear coordinates we need to solve a slightly different form of advection equation, e.g a one-dimensional curvilinear advection equation might read:

\begin{equation}

\frac{\partial \rho_s}{\partial t} + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( h_2 h_3 \rho_s u_1 \right) = 0

\end{equation}

Note that the form will be different depending on what coordinate we are advecting with respect to, as well as what type of quantity is being advected (e.g. a scalar vs. rank-1 tensor). Nonetheless the curvilinear advection equation will alway be of the form:

\begin{equation}

\frac{\partial f}{\partial t} + \frac{1}{h_\alpha^{(1)}} \frac{\partial}{\partial x_\alpha} \left( h_\alpha^{(2)} f u \right) = 0

\end{equation}

for advection with respect to the $\alpha$ coordinate (denoted $x_\alpha$). Differencing this general equation yields the update formula:

\begin{equation}

f^{n+1}_i = f^n_i + \frac{\Delta t}{ \left[ h_\alpha^{(1)} \right]_i ~ \Delta z_{i,1/2}} \left( \left[ h_\alpha^{(2)} \right]_{i+1/2} \varphi_{i+1/2}^n - \left[ h_\alpha^{(2)} \right]_{i-1/2} \varphi_{i-1/2}^n \right)

\end{equation}

In addition to this alteration the slope $\sigma$, which is measured with respect to the coordinate of advection needs to also be corrected. As computed in the model, $\sigma$ has units of [a/b], where a is the units of the variable being advected and b is the units of the coordinate being advected with respect to. In order for the flux values of Equation \ref{eqn:fluxes} to have the correct units, the velocity which multiplies the slope must be scaled using the metric factor of the dimension being advected, i.e.:

\begin{equation}

\varphi_{i-1/2} = \left\{ \begin{array}{c}

u_{i-1/2} f_{i}^n - \frac{1}{2} u_{i-1/2} \left( \Delta x_{1,i,b} + \frac{u_{1,i-1/2}}{\left[ h_1 \right]_{i-1/2}} \Delta t \right) \sigma_i^n \quad (u_{i-1/2} < 0) \\

u_{i-1/2} f_{i-1}^n + \frac{1}{2} u_{i-1/2} \left( \Delta x_{1,i,b} - \frac{u_{1,i-1/2}}{\left[ h_1 \right]_{i-1/2}} \Delta t \right) \sigma_{i-1}^n \quad (u_{i-1/2} \ge 0)

\end{array} \right. \label{eqn:fluxesfinal}

\end{equation}

Note that in Equation \ref{eqn:fluxesfinal}, the quantity $\frac{u_{1,i-1/2}}{\left[ h_1 \right]_{i-1/2} } \Delta t$ has the same units as $\Delta x_{1,i,b}$, as necessary to have a sensible update formula.

Advection in multiple dimensions is achieved by dimensional splitting, i.e. as a sequence of 1D advection problems. For example the continuity advection substep is composed of the following subsubsteps:

\begin{equation}

\frac{\partial \rho_s}{\partial t} + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_1} \left( h_2 h_3 \rho_s u_1 \right) = 0

\end{equation}

\begin{equation}

\frac{\partial \rho_s}{\partial t} + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_2} \left( h_1 h_3 \rho_s u_2 \right) = 0

\end{equation}

\begin{equation}

\frac{\partial \rho_s}{\partial t} + \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial x_3} \left( h_1 h_2 \rho_s u_3 \right) = 0

\end{equation}

where each substep update is applied using the previous step results as input.

\subsubsection{Boundary conditions for advection equations}

The boundary conditions for advection equations are implemented through the use of ghost cells \citep{Leveque:2002}. Generally, values of the plasma fluid state variables in the ghost cells are set to those on the corresponding outermost point of the domain, i.e. free-flow boundaries. Some exceptions to this exist at the top boundary ghost cells, which are often set to prevent excessive inflow into the simulation domain (see source code for details).