Numerical simulation of bore generation and morphology in thermal and Doppler ducts

We perform a series of numerical simulations employing a model describing nonlinear incompressible dynamics to explore bore generation and morphology for several ducting geometries. These include idealized thermal ducts in close proximity to, and more remote from, reflecting boundaries, an idealized Doppler duct remote from boundaries, and a combination of thermal and Doppler ducts at nearby altitudes. Based on observed bore environments, we assume initial long-wave perturbations having sinusoidal form and large amplitudes. Results indicate that close proximity of a duct to reflecting boundaries enhances bore development, that both thermal and Doppler ducts support bore development for long-wave perturbations of sufficiently large amplitudes, and that bore structures and evolutions exhibit both features anticipated by idealized bore theory and departures that remain to be explored in detail.


Introduction
The nonlinear phenomena known as bores and solitary waves have long been studied and are fairly well understood theoretically for idealized flows.As early as the late 1800s Lord Rayleigh considered conservation laws that captured some of the essential characteristics of the hydraulic jump (essentially the same as a bore in hydraulic systems), including the velocities and energetics involved.Thereafter, Korteweg and de-Vries (1895) presented the now-famous KdV equation, which was able to describe the dynamics of propagating solitary waves, or solitary wave trains, of unchang-ing shape.Later work by Benjamin et al. (Benjamin, 1967;Davies, 1967;Ono, 1975) extended the KdV theory and provided a simplified model of the equations of motion valid for continuously varying density profiles.An important result from this work is that internal fluids can support solitary waves and bores, just as hydraulic channel flow does.Essential to these solutions is the balance between nonlinear steepening and dispersion.Solitons represent a balance of these two effects, and when a disturbance of sufficient amplitude is present in a suitable environment, the feature will steepen and evolve into a soliton or a solitary wave train.The other essential feature of BDO theory is the existence of a horizontal duct for the disturbance to propagate along and vertically away from which the disturbance decays algebraically.For an in-depth numerical treatment of the BDO equation for a variety of initializations the reader is referred to Christie (1989).
Theory pertaining to nonlinear waves in continuouslystratified fluids was initially developed in applications to the oceanic thermocline.Thus, the initial theoretical framework was already in place when the first apparent bores were observed in the mesosphere and lower thermosphere (MLT) during the ALOHA-93 campaign (Taylor et al., 1995).The features seen in this and a series of subsequent observations (Medieros et al., 2001;Batista et al., 2002;Chung et al., 2003;Smith et al., 2003Smith et al., , 2005;;Brown et al., 2004;She et al., 2004;Fechine et al., 2005Fechine et al., , 2009;;Shiokawa et al., 2006;Snively et al., 2006) all share similar characteristics: their spanwise extents are much greater than their streamwise scales, they travel for hundreds of kilometers at the same altitude without significant dissipation, they exhibit a sharp "front" or onset, and they demonstrate crest creation.
The resemblance between the airglow observations during the ALOHA-93 campaign and an undular tidal bore described by Tricker (1965) led to the first theories attempting to explain these phenomena.These efforts, by Dewan andPicard (1998, 2001), explored extending insights from Published by Copernicus Publications on behalf of the European Geosciences Union.
B. Laughman et al.: Numerical simulation of bore generation and morphology hydraulic theory to the atmosphere, the most notable of these being the need for a duct and the nonlinear character of the event, including crest creation.
Further work in this area was undertaken by Seyler (2005).In this work, stratified fluid equations describing incompressible nonlinear dynamics were utilized, along with the ducting environment proposed by Dewan and Picard, to examine the temporal evolution of a long wavelength wave confined in the vicinity of a thermal duct.These results demonstrated steepening exhibiting similar behavior to that predicted by the BDO equation, and suggested that mesospheric bores may be excited in this manner.This work is especially noteworthy for two reasons.First, it is both quantitative and predictive in nature, and can evolve the model equations for a number of ducting environments and initializations.Second, unlike the BDO equation, which is an approximation to these dynamics, Seyler modeled the incompressible Navier-Stokes equations, thus providing more quantitative and realistic results to guide interpretation of observational data.
The goal of this paper is to extend these previous numerical studies to more realistic and more general atmospheric ducting environments.The dispersion relation (see below) that governs linear ducted wave responses in the atmosphere is a function of both the local mean temperature and velocity profiles.Indeed, both fields have been shown in theoretical and observational studies to contribute to ducting of linear gravity waves (Chimonas and Hines, 1986;Fritts and Yuan, 1989;Isler et al., 1997;Walterscheid et al., 1999;Snively and Pasko, 2005;Fritts and Janches, 2008).The excitation of these motions, however, is believed to occur through both linear and nonlinear processes, including local body forcing via momentum deposition, tunneling from adjacent regions of vertical propagation, and wave-wave interactions leading to frequency or wavenumber doubling (Chimonas et al., 1996;Hecht et al., 2001;Vadas andFritts, 2001, 2002;Walterscheid et al., 2001;Snively andPasko, 2003, 2008).Thus, both thermal and Doppler ducting will be explored.Additionally, our modeling will employ environments with a duct both in close proximity to, and more isolated from, upper and lower boundaries in order to explore their impacts on bore evolution and structure.
Our focus in this paper is on the ducting environments that can support bore-like development and responses rather than on the potential geophysical sources of these perturbations.As noted above, there appear to be a number of potential sources of linear gravity waves occurring at thermal and Doppler ducts, and we anticipate that this will also be true for nonlinear responses.Since an essential characteristic of undular bores is the creation of solitons and/or soliton trains from an initial disturbance, and these are often seen to occur in environments suggestive of the presence of much longer waves, we explore the responses to initial sinusoidal, but large-amplitude, perturbations in various idealized environments having thermal, Doppler, or both thermal and Doppler, ducts.The numerical model employed for these studies is described in Sect. 2. In Sect. 3 we present the results of this model for different ducting profiles.Section 3.1 is focused on thermal ducting, with subsections examining the role of the vertical boundaries as well as variations in the thermal structure.Section 3.2 presents the results of an isolated Doppler duct, and we present in Sect.3.3 some preliminary results of a dual-ducting environment containing both a thermal and a Doppler duct.We conclude in Sect. 4 with a few summary remarks and possible direction for future work.

Numerical model and boundary and initial conditions
Our current studies employ a numerical code solving the incompressible Boussinesq Navier-Stokes equations for direct numerical simulations (DNS) of bore generation and morphology assumed to occur in two spatial dimensions.Previous applications of this code have described threedimensional (3-D) studies of nonlinear dynamics, turbulence transitions, and turbulence evolutions and statistics for both stratified shear flow (Kelvin-Helmholtz) instability, or KHI, and gravity wave breaking (see Werne andFritts, 1999, 2001;Fritts et al., 2003Fritts et al., , 2006Fritts et al., , 2009a, b), b).The code employs a pseudo-spectral solver that computes the nonlinear advection terms in physical space.Linear terms and derivatives are handled in Fourier space, and transformations between physical and Fourier space are performed by high-radix FFTs.Potential temperature, vertical velocity, and the vertical component of vorticity are advanced in time with a 3rd-order Runge-Kutta scheme (Spalart et al., 1991).The horizontal velocity field is computed from these quantities after each time step, which is computed dynamically during the simulations to satisfy a CFL condition of 0.68.
The equations being solved are the incompressible Navier-Stokes equations subject to the Boussinesq approximation which support waves governed by the dispersion relation, Here u=(U 0 +u, v, w) is the full 3-D velocity field, P is the pressure field, ρ is the density, θ is the potential temperature, and g is gravity.Subscripts denote mean values except in the time derivative terms, N is the buoyancy frequency defined as N 2 =(g/θ 0 )(dθ 0 /dz), k=2π/λ x and m=2π/λ z , λ x and λ z are horizontal and vertical wavelengths, and ν, K, and α are kinematic viscosity, thermal diffusivity, and the coefficient of thermal expansion equal to 1/θ 0 .
A non-dimensional length scale is chosen to be the full width at half maximum (FWHM) of our stability duct, and our time scale is chosen to be the bouancy period at the maximum of the stability duct.By choosing a value of 9.55 m/s 2 for g and a value of 8500 K for θ 0 we define our temperature scale.With our velocity scale defined by our length and time scales we can select a Reynolds number appropriate for the mesopause.With a length scale of 6000 m and a stability maximum N=0.0628 s −1 we select a Reynolds number of ∼4000.This corresponds to a kinematic viscosity of ∼500 m 2 /s, which is unrealistically large for mesopause altitudes.This was done to suppress smaller scale motions in the simulation which would have dramatically increased resolution requirements.To compensate for the affects of this increased viscosity on the evolution of the means simulations were performed with the diffusion terms of the Navier-Stokes equations set to zero for the mean components of potential temperature and velocity.Finally we set ν=K resulting in a Prandtl number of 1.
All simulations presented in this paper assumed 2-D nonlinear dynamics and a domain having periodic horizontal boundary conditions and reflective upper and lower boundaries.All simulations also employ 1000 spectral modes in the horizontal and 501 modes in the vertical.The odd mode in the vertical representation allows the imposition of reflective (or rigid) upper and lower boundaries.This influence is omitted from post-processing and leads, on occasion, to a small apparent asymmetry about the midpoint in the vertical.The code also requires specification of the initial mean potential temperature and horizontal wind profiles, from which the stability profile, N 2 (z), and mean Richardson number, Ri=N 2 (z)/U 2 z (z), where N and U are the atmospheric stability and mean horizontal wind and the subscript denotes a vertical derivative, are obtained.
All our simulations are initialized with a sinusoidal waveform imposing the initial vertical velocity and potential temperature perturbations.The vertical wavelength is chosen to be comparable to the thickness of the duct in each case and the perturbation is confined in the vertical at the duct by specifying an amplitude envelope function for all runs except those described in Sect.3.1.1where the width of the domain is sufficient to confine the perturbation to the duct.The initial vertical velocity and potential temperature perturbations and the envelope function are written as All runs use a value of 0.2 for the non-dimensional background stability, N 0 .The thermal ducts in Sect.3.1.1use a maximum stability of ∼1 and zero background wind.The vertical wavelength remains unchanged for all runs with a non-dimensional value of 2. The non-dimensional horizontal wavelength is 20 for all runs except when noted in Sect.3.2, where it has a non-dimensional value of 40.Finally, the nondimensional perturbation amplitude, A, has a value of 0.14 for all runs except when noted in Sect.3.1.2and for the ducting structures present in Sect.3.3.

Results
We describe here the results of numerical experiments of bore excitation and evolution in both thermal and Doppler ducting environments.These results explore the influences of boundaries on the evolution of bore-like responses, as well as bore evolution for both idealized thermal and Doppler ducts, and in two cases where both duct types occur in close proximity.Our results exhibit both significant similarities to, and departures from, expectations based on the BDO equation examined extensively by previous authors.In particular, our results to be discussed in greater detail below indicate that close proximity of a duct to reflecting boundaries enhances bore development, that thermal and Doppler ducts support bore development for long-wave perturbations of sufficiently large amplitudes, and that bore structures and evolutions exhibit some features anticipated by idealized bore theory and significant departures from the idealized BDO solutions that remain to be explored in detail.Results are presented in a frame moving with the leading edge of the steepening in order to study the evolution more clearly.

Thermal duct -shallow domain
Our first simulation addressed the ducting environment posed by Dewan and Picard (1998) and examined numerically by Seyler (2005).The specific stability profile employed by Seyler (2005) was a hyper-Gaussian yielding a sharp N 2 (z) maximum centered in a domain with weak stability on either side.Because our model requires initialization with a mean temperature profile rather than a stability profile, it proved difficult to specify the hyper-Gaussian stability profile employed by Seyler (2005), N 2 (z)=δ 2 +exp(−z 4 ).Instead, we employed a stability profile composed of hyperbolic tangents that yielded an N 2 (z) profile that conforms closely to that of Seyler (2005).This stability profile was then integrated to provide an analytic expression for mean temperature composed of logarithms of hyperbolic cosines.These two stability profiles are shown together in Fig. 1a and exhibit close agreement.Hereafter this profile will be referred to as our composite profile.
The vertical extent of our model domain in this simulation is 2, with a horizontal extent (and initial perturbation wavelength) of 20 and a proximity to the rigid upper and lower boundaries displayed in Fig. 1a.The initial perturbation (Eq. 3) was chosen to match that of Seyler (2005), and its nonlinear evolution proceeds as expected.The envelope   , 34, 45, 88, 113, 191.function was not necessary for this simulation due to the narrow vertical domain and was not used.The initial perturbation steepens and evolves into a soliton train as anticipated by the BDO equation and modeled by Seyler (2005).Contour plots and horizontal traces of potential temperature displaced from the center of the duct are shown in Figs. 2 and 3.These fields demonstrate clearly the steepening of the initial sinusoidal long-wave perturbation and the creation of an increasing number of distinct solitons as the waveform evolves.
In order to address less idealized environments, we also considered the same evolution, but for a thermal duct described by a cosine variation in N 2 (z) given by where N 0 retains its value of 0.2.Our expectation that quasisinusoidal large-scale, low-frequency gravity waves (GWs) contribute a majority of the temperature variance throughout the MLT (Fritts and Alexander, 2003) suggests the use of such a sinusoidal ducting structure.The stability profile and its associated temperature profile are depicted by the dashed lines in panels (b) and (c) of Fig. 1 for comparison with the profile employed for the simulation described above.
The most significant difference from the composite profile is that the cosine variation of N 2 (z) leads to a somewhat narrower ducting structure.The results of this simulation are displayed with horizontal traces of potential temperature displaced from the center of the duct in Fig. 4 at the position and the same times as in Fig. 3.As expected, these results differ very little from those of Seyler (2005) and shown in Fig. 3.The major quantitative difference is the apparent more rapid evolution of solitary wave responses for the sharper and narrower cosine thermal duct.Referring to our chosen nondimensionalization, these shallow domain simulations show that for a duct with a FWHM of ∼6 km thickness and maximum buoyancy frequency of ∼0.0628 s −1 crests arise with an effective wavelength of ∼8 km and a propagation speed of ∼80 m/s in reasonable agreement with observed bore characteristics.

Thermal duct -deep domain
Much of the response seen in Sect.3.1.1(and Fig. 2) lie in close proximity to the rigid upper and lower boundaries.In an attempt to isolate the effects of the boundaries from the dynamics at the duct, we increased the vertical domain by a factor of 9 from 2 to 18.The initial perturbation used was similar to that in Sect.3.1.1except that it is confined by the envelope function (Eq.4).Aside from the change in the vertical domain size and the addition of the envelope function, all other model parameters remain unchanged.Figure 5 shows the potential temperature profiles for these deeper domain simulations.
In the context of BDO theory and the results of Sect.3.1.1,the observed evolutions were unexpected.Neither our composite duct nor the sharper cosine duct (Eq.5) led to steepening and soliton evolution in the deeper domains.Rather, the response was only weakly nonlinear, and despite the reduced  stability away from the thermal duct, wave energy penetrated throughout the vertical domain.Figure 6 depicts potential temperature contours in the center fifth of the domain initialized with the composite profile.The time evolution depicted in panels (a-f) clearly demonstrates a weakly nonlinear response, with the scales of motion remaining largely on the order of the initial perturbation.Figure 7 displays horizontal traces through the contours in Fig. 6 and exhibits a sharp contrast to the evolution and steepening seen in Fig. 2. While the behavior is clearly nonlinear, it exhibits no systematic steepening and no emerging bore-like behavior at later times.Similar results were obtained for the cosine duct.
Possible explanations for the failure to develop a bore were explored.Since steepening is a nonlinear affect, we increased the amplitude A of the perturbation by a factor of 2.5 to a non-dimensional value of 0.35 in order to increase the influence of nonlinearity.This was done with both the cosine and composite ducts with mixed results.
Figure 8 displays the temporal evolution of the potential temperature field when the model was initialized with the  cosine duct in our deep domain with a vertical extent of 18. Figure 9 consists of traces through the contours at z=−0.18.It is clear from both figures that the initial perturbation experiences nonlinear steepening; the apparent soliton-like feature located at x=14 at early times and the feature that arises at later times at x=10.5 argue strongly for this.Clearly, however, this is not a simple BDO-like response as seen in Sect.3.1.1,and the morphology of the response requires more analysis to be fully understood.Factors contributing to the departure from the BDO response likely include radiation away from the duct, which we believe to be the cause for the non-bore-like response seen in Figs. 6 and 7.As previously mentioned, the role of viscosity was briefly explored, and it appears that it does play a role in removing energy from the response.Another way in which viscosity could affect the solution is through diffusion (molecular or turbulent) of the mean ducting structure.A ducting structure changing in thickness and strength will affect the evolution of the perturbation on longer time scales.The model runs presented in this paper were rerun with no diffusion of the mean profiles, and the results were quantitatively very similar and thus are mentioned for completeness but not presented here.
An additional factor contributing to the departure from the expected BDO response likely lies with the exact structure of the duct.When the perturbation amplitude was increased for the composite duct, the response failed to achieve a borelike structure, while that for the cosine duct did yield a bore.Given that both stability profiles have the same maximum and minimum values, it seems that the exact form (or depth) of the duct plays a role in defining the response, which is consistent with the predictions of KdV and BDO theory.
A measure of nonlinearity for surface waves described by KdV theory is the parameter aλ 2 /h 3 , where a is the elevation of the perturbation above the free, unperturbed surface, λ is the wavelength, and h is the thickness of the duct (Lighthill, 1978).When this parameter is close to 1 the perturbation evolution is nearly sinusoidal and linear.As this parameter increases to about 16, the evolution becomes steeper and more strongly non-linear and cnoidal waves develop.When this parameter equals 16 solitary waves develop, and for values above 16 wave breaking occurs.The analogue of this parameter from BDO theory is a λ/h 2 (Christie, 1988;Koop and Butler, 1981).As is evident in Fig. 1 the cosine duct is narrower than the composite duct.Since the degree of nonlinearity depends inversely on the thickness of the duct, this appears to offer one explanation why nonlinear features were observed for the cosine duct, but not for the composite duct in the deep domain.Narrowing the width of the composite duct resulted in behavior that appeared slightly more nonlinear, but not to the extent exhibited by the cosine duct, indicating that the gradient of the stability profile may play a role in supporting this nonlinear behavior.Estimates of the nonlinearity from the parameter aλ/h 2 were ∼5 for the shallow domain of Sect.3.1.1and support the apparent nonlinearity of these observed features.While this value would ideally be higher, it is difficult to estimate due in part to the periodicity of the domain and the presence of other wave components superposed with the nonlinear response.Additionally, the BDO equation assumes a neutral background stability above and below the ducting region while our simulations employed a weakly stable background in order to suppress undesired instability dynamics.
Another explanation for the failure of a nonlinear response to develop lies with the radiation of wave energy.In the shallow domain discussed in Sect.3.1.1,the rigid boundaries confined wave energy within the domain in the region of enhanced stability.To model a reduced amount of available   , 8, 15, 23, 30, 38.energy, the perturbation amplitude of Sect.3.1.1was reduced from 0.14 to 0.01.The reduction of initial perturbation amplitude resulted in an essentially linear response which is consistent with the interpretation that removing energy from the ducting region inhibits nonlinear evolution.In the deeper domain discussed in Sect.3.1.2,energy was observed to leak away from the duct.With less energy at the duct, the resulting non-bore-like response makes sense, especially considering the nonlinear response that does arise when more energy is available at the ducting level at initialization.Some numerical studies were performed which are not reproduced here due to both the unphysical nature of the sta-bility profile and the complexity of the dynamical response.However, it is worth reporting that when the duct was surrounded by small regions of negative N 2 and subjected to the same initialization of Sect.3.1.2,wave energy was confined to the ducting region, radiation away from the duct was precluded, and nonlinear steepening and soliton creation occurred.This result is mentioned only because it suggests one possible path of study that may yield insight into the ducting environments necessary for the formation of bores and because it supports the notion that energy confinement to the ducting layer is necessary for bore formation to occur.This view is further supported by research done on the Morning Glory in Australia, which is a tropospheric nonlinear bore phenomenon.Both observational and numerical studies have demonstrated that while trapping is required to prevent the degradation of perturbation amplitude, it is not a sufficient condition by itself for nonlinear bore generation (Menhofer et al., 1997;Goler and Reeder, 2004).

Velocity duct -deep domain
Thus far our bore generation studies have focused on thermal ducts.Ducting is not governed only by the N 2 profile, however.Rather, it is based on the dispersion relation for m 2 , of which N 2 is only one contribution.It is clear from Eq. ( 2) that large-scale winds can play a crucial role in GW ducting and likely also in bore formation.Additionally, recent observations by Fechine et al. (2009) have demonstrated the importance of velocity ducting for mesospheric bores for the first time.
In order to explore the role of Doppler ducting in bore formation, we employ a mean horizontal velocity profile given by Ann. Geophys., 27, 511-523, 2009 The background temperature gradient was chosen to be the same as for the thermal ducting case in an effort to preclude radiation of wave energy away from the duct.We then perturbed this ducting structure with the same initial perturbation used in Sect.3.1.2.As with the cosine thermal duct discussed in that section, the Doppler duct yields a nonlinear evolution have similarities to, and departures from, those anticipated from the BDO equation.
Figure 10 displays contours of horizontal velocity and potential temperature at two different times which clearly exhibit strong nonlinearity and the evolution of an undular bore structure.Nonlinear steepening and crest creation are displayed more clearly with traces of the horizontal velocity field through the center of the domain in Fig. 11.In addition, we see that the leading edge of the bore exhibits a significant kink.
To assess whether the kink in the bore structure seen at later stages in Fig. 11 is physical, we performed two additional simulations: one in a domain (and with an initial perturbation wavelength) twice as long; a second with a domain twice as long, but with the same initial perturbation wavelength (thus two waves in the computational domain).The results for the longer initial perturbation are shown in Fig. 12 for the center half of the domain for comparison with the results in Fig. 11. Figure 13 shows the full domain.The latter reveals that bores scales that arise for the longer initial perturbation are slightly smaller than those for the perturbation with half the initial wavelength, suggesting that the scales may be more dependent on the perturbation amplitude and character of the duct than on the initial wavelength, though  additional work is needed to define the dependence of these dynamics on environmental factors and initial perturbations more quantitatively.
Results for the longer computational domain having the same initial perturbation wavelength exhibit the same structure shown in Fig. 11 and discussed above, indicating that the kink in this case is not an artifact of the limited computational domain, thus likely a consequence of the small horizontal wavelength of the initial perturbation.

Dual ducting environment -deep domain
The likely ducting environment in the MLT is largely defined by low-frequency GWs at smaller vertical scales.These GW motions induce significant velocity and stability variations across the GW phase.These result in potential Doppler and thermal ducts (at the velocity and stability maxima) that have relative vertical displacements that depend on the orientation of a perturbation relative to the direction of propagation of the GW.To examine the possible influences of superposed Doppler and thermal ducts, we examine two cases having both types of ducts, but with different degrees of separation.The two cases considered are shown in Fig. 14.In each case a thermal duct is centered in the domain, and a Doppler duct is offset by either a full or half duct width.The perturbation in each case is centered on the thermal duct, and the non-dimensional perturbation amplitude A is 0.28, twice that used for most other runs.The horizontal domain and perturbation wavelength are 20 and the vertical domain is 18, as in Sect.3.1.2.
Figure 15 displays the resulting waveforms in the potential temperature field for the wider duct separation.The results exhibit a significant asymmetric response that varies in time, suggestive of energy exchange between the two ducts throughout the evolution.The response is also both complex and strongly nonlinear, with clear wave steepening and the evolution of structures at smaller horizontal scales.This suggests that bore evolution can also occur in such dualducting environments, though the structures and evolutions are harder to interpret in terms of the predictions of the BDO equation or the single duct responses discussed above.
Figure 16 displays the corresponding evolution when the thermal and Doppler ducts overlap to a greater degree.
Again, nonlinear steepening and a transition to smaller scales play an important role here.Additionally, the overlapping character of the Doppler and thermal ducts appears to allow much closer and stronger coupling of the responses at, and the flow of energy between, the two ducts.These evolutions depart in significant respects from the responses to simpler ducting structures observed in the MLT and described above.As such, they may not be as readily identified in MLT observations, even if they occur frequently.Thus, further numerical studies of such environments, and more quantitative characterization of the specific wind and stability profiles accompanying observed bore-like structures will be necessary to assess their more general roles in realistic environments.

Conclusions
We performed direct numerical simulations of bore generation and evolution employing a high-resolution spectral model solving the incompressible Boussinesq Navier-Stokes equations in two dimensions.Our solutions demonstrated the viability of numerical modeling in describing mesospheric bore behavior, both in examining their generation by finiteamplitude perturbations at longer initial wavelengths and in defining their responses to varying ducting environments.We found general agreement between our model and the predictions of BDO theory for a shallow duct with rigid boundaries close to our ducting structure, with an estimated value for aλ/h 2 of ∼5.We also found agreement with mesospheric observations, with a crest velocity of ∼80 m/s for a wavelength of ∼8 km and a duct of FWHM of 6 km.Simulations employing a deeper domain that isolated the duct from the boundaries resulted in a response that failed to lead to steepening and bore formation, indicating that caution is needed in interpreting results when artificial boundaries are in close proximity to the ducting level.Undular bore responses were nevertheless able to be excited with the deeper domain with sufficiently strong initial perturbation amplitudes and sufficiently narrow thermal and Doppler ducts.In particular, Doppler ducts were found to yield very similar bore evolutions to thermal ducts, and the occurrence of both Doppler and thermal ducts in close proximity led to responses that were both highly nonlinear and much more structured and complicated than seen to arise for a single duct of either type.Indeed, bore generation and evolution for adjacent or overlapping Doppler and thermal ducts was found to be sufficiently complex and variable in time that BDO theory provides no useful guidance and such events in the MLT may not be recognized as bore responses.Particularly when the Doppler and thermal ducts overlapped to a significant degree (as we expect to occur frequently in the MLT, given the highly structured environments due to lowfrequency GWs), responses included rapid energy exchanges between the ducts that may have interesting implications for observations and other small-scale MLT dynamics.
The results presented here have served to demonstrate that bores can occur easily in response to various MLT environments and long-wave initial perturbations.Our results have also provided tantalizing evidence of the possible range of nonlinear responses in more complex environments, including Doppler ducts, than explored previously.Nevertheless, our present results have only begun the exploration of the factors, including generation mechanisms and propagation environments, that likely influence bore generation and evolution in the MLT.Significant uncertainties accompany the need for trapping wave energy at a ducting level, the range of potential source mechanisms, the dependence of bore structure and behavior on environmental and initial conditions, and the potential complexity of responses arising from these influences.It is not known, for example, whether bores arise only from trapped long-wave perturbations or whether they can also result from sustained long-wave propagation from another altitude.Numerical studies of the Morning Glory (Goller and Reeder, 2004) demonstrated that simply having a duct of high stability is not a sufficient condition for nonlinear responses, and wave radiation must be taken into account.We also have yet to explore other possible generation mechanisms, among them nonlinear interactions accompanying vertically-propagating GWs that have been suggested to lead to linear ducting responses (Chimonas et al., 1996;Snively and Pasko, 2008), GW breaking resulting in the generation of additional GWs at other spatial scales and frequencies (Fritts et al., 2003(Fritts et al., , 2009a;;Snively and Pasko, 2003), and the local body forces accompanying GW breaking.Finally, there is likely a broad range of environmental factors that will impact bore evolution, structure, and complexity which remains almost entirely unexplored to date, but which may prove necessary to understanding the range of responses at MLT altitudes.

Fig. 1 .
Fig. 1.Thermal ducting profiles used in Sect.3.1.1.(a) Dashed line is a plot of the hyper-Gaussian function employed by Seylor (2005); solid line is for our profile composed of hyperbolic tangent derivatives, hereafter referred to as our composite profile.(b) The solid line is our composite profile and the dashed line is the cosine ducting stability profile.(c) The temperature profiles that yield the stability profiles in (b).

Fig. 3 .
Fig. 3. Line plots through the contour plots depicted in Fig. 2. All cuts were taken at z=−0.34.The lowest trace corresponds to panel (a) of Fig. 2. Each trace has a constant offset added to it for clarity ending with the top trace which is a cut through panel (f).

Fig. 5 .
Fig. 5. N 2 (z) (left) and potential temperature profiles (right) for the composite (solid) and cosine (dashed) ducts employed for bore simulations in a domain 9 times deeper than shown in Figs. 3 and 4.

Fig. 7 .
Fig. 7. Line plots of cuts taken through the contours of Fig. 6.Cuts are taken near z=−0.36though vary slightly from panel to panel to demonstrate maximum response.

Fig. 10 .
Fig. 10.Contour plots of both potential temperature and horizontal velocity for the deep domain initialized with the Doppler duct depicted in Fig. 10.Panels (a) and (c) depict potential temperature contours at two different times during the evolution and panels (b) and (d) depict the coincident horizontal velocity field.

Fig. 11 .
Fig. 11.Cuts taken at z=0 through the horizontal velocity field obtained from the deep domain initialized with a velocity duct.Note that 2 of these cuts are from panels (b) and (d) of Fig. 11 and the rest are from contours not shown.Times are 0, 8, 15, 23, 30, 38.

Fig. 12 .
Fig.12.Horizontal velocity plots taken from a deep domain run initialized with a velocity duct for which the horizontal domain and wavelength were increased.Note that this is a subsection of the horizontal computational domain in order that direct comparison with Fig.12may more easily be made.Times are38, 47, 56, 66,  75, 83.

Fig. 13 .
Fig. 13.The same data as in Fig. 12 except shown for the full horizontal computational domain.

Fig. 14 .
Fig. 14.Dual ducting structures utilized in Sect.3.3.Solid lines depict the stability profile while the dashed lines in each case indicate the velocity ducting structure.(a) Velocity duct displacement is equal to the width of the duct.(b) Velocity duct displacement is equal to the half width of the duct.

Fig. 15 .
Fig. 15.Contour plots of potential temperature for the evolution of our initial perturbation in the presence of the dual ducting structure presented in Fig. 15, panel (a).Times are 32, 45, 69, 77, 83, 90.

Fig. 16 .
Fig. 16.Contour plots of potential temperature for the evolution of our initial perturbation in the presence of the dual ducting structure presented in Fig. 15, panel (b).Time are 23, 32, 53, 58, 64, 69.