3-D Mesoscale MHD Simulations of a Cusp-Like Magnetic Configuration: Method and First Results

We present a local mesoscale model of the magnetospheric cusp region with high resolution (up to 300 km). We discuss the construction and implementation of the initial configuration and give a detailed description of the numerical simulation. An overview of simulation results for the case of strongly northward interplanetary magnetic field (IMF) is then presented and compared with data from Cluster 2 spacecraft from 14 February 2003. Results show a cusp diamagnetic cavity (CDC) with depth normal to the magnetospheric boundary on the order of 1–2 RE and a much larger extent of∼5–9RE tangential to the boundary, bounded by a gradual inner boundary with the magnetospheric lobe and a more distinct exterior boundary with the magnetosheath. These results are qualitatively consistent with observational data.


Introduction
The magnetospheric cusps play a key role in plasma transport from the magnetosheath into the Earth's magnetosphere (Heikkila and Winningham, 1971;Frank, 1971). The cusps are the regions extending from the magnetic poles out towards the magnetopause. These regions are very complex by nature, not simply in terms of their physical geometry, but also in terms of plasma dynamics and spatio-temporal variations. This complex nature has spawned a rich body of space science literature. Though these features were present in magnetospheric models as early as the papers of Chapman and Ferraro (Chapman and Ferraro, 1931a,b), there re-Correspondence to: E. Adamson (adamson@mps.mpg.de) main many outstanding questions with regard to these regions. Extensive research has focused on the low and mid altitude cusp (Newell et al., 1989;Lockwood and Smith, 1992;Yamauchi et al., 1996). More recently, much research has focused on the exterior cusp. As research has progressed in this area, so has the general picture of the region. The exterior cusps have been identified as having three boundaries (Lavraud et al., 2004a;Paschmann et al., 1976;Haerendel et al., 1978;Onsager et al., 2001): (1) the lobe/mantle, (2) the dayside plasma sheet, and (3) the magnetosheath. Due to the convergence (divergence) of the magnetic field into (out of) the cusp region in the Northern (Southern) Hemisphere, the cusps act as focal regions for plasma flow and externally generated perturbations. The mapping of the open/closed boundary to the cusp region has enabled estimation of reconnection sites based on in situ cusp particle data (Fuselier et al., 2000;Trattner et al., 2004). The magnetic field in the vicinity of the cusps rotates through a complete 360 • about the dipolar axis. This ensures the existence of null-points where the dipolar magnetic field is bounded by an anti-parallel magnetosheath field. The complexity of the magnetic field orientation in the region suggests the presence of extended regions of significant magnetic shear. A typical feature of the exterior cusp is a region of extremely low magnetic field strength with respect to the surrounding field. These features have been descriptively referred to as cusp diamagnetic cavities (CDC's) due to the diamagnetic depression of the field (Erlandson et al., 1988). These CDC's are also reported to exhibit enhanced plasma thermal pressure and density, as well as turbulent magnetic fluctuations (Savin et al., 1998;Le et al., 2001;Savin et al., 2002Savin et al., , 2004Nykyri et al., 2004Nykyri et al., , 2006. Cusp energetic particles (CEP's) of up to at least 2 MeV were first reported in the high altitude polar cusp region by Chen et al. in 1996(Chen et al., 1997. The energization mechanism however, remains a topic of significant debate Chang et al., 1998). Modeling of the geomagnetic cusps is notoriously challenging. Global magnetospheric models have proven indispensable in the study of the interaction of the solar wind plasma with the Earth's magnetosphere (Fedder and Lyon, 1995;Raeder et al., 1995;Gombosi et al., 1998;White et al., 1998;Ogino, 1986;Tanaka, 1995). However, the exterior cusp region poses a significant challenge for these models. While global models do reproduce the high altitude magnetospheric cusps (identified by the typical enhanced density and diamagnetic field regions), the resolution is too coarse to accurately model the small scale features. The exterior cusp is composed of the diamagnetic region extending 1-2 R E normal to the magnetosheath boundary but much more expansive (∼6 R E ) in the transverse direction (Niehof, 2010), populated by largely stagnant plasma and bounded by a number of small scale features which vary depending on magnetosheath plasma conditions (Zhang et al., 2007). For northward (southward) IMF, field aligned flows have been reported along the CDC-lobe (Lavraud et al., 2004b) (CDCdayside plasma sheet) (Lavraud et al., 2005) boundary, consistent with outflow from a high latitude (sub-solar) reconnection site. The magnetopause in the vicinity of the cusp region has been reported to range in thickness from ∼50-6000 km, with an average thickness of ∼1600 km and a median thickness of only 800 km (Panov et al., 2008). Lavraud et al. (Lavraud et al., 2005) have associated sub-Alfvénic magnetosheath flows above the cusp with the existence of a plasma depletion layer adjacent to the cusp-magnetosheath interface. Global magnetospheric models typically consider a domain of some several 100 R E in the sun-tail direction and some 50 R E in the transverse direction. Even with the nonuniform grids commonly used in order to reduce the required grid densities in regions where lower resolution is possible, the simulations still require something on the order of 10 8 grid points, realizing a minimum resolution of ∼.1-.2 R E . At these resolutions, the cusp then poses a significant challenge for global models. It has also been noted by Siscoe et al. (Siscoe et al., 2005) that magnetic field depressions in the exterior cusp tend to be significantly greater than predicted by MHD models. This implies that the limited resolution of global models can not sufficiently reproduce the magnetic field gradients which define the diamagnetic currents within the cusp. Considering these small-scale features which add to the complexity of the cusp, it is clear that global magnetospheric models are presented with a significant challenge in modeling the region.
We have developed a mesoscale cusp-like magnetic field model in order to provide a better resolution of the entire cusp region than is possible in global models. The high resolution of the local model allows us to better reproduce the large gradients in the various plasma parameters which define the small scale features outlined above. The local model provides the additional benefit of allowing for the execution of controlled experiments. For instance, it provides a relatively easy means by which to examine the effects of varia-tions in local plasma parameters. This paper introduces this model and presents first results in order to illustrate some model properties. In Sect. 2 we discuss in detail the specifics of the model as well as the numerical methods implemented. Section 3 focuses on the method by which we generate an acceptable equilibrium configuration as an initial state for the dynamic simulation. In Sect. 4, we present initial results for strongly northward IMF and Sect. 5 consists of a summary of our results.

Numerical method
The evolution of a magnetized fluid is described by the full set of MHD equations. We implement these equations in the following form: where h = (p/2) 1/γ and ρ, p, η, u, and B represent density, pressure, resistivity, velocity, and magnetic field. The ratio of specific heats (γ ) is taken to be 5/3. All variables are normalized to typical system values, and are thus dimensionless. Length scale, density, and magnetic field strength are normalized to typical system values (L 0 = 1 R E , ρ 0 = 1 cm 3 , and B 0 = 50 nT). The pressure, velocity and time scale are then measured in units of normalized pressure p 0 = B 2 0 /8π = 1 nPa, Alfvén speed u A = B 0 / √ 4πρ 0 = 1090 km s −1 , and Alfvén time τ A = L 0 /u A = 6 s, respectively.
We utilize a leapfrog discretization method for the integration of the MHD equations. This is a finite difference method of second order in both time and space. The method is relatively easy to implement and affords the possibility of using a nonuniform grid. This reduces the required number of grid points in the current configuration at this resolution by a factor of three. The simulations presented herein are comprised of 203 × 203 × 203 grid points in x,y, and z, respectively.
The simulation code used here has been tested for conservation properties, varying grid and system size and compared with results from other methods (e.g., Schindler and Otto, 1989;Otto, 1990Otto, , 1999. We are confident from the testing and various applications of the code that the results are reliable.

Model
We construct a local model of the magnetospheric cusp region by placing a dipole at −10 R E along the z-axis with the dipole axis anti-parallel to z (no dipole tilt). The coordinates are such that the x-axis is directed toward the sun and the y-axis completes the right-handed system. The simulation domain ranges from −12 R E to 12 R E in both the xand y-directions, while the z-coordinate ranges from −5 R E to 6 R E , such that the dipole is 5 R E below the lower zboundary.
With ∇ · B = 0, the magnetic field may be represented by a vector potential A, such that B = ∇ × A. The dipole magnetic field is independent of the azimuthal angle, φ, lending a magnetic vector potential in cylindrical coordinates A = −κ(sinθ/r 2 )φ, where κ represents the dipole strength. By defining the original magnetic field in terms of a vector potential, it becomes straightforward to introduce a shielding current which effectively switches-off the magnetic field in the region z > 0, thus defining a magnetosheath region. This is accomplished through the inclusion of a hyperbolic tangent term in the vector potential. Modification of this dipolar magnetic vector potential to include the shielding current results in a vector potential of the form A = −κ(sinθ/r 2 )[1 − tanh(r cosθ/L)]φ, where L controls the width of the transition region from dipolar field to no field. This current acts to magnetically shield the dipole magnetic field along the z direction. In Cartesian components, this yields a magnetic field with components: A constant "draped" (no z-component) IMF may then be superposed in the magnetosheath region. The superposed model IMF is expressed as where φ is the IMF clock angle. Thus the model IMF lies purely in the x-y-plane. The initial magnetic field configuration for the case of strongly northward IMF may be seen in Fig. 1. It remains to complete a prescription of the plasma density and pressure. We note here, that the superposition of these magnetic fields does not result in a configuration in force balance, though this is the desired state of the initial configuration in order to insure that the plasma dynamics evolve as a result of the driving magnetosheath plasma flow and not due to the imbalance between forces in the initial construction of the fields. This initial non-equilibrium will be discussed later in this paper. This desire for force balance motivates our prescription of the pressure distribution.
Consider the MHD momentum Eq. (1). Static equilibrium requires: With these assumptions, it is clear that the requirement of force-free static equilibrium results in the magnetohydrostatic equation: As the magnetic field has already been prescribed, one is free to choose the pressure distribution. Consider the magnetohydrostatic Eq.
(2). Taking the scalar product of this equation with B: it follows that for a static equilibrium state, the pressure is constant on magnetic field lines. Thus, by expressing the initial pressure as a function of the field line equation, we ensure it's constancy on field lines. In three-dimensional spherical coordinates, field lines are given by r sinθ A φ = constant. We therefore define the initial pressure distribution as where p 1 = typical magnetospheric thermal pressure, p 2 = typical magnetospheric magnetic pressure, and α is a parameter which controls the cross-field line gradient. This gives a pressure distribution with low pressure in the dayside plasma sheet and the lobe, but an enhanced pressure filling the throat of the cusp and into the magnetosheath.
Turning now to the treatment of the density distribution, we note that the plasma densities in the magnetosheath and the cusp proper are typically a factor of 10 2 larger than those in the magnetospheric lobes (Lavraud et al., 2002;Panov et al., 2008). Thus, we choose a density distribution similar to that utilized for the pressure distribution: where ρ 0 and δρ represent the background and maximally enhanced densities, respectively. Similar to the pressure distribution, the result is a system with enhanced density along field lines having foot-points in the throat of the cusp which map into the magnetosheath. The third term populates the magnetosheath with an increased density. The initial distribution of density and pressure is illustrated in Fig. 2.
The largest gradients in the configuration discussed above are located in the center of the lower plane of the simulation domain, where the dipolar magnetic field lines converge. We reduce computational costs by using a nonuniform grid with highest grid resolution in this region of convergence. We define a grid in which the x-and y-resolutions are maximal ( x, y = 0.07) in the center of the simulation domain, and lowest ( x, y = 0.168) at the edges. The resolution in the z-direction is more uniform, with a maximum resolution of z = 0.05 at the lower boundary and minimum resolution of z = 0.069 at the top boundary. The grid is chosen to be more uniform in the z-direction in order to resolve not only the strong gradients toward the bottom of the domain where the field lines converge, but also to resolve the gra-dients within the magnetopause. The computational costs could be further reduced with the consideration of a more sophisticated grid with a high z-resolution within the magnetosphere and lower in the magnetosheath, for example.
Our initial configuration is not in equilibrium. We rely on a relaxation method to move the system towards force balance. The boundary conditions differ significantly between the relaxation and the dynamically evolving configurations. These differences will be discussed further in the respective sections of this paper. Here, it will suffice to mention the similarities in boundary conditions between the two configurations. Except where otherwise stated, Neumann boundary conditions are imposed, such that ∂f/∂x n = 0 (where f represents the general set of plasma parameters). The normal components of the magnetic field at the boundaries are calculated such that ∇ · B = 0, while a simple extrapolation method is applied to the transverse components. At the lower z-boundary, where the dipolar geomagnetic flux converges, this extrapolation does not strictly preserve the dipolar nature of the field. Rather, the approximation at the boundary results in a current sheet which generates large j ×B forces tending to evacuate plasma from the system on relatively short time scales. We resolve this issue by freezing the plasma at the lower z-boundary.

Relaxation method
In this section, we discuss the method by which we relax the initial configuration towards force balance, resulting in a sufficient equilibrium state.
In order to investigate the plasma dynamics of the region, it is desired to begin the dynamic simulation with a force free equilibrium initial configuration such that any unbalanced forces present in the initial configuration do not contribute in any significant way to the ensuing dynamics. A simple superposition of magnetic fields however, does not in general result in an equilibrium configuration. The above ansatz is an effective means of producing the desired MHD model of the region under consideration (i.e., a magnetospheric cusp) however, the system which results is not in equilibrium. Such a prescription of the pressure distribution guarantees that the thermal pressure is constant on magnetic field lines. However, the magnetohydrostatic Eq. (2) requires that the pressure gradient be balanced by the Lorenz force. The net force F ∇p + F j ×B is shown in the noon-midnight meridional plane for the initial state of the system in Fig. 3.
In order to reduce these forces and move toward an equilibrium state, a numerical relaxation is utilized, similar to that of Hesse and Birn (Hesse and Birn, 1993). In their approach, the system was frozen when the kinetic energy reached a maximum. We do not employ this energy based method. A complicated three-dimensional system such as our cusp model has the tendency to realize a kinetic energy maximum far too frequently to make such an approach efficient. During the relaxation, the system is closed (no normal flow at the boundaries). As the system evolves, the forces accelerate the fluid, converting potential energy into kinetic energy. By periodically freezing the system, thus removing the kinetic energy, the system moves toward an equilibrium state. The initial system exhibits unbalanced forces which act over various scales. There are forces which act over relatively short scales (i.e., across a current sheet) and forces which act over larger scales (i.e., imbalance between the day and night sides of the magnetosphere). The various ranges over which these forces act necessitate freezing the system at different time intervals. Initially, the system is frozen on short time scales (every 100 time steps or 0.05 s). This is necessary due to the large forces present in the initial prescription of the system. If the system is allowed to evolve over longer periods initially (i.e., utilizing a lower frequency for the relaxation), the initial force imbalance drives the system so hard that the final relaxed state no longer represents an acceptable magnetic cusp model. A high frequency relaxation efficiently damps these initial local forces, but requires far more computation time in order to sufficiently relax the forces which act over larger distances. Thus, it is preferred to strike a balance between a short relaxation interval, which necessitates a long evolution in order to realize a sufficiently relaxed state, and a longer relaxation interval which may result in a significantly different final state. The cusp model must be relaxed on short timescales in order to reduce the larger local forces. Once this is accomplished, the system may be allowed to evolve on larger timescales to reduce the long-distance force imbalance and to realize the desired quasi-equilibrium more rapidly. The actual relaxation schedule employed has been devised through extensive testing. It was determined through this testing that the system requires a short initial relaxation interval, but after the strong small-scale forces have been relaxed, the system is relatively insensitive to the chosen relaxation interval, eventually realizing the same force balance. The decisive factor then becomes more a question of computation time. We employ a relaxation schedule which removes these initial small-scale forces through a short relaxation interval and then alternate between short intervals (400 time steps) and larger intervals (4000 time steps), with preference given to the larger intervals in order to reduce computation time. Periodically reverting to short relaxation intervals ensures that small-scale forces remain damped. Defining the density of unbalanced forces in the system as F u = ∂(ρu)/∂t + ∇ · (ρuu), the momentum Eq. (1) may be rewritten as F u = −∇p + j × B. The square of this unbalanced force density may then be integrated over the simulation domain to indicate the proximity of the system as a whole to force balance at a given time: F 2 n = v F 2 u dv. It is also useful to monitor the maximum force in the system as this allows for a comparison between remnant forces in the final relaxed state and initial forces present in the dynamic system.
The results of this relaxation are shown in Fig. 4 for the case of 10 • IMF (strongly northward with a small positive B y component). Both the maximum (blue) and total integrated (green) forces are plotted. The system is being relaxed only up to a time of 60 Alfvén times. The continuation of the line plot beyond this time represents the dynamic evolution of the system and will be addressed in the next section. Note that the previously mentioned j × B forces at the boundaries are not included in this calculation. The relaxation reduces the total force in the system by roughly a factor of 500 and the maximum force by a factor of 250.
The initial configuration places a minimum in the magnetic field (which serves as a proxy for the magnetic nullpoint) at (−11.74,2.03,−0.06). This is the location within the magnetopause at which the magnetospheric and magnetosheath magnetic fields are most strongly anti-parallel. Through the relaxation, the magnetopause evolves from an initially flat current sheet in the x-y-plane, to a more hemispherical structure, bulging outward towards the center, and containing a slight dimple at the peak. The magnetic field in this relaxed state exhibits a minimum which has migrated sunward and slightly towards the noon-midnight meridian to (−8.8,1.6,0.3). The converging dipolar field defining the "throat" of the cusp also begins to tilt slightly in the direction of the anti-parallel field. Figure 5 shows the relaxed magnetic field distribution and current density in the noonmidnight meridian plane. This plane is shown here in order to illustrate the general features of the magnetic field just described, though the magnetic null-point is not contained in this plane.
The resulting distributions of density and pressure also exhibit significant evolution (Fig. 6). Both parameters have become enhanced in the throat and exterior of the cusp, while the range of values of each has roughly doubled. The xgradient of both parameters in the magnetosheath is a direct consequence of the closed system relaxing toward forcebalance. Note that the tilting of the cusp throat is also evidenced in the structure of the density and pressure.

Dynamic simulation
After the configuration has been sufficiently relaxed, a dynamic simulation may be initiated by imposing a homogeneous magnetosheath flow and modifying the boundary conditions. The boundary conditions are changed to allow inflow at the dayside boundary (maximum x-boundary) and outflow through the night-side boundary (minimum x-boundary). A position dependent (x,z) flow profile is constructed such that the streaming magnetosheath plasma flows over the top of the magnetopause current sheet. This initial flow profile is then used as the boundary condition for the flow at the dayside boundary, allowing for inflow in the magnetosheath while maintaining no flow normal to the boundary within the magnetosphere. The night-side boundary differs in that the entire boundary is open to allow for motion of the magnetopause within this boundary. Note that these are initial results, meant only to give an overview of the general features present in the simulation, as such the boundary conditions are relatively simple, though more elaborate boundary conditions are planned for future work. The magnetosheath flow is initialized with a magnitude of ∼220 km s −1 in the x-direction. The flow is sub-Alfvénic with an Alfvén Mach number of 0.7. For the following example, the resistivity is also increased by fifty percent over the value used during the relaxation in order to allow for a higher reconnection rate. The following results are from simulations with a 10 • IMF (strongly northward with a small positive B y component).
It is interesting in considering the previous plot of net forces throughout the relaxation to extend the plot so as to include the forces during the dynamic evolution of the system. This gives a means by which the remnant forces in the relaxed state may be compared with those during the dynamic evolution (see Fig. 4). The line plots in the figure for times after 60 Alfvén times represent the maximum and total forces during the dynamics. The maximum force steadily increases after the magnetosheath flow is initiated, while the total integrated force gradually decreases over the first 10 Alfvén times (between 60 and 70 Alfvén times in the figure) and subsequently remains relatively constant. This suggests that, although the maximum force is increasing, the volume containing these contributive forces is actually decreasing, such that the total force decreases. This behavior suggests that the forces in the relaxed configuration are indeed sufficiently balanced such that they result in no significant contribution to the dynamics.
Typical characteristics of the external cusps are evident in Figs. 5 and 6. The magnetic null-point is located on the night-side of the cusp, slightly offset from the noon-midnight meridian in the positive y-direction (eastward). This is a direct consequence of the IMF orientation as this is where the IMF is most strongly anti-parallel to the Magnetospheric field. This region corresponds to the largest magnetopause current due to the strong magnetic shear. The current sheet stretches out to both sides of this point around the surface of the magnetopause toward the dayside due to the anti-parallel projection of the field.

Comparison with Cluster 2 data
Plots of Cluster 2 data are shown in Fig. 7 for a cusp crossing on 14 February 2003. This event is quite interesting for a number of reasons and is discussed at length elsewhere (Nykyri et al., 2011a,b,c). Here, we discuss only the period from 18:30-19:30 for the purpose of comparison with our results from the local cusp simulation. The lagged IMF and solar wind dynamic pressure from ACE are shown in the bottom pane of Fig. 7. The IMF is northward and relatively constant for the period shown. For simplicity, we present data only from spacecraft 4, as it traverses the CDC while encountering both the lobe-CDC boundary and the CDC-magnetosheath boundary once with the exception of a brief secondary traversal of the CDC-magnetosheath boundary during this interval.
Between 18:30 and 19:30, SC4 travels from the magnetospheric lobe into the CDC and finally exits into the magnetosheath. SC4 detects a gradual inner boundary between the magnetospheric lobe and the CDC. Densities increase from magnetospheric values of roughly 0.2 particles cm −3 to values comparable to the magnetosheath (20 particles cm −3 ) by 18:53. This transition is accompanied by increasing plasma flows (up to 550 km s −1 ) with proximity to the CDC-lobe boundary. The enhanced flows are largely field-aligned with a small convective component (velocity components are not shown here, see Nykyri et al., 2011a), consistent with outflow from a high-latitude reconnection site tailward of the cusp. The CDC is indicated by strongly depressed and fluctuating magnetic field accompanied by an increase in plasma density. It is difficult, however, to identify a clear boundary between the lobe and CDC, as all parameters converge gradually towards CDC values. It appears that SC4 enters the CDC proper by at least 18:50, corresponding to the inner edge of the enhanced flow region. Plasma flow in this region is only 50 km s −1 . At 18:53 the spacecraft encounters brief but significant fluctuations in all plasma parameters nearly simultaneously with an abrupt increase of the solar wind dynamic pressure. These fluctuations are consistent with magnetosheath values and likely represent an abrupt encounter with the CDC-magnetosheath boundary caused by a compression of the magnetosphere due to enhanced solar wind dynamic pressure or possibly an encounter with a flux transfer event. Immediately following this encounter with the magnetosheath, the solar wind dynamic pressure decreases and SC4 again measures flows consistent with the inner edge of the CDC. These flows persist until 19:05 at which point SC4 reenters the CDC and corresponding region of stagnant flow. SC4 appears to remain within the CDC until 19:19 at which point it clearly exits into the magnetosheath , velocity components (x -solid black, y -dashed blue, zdash-dot red), magnetic field components, pressure (thermal pressure -black, magnetic pressure -blue, total pressure -red), current density components. The right pane shows (from top to bottom) the electric field components, plasma β, and the parallel current density. Shading represents the region of large plasma β indicating the CDC.
x = −4 R E after roughly two minutes of physical time. From top to bottom in the left pane of Fig. 8, the plots represent density, velocity components, magnetic field components, Pressure, and current density, while the right pane shows electric field components, plasma β, and parallel current density (all in simulation units).
Towards the lower boundary of the simulation domain (between −5 and −2 R E ), field-aligned flows up to ∼100 km s −1 are encountered, corresponding to the region of low density and high magnetic field strength. Adjacent to this region of enhanced flow is a region of relatively stagnant plasma (20-30 km s −1 ) convecting sunward (Fig. 8). Within this region of slow convection the magnetic field is reduced to a minimum of ∼9 nT. The transition from magnetospheric lobe plasma to CDC plasma is very gradual making it difficult to identify a specific boundary (similar to the Cluster data presented previously).
The shaded section in the plots represents the CDC (indicated by a normalized plasma β > 6). The CDC spans a distance of at most 2 R E in the z-direction (nearly normal to the magnetopause), but extends much farther in the transverse direction. This oblate geometry of the CDC may be seen in Fig. 9. These figures are cuts in the noon-midnight meridian plane, roughly where the CDC realizes its maximal extent tailward. For northward IMF, the cavity is identified as the region of high plasma β extending from the reconnection line (indicated by the maximum current density in Fig. 9), sunward. Plasma β (β = p/(B 2 /2µ o )) is a useful indicator of the CDC due to the strongly depressed magnetic field and enhanced plasma pressure within.
The external boundary of the CDC with the magnetosheath is characterized by a much more abrupt transition than the inner boundary with the magnetosheath. The sunward directed outflow from the reconnection site bifurcates into two layers, one flowing parallel to the magnetic field earthward into the cusp, the other convecting sunward along the inner edge of the cusp-magnetosheath boundary (Fig. 9). These flows are qualitatively consistent with those reported by Lavraud et al. (Lavraud et al., 2005) for instance. Enhanced sunward convection (∼60 km s −1 ) is evident along the inner edge of the CDC-magnetosheath boundary. Adjacent to this enhanced sunward convection, plasma parameters change relatively abruptly to magnetosheath values.

Discussion and conclusions
We have presented a new mesoscale local simulation of a magnetospheric cusp-like configuration. This model has been developed in order to address the plasma dynamics in the cusp region which develop on scales which are too small to be sufficiently resolved in global magnetospheric models. A detailed discussion of the construction of the model has been given. The initial model requires the use of a relaxation method in order to realize a sufficiently force-balanced equi-librium state. The details of this method are presented herein, along with the specifics of the numerical method, grid and boundary conditions.
We have presented our model results in comparison with an overview of Cluster 2 data for an event on 14 February 2003. The general features modeled by the simulation are consistent with the overall view of the cusp represented by the Cluster 2 data for northward IMF. The simulation results are for the case of 10 • IMF (strongly northward with a small positive B y component). This places the null point in the tailward side of the cusp, slightly offset from the noon-midnight meridian in the positive y-direction.
The inner boundary between the magnetospheric lobe and the CDC is a gradual transition identified by increasing plasma density, decreasing magnetic field strength, enhanced (largely field-aligned) flow consistent with lobe reconnection, and increased temperature. The CDC is the region of depressed magnetic field strength, relatively stagnant plasma, and densities comparable to magnetosheath values. The cavity has a maximum depth of ∼2 R E along CDCmagnetosheath boundary normal direction, decreasing with distance from the dipolar axis, consistent with the Cluster 2 data. The CDC has a much larger transverse extent of ∼5-9 R E . These dimensions and the shape of a relatively flat, but extended boundary agree well with observations (Nykyri et al., 2011a;Niehof, 2010). The transition from CDC to magnetosheath is much more abrupt than the inner boundary, and is characterized by an increase in magnetic field strength, decrease in temperature, and increase in plasma flow (largely −v x ). Observations show a much stronger gradient at the cavity boundary to the magnetosphere. Currently, it is not clear whether the soft transition in the simulation is caused by the limited resolution (note that the simulation has already a much higher resolution than, for instance, global models) or by physics not contained within our model. The simulation also does not show the large amplitude fluctuations of observed cavities. The nature and cause of these fluctuations is also not clear in the observations and there are indications that several processes contribute such as rapid motion of cavity boundaries, transient events, such as magnetic flux bundles moving along the cavity, or waves within the cavity or at the cavity boundary. Several of these mechanisms can be explored in our model and it is conceivable that such inputs may indeed explain the observed fluctuations.