Annales Geophysicae Numerical simulations of thermospheric dynamics : divergence as a proxy for vertical winds

A local scale, time dependent three-dimensional model of the neutral thermosphere was used to test the applicability of two previously published empirical relations between thermospheric vertical wind and velocity divergence, i.e., those due toBurnside et al.(1981) andBrekke(1997). The model self-consistently solves for vertical winds driven by heat and momentum deposited into the neutral atmosphere by high latitude ion convection. The Brekke condition accurately mimicked the overall “shape” of the three-dimensional model vertical wind field although, as written, it consistently overestimated the vertical wind magnitude by a factor of approximately 5/3, for the heating scenarios that we considered. This same general behavior was observed regardless of whether the forcing was static or rapidly changing with time. We discuss the likely reason for the Brekke condition overestimating the magnitude of our vertical winds, and suggest an alternative condition that should better describe vertical winds that are driven by local heating. The applicability of the Burnside condition was, by contrast, quite variable. During static heating, both the magnitude and the sign of the model vertical winds were predicted reliably at heights above those of maximum energy and momentum deposition per unit mass. However, below the thermal forcing, the Burnside condition predicted vertical winds of the wrong sign. It also introduced significant artefacts into the predicted vertical wind field when the forcing changed suddenly with time. If these results are of general applicability (which seems likely, given the way these relations are derived) then the Burnside condition could usually be used safely at altitudes above hmF2. But it should be avoided below this height at all times, and even at high altitudes during periods of dynamic forcing. While the Brekke condition (or our modified version of it) could likely be used in all circumCorrespondence to: S. L. Cooper (sl7cooper@students.latrobe.edu.au) stances, there are few experimental scenarios for which this would be useful. This is because evaluation of the Brekke condition would not usually be possible unless the vertical wind was already known in advance.


Introduction
The most difficult wind component to measure in Earth's thermosphere is the vertical, for two reasons.First, vertical winds have smaller speeds than horizontal winds.Second, ground-based remote sensing of thermospheric vertical wind is usually restricted to viewing solely in the zenith, which yields limited geographic coverage.These difficulties have motivated various attempts to infer thermospheric vertical winds indirectly from other data.For example, Burnside et al. (1981) describe a simple relation between vertical wind and the divergence of the two-component horizontal wind field.However, mathematical derivation of this relation requires a number of assumptions that may not hold in the real atmosphere.Brekke (1997) presents an alternative relation, in this case between the total divergence of the three-component wind field and its vertical component.
While a number of techniques exist for actually measuring thermospheric winds, Fabry-Perot Interferometers (FPIs) are the primary remote sensing tool, especially for altitudes above the E-region.These instruments have been utilised in numerous studies at both high and mid latitudes - Rees et al. (1984b), Crickmore et al. (1991), Aruliah (1995), Price et al. (1995), Innis et al. (1996), Dyson et al. (1997), Conde and Smith (1998), Ishii et al. (1999) and Greet et al. (2002) for example.Since they only measure the component of the Published by Copernicus Publications on behalf of the European Geosciences Union.S. L. Cooper et al.: Numerical test of Burnside condition wind parallel to their look direction, a single ground-based Doppler spectrometer can unambiguously measure vertical wind only when observing directly overhead.Furthermore, to derive all three wind components, an FPI must observe in at least two additional linearly independent off-zenith directions.In a common scenario, narrow-field observations are taken in a total of five directions: zenith, plus north, east, south and west at a 20 • to 30 • elevation angle.Subject to some assumptions, all three wind components can then be derived, together with an estimate of horizontal divergence.Alternatively, a number of all-sky imaging FPIs have been developed (Rees and Greenaway, 1983;Biondi et al., 1995;Conde and Smith, 1997;Ishii et al., 2002), however most of these cannot measure vertical winds directly in the zenith because the dispersion of wavelength with zenith angle is infinite at the center of the Fabry-Perot ring pattern.The exception to this is the all-sky instrument described by Conde and Smith (1997), which measures spectra from all portions of the sky by scanning the etalon gap over time.
When observing the 630 nm emission from atomic oxygen at around 240 km altitude, a five-direction observing scheme would sample thermospheric regions that are horizontally separated by around 500 km.Despite the thermosphere's very high kinematic viscosity, all-sky imaging FPI observations (Conde and Smith, 1998, for example) have observed relatively small scale structure in the wind field.Using a Doppler Imaging Fabry-Perot spectrometer Batten and Rees (1990) reported observations of structures in the F-region wind field at spatial scales as short as ∼50 km, with temporal scales as short as ∼10 min.They concluded that a spatial resolution much better than 500 horizontal km is required to resolve the smallest scale wind features that occur.The goal of the present study is to model thermospheric dynamics with adequate resolution to study features at these spatial and temporal scales which, until now, has not usually been feasible.
The vertical component of the thermospheric wind field is the one expected to exhibit most structure over small spatial scales.Smith (1998) reviewed the recent understanding of such vertical winds.Whilst he noted that thermospheric vertical wind speeds in excess of 100 m s −1 have been observed by a number of groups, the mechanisms that produce such large speeds are not well understood.Deng et al. (2008) have completed some modelling regarding the generation of such vertical winds and present a possible mechanism.However, the current study is concerned only with more common vertical wind speeds up to about 10 ms −1 , which are an order of magnitude smaller than typical horizontal winds, and are well within the speed range for which possible generation mechanisms are known.The typical 1 − σ uncertainty of FPI wind measurements is around 5 ms −1 which is 50% of the vertical wind speeds of interest here.Such large relative errors are obviously undesirable for observations; an alternative experimental procedure for measuring vertical winds would be very appealing.Two possible alternative methods have been tested in this study, the "Burnside" and "Brekke" conditions, which are defined below.
There are two components of vertical winds that persist in the thermosphere and contribute to the total vertical wind.The first component is the "barometric" component which is driven by the expansion and contraction of the thermosphere due to solar heating, effectively the atmosphere's "breathing" motion.The second is everything else that drives vertical winds, this includes the "divergence" component that arises from diverging and converging horizontal winds.Due to the solar heating pressure gradients will drive horizontal winds and to preserve mass continuity vertical winds will arise (Rishbeth et al., 1969).Upwelling can also be driven by intense local heating, such as Joule heating or particle precipitation.The present study is concerned with Joule heating and the subsequently driven winds, in fact tidal motions are not included in the model.The atmosphere's horizontal divergence and vertical wind are coupled through the requirement to conserve mass, a condition described by the continuity equation: where t is time, ρ is the mass density and u is the 3component neutral wind velocity vector.Using this and assuming hydrostatic equilibrium, it is straightforward to derive the "tendency equation" (Hess, 1959): where p(z) is the pressure at height z, ξ is a dummy variable for integration over altitude, u h is the horizontal wind, and ∂y y is the horizontal del operator.Equation (2) does not in general produce a simple and universal relationship between vertical wind and horizontal divergence.This is partly because the time derivative of pressure and the horizontal gradient of density are also involved in the coupling.But even if we simplify Eq. ( 2) by assuming that ∂p ∂t = 0 and ∇ h • ρ = 0, knowledge of horizontal divergence at one location is not in general adequate to predict the vertical wind there.However, by introducing two further assumptions, Burnside et al. (1981) obtained the following simple approximation, referred to here as the "Burnside condition": where H is the scale height.
The two additional assumptions used by Burnside et al. (1981) in deriving Eq. (3) from Eq. ( 2) are: -The horizontal divergence of the horizontal wind is independent of height: ∂ ∂z (∇ h • u h ) = 0 -The atmosphere is isothermal at heights above z: ∂T ∂z = 0 The Burnside condition has been tested in several observational studies -with differing conclusions.These include comparisons with observational data at mid-latitudes (Biondi, 1984;Sipler et al., 1995) and at high latitudes (Crickmore, 1993;Smith and Hernandez, 1995).Biondi (1984) used an FPI to observe the two-dimensional horizontal wind field and the corresponding overhead vertical wind.It was found that vertical winds were indeed associated with horizontally divergent flows but the Burnside condition overestimated the vertical winds.By contrast, Sipler et al. (1995) determined that average vertical winds during geomagnetically quiet times were larger than the Burnside condition would suggest.At high latitudes Crickmore (1993) observed that the Burnside condition underestimated the magnitudes of the measured winds by a factor of approximately two.Smith and Hernandez (1995) compared the vertical winds observed within the polar cap to those predicted by the Burnside condition and also found that the magnitude of the actual vertical winds were double those predicted.However, they found that the sign of the observed vertical wind was opposite to that predicted by applying the Burnside condition to the observed horizontal divergence -which is a much more serious discrepancy.Nevertheless, there was at least a correlation between the two quantities.However, it was concluded by Smith (1998) that the Burnside condition appears inappropriate at high latitudes.
An alternative relation between vertical wind and velocity divergence is given by Brekke (1997, p. 260) and referred to here as the "Brekke condition": where γ is the adiabatic constant, c p /c v .Unlike the Burnside condition, the Brekke condition uses the total divergence in the neutral wind field to infer vertical wind speed.The Brekke condition can be derived from the mass conservation equation (Eq.1), subject to some approximations.First, it is assumed that the atmosphere is quasi-static allowing the time derivative of density to be removed.Further approximations used include assuming the horizontal gradients in density are zero, the ideal gas law holds, hydrostatic equilibrium holds, and that the vertical wind displaces air parcels adiabatically.All but the last of these approximations are reasonable over most spatial and temporal scales of interest.The adiabatic assumption may or may not be reasonabledepending of course on whether the vertical motion is driven by local heating.
The existence of some relationship between divergence and vertical wind has been observationally verified by the studies cited above, but the best form for this relation remains unclear.There is thus considerable practical motivation to determine which condition best describes the actual behavior of the real atmosphere.This was the objective of the present study, in which we used a high spatial resolution and fully self-consistent numerical model of the neutral thermo-sphere to compare model vertical winds with those predicted by both the Burnside and Brekke conditions.
Numerous global-scale general circulation models of the atmosphere have been developed, Fuller-Rowell and Rees (1980), Dickinson et al. (1981), Richmond et al. (1992), Harris (2001) and Ridley et al. (2006) for example.In one such study, Rees et al. (1984a) applied a global scale model to simulate the thermospheric vertical wind field on specific days, and compared the results to corresponding observational data (Rees et al., 1984b).A major problem with this approach is that the global nature of these models precludes them from resolving dynamics on the scales that are likely to be important for driving vertical winds that are the focus of the work reported here.One global model that can potentially address dynamics at high resolution is the nested grid model of Wang et al. (1999).This model allows a region of higher spatial resolution to be nested into the global-scale model.The initial Thermosphere-Ionosphere Nested Grid (TING) model of Wang et al. (1999) contained only a one way coupling between the global and nested model domains, with the coarse grid model affecting the nested model.The TING model now incorporates two-way coupling between the nested grid and the global grids (Wang et al., 2005).However, the externally specified fields usually do not have adequate resolution to match the scales being studied here.
More recently, Ridley et al. ( 2006) reported the development of the Global Ionosphere-Thermosphere Model (GITM) which can have very flexible spatial resolution, allowing for higher spatial resolution than previous global scale models.This model solves on an altitude based grid rather than a pressure based grid, which means that hydrostatic equilibrium need not be enforced.This is an important feature of both the GITM and our model, however the GITM remains a global scale model.Deng and Ridley (2006) used the GITM with a spatial resolution of 2.5 • in latitude, 5 • in longitude and 0.25 of a scale height vertically to investigate various influences on the high latitude neutral winds.Although a similar region of the atmosphere was considered the neutral vertical winds were not the focus of this study.
The authors are unaware of any other 3-D thermospheric model with the combination of spatial resolution, temporal resolution, and domain size to match or exceed that used in this study.Furthermore, most global scale atmospheric models simplify the problem by introducing approximations into the governing equations -such as enforcing hydrostatic equilibrium (not the case in the GITM).Such approximations are valid over large spatial and temporal scales.Not only do they simplify the problem, they can also remove instabilities that would occur if features in the exact solutions appear on spatial and temporal scales smaller than the resolution of the discrete numerical solver.However, structures smaller than the typical global grid's cell size almost certainly do occur in the real atmosphere, and it is important to understand what cumulative role such features may play over larger scales.Omission of these terms from the models may be quite inappropriate.Deng and Ridley (2007) ran models with varying spatial resolution in a study investigating possible reasons for the underestimation of Joule heating in global scale models.Sub-grid dynamics affected the strength and the positioning of the Joule heating in the low resolution models.
Some modelling of smaller scale structures has also been undertaken using 2-D models (Fuller-Rowell, 1985;Walterscheid et al., 1985;Walterscheid and Lyons, 1992).Fuller-Rowell (1985) modelled the neutral wind response to forcing from electrodynamic features, both narrow and broad.The resolution of the nested regions of this model were similar to the resolution used in our current work.However, being only a 2-D transient model most likely limited this study from capturing all the processes included in the present work.It is worth noting that the neutral wind and temperature responses of both models do appear reasonably consistent.The temperature changes seen in both models are a few Kelvin in amplitude, and in both cases the vertical wind response is strongest about 10 min after a change in the forcing.Similar responses were also observed in 2-D modelling work by Walterscheid et al. (1985), Walterscheid and Lyons (1992) and Chang and St.-Maurice (1991).Chang and St.-Maurice (1991) determined that advection of momentum and energy can dominate over local scales near the small scale forcing features.
The neutral wind response to a stable auroral arc was simulated by Walterscheid et al. (1985).One interesting result was a region of cooling that arose near the arc soon after its appearance.Chang and St.-Maurice (1991) used a 2-D model of the upper atmosphere and compared results between an extreme case of forcing and a moderate case.In each case wave motions in vertical winds were observed, with the response to extreme forcing being much more pronounced, as expected.Non-linear terms in the governing equations were found to dominate in the extreme case, but play a lesser role in the moderate case.
More recently Sun et al. (1995) used a 3-D time dependant model to simulate the thermospheric response to diffuse aurora.The resolution of the model was 170 km longitudinally, 75 km latitudinally and 5 km vertically.A 2-D version of the model was also run and the results showed that the 3-D model was more accurate in simulating the neutral wind response.The 2-D models were shown to generate zonal winds that were about twice the strength of the 3-D zonal winds, with the 3-D zonal winds being more representative of the realistic atmospheric response.

The model
The model used in this study employs the finite volume method (FVM) to solve a set of time dependent, nonlinear coupled three-dimensional partial differential equations.The model domain extends 1000 km in the zonal direc-tion, 600 km meridionally and 300 km vertically.(The actual height range represented is 100 km to 400 km.)The resolution of the model is 20 km zonally, 15 km meridionally and 750 m vertically.The reason for the high vertical resolution is that the model does not enforce hydrostatic equilibrium, which allows small scale waves to develop with short vertical wavelengths.Such waves are often observed, especially during the early time steps of a model run, before minor numerical inconsistencies associated with the boundaries have damped out.
Although the spatial resolution of the model is higher than previous 3-D studies, the authors cannot exclude the possible existence of structures at even smaller spatial scales in the real atmosphere.The structure that we actually observe in our simulations does not approach the model's grid size, which suggests to us that the spatial resolution currently used is adequate for the phenomena that we are considering.Higher spatial resolution may well be required to model the response to even more localised forcing fields, but we have not yet studied behaviour at such scales.
There are five main equations that define the model: the three components of the momentum equation, the energy equation and the mass continuity equation.These are the familiar Navier-Stokes equations; the versions that are used in the model are listed below.Earth's curvature is assumed to be unimportant over the limited spatial extent of the model domain, which allows a simple Cartesian co-ordinate system to be used.Effects of Earth's rotation are included, by retention of the Coriolis term in our momentum equation.

Momentum equation
The momentum equation used is: with t representing time, ρ being the neutral mass density, u is the neutral velocity vector, µ is the co-efficient of viscosity, ν ni the neutral-ion collisional frequency, v the ion convection velocity vector, g the gravity vector, is Earth's angular velocity vector and p is pressure.
The momentum equation (Eq.5) describes the time rate of change of velocity, via terms that describe advection, viscosity, ion drag, gravity, Coriolis and the pressure gradient.

Energy equation
The Energy equation is: where e is internal energy per unit mass, κ is the thermal diffusivity co-efficient and φ is a viscous dissipation term.
Ann.The energy equation (Eq.6) describes the time rate of change of energy, using terms describing advection, thermal diffusivity, Joule heating, work and viscous heating.Temperature is retrieved by dividing the energy e by the specific heat capacity at constant volume c v .

Mass continuity equation
The requirement for mass conservation is used to close the set of equations (Eq.1).Terms in this equation correspond to the time rate of change of density and divergence of mass flux.

Predefined fields
Initial temperature and pressure fields are implemented from empirical fits to MSIS model fields (Hedin, 1991).This version of the model was used as we already had access to a working version of this code.While newer versions of the model are available this model is adequate for our purposes since our interest is in the changes produced by additional heating rather than trying to exactly match the atmospheric conditions at a particular time or location.The density field is obtained from these using the ideal gas law.These fields then evolve self-consistently as the model proceeds.The ion density and velocity fields are externally specified at all time steps.For the runs shown here, ion density peaked at 2×10 11 m −3 at an altitude of 300 km, and was specified with a Chapman-like profile for convenience.The ion velocity (convection) field drives the neutral dynamics in the model, and is varied according to the problem of interest but there is only a one-way coupling of the ion convection with the neutral atmosphere.Ion velocity is purely horizontal, and is constrained to be (predominantly, due to discrete sampling) non-divergent.Thayer and Killeen (1993) showed that the ion convection is non-divergent to first order, although it is acknowledged that observations of divergent ion convection do exist (André et al., 2003).The ion convection used is a narrow channel flowing predominantly in the zonal (westward) direction and peaking at a speed of 800 ms −1 , with a kink partway along the channel.It has a similar geometry to the convection pattern used in an earlier model (e.g.Fig. 1 of Cooper and Conde, 2006), however it now extends up to 400 km in altitude and tails off below an altitude 130 km.The meridional width of the channel is approximately 150 km.
The model is based on differential equations whose solution only becomes unique when constrained by boundary conditions.However, inconsistent boundary conditions can introduce artefacts or worse, they can cause the model to fail.The boundary conditions used in the current study are as follows: the normal derivatives of all parameters are set to zero on the side boundaries.On the bottom boundary the pressure and temperature are kept constant at their MSIS values, whilst the normal derivatives of the velocity components are each set to zero.On the top boundary the normal deriva-tives of all quantities are set to zero except for pressure, the gradient of which is set to ensure hydrostatic equilibrium.Inconsistencies (such as discontinuous second spatial derivatives) initially occur at the boundaries, launching waves and other small artefacts.Because of this, the model is initialised by running it for a few hours of model time with no forcing, which allows the artefacts to dampen out.The resulting fields, with neutral velocities of only fractions of ms −1 , are then used as starting conditions for subsequent model runs, with the forcing of interest being turned on a few time steps later.
We have investigated the performance of the Burnside and Brekke conditions, with high spatial and temporal resolution, by driving the neutral dynamics with ion convection and comparing the resulting vertical winds with those predicted by the two conditions.Specifically, the ion convection pattern was initially started with a kink 300 km from the eastern end of the model where it was maintained for 5400 model seconds after which time it was moved (instantaneously) some 400 km westward, where it remained for all of the following 5400 model seconds.The instantaneous movement of the ion convection kink was implemented to allow an assessment of the impulse reaction of the model.This allows a comparison between a "static" convection pattern and a time-varying, disturbed pattern.
The time step used in the model was 1 model second, however the fields were written out every 60 model seconds (the current temporal resolution of the model output).It was found that time steps larger than 1 model second could yield unacceptably large acoustic Courant numbers, which would cause the model to become unstable and fail.A similar process was identified by Chang and St.-Maurice (1991); the time step taken by their model was half a model second.

Results and discussion
Ion convection is the only external momentum source in the model.It also deposits energy by Joule heating.For the current model runs, the ion convection was simply turned on at time zero.Since the initial fields were obtained from a steady state model run with no convection, this meant the model essentially experienced an instantaneous step change in convection, rather than a gradual build up in speed.

Zonal wind and temperature responses
Neutral wind speeds near the convection channel began increasing as soon as the ion-drag forcing was turned on.As expected, the largest wind response appeared in the zonal component, because this is the direction parallel to the major ion-drag forcing.It is interesting to note that the maximum zonal velocity, after 3 model hours, was 291 ms −1 westward which is ∼36% of the ion convection speed.An example of the evolution of the zonal wind can be seen in Fig. 1  similar to the results obtained by Walterscheid et al. (1985) who observed the maximum velocity was some ∼40% of the zonal electric field drift velocity.The winds are negative here because the ion convection is moving in the westward direction.It must be noted that the zonal winds may increase further but the model run was concluded after 3 model hours, however modelling by Deng and Ridley (2006) indicated that it took 3 h for the neutral winds at an altitude of 300 km to approach their final values or 4 h at 200 km.Our simulations show that it is the viscosity term that mostly limits the steady-state neutral winds at ∼40% of the ion speed.
After 5400 model seconds the temperature change at an altitude of 300 km had a magnitude of ∼19 K, as shown in Fig. 2. The "convection channel" can be seen and there is some heating in the convection channel prior to the kink and interestingly there is some cooling in the channel after the kink.Large temperature changes occurred near side boundaries of the model domain, but these are presumed to be a consequence of imperfect boundary conditions, rather than real physical processes in the atmosphere.Of specific interest to this study is the change in the temperature in or near the ion convection.Heating occurred in the channel upstream of the convection kink.This is expected, as Joule heating should be strongest at the location of greatest ion-neutral velocity difference.However, just downstream of the kink the temperature was actually lower than it was at the start of the model.This region of cooling coincided with a region of large (total) divergence in neutral wind, which is the likely cause of the cooling.If the work term in the energy equation (Eq.6) was large enough this would cause the region to cool, as appears to be the case in this region.Regions of cooling have also been observed by Walterscheid et al. (1985).In their case these regions were thought to be due to either ver- tical displacements of air parcels or from dynamics resulting from acoustic disturbances.

Vertical wind response
We focus here on three altitudes at which we examine the applicability of the Burnside and Brekke conditions.These are 225 km, 275 km and 325 km.These three altitudes span the height range over which momentum forcing maximized.Upward wind commenced in our model as soon as forcing was initiated.Vertical speeds maximized 6-7 min after the start of the forcing, before settling down.This rapid response has been noted previously by others.Walterscheid et al. (1985) observed that the vertical winds maximised after 8 min of forcing, while plots from Fuller-Rowell (1985) show that the vertical wind was largest after 10 min (when compared to plots after 60 and 180 min).Both of these studies used two-dimensional models.This characteristic temporal response is understood to be an indication of the time scale over which the atmosphere reasserts hydrostatic equilibrium following the initial shock from the forcing.There was another spike in vertical wind following the move of the convection kink at 5400 model seconds.Vertical winds once again maximised 6-7 min after this perturbation, then gradually settled down.

The Burnside and Brekke conditions
Line plots of the model vertical winds (blue line) as well as the vertical winds predicted by the Brekke (green line) and Burnside (red line) conditions are presented in Fig. 3. Vertical winds shown in panel (a) are not large compared to typical observed vertical speeds at these heights.Other model runs (not presented here) with larger ion densities and convection speeds have generated substantially stronger vertical winds, but the conditions simulated in the present model run were intentionally modest.The key point in panel (a) is that the Burnside condition predicts downward winds at 225 km altitude, when the model vertical wind was flowing upward.This discrepancy is not too surprising; the Burnside condition only allows upward vertical wind to be associated with divergent horizontal flow, whereas a convective upwelling cell would include horizontal divergence of both signs -convergent at the base and divergent at the top.As described earlier, instances of anti-correlation between observed vertical winds and the corresponding Burnside prediction have been verified by observation (Smith and Hernandez, 1995).
Panels (b) and (c) (of Fig. 3) show that the Burnside vertical winds above 225 km are, for the most part, quite accurate.They mimic the variations in the vertical wind field with a slightly reduced magnitude during the period of forcing.This is in contrast to panel (a) in which the Burnside and the model winds were in opposite directions.Panel (c) shows that the Burnside winds have a slightly larger magnitude than the vertical winds at this height.
Considering the Brekke winds in these plots it can be seen that in each case the Brekke condition overestimates the wind magnitude relative to that actually occurring.For example, at 225 km the Brekke winds are 75% larger, 62% larger at 275 km and 69% larger at 325 km.It is important to note that it is only the magnitude of the Brekke winds that are inconsistent with the model vertical winds; the Brekke condition predicts the sign of the vertical winds correctly in all cases.Unfortunately, the Brekke condition is of little use in observational studies, if one had adequate data to determine all three divergent gradients the vertical wind itself would almost certainly be known already.
However, the model results do illustrate an important principle: the vertical winds are more directly related to the total divergence in the neutral wind field than just the horizontal divergence.It is also evident that the Burnside vertical winds are a smoother function of time than the model vertical wind field; many of the small scale oscillations in the model vertical winds do not reveal themselves in the Burnside prediction.The Brekke condition does seem to show these small oscillations more clearly.This further illustrates the better reliability of the Brekke condition, with the overestimation of the wind speed being the only issue.A possible reason for this error is discussed later.

Static and dynamic responses
The results at three points within the modelling domain show how the vertical wind and the Burnside and Brekke conditions vary with time.Figures 4 and 5 show how the two conditions behave spatially at two separate times.These two figures have the following format: The vertical winds are plotted in the top panel, the middle and bottom panels are  plots of the vertical winds minus the Brekke and Burnside predictions respectively.Figure 4 is a slice plot at an altitude of 325 km after 5400 model seconds.At this stage the ion convection has been in the same place for the entirety of the model run (apart from the initialisation period).This shows how the two conditions behave as a result of a "static" forcing regime.The top panel clearly shows the position of the convection channel and the vertical wind response to this.The middle and lower panels show the differences between the vertical winds and the predictions.It should be noted that although the differences between the magnitude of the vertical winds and the predictions are small compared to typical vertical winds speeds observed in the real atmosphere, they are at times a large percentage (up to ∼40%) of the prevailing model vertical winds.The middle panel shows that although the Brekke condition winds are too large in magnitude, it does predict the correct location for these winds relative to the flow channel.The bottom panel shows the differences between the vertical winds and the Burnside condition.This plot shows that although the Burnside condition is on average closer in magnitude to the model vertical wind, the differences are not as systematic as the Brekke condition.The important errors in the Brekke condition at 325 km altitude are due to an overestimation of the winds whereas the serious errors in the Burnside condition are not the magnitude, but in the position and sign of the winds.Of course, these are more serious concerns.
Figure 5 presents the same quantities 360 model seconds after the ion convection has changed position, i.e., at 5760 model seconds into the simulation.These panels show the more "dynamic" forcing scenario.As before, the difference from the Brekke condition (middle panel) is larger than the Burnside condition (bottom panel) however the differences are more systematic.The Brekke winds are in the appropriate positions, but are overestimated in magnitude.The Burnside winds are quite accurate within the convection pattern, but not in some regions bordering the ion convection.This is a serious concern as it means the Burnside condition will predict vertical winds where none exist.
It must also be remembered that the Burnside condition does not perform as well at the lower heights.This can be seen in Fig. 3, showing instances where the Burnside condition predicts winds in the opposite direction.The results also indicate that both conditions are more accurate when the forcing has been in place for an extended period of time.

The Brekke condition overestimation problem
The plots show that the Brekke condition better describes the shape of the vertical wind field than the Burnside condition.The remaining concern with the Brekke condition is the vertical wind magnitude.It is interesting to note that at all three heights presented in Fig. 3 the Brekke condition overestimates the vertical winds by roughly ∼67% -which is suspiciously similar to the factor of γ that appears in the condition.Indeed, we propose here a variation on the Brekke condition (that we will refer to as the "La Trobe condition") in which the adiabatic gas constant γ is simply dropped.Since the Brekke condition otherwise behaves so well spatially, it is worth considering whether there is theoretical justification for this modification.
The approximations made in deriving the Brekke condition were listed earlier.The key problem for the winds studied here is the assumption that vertical motion occurs adiabatically -which is obviously inappropriate for our model in which vertical winds are driven by Joule heating.It has been seen from other model runs that upwelling dampens out quickly when the heating source is removed.This indicates that the upwelling can only continue if heat is constantly being added to the air parcel which, obviously, is not an adiabatic process.The parcels then rise at a rate such that they remain slightly positively buoyant with temperature just slightly higher than the background atmosphere at each altitude.If there is little or zero vertical temperature gradient in the background atmosphere then, the rising gas parcels would represent an isothermal process rather than an adiabatic one.Rederiving the Brekke condition for an isothermal process leads to the La Trobe condition that: The only difference between the Brekke condition and the La Trobe condition in their final forms is the absence of the adiabatic gas constant in the La Trobe condition.Line plots of the La Trobe condition and vertical winds are presented in Fig. 6 which has a similar format to that of Fig. 3.This figure allows a comparison of the vertical winds (blue line) and the La Trobe condition (green line) for the entire time history of the model although we are going to focus on the period after 5400 model seconds.It shows that the La Trobe condition does still over estimate the vertical wind a little.Of course it retains the Brekke condition's ability to predict the spatial structure of the model vertical wind field, because the two relations differ only by a constant.The La Trobe condition is 17% larger than the model wind at an altitude of 225 km, 8% larger at 275 km and 13% at 325 km.Perfect agreement is of course not expected, as the actual atmosphere is not isothermal.But we do obtain a vast improvement on the overestimations of the Brekke condition mentioned earlier.
Figure 7 presents plots of the differences between the vertical winds and the La Trobe condition at an altitude of 325 km after 5400 model seconds (top panel) and at 5760 model seconds (bottom panel).These plots are of the same format as the difference plots presented in Figs. 4 and 5, respectively.The color scales used in Fig. 7 are different for the top and bottom panels.The color scale used in the top panel matches that of the middle and bottom panels of Fig. 4 whilst the scale of the bottom panel of Fig. 7 matches the lower two panes of Fig. 5 allowing an easier comparison of the La Trobe condition with the Burnside and Brekke conditions.A comparison of the top panel of Fig. 7 with the middle and bottom panels of Fig. 4 shows that the La Trobe condition is more accurate than the Brekke and Burnside conditions in the "static" forcing regime.A similar comparison of the bottom panel of Fig. 7 and the middle and bottom panels of Fig. 5 illustrates how the La Trobe condition compares against the Brekke and Burnside conditions in the "dynamic" forcing regime.Once again the La Trobe condition is more accurate than the other two conditions.The La Trobe condition retains the (correct) spatial form of the Brekke condition, while largely avoiding the overestimation of vertical speeds.ing regime.Here it has been used to study the effectiveness at high latitudes of using two simple formulae for predicting vertical winds: the Burnside condition, which uses the horizontal divergence, and the Brekke condition, which uses the total divergence.Access to all the fluid fields within the model has allowed comparison between the validity of using both horizontal velocity divergence and total velocity divergence of the neutral wind field as proxies for vertical winds.The reaction of the model to the imposed forcing is similar to the response seen in other models.Whilst observational data of the same spatial resolution is not available some of the larger scale responses are consistent with observations.Furthermore, the time constant observed in the model compares favourably with observational data.The simulated response of the model seems plausible with the main results of the study being: 1.There is a relationship between divergence in the neutral wind field and vertical winds.However, it was found that the total divergence was generally a better proxy than just the horizontal divergence.
2. For this reason, the Brekke condition was able to predict the correct position and direction of vertical winds, whereas the Burnside condition only predicted vertical winds of the correct sign above the height of maximum heating per unit mass (about 280 km in this case).Below this height, the Burnside condition predicted vertical winds with the opposite direction to those that actually prevailed.
3. A new prediction, the La Trobe condition, retains the spatial accuracy of the Brekke condition whilst going a long way to correct the overestimation error.
4. Both the Brekke and Burnside conditions are more accurate under the "static" forcing regime.
5. Most importantly, this study shows that the Burnside condition is not applicable for predicting vertical winds at high latitudes below the F-Region peak, under the forcing regimes considered here.However, providing the forcing is not changing rapidly with time, it is likely to be useable above the F-region peak.Since the Burnside condition can be evaluated from commonly available experimental data, this latter result is likely the most useful practical outcome of this work.

Fig. 1 .
Fig. 1.The time evolution of the zonal component of the neutral wind field at an altitude of 300 km.The position of the data used in this figure is shown as a circle in Fig. 2.

Fig. 2 .
Fig. 2. A plot of the difference (from time equals zero) in temperature at an altitude of 300 km after 5400 model seconds.The zonal direction spans 1000 km and the meridional direction spans 600 km as shown.The black circle represents the position of the data used in Fig. 1 whilst the convection channel is outlined by the black lines.

Fig. 3 .
Fig. 3. Plots of the vertical wind speed against time.The blue, green and red lines represent the model vertical wind, Brekke and Burnside conditions predictions respectively.Panel (A) is from a point at an altitude of 225 km that was within the ion convection for the first 5400 model seconds and then outside of it for the last 5400 model seconds.Panels (B) and (C) are at the same zonal and meridional positions but at altitudes of 275 km and 325 km respectively.It is important to note that the vertical scales are different for each plot.The latitudinal and longitudinal position of the point used in this figure is marked by circles in Figs. 4 and 5.

Fig. 4 .
Fig. 4. The top panel shows the vertical wind from the model at a height of 325 km after 5400 model seconds.The model extends 1000 km zonally and 600 km meridionally, the length and width of the plot respectively.The differences between the Brekke and Burnside conditions predictions and the model vertical wind (vertical wind less prediction) are shown in the middle and lower panels.The latitudinal and longitudinal position of the data used in Figs. 3 and 6 is represented by the circles.

Fig. 5 .
Fig. 5. Same format as Fig. 4 but after 5760 model seconds.The ion convection has moved to its new position and has been there for 360 model seconds.The latitudinal and longitudinal position of the data used in Figs. 3 and 6 is represented by the circles.

Fig. 6 .
Fig. 6.The model vertical wind is plotted (blue lines) as is the La Trobe condition predictions (green lines).As with Fig. 3 the line plots present the conditions after 5400 model seconds at altitudes of 225 km (A), 275 km (B) and 325 km (C), respectively.Again, the vertical scales are different for each plot.Circles in Figs. 4 and 5 show where data were taken from.

Fig. 7 .
Fig. 7.The top panel shows the difference between the vertical wind and the La Trobe condition prediction at an altitude of 325 km after 5400 model seconds.The same scale as the difference plots in Fig. 4 has been used.The bottom panel has the same format as the difference slice plots from Fig. 5 presenting the differences, after 5760 model seconds, between the model winds and the predictions using the La Trobe condition.Note that a different color scale is applied in the top and bottom panels, see text for explanation.