Simulation of high-energy radiation belt electron fluxes using NARMAX-VERB coupled codes

This study presents a fusion of data-driven and physics-driven methodologies of energetic electron flux forecasting in the outer radiation belt. Data-driven NARMAX (Nonlinear AutoRegressive Moving Averages with eXogenous inputs) model predictions for geosynchronous orbit fluxes have been used as an outer boundary condition to drive the physics-based Versatile Electron Radiation Belt (VERB) code, to simulate energetic electron fluxes in the outer radiation belt environment. The coupled system has been tested for three extended time periods totalling several weeks of observations. The time periods involved periods of quiet, moderate, and strong geomagnetic activity and captured a range of dynamics typical of the radiation belts. The model has successfully simulated energetic electron fluxes for various magnetospheric conditions. Physical mechanisms that may be responsible for the discrepancies between the model results and observations are discussed.


Introduction
The Van Allen radiation belts are composed of energetic charged particles that are trapped by the Earth's magnetic field. The electron belts exhibit a two-zone structure, with the inner and outer belts separated by an area of reduced flux known as the slot region. The relatively stable inner belt is located below 2 Earth radii, while the highly dynamic electrons of the outer belt can be found at distances greater than ≈ 3 R E . These energetic particles can affect geostationary satellites and damage spacecraft electronics [Baker, 1996]. As such, much attention has been devoted to developing accurate space weather prediction models to mitigate hazards caused by varying geomagnetic conditions. Particles trapped in the radiation belts undergo periodic motion described by adiabatic invariants. Interactions that occur on timescales faster than the periodic motion in question lead to violations of the corresponding adiabatic invariant and to diffusion of the particle in energy or pitch angle. This can happen through wave-particle interactions, and the effect of these violations on the bulk particle distribution is described by quasi-linear theory [Vedenov et al., 1961;Kennel and Engelmann, 1966]. Friedel et al. [2002] identified several mechanisms for the buildup of relativistic electrons in radiation belts. The most significant ones are believed to be radial diffusion and local acceleration [e.g., Horne and Thorne, 1998;Summers et al., 1998;Elkington et al., 1999;Green and Kivelson, 2004]. Radial diffusion is caused by gradients in the phase space density (PSD) of particles and can be enhanced by ultralow frequency (ULF) waves violating the third adiabatic invariant. Particles will diffuse radially inward or outward depending on the sign of the PSD gradient.
Local acceleration is caused by resonant wave-particle interactions with collective plasma instabilities found in the magnetosphere. These can change a particle's pitch angle and energy via violation of the first and second invariants. In particular, energy diffusion by chorus waves creates a peak in phase space density at L ≈ 4 [e.g., Summers et al., 1998;Brautigam and Albert, 2000;O'Brien et al., 2003;Chen et al., 2006;Iles et al., 2006;Shprits et al., 2007;Xiao et al., 2009aXiao et al., , 2010Reeves et al., 2013;Thorne et al., 2013a;Xiao et al., 2014]. It is postulated that the accelerated electrons can then be transported via radial diffusion to geostationary orbit, providing the bulk of high-energy (>2 MeV) particles encountered at these radial distances [e.g., Varotsou et al., 2005;Shprits et al., 2009]. Meanwhile, chorus, hiss, and electromagnetic ion cyclotron (EMIC) waves contribute to pitch angle scattering of trapped particles and their demise from the system via the loss PAKHOTIN ET AL.

10.1002/2014JA020238
cone [Xiao et al., 2009b[Xiao et al., , 2012. Another loss process is magnetopause shadowing [e.g., Desorgher et al., 2000], where drifting particles come into contact with the magnetopause on the dayside and escape to interplanetary space. The loss to the magnetopause can cause outward radial diffusion and extend losses down to the heart of the outer belt (below L ≈ 4) [Shprits et al., 2006a;Turner et al., 2012].
The physics behind radiation belt storm dynamics are not yet fully understood, and physics-based models, in particular, the wave models, are in a constant state of development. Additionally, they need a measure of the electron population entering the inner magnetosphere from the tail [e.g., Ganushkina et al., 2013Ganushkina et al., , 2014. Such a quantity is quite hard to measure empirically. Therefore, at present some of the most accurate models of radiation belt flux forecasting are data driven [e.g., Baker et al., 1990;Li et al., 2001;Kellerman et al., 2013]. In particular, the NARMAX (Nonlinear AutoRegressive Moving Averages with eXogenous inputs) model uses system identification to deduce a mathematical model from recorded data sets by identifying important model terms [Boynton et al., 2013a]. In that sense it is similar to the nonlinear equivalent of a correlation function. The input data is integral high-energy electron fluxes at geostationary orbit using data from the GOES-13 satellite, as well as solar wind parameters from L1. However, data-driven models need continuous data to generate predictions, and the NARMAX approach using geosynchronous orbit (GSO) data gives no indication of what happens at different radial distances.
This study presents a system that merges the physical and data-driven approaches, combining the accuracy of the data-driven NARMAX algorithm with the predictive range and resolution of the physics-based VERB code [Shprits et al., 2008a;Subbotin and Shprits, 2009;Shprits et al., 2009]. The VERB code is used to produce an electron phase space density profile as a function of energy, pitch angle, time, and L * , where L * refers to the dimensionless Roederer L * [Roederer, 1970] and is related to the third adiabatic invariant through where M is the Earth's dipole magnetic moment, Φ is magnetic flux enclosed by the guiding electron drift shell, and R E is radius of Earth. The VERB code takes two inputs: the geomagnetic Kp index and a measure of electrons entering the inner magnetosphere from the tail. While the Kp index is in public access and can be forecast with existing models, seed electron fluxes can be quite difficult to estimate. Subbotin et al. [2011a] used outputs from the Rice Convection Model (RCM) [Toffoletto et al., 2003, and references therein] as an input to the outer boundary. Combining the data-driven NARMAX model with the VERB code is an alternate solution and is employed in this study. The resulting scheme is called VNC-VERB-NARMAX Coupling-and generates accurate predictions of the entire outer belt region for a range of magnetospheric conditions. The rest of the paper is split into several sections. Section 2 summarizes the satellite instrumentation and orbital parameters, section 3 provides detailed descriptions of both models, and the coupling method is explained in section 4. Section 5 presents the results and section 6 summarizes the conclusions and caveats and outlines future work. f here, f is the electron phase space density, t is time, p is relativistic momentum, 0 is equatorial pitch angle, L is McIlwain's magnetic shell parameter (radial distance from the center of the Earth to the equatorial field line point) and L * is the Roederer parameter.
The , J, and Φ are the three adiabatic invariants that may be expressed as (1) is a function related to the bounce frequency and is approximated by Lenchek et al. [1961] T( 0 ) = 1.3802 − 0.3198(sin 0 + sin 1∕2 0 ) Radial diffusion is represented by the first term on the right-hand side of the equation, while diffusion in energy (momentum), pitch angle, and mixed terms are represented by the next four. Losses to the PAKHOTIN ET AL.
©2014. The Authors. atmosphere are represented by the last term f ∕ where is a characteristic electron lifetime assumed to be infinite outside the loss cone and equal to a quarter bounce time inside it. All pitch angle and energy diffusion coefficients are parameterized by Kp index either explicitly or implicitly through plasmapause location. They are bounce-and MLT-averaged and data on wave-particle interactions are drawn from statistical studies.
The first term in (1) is the radial diffusion term written in terms of , J, and L * and describes the radial diffusion of PSD with the radial diffusion coefficient D L * L * . Radial diffusion violates the third invariant while conserving the first and second. Meanwhile, wave-particle interactions violate the first and second invariants, leading to diffusion in energy and pitch angle. Radial diffusion is assumed to be driven by resonant interactions with ULF waves, with diffusion coefficient due to the magnetic component of ULF adopted from Brautigam and Albert [2000]: The radial diffusion coefficient was parameterized for the Kp index below 6, however, in the VERB model the same parameterization is used for all Kp indices.
The last four terms describe diffusion in momentum (D pp ), pitch angle (D 0 0 ), and mixed diffusion terms (D p 0 = D 0 p ). Energy and relativistic momentum are related via E = √ p 2 c 2 + m 2 c 4 − mc 2 where c is the speed of light.
The waves used in the simulation are dayside chorus, nightside chorus, and plasmaspheric hiss. The assumed wave model for the chorus waves is reproduced on Table 1 from Subbotin et al. [2011b]. Plasmaspheric hiss data are used from Orlova et al. [2014] where the wave normal angle distribution was taken from Table S3 in the supporting information of Thorne et al. [2013b]. The plasmapause location is calculated using Carpenter and Anderson [1992]: where Kp max24 is the maximum Kp values in the preceding 24 h. The diffusion coefficients are parameterized by the Kp index. VERB assumes a dipole field geometry for the computation and bounce averaging of pitch angle and energy diffusion coefficients.
VERB solves equation (1) numerically for L * from 1 to 7 where 7 is assumed to be the outer simulation boundary. An implicit finite differences method is used on an orthogonal grid for radial, energy, and equatorial pitch angle diffusion. Six boundary conditions are specified, two for each variable in the equation. The two pitch angle boundary conditions are as follows: zero PSD gradient at 0 • and 90 • . The PSD for low-energy electrons is constant to represent a balance of convective source and losses, while the PSD for high-energy electrons is zero (very few electrons are found at very high energies). Finally, there are assumed to be no PAKHOTIN ET AL.
electrons present at the low L * boundary to represent atmospheric losses. The high L * boundary is time dependent and is based on CRRES data. Following Brautigam and Albert [2000], an energy-independent flux B f (t) is used to reproduce the variation in convected outer boundary electrons. It is this parameter that requires NARMAX. Fluxes predicted by NARMAX at GSO are translated into expected differential fluxes at the outer boundary using a procedure outlined in the next section. These are translated from differential fluxes in standard units to fluxes as a function of B f (t) and input into the VERB code. This, along with the Kp index, provides all the necessary inputs to run the simulation.

Coupled Simulation
The VERB code was initialized with Kp index values from the German Research Centre for Geosciences Web site (http://www-app3.gfz-potsdam.de/kp_index/). The boundary flux-the second input to the system-was taken from NARMAX predictions for fluxes of > 800 keV and > 2 MeV. Differential flux was calculated from these two integral channels assuming exponential energy distributions for an energy of 1 MeV. The reason for this choice was to avoid extraneous dynamics that begin to predominate at higher energies . Because NARMAX predicts for geostationary orbit, it was necessary to recalculate fluxes for the outer boundary at L * = 7. A representative value of the average L * location of geostationary orbit was taken at L * = 6.2. Fluxes were recalculated using the assumptions that PSD remains constant at fixed first invariant between L * = 6.2 and L * = 7. The recalculated fluxes formed the outer boundary condition for VERB. A profile of particle flux as a function of time and L * was generated for an energy of 892 keV so as to compare it directly to output from the Van Allen probes mission MagEIS particle detector channel that has the same energy. This particular channel was chosen because it registers relatively high-energy particles that are nevertheless below 1 MeV.
From the profile generated by VERB, the predicted flux was extracted for each hour at geostationary orbit. Geostationary orbit location in L * is calculated using the ONERA library and the Kp-driven T89 magnetic field model [Tsyganenko, 1989]. The differential flux taken from the model at 892 keV and 1000 keV was used to calculate flux from the integral energy channel > 800 keV to compare with GOES-13 data.

Results
Electron flux observations from the GOES-13 and Van Allen probe spacecraft of electron energies above 800 keV were compared with fluxes predicted by the coupled model. This was done for several time periods. The first event, from 1 to 25 August 2013, featured quiet periods punctuated by extended periods of relatively high geomagnetic activity with storm sudden commencements (SSCs). The second time period, 1 to 23 September 2013, is largely geomagnetically quiet with Kp index never rising above 3+. The third time period, 1 to 20 October 2013, featured very rapid changes of Kp indices and instances of high geomagnetic activity. Figure 1 shows, for the time period 1-25 August 2013 (a) the Dst and Kp indices, (b) the differential daily averaged flux from GOES-13, NARMAX, and VNC at geostationary orbit, (c) data from the RBSP-B probe for an energy channel corresponding to 892 keV, and (d) the VERB simulation with a NARMAX-driven outer boundary. It can be seen that there are three periods of electron flux enhancements: 1-4 August, 5-16 August, and 17-26 August with a small decrease in flux around 22-23 August. All these are reproduced by NARMAX (Figure 1b) and as a consequence reflect conditions on the outer boundary in the VNC simulation ( Figure 1d). The intensity of dropouts is somewhat underestimated by the model, while fluxes on 4 August exceeded NARMAX predictions and, consequently, VNC reported smaller fluxes than were measured by the Van Allen probe. Otherwise, the model successfully captured the dynamics of this time period, with all three periods of activity visible on Figure 1d.

First Period: 1-25 August 2013
The reason for the underestimation of fluxes on 4 August can be deduced from observing the input to the outer boundary (Figure 1b). On 5 August, increased flux was observed on GOES-13 (red) but not NAR-MAX (blue) which in this case registered the increase with a delay. The injection can be seen on the model (Figure 1d) around 6 August; the Kp index drops at this point. Since acceleration wave models in VERB are parameterized by the Kp index, the magnitude of these processes was too low to reproduce the quite intense flux profile seen on MagEIS data for that time. Thus the variation in this case is a consequence of the fact that VNC is driven not by real GOES-13 data but by the NARMAX predictions of it. Allen probe B MagEIS data for 892 keV, and (d) VERB simulation of the outer radiation belt using the NARMAX-derived differential flux as the outer boundary. The black curve in Figure 1d represents geostationary orbit.

10.1002/2014JA020238
The GOES-13 geosynchronous satellite orbit is represented on Figure 1d as a black sinusoidal curve between L * = 5 − 7. The location of GSO in L * space varies throughout the day. Figure 2 shows the same information as Figure 1 for the time period 1-23 September 2013. It can be seen here that under these somewhat more quiet conditions, energetic electron fluxes are lower and have two main enhancements, 1-11 September and 19-24 September. In particular, the time periods 1-11 September features a long, steady period of high fluxes around L * = 5, which subsequently falls off by 13 September.

Second Period: 1-23 September 2013
There is another small enhancement barely visible on Figure 2c on 15-16 September on the edge of the satellite's orbital coverage. It can be seen more clearly in Figure 2b from both GOES and NARMAX, suggesting that the particles came from the outer boundary through geostationary orbit. The long buildup of electron flux can be seen on Figure

Third Period: 1-20 October 2013
This period features three strong flux enhancement periods that can be seen on Figure 3: 2-9 October, 9-15 October, and 15-20 October. Sharp dropouts are evident between these time periods (Figure 3c) that can be characterized as main phase depletions from the Dst index on Figure 3a. The three events here appear to be caused by strong injections from the outer boundary, captured by GOES-13 at GSO and predicted successfully by NARMAX (Figure 3b). Meanwhile, the profile of the Kp index is unusual in that it contains periods of very low geomagnetic activity, then rapid rise to high indices (as high as Kp = 7 on 2 October) and shortly afterward, a drop back to very low activity.
Under these conditions, the simulation results begin to deviate from the observations. A small enhancement (1-2 October) from above L * ≈ 5.5 down to L * = 4.5 is observed on both the model and observations, followed by a sharp dropout on 2 October. However, the enhancement afterward (2-9 October), although visible on the model, is over an order of magnitude smaller than in observations. There are several possible reasons for this discrepancy. First, NARMAX underestimates GSO flux on 2-3 October (as is evident by the second blue and red data points on Figure 3b). As a result, flux at the outer boundary is underestimated, and this carries through into lower L shells. Second, local acceleration and radial diffusion rapidly drop off in the simulation after 3 September because they are parameterized by the Kp index, which drops to near zero. Since the convected population is underestimated, and little local acceleration occurs, the result is an underestimate in fluxes across the whole outer belt for that time period. The peak in fluxes in the simulation around 2-5 October is also around L * = 3.5 while the observations suggest that it is around L * = 4 − 4.5. The peak in flux in the simulation is most likely as a result of electrons adiabatically accelerated because of high geomagnetic activity. The observations do show an enhancement inward of L * = 4, but electron fluxes below L * = 3.5 fall off rapidly, likely due to scattering by plasmaspheric hiss. It is possible that the model underestimates the effect of plasmaspheric hiss during this time period, perhaps due to the unusual Kp index dynamics.
The other two enhancements (9-15 October and 15-20 October) are predicted more accurately, with fluxes peaking around L * 4 − 4.5, consistently with observations. Flux drops around 8-9 October and 14-15 October are also seen in the model and observations. In both enhancements, the magnitude of fluxes is again underestimated by the model. Meanwhile, the effect of hiss scattering below L * = 4 appears to be underestimated by the model as well, resulting in a less clear outer belt edge than observed by MagEIS.
Since the effect of all waves appears to be underestimated, it is possible to suggest that the source of the discrepancy is the wave model because of very rapid changes in Kp index. Any inaccuracy with outer boundary flux will also carry through into local acceleration, since the outer boundary provides the seed and source population for the outer belt.
©2014. The Authors. Figure 4 shows GOES-13 fluxes for the > 800 keV channel compared with fluxes predicted by VNC. Differential fluxes were obtained for two energies-891 keV and 1 MeV-and used to calculate > 800 keV integral fluxes assuming exponential energy distributions. It can be seen that, absent some underestimations at periods of rapid flux decrease, the model follows the data reasonably well. The same analysis was performed on the 1-23 September 2013 time interval, with the results visible on Figure 5. As noted above, the sinusoidal flux variation is a diurnal effect due to the spacecraft traversing different L * regions at different magnetic local times. Black crosses on the Figures 4 and 5 show NARMAX predictions for those time periods. It can be seen that NARMAX predicts variations quite well at geostationary orbit and can be successfully used for outer boundary information as an input to the VERB system. Figure 6 shows geosynchronous orbit data for 1-20 October 2013. Under the extreme geomagnetic conditions presented, the model (red) differ from the data (blue). In particular, the model estimate of the first enhancement around 3-4 October is over an order of magnitude lower than that observed. However, NARMAX predicts the enhancement effectively, as can be seen by the third black cross for 3 October. This suggests that the increased flux simply does not propagate from the outer model boundary to lower L * values. The reason for this could be the rapid drop in all diffusion coefficients as the Kp index suddenly drops on 3 October. This would mean that any increases in outer boundary flux do not diffuse below the outer boundary, and there are few particles inside the belt to accelerate. Further, as mentioned above, local acceleration and pitch angle scattering would also be reduced under low Kp indices. Finally, the choice of magnetic field model can play a part since solar wind data and Dst index are not explicitly included in the T89 model.

Geostationary Orbit
The other two enhancements are captured more accurately, though it can be seen that model fluxes on 9-15 October and 15-20 October are still somewhat lower than those observed on GOES-13. Flux depletions occur too rapidly (on the scale of a few hours) for the model to react to them.

NARMAX as an Outer Boundary Driver
This study presents a first attempt to use predictions for GOES-13 data to drive a radiation belt diffusion code. It is believed that the NARMAX segment of the model would benefit from a higher time resolution. Rapid changes in GSO flux should then be captured sooner. In an environment of rapidly changing geomagnetic activity, a situation can arise where a rapid flux buildup or drop is only registered by the model 24 h later, by which point the geomagnetic indices driving the wave model may have PAKHOTIN ET AL.
©2014. The Authors. changed completely. This propagates the error from the outer boundary into the heart of the belt. Further, at high geomagnetic activity the location of GSO in L * may drop below L * = 6.2, meaning the flux transfer operation may no longer be valid. One solution is to set the outer boundary lower. However, this may affect GSO prediction capacity, as a low boundary can take GSO out of the simulation domain.
High-energy channels (> 800 keV and > 2 MeV) are used for outer boundary flux determination. Meanwhile, it is believed that convected lower energy electrons form the source and seed populations for the radiation belt. A lower energy channel may thus give better predictions, in particular with a higher time resolution that will allow to isolate fluxes from nightside MLT sectors only. Work to this end is currently ongoing.

Magnetopause Shadowing
Magnetopause shadowing is a phenomenon caused by a pressure pulse from the solar wind, for example, a fast CME, compressing the magnetosphere. If the boundary of the magnetosphere is near geosynchronous orbit, the drifting radiation belt electrons impact the magnetopause and are lost from the system to the solar wind. This causes a rapid drop in particle counts [e.g., Shprits et al., 2006aShprits et al., , 2006bHorne et al., 2013], affecting all energies [e.g., Bortnik et al., 2006] and potentially accounting for over 90% of electron loss outward of L * = 5 [Yu et al., 2013]. Since the pressure pulse in the solar wind can only be detected at L1 and takes about 40 min to reach the Earth, it is hard to predict for such an event several hours in advance. However, since dropouts cause depletions, rather than buildups, of relativistic electrons, it will lead to an over-prediction of flux levels and should not pose any real risk. Improving the coupling between the model and the solar wind together with a model for magnetopause location [Shue et al., 1997] is a subject of current work.

Kp Index
The Kp index is used to drive both the VERB model (where it influences the location and effect of the local diffusion region) and also the T89 model for calculating the magnetic shell and constructing the diurnal variation dynamic. It is well known [e.g., Roederer, 1970] that MLT location accounts for most of the variation of energetic electron flux, and it is seen again in the observations above. Thus having an accurate measure of Kp index is crucial for the success of the prediction. To move toward forecasting, a model that predicts Kp index should be employed, such as the Wing Kp forecast model [Wing et al., 2005]. Research to this end is currently ongoing.

Computing Time Considerations
As outlined in Subbotin et al. [2010], including the full Fokker-Planck diffusion equation involves including mixed energy and pitch angle diffusion coefficients. This can place computational strain on the model and slow down the simulation. One effect of including mixed diffusion terms is an effective inhibition of local acceleration; the peak in PSD associated with local acceleration forms later and is not as large as when mixed diffusion terms are omitted. Changes to the dynamics of the source region in turn affect simulated fluxes at geostationary orbit, and it was found that the VERB-NARMAX model performance was improved by inclusion of the mixed terms. All the above simulations were performed with mixed diffusion turned on. In practical terms, if the aim is to create a system that would generate forecasts every hour, the system has to physically have time to run all the necessary calculations. This can be done by assigning the simulations as a series of batch jobs to a supercomputer or by running the simulation with a reduced grid size or for smaller time periods. Omission of the mixed diffusion terms has been found to produce a satisfactory fit to observations, although with reduced accuracy.

Summary
A model fusing the data-driven NARMAX and physics-driven VERB codes has been developed and tested. The model combines the accuracy of the data-driven approach with the predictive range and resolution of physical methodology. The coupled simulation was tested on three separate time periods, simulating 68 days in total. The aim was to test the long-term performance of the model at predicting energetic electron fluxes at and inside of geostationary orbit in conditions typically encountered in the magnetosphere. The time periods studied contain all typical features of radiation belt dynamics: SSCs, dropouts, refilling, and rapid changes in geomagnetic indices. The combined system was found to reasonably accurately model fluxes during quiet and disturbed conditions, capturing all the essential dynamics. Sources of discrepancies, as well as caveats for the use of the model in real-time radiation belt forecasting, have been outlined in the manuscript.
These simulations show that the dynamic evolution of the outer regions of the radiation belts presents a strong driver of the dynamics of the radiation belts. Kim et al. [2012] further showed that losses to the magnetopause can extend to L shells much lower than VERB typical outer boundary and may explain overestimation of fluxes during dropout events. Accurate simulations of the dropouts will include accounting for the MLT-dependent loss to the magnetopause and accurate simulation of adiabatic effects during storms [Kim et al., 2010;Saito et al., 2010;Su et al., 2011;Kim and Chan, 1997]. The model has been found to slightly overestimate energetic electron fluxes, especially during periods of rapid radiation belt depletion, but simulates radiation belt rebuilding with good accuracy. Data from the Van Allen probes showed good agreement and allowed to validate the system for the entire outer radiation belt region.
Our simulations also confirm suggestions of Kim et al. [2012] that losses to the magnetopause can extend to the heart of the inner belt and cannot be modeled simply by the variations of the outer boundary. Future work should include MLT dependent simulations that account for convection, adiabatic effects and MLT-and pitch-angle-dependent loss to the magnetopause.
The study shows that high-energy electron fluxes near the outer boundary can be used to predict radiation belt dynamics. Since convection mainly brings low-energy particles to the radiation belts, this is somewhat surprising, particularly considering the low time resolution of NARMAX. In essence, NARMAX in its current form captures all high-energy particles at GSO, including those locally accelerated inside the belt. An intriguing question emerges-what proportion of high-energy particles observed at GSO are the high-energy tail of an injected population, and what proportion are locally accelerated electrons that have then radially diffused outward to GSO? A comparative study of the effectiveness of low-energy and high-energy channels from GSO satellites would shed light on this matter. acknowledge the efforts of the Van Allen probes team in creating the data products for the MagEIS instrument and also the work done by the GOES team. Contribution toward this work from Adam Kellerman has been greatly appreciated. The work at UCLA was supported by the Lab Fee grant, NASA grants NNX10AK99G and NNX13AE34G, and NSF grant AGS-1203747.
Michael Liemohn thanks two anonymous reviewers for their assistance in evaluating this paper.