Modeling formation of turbulent Sporadic-E clouds using realistic wind data

A high resolution two-dimensional multi-fluid model of sporadic-E layers was developed and driven with physically realistic mesosphere, lower thermosphere (MLT) winds measured over Albuquerque, New Mexico. The realistic E-region winds are produced by the HYdrodynamic Point-wise Environment Reconstructor (HYPER) model that ingests meteor derived wind observations from a Spread-spectrum Interferometric Multistatic meteor radar Observing Network (SIMONe) system combined with the Navier-Stokes equations to provide high resolution three-dimensional wind fields over time. Sporadic-E dynamics are simulated using both realistic winds from HYPER as well as idealized hyperbolic tangent windshears to compare and contrast. Overall, the model shows greater inhomogeneity and irregularity using realistic winds with no clear peaks in the spectra, unlike the periodic density structures from the idealized windshears. Furthermore, range-time-frequency (RTF) observations from a local ionosonde were used to compare sporadic-E observations with the simulations. In general, the simulations show Kelvin-Helmholtz billow formation during the periods with range-spread sporadic-E in ionosonde observations, but the ionosonde virtual heights are 5–10 km above the simulated peak densities, likely due to altitude limitations from meteor radar observations. Ultimately, the use of realistic winds to drive sporadic-E models provides more insight to study complex dynamics and decipher observations of turbulent layers.

Sporadic-E (denoted Es) is a strong enhancement of metallic ions in the E-layer of the ionosphere that is known to disrupt High Frequency (HF) skywave propagation (McNamara, 1991; Fabrizio, 2013), Very High Frequency (VHF) amateur radio (Neubeck, 1996), and cause Global Navigation Satellite System (GNSS) loss-of-lock and navigational errors up to meters in magnitude (Dubey et al., 2006; Yue et al., 2016). Populations of long lived metallic ions, deposited by meteor ablation, are convected by upper atmospheric winds and concentrated into denser structures by windshears combined with Earth’s magnetic field induced Lorentz forces, forming ion clouds with densities often greater than the background E-layer (Whitehead, 1989). The windshear model accounts for the tidal winds in the mesosphere and lower thermosphere (MLT) that provide vertical ion convergence (VIC) for zonal and meridional windshears (Mathews, 1998), which must outpace diffusion for Es layers to increase in density (Haldoupis, 2011). This model of Es formation is capable of explaining many global trends when combined with knowledge of diurnal and semi-diurnal tides (Chu et al., 2014; Yamazaki et al., 2022; Hodos et al., 2022), but does not account for small-scale structure or turbulence that is known to form in Es layers due to instabilities (Hysell et al., 2009, Hysell et al., 2012; Hysell et al., 2016).

MLT instabilities have been observed with a wide range of instruments over the last several decades. Quasiperiodic (QP) radar echoes from the Middle and Upper Atmosphere (MU) radar in Japan were believed to be caused by gradient drift instabilities related to sporadic-E layers (Yamamoto et al., 1991). These QP echoes have since been observed at a variety of other locations (e.g., Hysell and Burcham (1999) in the eastern United States, Chau and Woodman (1999) over Piura Peru), indicating that the underlying structures are common features of mid-latitude sporadic-E and likely due to shear instabilities such as the Kelvin-Helmholtz (K-H) instability (Larsen, 2000). In addition to radar observations, these instabilities have been observed using ionosondes (Hussey et al., 1998), airglow imaging (Hecht, 2004), rockets (Larsen et al., 2005), lidar (Ejiri et al., 2019; Hecht et al., 2021), and Long Wavelength Array (LWA) radio telescopes (Obenberger et al., 2025). Models of the instability formation in neutral gasses have been developed to provide insight into driving conditions, structure, and temporal dynamics. A Richardson number less than 1/4 indicates an unstable condition where K-H billows can grow, and the Reynolds number also affects the strength and depth of K-H billows (Fritts et al., 2014). Three-dimensional high resolution simulations were able to capture complex dynamics observed in airglow imagers, including tubes and knots within K-H billows (Fritts et al., 2021). The K-H billow simulations were shown to produce a Kolmogorov spectrum with a broad range of scale sizes (Werne and Fritts, 1999).

Specifically for sporadic-E layers, plasma dynamics for instabilities have been modeled using idealized windshears by several previous studies. In Bernhardt (2002), a two dimensional linearized perturbation model was developed for the growth of a K-H eigenmode using a hyperbolic tangent (tanh) neutral windshear. This time-dependent quasi-neutral plasma transport model showed K-H billow formation with ions forming concentrated vortex centers and hollow cells at low and high E-region altitudes respectively. A two-fluid model was demonstrated in Knisely and Emmons (2023) using similar conditions and windshears as Bernhardt (2002), resulting in K-H billow formation that was shown to produce a time-dependent perturbation on GNSS Radio Occultation (RO) signals as the K-H billows developed and evolved. A three dimensional spectral model using incompressible Navier-Stokes for neutral dynamics and extended Magnetohydrodynamic (MHD) for plasma dynamics was developed in Hysell et al. (2024) to study Ekman spiral wind patterns, which contain both zonal and meridional components that rotate as a function of altitude. Driving the model with an idealized Ekman spiral neutral wind pattern for initial conditions, the simulations were shown to match well with mid-latitude patchy sporadic-E observations from coherent scatter radars, which show QP echoes with secondary irregularities (braids).

While the K-H instability modeling efforts for sporadic-E layers have shown success using idealized neutral wind shears (hyperbolic tangents), recent advances in multistatic meteor radars (MMR) have greatly improved the spatiotemporal coverage and resolution of MLT wind observations at sporadic-E altitudes (Chau et al., 2019). Spread-spectrum Interferometric Multistatic meteor radar Observing Network (SIMONe) systems are capable of measuring over 100 k meteor counts daily, and the derived wind patterns provide insight into dynamics behind nonlinear responses such as mesospheric bores (Ramachandran et al., 2023; Ramachandran, 2023). While the MMR observations drastically improve coverage and resolution, they are not able to fully resolve submesoscale dynamics and require regularization. Recently, a physics-informed neural network approach was implemented to develop HYPER (HYdrodynamic Point-wise Environment Reconstructor), which uses MMR data along with the Navier-Stokes equations to provide physically realistic high resolution 3-D wind fields over time (Urco et al., 2024).

This work uses HYPER wind fields from SIMONe measurements in Albuquerque, New Mexico to simulate K-H billow formation in sporadic-E layers using realistic wind data during sporadic-E events. The model is also driven with idealized hyperbolic tangent wind shears to compare and contrast the predictions and power spectra. Further, observations from a local ionosonde are used to provide insight into the actual sporadic-E dynamics over time as an initial validation of model predictions of K-H instability formation.

The model in this study calculates the advection of ions within a dynamic background neutral wind. It uses compressible Navier Stokes for the dynamics of the neutrals, and an equilibrium drift-diffusion model for the ion advection, incorporating drag due to neutral collisions and Lorenz forces due to the magnetic field in a conductivity tensor (Anderson, 1995; Kelley, 2009). To improve readability, the details of the model are provided in Section 6, in the Supplementary Material. Here, we present the physical parameters and realistic wind data used in the simulations.

E-region properties

The specific physical properties used in the simulations at an altitude of 100 km are given in Table 1. Magnetic field parameters are estimated using the International Geomagnetic Reference Field (IGRF; Alken et al., 2021), and atmospheric neutral densities and temperatures are from the US Standard Atmosphere 1976 model, in the altitude segment starting from 86 km geometric altitude (NOAA, 1976). At 100 km, the Hall parameter for ions, the ratio of their gyrofrequency to the collision rate with the background neutrals, is much less than one, indicating that they are mostly unmagnetized and will tend to drift with the background neutral flow. Electrons are already magnetized in the E-Region, while ions require higher altitudes to become magnetized (Rishbeth and Garriott, 1969).

In this iteration of the model, iron ions (F e +) are advected by collisions with the background neutrals and also experience a Lorenz force expressed in the ion mobility tensor. A mostly zonal windshear that crosses the mostly meridional magnetic field leads to a force that compresses layers of ions into denser sheets. Turbulent structures in neutral winds then wrap the ion sheets into denser structures such as the clouds and quasi-periodic echo “rod” structures. A similar mechanism was also reported by Larsen (2000) and Pan and Larsen (2000).

Initially, the model was used without reference to chemistry, especially with regards to recombination and charge exchange with background ions. Initial ion/electron layers formed and concentrated at the altitude of the largest windshears, which were often significantly below layers observed in ionosonde data. However, ion lifetimes decrease exponentially at lower altitudes, such that chemistry plays a vital role in sporadic-E formation. From the iron ion lifetime dependence as a function of altitude calculated by MacDougall et al. (2000), a simple lifetime model for ions k recomb,s was developed to account for recombination at lower altitudes. In this model k r e c o m b, F e + is the 1 / e recombination lifetime in hours (hr) for iron ions (F e +) and z is the altitude in km. The two parameters of Equation 1 are linearly fit to the 1 0 5 c m − 3 series in MacDougall et al. (2000), Figure 3. This series corresponds to a plasma frequency of 2.8 MHz, which is within the sensitivity of ionosonde observations for sporadic-E. In the model, the ion density decays as a function of altitude at the modeled rate in the ion density conservation equation. F e + (as opposed to Mg + or Na +) were the predominant focus of this paper due to the known role in sporadic-E at lower altitudes (Roddy et al., 2004) and availability of the ion lifetimes as a function of altitude.

The addition of the ion lifetime corrects model behavior and limits the appearance of significant sporadic-E concentrations to altitudes above 90 km.

Realistic wind and ionosonde datasets

One challenge in modeling the MLT region is the sparseness of wind data available for altitudes above 90 km. While previous studies have combined ionosonde with limited satellite mission and sounding rocket observations, fully 4-D wind-fields (three spatial and time) coordinated with ionospheric data have not generally been available (Durre et al., 2006; Swinbank and Ortland, 2003). Incoherent scattering radars such as Aricebo can provide ion velocities, which are about the same as neutral wind velocities below 110 km, but these datasets are from relatively rare ISR installations (Gong and Zhou, 2011). With the recent installation of a SIMONe system in Albuquerque, New Mexico, which also contains a local ionosonde, a combination of realistic wind observations and ionospheric conditions is now available. The ionosonde operated by the Air Force Research Laboratory at Kirtland Air Force Base (35.0 ° N, 253.5 ° E) is a Lowell Digisonde Portable Sounder 4D (DPS4D) (Reinisch et al., 2008) operating at a cadence of one ionogram every 2 minutes. The SIMONe radar is a Spread Spectrum Interferometric Multi-Static Meteor Radar, located in central New Mexico (SIMONe-NM). It operates between 28 and 50 MHz and uses phase-coded spread-spectrum continuous wave (CW) signals to measure the Doppler signal of meteors. Currently, SIMONe-NM uses six transmitters and four receivers spread throughout New Mexico (Obenberger et al., 2026).

Two days from June 2024 were selected for analysis due to strong sporadic-E activity indicated by manual inspection of the Kirtland ionosonde observations and the availability of quality SIMONe data. Since only 2 days are analyzed here, the impacts of seasonal, geomagnetic, and solar cycle variations cannot be determined. The SIMONe observations of meteor velocity were analyzed with HYPER reconstruct dynamically consistent and physically realistic wind fields by combining the measurements with Navier-Stokes constraints. HYPER produces high resolution physically realistic wind estimates in a 300 × 300 km 2 region with an altitude range of 75–107 km. The data have a horizontal sampling of 20 km, a resolution of ∼ 40 km, a vertical resolution of 0.5 km, and a time resolution of 30 min. Figure 1 shows the three components of the zenith-pointing wind over time for 12 June 2024 above Albuquerque.

Ionograms for the month of June were scanned for dates and times of interest where sporadic-E layers formed. Sporadic-E is indicated on the ionogram as an unusually flat return at E-region altitudes, generally exceeding the background peak E-region frequency, foE (Piggott and Rawer, 1961). Figure 2 shows a quiescent period and one of the sporadic-E periods on 12 June 2025. Range-time-frequency (RTF) graphs are constructed from ordinary-mode returns of the collected ionograms, which closely correspond to the electron density as a function of altitude, and are shown in Figure 3 for 12 and 14 June 2024. For each ionogram RTF datafile, pynasonde, an open source python library (Chakraborty et al., 2026), was used to read in the O-mode frequency data. The maximum reflected frequency was binned as a function of virtual altitude (1/2 the group path) from the O-mode channel and plotted as a function of altitude and the ionogram time. As observed in the RTFs, there are several occurrences of “V” shaped sporadic-E returns over these 2 days, indicating the presence of billowy sporadic-E layers. The V-shaped layers are produced by off-zenith backscatter from broken ion clouds drifting over the ionosonde site (see Obenberger et al. (2025) for more details). Large range spreads are also indicative of billowy turbulent layers, whereas flat layers (e.g., 1800–2,300 on 12 June 2024 in Figure 3) imply the absence of billows.

Wind data for periods containing sporadic-E were searched for zonal sections containing high windspeeds (greater than 50 m/s) and high vertical windshear gradients (greater than 40 m/s/km). Winds from these times having the simultaneous occurrence of sporadic-E, high windspeeds, and shears were used as initial conditions for the simulations. Details of selected regions and times are provided in Table 2. Figure 4 shows the Richardson numbers from the 3-D wind data for 12 June 2024. The minimum Richardson number indicates the strongest shear layer present in the wind gradients above 83 km altitude for the entire 300 × 300 km 2 horizontal domain. Richardson numbers consistent with buoyancy opposing instabilities are present simultaneously with observed sporadic-E layers, but the magnitude is not proportional to the observed turbulence of the sporadic-E layer. More discussion of this will be given in Section 4.

Idealized hyperbolic tangent windshear

An idealized hyperbolic tangent (tanh) windshear with small initial density perturbations was used as a point of comparison with prior modeling efforts. This symmetric initial condition produces regular Kelvin-Helmholtz billows which compress the ions into sheets, then wrap into “rod” structures similar to the quasi-periodic echoes reported in Urbina et al. (2000).

For the geometry outlined here, x ̂ points to the east (right), y ̂ points north (into the page), and z ̂ points vertically upward (up). The initial conditions were chosen to be similar to the initial conditions in the 100 km scenario of Bernhardt (2002) shown in Equations 2–5 below: where d 1 = 250 m, u 0 = 100 m/sec, ρ 0 = 6.20 × 1 0 − 7 kg/ m 3, d 2 = 5688.83 m, L x = 9000 m, L z = 500 m, z mag = 250 m, T = 182 K, and | B | = 3.7 × 1 0 − 5 T. In the simulated tanh windshear condition, the magnetic field is pointed out of the page (− y ̂). These conditions correspond to a Richardson number of 0.01, where anything less than 1/4 is an unstable condition. The wavelength of the most unstable mode should be 8 times the width of the shear layer (Larsen, 2000), or 4 km, although in this case the initial 9 km perturbation appears to grow to dominate the flow.

The initial metallic ion layer was assumed to start as a gaussian sheet with a given peak background density and thickness. Wind velocities crossed with the magnetic field lead to compression or dispersion of the ion layer. The relative orientation of the windshear to the magnetic field vector governs whether a given shear will compress ions into rods in the center of the K-H billows, or disperse them and possibly concentrate them in layers outside the vortices (Figure 5). The neutral wind develops the Kelvin-Helmholtz instability, producing vortices (Figure 6). The compressed ion sheet wraps up in the vortices and compresses further into denser rod shapes. These rod shapes resemble the quasi-periodic echo behavior seen in Figure 4 of Urbina et al. (2000). Werne (2025) has described these features as due to the Lorentz Force, and has performed similar modeling with similar results using a Direct Numerical Simulation (DNS) spectral method code.

Bernhardt (2002) modeled the advection of ions and electrons in a linearized perturbative model of a growing Kelvin-Helmholtz instability. Periodic rod-structures were produced at 100 km altitude (Figure 15 in Bernhardt (2002)), however, the different physical assumptions also result in different behavior. For example, the enforced quasi-neutrality of the plasma during evolution also produced some two-layer structures near 120 km altitude, where electron magnetization plays a key role in plasma dynamics. The physics in the model developed for this study follows ion advection under neutral wind and Lorenz forces, and assumes that the electrons will follow ion populations.

When the recombination model is implemented for the ideal windshear case, the total ion density attenuates over time, leading to a lower ion density (note the different scales in Figures 5, 7). The dispersing case is more strongly affected by differential recombination rates because of the larger altitude spread of the ion concentrations. Surprisingly, because of the recombination timescales, the layers are denser as they pass beneath the K-H billows and more attenuated above the K-H billows, attenuating incompletely in their passage below the billow. The recombination 1 / e timescales are 378, 859, and 1953 s at 93, 95, and 97 km, compared with a K-H billow period of ∼ 150 s.

Realistic winds initial condition

In addition to the idealized case presented in Section 3.1, the real 4-D wind datasets were used as initial conditions for the model. Realistic windshears can resemble the simple case of the tanh windshear, but can also have additional structures that include meridional components of the shear, some vertical motion, and multiple layers. Vertical ionograms for the month of June 2024 were scanned for sporadic-E layers, and the 3-D wind data for dates, times, and locations with zonal winds and wind-shear magnitudes exceeding 50 m/s, and 40 m/s/km. Regions where both occur are tabulated and described in Table 2.

For these cases, the declination of the magnetic field was taken into account. 2-D regions were selected such that the y-axis points along the magnetic heading, so that a vector parallel to the surface with the same declination is perpendicular to the plane, pointing into the page. The magnetic field points primarily towards the Earth’s surface over Albuquerque, with an inclination angle of 62° that is taken into account in the simulations. Additionally, recombination is taken into account via an effective ion loss term as a function of altitude (Equation 1).

As displayed in Figure 3 on 14 June at 0530 UTC (Sim14), a billowy sporadic-E pattern is present in the ionosonde observations with a large range spread extending from ∼ 105–150 km. When searching the 3-D wind dataset, a complex multi-layer wind pattern was observed, with multiple layers (Figure 8a). A gaussian ion layer was initialized (biased towards the top of the domain by 4 km), and a simulation was run with these initial conditions. A region 80 km in length, 85–107 km in altitude was selected from the 3-D wind data. The behavior of the modeled ion density over time is shown in Figures 8b–d.

On 12 June at 1500 UTC (Sim12pm), a relatively flat sporadic-E layer was present in the ionosonde returns (Figure 3) that was dense and long-lived. As shown in (Figure 9a), the 3-D wind data indicates two strong windshears, one at 85 km and another at 100 km. In simulations conducted prior to the introduction of the recombination chemistry (not shown in this paper), the 85 km windshear entrained ions causing peak sporadic-E density at unusually low altitudes. Accounting for Fe + lifetimes (Equation 1), metallic ions convected below ∼ 90 km rapidly recombine. The ion simulations over time (Figures 9b–d) show a set of K-H billows at 100 km that appear to convect the ions, which become dense in the middle of the 95 km, 90 m/s wind layer.

On 12 June at 0710 UTC (Sim12am), another complex, layered wind formation is seen in the 3-D wind data shown in panel a of Figure 10. A sporadic-E layer formed with a lower density than the 14 June 0530 UTC scenario (Figure 3), although the layer still appeared to have a similar billow-induced range spread. The simulated temporal ion dynamics for this scenario using realistic wind inputs are displayed in Figure 10, showing another relatively turbulent neutral wind pattern that advected the ions almost through the accessible domain, if not for the recombination.

In comparison, the simulated ion densities using realistic winds show more inhomogeneities and irregularities than the idealized tanh windshears. Although both sets of simulations produce K-H billows, the idealized windshear produces billows at regular spatial intervals of about 8 km while the realistic wind simulations do not show regular intervals between the billow structures. Furthermore, realistic wind simulations no longer show roughly symmetric spiral shapes, but instead have a more chaotic and turbulent structure. The scale of the turbulent structures present in the simulations is governed by instabilities growing from the initial conditions and the simulation grid sizes on the order of meters. The initial data has a horizontal resolution of 20 km, and is only the seed of growth for this finer structure. Smaller scale structure develops from the initial more uniform flow as momentum and viscosity move kinetic energy from larger to smaller scales, eventually terminating in the viscosity dominated Kolmogorov lengthscale. Qualitatively, the asymmetric structures from realistic wind simulations align closer with QP echoes that show more irregularity than the constant period predicted by the idealized windshear scenario (Chau and Woodman, 1999; Hysell and Burcham, 1999).

A quantitative comparison of length-scales present in each simulated ion cloud is given by the scaled horizontal Fourier transform of the ion density shown in Figure 11. The Fourier transform of the ion density within the model domain is taken along the x-axis and the absolute value is averaged across the middle third of altitudes to show lengthscale detail. For ideal tanh windshears, the ion density pattern is periodic such that every other Fourier component contributes significantly, the most important being the horizontal length scale of K-H billows at 9 km. Significant details are present down to about 2 km. The realistic wind scenarios are more chaotic and show less periodicity at a single length-scale. For Sim14, the largest length-scale is 20 km, with detail down to about 3 km (where the FFT components fall below about 0.2 on this normalized comparison scale). For Sim12pm, the length-scale range is about 20 km down to 8 km; for Sim12am, 11 km down to about 5 km. The second plot shows a power spectrum averaged across all altitudes at several times. The power spectra show a scaling law that appears to be common across scenarios: ion density features appear to decay with a k − 1 dependence, where k is the spatial frequency.

Comparison with ionosonde observations

Ionosonde RTF observations are compared with the Sim14 (0530 UTC, 2330 Mountain Daylight Time, MDT), Sim12pm, and Sim12am cases depicted in Figures 8–10. Closeup comparisons are given in Figures 12–14. The vertical O-mode ionosonde returns are closely related to the plasma frequency due to reflections that occur at altitudes where the signal frequency equals the plasma frequency (Davies, 1990). Therefore, the ionosonde RTF observations provide an estimate of ion dynamics as a function of altitude (range) and time to compare with the simulations. Given the relatively low vertical resolution of the ionosonde observations, the comparison is based on the altitude range spread in the ionosonde data to the range spread of the features in the simulated ion clouds. The two early morning initial winds develop turbulence which advects a simulated ion layer into broken concentrations spread out widely in altitude. The corresponding ionosonde data shows a V-shape, which indicates the approach and recession of a concentrated ion cloud, and a “puffier” larger range-spread in altitude. The sim12pm case develops a broader flatter ion layer that is concentrated on the underside of a shear layer or less-turbulent billow structure. The corresponding ionosonde data is of a longer lived, unbroken and mostly flat ion density profile. Our wind data only exists where meteors begin to ablate up to about 105 km, and in addition, the resolution of both wind data and ionosonde data is coarser than the features we would like to simulate. Despite the limitations of the data, we suggest the behavior of the simulations may correspond to some of the features observed in the ion density patterns.

Note that power spectra were not calculated for ionosonde RTF observations due to the limited vertical resolution of ∼ 5 km that cannot resolve the small spatial scales of interest. However, range-spread in ionosonde sporadic-E measurements implies the presence of small-scale structures due to signal reflections from a variety of locations and altitudes for a single transmit frequency, and V-shaped returns (as observed in Figure 12) indicate off-zenith backscatter from layers drifting over the ionosonde site (Obenberger et al., 2025).

The relatively complicated initial windshears present for Sim14 and Sim12am form highly turbulent structures that convect ions throughout the simulated volume. Recombination reduces ion lifetimes at low altitudes, although some short-lived concentrations appear to convect below 90 km. In the ionosonde RTF observations, a strong return over a relatively wide range spread appears over the course of 30 min and lasts for another 30–60 min. Our model implies a connection between the complicated multi-layer initial winds, leading to highly turbulent flow that advects the ion density over a wider altitude range even as it concentrates them into denser more strongly reflective features. In Obenberger et al. (2025), the V-shaped returns in RTF ionosonde measurements are described in detail and compared against Long Wavelength Array (LWA) observations. These V-shapes appear somewhat similar to the puffy altitude range-spread observed in the 12 and 14 June early morning cases (sim12am and sim14). Obenberger et al. (2025) describes these observations as indicating a concentrated density that approaches and recedes in apparent range as the cloud passes over the Digisonde, which would be characteristic of a more broken density pattern. Either explanation, varying range to receiver or physical altitude range-spread, may be responsible for the apparent altitude spread of these early morning cases. Both likely play a role given the broken ionosonde traces.

For Sim12pm, a relatively dense, long lived, and flat layer appears in the ionosonde RTF at an altitude of 105 km, with only about a 5 km thickness (the resolution limit for the ionosonde range reporting). In the simulations, a series of K-H billows from 95 to 105 km wrap and concentrate ions into sheets at about 95–100 km altitude. The simulated sheets appear to have a thickness of about 2.5 km. The model has a single strong initial shear layer that develops less complicated turbulence. The relatively flat Sim12pm virtual height likely indicates the presence of a broad layer concentrated by a shear that is either not turbulent or does not strongly advect the ions away from the principal layer as it concentrates them.

The mean altitude indicated by the ionosonde RTF observations is ∼ 110 km in the Sim14 scenario, ∼ 105 km for Sim12pm, and ∼ 100 km for Sim12am. In the simulations, the mean altitudes associated with the enhanced density regions are approximately 95 km (Sim14), 100 km (Sim12pm), and 100 km (Sim12am). Part of the altitude bias is due to the lack of meteor observations above ∼ 105 km which results in no reconstructed winds above 107 km. Important wind structures are likely present above this altitude, but are not measurable with current meteor radar sensitivities. However, if the multi-layered character of the winds persists above the range of sensitivity, the multi-layered initial conditions that evolve in turbulence may explain the characteristics of these high-altitude-spread clouds.

Furthermore, ionosonde virtual heights are greater than actual heights due to the group delay of high frequency signals through the ionospheric plasma (Rishbeth and Garriott, 1969). For sporadic-E layers, this positive virtual height bias varies as a function of altitude, ranging from zero bias for actual layer altitudes at 100 km to ∼ 5 km for layers at 110 km (see Figure 5 of Haldoupis et al. (2024)). However, the virtual height bias also depends on the transmitted signal frequency, where higher radio frequencies relative to the background E-layer plasma frequencies result in closer agreement between sporadic-E virtual and actual heights due to indices of refraction closer to unity. From this, it appears that the main cause for altitude discrepancies between the simulations and the ionosonde observations is due to the lack of wind observations above ∼ 105 km.

This two dimensional model has been used to simulate ion advection within a background neutral wind, using assumptions of equilibrium drift, with the dominant force being from collisional drag on the ions and the Lorenz force, expressed in the conductivity tensor. These physics lead to behavior for ion layers where zonal velocity shears in the background neutral flow, oriented correctly relative to the meridional magnetic field tend to concentrate initially diffuse metallic ion layers.

The presence of simple shears leads to K-H instabilities when the shear magnitude exceeds 80 m/s/km (corresponding to a Richardson number below 1/4 in this region of the atmosphere). These billows, depending on direction, concentrate the ions further into one-dimensional rod structures, resembling the quasi-periodic echoes observed by Urbina et al. (2000), Larsen (2000), and Wang et al. (2025). The idealized tanh windshear model produces regular rod concentrations in the vortices of the Kelvin-Helmholtz instability. The initial ion layer compresses and then rolls into dense structures in the center of the vortices. This behavior matches the behavior found by Werne (2025) with a Direct Numerical Simulation (DNS) simulation using similar physical assumptions. Bernhardt (2002) also predicted similar behavior in the low altitude case from a linearized perturbation model using different physical assumptions (enforced quasi-neutrality of drifting plasma), resulting in the formation of hollow cells in his 120 km altitude scenario. In contrast, our idealized windshear simulations predict enhanced densities in the billow centers.

In our realistic meteor-radar derived wind data, more complexity is often present (Urco et al., 2024). Multiple shears at different altitudes have been observed, as well as turning shears with zonal and meridional winds layered on top of each other can be seen in Figure 1. Even with the limitations of the 2-D model, the turbulence developing from multiple thin windshear layers advect ions into more irregular structures than the simple ideal windshear scenario. Broken and isolated ion densities are concentrated and advected over large altitude ranges forming billowy structures, such as those shown in sim14 and sim12am. Likewise, the large altitude range spread and V-shapes in the corresponding ionosonde data indicate isolated and broken clouds approaching and receding from the detector, and possible billowy structure. Simpler initial conditions with one or two dominant windshears, such as the sim12pm case, do not result in as much turbulence and the ions are concentrated into unbroken sheets. This agrees with the ionosonde measurements of a flat unbroken layer on 12 June 1530 UTC. Because instability growth depends on vertical windshear, uncertainties in the HYPER-reconstructed shear profiles (on the order of ± 5 m/s) may influence the timing and structure of the simulated features.

The minimum Richardson number, shown in Figure 4 indicates the strongest | ∂ u ⃗ ∂ z | 2 shear layer present in the 3-D wind dataset at an altitude above 83 km. Richardson numbers below 1/4 indicate the potential for instability but the relative magnitude may not be telling the entire story: One key difference between the 0700 and 1530 wind-fields is that at 0700 (sim12am) there are many layered shears, even though they are each individually weaker. At 1530 (sim12pm) there is only one strong shear layer. This could explain why we get puffy broken sporadic-E at 0700, indicative of turbulence, and a flatter layer at 1530. 1530 may be unstable but may be producing a more organized unbroken K-H pattern, and the blanketing layer is a sheet on the underside of the K-H billows.

Although the simplified 2-D model of 3-D winds certainly lacks important contributions from the ignored spatial dimension, the model results using realistic winds provide helpful insight into both spatial and temporal variations that cannot be captured by idealized tanh windshears set at arbitrary altitudes and locations. In addition to improving agreement with QP echoes by predicting more irregularity than idealized windshear simulations, using realistic wind profiles to form sporadic-E layers at the correct altitudes and locations is vital for comparison against observational techniques capable of tracking horizontal motions such as the Long Wavelength Array (Obenberger et al., 2025) or ground-based GNSS receiver networks (Maeda and Heki, 2015). Additionally, GNSS radio occultation (RO) observations are highly sensitive to altitude, length, and density variations, as well as small-scale structure due to turbulence (Emmons et al., 2022; Knisely and Emmons, 2023). Because of this, improvements in modeling physically realistic turbulence developing from initial wind-shear structures can improve our ability to understand and interpret the sort of disruptions we see in RO data. Ultimately, improvements in modeling sporadic-E layers using realistic winds can help to improve interpretations and characterizations of the observations of and disruptions by sporadic-E layers.

Future work to improve the model will proceed along several fronts. The ability to model certain phenomena, such as the gradient drift instabilities, believed by Hysell et al. (2009) and Larsen (2000) to be important in shaping ion clouds, will require a multi-fluid approach where ions and electrons can drift separately and build up resolvable charge concentrations. Taking the lessons learned with this 2-D model and returning to the 3-D GPU model will enable the simulation of inherently 3-D wind situations, such as the one depicted at times in Figure 1 where layered zonal and meridional wind components resulting in a turning shear are present. Conditions such as these were modeled in Hysell et al. (2024) using a quasi-neutral MHD spectral method model of an Ekman spiral. His work produced braided structures that will require a 3-D model to resolve. Finally, the inclusion of additional ion species and a robust plasma kinetics approach will allow modeling of ion chemistry in a more complete manner than the Fe + recombination model used in this paper.

A 2-D model of ion advection in an evolving background neutral wind has been introduced. The combination of a co-located meteor-radar and ionosonde in Albuquerque, New Mexico provides upper atmospheric 3-D winds that can be used as initial conditions for modeling along with local plasma density observations to compare against the simulations. The model predicts physically reasonable behavior consistent with prior work using an idealized hyperbolic tangent initial windshear which develops Kelvin-Helmholtz billows (Bernhardt, 2002; Werne, 2025). Using realistic winds as input, additional structures such as chaotic large-range-spread broken concentrations, develop that are not observed in the idealized windshear simulations.

Comparing against ionosonde observations, windshears that resulted in more turbulent wind-fields can possibly be associated with the greater altitude range-spread observed in ionosonde range-time-frequency (RTF) observations. RTF V-shapes are associated with turbulent and billowy clouds approaching and receding from the ionosonde, and the simulations predict significant small-scale structure from turbulence during these periods. A windshear that resulted in a flatter, less turbulent layer in RTF observations could explain a flatter, broader ion sheet observed in the afternoon of 12 June 2024. The simulated peak densities were 5–10 km below the ionosonde indicated sporadic-E virtual heights. While the virtual height to real height bias may account for a small difference in E-region altitudes, the primary cause of the difference is likely due to the lack of meteor-radar wind observations above 107 km. Further, given the relatively coarse altitude resolution of ionosonde observations, a direct comparison of small-scale structure to the simulations was not possible. However, this type of comparison may provide more useful information in future analyses focused on lower altitude billowy sporadic-E layers.

For future work, the model will be extended into 3-D to include rotational effects as a function of altitude. Selecting a 2-D plane is an oversimplification of the complex 3-D wind patterns, and the impacts of including the third dimension should be explored. Further, the simplified chemistry approach that only includes estimated Fe + lifetimes will be refined to include a more complete plasma kinetics package with additional metallic ion species. Finally, the next comparison will either focus on lower altitude sporadic-E layers well below the 107 km cutoff for wind observations, or will use incoherent scatter radar observations for comparison and validation.