Accounting for Variability in ULF Wave Radial Diffusion Models

Many modern outer radiation belt models simulate the long‐time behavior of high‐energy electrons by solving a three‐dimensional Fokker‐Planck equation for the drift‐ and bounce‐averaged electron phase space density that includes radial, pitch‐angle, and energy diffusion. Radial diffusion is an important process, often characterized by a deterministic diffusion coefficient. One widely used parameterization is based on the median of statistical ultralow frequency (ULF) wave power for a particular geomagnetic index Kp. We perform idealized numerical ensemble experiments on radial diffusion, introducing temporal and spatial variability to the diffusion coefficient through stochastic parameterization, constrained by statistical properties of its underlying observations. Our results demonstrate the sensitivity of radial diffusion over a long time period to the full distribution of the radial diffusion coefficient, highlighting that information is lost when only using median ULF wave power. When temporal variability is included, ensembles exhibit greater diffusion with more rapidly varying diffusion coefficients, larger variance of the diffusion coefficients and for distributions with heavier tails. When we introduce spatial variability, the variance in the set of all ensemble solutions increases with larger spatial scales of variability. Our results demonstrate that the variability of diffusion affects the temporal evolution of phase space density in the outer radiation belt. We discuss the need to identify important temporal and length scales to constrain variability in diffusion models. We suggest that the application of stochastic parameterization techniques in the diffusion equation may allow the inclusion of natural variability and uncertainty in modeling of wave‐particle interactions in the inner magnetosphere.


Introduction
The Van Allen outer radiation belt is a typically quiescent torus-shaped region in near-Earth space between 13,000 and 40,000 km radial distance consisting mainly of electrons between 100s of keV and multiple MeV trapped by the Earth's geomagnetic field.Protons are also present and modeled in the radiation belts (Vacaresse et al., 1999), but here we focus on the high-energy electron population.The behavior of electrons in the outer radiation belt is affected by multiple processes, some of which are immediate responses to solar wind forcing, whereas some are more indirect energy pathways involving energy stored in the substorm cycle.Numerical modeling is a powerful tool to provide deep understanding of the behavior of the outer radiation belt, allowing us to quantify the effects of different processes (e.g., Glauert et al., 2014;Reeves et al., 2012;Shprits et al., 2008).
From a more practical standpoint, the ability to model these physical processes is becoming increasingly important as Earth becomes more dependent on space-based technologies.As of 31 March 2020 there were 135 satellites operating in medium Earth orbit (MEO; 2,000-35,786 km) and 554 in geostationary orbit (GEO;35,786 km), therefore operating in the heart of the belt (https://www.ucsusa.org).Outer radiation belt electrons can be hazardous to these spacecraft, but there are insufficient in situ measurements available to monitor the radiation environment directly.There remains a pressing need to develop accurate models of the outer radiation belt for operational purposes in addition to promoting further physical understanding.
One effective method to study the dynamics of the outer belt electrons is to model the evolution of electron phase space density (PSD) f(M, J, Φ; t) by a Fokker-Planck equation as a function of the three adiabatic invariants and time (Schulz & Lanzerotti, 1974).Here M, J, and Φ are the first, second, and third adiabatic invariants, respectively.It is helpful to consider Φ in terms of the adiabatic reference parameter L*, defined by (Roederer, 1970).Since a first-principles model of wave-particle interactions in the outer radiation belt is intractable across its large volume and long timescales, all the physics within the outer radiation belt can be effectively described by diffusive processes.Each type of diffusion-pitch angle, energy, and radial-by each wave mode is described in the Fokker-Planck equation by a diffusion coefficient D ij .A myriad of different wave-particle interactions is important for the radiation belts.For example, very low frequency (VLF) whistler mode chorus mediate energy diffusion (Thorne et al., 2013), whereas VLF whistler mode hiss (Lyons & Thorne, 1973;Meredith et al., 2007) and ULF electromagnetic ion cyclotron (EMIC) waves (Kersten et al., 2014) predominantly diffuse in pitch-angle and therefore contribute to loss.ULF wave-driven radial diffusion at Pc-5 frequencies is considered to be an important and effective mechanism to transport and accelerate relativistic electrons in the outer radiation belt (Elkington et al., 2003;Mann et al., 2013;Ozeke et al., 2017Ozeke et al., , 2018;;Shprits et al., 2008).
In this paper we focus on radial diffusion as a result of ULF waves, which in the diffusion framework can be modeled as a straightforward one-dimensional problem.All of the physics is contained in the radial diffusion coefficient D LL , which is proportional to ULF wave power.A wealth of data exists both on the ground and in space to calculate ULF wave power and construct D LL (Dimitrakoudis et al., 2015;Li et al., 2017;Liu et al., 2016;Ozeke et al., 2012Ozeke et al., , 2014;;Ukhorskiy et al., 2009).Empirical models formulate analytic expressions for D LL from ULF wave power data over long timescales, aiming to capture the spatiotemporal evolution of D LL in such a way that although rapid changes cannot be accurately captured, the long timescale behavior of the outer radiation belt may be adequately described (e.g., Ozeke et al., 2018).In this paper, we wish to highlight the numerical consequences of using different methods for modeling the temporal and spatial variability of D LL with more realistic values that represent the underlying probability distribution of ULF wave power.
Many theoretical approximations exist for the radial diffusion coefficient D LL based on a variety of assumptions and approximations (Ali et al., 2016;Birmingham, 1969;Cornwall, 1968;Elkington et al., 2003;Fälthammar, 1966Fälthammar, , 1968;;Fei et al., 2006;Lejosne et al., 2013;Liu et al., 2016;Schulz & Lanzerotti, 1974).All of these approximations are constrained by some statistical parameterization of ULF wave power obtained from many years of space or ground-based observations.The most widely used D LL parameterizations in radiation belt models parameterize by the geomagnetic index K p (Brautigam & Albert, 2000;Ozeke et al., 2012Ozeke et al., , 2014)).These parameterizations are deterministic with a single output for each value of Kp.
Typical approaches in radiation belt modeling follow a classical parameterization approach whereby average or median D LL values are used.These values only change when the fit parameters change, and therefore, there is a chance that the full range of variability of D LL is not captured in this classical approach.In numerical weather prediction and climate modeling, classical parameterizations have proven to be insufficient.Instead, stochastic parameterizations are used to capture the whole distribution of behavior in underlying physical processes to yield improved results.Note that previous attempts to capture more realistic variability in ULF-mediated radial diffusion have used observations to recreate event-specific models of diffusion (Perry et al., 2005;Riley & Wolf, 1992;Tu et al., 2012).These types of study, although potentially more accurate, are limited to test cases with available data in space and time.We propose that in cases where direct data is lacking, it is still possible to capture the full range of behavior in the problem using stochastic parameterizations (e.g., Watt et al., 2017), and we demonstrate a simple implementation of this technique in this paper.
Here we present a series of idealized numerical experiments of radial diffusion over a hypothetical period of constant geomagnetic activity.These experiments offer a proof of concept intended to explore the spatiotemporal impacts of including stochastic variability in comparison with the (Ozeke et al., 2014) ULF radial diffusion coefficients in the radial diffusion equation and highlight current deterministic model limitations.Any significant discrepancies between the deterministic and stochastic models should motivate further research questions to better understand the physical processes underlying ULF wave-driven radial diffusion to include in our models for improved accuracy.The remainder of this paper is structured as follows.Sections 2-4 describe the radial diffusion problem, implementation of stochastic parameterization, and setup and description of the idealized experiments, respectively.Section 3 presents the results from the numerical experiments.Section 4 discusses the impact of the results in the wider context of the outer radiation belt.Section 5 describes conclusions and remarks from this paper.

Modeling the Radial Diffusion Equation
We focus on the radial diffusion equation as a simplified approximate model of electron behavior in the outer radiation belt.Although the one-dimensional description of radial diffusion has successfully reproduced electron behavior during some events (e.g., Ozeke et al., 2018;Shprits et al., 2005), the diffusion framework itself is not always accurate.Previous studies have calculated radial diffusion coefficients directly in "event-specific" analysis (e.g., Ukhorskiy et al., 2009) and demonstrate that diffusion-based models can have difficulty accurately rendering event-specific dynamics (Ukhorskiy et al., 2009).Here, we intend these numerical experiments as a straightforward demonstration of the concept of stochastic parameterization.Radial diffusion is also a valid and important part of more complicated outer radiation belt models, where it is joined by diffusion processes in velocity space due to other wave modes.Over the long timescales studied in diffusion models, we observe that empirical models for D LL , in whichever theoretical framework they are constructed, naturally have some uncertainty.Investigating the consequences of that uncertainty is our aim in this work.
In this demonstration we simplify the behavior of high-energy electrons in the outer radiation belt and focus on radial diffusion across Roederer L* (Roederer, 1970), hereon denoted L. Here, the first and second adiabatic invariants, M and J, are conserved.The evolution of the distribution function of trapped particles f(M, J, Φ; t) can be related to the distribution function at time t + Δt (without sources or sinks) (1) where Π(Φ − ϕ, ϕ, t) is the probability that a particle with an invariant shell coordinate Φ − ϕ at time t will end up with coordinate Φ at time t + Δt.By Taylor expanding f, Π to first order in t on the left and second order in Φ in the integral, we obtain the one-dimensional Fokker-Planck equation Here D Φ and D ΦΦ are the firstand second-order Fokker-Planck diffusion coefficients, respectively.If we assume the following relation for D Φ , the average change of Φ per unit time for one particle on the shell Φ during that time interval For radial diffusion to be effective, a radial gradient in the PSD is required, which we assume here.A precipitation loss term is often also added to Equation 4, which is ignored here in the idealized case.Radial diffusion is considered across L = 2.5-6.Dirichlet and Neumann boundaries are imposed on the inner and outer boundaries, respectively: ∀t; (5) In reality the gradient across the outer boundary will not be 0, and many radiation belt models either determine the outer boundary from electron flux data observed by spacecraft (e.g., Drozdov et al., 2017;Glauert et al., 2018;Shin & Lee, 2013) or use plasmasheet characteristics (Christon et al., 1988(Christon et al., , 1991) ) and magnetic activity dependencies (Bourdarie & Maget, 2012) for analytic fits (Maget et al., 2015).
In Equation 4, D LL represents the ULF wave radial diffusion coefficient.Constructed through a coordinate transform of the flux invariant diffusion coefficient, D ΦΦ , D LL is formally defined by (Roederer & Zhang, 2014) where R s , ΔR s /R s , and τ d are the dipole-distortion parameter, its relative fluctuation, and the drift period, respectively.Here, <> denotes the drift-average operator.In a realistic setting, R s would be represented by a parameter that globally describes magnetospheric activity, such as Kp or ULF wave power.Application of different frameworks to describe large-scale fluctuations of electric and magnetic fields (e.g., Brautigam & Albert, 2000;Brautigam et al., 2005;Lejosne et al., 2013;Ozeke et al., 2012Ozeke et al., , 2014) ) employ different assumptions, but many ultimately require some estimate of the power spectral density of ULF fluctuations in electric and/or magnetic fields.We note that from Equation 7 and from theoretical estimates of D LL , there are inherent minimum temporal scales on which D LL is constructed: by definition D LL is constructed for timescales longer than the drift period of the electrons, longer than a few periods of the ULF wave fluctuations, and of the same order or longer than the solar wind driving processes that induce the ULF fluctuations.In many cases, ULF power spectral density is estimated from observations over a period of at least an hour (see Ozeke et al., 2014), and so we employ this as the smallest timescale of variability in our study.
We consider as a deterministic reference model the empirical L and K p parameterized D LL presented by Ozeke et al. (2012Ozeke et al. ( , 2014)).This model is a simplification of the theoretical analysis presented by Fei et al. (2006) and assumes that median ULF wave power is representative of expected ULF wave power.The most notable feature of this model is that the uncertainty in the statistical representation of ULF power spectral density has been quantified, allowing us to perform this demonstration using observationally derived constraints.Other models exist, which are similarly parameterized by Kp activity, with some following the same theoretical framework as Fei et al. ( 2006) (e.g., Brautigam et al., 2005) and others pursuing other frameworks (e.g., Lejosne et al., 2013), but all do not explicitly state and characterize the uncertainty in their models as in Ozeke et al. (2012Ozeke et al. ( , 2014)).We note that the accuracy of the theoretical framework used to estimate D LL is beyond the scope of this paper and direct the interested reader toward Lejosne (2019) for a thorough review of such frameworks.We reiterate that since the (Ozeke et al., 2014) empirical D LL model contains explicit estimates of uncertainty, that makes it appropriate for use in our demonstration.
Since the azimuthal electric field radial diffusion coefficient, D E LL , typically dominates, in these idealized experiments we omit the compressional magnetic component and base our stochastic parameterization around the model for D LL ¼ D E LL , expressed per day by We describe in the following section how we implement our estimates of D E L L ðtÞ, by perturbing Equation 8 in such a way as to recover a better representation of the underlying distribution of D E L L across a period of time.

10.1029/2019JA027254
Journal of Geophysical Research: Space Physics We solve the radial diffusion equation using a modified Crank-Nicolson second-order finite difference scheme presented by Welling et al. (2011), which is semi-implicit and unconditionally stable: where L 2 for modeling simplicity.The chosen grid and time steps for our numerical experiments are 0.1L and 1 s, respectively, following extensive model verification of the numerical scheme to determine a suitable trade off between numerical error and computational cost for the experiments (see the supporting information).

Stochastic Parameterization
We suggest that the most physically intuitive method to implement stochastic parameterization is to focus efforts on the representation of the diffusion coefficient, since it is the variable that contains all the information about the wave-particle interaction.The diffusion coefficient parameterization has been shown to result in a large amount of variability, especially during storm times (Murphy et al., 2016).In this work, we choose a straightforward method to model D LL (L,t) that involves constructing a noisy temporal or spatial series that retains the key known properties of the distribution of D LL .More sophisticated techniques, such as autoregressive moving average (ARMA) models, can be used to create spatiotemporal series of the diffusion coefficients with the appropriate autocorrelative properties.However, these rely on important characteristic scales of spatial and temporal variability that are not yet known.
We do, however, have access to some information constraining the expected distribution of D LL .Bentley et al. (2018) found that the probability distribution of ground-based ULF wave power appears log-normal (LN).We infer from this that D LL is also likely to be approximately LN; indeed, Ozeke et al. (2014) confirm that the distribution of D LL in space is not Gaussian and is log-symmetric, since the interquartile range (IQR) is reported between one third and three times the median.Hence, it is appropriate to construct a noisy time series for D LL by multiplying the median D LL by a random LN noise factor ϵ, resulting in a time series that, when aggregated over a long period of time, reproduces the required LN distribution.If we constructed a noisy temporal or spatial series by adding Gaussian noise to the median D LL , the resulting distribution of D LL cannot be LN since it has the potential to include negative values of diffusion, which would also be difficult to interpret in this context.
To investigate the consequences of variability, we consider ensembles of numerical experiments.In each case we compute the solutions of the radial diffusion equation using Equation 9, where D LL (t) is separately constructed each time using the methods described below.Our recreations of D LL (t) do not alter the underlying Fokker-Planck diffusion theory but produce realizations of D LL that better recover the underlying distribution of ULF power spectral density.Future work will seek to identify the most appropriate methods to model both the diffusion coefficient and its variability, but the straightforward methods we adopt here serve to illustrate the behavior of the radial diffusion equation when stochastic parameterization is adopted using known constraints.
In each numerical experiment we run an ensemble with 250 ensemble members, providing a span of possible realizations of 48 hr D LL time series resulting from the inclusion of a stochastic variability.Convergence testing of our numerical experiments (see the supporting information) demonstrates that 250 ensemble members is sufficient to realize the behavior of the experiment.
In all experiments we choose K p = 3, corresponding to "unsettled" geomagnetic activity.Unsettled geomagnetic activity allows us to explore stochastic variabilities during periods where the radial diffusion coefficients are large enough to see changes after 48 hr.We also wish to avoid the illogical situation of having a very high level of geomagnetic activity while enforcing a constant outer boundary.For the demonstrations approximated in this paper, a compromise of K p = 3 was felt to be appropriate.The initial PSD is chosen to provide a peak inside the computational domain as expected in the outer radiation belt, and a zero gradient at the outer boundary, for ease of computation in these illustrative experiments where we have chosen A = 9 × 10 4 , μ = 4, σ = 0.38, B = 0.05, and γ = 5 and erf is the error function.Such a profile is reasonable when compared to satellite observations (e.g., see Figures 1 and 2 in Boyd et al., 2018).
If one wanted to do the equivalent in L space (with a transformed diffusion equation), it suffices to use (Roederer & Zhang, 2014) The initial PSD profile and proposed boundary conditions result in the expected radial diffusion process drawing PSD from central L toward both boundaries.

Experiment 1: Temporal Variability of D LL
Our first experiment focuses on the temporal variation of D LL across a range of timescales.We employ a simple method, where the D LL in Equation 8 is multiplied by a random factor ϵ, which changes every Δt.The same factor ϵ is applied at each value of L in the model.The choice of distribution of ϵ is guided by the statistical analysis presented by Ozeke et al. (2014), who found that the IQR of observed wave power implies that D LL lies between a third of and three times the model value 50% of the time.We use this information to control the variance of the noise.Combined with recent studies that suggest that ULF wave power spectral densities appear LN (Bentley et al., 2018), we construct a log-normally distributed variability with the following parameters: where

Experiment 2: Spatial Variability of D LL
In Experiment 1, D LL was constructed with perfect correlation across all L, with the same ϵ applied to all Lshells.This is one extreme of L spatial correlation, with the (Ozeke et al., 2014) D L L scaling as a smooth, monotonically increasing profile.We hereon refer to this approach as global variability.However, we must consider that although the statistical profile of D L L (L) is smooth, individual cases of D L L (L,t) may be less smooth.In this experiment, we investigate how radial diffusion responds to a realized D L L , which may vary on local spatial scales, and not necessarily be a smooth monotonically increasing function of L.
We now consider the log-normally distributed variability applied every 3 hr, comparing the global variability with local spatial correlation scales.We consider cases where D L L varies independently on spatial scales of 1L,0.5L, and 0.1L.Example ensemble members for each of these cases are shown in Figure 2. The final case denotes the other extreme where measures of D L L (L,t) are independent at all grid points, that is, that independent ϵ is applied at each grid point in L to create an ensemble of D L L both spatially and temporally.
We have retained temporal variability in this experiment to maintain our goal of creating D L L time series that represent realistic values.Ground magnetometer ULF wave power measurements, and consequently D L L , do not typically remain constant over 2 days (e.g., Olifer et al., 2019).Results from differing spatial variability scales can therefore be interpreted in conjunction with the 3-hourly temporal variability.
In a more physical realization, we would expect spatial correlations across L to be less crude and abrupt, and are likely to exhibit smoother variations with appropriate length scales.However, for the purpose of this demonstration, we have chosen the simplest way to apply spatial variability in the model to motivate the importance of understanding the spatial structure of radial diffusion across L. .Larger variances may be necessary if the variability of D LL is not simply due to the variability in observed ground-based ULF power spectral density.Smaller variances have been considered to see the effect of an "improved" parameterization (i.e., one where the parameters are chosen in a way that minimizes the variance).In each of these cases, ensemble D LL time series are formulated by applying variability globally across L every 3 hr, with the distribution of the variability LN.

Experiment 4: Shape of the D LL Probability Distribution
Each experiment (1-3) utilized a log-normally distributed variability, chosen based on statistical studies of ULF wave power spectral densities parameterized by solar wind variables (Bentley et al., 2018).The IQR presented by Ozeke et al. (2014) describes the uncertainty in the deterministic parameterization, but we do not know how the D LL s are distributed in a K p -based model.Adopting the values and log-symmetric nature of the (Ozeke et al., 2014) IQR in order to preserve statistical averages (a zero mean and median in the logarithm), a range of log-symmetric distributions for the variability are tested.We consider log-uniform (LU), LN, log-Laplace (LL) and log-Cauchy (LC) distributions, which provides a set of distributions ranging from bounded to heavy tailed (for further information about each of these distributions, please see the supporting information).Since the heavy tailed distributions can easily produce variabilities resulting in a D LL which is unrealistically many orders of magnitude larger than the deterministic solution, for this experiment we bound the variability by 3 orders of magnitude (i.e., the variability can increase/decrease D LL up to a maximum/minimum of 3 orders of magnitude compared to the reference value).The respective probability density functions (PDFs) of the variability distributions are as follows: for x > 0, where I [,] is the characteristic function.Here the quantities a, b, σ N , σ L , and σ C are the parameters of the underlying uniform, normal, Laplace, and Cauchy distributions, respectively.The parameters were calculated from their corresponding cumulative density functions in order to preserve the IQR specified by Ozeke et al. (2014) (see the supporting information).

Results
The figures showcasing results for each experiment generally follow the same format.The initial PSD and resulting PSD from the constant deterministic D LL are shown.By the log-symmetric nature of the D LL probability distributions in each experiment, the constant deterministic D LL is precisely the median diffusion coefficient from the ensemble and a natural reference for comparison.The mean diffusion coefficient is deliberated in section 6.There is no convention regarding which statistical measure is most appropriate in ensemble modeling (Knutti et al., 2010), and we have therefore shown two natural measures, the ensemble mean and median.By ensemble mean (median) PSDs, we imply the PSD profile resulting from taking the mean (median) across all ensemble members at each L, and not representing a specific member of the ensemble.The kernel density estimates (KDEs) of the ensembles are also shown.Kernel density estimation is a mathematical process of finding an estimate PDF of a random variable, inferring attributes of a population based on a finite data set.In the case of our ensembles, the contribution of each ensemble member value in L-PSD space is smoothed out into a region of space surrounding it.Aggregating each of these smoothed points provides an image of the overall ensemble structure and density function.Ensemble modes, another useful measure of the ensemble result, can be estimated from this density function (Kourentzes et al., 2014).
In our figures KDEs shown are relative to each column, meaning that if a single L column were extracted, the result would be a PDF estimate of the PSD at that particular L. KDEs are therefore useful in an ensemble setting since they allow us to see where ensemble member solutions cluster in the phase space.In our estimates the KDEs are calculated over 100 bins.

Experiment 1-Temporal Scales
Results of the ensembles for the variety of temporal variability scales are shown in Figure 3.For ensemble medians, inclusion of a LN variability results in more diffusion than the constant deterministic D LL at all variability temporal scales less than 24 hr, with the magnitude of diffusion increasing as the temporal scale decreases.The ensemble median for a temporal variability of 24 hr is identical to the deterministic solution, suggesting that on long timescales, a deterministic parameterization of D LL is sensible for a D LL with daily variation.Results for the ensemble mean are similar, except we observe more diffusion than the constant D LL at all temporal scales.This is unsurprising since the (Ozeke et al., 2014) D LL is based on the median of log-symmetric distributions, where means are larger than medians.Therefore, the ensemble D LL time series at all temporal scales will have a mean larger than both the deterministic approximation and ensemble median, resulting in more diffusion.An interesting result lies in the comparison of ensemble medians and means.On the most rapid temporal D LL variability of 1 hr, results from the ensemble mean and median are identical.As the temporal variability becomes less rapid, both exhibit less diffusion, but the profiles separate with the ensemble median displaying increasingly less diffusion than the mean as it approaches the deterministic solution at daily variability.
Over all temporal variability scales, the occurrence of possible states in the set of all ensemble solutions spans similar regions.For the rapid 1 hr variability, the set of all solutions is more diffusive than the deterministic case.The deterministic solution becomes increasingly closer to the denser region of ensemble solutions with larger temporal scales, falling exactly in the region of highest probability for daily variation.We see that increasing the frequency of D LL variability tends to a single mode solution in density, which is more diffusive than that produced by the deterministic model.Inclusion of the variability expressed by Ozeke et al. (2014) in their 3-hourly deterministic model produces a span of solutions, which vary greatly from the deterministic case at all L, most of which are more diffusive.The use of the median-based deterministic parameterization may therefore not be robust.When we allow the stochastic D LL to vary daily, however, the deterministic solution fell exactly in the regions of highest probability, emphasizing again that the deterministic approximation is more suitable for a daily varying D LL .When including variability, the deterministic parameterization frequently produces lower estimates of radial diffusion, so understanding the temporal variability of ULF wave power spectral density is important to know the extent of potential underestimation.

Experiment 2-Spatial Scales
Ensemble results for Experiment 2 are shown in Figure 4. We find that on average all spatial scales of variability result in similar levels of diffusion, but all exhibit more diffusion than the deterministic solution.In each case the ensemble means and medians are almost identical.Most importantly, we observe variance reduction in the set of ensemble solutions as independence of D LL measurements occurs on increasingly smaller spatial scales, with the distributions tending toward a single mode solution of diffusion similar to those exhibited by the ensemble median and mean.A smaller variance implies possibility of a stronger parameterization with reduced uncertainty.It is important to investigate instantaneous observations of ULF wave power across multiple latitudes to better understand spatial correlations and coherence across L*, since regions of independent power measurements could allow for better parameterizations of D LL .

Experiment 3-Variance
Figure 5 shows the ensemble results for Experiment 3, with each variance expressed in terms of the variability IQR.It is evident that radial diffusion is very sensitive to the width of the variability distribution.Just doubling the multiplicative scaling of the IQR suggested by Ozeke et al. (2014) results in significantly more diffusion in both ensemble averages, reducing the peak in PSD by around 20,000.The shape of the distribution for the set of all ensemble solutions also drastically changes, with a large density of solutions tending to the asymptotic result controlled by the boundary conditions.Although a wider variability distribution equally allows for both significantly larger and smaller vales of D LL , the radial diffusion equation is clearly heavily sensitive to the larger values that drive radial diffusion to significant levels beyond the deterministic approximation.
As seen in the other experiments, introduction of any variability regardless of its width results in more diffusion than the deterministic solution, when considering ensemble averages.However, if the uncertainty in the deterministic model were to have a slightly smaller multiplicative IQR of ±2 the (Ozeke et al., 2014) D LL , the variance of all ensemble solutions decreases significantly.With this smaller variance, the ensemble mean and median PSDs are closer to the deterministic model, which also falls within the set of ensemble solutions.This suggests that parameterization of ULF radial diffusion coefficients should prioritize variance reduction in order to be better representative of the underlying physical process, which draws upon the efficiency of binning by geomagnetic index K p , from which most of the uncertainty arises (Ozeke et al., 2014).

Experiment 4-Underlying Distribution
Ensemble results for Experiment 4 are shown in Figure 6.Differences between the heavy and nonheavy tailed distributions are apparent in the ensemble medians.Although studies suggest that ground-based ULF power spectral density is LN when parameterized by solar wind variables (Bentley et al., 2018), the distribution of uncertainty in the K p -based (Ozeke et al., 2014) model is not disclosed.If the distribution were to be heavy tailed or LU (which may be considered to have the heaviest tail as all values in the uniformly distributed component have equal chance of being sampled), we see more than double the median diffusion than for a log-normally distributed variability.For scenarios where the expected ULF wave power is not a statistical average, the assumed LN variability can exhibit as much diffusion as some of the heavy tailed variabilities, but this is more unlikely as shown in the KDEs.In any case, with the inclusion of variability in D LL for all probability distributions, we see significantly more diffusion than the deterministic solution, with notable variance in ensemble solutions for all variability distributions.The heavier tailed variabilities have denser regions approaching that of the asymptotic solution, and the shape of the KDEs across L-shells is quite distorted contrary to the smoothness seen for a LN D LL .Since there are multiple components of interest in the ensemble results, studies investigating the true underlying probability distribution of ULF wave power are vital to quantifying the shortfall and uncertainty introduced by a deterministic empirical D LL based upon statistical averages.

Discussion
In the outer radiation belt, radial diffusion has the ability to both accelerate electrons to relativistic energies and produce fast losses, where the efficiency of the acceleration increases with increasing ULF wave activity (Elkington et al., 2003;Shprits et al., 2008).Many models use an empirical deterministic radial diffusion coefficient dependent on L and K p , which may sacrifice accuracy (Brautigam & Albert, 2000;Brautigam et al., 2005;Ozeke et al., 2012Ozeke et al., , 2014)).In this paper we present idealized numerical experiments, which investigate the impact of including variability in the radial diffusion equation.Our experiments reintroduce the variability into a parameterized model, where D LL has been binned by Kp.We use the observationally constrained variability in the model to model a variable D LL that reproduces a realistic distribution of values and compare against the constant parameterized value.We employ constant boundary conditions and only study one value of the controlling parameter Kp.In this way, we isolate only the variability of D LL due to its parameterization by Kp.
In all experiments we found that the mean and median of the ensembles exhibit increased diffusion above that for the deterministic approximation.One way to interpret these results is that when the likelihood of strong radial diffusion is large over a particular period (either because the variance in the parameterization is large or because the underlying distribution has a heavy tail), then the diffusion exceeds what one would expect from using a constant diffusion coefficient.It is important to bear in mind that the times where diffusion is weak will not counteract the times when diffusion is strong because there is no means of reversing the diffusion; hence, the periods when diffusion is much stronger than the median will dominate the temporal evolution of the experiment.When the diffusion varies more rapidly, then each member of the ensemble is more likely to contain a period of strong diffusion over the fixed 48-hr experiment length, thus contributing to a stronger diffusion in the mean/median of the ensemble.The ensembles are also sensitive to the size of the variance (see Experiment 3), again suggesting that it is the likelihood of ensemble members containing periods of very strong diffusion that dominates the ensemble results.
The collected range of numerical experiments suggests that over extended time periods, infrequent instances of very efficient ULF wave-particle interactions make important contributions to radial diffusion and should be included in models in some way.We also note that by using an ensemble framework, the uncertainty in the PSD is explicitly quantified, providing the means to provide a range of confidence in the model for more accurate radiation belt modeling.The quantification of uncertainty in D LL is also important for future data assimilation methods.
Experiment 1 indicates that the amount of diffusion depends upon how rapidly the diffusion coefficient varies.Hence, it is important to understand the timescales of variability.ULF wave power can vary on a range of timescales, which would ideally be accounted for in the radial diffusion coefficient.For example, ULF wave power can increase and persist on the order of tens of minutes during an auroral activation due to substorms (Rae et al., 2011), while decaying on hourly timescales during strong poloidal wave events (Liu et al., 2011).Parameterization of D LL with K p may therefore not be optimal, since it may not vary quickly enough.
We found that variation of D LL with the added inclusion of local spatial variabilities on a range of length scales resulted in more diffusion that the deterministic solution (see Experiment 2).However, when considering the ensemble averages, all levels of spatial coherence across L* performed similarly.Since applying variability to subglobal spatial scales still allows for an enhanced D LL at several L, this result is somewhat counterintuitive to those found in the other experiments.While it was found that instances of weaker diffusion cannot counteract the temporal evolution imposed by instances of stronger diffusion, counteractions can occur across spatial scales, creating a net diffusion that seems to follow that observed by a globally applied variability.More interestingly, we found that the variance of the possible states in the set of all ensemble solutions decreases significantly with variability applied to increasingly smaller subglobal spatial scales.It is important to understand and quantify these spatial scales.Rae et al. (2019) showed the evolution of ground-based ULF wave power during geomagnetic storms.ULF wave power can exhibit spatial coherence across ranges of L but does not rise and fall everywhere simultaneously due to the complicated evolution of cold plasma density and magnetic field strength in the inner magnetosphere.They also present evidence that the temporal variability of ULF wave power may vary with L. It may also be that spatial coherence varies with time and geomagnetic activity.The spatial variability (in the radial direction) of drift-averaged diffusion due to ULF waves throughout the outer radiation belt promises a rich vein of future work.
Sensitivity of radial diffusion to the variance of the full probabilistic distribution of D LL was explored in Experiment 3.For small variances, the diffusion results approach those of the deterministic model, as expected.But as the variance is increased, the diffusion results rapidly diverge.These results suggest that it is worth seeking alternative parameterizations that focus on variance reduction in the construction of the diffusion model.Another way to reduce the variance in the parameterization may be to focus on the calculation of D LL itself.For example, D E LL in the Ozeke et al. (2014) model was constructed via a mapping technique that utilized several assumptions: constant (low) wave number m = 1, constant width of the wave activity in latitude, and constant ionospheric conductance parameters (Ozeke et al., 2009).These quantities are typically not constant and contribute to the uncertainty in the deterministic model and should be included in the stochastic parameterization.The theoretical background from which D LL is based may also produce uncertainties.Several analytical diffusion rates based on magnetic and electric field assumptions exist, with L dependence ranging from L 6 -L 11 and frequency dependence on a range of wave modes (e.g., Birmingham, 1969;Cornwall, 1968;Elkington et al., 2003;Fälthammar, 1966Fälthammar, , 1968;;Fei et al., 2006;Schulz & Lanzerotti, 1974).If enough of the underlying variability in the deterministic model is known, the better the variability in the stochastic models can be characterized or accounted for.It should be mentioned however that natural variability might exist, which cannot be parameterized by any means.Deducing levels of natural variability in ULF wave-driven radial diffusion is necessary in understanding information always lost by a deterministic model.If these levels are substantial, our results suggest that a stochastic approach to modeling radial diffusion may be more robust.
The response of radial diffusion to higher likelihoods of an enhanced D LL , which dominates temporal evolutions, was explored in Experiment 4. It is evident that significantly more radial diffusion occurs for heavier tailed variabilities, indicating that the amount of diffusion is controlled by the relative importance of the large values of D LL in the distribution.A global upper bound for possible ULF wave power is justified since it is counterintuitive for ULF waves to have infinitely large power in a finite-sized magnetosphere.The shape of the distribution is therefore important.It may also be that the shape of the distribution of D LL is not constant.During quiet times when the outer radiation belt is relatively quiescent, the variability might be better represented heavily skewed to the left with a single small upper bound on ULF wave power.In a storm-time model where ULF wave activity is enhanced during the main and recovery phase (Murphy et al., 2011;Murphy et al., 2015;Rae et al., 2011), a right skewed ULF wave power distribution that favors larger ULF wave powers might be more suitable.Further research into tail values of the distribution of ULF wave power is important to constrain the physical upper bound of power variability to include in stochastic models.
In each of our experiments, ensemble averages and KDEs were compared to a (Ozeke et al., 2014) constant deterministic solution, which is based on the median of statistical ULF wave power.However, it may be more fair to compare the evolution of our numerical ensembles with an experiment where D LL is kept constant, but at the mean value of the distribution, especially since the ethos of constructing a diffusion coefficient is to consider the average behavior of the waves.Figure 7 indicates the results of a number of numerical experiments with constant D LL (mean, solid pink; upper quartile, dashed pink; and lower quartile, dash-dot pink) compared with the ensemble result using a LN distribution with Δt = 1 hr.We observe that the mean-based D LL only causes slightly more diffusion than the median based and is also significantly less diffusive than the ensemble averages.While inclusion of the LQ-and UQ-based D LL does result in a broad span of possible PSD solutions, the UQ produces diffusion only as strong as the ensemble averages, falling short of the regions of highest density seen in the ensemble solutions.It is apparent that having a deterministic representation of D LL fails to represent the underlying distribution of radial diffusion solutions found from the stochastic D LL time series, which better represent the true underlying distribution of ULF wave power.Our ensemble modeling highlights where efforts should be placed to get a better description of D LL , so that we can aim for a parameterization with a quantified uncertainty that truly represents the underlying distribution of possible solutions of the radial diffusion equation.
Diffusion due to other types of wave-particle interactions is important in the outer radiation belt, and similar modeling strategies may be required.Diffusion in pitch angle and energy due to higher-frequency waves is also highly variable (Watt et al., 2019), potentially with different time and length scales depending on location in the magnetosphere.It will be necessary to repeat similar numerical experiments to determine the stochastic parameters necessary to use in stochastic parameterizations of pitch angle and energy diffusion and then design observational analyses that can best constrain those parameters.

Conclusions
Our idealized experiments highlight the spatiotemporal impacts of including stochastic parameterizations in the ULF wave-driven radial diffusion.We have shown that diffusion is increased above the deterministic model when the diffusion coefficients vary more rapidly, when the spatial correlation of the diffusion across L-shells ranges from fully coherent to completely independent, and when the variance of the distribution is increased, or a more heavy-tailed distribution is used.We have demonstrated that future research should focus on the temporal evolution of ULF wave power, the spatial correlations of diffusion across L-shells, and the underlying distribution and variance of the radial diffusion coefficients.The successful implementation of a stochastic radial diffusion model requires variability parameters that are derived appropriately; that is, spatial and temporal scales of the variability may themselves vary in time and space.Our research motivates further investigation of stochastic methods for use in radiation belt diffusion models as a method to include the variability of wave-particle interactions in the inner magnetosphere.

Figure 1 .
Figure 1.Example ensemble member D LL time series shown for a range of temporal variability scales.In each case, the constant (Ozeke et al., 2014) deterministic D LL is multiplied by a log-normal variability at the relevant hour of variability, constrained by the empirical model and ULF wave power observations, and persists until to the next hour of variability where the process is repeated.Examples are shown for variability temporal scales of 1, 3, 6, 12, and 24 hr, along with the constant D LL with no variability.D LL shown here has units s −1 in line with the 1 s time step used in our numerical scheme.
1:34896 Þ are the parameters of the normally distributed logðϵÞ.Note that for a normally distributed random variable, the IQR is approximately 1.34869 multiplied by the standard deviation.We consider variability Δt = 1, 3, 6, 12, and 24 hr, and example ensemble members for each of these cases are shown in Figure1.They are effectively artificial representations of what might be observed in situ.

Figure 2 .
Figure 2. Example ensemble member D LL time series shown for a range of spatial variability scales.In each case, every 3 hr the constant (Ozeke et al., 2014) deterministic D LL is multiplied by log-normal variabilities on a variety of local spatial variability scales, constrained by the empirical model and ULF wave power observations, and persists for 3 hr where the process is then repeated.Examples are shown for variability spatial scales of 1L, 0.5L, and 0.1L, along with the global variability case and constant D LL with no variability.D LL shown here have units s −1 in line with the 1 s time step used in our numerical scheme.

4. 3 .
Experiment 3: Width of the D LL Probability Distribution The empirical (Ozeke et al., 2014) D LL parameterization is based on the median of statistical ULF wave power, and uncertainty in the parameterization has the multiplicative IQR 1 3 D L L ; 3D L L mentioned previously.We compare the IQR suggested by Ozeke et al. (2014) with larger and smaller IQRs, namely,

Figure 3 .
Figure 3. Ensemble results for the final PSD at the end of Experiment 1 for a range of temporal variability scales (1, 3, 6, 12, and 24 hr, respectively).The median (dashed), mean (dash-dot) ensemble profiles are shown, as well as the initial PSD profile (dotted) and the deterministic solution with constant deterministic D LL (solid).Ensemble kernel density estimates of the resulting electron PSD are also shown.

Figure 4 .
Figure 4. Ensemble results for the final PSD at the end of Experiment 2 for a range of spatial variability scales (global, 1L, 0.5L, and 0.1L, respectively).The description of lines and KDEs are as in Figure 3.

Figure 5 .
Figure 5. Ensemble results for the final PSD at the end of Experiment 3 for a range of log-normal variability IQRs (±2, ±3, ±6, and ±10 of the deterministic D LL , respectively).The description of lines and KDEs are as in Figure 3.

Figure 6 .
Figure 6.Ensemble results for the final PSD at the end of Experiment 4 for a range of variability probability distributions (Log-Normal, Log-Laplace, Log-Uniform, and Log-cauchy, respectively).The description of lines and KDEs are as in Figure 3.

Figure 7 .
Figure 7. PSD resulting from the radial diffusion equation after 2 days with constant Kp = 3, shown for a constant deterministic D LL based on the mean pink), LQ (dash-dot pink) and UQ (dash pink) of ULF wave power.These plots are laid over the first subplot in Figure 3.