Bayesian inference of quasi-linear radial diffusion parameters using Van Allen Probes
Rakesh Sarma, Mandar Chandorkar, Irina Zhelavskaya, Yuri Shprits, Alexander Drozdov, Enrico Camporeale
mmanuscript submitted to
JGR: Space Physics
Bayesian inference of quasi-linear radial diffusionparameters using Van Allen Probes
Rakesh Sarma , Mandar Chandorkar , Irina Zhelavskaya , , Yuri Shprits , ,Alexander Drozdov , Enrico Camporeale , , Centrum Wiskunde & Informatica, Amsterdam, The Netherlands GFZ German Research Centre For Geosciences, Potsdam, Germany Institute of Physics and Astronomy, University of Potsdam, Germany University of California, Los Angeles, CA, USA CIRES, University of Colorado, Boulder, CO,USA NOAA, Space Weather Prediction Center, Boulder, CO, USA
Key Points: • We present the first application of Bayesian parameter estimation to the problem ofquasi-linear diffusion in the radiation belts • The Bayesian approach allows the problem to be cast in probabilistic terms and forensemble simulations to be run • An improved accuracy is demonstrated when compared against standard deterministicmodels
Corresponding author: Rakesh Sarma, [email protected] –1– a r X i v : . [ phy s i c s . s p ace - ph ] F e b anuscript submitted to JGR: Space Physics
Abstract
The Van Allen radiation belts in the magnetosphere have been extensively studied usingmodels based on radial diffusion theory, which is based on a quasi-linear approach with pre-scribed inner and outer boundary conditions. The 1-d diffusion model requires the knowledgeof a diffusion coefficient and an electron loss timescale, which are typically parameterizedin terms of various quantities such as the spatial ( L ) coordinate or a geomagnetic index(for example, Kp ). These terms are empirically derived, not directly measurable, and henceare not known precisely, due to the inherent non-linearity of the process and the variableboundary conditions. In this work, we demonstrate a probabilistic approach by inferring thevalues of the diffusion and loss term parameters, along with their uncertainty, in a Bayesianframework, where identification is obtained using the Van Allen Probe measurements. Ourresults show that the probabilistic approach statistically improves the performance of themodel, compared to the parameterization employed in the literature. The Van Allen radiation belts consist of energetic electrons and protons originating fromthe solar wind, which are trapped by the Earth’s magnetic field. The dynamics of theseparticles are affected by an interplay of different mechanisms, including local accelerationand losses due to wave-particle interactions and external injections (Shprits, Elkington, etal., 2008; Shprits, Subbotin, et al., 2008; G. Reeves et al., 2013; Ukhorskiy & Sitnov, 2012).Understanding and forecasting the radiation belt particles is an essential part of SpaceWeather, since these particles can interfere with the satellites orbiting around the Earthand cause data loss. Hence, the estimation of the flux distribution in the space-time domainis a long-standing challenging problem in the space physics community (G. D. Reeves et al.,2003; Shprits, Elkington, et al., 2008; Shprits, Subbotin, et al., 2008).The standard way of studying radiation belt particles dynamics is through a quasi-linear approach, that dates back to the seminal paper of Kennel and Engelmann (1966) andhas been routinely applied to magnetospheric particle since the 1970’s (Schulz & Eviatar,1969; Lanzerotti et al., 1970; Schulz & Lanzerotti, 2012).The quasi-linear approach studies the particles’ dynamics based on their quasi-periodicorbits with respect to the field lines of the Earth’s magnetic field. Specifically, the cyclotronmotion (or gyration) is the circular motion of particles around the field line, the bouncemotion is due to the mirror force and along the field lines from one hemisphere to the other,and the drift motion is the orbit in the east or westward direction. These orbits are associatedto the conservation of adiabatic invariants, respectively named the first ( µ ), second ( K ) andthird ( L ∗ ) invariant (Roederer, 1970). The evolution of the particle density is then studiedas a diffusive process in the adiabatic invariant space, where the effect of resonant wave-particle interactions is described through diffusion coefficients. The third adiabatic invariantis associated with the longest timescale, and it can be violated by wave-particle interactionswith Ultra-Low Frequency (ULF) waves. Therefore, as a first approximation, the systemcan be studied as radially diffusive, assuming the other two invariants to be conserved. Inthis study, we estimate the phase space density (PSD) based on the case of pure radialdiffusion, which has been extensively investigated in the literature.The radial diffusion formulation involves a parametric representation of the diffusioncoefficient and the electron loss timescale, which are generally formulated as varying in L (the spatial coordinate) and Kp (a geomagnetic index used as a proxy for the amplitude ofgeomagnetic perturbations).The most extensively employed diffusion rate parameterization is obtained by Brautigamand Albert (2000) based on the October, 1990 storm. There have been other developmentswhere the radial diffusion coefficient is evaluated using various approaches. Fei et al. (2006)constructed time dependent diffusion terms from MHD simulations of September 1998 storm –2–anuscript submitted to JGR: Space Physics using ULF electric and magnetic field power spectral density. In Ozeke et al. (2012, 2014),the diffusion coefficient is expressed as a sum of two terms due to azimuthal electric andcompressional magnetic fields and expressions for these two terms are derived. A data-drivenapproach is employed to determine the radial diffusion coefficients in Su et al. (2015). Theimportance of the electron loss timescale in the radial diffusion modeling during storm-timeis demonstrated in Shprits and Thorne (2004), where the interplay between inward radialdiffusion and loss terms on the PSD is investigated. Summers et al. (2007) obtained theloss timescale due to the combined effects of chorus, plasmaspheric hiss, and EMIC waves.Shprits et al. (2007) has parameterized the loss timescale due to chorus waves as a functionof the geomagnetic index, AE . Further extension to include the effect of multiple storms isinvestigated in Tu et al. (2009) with an internal heating term. In Ali et al. (2016), the mag-netic and electric field measurements from the Van Allen Probes are used to compute powerspectral densities of both components, which are used with the Fei et al. (2006) formulationto obtain magnetic and electric components of the radial diffusion components.It is important to emphasize that even in three-dimensional studies of the radiationbelt that solve the diffusion equation in the whole adiabatic invariant space, hence takinginto account energy and pitch-angle scattering, the radial diffusion coefficient is still oftenparameterized using the Brautigam and Albert formula (see, e.g., Subbotin and Shprits(2009); Su et al. (2010); Tu et al. (2013); Bourdarie and Maget (2012); Welling et al.(2012))However, a single parameterization might not generalise well for different geomagneticconditions. Moreover, the diffusion problem is significantly influenced by the boundaryconditions, hence a deterministic parameterization might not be adequate to represent theuncertainties due to variable particle injections at the boundary. The aim of the currentinvestigation is to obtain a probabilistic representation of the PSD by introducing uncer-tainties in the parametric representation of coefficients in the underlying partial differentialequation. This is one of the first applications of the Bayesian framework approach to pa-rameter estimation of the radial diffusion equation. In the literature, data assimilation withan Extended Kalman filter technique is employed for determining the lifetime of electronsby Kondrashov et al. (2007).In this study, we perform Bayesian parameter identification approach for all the termsdefining the 1-d diffusion equation. Once the diffusion coefficient and the electron timeloss are defined as probability density functions, one can run an ensemble of simulations bysampling different values of the parameters and hence being able to estimate the uncertaintyof the output PSD (Camporeale et al., 2016; Camporeale, 2019). We use a data-drivenrepresentation of the input parameters, by employing Van Allen Probes data to identify theparameter distribution in a Bayesian setting (Spence et al., 2013).The manuscript is structured in the following parts: Section 2 provides an introductionto the radial diffusion model and the chosen parameterization of coefficients. In Section3, the framework for uncertainty propagation and the Bayesian identification is discussed.The results, discussion and comparison to the Van Allen measurements and other parame-terization in literature are discussed in Section 4. Finally, the findings and future researchdirections are discussed in Section 5. Several radiation belt models (Li, 2004; Subbotin & Shprits, 2009; G. D. Reeves etal., 2012; Albert et al., 2009) have been developed to quantify the radial transport. Thetrapped particles are quantified for given adiabatic invariants ( µ, K, L ) at time t with thePSD, f ( µ, K, L, t ). Under the assumption that the invariants ( µ, K ) are conserved, the radialdiffusion model is a one-dimensional model based on the modified Fokker-Planck equation –3–anuscript submitted to JGR: Space Physics (Walt, 1970), and is given by: ∂f∂t = L ∂∂L (cid:20) D LL L − ∂f∂L (cid:21) − fτ , (1)where, D LL is the diffusion coefficient and τ is the loss timescale, which is essentially a cor-rection term for unaccounted dynamics (such as pitch-angle and energy scattering). Variousparameterizations for D LL and τ are available in literature. In particular, in this study, werefer to the seminal work on the statistical diffusion rate of D LL proposed by Brautigamand Albert (2000), which models the geomagnetic storm of October 9 1990, given by: D LL ( L, t ) = α D L β D b D Kp ( t ) , (2)where they use the values α D = 4 . × − , β D = 10 . b D = 0 . electromagnetic part. Use of additional electrostatic part leads to incorrect simulation resultsas was shown by Kim et al. (2011). The choice of the L and Kp -dependent parameterizationfor the diffusion coefficient is motivated by the practise in literature. In the later part ofthis manuscript, we show that this choice leads to close agreement of PSD with respect toVan Allen Probes measurements. It has also been reported in earlier investigations (Aliet al., 2016) that D LL exhibited a weak energy dependence in the range of µ between 500MeV/G and 5000 MeV/G. Moreover, in this study we seek to demonstrate the applicationof Bayesian calibration to the diffusion problem, which can be applied to other choices of D LL parameterization in future investigations.For the electron lifetime τ , we employ Kp and L − dependent parameterization based onGu et al. (2012) (without energy dependence) inside the plasmasphere, while a model basedon Ozeke et al. (2014) is used outside the plasmasphere. The plasmapause position L pp is estimated using a recently developed Plasma density in the Inner magnetosphere Neuralnetwork-based Empirical (PINE) model (Zhelavskaya et al., 2017). The PINE density modelis developed using neural networks and is trained on the electron density data set from theVan Allen Probes Electric and Magnetic Field Instrument Suite and Integrated Science(EMFISIS) (Kletzing et al., 2013). The model reconstructs the plasmasphere dynamicswell (with a cross-correlation of 0.95 on the test set), and its global reconstructions ofplasma density are in good agreement with the IMAGE EUV images of global distributionof He+. The MLT-averaged plasmapause position is calculated using the output of thePINE model by applying a density threshold of 40 cm − to separate the plasmasphere fromthe outside of the plasmasphere. The Kp index is obtained from the OMNIWeb database.The parameterization for τ that we employ in (1) is given by: τ ( L, t ) = ( α τ + β τ L + b τ L ) /Kp ( t ) for L ≤ L pp = c τ /Kp ( t ) for L > L pp . (3)We choose a 1-year period from October 2012 to September 2013 for the purpose of analysis.The initial and outer boundary conditions are interpolated from the Van Allen Probes data.As mentioned, (1) is obtained for a constant value of ( µ, K ). In this study, we compare allthe results with Van Allen measurements for µ = 700 MeV/G and K = 0 . . · Re.In order to estimate the accuracy of our predictions, we use the relative error as a metric,given by: (cid:15) = | f va − f | f va , (4)where f va is the PSD value obtained from the Van Allen Probes and f is the PSD-estimatefrom (1). The absolute value of the discrepancy is used here, since we are interested inestimating the overall performance of the solver throughout the domain, which will beobtained by integrating this error across the domain. Now, as defined in Section 1, we are –4–anuscript submitted to JGR: Space Physics interested in obtaining an informed estimate on the parameters defining the coefficients D LL and τ in (1). In the next section, we introduce the Bayesian framework which is used toidentify these parameters. The Bayesian approach to the calibration of computer models was introduced in Kennedyand O’Hagan (2001), where the term calibration refers to adjusting the free parameters inorder for the model output to fit the observations. In our case, the forward model is rep-resented by (1), that is solved numerically with a standard finite-difference scheme on auniform grid in (
L, t ). Because solving this equation numerically is relatively fast, we optfor a standard Markov-Chain Monte Carlo (MCMC, Brooks et al. (2011)) procedure toexplore the space of unknown free parameters in (2) and (3), that are collected in a multi-dimensional vector Λ defined as: Λ := ( α D , β D , b D , α τ , β τ , b τ , c τ ) . In other words, Λ is the set of uncertain parameters that are to be identified from thisinvestigation. The ground truth for the PSD is taken from Van Allen Probes measurements,and we derive a Bayesian model trained over a data-set of 30 days from 01-Oct-2012 to 30-Oct-2012. The objective is to demonstrate that the method is generalised for time-periodson which it is not trained, hence the model will be tested for the rest of the year (November2012 to September 2013). The vector of parameters Λ is treated as a random variable,meaning that it is associated to an (unknown) probability density. The scope of the Bayesianinference is to estimate the probability of Λ , given the PSD observations, that we denotewith f + . Hence, we can use the classical Bayes’ rule (Gelman et al., 2004): P ( Λ | f + ) ∝ P ( f + | Λ ) · P ( Λ ) . (5)where, P ( f + | Λ ) is the likelihood that defines the discrepancy between the model estimateand the Van Allen Probes measurements, for a given realization of Λ . The term P ( Λ ) is the prior distribution of Λ , which encodes all prior physical information one might have aboutthe parameters. P ( Λ | f + ) is called the posterior distribution, which in general cannotbe expressed in closed form, but can be sampled through a Monte Carlo procedure (seebelow). Finally, once a sufficient number of samples from the posterior distribution hasbeen collected, an ensemble simulation can be obtained by propagating each realization of Λ through (1) to obtain the posterior predictive distribution of the PSD. As a first step, the priors on the parameter set Λ are defined. In the present investigation,we assume uniform priors, given by: Λ ∼ U ( Λ min , Λ max ) (6)where, Λ min and Λ max are chosen such that a wide domain is defined. The bound foreach parameter of Λ is shown in Table 1. Existing parameter values identified in liter-ature (Brautigam & Albert, 2000; Ozeke et al., 2014) are used as a reference, such thatparameterizations obtained in these investigations belong to the set of the defined bounds.We employ a Gaussian likelihood, by introducing an additive random variable η , speci-fied as an unbiased normal distribution with a covariance matrix Σ. Under the assumptionthat there are no modeling errors, the statistical model can then be written as: f + = f ( Λ ) + ηη ∼ N (0 , Σ) , (7) –5–anuscript submitted to JGR: Space Physics
Table 1.
Upper and lower bounds of prior defined in (6) for seven parameters in Λ Parameter Lower bound Upper bound α D . . × − β D . . b D . . α τ /Kp . /Kp . /Kpβ τ . . b τ . . c τ . . f ( Λ ) is the output of the forward model (1), which gives the likelihood as: P ( f + | Λ ) := P η ( f + − f ( Λ )) = | Σ | − n/ exp (cid:20) − (cid:0) f + − f ( Λ ) (cid:1) T Σ − (cid:0) f + − f ( Λ ) (cid:1)(cid:21) , (8)where | Σ | is the determinant of the covariance matrix and n is the number of Van Allenprobe measurements.For sampling the posterior, we use the Metropolis-Hastings Markov Chain Monte Carlo(MCMC) algorithm (Hastings, 1970), with a reversible Markov random-walk. The algorithmperforms a random walk by sampling from a proposal distribution with a prescribed variance.Each new sample distribution is conditional only on the current sample, hence the generatedsequence of samples resembles a Markov chain. Once a sample is chosen, a ratio of thedensities for the two consecutive samples is computed, which informs the iterator if thejump results in increase/decrease of the posterior density. If the jump results in a probabilityhigher than an acceptance limit, the sample is accepted with that probability or it is rejected.The proposals are selected from a normal distribution, with tuned variance and a total of20000 samples, where initial 1000 samples discarded as burn-in period, and a thinning factorof 2 (Gelman et al., 2004) is applied for the remaining samples to improve independence. It isto be noted here that choice for the burn-in period and thinning factor is problem-specific andit is decided based on the convergence of the parameters. The standard checks of convergenceand auto-correlation for MCMC were applied (Hastings, 1970; Chib & Greenberg, 1995). In this section, we show the application of the Bayesian inference for identifying theparameters in (1), which are updated with the Van Allen Probes measurements. The 7-parameter uncertain space given by Λ is propagated with the forward solver and then theBayesian framework is applied. Figure 1 shows the posterior distributions of the identified parameters Λ of D LL and τ . As mentioned, uniform priors are specified with wide domains, where for example, theprior of β D ∼ U (1 , c τ , whose uncertainty is still high after identification,which shows that the data is not informative enough to infer c τ . It is also interesting tonote that c τ in (3) is employed above the plasmapause location L pp , which implies that thisterm is mostly employed at higher values of L , where the flux from the boundary dominatesthe density predictions. As a result, the inferred values of c τ strongly depend on the influx ofparticles at the outer boundary, which is highly uncertain and a deterministic representationis not possible. –6–anuscript submitted to JGR: Space Physics
Furthermore, parameters β τ and b τ have values close to zero, which are coefficients for L and L -dependence terms of τ in (3) respectively. This shows the variation of τ withrespect to L is minimal. It should also be mentioned here that other investigations of higherorder variation of τ with L also showed negligible dependence.In Figure 2, 2-d marginal distributions of each of the parameter-pairs for five of theidentified parameters in Λ are plotted along with the 1-d probability densities. The param-eters β τ and b τ are not plotted here since their identified samples are close to zero. It isinteresting to observe that clear correlations are identified between many parameter-pairs.Parameters α D and β D have a clear negative correlation with an asymptotic trend. Further-more, parameter-pairs α D and α τ , β D and c τ , α τ and c τ have negative correlation, while α D and c τ seem to have a positive correlation. Other parameter pairs do not seem to showa clearly visible correlation behaviour.For the identified posteriors, we determine the maximum a posteriori (MAP) estimateof the probability distribution, which is the parameter-value corresponding to the point ofmaximum probability. Table 2 compares the MAP estimate of the posteriors to the parame-terization employed in the 1-d model of Drozdov et al. (2017), where two parameterizationsfor the 1-d diffusion problem are presented. In this study, we compare our estimates to the1-d model based on Brautigam and Albert (2000) diffusion parameterization in Drozdov etal. (2017). Analytical Kp -dependent model based on Shprits et al. (2005) was employedoutside the plasmasphere, where τ = 3 /Kp days is used. In order to take into account theeffect of the variable outer boundary, Shprits et al. (2006) used τ = 5 /Kp as some of the losswas provided by the outward diffusion and longer lifetimes resulted in a better comparisonwith observations. Since we have a variable boundary in this investigation, we compareour estimates to a parameterization with a longer lifetime similar to (Drozdov et al., 2017;Ozeke et al., 2014), given by τ = 6 /Kp days. It is to be noted that for the subsequentsections of this manuscript, we use the same choice of parameterization as reference.Figure 3 compares the probabilistic D LL obtained with the identified posteriors shownin Figure 1 to the D LL employed in Brautigam and Albert (2000). At low values of Kp , thetwo predictions are close to each other, while the discrepancy between them progressivelyincreases as Kp is increased (note the vertical logarithmic scale). With Brautigam andAlbert (2000) parameterization, D LL has higher values at higher L with increase in Kp compared to the probabilistic treatment. This could be related to the Kp -dependent τ parameterization, both above and below the plasmapause location.It is observed that the MAP estimates for the parameters of D LL namely α D , β D and b D , have values close to the reference parameterizations. For β τ and b τ , reference valuesdo not exist and they both are close to zero. In (3), α τ is Kp -dependent, while in thereference, it has a constant value of 10 . . /Kp above the plasmapause location.With the current choice of parameterization, this term would be significantly smaller forhigh Kp while being comparable to the reference value for low Kp . Furthermore, the term c τ also significantly varies with respect to the reference value. As discussed already, c τ hashigh variance, and the MAP estimate is higher than the reference. Hence it is observedthat in terms of the MAP comparison, the lifetime estimates are giving longer time scalesthan that estimated in literature. This is clearly due to omission of local acceleration inthe 1-d diffusion model. Similarly, the 1-d model in Drozdov et al. (2017) underestimatedthe observations. Hence, the Bayesian parameter estimation is trying to compensate for themissing physical processes. Here we investigate the effect of the predicted uncertainties of the diffusion parameterson the evolution of the PSD. In practice, because the posterior distributions are derived by ∼ –7–anuscript submitted to JGR: Space Physics
Table 2.
Comparison of parameterization based on the MAP estimate of the identified posteriorto 1-d model in Drozdov et al. (2017)
Parameter MAP estimate Reference α D . × − . × − β D .
43 10 . b D .
47 0 . α τ /Kp . /Kp . β τ .
01 - b τ . c τ .
22 6 . Probabilistic variation of density
In order to visualise the probabilistic field of density in the (
L, t ) domain, we plot thePSD at various time instances. Figure 4 shows the PSD variation at multiple time snapshots,advancing in time from Figure 4a to 4j. The solid line represents the mean of the posteriorpredictions, while the shaded area shows the 1 σ confidence interval. At the beginning of thesimulation near the initial condition (Figure 4a), it can be observed that the uncertainty islow. Also there is a good agreement with respect to the Van Allen Probe measurements (reddots), which is expected since the problem is initialised with the data. The low uncertaintyis also expected since the identified posteriors have low variance.The solution is further time-marched and the uncertainty progressively increases (Fig-ures 4b-4e). It can be observed that the uncertainty near the upper boundary is smaller,while it is higher around the lower boundary at all time instances. This is due to the in-fluence of the upper boundary condition, which is interpolated from the Van Allen Probedata. It is interesting to observe the influence of the injection of particles from the outerboundary. In Figure 4f, the uncertainty around L = 4 is observable, while after the influxof particles in Figure 4g, the uncertainty reduces significantly. Also the Van Allen Probedata has significant noise and discontinuities in Figure 4g, since higher density gradientsdiminishes the smoothness of interpolation. The uncertainty again starts to increase inFigure 4h and a similar effect with injection of particles can be observed in Figures 4i - 4j.The red dots indicate the interpolated value of the Van Allen Probes data at a given time.Due to the effect of interpolation, it is not expected that the simulation outputs and thered dots would coincide perfectly. Overall the model performs well in terms of agreementwith the Van Allen Probes measurements, often predicting within one standard deviation.Furthermore, the lack of agreement at lower values of L is accounted by the confidenceintervals. It is observed that the model discrepancy is higher when there is sudden influxof particles, which cannot be properly accounted for in a diffusive process. However, withour probabilistic approach, we are able to obtain overlap with respect to the data at mostvalues of L . Continuous Rank Probability Score (CRPS)
In order to measure the performance of the probability forecasts, we employ the Con-tinuous Rank Probability Score (CRPS), which is a metric to evaluate the quadratic dis-crepancy between the cumulative distribution function (CDF) F of the forecasts and theempirical CDF of the data or observations F d , which is given by: CRP S = (cid:90) + ∞−∞ (cid:0) F ( f ) − F d ( f ) (cid:1) df, (9) –8–anuscript submitted to JGR: Space Physics where F d is a Heaviside function for a scalar observation. Thus, CRPS is computed byintegrating the area between the two distribution functions. If the computed area is low,the forecasts are close to the observations. Also, CRPS is a convenient measure to comparethe performance of deterministic and probabilistic forecasts (Camporeale et al., 2019). Ifthe forecast is deterministic, CRPS reduces to a mean absolute error. Because we havea probability density of the PSD for any point of the domain ( L, t ), we can numericallycalculate the CRPS over the whole domain. Figure 5 compares the CRPS obtained fromthe probabilistic predictions of PSD with the identified posteriors shown in Figure 1 andthe deterministic predictions with the parameterization referenced in Table 2. With respectto assessing the quality of the predictions, values closer to zero show that the predictor isperforming better. A qualitative comparison shows that the probabilistic forecast performsbetter at most of the locations in the (
L, t ) domain. Figure 6 shows the comparison of thedensity of CRPS in the whole domain with the probabilistic and deterministic approaches.The mean of the CRPS is 0 .
429 and 0 . . . The advantage of the proposed method clearly stays in the ability of deriving a prob-abilistic estimate of the PSD and the consequent ability of identifying cases of large uncer-tainty. However nothing precludes to use the information on the posterior distribution ofthe parameters in a deterministic way, by running a simulation corresponding to the MAPestimate of the parameters. In this section, we compare the PSD at µ = 700 MeV/G and K = 0 . . · Re, predicted by the MAP estimate of the identified posteriors of Λ to theVan Allen Probe measurements. As a reference we again employ the 1-d model comparisonprovided in Drozdov et al. (2017) (see Table 2). Additionally, we compare the predictions toother parameterizations existing in literature. We define D LL based on parameterizationsobtained in Ozeke et al. (2014) and Ali et al. (2016). In Ozeke et al. (2014), D LL is definedas the sum of azimuthal electric field D ELL and compressional magnetic field D MLL , given by: D ELL = L . · − (0 . L +0 . Kp ) ,D MLL = L . · − ( − . L +0 . L − . Kp +0 . Kp ) . (10)In Ali et al. (2016), the drift-averaged power spectral densities are used to obtain the electricand magnetic components of the radial diffusion coefficient, and a genetic algorithm is usedto derive a simple model for each component in the least squares sense, which is given by: D ELL = exp ( − .
951 + 0 . · Kp · L + 1 . · L ) ,D MLL = exp ( − .
253 + 0 . · Kp · L + L ) . (11)Figure 7 shows the PSD in the form of scatter plots and compares five different realisations.It can be observed that the forward model estimations have good qualitative agreementwith the Van Allen measurements. The model predictions by (1) with the MAP-based andthe reference parameterizations are in close agreement. This is expected since the MAPestimate of the posteriors, in particular the diffusion parameters, match closely with thereference values, which can be seen in Table 2. The PSD estimates from the diffusion modelbased on (10) (Ozeke et al., 2014) are also close to the MAP-based model predictions. Forthe diffusion model based on (11) (Ali et al., 2016), the PSD predictions at higher valuesof L are in qualitative agreement, however there is higher discrepancy at lower L values. Ithas to be mentioned here that the model (11) was calibrated in the range 3 . ≤ L ≤ . ≤ Kp ≤
5. Hence in this case we are performing extrapolation, which may lead tolower performance of the model.We also compare the relative error given by (4) in estimation of PSD with scatterplots shown in Figure 8. As expected, predictions from the MAP-based and reference –9–anuscript submitted to
JGR: Space Physics parameterization match closely. The MAP-based model performs better at lower valuesof L . Similarly, the model from Ozeke et al. (2014) also exhibits similar error behaviour.As expected, the model from Ali et al. (2016) does not perform well at lower values of L .Overall, in terms of the average relative error (cid:15) avg = (cid:80) i (cid:15) i /m , where m is the number of VanAllen Probes measurements, it is found that (cid:15) avg = 0 . . .
078 and 0 .
326 with theMAP-based, reference parameterization, based on (Ozeke et al., 2014), and (Ali et al., 2016)respectively. It can thus be concluded that the MAP-based model performs better than allthe compared models and a small improvement is obtained if the MAP-based parametervalues are employed in a deterministic setting.
In this paper, we demonstrated a probabilistic framework for estimating the PSD ofenergetic particles in the magnetosphere, following a standard quasi-linear radial diffusionequation. The probabilistic treatment is based on the Bayesian framework, where, initially,uncertainties are assumed on the parameterization of the diffusion and loss terms, whichare iteratively updated using the Van Allen Probes measurements as ground truth. Finally,the informed posteriors are propagated to obtain a probabilistic estimate of particles’ phasespace density. It is shown that the data is able to reduce the uncertainties considerably, andthe probabilistic treatment improves the quality of the predictions. This approach providesan alternate approach to understanding the 1-d radial diffusion, since a deterministic formof parameterization may not be able to represent the non-linear behaviour of this boundary-dependent problem.In terms of the choice of the prior, Ali et al. (2016) obtained uncertainty bounds andattempted to define the probability distribution of the total diffusion coefficient. This couldbe useful to replace the uniform prior in future investigations, in order to define the domainand improve the convergence of the Bayesian model. The Bayesian estimator predictedlonger time scales for the lifetime parameter, to take into account the local acceleration,which is omitted in the 1-d model. In the future, a similar approach can be extended to3-D modeling that includes calculation of loss and local acceleration from physics. Thatwould help better quantify the relatively unknown radial diffusion rates. Inferior radialdiffusion rate is a very challenging task. In-situ measurements do not allow for inferringthe m-numbers of waves or observing global distribution of waves in MLT and L. Groundmeasurements allow to find more global maps but observations on the ground may bevery inaccurate during disturbed conditions, as ionosphere may shield the waves and theymay also be guided to another location, and it may be difficult to trace back the groundobservations into space. The Bayesian parameter estimation certainly has a great potentialto shed light on this important problem. –10–anuscript submitted to
JGR: Space Physics
Figure 1.
Posterior distributions of identified parameters Λ of D LL and τ .–11–anuscript submitted to JGR: Space Physics
Figure 2.
Identified 1-d posterior densities and 2-d marginal distributions for five parametersdefining the diffusion and loss terms in (1).
Figure 3.
Probabilistic representation of D LL (in log scale) with the identified posteriors, with Kp = 2, 4, 6 and 8 from left to right. The red dots represent D LL obtained with Brautigam andAlbert (2000) parameterization, while the blue solid lines represent the mean of the probabilistic D LL obtained from the posteriors. –12–anuscript submitted to JGR: Space Physics
Figure 4.
Snapshots of PSD at µ = 700 MeV/G and K = 0 . . · Re for different timeinstances. Top panel shows the scatter plot of Van Allen Probes measurements with red linesindicating the time instances at which the probabilistic response is plotted, advancing in time from(a) to (j). Each subfigure from (a) to (j) shows the posterior-propagated PSD with mean (solid blueline), 1 σ confidence intervals (blue shade) and the Van Allen measurements (red dots) correspondingto the time instance t to t in the top panel. –13–anuscript submitted to JGR: Space Physics
Figure 5.
Scatter plot of CRPS for the time period from October 2012 to September 2013 fromVan Allen Probes. The top panel correspond to the probabilistic predictions of PSD by (1) withthe identified posteriors, while the bottom panel is obtained from the deterministic predictions of(1) with the reference parameterization (Drozdov et al., 2017).
Figure 6.
Density plot of CRPS samples shown in Figure 5. The black solid line correspond tothe probabilistic predictions of PSD by (1), while the red dotted line correspond to predictions bythe reference parameterization (Drozdov et al., 2017).–14–anuscript submitted to
JGR: Space Physics
Figure 7.
Scatter plot of PSD (in logarithms of base 10) at µ = 700 MeV/G and K = 0 . . · Re in L − t domain for the time period from October 2012 to September 2013 obtained fromVan Allen Probes measurements (top panel), model predictions by (1) and parameterized by: theMAP-estimate of posteriors (second panel), reference (Drozdov et al., 2017) (third panel), (10)from Ozeke et al. (2014) (fourth panel), (11) from Ali et al. (2016) (fifth panel), and Kp index fromOMNIWeb database (bottom panel). Figure 8.
Scatter plot of relative error (cid:15) in L − t domain given by (4), in estimation of PSDfor the time period from October 2012 to September 2013 based on model prediction by (1) andparameterized by: the MAP-estimate of posteriors (top panel), reference (Drozdov et al., 2017)(second panel), (10) from Ozeke et al. (2014) (third panel), and (11) from Ali et al. (2016) (bottompanel). –15–anuscript submitted to JGR: Space Physics
Acknowledgments
The Kp dataset is obtained from the Omniweb database ( https://omniweb.gsfc.nasa.gov ). The data generated for this research is available online at: https://doi.org/10.4121/uuid:399b9a42-f07e-4b30-a4aa-86abfeeddac0 . References
Albert, J. M., Meredith, N. P., & Horne, R. B. (2009). Three-dimensional diffusion sim-ulation of outer radiation belt electrons during the 9 october 1990 magnetic storm.
Journal of Geophysical Research: Space Physics , (A9). Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2009JA014336 doi: 10.1029/2009JA014336Ali, A. F., Malaspina, D. M., Elkington, S. R., Jaynes, A. N., Chan, A. A., Wygant, J.,& Kletzing, C. A. (2016, 10). Electric and magnetic radial diffusion coefficientsusing the Van Allen probes data. Journal of Geophysical Research: Space Physics , (10), 9586–9607. Retrieved from https://doi.org/10.1002/2016JA023002 doi:10.1002/2016JA023002Bourdarie, S., & Maget, V. (2012). Electron radiation belt data assimilation with anensemble kalman filter relying on the salammbˆo code. In Annales geophysicae (Vol. 30,p. 929).Brautigam, D. H., & Albert, J. M. (2000, 1). Radial diffusion analysis of outer radiationbelt electrons during the october 9, 1990, magnetic storm.
Journal of GeophysicalResearch: Space Physics , (A1), 291–309. Retrieved from https://doi.org/10.1029/1999JA900344 doi: 10.1029/1999JA900344Brooks, S., Gelman, A., Jones, G., & Meng, X.-L. (2011). Handbook of markov chain montecarlo . CRC press.Camporeale, E. (2019). The challenge of machine learning in space weather: Nowcastingand forecasting.
Space Weather , (8), 1166–1207.Camporeale, E., Chu, X., Agapitov, O., & Bortnik, J. (2019). On the generation of proba-bilistic forecasts from deterministic models. Space Weather , (3), 455–475.Camporeale, E., Shprits, Y., Chandorkar, M., Drozdov, A., & Wing, S. (2016). On thepropagation of uncertainties in radiation belt simulations. Space Weather , (11),982–992.Chib, S., & Greenberg, E. (1995). Understanding the metropolis-hastings algorithm. The American Statistician , (4), 327–335. Retrieved from Drozdov, A. Y., Shprits, Y. Y., Aseev, N. A., Kellerman, A. C., & Reeves, G. D. (2017, 1).Dependence of radiation belt simulations to assumed radial diffusion rates tested fortwo empirical models of radial transport.
Space Weather , (1), 150–162. Retrievedfrom https://doi.org/10.1002/2016SW001426 doi: 10.1002/2016SW001426Fei, Y., Chan, A. A., Elkington, S. R., & Wiltberger, M. J. (2006, 12). Radial diffusionand mhd particle simulations of relativistic electron transport by ulf waves in theseptember 1998 storm. Journal of Geophysical Research: Space Physics (19782012) , (A12). Retrieved from https://doi.org/10.1029/2005JA011211 doi: 10.1029/2005JA011211Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. (2004). Bayesian data analysis (2nded.). Chapman and Hall/CRC.Gu, X., Shprits, Y. Y., & Ni, B. (2012). Parameterized lifetime of radiation belt electrons in-teracting with lower-band and upper-band oblique chorus waves.
Geophysical ResearchLetters , (15). Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2012GL052519 doi: 10.1029/2012GL052519Hastings, W. K. (1970). Monte carlo sampling methods using markov chains and theirapplications. Biometrika , (1), 97–109. Retrieved from –16–anuscript submitted to JGR: Space Physics
Kennedy, M. C., & O’Hagan, A. (2001). Bayesian calibration of computer models.
Journalof the Royal Statistical Society: Series B (Statistical Methodology) , (3), 425–464.Kennel, C., & Engelmann, F. (1966). Velocity space diffusion from weak plasma turbulencein a magnetic field. The Physics of Fluids , (12), 2377–2388.Kim, K., Shprits, Y., Subbotin, D., & Ni, B. (2011, 10). Understanding the dynamic evo-lution of the relativistic electron slot region including radial and pitch angle diffusion. Journal of Geophysical Research: Space Physics (19782012) , (A10). Retrievedfrom https://doi.org/10.1029/2011JA016684 doi: 10.1029/2011JA016684Kletzing, C. A., Kurth, W. S., Acuna, M., MacDowall, R. J., Torbert, R. B., Averkamp,T., . . . Tyler, J. (2013, Nov 01). The electric and magnetic field instrumentsuite and integrated science (emfisis) on rbsp. Space Science Reviews , (1), 127–181. Retrieved from https://doi.org/10.1007/s11214-013-9993-6 doi: 10.1007/s11214-013-9993-6Kondrashov, D., Shprits, Y., Ghil, M., & Thorne, R. (2007, 10). A kalman filter tech-nique to estimate relativistic electron lifetimes in the outer radiation belt. Jour-nal of Geophysical Research: Space Physics (19782012) , (A10). Retrieved from https://doi.org/10.1029/2007JA012583 doi: 10.1029/2007JA012583Lanzerotti, L., Maclennan, C., & Schulz, M. (1970). Radial diffusion of outer-zone electrons:An empirical approach to third-invariant violation. Journal of Geophysical Research , (28), 5351–5371.Li, X. (2004, March). Variations of 0.76.0 mev electrons at geosynchronous orbit as afunction of solar wind. Space Weather , (3), 1-10. doi: 10.1029/2003SW000017Ozeke, L. G., Mann, I. R., Murphy, K. R., Rae, I. J., & Milling, D. K. (2014, 3).Analytic expressions for ulf wave radiation belt radial diffusion coefficients. Jour-nal of Geophysical Research: Space Physics , (3), 1587–1605. Retrieved from https://doi.org/10.1002/2013JA019204 doi: 10.1002/2013JA019204Ozeke, L. G., Mann, I. R., Murphy, K. R., Rae, I. J., Milling, D. K., Elkington, S. R., . . .Singer, H. J. (2012, 4). Ulf wave derived radiation belt radial diffusion coefficients. Journal of Geophysical Research: Space Physics (19782012) , (A4). Retrieved from https://doi.org/10.1029/2011JA017463 doi: 10.1029/2011JA017463Reeves, G., Spence, H. E., Henderson, M., Morley, S., Friedel, R., Funsten, H., . . . others(2013). Electron acceleration in the heart of the van allen radiation belts. Science , (6149), 991–994.Reeves, G. D., Chen, Y., Cunningham, G. S., Friedel, R. W. H., Henderson, M. G.,Jordanova, V. K., . . . Zaharia, S. (2012). Dynamic radiation environment as-similation model: Dream. Space Weather , (3). Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2011SW000729 doi: 10.1029/2011SW000729Reeves, G. D., McAdams, K. L., Friedel, R. H. W., & O’Brien, T. P. (2003). Accelerationand loss of relativistic electrons during geomagnetic storms. Geophysical ResearchLetters , (10). Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2002GL016513 doi: 10.1029/2002GL016513Roederer, J. G. (1970). Dynamics of geomagnetically trapped radiation (Vol. 2). Springer-Verlag.Schulz, M., & Eviatar, A. (1969). Diffusion of equatorial particles in the outer radiationzone.
Journal of Geophysical Research , (9), 2182–2192.Schulz, M., & Lanzerotti, L. J. (2012). Particle diffusion in the radiation belts (Vol. 7).Springer Science & Business Media.Shprits, Y. Y., Elkington, S. R., Meredith, N. P., & Subbotin, D. A. (2008). Reviewof modeling of losses and sources of relativistic electrons in the outer radiation beltI: Radial transport.
Journal of Atmospheric and Solar-Terrestrial Physics , (14),1679 - 1693. Retrieved from (Dynamic Variability of Earth’s Radiation Belts) doi:https://doi.org/10.1016/j.jastp.2008.06.008Shprits, Y. Y., Meredith, N. P., & Thorne, R. M. (2007, 6). Parameterization of radia- –17–anuscript submitted to JGR: Space Physics tion belt electron loss timescales due to interactions with chorus waves.
GeophysicalResearch Letters , (11). Retrieved from https://doi.org/10.1029/2006GL029050 doi: 10.1029/2006GL029050Shprits, Y. Y., Subbotin, D. A., Meredith, N. P., & Elkington, S. R. (2008). Review of mod-eling of losses and sources of relativistic electrons in the outer radiation belt II: Localacceleration and loss. Journal of Atmospheric and Solar-Terrestrial Physics , (14),1694 - 1713. Retrieved from (Dynamic Variability of Earth’s Radiation Belts) doi:https://doi.org/10.1016/j.jastp.2008.06.014Shprits, Y. Y., & Thorne, R. M. (2004, 4). Time dependent radial diffusion model-ing of relativistic electrons with realistic loss rates. Geophysical Research Letters , (8). Retrieved from https://doi.org/10.1029/2004GL019591 doi: 10.1029/2004GL019591Shprits, Y. Y., Thorne, R. M., Friedel, R., Reeves, G. D., Fennell, J., Baker, D. N., &Kanekal, S. G. (2006, 11). Outward radial diffusion driven by losses at magnetopause. Journal of Geophysical Research: Space Physics (19782012) , (A11). Retrievedfrom https://doi.org/10.1029/2006JA011657 doi: 10.1029/2006JA011657Shprits, Y. Y., Thorne, R. M., Reeves, G. D., & Friedel, R. (2005, June). Radial diffusionmodeling with empirical lifetimes: comparison with CRRES observations. AnnalesGeophysicae , (4), 1467-1471. Retrieved from https://hal.archives-ouvertes.fr/hal-00317747 Spence, H. E., Reeves, G., Baker, D., Blake, J., Bolton, M., Bourdarie, S., . . . others (2013).Science goals and overview of the radiation belt storm probes (rbsp) energetic particle,composition, and thermal plasma (ect) suite on nasas van allen probes mission.
SpaceScience Reviews , (1-4), 311–336.Su, Z., Xiao, F., Zheng, H., & Wang, S. (2010). Steerb: A three-dimensional code forstorm-time evolution of electron radiation belt. Journal of Geophysical Research:Space Physics , (A9).Su, Z., Zhu, H., Xiao, F., Zong, Q. G., Zhou, X. Z., Zheng, H., . . . Wygant, J. R. (2015). Nature Communications , (1).Subbotin, D. A., & Shprits, Y. Y. (2009). Three-dimensional modeling of the radia-tion belts using the versatile electron radiation belt (verb) code. Space Weather , (10). Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2008SW000452 doi: 10.1029/2008SW000452Summers, D., Ni, B., & Meredith, N. P. (2007, 4). Timescales for radiation belt electronacceleration and loss due to resonant waveparticle interactions: 2. evaluation for vlfchorus, elf hiss, and electromagnetic ion cyclotron waves. Journal of GeophysicalResearch: Space Physics (19782012) , (A4). Retrieved from https://doi.org/10.1029/2006JA011993 doi: 10.1029/2006JA011993Tu, W., Cunningham, G., Chen, Y., Henderson, M., Camporeale, E., & Reeves, G. (2013).Modeling radiation belt electron dynamics during gem challenge intervals with thedream3d diffusion model. Journal of Geophysical Research: Space Physics , (10),6197–6211.Tu, W., Li, X., Chen, Y., Reeves, G. D., & Temerin, M. (2009, 2). Stormdependent radiationbelt electron dynamics. Journal of Geophysical Research: Space Physics (19782012) , (A2). Retrieved from https://doi.org/10.1029/2008JA013480 doi: 10.1029/2008JA013480Ukhorskiy, A., & Sitnov, M. (2012). Dynamics of radiation belt particles. In The van allenprobes mission (pp. 545–578). Springer.Walt, M. (1970). Radial diffusion of trapped particles. In B. M. McCormac (Ed.),
Particlesand fields in the magnetosphere (pp. 410–415). Dordrecht: Springer Netherlands.Welling, D., Koller, J., & Camporeale, E. (2012). Verification of spacepy’s radial diffusionradiation belt model.
Geoscientific Model Development , (2), 277.Zhelavskaya, I. S., Shprits, Y. Y., & Spasojevi, M. (2017, 11). Empirical modeling ofthe plasmasphere dynamics using neural networks. Journal of Geophysical Research: –18–anuscript submitted to
JGR: Space Physics
Space Physics , (11), 11,227–11,244. Retrieved from https://doi.org/10.1002/2017JA024406 doi: 10.1002/2017JA024406doi: 10.1002/2017JA024406