## Abstract

We argue the merits of having a simplified method to simulate equatorial plasma bubble (EPB) plume structures for practical usage. The capability to realistically model EPB plume structures in simulations would be advantageous when assessing the severity of ionospheric threats. Such advantages would arise because a realistic model of EPB plume structures could allow nonstationary scintillation signals to be simulated. Although EPB plume structures can be modeled via first-principle physics-based models, these models tend to be computationally demanding. High-performance computing facilities might be able to offer some remedy, but serious handicaps would remain for those without access to such advanced facilities. We investigated multiple options that utilize the diffusion-limited aggregation (DLA) fractal process to generate bifurcating structures that resemble typical EPB plume structures. We combined the DLA algorithm with the International Reference Ionosphere model to simulate EPBs in three dimensions. Initial tests of this modeling approach indicate promising results.

- equatorial plasma bubbles
- fractal structure
- GNSS
- low-latitude ionosphere
- numerical simulation
- scintillations
- total electron content

## 1 INTRODUCTION

Equatorial plasma bubbles (EPBs) are low-latitude ionospheric phenomena that occur during nighttime hours, primarily during the post-sunset period, driven by the Rayleigh–Taylor instability (see, e.g., Chaturvedi and Ossakow, 1977; Sultan, 1996) at the bottom of the ionospheric F-layer around dusk. The occurrence of EPBs may be observed as equatorial spread-F in ionogram traces and radar echoes, as total electron content (TEC) depletions in global navigation satellite system (GNSS) measurements from ground receiver stations, and as dark bifurcating bands in all-sky airglow observations (e.g., Booker & Wells, 1938; Woodman & La Hoz, 1976; Makela et al., 2006; Pradipta et al., 2015; Hosokawa et al., 2020). Within the EPB-associated TEC depletions, there are ionospheric plasma density irregularities that can cause amplitude scintillations in GNSS as well as in communication satellite (SATCOM) signals. Severe ionospheric scintillations due to EPBs may reduce the availability of GNSS and SATCOM services over low-latitude regions during nighttime hours. Further, steep TEC gradients associated with EPBs may also adversely affect the integrity of space-based augmentation systems (SBASs) and ground-based augmentation systems (GBASs) in civil aviation (see, e.g., Saito et al., 2017; Yoon et al., 2017). The occurrence of EPBs follows a regular seasonality (Burke et al., 2004; Gentile et al., 2006a, 2006b), but there is significant day-to-day variability in EPB occurrence as well (Ogawa et al., 2006, 2009; Joshi et al., 2019), which is far more difficult to forecast (e.g., Li et al., 2021).

Amplitude scintillations due to plasma density irregularities that reside within EPBs depend on several factors, including the background ionospheric density and the wavenumber spectrum of plasma density irregularities (see, e.g., Groves and Carrano, 2016). When the background ionospheric density is high, most EPBs are likely to have deeper TEC depletions and stronger plasma density irregularities. For this reason, low-latitude scintillations tend to be stronger during solar maximum years and significantly weaker during solar minimum years, even if the seasonal occurrence of EPBs remains relatively stable throughout the solar cycle. Meanwhile, the wavenumber spectrum of plasma density irregularities within EPBs also influences the strength of ionospheric scintillations. The wavenumber spectrum evolves over time from the initial growth of the EPB structures until the structures start to decay in the early morning hours. Freshly generated EPBs tend to contain a greater amount of Fresnel-scale density irregularities that can cause scintillations in GNSS and ultra-high-frequency SATCOM signals, whereas decaying EPBs tend to be devoid of such Fresnel-scale irregularities. For practical purposes, scintillations in GNSS and SATCOM signals are often simulated by a phase-screen model (Rino, 1979a, 1979b; Beniguel, 2002; Beniguel & Hamel, 2011), which replaces an extended ionospheric slab of finite thickness with a thin equivalent diffraction screen.

In assessing potential threats posed by EPBs and scintillations, it is also important to note the various stages of the EPB life cycle throughout the nighttime hours, which control the precise nature of the threats. Once the initial development of EPBs successfully passes the initial linear growth phase, the subsequent nonlinear growth of EPBs typically leads to the development of tall plume structures that protrude into the upper region of the ionospheric F-layer. The complete geometry of EPB plume structures can be simulated numerically via first-principle physics-based models such as physics-based model (PBMOD), Sami3 is a model of the ionosphere (SAMI3), or high-resolution bubble (HIRB) (see, e.g., Retterer, 2010; Huba et al., 2015; Yokoyama et al., 2019). Alternatively, semi-analytical fluid models of EPB growth may also be used to simulate the nonlinear development stage (Bernhardt, 2007a, 2007b). In these physics-based models, the introduction of some seed perturbations around the sunset period is typically required to initiate the instability process. Here, one may be forced to apply randomly selected seed perturbations at the bottomside F-layer because of a limited knowledge of the state of the ionosphere at these length scales, which is not ideal. Yet another technical challenge is the computationally demanding nature of numerical partial differential equation (PDE) solvers employed for running these physics-based models.

Despite the aforementioned difficulties, there are distinct advantages that would be gained by modeling the EPB plume structures in simulations. The inclusion of realistic EPB plume structures in scintillation modeling would allow us to simulate nonstationary scintillation signals as the line-of-sight (LOS) between ground-based receivers and orbiting satellites intersects different parts of the EPB plume structures at different times. Further, in the context of technological impacts on GNSS constellations, the fact that Fresnel-scale irregularities are clustered in localized patches also implies that not all GNSS satellites will be affected by scintillations simultaneously. An illustrative example is a study by Rino and Carrano (2011) in which EPB plume structures generated in a PBMOD simulation were used in scintillation modeling. There have also been examples involving simple analytical square-shaped vertical EPB structures (Saito et al., 2009a, 2009b; Harris et al., 2011; Fujiwara & Tsujii, 2016; Supriadi & Saito, 2019), which helped evaluate how potential hazards due to EPBs on a regional ground-based facility may materialize.

Given the advantages of using simulated EPB plume structures in scintillation modeling, it is therefore desirable to address the challenge of a computational bottleneck. Within the context of first-principle physics-based models, there are a number of possibilities that may be pursued, including the use of high-performance computing facilities and the development of better numerical PDE solvers. High-performance computing facilities may be able to offer some remedy to the situation, but a handicap would remain for those without access to such facilities. Likewise, any breakthrough in improving numerical PDE solvers is difficult to predict and cannot be relied upon in the short term. Moreover, with respect to detailed spatial scales, even physics-based models have limitations, as small-scale structures generated in first-principle physics-based models may not necessarily reproduce small-scale structures of real-world EPBs. Considering these constraints, perhaps a more practical possibility in this case is to formulate simplified ways of generating sufficiently realistic EPB plume structures in numerical simulations, achieving a lower computational cost compared with full physics-based models. Deprioritizing or abandoning self-consistent first-principle physics may be considered an acceptable sacrifice as long as the EPB plume structures generated by the simplified model still meet a minimum geometrical resemblance to experimental observations (see, e.g., Kil, 2015; Balan et al., 2018). If this option is indeed feasible, the availability of such a simplified-yet-realistic model of EPB plume structures could potentially open the door to ubiquitous implementation of EPB and scintillation modeling, including by those who have access to only modest computing facilities.

Figure 1 illustrates a potential merit in the implementation of such a simplified-yet-realistic model of EPB plume structures. Because the simplified EPB plume structure model may no longer contain any first-principle physics, the utilization strategy differs considerably from that of full physics-based models. In this case, the simplified model would be used to simulate numerous scenarios of EPB plume structure geometry, which would be highly impractical for full physics-based models. Combined with satellite ephemeris data, the ensemble of modeled EPB plume structures can be used to assess the expected severity of scintillation for different satellite links in each scenario. Such statistical ensemble simulations could be used for characterizing the full range of EPB- and scintillation-related hazards to technological systems that operate in low-latitude regions and depend on links to satellites in orbit. Given the rapid adoption of technology by numerous developing economies in low-latitude regions (Poushter, 2014; Elmi, 2020) and the steadily increasing roles of space infrastructure in our modern digital economy (OECD Space Forum, 2020), a comprehensive characterization of EPB- and scintillation-related hazards will be needed. Thus, there is a good opportunity to fulfill this particular need.

In this paper, we present the construction of one such simplified EPB plume structure model. We begin by discussing the basic geometrical structures of EPB plumes and comparing them with geometrical structures that share similar characteristics yet can be simulated via considerably simpler mathematical methods. Finally, we outline the basic logic/scheme of our simplified EPB plume structure model and present some test results as a case demonstration.

## 2 BIFURCATING PLUME STRUCTURES OF EPBS

Figure 2 depicts an example of EPB plume structures in comparison with two similar geometrical shapes. Figure 2(a) shows a tilted EPB plume structure with some branching patterns, based on airglow observations conducted at Haleakala, Hawaii (Makela & Kelley, 2003). Here, the airglow emissions originated from a horizontal layer at an altitude of approximately 250 km, and they were subsequently mapped onto a vertical plane at the geomagnetic equator following the geomagnetic field lines. Figures 2(b) and 2(c) show other geometrical shapes that can be compared with the EPB plume structures. Figure 2(b) shows a geometrical structure known as viscous fingering patterns, obtained by air being pushed through thick oil while everything is confined in a flat two-dimensional (2-D) space. Figure 2(c) shows a geometrical structure known as diffusion-limited aggregation (DLA) fractal patterns, generated numerically through the random walk (Brownian motion) of many test particles. We can clearly see some similarities between these three geometrical structures, especially in terms of their branching patterns. The geometric similarities point to a realistic possibility of approximating the shape of EPB plume structures with either viscous fingering patterns or DLA fractal patterns.

Whereas the growth and development of EPBs in low-latitude ionosphere are controlled by the Rayleigh–Taylor instability, the formation of viscous fingering patterns in a Hele-Shaw cell is caused by the Saffman–Taylor instability (Saffman, 1986). In contrast, the DLA fractal patterns mimic the growth of dendritic geometrical shapes such as coral reefs and electrical discharges (Witten & Sander, 1981). The DLA mechanism works by having the structure’s elementary building blocks float around as point particles in random walks until they eventually hit the main structure and stop moving. Both viscous fingering and DLA naturally lead to the formation of branching patterns, and the fractal dimensions of the resulting patterns from either mechanism are comparable (e.g., Daccord et al., 1986). To the best of our knowledge, the fractal dimensions of EPB plume structures in the ionosphere have not been investigated in detail. However, by visual inspection of the structures (cf. Figure 2), we may suspect that the fractal dimensions of EPB plume structures could be similar those of viscous fingering and DLA.

To verify and confirm this premise, we investigated the fractal dimensions of the branching structures depicted in Figure 2 by employing computational image analysis software, including ImageJ (https://imagej.nih.gov/ij/) from the National Institute of Health, FDEstimator (http://www.fractal-lab.org/Downloads/FDEstimator.html) from the Fractal Lab, and FDIM (http://reuter.mit.edu/soft-ware/fdim/) from the Laboratory for Computational Longitudinal Neuroimaging at the Massachusetts Institute of Technology. It was found that the EPB plume (cf. Figure 2(a)) has a fractal dimension in the range of 1.71–1.80. Meanwhile, the other branching structures (cf. Figures 2(b) and 2(c)) have fractal dimensions in the range of 1.67–1.85, in general agreement with the underlying literature on DLA (Witten & Sander, 1981, 1983).

Another common morphological feature of EPBs is that the western walls are typically more structured (Tsunoda & White, 1981; Tsunoda, 1981, 1983; Mendillo & Baumgardner, 1982; Yokoyama et al., 2015). Although a regular DLA fractal mechanism is generally isotropic, modified variants of the DLA mechanism also exist whereby directional asymmetry in the branching complexity can be induced (Kondoh & Matsushita, 1986; Meakin et al., 1987; Dupraz et al., 2006; Duarte-Neto et al., 2014). This offers additional flexibility in possible ways of matching real-world EPB morphology with simulated fractal geometry through fine-tuning.

At this point, an important question must be addressed regarding the practicality of modeling viscous fingering patterns and DLA fractal patterns. Modeling of viscous fingering patterns requires computational fluid dynamics (CFD). Although viscous fingering simulations will most likely be less complex than numerical simulations of collisional and magnetized plasma in the Earth’s ionosphere, the CFD calculations can still be computationally demanding if we are required to solve coupled PDEs on a fine grid in the process. Hence, using simulated viscous fingering patterns as a proxy for EPB plume structures may not be effective in practice. In contrast, DLA fractal patterns are relatively straightforward to generate in numerical simulations. The simplicity of generating DLA fractal patterns arises from the fact that the bulk of the computational process only involves the use of a random walk (Brownian motion). Thus, it is more reasonable to adopt DLA fractal patterns as a proxy for EPB plume structures in our simplified model.

The inherent role of the random walk (Brownian motion) in the DLA fractal algorithm, which calls for the use of a random number generator, serendipitously fits the intended utilization of the simplified EPB plume structure model as a tool for statistical ensemble simulations. Realizations from two separate runs of the DLA fractal algorithm, although fundamentally equivalent in terms of their underlying statistical properties, could appear very different in terms of their geometrical branching patterns. Each realization would represent a scenario that can be used to build a larger ensemble. By adopting the DLA approach, the deterministic nature of the EPB plume structure model is sacrificed. However, this loss of the deterministic nature may not necessarily be a cause for deep regrets, as even first-principle physics-based models, whose major strength is their deterministic nature, often require randomly generated seed perturbations for simulating EPB plume structures.

## 3 A SIMPLIFIED OPTION FOR AN EPB PLUME STRUCTURE MODEL

As outlined in the previous section, here we utilize the DLA fractal algorithm to generate branching patterns in the EPB plume structure model. Although the branching patterns generated by the DLA algorithm are almost guaranteed to resemble typical EPB plumes, there are additional technical requirements that must be satisfied by the EPB plume structure model. These requirements regarding EPB plume features are listed in Table 1.

Although these additional requirements can be satisfied in a number of ways, we apply the simplest possible route in this work. Requirements regarding the initial and maximum altitude of EPBs can be satisfied by strictly defining the spatial domain of the DLA simulation. Here, the absolute ceiling altitude limiting the EPB plume height is given by an empirical formula, *h*_{apex}[km] = 497.12 + 6.09 · F10.7[sfu], which is based on phenomenological results reported by Joshi et al. (2022). In addition to raising the *h*_{apex} ceiling, in order to actually make the EPB plumes larger/taller, the number of Brownian random-walking virtual particles (i.e., the number of DLA building blocks *n*_{virtdust}) must also be increased in proportion. Temporal requirements regarding the start time and decay process of EPBs can be satisfied by setting a proper start switch and adjustable fading envelope according to solar local time (SLT). Dynamic requirements regarding the EPB drift motion and tilt angle can be satisfied by adopting an appropriate zonal drift velocity model. Finally, the geomagnetic conjugate mapping requirement can be satisfied by employing a relatively simple geomagnetic field line model. Here, we used a simple quasi-dipole model (see, e.g., VanZandt et al., 1972), with the geomagnetic equatorial line and the declination angles pre-computed using the International Geomagnetic Reference Field model (Thebault et al., 2015; Alken et al., 2021).

Figure 3 shows a diagram of the basic scheme that we adopted for this simplified EPB plume structure model. The scheme consists of two separate components: (1) plume geometry calculations performed using the DLA fractal algorithm and (2) the background ionospheric plasma density model. Figure 3(a) shows a snapshot of the output from the plume geometry calculations, and Figure 3(b) shows a corresponding snapshot of the background ionospheric plasma density model at the same time epoch. The colorbar in Figure 3(a) delineates the relative density depletion (*rdd*) values ranging from −1 to 0, where 0 indicates the unperturbed condition and −1 indicates fully depleted plasma. The plumes are initially assigned a 90% depletion level, i.e., *rdd* = −0.90, which is a moderation of the extreme 99.9% depletion level reported by Tsunoda et al. (1982). Later, the *rdd* values gradually decay toward zero following a prescribed fading envelope. The colorbar in Figure 3(b) represents the plasma density values (in units of m^{−3}). Figure 3(c) shows the final model output, where the simulated plume geometry is incorporated into the ionospheric plasma density distribution. Mathematically, we may express this operation as *N _{e}* =

*N*

_{e0}× (1–|

*rdd*|), where

*N*

_{e0}is the background ionospheric density and

*N*is the resulting ionospheric density. All of these calculations are performed along the vertically oriented geomagnetic equatorial plane. For the background ionospheric density, we used outputs from the international reference ionosphere (IRI) model, onto which the DLA-based EPB plume structures are subsequently stamped. In principle, the background ionospheric density from other models, such as NeQuick, may also be used in this scheme.

_{e}The longitudinal (i.e., east–west) separation between adjacent plumes in the simulation is controlled by making the initial aggregation site (at an altitude of ~150 km) undulated with a sinusoidal shape. This arrangement would induce the DLA aggregates to start growing preferentially at the tips of the undulation. Here, the undulated initiation site is set to have a spatial periodicity of 5° longitude (~550 km) horizontally and an amplitude of 20 km (i.e., 40 km peak-to-trough) vertically, resembling large-scale wave structures (LSWSs) or periodic upwelling in the bottom ionospheric F-region (Tsunoda, 2015 and references therein).

There is some flexibility in how the altitude of the initial aggregation site may be selected, but it generally should be slightly below the bottomside ionospheric F-layer, where EPB plumes typically begin developing. If the altitude of the initial aggregation site is set either too low or too high, we could lose the realism of this DLA-based approach for modeling EPB plumes. If it is set extremely low (e.g., below 60 km), the resulting EPB plumes could develop too many fingers/branches before reaching the topside ionosphere. If it is set too high, (e.g., at ~350 km near the F-peak), the resulting EPB plumes would appear as if they emerge suddenly out of nowhere.

The precise timing for when these DLA-based structured plasma density depletions start to grow and begin to decay was coded as adjustable parameters, set to local dusk and 1 hr after local midnight, respectively. More precisely, the DLA process itself is initiated at *t*_{1} = 18.5 SLT. Meanwhile, the fading envelope controlling the decay of *rdd* as a function of SLT is given by the following:

1

where *t*_{2} = 26.0 SLT (≡ 02:00 SLT) is a time offset parameter and *τ* = 0.3 hr is the decay time constant. In the model, we also allowed the DLA-based EPB plume structures to drift in the eastward direction after they had started to grow. In particular, we adopted an F10.7- and Kp-dependent zonal drift velocity profile that also contains an altitude shear, based on empirical findings published in the scientific literature (Sheehan & Valladares, 2004; Kil et al., 2009; Barros et al., 2018). For a given day of the year, the nighttime zonal drift velocity at an altitude of 350 km as a function of SLT is calculated using the empirical F-region drift model of Sheehan & Valladares (2004), which is based on a set of biquadratic fit coefficients (with the day of the year, Kp, and F10.7 as independent variables). Given the zonal drift velocity at 350-km altitude, the zonal drift velocities at other altitudes are calculated by assuming a shear (i.e., vertical gradient) in the zonal drift velocity. For simplicity of the shear model, we nominally assume that the zonal drift velocity (positive eastward) changes by −6 m/s for every 100-km increase in height, based on the observations of Kil et al. (2009) and Barros et al. (2018). For the same reason, we assume that the zonal drift velocity changes by +6 m/s for every 100-km decrease in height. The minimum height that we consider for this sheared velocity profile is 100 km, and the maximum height depends on the altitude of the tallest bubble being simulated. This shear is only activated when the magnitude of the zonal drift velocity at 350-km altitude exceeds 10 m/s, in order to ensure that the shear diminishes at low speeds. Using these settings in the model, we were able to match many essential physical characteristics of EPBs that have been observed empirically. However, we note again that although the model displays a realistic geometrical appearance, the output of this simplified EPB plume structure model is not guaranteed to satisfy certain mathematical properties (e.g., incompressibility) that would be inherent in self-consistent solutions generated by full physics-based models.

We also note that because the DLA plume geometry and background ionospheric plasma density come from two separate pipelines (cf. Figures 3(a) and 3(b)), the computational mechanism by which the DLA-based EPB plume structures are initiated and grow is independent of the precise configuration of the background ionospheric density. This basic architecture sets the DLA-based approach (at least in its current version) apart from physics-based models.

Figure 4 shows a few alternative perspectives of the simulated EPB plume structures from the simplified model to demonstrate the realism of this model. Figure 4(a) shows a three-dimensional (3-D) view of the simulated EPB plume structures, using two representative orthogonal planes to display the ionospheric plasma density distribution. Meanwhile, Figure 4(b) shows the vertical TEC representation of the simulated EPB plume structures, after the ionospheric plasma density values have been integrated vertically at several altitude layers from 0 to 9900 km. The vertical TEC representation provides a bird’s eye view of the EPB plume configuration. Geomagnetic field mapping based on the quasi-dipole model (e.g., VanZandt et al., 1972) was used to create the full 3-D extension and the vertical TEC representation of the EPB plume structures.

In Figure 4(b), the solid magenta line represents the geomagnetic equator, and the dashed magenta lines represent ±20° magnetic latitude (MLAT) lines. After the height integration, the branching structures have largely disappeared, especially in the area around the geomagnetic equator. However, we can clearly observe branching patterns at the northernmost and southernmost tips of the EPB plumes away from the geomagnetic equator line. From this bird’s eye perspective, we can also recognize similarities between the simulated EPB plume structures and actual EPB structures observed over the Atlantic sector in optical imagery data from the Global-scale Observations of the Limb and Disk (GOLD) satellite mission (Karan et al., 2020; Rodriguez-Zuluaga et al., 2021). These similarities suggest the possibility of cross-comparison between statistical ensembles produced by this simplified EPB plume structure model and empirical ensembles of EPB geometrical morphology obtained via various observation methods.

In the following section, we present a comparative analysis and characterization of the DLA-based simulated EPB plume structures, with respect to observation data and other models.

## 4 COMPARATIVE ANALYSIS AND CHARACTERIZATION

Figure 5 shows a comparative analysis of basic EPB morphology between DLA-based simulation results and a representative set of real-world observations. Figure 5(a) presents the DLA-based simulated PEB plume structures in terms of the IRI model TEC (cf. Figure 4(b)), displayed using geomagnetic coordinate representation. Each column in the color plot depicts a slice of vertical TEC along a magnetic meridian that crosses the geomagnetic equator at the specified geographic longitude. Note also some numerical artifacts between ±5° MLAT in the form of horizontal stripes, which correspond to pixels that were missed by the flat (near-zero inclination) magnetic field lines during the conjugate mapping process north and south from the geomagnetic equator. The EPB plume structures can be recognized as regions of depleted vertical TEC. The SLT for various longitudes at this epoch is indicated just above the horizontal axis. Older bubbles that have already begun decaying are located toward the right-hand side of the plot (greater SLT values), whereas newer bubbles that have just started growing are located toward the left-hand side of the plot (smaller SLT values). Figures 5(b)–(d) display several examples of EPB plume structures captured in measurements of optical airglow emissions from the upper atmosphere, represented in a mix of geographic and geomagnetic coordinate systems. These examples are taken from works reported by Kelley et al. (2003), Kil et al. (2009), and Comberiate & Paxton (2010), respectively.

A few common features can be discerned between the simulation result and real-world observations. The first specific commonality is observed in Figures 5(a) and 5(b), where we find prominent florettes at the tip of the EPB plumes. This feature is generally expected because the tip is the freshest part of the EPB plume, where new small-scale branching structures may still be developing. In contrast, at the root/stem of the EPB plume, some of the small-scale structures may have started to decay/dissipate. The second set of commonality between the simulation result (Figure 5(a)) and observations (Figures 5(b)–(d)) is the inverse-C bending shape of the EPB plume. Such an inverse-C shape along the horizontal (longitude–latitude) plane is associated with the westward tilt of EPB plumes when viewed on the vertical equatorial plane. For real-world EPBs, the westward tilt is caused by the altitude-dependent zonal drift velocity of the EPBs, due to altitude variations in plasma conductivity. In the DLA-based simulation, this westward tilt is induced by prescribing a zonal velocity profile that contains a shear, as described in Section 3. As shown in Figure 5(a), the inverse-C bending is more prominent for older bubbles because of the ample amount of time that has elapsed for the shear to induce full tilt/bending. At approximately 20°W longitude, the simulated EPB plume exhibits a westward bending of ~3° longitude over a ~12° latitude extent from the geomagnetic equator. Meanwhile, the sample observations (cf. Figures 5(b)–(d)) indicate a westward bending of 4°–6° longitude over a ~20° latitude extent from the geomagnetic equator. Overall, the amount of inverse-C bending in the simulation and representative observations shows reasonable agreement.

Next, we examine the rate of EPB plume growth in the DLA-based simulation, with comparison to a few representative cases of physics-based model runs drawn from the literature. This analysis can characterize the typical time for a fully developed EPB plume to form, at which point its apex height has reached the topside ionosphere.

Figure 6 shows sequential snapshots of DLA-based EPB plume simulation results around 0°E longitude from 19:57 Coordinated Universal Time (UTC) (*t* = 19.95 hr) to 21:42 UTC (*t* = 21.70 hr) with 3-min (= 0.05 hr) intervals. The growth progression of two EPB plumes near 0°E longitude was tracked to determine the time evolution of their apex heights. In each frame, we mark the altitude of the plume base at the bottomside F-region, as well as the altitude of the plume top. The difference between the two altitudes (i.e., from base to top of a given plume) is designated the net plume height. An epoch shortly before the first snapshot, determined by a backward linear extrapolation of the first and second snapshots down to zero net plume height, was taken as the initial time (*t*_{0}) for assigning the elapsed time. This initial time was numerically determined to be *t*_{0} = 19.82 hr. We performed the same procedure on a few cases of physics-based model simulations drawn from the literature. Represented in the samples were PBMOD simulation runs from Retterer (2010) and Huang et al. (2012), a HIRB simulation run from Yokoyama et al. (2014), and an equatorial bubble simulation code in two dimensions (PBM2D) simulation run from Carrasco et al. (2020).

In these sampled physics-based simulation runs, the initial amplitude perturbations ranged from ~1% (Retterer, 2010) and ~5% in plasma density (Huang et al., 2012; Carrasco et al., 2020) to ~17% in integrated Pedersen conductivity (Yokoyama et al., 2014). The initial perturbations were sinusoidal, with wavelengths of 200–1200 km. As indicated by the listed choice of initial perturbation parameters, these samples represent a diverse range of simulation conditions. The diverse sampled cases are actually good for comparison with the DLA model, because the DLA model works very differently from physics-based models in terms of how the EPB plumes are initiated, driven, and/or controlled. Thus, the comparison strategy in this case is to estimate a range of EPB plume growths from several known samples that had previously been realized in published modeling works and to subsequently determine whether the EPB plume growth progression from the DLA-based model generally falls within the range of these samples.

Figure 7 shows a composite plot of net plume height versus elapsed time for the DLA-based simulation and the sampled cases of physics-based simulations. In the graph, the two blue curves represent the two tracked plumes from the DLA-based simulation (cf. Figure 6). The physics-based simulations are represented by green curves. One curve each represents HIRB (Yokoyama et al., 2014) and PBM2D (Carrasco et al., 2020) simulation runs. Three curves represent PBMOD simulations: PBMOD case #1 is from Huang et al. (2012), and PBMOD cases #2 and #3 are from Retterer (2010). During the first 30 min, the growth of the DLA-based EPB plumes was comparable to the highest two cases of EPB plume growth from physics-based simulations. Past the first 30 min and approaching the 1-hr mark, the growth of the DLA-based EPB plumes corresponded to the lower quartiles of the samples. Past the 1-hr mark, the EPB plume growth began to taper off for both the DLA- and physics-based simulations, as the net height of the simulation EPB plumes reached at least 300 km relative to the bottomside F-region (i.e., reaching an actual altitude of at least 500 km). Overall, the growth of DLA-based EPB plumes was found to be on par with that of the sampled physics-based simulation runs, taking roughly 1 hr for a plume to become fully developed by reaching topside ionospheric altitudes. In theory, the growth rate of EPB plumes in a DLA-based simulation can be adjusted by manipulating the number of Brownian motion cycles performed in each time step. The ability to achieve precise control over the temporal growth of EPB plumes using DLA will be the subject of future research investigations.

Another set of examination that can be performed is a power spectrum characterization of the simulated EPB plume structures. Here, we are generally limited to structures larger than the outer scale size of individual plumes, as the DLA-based model does not produce structures much smaller than the width of a plume. Nevertheless, such a characterization would still provide useful information. Two approaches were taken: one using spatial profiles of density and vertical TEC at a fixed time and one using the temporal profile of the vertical TEC along a fixed LOS.

Figure 8 shows the setup and results of power spectrum characterization obtained via the first approach, i.e., using the ionospheric plasma density and TEC at a fixed time. Figure 8(a) presents a colormap plot of ionospheric density (along the geomagnetic equatorial plane) at 03:00 UTC as a function of longitude and altitude. A dashed line at 350-km altitude marks the cross-sectional cut line on which we examine the ionospheric density variation. Figures 8(b) and 8(c) show the spatial profile of vertical TEC and ionospheric density (on the specified cut line), respectively, as a function of longitude at this time epoch. Figures 8(d) and 8(e) show the computed power spectra (on a logarithmic scale) of vertical TEC and ionospheric density, as a function of spatial frequency. Fit lines to the power spectra (with confidence bounds) are also shown in the graphs. These fit lines reveal a constant slope in the log–log scale, which demonstrates power-law behavior in both power spectra. The power-law exponents (derived from the slopes) were −1.88 and −1.30 for the vertical TEC and ionospheric density power spectra, respectively. These power-law exponents (or spectral slopes) from the fits are considerably shallower than the spectral slopes of irregularities in the interior of EPB plumes, which are typically −2.50 to −3.00 (e.g., Livingston et al., 1981; Rino et al., 1981, 2016, 2018; van de Kamp & Cannon, 2009). We did not obtain these steeper spectral slopes because we are characterizing only the power spectra of vertical TEC and density structures with scale sizes larger than the outer scale size of an EPB plume. It is known that ionospheric irregularities in the two regimes have different spectral slopes, with irregularities smaller than the outer scale size having a steeper spectral slope (e.g., Nickisch, 2004; Carrano & Rino, 2016). For irregularities larger than the outer scale size, spectral slopes between −1.00 and −1.99 have been experimentally observed (Basu et al., 1983; Kelley & Hysell, 1991; Vijayakumar & Pasricha, 1997; Kil & Heelis, 1998).

Figure 9 shows the setup and results of power spectrum characterization obtained via the second approach, i.e., using the temporal profile of the vertical TEC from a fixed LOS. Figure 9(a) displays a colormap plot of ionospheric density (along the geomagnetic equatorial plane) at 03:00 UTC as a function of longitude and altitude. The SLT values are indicated just above the horizontal axis. The pink line indicates the designated vertical LOS (at 60°W longitude) used for the characterization. Over time, we allowed the EPB plumes to pass through the LOS, and we recorded the vertical TEC as a function of time. Figure 9(b) shows a keogram of vertical TEC values for various longitudes as a function of UTC. In general, a keogram is an assembly of sequential strips along a certain line, which was originally invented and named by Eather et al. (1976). The aforementioned procedure is equivalent to taking a slice of this TEC keogram at 60°W longitude, indicated by the dashed line. Figure 9(c) shows the time series plot of vertical TEC obtained using this procedure, with SLT values for the time axis. The dashed line indicates the baseline (unperturbed) TEC that would be observed if the EPB plume structures were absent, and the purple line segment on the time axis indicates the subset of the vertical TEC time series that was used for computing the power spectrum. The EPB plumes drifted eastward, crossing the designated LOS at a velocity of ~110 m/s, and we used the velocity value to convert from temporal frequency to spatial frequency. Figure 9(d) shows the computed power spectrum (on a logarithmic scale) with a fit line and confidence bounds. Here, a spectral slope of −1.79 was found from the fit, which is comparable in magnitude to the spectral slope in Figure 8(d).

To complement the power spectrum characterization, additional data graphs are presented in Figures 9(e) and 9(f). Figure 9(e) presents a time series plot of the net vertical ΔTEC, obtained by subtracting the baseline from the original vertical TEC time series (cf. Figure 9(c)). This plot highlights the depth of the EPBs that appear in the simulation, in this case reaching almost −10 TEC units (TECU) at the greatest magnitude. Figure 9(f) displays a time series plot of the apparent TEC gradient, calculated using a method equivalent to the time-step method. The EPB drift velocity of ~110 m/s was used to convert *d*TEC/*dt* into ∇TEC in this case. Here, the greatest ∇TEC magnitude is just under 0.2 TECU/km (= 32 mm/km at the GPS L1 frequency). This peak value appears rather small for TEC gradients associated with EPBs, which is partly due to the size of the timegrid for the simulation, which is 0.05 hr = 3 min = 180 s. One may compare this to the time interval for GNSS observation data from a continuously operating receiver station network, which is typically 30 s. The relatively coarse timegrid for the simulation might effectively smear the EPB side wall(s) as perceived along the designated vertical LOS (cf. Figure 9(a)). Thus, the apparent ∇TEC magnitude might be underestimated. The underestimation factor depends on the spatial size of the transition region between the EPB interior and ambient ionosphere. This edge of an EPB is often characterized as a major source of TEC gradients in low-latitude GBAS threat models (Yoon et al., 2017, 2019; Wang et al., 2021). In a case example analyzed by Saito & Yoshihara (2017), the steepest part of the EPB side edge/transition region was estimated to be ~5.3 km in size, which we shall adopt for our present assessment. Given the EPB drift velocity of ~110 m/s in our case, a short time step of Δ*t* = 30 s will give a spatial step size of Δ*x* = 3.3 km. Meanwhile, a longer time step of Δ*t*′ = 3 min will give a spatial step size of Δ*x*′ = 20 km. Because Δ*x*′ = 20 km > 5.3 km, the steepest part of the EPB side edge will be undersampled, leading to an underestimated ∇TEC magnitude. The underestimation factor is 20 km/5.3 km = 3.8 (reducing the spatial step size below 5.3 km will not change the edge steepness any further in the sampling). Adjusted for this underestimation, the EPB-associated peak ∇TEC in this simulation would be corrected to 3.8 × 32 mm/km = 122 mm/km at the GPS L1 frequency.

Although not reaching the |∇TEC| upper bound in existing SBAS and GBAS threat models (Saito et al., 2017; Yoon et al., 2017), this compensated value is generally on par with typical magnitudes of TEC gradients observed experimentally at low-latitude locations (Pradipta & Doherty, 2016; Saito & Yoshihara, 2017; Chang et al., 2021; Affonso et al., 2022). Regardless, the objective here is not to match the absolute upper bound ∇TEC. Rather, the objective is to have an ability to represent instances of anomalous ∇TEC with enough contrast relative to nominal/unperturbed conditions over low-latitude regions. In terms of contrast between EPB-related anomalous ∇TEC and the envelope of nominal ∇TEC, as indicated in Figure 9(f), the anomalous-to-nominal |∇TEC| ratio is approximately 4.1 in this case. In measurement data from Rungraengwajiake et al. (2015) and Budtho et al. (2018), |∇TEC| values of 20–100 mm/km were frequently observed during EPB occurrences, and the typical size of the ∇TEC envelope under nominal condition was ~15 mm/km. These findings imply that anomalous-to-nominal |∇TEC| ratios in the range of 2.0–7.0 can be roughly considered as sufficient contrast for EPB-related anomalous ∇TEC.

## 5 BASIC COMPUTATIONAL PERFORMANCE ASSESSMENT

In this section, we present a basic runtime characterization of various computational modules. Figure 10 illustrates the sequence of computational modules that are involved in a batch run. The simulation domain is over a longitude span of 130°, geographic latitude range of −40° to +40°, and altitude range of 0 to 1100 km. The spatial grid spacings are `dlon` = 0.1°, `dlat` = 0.1°, and `dh` = 10 km. The smooth background ionospheric density from the IRI model is computed on a coarser grid with `dlon`’ = 7.5° and `dlat`’ = 7.5°, which is subsequently interpolated onto the finer grid. The time step between consecutive epochs is `dt` = 0.05 hr (= 3 min). A more detailed list of the initial simulation parameters and their basic settings is given in Table 2. For a simulation covering a 24-hr period, there is a total of 481 data epochs or time frames. Test simulations were performed using MATLAB version 7 on a Windows 7 laptop computer (manufacture year: 2013) with Intel Core i5-3337U 1.80 GHz CPU and 4.00 GB RAM. The recorded model runtimes for various modules are summarized in Table 3, and additional details are provided below.

Parameter initialization takes 0.001 s, which is negligible. Plume geometry construction using DLA takes 683 s = 11.4 min. The rendering of various 2-D and 3-D graphics of *rdd* values takes 14092 s = 235 min = 3.9 hr. The IRI model run to generate the background electron density configuration along the vertical geomagnetic equatorial plane takes 373 s = 6.2 min. The IRI model run along a horizontal hmF2 surface takes 946 s = 15.8 min. The IRI model run along multiple flat horizontal planes from 150-km to 750-km altitude (100-km spacing) takes 6408 s = 107 min = 1.8 hr. Graphics rendering of the resultant electron density configuration along the vertical geomagnetic equatorial plane takes 1461 s = 24.4 min. Rendering of 3-D graphics of the resultant electron density configuration (along both the vertical geomagnetic equatorial plane and horizontal hmF2 plane) takes 8168 s = 136 min = 2.3 hr. Graphics rendering of the resultant TEC configuration takes 5249 s = 87.5 min = 1.5 hr.

For a given calendar date of the simulation, one only needs to perform the IRI model runs once, which, in this case, would take a fixed combined runtime of 6.2 min + 15.8 min + 1.8 hr = 2.2 hr. Some of the graphics rendering modules (especially those for the *rdd* values) are not essential and can be skipped, which saves 3.9 hr of computational time in each simulation run. The module that would be most beneficial to repeat multiple times is the plume geometry construction, which has a runtime of slightly over 10 min for each simulation. Hence, in principle, it would be possible to generate a statistical ensemble of several hundreds of separate scenarios of EPB growth and evolution after only a few days of computational work.

It should be noted here that the runtimes of some modules are rather long because of inefficiency in the current implementation of the program. In particular, the IRI model computation had an unreasonably long runtime in the tests because of a mixed implementation involving repeated system calls from MATLAB to the FORTRAN executables for generating individual profiles. In the future, the implementation can be optimized by letting the bulk of the IRI model computation run via the FORTRAN code alone, without multiple unnecessary system calls from MATLAB. This may reduce runtimes in the future.

## 6 POTENTIAL APPLICATIONS

When we consider the steady progress made in the history of EPB modeling efforts, the main qualities most often demanded from numerical EPB models include the capability of “*reproducing EPBs that grow nonlinearly from the crest of a large-scale upwelling, bifurcate at the top of EPBs, then become turbulent at the topside of the F-region*” (Yokoyama, 2017). In general terms, the DLA-based simulated EPB plume structures possess these essential qualities. Of course, there is also the higher goal of “*a full high-resolution global simulation [that] could model all phenomena at all scales*” (Yokoyama, 2017), which would be ideal for ionospheric scintillation simulation purposes. Strictly speaking, the DLA-based EPB modeling would not achieve this aspiration directly. Nevertheless, with some modification and augmentation, this approach can still serve as a useful platform for simulating low-latitude scintillations, including the nonstationary nature of scintillation signals.

Conventionally, simulations of scintillation due to ionospheric density irregularities are performed by using the phase-screen model, which replaces a slab ionosphere of finite thickness *L* with a single thin diffracting phase screen (see, e.g., Rino, 1979a; Cervera & Thomas, 2006; Vasylyev et al., 2022). In a phase-screen model, ionospheric irregularities are represented by an irregularity strength parameter *C _{s}* (or

*C*≡ 1000

_{k}^{2v+1}

*C*) defined as follows:

_{s}2

with *v* = *p*/2, where *p* is the phase spectral slope, *q*_{0} is the outer scale wavenumber, Γ is the gamma function, and is the variance of the plasma density fluctuations. It is implicitly assumed that the irregularities uniformly fill the entire ionospheric slab. The difference between *C _{s}* and

*C*is that

_{k}*C*represents the 3-D spectrum of density irregularities at a spatial wavenumber of 1 rad m

_{s}^{−1}, whereas

*C*represents the strength of the irregularity spectrum at a spatial wavelength of 1 km. The weak-scatter amplitude scintillation index may then be computed using the following expression:

_{k}3

where *r _{e}* is the classical electron radius,

*λ*is the wavelength of the radio waves,

*θ*is the angle relative to normal incidence to the phase screen,

*F*is the Fresnel filter factor, and

*Z*is the Fresnel zone parameter. The Fresnel zone parameter is given by

*Z*=

*λz*sec

_{r}*θ*/4

*π*, where is the reduced height of the phase screen (

*z*

_{1}and

*z*

_{2}are the distances from the screen to the radio wave source and to the observer, respectively). The amplitude scintillation index

*S*

_{4}can finally be estimated by letting

*S*

_{4w}undergo saturation via .

For more general cases in which we do not have a uniform slab ionosphere, we will at least need to replace *C _{s}L* with an integral ∫

*C*in the formulation (likewise for

_{s}dz*C*), as we no longer have a constant-value

_{k}L*C*(or

_{s}*C*) that fills the entire ionospheric layer. In such situations, the EPB plume model(s) may be invoked. For example, one could use the DLA-based simulated EPB plume structures to set up a temporally and spatially varying irregularity strength parameter

_{k}*C*that is nonzero only within the interior of the plumes. This approach has been previously utilized in the context of scintillation due to ionospheric irregularities inside EPBs along radio occultation ray paths, but with the use of a simple Gaussian-shaped spatial envelope representing the EPB plume (Carrano et al., 2011; Ludwig-Barbossa et al., 2023). In such cases, one would have “direction-dependent”

_{k}*C*values that are highly sensitive to the precise geometrical orientation of the LOS to GNSS satellites and how they intersect the EPB plumes, thus emulating major factors governing the nonstationarity and directional dependence of low-latitude scintillations (Anderson & Strauss, 2005; Abadi et al., 2014; Rino et al., 2023). The formulation and construction of such a spatiotemporally varying

_{k}L*C*suitable for this purpose will be the subject of future work.

_{k}Another potential application of the simplified EPB plume model is to simulate the presence of anomalous TEC gradients over a certain area at low latitudes, which pose considerable threats to the availability and integrity of SBAS/GBAS services in civil aviation. We alluded to this particular application at the end of Section 4, with an estimated TEC gradient magnitude that can be expected in the simulation runs. Although amounting to only a fraction of the |∇TEC| upper bound in existing SBAS/GBAS threat models (Saito et al., 2017; Yoon et al., 2017), the estimated TEC gradient in the model example has a clear contrast relative to nominal/unperturbed conditions and thus may represent instances in which anomalous ionospheric gradients are present over or near SBAS/GBAS facilities. The simplified model could be used to construct various simulated scenarios in which EPB plume structures (containing adverse gradients) grow/develop and drift into place over an SBAS/GBAS service area. Assessing how SBAS/GBAS service is affected in specific simulation scenarios (or in statistical ensembles of multiple simulation scenarios) may be useful in the design and planning of effective early warning systems or other mitigation methods against ionospheric threats to SBAS/GBAS facilities at low latitudes. Here, the basic philosophy is aligned with that put forward previously by Saito et al. (2009a, 2009b), Saito & Fujii (2010), and Yoon et al. (2019).

## 7 FURTHER IMPROVEMENTS AND OTHER POTENTIAL ROLES

As a case demonstration of the simplified EPB model, we used the IRI model as the background ionospheric density distribution in this paper. However, there is a large degree of flexibility in the choice for background ionospheric models, such as the NeQuick model (Nava et al., 2008; European Commission, 2016), if another empirical ionospheric model is desired. Alternatively, a background ionospheric density distribution obtained from a physics-based model, such as a thermosphere-ionosphere-electrodynamics general circulation model (TIEGCM) (Richmond & Maute, 2014; Qian, 2014), can also be utilized with this simplified EPB model. Such a “hybrid” pairing option predicates that the DLA-based simplified EPB plume model described in this paper does not have to be completely detached from other physics-based models. In fact, additional information available from other physics-based models can be utilized to make further improvements to the DLA-based simplified EPB plume model. For example, the Rayleigh–Taylor instability growth rate calculated via the TIEGCM model (Carter et al., 2014) can be used to adjust the number of test particles in the DLA simulations, thus controlling the success of the EPB plumes in their growth and development. Such options would pave the way for making the simplified EPB model more versatile in the future.

Putting the issue more broadly, we will be facing a larger question regarding what would be considered comparable inputs between physics-based and DLA-based models. Answering this question means determining a precise relationship between initial conditions in physics-based models and suitable/equivalent parameter settings in a DLA-based model, which essentially constitutes a “parameter interfacing” problem. Having a well-calibrated “input parameter interface” between these two classes of EPB plume structure models would be ideal, but perhaps can only be successfully accomplished over the long term. Such characterization would be an important part of future research investigations.

With the possibility of a simplified-yet-realistic 3-D EPB plume structure model, there may also be opportunities for linkage with other technology categories. Potential opportunities could emerge in the area of graphics visualization technology, such as virtual reality, augmented reality, and 3-D volumetric displays. With current trends in the commercialization of space activities, there may be growing interest in displaying various space phenomena and the orbit of space assets using advanced graphics visualization techniques. A comprehensive 3-D visualization of risks and hazards posed by EPB/scintillation using these novel technologies may help to achieve some of the main goals of space situational awareness or space domain awareness for both public and private aerospace sectors. In addition, these tools may also be useful in the context of educational settings.

If we consider a more formal institutional framework, the simplified EPB plume model may have a potential role as a component in observing system simulation experiments (OSSEs) for low-latitude ionospheric space environment monitoring. An OSSE is a modeling experiment often used to evaluate the value of a new observing system (e.g., Arnold & Dey, 1986; Hoffman & Atlas, 2016; Zeng et al., 2020). Although the context may differ, the basic motivation behind OSSEs is related to that of crisis simulation exercises used in the global finance/banking sector (Lee, 2023) or space weather reasonable worst-case scenario exercises (Hapgood et al., 2021) and realistic space weather scenario response and recovery exercises (National Science and Technology Council, 2015) in space weather impact preparedness efforts. However, a significant amount of refinement and optimization in the DLA-based EPB plume model would be needed before adoption into a potential OSSE role.

## 8 SUMMARY AND CONCLUSION

We have discussed the theoretical possibility of formulating a simplified EPB plume structure model for practical applications. In foundational terms, we presented a justification of why such a simplified model would not be a backward step among the current availability of physics-based models. In specific/technical terms, we outlined a concrete formulation of such a simplified EPB plume structure model based on the DLA fractal algorithm. This simplified model operates using two separate main modules: one module for generating the bifurcating patterns of EPB plumes using the DLA algorithm and one module for generating the background ionospheric plasma density distribution. The outputs from these two modules are subsequently combined to generate the final model output. In addition to these two main modules, several other additional settings were put forth to make the simulated EPB plume structures resemble as many known physical characteristics of EPBs as possible. The results of initial tests qualitatively indicate a reasonably high degree of realism when compared with actual EPB plume structures observed in empirical data. However, we note that specific structures generated in the model simulation simply represent possible realizations, which necessitates statistical ensemble runs of the model for us to draw a deeper meaning.

Ideally, the usage of this simplified model of EPB plume structures should be combined with adequate reference to a long-term EPB and scintillation occurrence climatology (e.g., Pradipta et al., 2021) in order to ensure a proper context for this low-latitude ionospheric phenomenon. Combined with long-term climatology, this simplified EPB model can help to establish a better collective awareness of the hazards posed by EPBs/scintillations over low-latitude regions. Finally, it is hoped that the implementation of EPB/scintillation modeling as a regional space weather hazard assessment tool would become ubiquitous, enabling a much lower dependence on high-performance computing facilities.

## HOW TO CITE THIS ARTICLE

Pradipta, R., Carrano, C. S., Groves, K. M., & Doherty, P. H. (2024). Can numerical simulations of equatorial plasma bubble plume structures be simplified for operational and practical usage? *NAVIGATION, 71*(2). https://doi.org/10.33012/navi.645

## CONFLICT OF INTEREST

The authors declare no potential conflicts of interest.

## ACKNOWLEDGMENTS

This work was supported by the Federal Aviation Administration through Memorandum of Agreement DTFAWA-17-X-80005 with Boston College.

## Footnotes

↵2 We dedicate this paper to the memory of Patricia Doherty (d. 14 July 2022), who continues to inspire us.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.