Quantifying Uncertainties in the Quiet‐Time Ionosphere‐Thermosphere Using WAM‐IPE

This study presents a data‐driven approach to quantify uncertainties in the ionosphere‐thermosphere (IT) system due to varying solar wind parameters (drivers) during quiet conditions (Kp < 4) and fixed solar radiation and lower atmospheric conditions representative of 16 March 2013. Ensemble simulations of the coupled Whole Atmosphere Model with Ionosphere Plasmasphere Electrodynamics (WAM‐IPE) driven by synthetic solar wind drivers generated through a multi‐channel variational autoencoder (MCVAE) model are obtained. Applying the polynomial chaos expansion (PCE) technique, it is possible to estimate the means and variances of the QoIs as well as the sensitivities of the QoIs with regard to the drivers. Our results highlight unique features of the IT system's uncertainty: (a) the uncertainty of the IT system is larger during nighttime; (b) the spatial distributions of the uncertainty for electron density and zonal drift at fixed local times present 4 peaks in the evening sector, which are associated with the low‐density regions of longitude structure of electron density; (c) the uncertainty of the equatorial electron density is highly correlated with the uncertainty of the zonal drift, especially in the evening sector, while it is weakly correlated with the vertical drift. A variance‐based global sensitivity analysis suggests that the IMF Bz plays a dominant role in the uncertainty of electron density. A further discussion shows that the uncertainty of the IT system is determined by the magnitudes and universal time variations of solar wind drivers. Its temporal and spatial distribution can be modulated by the average state of the IT system.

ZHAN ET AL.
The sources of the day-to-day variability of the IT system mainly come from the varying solar radiation, geomagnetic activity, and meteorological forces.These resources have different impacts on the IT system at different locations and local times.Previous studies (Fang et al., 2018;Materassi et al., 2020) show that geomagnetic activity is the main contributor to the NmF2 variability and TEC variability globally, and the contributions from solar radiation, solar wind, and meteorological force at low latitudes are distinct from those at mid-and high-latitude regions.The day-to-day variability of maximum electron density (N max ) represented by its standard deviation shows the smallest contribution from solar irradiance and equal contributions by the geomagnetic activity and meteorological forces from below (Materassi et al., 2020).Materassi et al. (2020) also shows a different variability across the E-and F-regions.
Previous studies have shown pronounced day-to-day variability of plasma drifts, electron density, neutral winds, and Rayleigh-Taylor instability growth rates in the equatorial F region due to tidal forcing from the lower atmosphere (meteorological forcing), solar radiation, solar wind, and geomagnetic storms (Wang et al., 2021).In previous observational studies (B.Fejer & Scherliess, 2001;B. G. Fejer et al., 2005), authors have usually used K p = 3 or 4 as a threshold to divide the data set into quiet and active groups and attribute the characteristics of quiettime day-to-day variability of PRE or ESF to the tidal forcing propagating from the lower atmosphere.However, during this relatively quiet condition, the solar wind impacts may still exist and, therefore, need to be quantified.While it is difficult to specify the uncertainties from all these sources in one study, we will specifically investigate the impact of the uncertainty of solar wind on the IT system during quiet to moderately disturbed conditions (the maximum K p during the current and previous days is smaller than 4).This study quantifies the degree to which relatively quiet solar wind conditions can impact the equatorial and low-latitude IT system.
The solar wind parameters directly interact with the high-latitude IT system and have an impact on the low-latitude and equatorial ionosphere mainly through large-scale traveling ionospheric disturbance (LSTID) induced by Joule heating, prompt penetration electric field (PPEF), and disturbance dynamo electric field (DDEF) (i.e., Bagiya et al., 2011, and references therein).These mechanisms have different local time and longitudinal dependencies (Huang et al., 2010;Pavlov & Pavlova, 2007;Pavlov et al., 2008), and which one will dominate the variation of the equatorial ionosphere is still an open question (B.G. Fejer, 2011).
Uncertainty quantification (UQ) is an emerging field that aims to compute the amount of uncertainty associated with particular quantities of interest (QoIs) in a physical system (i.e., Licata et al., 2022;Jivani et al., 2023).One technique widely used in UQ is to approximate the QoI in a series of orthogonal polynomials known as polynomial chaos expansion (PCE) (Doostan & Owhadi, 2011;Ghanem & Spanos, 2002;Xiu, 2010) and estimate the QoI low-order statistics or the sensitivity of QoI to the input through the expansion coefficients (Crestaux et al., 2009).PCE has been applied in a recent space weather study by Jivani et al. (2023) to identify the important model parameters that affect the solar wind prediction with the Alfven Wave Solar atmosphere Model (AWSoM).
The remainder of this paper is organized as follows: In Section 2, we will introduce the WAM-IPE model, the selection of QoIs and uncertain input parameters, the PCE-based UQ and sensitivity analysis method, and the numerical experiments we conduct.In Section 3, the results from simulations with WAM-IPE focus on the temporal and spatial variation of the uncertainties of QoI, and the sensitivity analysis results regarding the uncertain input parameters are presented.We demonstrate that even during quiet conditions, the solar wind parameters can vary within a relatively large range and affect the equatorial and low-latitude IT system.In Section 4, we will discuss the implications of the temporal and spatial distribution of the uncertainties in the IT system and the sensitivity analysis results.The main findings and conclusion are included in Section 5.

Models and Methods
Figure 1 shows the diagram of the whole workflow in this study.We will quantify the uncertainty of WAM-IPE outputs (electron density N e , plasma flow V i , and neutral wind U n ) corresponding to the varying solar wind drivers (interplanetary magnetic field Bz, solar wind speed, and solar wind density).UQ and sensitivity analyses of these QoIs using a data-driven PCE method are conducted.Instead of using historical solar wind measurements as drivers, we will train and generate random samples via a multi-channel variational auto-encoder (MCVAE), which captures the statistical features from historical measurements of multiple solar wind parameters simultaneously ZHAN ET AL.

10.1029/2023SW003665
3 of 21 through the encoder and generates synthetic realizations compatible with the data through the decoder.The MCVAE method is an example to show how we deal with large dimensions of input parameters.PCE-based uncertainty propagation method has been used in other areas but is very new in Space Physics/Space Weather community.This study is also a first step of our project which intends to investigate how the generation of equatorial plasma bubbles varies with different input parameters which include solar wind, forces from the lower atmosphere, and so on.

WAM-IPE and QoIs
WAM-IPE is a physics-based model developed by researchers at the University of Colorado Boulder and the National Oceanic and Atmospheric Administration Space Weather Prediction Center (NOAA SWPC) to provide operational Ionosphere-Thermosphere-Mesosphere (ITM) forecasts.The current version of WAM was built on the Global Spectral Model (GSM) of the NOAA Global Forecasting System (GFS) weather model (GSM/WAM for short).The current spatial resolution is T62 (roughly 1.8° × 1.8° in latitude-longitude) with 150 vertical levels spanning from sea level through the upper thermosphere to a pressure level of 3 × 10 −7 Pa (between 400 and 600 km height, depending on solar activity).The IPE part uses a minimum of 13,600 flux tubes, with 170 in latitude and 80 in longitude directions, resulting in a longitudinal resolution of 4.5°.The number of flux tubes is adjustable, and a non-uniform grid is possible.
The current version of WAM-IPE includes physical and chemical processes relevant to the upper atmosphere (Akmaev et al., 2008), as well as electrodynamic and plasma processes simulated by the IPE component (Maruyama et al., 2016;Sun et al., 2015).IPE is a time-dependent, 3-D model of the global ionosphere and plasmasphere in magnetic flux tube coordinates.Within a given IPE flux tube, the solver based on the Field Line Interhemispheric Plasma (FLIP) model (Richards, Fennelly, & Torr, 1994) simulates the plasma density, composition, and temperature.An electrodynamics solver (Richmond, 1992) calculates the electric field and feeds it back to the plasma transport algorithm.The high-latitude potential is specified by the Weimer empirical model (Weimer, 2005), and the model uses an arithmetic mean of the daily F10.7 proxy and the 81-day average of the F10.7 proxy to feed a combination of the EUVAC/HEUVAC solar irradiance models (Richards, Torr, et al., 1994).
WAM fields, including winds, temperature, and molecular and atomic atmospheric composition, are fed into IPE to enable the plasma to respond to changes driven by the neutral atmosphere.The coupling is based on time-dependent 3D re-gridding carried out by the Earth System Modeling Framework (ESMF).The model ingests solar wind and geomagnetic inputs, provided both by direct observation and by the Space Weather Forecast Office (SWFO) forecast.The WAM-IPE Forecast System (WFS) that runs in operation at SWPC currently provides 2-day forecast products utilizing the forecast solar wind and geomagnetic indices provided by SWFO at SWPC.
In this study, we drive WAM-IPE with the MCVAE-generated solar wind data and define QoIs as the ionospheric electron density, plasma drift, and neutral wind to evaluate their uncertainties.

Drivers and Their Data-Driven Uncertainty Representation
To create the solar wind drivers using MCVAE, 5-min time resolution solar wind measurements during 1981-2021 were obtained from NASA OMNI-Web (https://omniweb.gsfc.nasa.gov/).The parameters, including solar wind proton density, speed, and IMF Bz, that are used to drive WAM-IPE are selected.We only use the solar wind data under geomagnetic quiet conditions, which are defined by the days when the maximum K p is smaller than 4 on the current day and the previous day.Otherwise, this day is defined as a disturbed period.This criterion is used because the geomagnetic storm effects may last for 1-2 days (Scherliess & Fejer, 1997).
As each solar wind parameter has 288 values per day, it can lead to significant complexity for the UQ procedure.Thus, it is important to first reduce the dimensions of these parameters.It is also critical to model the joint probability density function of the data from their historical measurements in order to generate new samples from the joint density for our UQ and sensitivity analysis.In the present study, the VAE technique, an approach that meets these two requirements, was chosen.In particular, we use the multi-channel version of VAE (MCVAE) (Antelmi et al., 2019b) as it allows us to preserve the statistical dependence among the drivers.
For simplicity, the encoder of a three-layer neural network in one channel VAE will map the solar wind measurements to the latent variables   ∈ ℝ  , and the decoder consisting of another 3-layer neural network will reconstruct the measurements from the latent space.The latent variables z are formed by the linear transformation of a standard Gaussian vector ξ via learnable mean μ and standard deviation σ vectors, where ⊙ denotes entry-wise product and I is the identity matrix of size d.The goal of the VAE training is to obtain the best reconstruction of the data while ensuring   (, ) distribution on z.The setup of MCVAE is adjusted from that used in Antelmi et al. (2019b) and more details can be found in the paper.A brief description can be found in Appendix A.
After training and testing the MCVAE model, we find that a latent space of size 30 leads to the smallest validation error in the MCVAE training.As such, the Gaussian random vector ξ in Equation 1 describing the uncertainty in the drivers (via the decoder) is of dimension d = 30 and is used to generate 500 synthetic solar wind realizations.As we shall explain in Section 2.3, ξ plays a key role in the construction of our uncertainty propagation scheme.Figure 2 shows the 500 samples (gray lines) of IMF Bz, solar wind density, and solar wind speed used to drive WAM-IPE.The corresponding means and standard deviations are shown in red curves and blue error bars, respectively.The majority of the IMF Bz samples vary between −10 and 10 nT with a mean of around 0, and the average shows slight universal time variation.Likewise, the majority of the solar wind density samples vary below 10 m −3 , with a mean of around 5 m −3 .The majority of the solar wind speed samples vary between 300 and 450 km/s, with a mean of around 380 km/s.These values correspond to typical geomagnetic quiet conditions.

Non-Intrusive Polynomial Chaos Expansion
PCE provides a framework to approximate the solution of a stochastic system by projecting it onto a basis of orthonormal polynomials of random inputs.We summarize the key aspects of PCE below and refer the interested reader to Ghanem and Spanos (2002), Xiu (2010), and Doostan and Owhadi (2011   () is the orthogonal polynomial of degree i k in the kth variable.The choice of    () depends on the distribution of ξ k and follows the so-called Askey scheme (Xiu & Karniadakis, 2003).In the present study, the input drivers are described by the non-linear transformation of the VAE's latent variables   ∈ ℝ  (with d = 30) via the decoder, as discussed in Section 2.2.The variables z are in turn linear transformations of the standard Gaussian variables ξ in Equation 1.Therefore, the random variables ξ k , k = 1, …, 30, in constructing the 30-dimensional PCE are independent, standard Gaussian variables, and the basis functions    () are of Hermite type.
In practice, the expansion in Equation 2 shall be truncated for computational purposes.Considering all d-dimensional Hermite polynomials of total degree not exceeding p, u(ξ) can be approximated by where In this case the size of   , hence the number of terms in Equation 3, is given by P = (p + d)!/(p!d!).For the interest of notation, and up to a choice of arrangement, Equation 3 can be rewritten as (4) The main task in PCE is to find the coefficients α i for which several approaches are available; see, for example, (Doostan & Owhadi, 2011;Ghanem & Spanos, 2002;Hampton & Doostan, 2015;Xiu, 2010).In the present study, we employ the compressive sampling approach via ℓ 1 -minimization (Doostan & Owhadi, 2011;Hampton & Doostan, 2015), which requires a set of corresponding samples of inputs ξ and the solution of interest u(ξ), denoted by , obtained via Monte Carlo sampling, for instance.Evaluating Equation 4 at these samples leads to the linear system α α α ≈ , (5) where α = (α 1 , …, α P ) is the unknown coefficient vector, u = (u(ξ (1) ), …, u(ξ (N) )) contains the samples of the quantity of interest at each time or location, and Ψ(i, j)≔ψ j (ξ (i) ), for i = 1, …, N and j = 1, …, P, is the so-called measurement matrix.To solve the system of Equation 5for α we utilize the SPGL1 solver (van den Berg & Friedlander, 2008) via the MATLAB implementation in van den Berg and Friedlander (2019).Following Doostan and Owhadi (2011) and Hampton and Doostan (2015), an accurate solution to Equation 5 requires a number of solution realizations N that linearly depends on the number of the dominant coefficients, possibly much smaller than P.
The accuracy of the PCE model Equation 4depends on the choice of p, N, and the smoothness of u(ξ) with respect to ξ.Assuming higher order derivatives of u(ξ) with respect to ξ are bounded and large enough samples of u(ξ) are available, increasing p results in more accurate PCE approximations of u(ξ).However, in practice, a certain number of samples N can be afforded, in which case increasing p may lead to overfitting.To avoid this and to set a suitable tolerance on ‖Ψα − u‖ 2 needed by SPGL1, we utilize the cross validation algorithm used by Doostan andOwhadi (2011), andPeng et al. (2014) with 75% of the samples of u as a training set and 25% as a validation set.To find the best order of the PCE, we test p = 1, 2, 3 and find that p = 2 leads to a smaller validation error.
As described in more detail in Section 2.4, N = 500 samples of electron density, plasma drifts, and neutral winds are used to generate the PCEs of this study.
ZHAN ET AL. 10.1029/2023SW003665 6 of 21 Once computed, the expansion coefficients α can be used to estimate the statistics of u, construct a surrogate model that maps ξ to u(ξ), or perform sensitivity analysis, as discussed in the following section.For example, we obtain the mean and variance of u, respectively, by

Sensitivity Analysis
To identify the drivers ξ k that contribute the most to the variability of the IT system, we adopt a variance-based global sensitivity analysis method to decompose the QoI variance into parts corresponding to each input variable as well as their combinations.The Sobol' index (Sobol', 2003) is a widely used variance-base method that can be easily obtained through the PCE expansion coefficients.Based on the PCE results, the expressions for the first-order and total effect Sobol' indices based on PCE (Sudret, 2008) are, respectively, given by where , and where is the total variance of the QoI and obtained as described in Section 2.3.1.The first order Sobol' index quantifies the expected reduction in the variance of the QoI when we fix the input ξ k to its mean, and the total effect index describes the contributions from the PCE terms involving ξ k in which ξ k appears either individually or in combination with other inputs.Larger Sobol' indices imply a larger contribution from an input.For this study, we select the total effect Sobol' index     to present the total contribution from each input parameter to the variance of QoIs and sensitivity analysis.As each solar wind parameter corresponds to 10 values in the latent space, we will sum up the corresponding 10 values of Sobol' indices at each time and location to represent the contribution from the uncertainty of each solar wind parameter to the uncertainty of the QoIs.

Design of Experiments
To investigate the uncertainties associated with the solar wind drivers, we fix the F10.7 (120 solar flux units, sfu), average Kp index Equation 2, and other parameters in the model.As mentioned above, N = 500 samples of IMF Bz, solar wind speed, and solar wind density are used to drive the WAM-IPE and obtain realizations of QoIs (electron density, vertical and zonal plasma drifts, and meridional and zonal neutral winds in this study).The solar wind drivers are ingested into WAM-IPE through the high latitude electric potential model Weimer.
All the simulations are made on the same day (16 March 2013) to ensure the impact from the lower atmosphere perturbations in WAM-IPE stays constant and all the variability mainly comes from changes in the solar wind parameters.While we only focus on the quiet time variability of the IT system during varying solar wind conditions, the approach and methods utilized in this study can be applied to other problems (i.e., different drivers, geophysical conditions, and QoIs) depending on the objectives we aim at.

Results
We present the UQ and sensitivity results for electron density, plasma drifts, and neutral winds in this section.The universal time (UT) and local time (LT) evolutions of the global distribution of uncertainties in QoIs are discussed in detail, with a particular focus given to the equatorial and low-latitude regions.In general, the uncertainty in the equatorial and low-latitude regions (the white line-bounded region) is larger during the nighttime.The uncertainty at 0500 UT has large values in a limited region compared to the other UTs.At 1100 UT, larger uncertainty appears in the evening sector between 2000 and 0030 LT at much broader regions at all 3 altitudes.We also see that the regions with large uncertainty at 250 km are closer to the magnetic 20°N and 20°S and emerge at the magnetic equator at 350 km.This feature can be associated with the equatorial ionization anomaly (EIA).This will be discussed in detail in Section 4.1.In the post-midnight, at 250 km, we see larger uncertainty along the magnetic 20°N and 20°S during 0000-0600 LT and smaller uncertainty during 0400-0600 LT at the equator.At higher altitudes, the uncertainty along the magnetic 20°N and 20°S becomes much smaller, while the uncertainty at the equator region between 0400 and 0600 LT still exists but with slightly smaller magnitudes.

Universal Time Evolution
At 1700 UT, we see evening sector features at all three altitudes similar to those occurring at 1100 UT but in a much narrower latitude region at the three altitudes.We see smaller uncertainty along the magnetic 20°N and 20°S at the morning terminator between 0400 and 0800 LT at 300 and 350 km.We also see great uncertainty in the northern hemisphere between 0400 and 0600 LT at 250 and 300 km.At 2300 UT, similar to 1700 UT, we see great uncertainty in the evening sector in the equatorial and low-latitude regions and smaller uncertainty during 0600-0800 LT along magnetic 20°N and 20°S.
From the description of the global distribution of uncertainty in the electron density, the dominant locations for large uncertainties occur in the evening sector and in the low-latitude and equatorial regions, with some differences with altitude.The altitudinal dependence of electron density has been reported in Chen et al. (2016) with maximum electron densities located away from the magnetic equator at lower altitudes and at the magnetic equator at higher altitudes.The larger uncertainty in the evening sector and the altitudinal difference could be associated with EIA.Therefore, the universal time evolution of uncertainty is associated with not only the variation of the solar wind driver (0500 UT vs. 1100 UT) but also the background electron density.

Local Time Evolution
We have shown that the uncertainty of electron density is larger at night.Thus, we further focus on the local time evolution of the uncertainty at night.In Figure 4, the global distribution of electron density uncertainties at 6 different local times (2000, 2200, 0000, 0200, 0400, and 0600 LT, fixed local time at all the longitudes) at 3 altitudes (250, 300, and 350 km) is presented.
In the postmidnight sector (0000-0600 LT), larger uncertainty appears around the magnetic 20°N and 20°S and at mid-to high latitudes at 250 km.The region and magnitude of large uncertainties become smaller at higher altitudes.At 2000 and 2200 LT, the distribution of large uncertainty presents a longitude structure with four peaks along the magnetic equator.The peaks are consistent at specific longitude sectors (30-90°E, 120-180°E, 240-300°E, and 330-30°E) at the 3 altitudes, with a shrinking trend when the altitude increases.This longitude structure could be associated with the longitude structure of the EIA (Chen et al., 2016).The peaks also present a slight shift in longitude from 2000 LT to 2200 LT.

Universal Time Evolution
We see larger electron density uncertainties in the low-latitude and equatorial regions, so we will only look at the uncertainties of plasma drifts and neutral winds in this and the next sections.Plasma drifts and neutrals at   (Scherliess & Fejer, 1997).We present the results for the uncertainty of vertical and zonal plasma drifts at 300 km and at 0500, 1100, 1700, and 2300 UT in Figure 5.
At 0500 UT, the uncertainties of zonal and vertical plasma drift in the low-latitude region are mostly small, with slightly larger uncertainties of zonal drift appearing during 2200-2400 LT.At 1100 UT, vertical drift shows large uncertainty between 2000 and 2200 LT and slightly smaller uncertainty during postmidnight from 0000-0600 LT to 0800-0900 LT at all latitudes.Zonal drift shows large uncertainty during 2000-0300 LT at all altitudes and smaller uncertainty around 0800 LT around the magnetic equator.At 1700 UT, vertical drift also shows large uncertainty during 2100-2400 LT and during postmidnight, mainly between 0300 and 0900 LT, with large magnitudes around 0600 LT.Zonal drift shows a similar distribution of uncertainty in the evening and
From the description above, we see a close correlation between the spatial and temporal distribution of the uncertainties of zonal drift and electron density, while the uncertainty of vertical drift tends to appear independently from that of zonal drift and electron density.

Universal Time Evolution
The universal time evolution of the uncertainty of neutral winds at vertical pressure level 140 (∼300 km) is presented in Figure 7 with the same format as Figure 5.For zonal wind, larger uncertainties occur at 1100 and 1700 UT in the nighttime sector (2000-0200 LT and 0400-0800 LT).At 2300 UT, zonal wind also shows large uncertainty during 0400-0800 LT.For meridional wind, the uncertainty shows a very different feature from that of zonal wind.At 0500 LT, larger uncertainty occurs in the morning sector (0200-1200 LT in the northern hemisphere and 0400-0800 LT in the southern hemisphere).At 1100 UT, large uncertainty occurs in the evening sector in the southern hemisphere (with a moving trend from high to low latitudes and across the geographic equator) and in the morning sector in the northern hemisphere.The uncertainty during the daytime also happens in the afternoon sector, with a moving trend from high to low latitudes.At 1700 UT, larger uncertainty appears in the evening sector, with a trend moving from the northern hemisphere high latitude to low latitude and across the equator to the southern hemisphere in the local postmidnight sector.At 2300 UT, larger uncertainty mainly appears in the postmidnight sector away from the equator in the two hemispheres.These local time evolution features of neutral wind uncertainty indicate the effects of the universal time variation of the solar wind drivers and the differences in the zonal and meridional directions.

Local Time Evolution
The local time evolution of neutral wind uncertainty is presented in Figure 8 with the same format as Figure 6.At the equatorial and low-latitude regions, larger uncertainty appears mainly at 2200 LT and 0000 LT between 0 and 180°E, at 0400 LT between 60 and 120°E, and at 0600 LT between 180 and 300°E.Uncertainty in the meridional direction is more scattered in longitudes.We also notice that the regions of larger uncertainty shift in longitude sectors with time.For example, the larger uncertainty at 2000 LT around 120°E in the southern hemisphere shifts to around 180°E at 2200 LT, to around 240°E at 0000 LT, and to around 300°E at 0200 LT.This is a signature of the longitudinal structure of electron density modulated by the tides in the lower atmosphere.

Sensitivity Analysis
In Figure 10, we present the sensitivity analysis results in the form of the product of the normalized standard ) and Sobol' indices for IMF Bz, solar wind density, and solar wind speed in figure.We do so to highlight the regions of large uncertainty and Sobol' index.We limit our presentation to results at 0500 and 1100 UT, as two examples.At 0500 UT, the sensitivity value for solar wind speed appears in a large region with larger magnitudes (especially at 250 km) than that for IMF Bz and solar wind density.At 1100 UT, we see the high-value region in the evening sector (2000-2400 LT), especially at 250 km, indicating the dominant contribution from IMF Bz to the uncertainty of electron density.We also see a region with great values in the dawn sector (0400-0600 LT) for solar wind speed, especially at 250 km.We find that the uncertainties of electron density in the F region are larger in the nighttime than in the daytime and are the largest during dusk and dawn sectors due to the variation of solar wind parameters.We also notice that the uncertainty shows longitudinal dependence with four peak structures.The local time variation of standard deviations has also been reported in previous numerical and observational analyses (Fang et al., 2018).Fang et al. (2018) showed that the variability of TEC and vertical plasma drift associated with the geomagnetic activity is stronger in the nighttime, especially in the postmidnight sector, than in the daytime at low latitudes.This could be associated with the lower background plasma density at night, and the disturbance dynamo electric field and prompt penetration electric field can change the plasma density and drift to a larger extent.The longitudinal  2020) showed that the quiet-time day-to-day variability of E × B drift is largest at dawn and is mainly driven by the E region winds.In this study, as the lower atmosphere is constant, the variation of the solar wind could impact the equatorial ionosphere through the penetration electric field in a short time (less than 3 hr), the disturbance dynamo electric field in the long term (more than 20 hr), and the traveling ionospheric disturbance (TID) caused by the traveling atmospheric disturbance (TAD) induced by Joule heating at high latitudes (see, e.g., Bagiya et al., 2011;Chakrabarty et al., 2015).
We found that the uncertainties of electron density and plasma drifts show a longitude structure with four peaks at fixed evening local times (2000 LT or 2200 LT).This indicates a close correlation with the longitude structure of the equatorial ionization anomaly (EIA) reported in previous studies (Chen et al., 2016;Lin et al., 2007).The vertical variation of uncertainties in electron density also shows the signature of EIA.We present the mean electron at 2200 LT at 250, 300, and 350 km in Figure 11.We see a clear correlation between the low electron density region and the large uncertainty regions shown in Figure 4. We also observe a high correlation between the uncertainties of the equatorial electron density and zonal plasma drift.The question is then: what controls the longitude dependence of the uncertainty of the zonal plasma drift, hence the longitude dependence of the uncertainty of the electron density.B. G. Fejer et al. (2008) study the longitudinal dependence of disturbance vertical plasma drifts and show that the disturbance dynamo downward drift is largest in the eastern hemisphere in the dusk sector, and the disturbance dynamo upward drift is largest in the late night in the western hemisphere.B. G. Fejer (2011) further shows that the longitude dependence of the disturbed plasma drift is associated with the longitude location where the high latitude energy deposition takes place.This is also revealed in the uncertainty of meridional winds in Figure 7.This again confirms the importance of the universal time variation of solar wind drivers.

Solar Wind Driver and Uncertainty of the IT System
The two examples for sensitivity analysis in Section 3.5 show that solar wind speed plays a dominant role in the electron density uncertainty at 0500 UT, while IMF Bz plays a dominant role at 1100 UT.This difference indicates the importance of the universal time variation of the solar wind drivers.When we look at the solar wind drivers in Figure 2, IMF Bz at 0500 UT is mainly positive (northward), while it is around 0 at 1100 UT.Northward IMF Bz indicates weaker geomagnetic activity than zero or southward IMF Bz, when the solar wind speed will become important.WAM-IPE incorporates the empirical high-latitude electric field model (Weimer, 2005) which takes solar wind drives as inputs.The solar wind speed and density will impact the IT system through the solar wind dynamic pressure term (ρv 2 , ρ, solar wind density, v, solar wind speed) in the model.Apparently, solar wind speed will play a larger role as apposed to solar wind density.The solar wind pressure will further affect the penetration electric field and therefore the low-latitude and equatorial IT system.We also see smaller uncertainties of plasma drifts and neutral winds at 0500 UT in Figures 5 and 7, respectively.Therefore, the universal time variation of solar wind drivers and the polarity of IMF Bz are of critical importance to the uncertainty of the IT system.A model simulation study made by Greer et al. (2017) showed that the strongest enhancements of TEC occur in the American and Pacific sectors when the storm onset happens during 1600-2400 UT.While the storm onset is

Correlation Between Plasma Drift and Electron Density
As mentioned above, the uncertainty of electron density shows a closer correlation with the uncertainty of zonal drift as apposed to vertical drift.We performed a correlation analysis via a linear fitting with the uncertainties of zonal and vertical plasma drifts as input and the uncertainty of electron density as output.The results corresponding to the three locations above (283.5°E,163.5°E, and 43.5°E) are presented in Figure 12.We ignored the data points with the standard deviation of electron density was smaller than 0.01 in the logarithmic scale.
From these results, we see that the electron density is generally correlated more with zonal drift than with vertical drift, especially when the electron density uncertainty is large and in the dusk sector.This indicates that the electron density uncertainty is more closely associated with the zonal plasma drift.Previous studies showed that the equatorial zonal plasma drift is closely associated with the radial penetration electric field (Immel et al., 2004) at high latitudes.
The role of equatorial zonal plasma drift has been comprehensively studied in Pavlov and Pavlova (2007), Pavlov and Pavlova (2008), Pavlov et al. (2008), andPavlov andPavlova (2013) under different conditions.Pavlov and Pavlova (2007) showed that the zonal drift could modify the nighttime equatorial ionosphere and lead to an up to 2.4-fold underestimation of NmF2 in equinox high solar activity conditions.Huang et al. (2010) further showed that the variation of the equatorial ionospheric ion density has an in-phase correspondence with the ion eastward velocity and an anti-phase correspondence with the ion upward velocity in the dusk sector.In Figure 12 of this study, we observe a high correlation between the electron density and zonal plasma drift and a smaller correlation between the electron density and vertical plasma drift.This might indicate that the zonal plasma drift has a local effect on the electron density, while the vertical plasma drift has a non-local effect on the electron density.The other possible reason might be that the local effect of vertical plasma drift on the equatorial electron density is countered by the effect of meridional wind, which may in turn increase or decrease the electron density depending on the direction of meridional wind.
We also note that the WAM-IPE model is based on one-way coupling in that the variability of plasma flow due to the direct impact of the prompt penetration electric field will not feed back to the neutrals.This indicates that the variability of the neutral winds might be underestimated and the variability of the plasma flow might be overestimated.Due to the model, we cannot address or quantify the overestimated or underestimated uncertainty in this study.

Conclusion
We have conducted a comprehensive analysis of the uncertainties of the electron density, plasma drifts, and neutral winds by applying advanced data-driven modeling, UQ, and sensitivity analysis methods to the WAM-IPE simulation in the presence of varying solar wind conditions.We provide insight on various attributes of such uncertainties in terms of local time, altitude, longitude, and latitude dependence due to the varying solar wind drivers.The key findings are summarized as follows: 1.The uncertainties of the equatorial electron density, plasma drifts, and neutral winds are larger at night.The uncertainties of zonal drift are mostly larger than those of vertical drift.The universal and local time evolution at different altitudes indicates a close correlation between the uncertainties of electron density and zonal drift.
Our analysis further indicates that the electron density is more strongly correlated with the zonal drift than with the vertical drift.2. The sensitivity analysis via Sobol' indices indicates the importance of the universal time variation of solar wind drivers and the dominant role of the IMF Bz polarity in the uncertainties of the equatorial and low-latitude IT system.When IMF Bz is northward, the solar wind speed will play a larger role, while when it is 0 or southward, IMF Bz will play a larger role.The combination of IMF Bz and solar wind speed can determine the cross-polar cap electric potential and, therefore, the penetration electric field, which directly impacts the ionospheric state.The Joule heating related to the electric field can further impact the thermospheric state through the propagation of neutral wind.3. The average state determined by the fixed F10.7 and lower atmospheric conditions can modulate the solar wind-associated variability of the IT system.The spatial (longitudinal, latitudinal, and vertical) distribution of the electron density uncertainty at fixed local times in the evening sectors shows modulation of the EIA, with larger uncertainty occurring at regions of low density (see Figures 4 and 11).The universal time evolution of the uncertainty of meridional wind also shows a propagating trend from one hemisphere to the opposite hemisphere.The local time evolution of the uncertainty of meridional winds shows a shift in the longitude direction, which is again associated with the longitudinal structure of the IT system modulated by tides propagated from the lower atmosphere.
In addition, we find that larger uncertainties in the IT system from variations in solar wind mostly occur after 2000 LT, and the uncertainty of the vertical drift around sunset hours is small.This indicates that the uncertainties of the IT system due to the quiet-time solar wind variation might not be the only factor responsible for the day-to-day variability of postsunset ESF occurrence.However, a separate analysis relating the magnitude of this variability to the generation of ESF would be needed in order to reach this conclusion.Further analysis, considering the variability from the lower atmospheric tidal forcing, is needed and will form the basis of a future study.

Appendix A: Multichannel Variational Auto-Encoder (MCVAE)
A brief description of the structure of MCVAE and how it works is given in this section.As we mentioned above, MCVAE is an extension of vanilla VAE, we will first introduce vanilla VAE and then how it is used to build MCVAE.
ZHAN ET AL.In Figure A2, MCVAE is an extension of VAE with multiple input data sources (x 1 , x 2 , x c ) and independent encoder and decoder for each data source.The difference is that the reconstruction loss is the average of all the channels which can force the correlation among the input data sources reserved.Interested readers are referred to Antelmi et al. (2019b) for the basics of MCVAE.In our case, we have 3 input data sources (IMF Bz, solar wind density, and solar wind speed), and each data source has 288 values per day and all the days with solar wind measurements available during 1981-2021.The encoder for each channel is a three-layer neural network with 288, 149, and 10 values in the first, second, and third layers, respectively.The decoder for each channel is a similar neural network but with the number of values in each layer reversed (10,149,288).When generating synthetic samples, 3 independent vectors of ξ following Gaussian distributions with length 10 are fed to the decoders.
The comparisons between historical measurements and generated drivers during quiet conditions are shown in Figure A3 and Figure A4. Figure A3 shows the histograms of measurements on the top and synthetic drives on the bottom.

Figure 1 .
Figure 1.Diagram for the whole workflow of this study.On the left, it shows the procedure to apply MCVAE to generate new solar wind drivers.On the right, it shows the key elements for UQ analysis by using the WAM-IPE simulations and PCE.
) for more detailed exposition.

Figure 2 .
Figure 2. (gray lines) Variations of the 500 samples of IMF Bz (top), solar wind density (middle), and solar wind speed (bottom) to drive WAM-IPE and the corresponding mean (red) and one standard deviation range (blue).

Figure 3
Figure3shows the uncertainty results for electron density from four chosen universal times (UTs) as a function of local time and latitude at 250, 300, and 350 km.The uncertainties are calculated as the standard deviation of the logarithmic-scale electron density.The white lines label the geomagnetic equator and 20°N and 20°S.

Figure 3 .
Figure 3. Universal time evolution of the electron density uncertainty as a function of local time and latitude from 4 different UTs (top: 0500, 1100 UT, bottom: 1700, 2300 UT) and at 3 different altitudes (250, 300, and 350 km, labeled in the upper left of each subplot).The white lines indicate the magnetic equator and 20°N and 20°S.
equatorial regions can be impacted by the solar wind-associated electric field and Joule heating associated disturbance dynamo electric field

Figure 4 .
Figure 4. Local time evolution of the electron density uncertainty as a function of longitude and latitude at different local times (LTs) from top to bottom (2000, 2200, 0000, 0200, 0400, and 0600 LT) at three different altitudes (250, 300, and 350 km) from left to right.The discontinuity of the data (sharp change of values at certain longitudes) is due to the simulation's start and end.
that at 1100 UT.The difference is that large uncertainty also appears during 0400-0500 LT.At 2300 UT, the vertical drift only shows small uncertainty during 0600-0800 LT and much smaller uncertainty during 2000-2200 LT.The zonal drift again shows large uncertainty during 2100-0100 LT and during 0400-0600 LT.

Figure 5 .
Figure 5. Universal time evolution of the uncertainties of vertical (left) and zonal (right) plasma drifts at 0500, 1100, 1700, and 2300 UT (from top to bottom) at 300 km.The white line indicates the geographic equator, and the white lines indicate the magnetic equator and 20°N and 20°S.

Figure 6 .
Figure 6.Local time evolution of the uncertainty of vertical (left) and zonal (right) plasma drifts at 2000, 2200, 0000, 0200, 0400, and 0600 LT (from top to bottom) at 300 km.The magenta line indicates the geographic equator, and the white lines indicate the magnetic equator and 20°N and 20°S.

Figure 9
Figure9presents the local time and vertical variation of electron density, plasma drift, and neutral wind uncertainties from three fixed locations along the magnetic equator.The first location is the location of Jicamarca (283.5°E) while the other two (163.5°Eand 43.5°E) are 120° away in longitude.The results are presented as a function of altitude (vertical pressure level for neutral winds) and local time.From the first row, we observe that the uncertainty of electron density at Jicamarca is large in the evening sector between

Figure 7 .
Figure 7. Universal time evolution of the uncertainty of zonal (left) and meridional (right) winds from 4 different UTs (from top to bottom).

Figure 8 .
Figure 8. Local time evolution of the uncertainty of zonal (left) and meridional (right) winds from 4 different LTs (from top to bottom).

Figure 9 .
Figure 9. Uncertainty of electron density in logarithmic scale (first row, std(log 10 N e )), zonal drift (second row), vertical drift (third row), zonal wind (fourth row), and meridional wind (bottom row) from three locations (283.5°E,163.5°E, and 43.5°E) along the magnetic equator.The vertical axes for zonal and meridional winds are pressure levels.

Figure 10 .
Figure 10.Global distribution of sensitivity analysis results for the total effect Sobol' index at 0500 UT (first 3 rows) and 1100 UT (bottom 3 rows) as the product of the normalized standard deviation and Sobol' index for IMF Bz (top), solar wind density (second row), and solar wind speed (bottom) at 250 km (left), 300 km (middle), and 350 km (right).The white lines indicate the magnetic equator and 20°N and 20°S.

Figure 12 .
Figure 12.Correlation between the electron density and plasma drifts (left: vertical, right: zonal) from the three locations (283.5°E,163.5°E, and 43.5°E) at 250, 300, 350, 400, and 500 km (from top to bottom) selected from Figure 9.The red lines correspond to the linearly fitted curves.The correlation coefficients, r 2 , and root mean squared error (RMSE) are labeled in each subplot.Only uncertainties larger than 0.01 are considered for correlation analysis.
FigureA1, VAE consists of an encoder to infer a lower dimensional representation (the mean μ and standard deviation σ of a resulting multivariate Gaussian random vector) of the data in a probabilistic way and a decoder to transform the latent representation back to the original observation space by sampling the resulting decoder output distribution.By enforcing a prior on the latent variable z  ( =  +  ⊙ ,  ∼  (, )) , the VAE can generate new samples of x by sampling the latent variable and evaluating it through the decoder.Interested readers are referred toKingma and Welling (2019) for the basics of VAE.
Figure A4 shows the scatter plots with Speed versus IMF Bz on the left, Density versus IMF Bz in the middle, and Speed versus Density on the right.The distributions indicate that the ranges and shapes of the distributions of the synthetic drivers are very similar to those of the measurements.The similar shapes of scatter plots indicate that the correlation among the solar wind parameters is reserved.

Figure A1 .
Figure A1.Structure of variational autoencoder.x is the input data and  x , g ϕ (x) and h θ (z) is the neural network for encoder with ϕ and θ as the model parameters, respectively.

Figure A2 .
Figure A2.Structure of multi-channel variational autoencoder.c represents the number of channels or data sources.z is the shared common latent space.

Figure A3 .
Figure A3.Distributions of solar wind parameters of historical measurements (top) and synthetic samples (bottom).

Figure A4 .
Figure A4.Scatter plots of solar wind parameter pairs (from left to right: solar wind speed vs. Bz, solar wind density vs. Bz, and solar wind speed vs. solar wind density) of historical measurements (top) and synthetic samples (bottom).