The effects of mesoscale regions of precipitation on the ionospheric dynamics , electrodynamics and electron density in the presence of strong ambient electric fields

We have developed a new high resolution twodimensional model of the high latitude ionosphere in which nonlinear advection terms are closely coupled with the electrodynamics. The model provides a self-consistent description of the ionospheric feedback on the electrodynamical perturbations produced by auroral arc-related particle precipitation in regions with strong ambient electric fields. We find in particular that a heretofore neglected ion Pedersen advection term can introduce considerable changes in the electron density profile, the current density distribution, the conductivities and the electron temperatures. We find that the convective effects can carry the ionisation more than 150 km outside the precipitation region in a few minutes, with attendant large changes in the current distribution and E-region densities that become enhanced outside the region of particle precipitation. The production of a tongue of ionisation that slowly decays outside the auroral boundaries contrasts with the sharp geometric cut-off and associated stronger current densities found in previous studies.


Introduction
A numerical investigation has been made into the dynamics of mesoscale auroral structures, using a new numerical tool that has been developed to study transport and electrodynamics in the upper atmosphere.The ultimate goal of the model is to provide a high resolution self-consistent description of Correspondence to: J. D. de Boer (john.deboer@rmc.ca)the coupling between the thermosphere and the ionosphere in the auroral region.The model, for this reason, also includes a neutral atmospheric component.However, for the study at hand, neutral atmospheric feedback effects need not be considered: for the time scale of a few minutes which is of interest for the present publication, the neutral response is simply minimal.We therefore leave the description of the neutral part of our code to a future publication.
Our model uses a two-dimensional (2-D) domain in the meridional plane, and assumes a dipole magnetic field.The model also couples the electric potential to ensuing changes in the ionospheric conductivity through temperature and density changes.The 2-D grid allows for latitudinal resolutions as small as 400 m in size.This provides a much more detailed spatial resolution than what a 3-D or global model can.For example, the TIME-GCM model (Roble and Ridley, 1994) has a latitudinal resolution of 5 • (≈570 km).The CTIM (Fuller-Rowell et al., 1996) and the CTIPE (Fuller-Rowell et al., 2002), which is based on the CTIM, have a resolution of 2 • of latitude (≈230 km).The GITM (Ridley et al., 2006) has a variable latitudinal resolution, with a minimum of 1.25 • (≈140 km) in the auroral ovals.Even though these models all offer self-consistent treatments of the coupled thermosphere and ionosphere we show here that some important smallerscale physics (in time and space) is missing from the large scale treatment.This being stated, no direct comparison between the global models and the present study is intended, given that we are really looking here at the different physics that takes place on smaller scales.east-west structure, such that gradients normal to the plane of the domain (east-west gradients) are unimportant compared to gradients in the plane (north-south and vertical).
The geometry of our mesh reflects a dipole magnetic field geometry.Our ionospheric model is thus inherently 2-D, rather than an array of 1-D models as was done in some of the earlier work.In the process of developing a more robust and realistic model that handles nonlinear terms more easily for an improved geometry, we have been able to uncover some aspects of the dynamics of auroral arcs which were not revealed by previous studies.Specifically, the convection of ions across magnetic field lines (as they carry the Pedersen currents) and the inclination of the magnetic field from vertical are essential features of the model that introduce important differences with earlier results.
The backbone of the physical processes that we discuss here have been described in Noël et al. (2000).Briefly stated, when an increase in conductivity is introduced after the onset of an elongated east-west precipitation event and an ambient north-south electric field is present, charges accumulate at the edge of the precipitation region in an attempt to maintain uniform Pedersen currents.However, the charge accumulation immediately sets up parallel electric fields and strong ensuing current densities.In that sense, parallel currents arise from horizontal gradients in the Pedersen conductivity, σ P , in the E-region: the Pedersen current density is proportional to the product of σ P and the perpendicular component of the electric field, E ⊥ , so that wherever a horizontal gradient of the conductivity is present there must also be a parallel current to conserve charge.These parallel currents are carried almost entirely by electrons, and they are closed within the magnetospheric dynamo, which is beyond our consideration here.The parallel currents enjoy a much lower resistance than the Pedersen currents, so that E is always much weaker than E ⊥ .
An important element of the physics to be discussed here is the magnetisation parameter for each ion species, κ i , defined as the ratio of the cyclotron frequency to the momentum transfer collision frequency.The Pedersen mobility of the species can be expressed as where B is the magnetic field strength.The Pedersen mobility is at a maximum value of 1/(2B) at the same altitude as that at which κ i is equal to unity.Under the conditions studied here, the µ P,i for O + , NO + and O + 2 all reach that maximum between 119 and 120 km.At that height and with the 100 mV/m field chosen for this study, the ions have a Pedersen-component drift velocity of v x =900 m/s.The effects of this drift will be examined in this article.

Description of the model
A brief description of the important characteristics of the numerical model is given in this section.

The computational domain
The model operates within a 2-D domain which is a meridional slice through an east-west oriented arc of precipitation.The upper and lower boundaries are curves of constant altitude.The latitudinal bounds of the domain can be adjusted, but are typically centred on 70 • N. The scenarios presented used a domain spanning only about 2.5 • of latitude.The exact boundaries follow magnetic field lines through latitudes set at the lower boundary.
The earth is assumed to be spherical and the magnetic field is a centred, non-tilted dipole.Since there is no tilt the model does not represent any real terrestrial meridian accurately.However this is not a problem as long as the present model is to be run for short (less than an hour or so) simulations of representative auroral activity.

Structure of the computational grid
The computational nodes are arranged on a discrete number of magnetic field lines.The spacing between the lines can be smoothly reduced in an area of interest, as it has been for some of the results presented here.The default spacing was set at 1.6 km, but it was reduced to as small as 400 m along the northern edge of the arc to obtain better resolution.
The population of computational nodes on each field line is related to the separation between the field lines.At the bottom of the domain, the vertical spacing is equal to the spacing between field lines, and it gradually increases, starting at 120 km altitude, to 1.5 times the spacing between field lines.
The model takes its name from this mesh scheme which is structured in one direction but not in the other: "Quoit" stands for Quasi-Unstructured treatment of the Oval's Ionosphere and Thermosphere.Some detail of the mesh connectivity is shown in Fig. 1.
The neighbours for each node are determined according to the Delaunay triangulation (Delaunay, 1934).Some iterations are then performed to smoothen the mesh by allowing each node to slide incrementally towards the mean position of its neighbours.Each node is constrained to slide only along its field line.The triangulation is re-computed after each iteration.
Spatial gradients of the velocity moments (the bulk transport properties) of the neutral species are determined using the neighbours of the Delaunay triangulation, while spatial gradients of the velocity moments of the charged species are determined in two ways: for the parallel gradients, each interior node has two unique neighbours on the same field line which can be used for central or forward differences.For  the perpendicular gradients, a virtual neighbour is created on the field line on each side of a node, north and south.The values of a quantity at a virtual neighbour are computed by interpolation between the nearest node on that field line and its upper or lower neighbour.This scheme is illustrated in Fig. 2.
This interpolation scheme for perpendicular gradients introduces numerical diffusion into computations of perpendicular transport, since forward (upwind) differences are calcu-Fig.3. Detail of the alternate mesh used to assess numerical diffusion of perpendicular transport.lated using interpolation.In order to quantify this numerical artefact, an alternate discretisation was created in which the nodes are arranged in a quasi-rectangular array, using discrete values of both the shell value L and a coordinate orthogonal to L. The equations for this dipole coordinate system can be found in Sect. 2 of Huba et al. (2000).Some detail of the grid is shown in Fig. 3.
The first type of mesh (Fig. 1) is desirable for several reasons (mostly related to modelling the neutrals) and was used for the results presented in Sect. 4. Nevertheless, the results were compared to results obtained with the alternate mesh (Fig. 3).The differences were small, and they did not affect the qualitative conclusions.

Ionospheric species and their moments
Six ion species (H + , N + , O + , N + 2 , NO + and O + 2 ) are treated with a subset of the 8-moment approximation: the standard 5-moment set plus the parallel component of heat flow.The parallel heat flow is calculated with the Fourier approximation.Each of the ion species has a number density and a unique velocity vector, but the ions are assumed to have a common temperature and heat flow, which is adequate for E-region studies.
Electrons are treated with a similar subset of the 8-moment approximation.However, in the case of the electrons the parallel heat flow is modelled in a time-dependent manner including the thermo-electric term, rather than via the Fourier law.Charge quasi-neutrality is imposed on the electron number density.
The partial differential equations governing the transport of the ions and electrons may be found in Blelly and Schunk (1993) and Blelly et al. (1996).These papers also give the closed expressions for parallel electron drift and the parallel polarisation field associated with quasi-neutrality.
The model also uses the Navier-Stokes (N-S) approximation for the neutral species.The number densities of eight neutral species (H, He, N, O, N 2 , NO, O 2 and Ar) are modelled.The neutral species are assumed to have a common drift velocity v n and temperature T n .The transport equations for the neutral moments were obtained from Prölss (2004).

Coupling and source terms
The moment equations for the various species, both charged and neutral, are coupled through collision frequencies.Expressions for the ion-neutral, electron-neutral and electronion momentum transfer collision frequencies were obtained from Schunk andNagy (1980, 2000).Besides the homogeneous transport terms and the terms describing transfers of momentum and energy between species through elastic collisions, the following source terms or coupling effects were included when the rates of change of the modelled moments were computed: 1. Continuity equations: -ionisation, dissociation, excitation and electron heating due to precipitating energetic electrons -ionisation due to solar EUV and night-time EUV -ion chemistry, including recombination 2. Energy equations: -Joule heating including wave heating (Farley-Buneman instability) -inelastic electron-neutral energy transfer -the thermo-electric effect The chemistry model uses the reaction rates for O + given by St.-Maurice and Laneville (1998) and the NO + dissociative recombination rate given by Noël et al. (2000).The quantitative results depend to some degree on these particular expressions, since they are very important reactions for the conditions studied.Beside those rates, the model also includes a number of other reaction rates obtained from Rees (1989), Blelly et al. (1996), Swaminathan et al. (1998), Huba et al. (2000) and Prölss (2004).But the precise values of those rate coefficients are not expected to be very important for reproducing or testing the results presented here.
The night-time photoionisation model was adopted from the one described by Huba et al. (2000) and the algorithms in their Sami2 release 0.98.This parameterisation describes re-emission from the geocorona and EUV in starlight.A daytime solar EUV ionisation model was also adopted from the same sources, but for this study it was switched off.
The expressions for the Farley-Buneman instability (threshold field and electron heating) are due to Dimant and Milikh (2003).The expressions developed by Robinson (1986) which are quoted by Noël et al. (2005) were also used for comparison, but were found to give results fairly similar to those of Dimant and Milikh.Electron cooling rates due to inelastic collisions with neutrals were obtained from Schunk and Nagy (2000).

Computing the electric potential
The cyclotron and collision frequencies of the charge carriers correspond to time scales that are very much shorter than any of the dynamics under study.It is true that the current density in aurorae can vary significantly over periods of a second or less when the effects of Alfvén waves are considered, as Zhu et al. (2001) have done.However on time scales longer than about one second, obtaining the current density is essentially a D.C. problem since the time-rate-of-change of the magnetic and electric fields becomes negligible.A polarisation field arising from space charge density is implicit in the electric field perturbation obtained as described below, but the imbalance between n e and n i is negligible when calculating the velocity moments of the electron gas.The magnetic field perturbation arising from the current density (e.g. from Hall currents) may also be ignored, to a first approximation, as a correction to the geomagnetic field itself.
Given these reasonable approximations, the current density distribution may be computed at any instant from the 2-D network of Birkeland and Pedersen conductivities.Under this assumption of D.C. currents, the electric potential field is constrained by the charge-conservation condition ∇ • J = 0, which may be expanded to where ρdeg is the rate of charge deposition from degraded primary electrons; this is a very small term, but we include it for consistency.This has a form similar to the Poisson equation, except for the different length scales between the perpendicular and parallel directions implied by the markedly different values of Pedersen and Birkeland conductivity, σ P and σ B .This form suggests that the electric potential can be solved for by successive numerical relaxation.
To obtain a numerical expression corresponding to ∇ •J = 0, we consider a quasi-rectangular cell around each node, aligned with the field lines.The four neighbours are the same as those used for calculating spatial gradients of the charged species' moments: one each above and below on the same field line, and an interpolated neighbour on each neighbouring field line.We may write a first-order expression for the current flowing into the cell through each of the four sides: where the conductances C i are the product of either σ P or σ B , as appropriate, and the width of the cell face, divided by the neighbour's distance.A is the area of the cell.Setting I net equal to zero and re-arranging for we get Obviously the net current will be zero for some value of which is an appropriately-weighted average of the neighbours' potentials, offset very slightly by the degraded precipitating charges.
Before each time step, the conductivities are computed.Then a number of iterations are performed to relax each node's potential to a value which contributes to a consistent field for the whole domain.Those values of potential are also used as the initial values for the next time step.The potentials of the upper, northern and southern edges are set to Dirichlet conditions which correspond to a perpendicular field strength of 100 mV/m.The lower boundary has a Neumann (∂ /∂z = 0, or J = 0) boundary condition.E and E ⊥ are computed from the potential field using first-order central differences.
As mentioned, the upper boundary has a Dirichlet condition on , although this is not 100% appropriate.Even though the parallel conductivity σ B is very high above the Eregion, some resistance to parallel currents remains above the upper boundary, wherever the boundary is set.However we have made runs with the upper boundary set at both 300 and 600 km, and we find no significant change to the E-region currents.

Coordinates
For the computation of mobility and conductivity, our code uses a coordinate system whose x-, y-and z-axes are defined by B, E ⊥ and B × E, respectively.Thus ambipolar drift is inherently in the sense of negative z while Hall current is in the positive z-direction, in this coordinate system.But the physical orientation of this system varies with the direction of E ⊥ , and therefore a geographically-tied coordinate system is required in which to represent E itself.
The components of the vector E are stored in a system whose x-component is parallel to B (quasi-downward), whose y-coordinate is perpendicular to B and in the plane of the domain (chosen as quasi-northward), and whose zcomponent is defined to make a right-handed system (eastward).For the results presented in this paper, E always lies in the meridional plane.Vector quantities related to charged constituents (v i , v e and J ) are also stored in this second system.
In Sect.2.5 above, and in the discussions below, however, the coordinates x, y and z are used in a sense which is typical for 2-D studies of the auroral oval, namely x to the right (southward), y into the page (eastward) and z vertical.Since we are treating a dipole field with approximately 80 • inclination in the area of interest, there is a potential for ambiguity in this definition, and indeed by x we really mean the sense Lastly, let us note that the words "vertical" and "horizontal" are used only in respect to the effective geopotential, and not the magnetic field topology.

Conditions studied
The conditions under investigation in this paper were largely the same as those studied by Noël et al. (2000Noël et al. ( , 2005)).A discrete arc of energetic electron precipitation was imposed at the top boundary, centred on a shell whose footpoint is 70.0 • latitude.The expression for the latitudinal pattern is the same as Eq. ( 35) in Noël et al. (2000) and is where x is the distance north or south of the centre of the arc.For the results presented here, the half-width was set at = 10 km, and the sharpness of the cut-off was set at a = 5 km −1 .
The spectral flux of the energetic electrons was the same as that used by Noël et al. (2000Noël et al. ( , 2005)).It is shown in Fig. 4, and had an isotropic pitch-angle distribution.The precipitation is "switched on" at t = 0 and left on.The modelled height range covers from 90 to 300 km for the results presented in this paper.The meridional boundaries follow magnetic field lines through 67.9 • and 70.5 • N at the lower boundary for the scenarios with a southwards E, and through 69.3 • and 72.0 • N for the scenario with northwards E, the domain being biased in each case towards the direction of Pedersen advection out of the arc.
The neutral number densities (except for NO; see below) and the temperature are initialised using NRLMSISE-00 (Picone et al., 2002), specifically Beta Release 2.0 in C. The neutral winds were initialised using HWM-93 (Hedin et al., 1996) which was ported into C for this purpose.The solar conditions passed to both packages were F10.7A=81, F10.7=100, and a daily AP index of 4. The local time was set at 6 a.m. on 1 January.This choice of conditions for the neutral parameters is arbitrary, and it does not situate the ionospheric study at any particular value of MLT.
The solar zenith angle used in the night-time photoionisation model was set to a constant 100 • for this study.Although the region under study would not actually be out of direct sunlight until the sun were more than about 108 • from the zenith, the fluxes from the night-time EUV model yielded initial electron densities comparable to those used in the earlier study, and they were used as a proxy for all EUV and soft precipitation fluxes.The ionosphere was initialised by applying the night-time EUV model fluxes to the thermosphere for 20 min without precipitation.A column of this ionosphere was parameterised to use as a default starting point.For each of the three scenarios presented in the results, this parameterisation was then used to initialise a further 20 min run with night-time EUV and the particular electric field applied, but still no precipitation.A column with all of the species not in MSIS (ions, electrons and NO) was saved and used to initialise the runs presented.Therefore there is a "background" night-time ionosphere outside of the arc, whose Eregion profile is in steady state until precipitation ionisation intrudes due to advection.The initial electron density profiles are shown in Fig. 16.

Cases studied
In Table 1 we present the three scenarios that were studied.The first two use a dipole magnetic field and therefore have a B field inclination of about 80 • in the auroral zone.Thus scenario A with southward electric field also has a downward component of E. Similarly, scenario B with a northward electric field has an upwardly tilted E. Scenario A could be representative of either 06:00 MLT with a dawn-to-dusk oriented convection field, or 18:00 MLT with the opposite magnetospheric convection; scenario B represents the complementary possibilities.Scenario C uses a strictly vertical B field, which permits us to separate the effects of advection and divergence.
Table 1 also shows the two distinct cases for which results are presented for scenario A. In case 1 we have removed the Pedersen drift term in order to compare the results with previous work.Case 2, as well as scenarios B and C, include the effect of Pedersen drift.

Effect of Pedersen drift at t=30 s (scenario A)
The most significant feature obtained with this model is in the E-region ionisation.We observe that it is spread out away from the arc of precipitation, in the same direction as the applied electric field, rather than being concentrated within the arc where it is produced.This means that the gradients of σ P which cause parallel currents become steady ramps rather than sharp steps, in contrast with earlier results.
To demonstrate the importance of the convection of ionisation on the current distribution, the results of two runs are shown in Figs. 5 through 7. The only difference between these two cases is that in the left-hand frame of each figure (case 1) the physical effect of ion convection across field lines has been neglected.
Figure 5 shows the electron density after 30 s for the two different runs.The results in the left-hand panel are comparable to those obtained by Noël et al. (2005), with the electron density found to be fairly uniform throughout the arc.The results in the right-hand panel are those obtained by including the effect of ion advection.The production rate of electrons (and ions) is essentially identical to the first case.However, the ions are subjected to an equatorward (and slightly downwards) drift due to their Pedersen mobility in the background southward pointing electric field.This mobility reaches a peak at about 120 km altitude.This effect breaks the northsouth symmetry of the results obtained in earlier studies.
There is also an accumulation of ionisation below the highmobility layer due to the inclination of the magnetic field lines away from vertical.As the Pedersen drift carries the ions to lower altitudes, the drift slows down due to more frequent collisions with neutrals, which have a scale height of about 7 km.The net effect is to reduce the ion density above 120 km, most markedly within one neutral scale height of 120 km, as observed in the results from our calculations.Furthermore, after maintaining a relatively strong downward drift until about one neutral scale height below 120 km, the ions begin to slow down, so that below about 113 km they "pile up" as a result of the converging motion.
After 30 s, the effect of the ion Pedersen convergence can be seen in the height of the E-region peak.The peak ionisation rates found with the model for O, N 2 and O 2 occurred at 120, 115 and 113 km altitude, respectively.(Total electron production was at a maximum at 115 km.)Without Pedersen drift the peak concentrations of NO + and O + 2 would be at 115 and 111 km, respectively (case 1).With the inclusion of the advection terms, the peaks appear at 113 and 109 km, respectively.
The effect of the ion Pedersen drift can also be seen above 120 km.In case 1, the electron number density n e has only one maximum at 114 km altitude, seen as a bright yellow band in the left-hand panel of Fig. 5.But in case 2 there is a local minimum in electron density at about 123 km along all of the field lines within the arc, and a secondary maximum between 140 and 160 km.In case 2 the minimum at 123 km ranges from n e = 2.0×10 10 m −3 at the northern edge of the arc to 1.3×10 11 m −3 at the southern edge, compared to 1.5×10 11 m −3 at 123 km throughout the interior of the arc in case 1.As scenario C will show, the minimum is attributable mainly to advection outside the arc, which is greatest near this height.However divergence of the drift also contributes to forming the minimum at 123 km.
The differences in the results between cases 1 and 2 are just as significant for current density as they are for electron (and ion) density.Figure 6 shows the current density J as quiver plots for the two cases, at t = 30 s. Again, the left-hand panel shows case 1 with Pedersen drift being neglected and the right-hand panel (case 2) has the drift included.In case 1, there are two sharp horizontal gradients in σ P which must generate concentrated parallel currents, one upwards and one downwards, with current densities related to the sharpness of the latitudinal cut-off of the precipitation.Case 2 had two broad horizontal gradients in σ P : the quantity ∂σ P /∂x is positive throughout the interior of the arc and negative over a broad area to the south of the arc.However, the gradients of σ P are much weaker than in case 1, although they exist over a much wider area.As a result, there is a downwards current throughout the precipitation arc, and a wider region of upwards current spread out over nearly 30 km on the southward side of the precipitation region.
Figure 7 shows the parallel component of the current density J at 30 s.The left and right panels correspond to the respective panels in the previous figure.Note the different colour scales used for each panel.Case 1 had two very narrow channels of current, comparable to earlier studies.The parallel current density is high at the edges of the arc, and www.ann-geophys.net/28/1345/2010/Ann.Geophys., 28, 1345-1360, 2010 low everywhere else.In case 2 we obtain considerably lower peak parallel current densities than in the past, with distributed currents throughout the interior of the arc and on its south side.The peak magnitude of parallel current density is about 10 times lower when advection through Pedersen drift is included (case 2) than without it (case 1).
Figure 8 shows the latitudinal profile of J in case 2 at an altitude of 140 km.This figure provides a better appreciation of the shape of the current profile than the colour scale in the previous figure, and shows the asymmetry of the current distribution.

Evolution from t=30 s to 2 min (case 2)
The results presented henceforth all include the effect of Pedersen advection (case 2).In the scenarios under study, a steady state was reached within about 5 min.The electron density within the arc below 150 km reached 80% of its steady-state value within about one minute.At higher altitudes the recombination time scale is longer so it took longer to reach equilibrium.However the electron density there does not significantly affect the E-region current distribution.Outside the arc, 50 km to the south, electron density takes approximately two minutes to reach 80% of its steady-state Fig. 8. Parallel current density J (black curve) after 30 s with ion advection (case 2).This is the same data as the right-hand side of Fig. 7, but it shows J across a section through 140 km altitude.A positive current is upwards.The conductances P (red) and H (blue) are integrated along field lines and shown in units of siemens on the right-hand axis.
value: there is very little change during the first minute while the ions are advancing southwards, and then during the second minute the electron (and ion) densities increase steadily towards equilibrium.
The evolution during the first several minutes can be described as a steadily growing tongue of ionisation moving equatorward of and slightly below the area of ion production.The leading edge of this tongue is determined by the peak Pedersen drift, which for NO + occurs at a height of 119 km and was about 900 m/s at that height, for the imposed 100 mV/m field.The resulting steady state solution is shown in the next section.
One curious aspect of the approach to equilibrium was unexpected and deserves mention, namely, a region of elevated electron temperatures protruding above the advancing tongue of ionisation, roughly between 140 and 180 km altitude.Figure 9 shows a time series of electron temperature from 15 s (top panel) to 120 s (bottom panel), in which one can see this transient feature.The reason for this phenomenon is that the growing region of ionisation extends underneath a region with very low electron density.Yet the advancing incursion of ionisation has a strong latitudinal gradient of σ P , so there must be upward parallel currents to close the circuit.These parallel currents are carried by electrons which, despite their relative scarcity, still offer a much lower resistance along field lines than the abundant ions at 120 km offer across field lines.Nevertheless, the electrons must move downwards at up to 40 km/s to provide the required current density, hence the high temperature which results from friction.The electron temperature reaches a peak of 5500 K about 25 s after  the start of precipitation at 160 km altitude, 10 km south of the arc.The peak flattens out rapidly after that, while it also moves downwards and farther equatorwards.As the tongue extends to its steady-state configuration, the Pedersen conductivity gradients soften and the parallel current is spread over a very wide area, with an accompanying reduction in the required electron velocity.In comparison, the narrow electron temperature enhancement on the northern edge of the arc maintains a steady and stationary shape with a maximum value of 3300 K near 140 km altitude.

Steady state with ion advection (case 2)
As mentioned earlier, the ionospheric and electrodynamic features of the system essentially reach an equilibrium within about two minutes of simulation time.The results presented in this section are for five minutes elapsed time, but are also fairly representative of any time after two minutes.
Figure 10 illustrates how the E-region ionisation created by the arc of precipitation becomes smeared out over more than 150 km to the south of the precipitation region.The Eregion peak tapers off from n e = 3.0×10 11 m −3 at the southern edge of the arc, to 1.4 × 10 11 m −3 at 50 km to the south, to 9 × 10 10 m −3 at 100 km south, and to 4.5 × 10 10 m −3 at 150 km south.Compared to the southern edge of the precipitating region, the peak has fallen off by 1/e about 80 km to the south.Within this tail of ionisation, the peak concentrations of NO + and O + 2 are found at about 108 and 107 km, respectively.These heights are even lower than the altitudes of 113 and 109 km found within the arc, since they are not biased upwards by the ion production terms.
Figure 9 shows a stationary region of elevated electron temperature just poleward of the arc.This feature reaches a steady state quickly.It occurs because the electron density drops off to a very low value outside the precipitation arc on its north side, whereas the area of downward current density extends slightly outside the arc.The reason for this is that σ P is elevated right up to the edge, so a downward current must close the circuit just outside the arc.The intensity of the J on this edge, and the accompanying rise in T e , are not as large as those found by Noël et al. (2005), since the Pedersen drift carries most of the ionisation away from the edge of the arc as soon as it is produced.We find that T e reaches a peak of 3300 K near 140 km altitude, and that the width of the peak (between half-maxima) was about 4 km.
In order to confirm that the hot and narrow electron layer poleward of the arc was adequately resolved, we increased the spatial resolution of the mesh on the northern edge of the arc.The latitudinal profile of T e at 140 km altitude is shown in Fig. 12 for four different mesh spacings: 1 km, 800, 600 and 400 m.The peak temperature of about 3300 K was consistent, and the position of the maximum at about 500 m outside the arc was consistent to within the grid spacing of each mesh.
Figure 9 also shows a layer of elevated T e around 110 km due to the Farley-Buneman instability (wave heating).The threshold field E Thr (notation of Dimant and Milikh, 2003) for the instability was about 19 mV/m at its lowest point near 105 km altitude under the initial conditions with T e = T n .
With T e elevated as a result of the instability occurring, the parameter E Thr rises in such a way that its minimum is approximately 36 mV/m at 101 km.
The perpendicular polarisation field set up by the changes in Pedersen conductivity modifies the perpendicular electric field.The field strength rises to 109 mV/m just outside the arc and drops to 94 mV/m just inside.The effect of this perturbation is visible at the lower-left of each panel in Fig. 9 because of the sensitivity of T e to E ⊥ in the areas Fig. 12. Electron temperature profile at 140 km altitude on the northern edge of the arc, at t = 5 min, for a series of mesh resolutions from 1 km to 400 m.The horizontal axis is lateral position, positive toward the south (into the arc).The peak value of T e changed by less than 50 K, and the location of the peak was consistent within the resolution of the mesh.
of occurrence of the wave heating.Incidentally, this change in field strength over a distance of about 2 km can be used to estimate the space charge density using Gauss' law.It is roughly 4 × 10 2 m −3 , or about one part in 10 8 of the electron density, which supports our earlier assertion that the space charge density can be neglected in the ionospheric transport equations.
The Hall current density for case 2 is shown in Fig. 11.One can see that its spatial distribution is very similar to that of the electron density below 120 km.We now show that the Hall current density is quite different for a northward, or poleward, pointing electric field.

Northward E field (scenario B)
Scenario B has an electric field oriented northwards, opposite in direction to the previous scenario.The precipitation ionisation pattern is the same as in the southward E scenario, but now the ions advect northwards from their point of creation.Because of the inclination of the B field, they also move upwards.Therefore the results in this scenario cannot be described as being approximately a mirror image of those in scenario A. We only present results at t = 5 min, which effectively represents steady state.
The most profound effect of the reversal of the electric field is the upward component of Pedersen drift, which at-  tempts to empty out the band of altitudes where µ P, i is high and piles the ionisation up above 130 km where this drift begins to wane. Figure 13 shows the electron density as a function of position at t = 5 min for a northward electric field.We observe a gradual slope in the altitude of the peak electron density: 50 km to the north of the arc the peak is at 125 km altitude; 100 km to the north it is at 131 km; and 150 km to the north it is at 134 km.
Within the precipitation region, there are two E-region maxima, one at 112 km and the other at 135 km, with a subtle minimum at 122 km.The minimum is comparable to those found in scenarios A (case 2) and C, and is therefore mostly attributable to transport out of the arc, rather than to divergence of the drift.That is, the minimum is due to the −v x (∂n i /∂x) term in the ion continuity equation rather than −n i (∂v x /∂x).
The effect of the divergence of Pedersen drift can be detected in the "background" ionisation outside of the arc.The background is too weak to be visible in Fig. 13, but it is shown in Fig. 16.There is a peak n e = 1.4 × 10 10 m −3 at 98 km altitude, a minimum of 3.5 × 10 9 m −3 at 119 km, and a second peak of 1.2 × 10 10 m −3 at 146 km.The background had no latitudinal variation, so that the double peak must have been generated by the divergence term in the ion continuity equation.
The Hall current density J H is shown in Fig. 14.Note that the Hall current has changed sign from scenario A (Fig. 11), and is now eastward due to the northward E. It has a peak magnitude comparable to the previous case, but is much more localised in spatial extent.

Horizontal E field (scenario C)
Figure 15 shows the steady-state profile of electron density as a function of position with a geometry of vertical B field lines.This scenario is intended to separate the effects of B field tilt and vertical Pedersen drift from those of advection.As one might expect, the plume of ionisation stays at a constant altitude in this geometry.
Just inside the right side of the arc, there is a maximum n e = 3.0 × 10 11 m −3 at 112 km, and a second maximum of 1.5 × 10 11 m −3 at 139 km, separated by a subtle minimum of 1.4 × 10 11 m −3 at 123 km.Since there is no vertical drift here, the minimum at 123 km is attributable only to depletion caused by advection out of the arc, which is strongest at Outside the arc, there is a single n e peak at 115 km, which is quite sharp on the upper side due to the sharp cut-off in wave heating.This abrupt vertical gradient in n e is an artefact of the vertical B field geometry.

Electron density in the absence of the precipitation arc
Figure 16 shows the "background" profiles of electron density used to initialise each scenario.As mentioned in Sect.3, these data were obtained by running the model with the scenario-specific electric field and the EUV model applied, but without precipitation, until steady-state was reached.These profiles are also the condition in which the ionosphere remained "upstream" of the Pedersen-advected precipitation ionisation.These profiles on their own provide some interesting results, although we have preferred to focus on investigating the effects of discrete precipitation arcs.Yet one can see in them how just the vertical component of the Pedersen drift, due to the inclination of the field lines, influences the height and strength of the E-region maxima.

Discussion
The effect of ion drift-velocity divergence on the height and thickness of the E-region and of sporadic E-layers has been anticipated and modelled by Nygrén (1990), Kirkwood and Nilsson (2000) and MacDougall and Jayachandran (2005), although the latter two studies were focused only on the context of metallic ions in sporadic E. We find in our 'background' ionisation a confirmation of the vertical-drift phenomenon with general application to all E-region ion species.Moreover, these studies were 1-D, so they looked only at vertical convergence with the assumption of planar uniformity, i.e. (∂/∂z) (∂/∂x), which is not the case when examining precipitation arcs.Nygrén et al. (2008) have noted the importance of horizontal transport in forming sporadic E-layers.For precipitation ionisation, we find that all terms in the vector drift have an important influence on the E-region electron density profile and on the coupling of the ionospheric dynamics with the electrodynamics.

Extent of meridional ion transport
The length scale of the southward plume of ionisation can be seen to be limited by the diagonal trajectory an ion takes through the region of fast Pedersen drift.Let the neutral scale height be H n , which is about 7 km at 120 km altitude.If we consider the band of altitudes from H n above 120 km to H n below it, this rapid part of the trajectory where µ P ,i is high has a length of about 2 H n tanI , where I is the inclination of the magnetic field lines, which is about 80 • at 70 • MLat.Therefore, the ions can convect a distance of order 80 km before reaching the bottom of their layer of highest mobility and beginning to decelerate due to a higher collision frequency.Such a trajectory takes about 110 s to complete with the E ⊥ that we applied, and at either end of it v x is 65% of its maximum value.Travelling from 2 H n above 120 km to 2 H n below it requires 330 s, and v x is 26% of its maximum value at the beginning and end.
The size of the southward plume might also be thought to be limited by the chemical lifetime of the ion species.The trajectory described in the previous paragraph takes about 110 s to complete, whereas at 108 km altitude NO + has a lifetime of about 80 to 120 s, getting longer the farther south one looks in the plume.The recombination rate of the major ion species is effectively proportional to the square of its concentration.Therefore, the lifetime does affect the meridional profile of the plume, but it does not limit how far south the ionisation can extend.
One implication of this advection is that, while auroral emissions are directly related to local production of ionelectron pairs, they would appear to be difficult to relate directly to electron concentration and conductivity.

Effect of the electric field polarity
The scenarios with northward and southward E are not symmetrical, but they nearly have a rotational symmetry.Note that if we assume ν in ≈ ν 0 exp(−z/H n ), then we can write κ i ≈ exp((z−z m )/H n ), where z m is the altitude at which κ i is equal to unity (the "magnetisation boundary").Using Eq. ( 1) and substituting this form for κ i , we can say From this relationship we can appreciate that the magnitude of the Pedersen drift does not change when the polarity of E changes, and that its vertical profile is more or less symmetrical in the vertical axis about the ion magnetisation boundary, to the extent that H n is comparable over several scale heights.
There is therefore a modest degree of rotational symmetry between scenarios A and B, although it is not complete since ion-neutral chemistry is height dependent and T e and H n are not constant.
The other significant asymmetry is that in the southward E scenario, the ions are driven lower towards altitudes where they are progressively demagnetised and their E × B motion ceases.But in the northward E scenario, the ions become steadily more magnetised, and the magnitude of their E × B motion approaches E ⊥ /B.Therefore, the former scenario has very little zonal advection of ions, and therefore significant Hall current density over the whole area of the plume, whereas the latter scenario has strong zonal advection of Eregion ions, but significant Hall currents only within or near the arc.This asymmetry is obvious when comparing the Hall current distributions in Figs.11 and 14.

Assumption of zonal symmetry
The distribution of Hall current density in the southward E scenario is shown in Fig. 11.The electrons, having a magnetisation κ e much greater than unity, are nearly unimpeded in their eastward E ×B drift, whereas a reduction in the ions' E × B drift due to collisions is the source of the Hall current density.The peak speed of this drift is about 1.9 km/s under the conditions studied.Under the assumption of negligible zonal gradients, the E × B drift does not affect the plasma number density in the 2-D domain, since the drift is at right angles to the domain.However a zonal gradient in either the spatial distribution or the energy spectrum of the precipitation will lead to effects which cannot be elucidated with a 2-D model.

Ionosphere-magnetosphere coupling considerations
A standard view of the magnetosphere-ionosphere (M-I) system in numerical models of the ionosphere consists of a magnetosphere that provides an ideal D.C. voltage source at the upper boundary, while the ionosphere acts as a resistive load.
An alternate view of the magnetosphere is one that provides www.ann-geophys.net/28/1345/2010/Ann.Geophys., 28, 1345-1360, 2010 charges to the ionosphere through Birkeland currents typically carried by energetic electrons, in which case the ionosphere provides a cross-field potential drop that is consistent with its integrated conductivity; that is, the magnetosphere acts as an ideal current source.The former idealisation is more common among ionospheric models, but the latter has also been used, e.g. by Richmond and Matsushita (1975).
In both cases some magnetosphere inputs (energetic particle precipitation, and either electric field or current density) are predicated without regard to the response of the ionosphere, i.e. the M-I system is uncoupled.
In the problem considered here, however, either idealised view appears to be quite inadequate.When we consider a region of hard precipitation where a strong ambient perpendicular electric field already exists, the precipitating current densities become considerably smaller than the Birkeland currents triggered by the ionosphere itself, i.e. by the E-region Pedersen conductivity gradients.
In the context of the electrodynamical calculations of the type carried out in the present paper, the ionospheric parallel electric fields and associated Birkeland currents can have two distinct origins.One source of ionospheric Birkeland currents lies at the centre of the precipitation region, where the perpendicular electric field set up by the decelerated precipitating particles is very small, but conversely, the parallel field cannot vanish.For the case studied in the present paper, this contribution is not any more important than the already small parallel current densities carried by the precipitating electrons.The second source of ionospheric currents is caused by the change in Pedersen conductivities introduced by the precipitating particles as they ionise the E-region.It is associated with the couplet of upward and downward parallel electron currents which must accompany perpendicular ion motion.This strong source of parallel currents has been shown here to be much more strongly antisymmetric than anticipated by earlier work, e.g.Noël et al. (2005), and to extend to considerable distances from the precipitating region, owing to the advection of plasma in response to E-region Pedersen drifts.Irrespective of this antisymmetry, we are looking here at the ionosphere as a source of parallel (Birkeland) currents which is not in phase with the precipitating currents and not even centred on the precipitating currents.This suggests that coupling of the M-I system may be essential to a complete understanding of the electrodynamics of precipitation arcs.
It is also possible that the sudden onset of Birkeland currents resulting from ionospheric dynamics could be a source of Alfvén waves launched from the ionosphere.That is, transients in the 2-D meridional current system might arise from temporal changes in the resistive load, as well as from changes in the voltage source.A study of this mechanism would, however, require a different temporal resolution than what we have used here.

The Cowling effect
A partial Cowling effect is associated with the enhancement of auroral electrojet currents.See Yasuhara et al. (1985) for an explanation of the application of the Cowling effect to auroral geometry and observational evidence.However this effect is operative when a dynamo of atmospheric origin is generating a westward electric field, and the ionosphere acts as a generator, rather than a load, in the 2-D (meridional) current pattern.In this situation, those investigators find a partially effective Cowling mechanism, where there is enough resistance to field aligned currents that a southward polarisation field partially blocks the northward Hall current.The electric field is of order 10 mV m −1 .
The situation we are studying is one in which the magnetosphere is driving a field aligned current pattern in the meridional plane.As Yasuhara et al. note, the Cowling effect would not be dominant in such a situation.But it has been suggested, for example by Amm et al. (2008), that intense horizontal currents associated with aurorae could be the result of a Cowling mechanism.Now, the magnetospheric dynamo in our geometry produces a zonal Hall current.This Hall current would have to be closed via field aligned currents into the magnetosphere at each end of the east-west arc.If there is resistance to this current closure, then the possibility of a partial Cowling effect with a different geometry is created.However the extra current which the magnetosphere is required to close in this case is of measure W/L compared to the currents which it is driving, where W is the north-south width of the arc and L is its east-west length, compared to a measure L/W when the zonal electrojet is driven by an atmospheric dynamo.Therefore we do not think such a Cowling-type effect will be important in cases where the magnetospheric dynamo generates a strong electric field at right angles to an arc of precipitation.
We hope to extend our study to include a zonal component of the electric field.In an extreme case, its value could be implicitly coupled so as to eliminate the zonal current integrated over the meridional plane.This would simulate the effect of an east-west polarisation field resulting from Hall currents parallel to the arc that cannot be closed within the magnetosphere.A smaller value of the zonal field is probably more realistic.Amm et al. (2008) note that a 3-D model should be required in order to study the Cowling effect at high latitudes, so there may be a limit to what we can obtain with a simple extension of the present 2-D study.However a zonal electric field component remains a priority for our future work.

Conclusions
As a result of the present investigation we conclude that Eregion ionospheric conductivities in the auroral zone cannot be derived directly from a local production-loss equilibrium Ann.Geophys., 28,[1345][1346][1347][1348][1349][1350][1351][1352][1353][1354][1355][1356][1357][1358][1359][1360]2010 www.ann-geophys.net/28/1345/2010/or even from local time histories of the same; advection of charge carriers plays an essential role on any scale less than about 150 km.Besides passing electric current, the conductivity also convects itself.The conductivity of the E-region is typically described by the three unique elements of its conductivity matrix.However, the full story is much more complex, requiring an appreciation of the trajectory the ions take after their production.Two minutes constitutes a long time scale for some auroral activity; therefore the steady state may not necessarily constitute a realistic way to study the electrodynamics of arcs.However, we must stress that, even over a time period less than 30 s, the advection of charge carriers has a strong effect on the parallel currents and electron temperature associated with an arc of precipitation.In particular, the parallel currents are more distributed than expected from a simple consideration of the precipitation boundaries, as they spread out in the direction of the E field over a characteristic distance of about 80 km.

Fig. 1 .
Fig. 1.Detail of the mesh at the lower boundary.This example shows a resolution of 600 m along a field line in an area of interest, gradually coarsening to 1.9 km in the rest of the domain.

Fig. 2 .
Fig. 2.The construction of a virtual neighbour C used to calculate a perpendicular gradient at the node x 0 .A is the nearest real node to C, and it is always a Delaunay neighbour of x 0 .B is the upper or lower parallel neighbour of A, chosen such that C is interpolated.

Fig. 4 .
Fig. 4. The energy spectrum of the precipitating electrons used for the runs shown.Reproduced from Noël et al. (2000).

Fig. 5 .
Fig. 5. Electron density after 30 s. Case 1 (left) suppresses ion advection.Case 2 (right) shows the results with the correct transport terms included.

Fig. 6 .Fig. 7 .
Fig.6.Quiver plots of current density J after 30 s without (case 1, left) and with (case 2, right) ion advection.The vector scale is the same on both sides, however on the left the field intensity on the edges of the arc is so strong that the arrows coalesce.(The arrow heads scale with the length of the arrows, and they appear foreshortened near 110 km because they have a component out of the page.)

Fig. 9 .
Fig.9.Electron temperature at 15, 30, 60 and 120 s.The spatial scales and colour scales are the same in all four panels.The precipitating region lies between the inner edges of the two enhanced electron temperature regions seen in the top panel.The pale yellow area in the first panel is off scale; the maximum value of T e in that panel is just under 5500 K.The scale was chosen to make the next three panels more readable.

Fig. 11 .
Fig. 11.Hall current density J H after 5 min.The Hall current is westward here, and is shown as a positive value.

Fig. 14 .
Fig.14.The Hall current density J H at t = 5 min with a northward E field (scenario B).Note that its direction has changed compared to Fig.11.The Hall current is eastward here (into the page), which happens to be negative in our computational sign convention.

Fig. 16 .
Fig. 16.These profiles of steady-state electron density were obtained with the background EUV ionisation only, and the electric field specific to each scenario.

Table 1 .
A key to the scenarios and cases for which results are presented.All scenarios used a 100 mV/m field and the same precipitation pattern and spectrum.