Search

Hollow cathode plasma contactors have an extensive history as spacecraft neutralizers, but questions remain on whether the technology can neutralize significant positive charging in tenuous space plasmas. Simulations and a representative semianalytical model (which we term the ion emission model) predict effective neutralization in these scenarios as ion current is emitted from the surface of the quasi-neutral contactor plasma according to the space-charge limit. This article focuses on experimental plasma measurements of a hollow cathode biased to several positive voltages with respect to a surrounding vacuum chamber in order to validate the ion emission model. Particular attention is paid to the emitted ions (those reaching the chamber wall) and to the ion-rich sheath region where the ion emission process takes place. Retarding potential analyzer (RPA) measurements were performed to understand ion flow velocity and direction as these parameters relate directly to the current predicted by the ion emission model. Planar probe measurements were performed at the chamber wall to compare local emission currents to model predictions. Evidence of collisions within the plasma (particularly charge exchange collisions) and simple models predicting ion energetics are presented. These experiments suggest that ion emission from the contactor plasma acts according to the ion emission model and adds to the physical understanding of spacecraft neutralization using hollow cathode plasma contactors as it may occur in tenuous space plasmas.

Ganguli and Crabtree have written a comment about a recent article by the authors listed above on radiation-belt remediation. They have objected to our evaluation of the Naval Research Laboratory’s chemical release concept which states that this concept may be impractical due to an apparently low overall efficiency. In their comment, they provide a scientific argument and refer to the published literature to counter our statement. Here, we provide more details on our numerical calculations and experimental results which led to this evaluation.

Energetic electrons can be trapped in Earth’s magnetic field, forming the radiation belts (also known as the Van Allen Belts). These electrons, which can originate from the solar wind or a high-altitude nuclear explosion (HANE), have the potential to damage satellites in low-Earth orbit (LEO). For example, in 1962, the U.S. detonated a nuclear device at an altitude of about 400 km in the Starfish experiment. The resulting enhancement of the radiation belts disabled several satellites within a few months and energetic electrons remained in the radiation belts for up to several years. In order to address this potential vulnerability, schemes have been proposed to drain electrons from the radiation belts, with the most promising approaches based on using high-power very-low-frequency (VLF) waves to scatter the electrons into more field-aligned trajectories, forcing them to precipitate into Earth’s atmosphere. This paper will provide an overview of enhanced electron distributions in the radiation belts as well as approaches to VLF wave belt remediation including the use of either antennas or relativistic electrons beams in space to generate the VLF waves.

When the Earth’s moon is in the supersonic solar wind, the darkside of the Moon and the lunar plasma wake can be very dangerous charging environments. In the absence of photoelectron emission (dark) and in the absence of cool plasma (wake), the emission or collection of charge to reduce electrical potentials is difficult. Unique extreme charging events may occur during impulsive solar-energetic-electron (SEE) events when the lunar wake is dominated by relativistic electrons, with the potential to charge and differentially charge objects on and above the lunar surface to very-high negative electrical potentials. In this report the geometry of the magnetic connections from the Sun to the lunar nightside are explored; these magnetic connections are the pathways for SEEs from the Sun. Rudimentary charging calculations for objects in the relativistic-electron environment of the lunar wake are performed. To enable these charging calculations, secondary-electron yields for impacts by relativistic electrons are derived. Needed lunar electrical-grounding precautions for SEE events are discussed. Calls are made 1) for future dynamic simulations of the plasma wake in the presence of time-varying SEE-event relativistic electrons and time-varying solar-wind magnetic-field orientations and 2) for future charging calculations in the relativistic-electron wake environment and on the darkside lunar surface.

Insufficiently accurate magnetic-field-line mapping between the aurora and the equatorial magnetosphere prevents us from determining the cause of many types of aurora. An important example is the longstanding question of how the magnetosphere drives low-latitude (growth-phase) auroral arcs: a large number of diverse generator mechanisms have been hypothesized but equatorial magnetospheric measurements cannot be unambiguously connected to arcs in the ionosphere, preventing the community from identifying the correct generator mechanisms. Here a mission concept is described to solve the magnetic-connection problem. From an equatorial instrumented spacecraft, a powerful energetic-electron beam is fired into the atmospheric loss cone resulting in an optical beam spot in the upper atmosphere that can be optically imaged from the ground, putting the magnetic connection of the equatorial spacecraft’s measurements into the context of the aurora. Multiple technical challenges that must be overcome for this mission concept are discussed: these include spacecraft charging, beam dynamics, beam stability, detection of the beam spot in the presence of aurora, and the safety of nearby spacecraft.

The 1970s and 1980s were heydays for using active electron beam experiments to probe some of the fundamental physical processes that occur throughout the heliosphere and in astrophysical contexts. Electron beam experiments were used to study spacecraft charging and spacecraft-plasma coupling; beam-plasma interaction physics; magnetic bounce and drift physics; auroral physics; wave generation; and military applications. While these experiments were enormously successful, they were also limited by the technologies that were available at that time. New advances in space instrumentation, data collection, and accelerator technologies enable a revolutionary new generation of active experiments using electron beams in space. In this paper we discuss such an experiment, the Beam Plasma Interactions Experiment (Beam PIE), a sounding rocket experiment designed to (a) advance high-electron mobility transistor-based radio frequency (RF) linear accelerator electron technology for space applications and (b) study the production of whistler and X-mode waves by modulated electron beams.

We discuss the development, analysis, implementation, and numerical assessment of a spectral method for the numerical simulation of the three-dimensional Vlasov–Maxwell equations. The method is based on a spectral expansion of the velocity space with the asymmetrically weighted Hermite functions . The resulting system of time-dependent nonlinear equations is discretized by the discontinuous Galerkin (DG) method in space and by the method of lines for the time integration using explicit Runge–Kutta integrators . The resulting code, called Spectral Plasma Solver (SPS-DG), is successfully applied to standard plasma physics benchmarks to demonstrate its accuracy, robustness, and parallel scalability.

Whistler wave-particle interactions play an important role in the Earth inner magnetospheric dynamics and have been the subject of numerous investigations. By running a global kinetic ring current model (RAM-SCB) in a storm event occurred on Oct 23–24 2002, we obtain the ring current electron distribution at a selected location at MLT of 9 and L of 6 where the electron distribution is composed of a warm population in the form of a partial ring in the velocity space (with energy around 15 keV) in addition to a cool population with a Maxwellian-like distribution. The warm population is likely from the injected plasma sheet electrons during substorm injections that supply fresh source to the inner magnetosphere. These electron distributions are then used as input in an implicit particle-in-cell code (iPIC3D) to study whistler-wave generation and the subsequent wave-particle interactions. We find that whistler waves are excited and propagate in the quasi-parallel direction along the background magnetic field. Several different wave modes are instantaneously generated with different growth rates and frequencies. The wave mode at the maximum growth rate has a frequency around 0.62 ω c e , which corresponds to a parallel resonant energy of 2.5 keV. Linear theory analysis of wave growth is in excellent agreement with the simulation results. These waves grow initially due to the injected warm electrons and are later damped due to cyclotron absorption by electrons whose energy is close to the resonant energy and can effectively attenuate waves. The warm electron population overall experiences net energy loss and anisotropy drop while moving along the diffusion surfaces towards regions of lower phase space density, while the cool electron population undergoes heating when the waves grow, suggesting the cross-population interactions.

We present an electrostatic Particle-In-Cell (PIC) code on multi-block, locally structured, curvilinear meshes called Curvilinear PIC (CPIC). Multi-block meshes are essential to capture complex geometries accurately and with good mesh quality, something that would not be possible with single-block structured meshes that are often used in PIC and for which CPIC was initially developed. Despite the structured nature of the individual blocks, multi-block meshes resemble unstructured meshes in a global sense and introduce several new challenges, such as the presence of discontinuities in the mesh properties and coordinate orientation changes across adjacent blocks, and polyjunction points where an arbitrary number of blocks meet. In CPIC, these challenges have been met by an approach that features: (1) a curvilinear formulation of the PIC method: each mesh block is mapped from the physical space, where the mesh is curvilinear and arbitrarily distorted, to the logical space, where the mesh is uniform and Cartesian on the unit cube; (2) a mimetic discretization of Poisson's equation suitable for multi-block meshes; and (3) a hybrid (logical-space position/physical-space velocity), asynchronous particle mover that mitigates the performance degradation created by the necessity to track particles as they move across blocks. The numerical accuracy of CPIC was verified using two standard plasma–material interaction tests, which demonstrate good agreement with the corresponding analytic solutions. Compared to PIC codes on unstructured meshes, which have also been used for their flexibility in handling complex geometries but whose performance suffers from issues associated with data locality and indirect data access patterns, PIC codes on multi-block structured meshes may offer the best compromise for capturing complex geometries while also maintaining solution accuracy and computational efficiency.

We prove the convergence of a spectral discretization of the Vlasov-Poisson system. The velocity term of the Vlasov equation is discretized using either Hermite functions on the infinite domain or Legendre polynomials on a bounded domain. The spatial term of the Vlasov and Poisson equations is discretized using periodic Fourier expansions. Boundary conditions are treated in weak form through a penalty type term, that can be applied also in the Hermite case. As a matter of fact, stability properties of the approximated scheme descend from this added term. The convergence analysis is carried out in details for the 1D-1V case, but results can be generalized to multidimensional domains, obtained as Cartesian product, in both space and velocity. The error estimates show the spectral convergence, under suitable regularity assumptions on the exact solution.

We investigate the properties of the ion sky maps produced by the Dual Ion Spectrometers (DIS) from the Fast Plasma Investigation (FPI). We have trained a convolutional neural network classifier to predict four regions crossed by the MMS on the dayside magnetosphere: solar wind, ion foreshock, magnetosheath, and magnetopause using solely DIS spectrograms. The accuracy of the classifier is >98%. We use the classifier to detect mixed plasma regions, in particular to find the bow shock regions. A similar approach can be used to identify the magnetopause crossings and reveal regions prone to magnetic reconnection. Data processing through the trained classifier is fast and efficient and thus can be used for classification for the whole MMS database., Accepted to JGR: Space Physics

Planned active space experiments and ideas for future active space experiments are reviewed. Three active experiments being readied are DSX (Demonstration and Space eXperiments), SMART (Space Measurement of Rocket-released Turbulence), and BeamPIE (Beam Plasma Interaction Experiment). Ideas for future experiments include relativistic-electron-beam experiments for magnetic-field-line tracing, relativistic-electron-beam experiments to probe the middle atmosphere, plasma-wave launching using superparamagnetic-nanoparticle amplification of magnetic fields, the heavy-ion mass loading of collisionless magnetic-field-line reconnection, the use of electrostatically charged tethers to pitch-angle scatter radiation-belt particles, cold plasma releases to modify magnetospheric plasma physics, and neutral-gas releases to enhance neutral-particle imaging of the magnetosphere. Technologies that are being developed to enable future space active experiments are reviewed: this includes the development of compact relativistic accelerators, superparamagnetic particle amplified antennae, CubeSats, and a new understanding of how to control dynamic spacecraft charging. New capabilities to use laboratory facilities to design space active experiments as well as new computer-simulation capabilities to design and understand space active experiments are reviewed.

Determining the magnetic connectivity between magnetospheric phenomena and ionospheric phenomena is an outstanding problem of magnetospheric and ionospheric physics. Accurately establishing this connectivity could answer a variety of long-standing questions. The most-viable option to solve this is by means of a high-power electron beam fired from a magnetospheric spacecraft and spotted at its magnetic footpoint in the ionosphere. This has technical difficulties. Progress has been made on mitigating the major issue of spacecraft charging. The remaining physics issues are identified, together with the need for a synergistic effort in modeling, laboratory experiments and, ultimately, testing in space. The goal of this commentary is to stimulate awareness and interest on the magnetosphere-ionosphere connectivity problem, and possibly accelerate progress towards its solution.

Based on discussions held at a workshop in Bra-Pollenzo, Italy, this paper elaborates upon 19 outstanding questions of plasma physics in the Earth's magnetosphere. The questions are grouped according to (a) driving processes, (b) radiation belt and ring current issues, (c) auroral physics, (d) internal plasma processes, and (e) magnetosphere-ionosphere mapping issues. Future needs for magnetospheric plasma physics (measurements, techniques, simulations, theories, studies) are outlined.

An active mapping mission is described that unambiguously connects measurements in the Earth's magnetosphere to visible aurora in the atmosphere. The core of the mission is an electron-beam source operated on a spacecraft in the equatorial magnetosphere, with the electron beam traveling along the Earth's magnetic-field lines to the atmosphere, depositing its energy to create an optical beam-spot in the atmosphere at the footpoint of the spacecraft's magnetic-field line. This optical spot can be imaged by ground-based cameras, putting the location of the spacecraft's magnetic footpoint into the context of the optical aurora. Scientific instruments carried on the spacecraft make critical measurements of the properties of the magnetosphere at the locations where the magnetosphere powers the aurora, allowing the determination of the plasma-physics mechanisms by which the magnetosphere drives the aurora, in particular answering the outstanding question of how the magnetosphere drives low-latitude auroral arcs. Long-standing questions in magnetosphere-ionosphere coupling that have not been answered because we could not unambiguously connect locations in the magnetosphere with their image in the ionosphere will finally be addressed. In this paper the properties of a “standard” growth-phase auroral arc are collected, theories of the magnetospheric generation of auroral arcs are reviewed, and critical magnetospheric measurements to discern the mechanisms that drive auroral arcs are determined. Further, the plasma physics of the experiment is investigated, including spacecraft-charging mitigation, beam stability, beam scattering, and electron orbit theory. Tradeoffs (keV versus MeV) concerning the energy of the electron beam are enumerated.

We study numerically small-scale reconnection events in kinetic, low-frequency, quasi-2D turbulence (termed kinetic-Alfv\'en turbulence). Using 2D particle-in-cell simulations, we demonstrate that such turbulence generates reconnection structures where the electron dynamics do not couple to the ions, similarly to the electron-only reconnection events recently detected in the Earth's magnetosheath by Phan et al. (2018). Electron-only reconnection is thus an inherent property of kinetic-Alfv\'en turbulence, where the electron current sheets have limited anisotropy and, as a result, their sizes are smaller than the ion inertial scale. The reconnection rate of such electron-only events is found to be close to $0.1$., Comment: 8 pages, 3 figures

Energetic electron beams have been proposed for tracing magnetic field lines from the magnetosphere down to the ionosphere, in active experiments aimed at diagnosing mechanisms at play in the coupling between magnetosphere and ionosphere. It is recognized however that in the absence of an efficient mitigation technique, this approach would lead to unacceptably large spacecraft charging and positive potential buildup, which would result in environmental hazard for the spacecraft. This problem would be particularly acute in low density regions of the magnetosphere of interest in the study of magnetic field reconnection and substorm dynamics. A solution to this predicament could consist of creating a plasma contactor whereby a gas puff would be ionized, leading to the evacuation of positive charges and collection of cold electrons, thus compensating for the charges lost in the electron beam. A possible alternative is presented here, which consists of attaching a large passive conducting surface to the spacecraft, a “tethered capacitor”, from which negative charges would be drawn to compensate for those lost from the beam. This capacitor would then charge to a large positive potential, leaving the spacecraft and electron gun at a lower, acceptable positive potential. The tethered capacitor could have a relatively small mass; consisting only of a thin conducting surface that would be “inflated” as a result of repulsive electrostatic forces. This charge mitigation concept, as applied to active electron beam experiments, is explored using three dimensional particle-in-cell (PIC) simulations from which scaling laws can be inferred for the spacecraft and tethered capacitor potentials under proposed electron beam operations.

Predicting variations in the near-Earth space environment that can lead to spacecraft damage and failure is one example of “space weather” and a big space physics challenge. A project recently funded through the Los Alamos National Laboratory (LANL) Directed Research and Development (LDRD) program aims at developing a new capability to understand, model, and predict Space Hazards Induced near Earth by Large Dynamic Storms, the SHIELDS framework. The project goals are to understand the dynamics of the surface charging environment (SCE), the hot (keV) electrons representing the source and seed populations for the radiation belts, on both macro- and micro-scale. Important physics questions related to particle injection and acceleration associated with magnetospheric storms and substorms, as well as plasma waves, are investigated. These challenging problems are addressed using a team of world-class experts in the fields of space science and computational plasma physics, and state-of-the-art models and computational facilities. A full two-way coupling of physics-based models across multiple scales, including a global MHD (BATS-R-US) embedding a particle-in-cell (iPIC3D) and an inner magnetosphere (RAM-SCB) codes, is achieved. New data assimilation techniques employing in situ satellite data are developed; these provide an order of magnitude improvement in the accuracy in the simulation of the SCE. SHIELDS also includes a post-processing tool designed to calculate the surface charging for specific spacecraft geometry using the Curvilinear Particle-In-Cell (CPIC) code that can be used for reanalysis of satellite failures or for satellite design.

Spectral (transform) methods for solution of Vlasov-Maxwell system have shown significant promise as numerical methods capable of efficiently treating fluid-kinetic coupling in magnetized plasmas. We discuss SpectralPlasmaSolver (SPS), an implementation of three-dimensional, fully electromagnetic algorithm based on a decomposition of the plasma distribution function in Hermite modes in velocity space and Fourier modes in physical space. A fully-implicit time discretization is adopted for numerical stability and to ensure exact conservation laws for total mass, momentum and energy. The SPS code is parallelized using Message Passing Interface for distributed memory architectures. Application of the method to analysis of kinetic range of scales in plasma turbulence under conditions typical of the solar wind is demonstrated. With only 4 Hermite modes per velocity dimension, the algorithm yields damping rates of kinetic Alfvén waves with accuracy of 50% or better, which is sufficient to obtain a model of kinetic scales capable of reproducing many of the expected statistical properties of turbulent fluctuations. With increasing number of Hermite modes, progressively more accurate values for collisionless damping rates are obtained. Fully nonlinear simulations of decaying turbulence are presented and successfully compared with similar simulations performed using Particle-In-Cell method.

Active Experiments in Space: Past, Present and Future; Santa Fe, New Mexico, 11–14 September 2017

A new model describing dust charging and heating in unmagnetized plasmas in the presence of large electron emission currents is presented. By accounting for the formation of a potential well due to trapped emitted electrons when the dust is positively charged, this model extends the so-called OML+ approach, thus far limited to thermionic emission, by including electron-induced emission processes, and in particular low-energy quasi-elastic electron reflection. Revised semi-analytical formulas for the current and heat fluxes associated with emitted electrons are successfully validated against particle-in-cell simulations and predict an overall reduction of dust heating by up to a factor of 2. When applied to tungsten dust heating in divertor-like plasmas, the new model predicts that the dust lifetime increases by up to 80%, as compared with standard orbital-motion-limited estimates.

A spectral method for the numerical solution of the multi-dimensional Vlasov-Maxwell equations is presented. The plasma distribution function is expanded in Fourier (for the spatial part) and Hermite (for the velocity part) basis functions, leading to a truncated system of ordinary differential equations for the expansion coefficients (moments) that is discretized with an implicit, second order accurate Crank-Nicolson time discretization. The discrete non-linear system is solved with a preconditioned Jacobian-Free Newton-Krylov method.It is shown analytically that the Fourier-Hermite method features exact conservation laws for total mass, momentum and energy in discrete form. Standard tests involving plasma waves and the whistler instability confirm the validity of the conservation laws numerically. The whistler instability test also shows that we can step over the fastest time scale in the system without incurring in numerical instabilities. Some preconditioning strategies are presented, showing that the number of linear iterations of the Krylov solver can be drastically reduced and a significant gain in performance can be obtained.

Experiments where a high-voltage electron beam emitted by a spacecraft in the low-density magnetosphere is used to probe the magnetospheric configuration could greatly enhance our understanding of the near-Earth environment. Their challenge, however, resides in the fact that the background magnetospheric plasma cannot provide a return current that balances the electron beam current without charging the spacecraft to such high potential that in practice prevents beam emission. In order to overcome this problem, a possible solution is based on the emission of a high-density contactor plasma by the spacecraft prior to and after the beam. We perform particle-in-cell simulations to investigate the conditions under which a high-voltage electron beam can be emitted from a magnetospheric spacecraft, comparing two possible routes that rely on the high-density contactor plasma. The first is an “electron collection” route, where the contactor has lower current than the electron beam and is used with the goal of connecting to the background plasma and collecting magnetospheric electrons over a much larger area than that allowed by the spacecraft alone. The second is an “ion emission” route, where the contactor has higher current than the electron beam. Ion emission is then enabled over the large quasi-spherical area of the contactor cloud, thus overcoming the space charge limits typical of ion beam emission. Our results indicate that the ion emission route offers a pathway for performing beam experiments in the low-density magnetosphere, while the electron collection route is not viable because the contactor fails to draw a large neutralizing current from the background.

The idea of using a high-voltage electron beam with substantial current to actively probe magnetic field line connectivity in space has been discussed since the 1970s. However, its experimental realization onboard a magnetospheric spacecraft has never been accomplished because the tenuous magnetospheric plasma cannot provide the return current necessary to keep spacecraft charging under control. In this work, we perform Particle-In-Cell simulations to investigate the conditions under which a high-voltage electron beam can be emitted from a spacecraft and explore solutions that can mitigate spacecraft charging. The electron beam cannot simply be compensated for by an ion beam of equal current, because the Child-Langmuir space charge limit is violated under conditions of interest. On the other hand, releasing a high-density neutral contactor plasma prior and during beam emission is critical in aiding beam emission. We show that after an initial transient controlled by the size of the contactor cloud where the spacecraft potential rises, the spacecraft potential can settle into conditions that allow for electron beam emission. A physical explanation of this result in terms of ion emission into spherical geometry from the surface of the plasma cloud is presented, together with scaling laws of the peak spacecraft potential varying the ion mass and beam current. These results suggest that a strategy where the contactor plasma and the electron beam operate simultaneously might offer a pathway to perform beam experiments in the magnetosphere.

Local acceleration via whistler wave and particle interaction plays a significant role in particle dynamics in the radiation belt. In this work we explore gyroresonant wave-particle interaction and quasi-linear diffusion in different magnetic field configurations related to the 17 March 2013 storm. We consider the Earth's magnetic dipole field as a reference and compare the results against nondipole field configurations corresponding to quiet and stormy conditions. The latter are obtained with the ring current-atmosphere interactions model with a self-consistent magnetic field (RAM-SCB), a code that models the Earth's ring current and provides a realistic modeling of the Earth's magnetic field. By applying quasi-linear theory, the bounce- and Magnetic Local Time (MLT)-averaged electron pitch angle, mixed-term, and energy diffusion coefficients are calculated for each magnetic field configuration. For radiation belt (∼1 MeV) and ring current (∼100 keV) electrons, it is shown that at some MLTs the bounce-averaged diffusion coefficients become rather insensitive to the details of the magnetic field configuration, while at other MLTs storm conditions can expand the range of equatorial pitch angles where gyroresonant diffusion occurs and significantly enhance the diffusion rates. When MLT average is performed at drift shell L=4.25 (a good approximation to drift average), the diffusion coefficients become quite independent of the magnetic field configuration for relativistic electrons, while the opposite is true for lower energy electrons. These results suggest that, at least for the 17 March 2013 storm and for L≲4.25, the commonly adopted dipole approximation of the Earth's magnetic field can be safely used for radiation belt electrons, while a realistic modeling of the magnetic field configuration is necessary to describe adequately the diffusion rates of ring current electrons.

A spectral method for kinetic plasma simulations based on the expansion of the velocity dis- tribution function in a variable number of Hermite polynomials is presented. The method is based on a set of non-linear equations that is solved to determine the coefficients of the Hermite expansion satisfying the Vlasov and Poisson equations. In this paper, we first show that this technique combines the fluid and kinetic approaches into one framework. Second, we present an adaptive strategy to increase and decrease the number of Hermite functions dynamically during the simulation. The technique is applied to the Landau damping and two-stream instability test problems. Performance results show 21% and 47% saving of total simulation time in the Landau and two-stream instability test cases, respectively.

A number of electron beam optics elements allow for substantial control over the phase space of the relativistic bunches. We study a scheme capable of generating a beam distribution consisting of several well-separated energy bands. Such a distribution is unstable and the two-stream instability driven by the space charge develops. We investigate the instability analytically and numerically.

[1] The particle dynamics in the Earth's radiation belt is generally modeled by means of a two-dimensional diffusion equation for the particle distribution function in energy and pitch angle. The goal of this paper is to survey and compare different numerical schemes for the solution of the diffusion equation, and to outline the optimal strategy from a numerical point of view. We focus on the general (and more computationally challenging) case where the mixed terms in the diffusion tensor are retained. In Part 1, we compare fully implicit and semi-implicit schemes. For the former, we have analyzed a direct solver based on a LU decomposition routine for sparse matrices, and an iterative incomplete LU preconditioned Generalized Minimal REsidual solver. For the semi-implicit scheme, we have studied an alternating direction implicit scheme. We present a convergence study for a realistic case that shows that the time step and grid size are strongly constrained by the desired accuracy of the solution. We show that the fully implicit scheme is to be preferred in most cases as the more computationally efficient.