A variational approach to Data Assimilation in the Solar Wind
CConfidential manuscript submitted to
Space Weather
A variational approach to Data Assimilation in theSolar Wind
Matthew Lang , Mathew J. Owens Le Laboratoire des Sciences du Climat et de l’Environnement, CEA-CNRS-UVSQ, 91191 Gif Sur Yvette,France Department of Meteorology, University of Reading, Reading, Berkshire, UK
Key Points: • This paper develops a variational data assimilation (DA) method and applies itto a simple solar wind propagation model. • Such a DA scheme enables the inner boundary of the solar wind model to beupdated using observations in near-Earth space. • Experiments performed with both synthetic and STEREO observations showthat the DA method is able to reduce errors in the solar wind speeds.
Corresponding author: Matthew Lang, [email protected] –1– a r X i v : . [ phy s i c s . s p ace - ph ] O c t onfidential manuscript submitted to Space Weather
Abstract
Variational Data Assimilation (DA) has enabled huge improvements in the skill ofoperational weather forecasting. In this study, we use a simple solar-wind propagationmodel to develop the first solar-wind variational DA scheme. This scheme enablessolar-wind observations far from the Sun, such as at 1 AU, to update and improvethe inner boundary conditions of the solar wind model (at 30 solar radii). In thisway, observational information can be used to improve estimates of the near-Earthsolar wind, even when the observations are not directly downstream of the Earth. Us-ing controlled experiments with synthetic observations we demonstrate this method’spotential to improve solar wind forecasts, though the best results are achieved in con-junction with accurate initial estimates of the solar wind. The variational DA schemeis also applied to STEREO in-situ observations using initial solar wind conditions sup-plied by a coronal model of the observed photospheric magnetic field. We considerthe period Oct 2010-Oct 2011, when the STEREO spacecraft were approximately 80 ◦ ahead/behind Earth in its orbit. For 12 of 13 Carrington Rotations, assimilation ofSTEREO data improves the near-Earth solar wind estimate over the non-assimilatedstate, with a 18 .
4% reduction in the root-mean-squared-error. The largest gains aremade by the DA during times when the steady-state assumption of the coronal modelsbreaks down. While applying this pure variational approach to complex solar-windmodels is technically challenging, we discuss hybrid DA approaches which are simplerto implement and may retain many of the advantages demonstrated here.
In meteorology, data assimilation has long been used to improve initial condi-tions for forecasting, leading to a reduction in the ‘butterfly effect’ and hence im-provements in forecasting skill. Furthermore, improvements in the implementation ofdata assimilation methods into numerical weather prediction models have led to hugeimprovements in the forecasting accuracy of longer lead-times over the past 20 − Kalnay , 2003]. However, space weather forecasting has yet to exploit the greatpotential available from implementing data assimilation methods into their forecastingmodels.For space-weather forecasting, data assimilation has been attempted in threemain domains: the photosphere, the solar wind and the ionosphere. Ionospheric dataassimilation is arguably the most mature and as the number of observations of the iono-sphere increases, so does its importance and effectiveness. Various data assimilationmethods have been applied to the ionosphere, such as 3DVar [
Bust and Mitchell , 2008],4DVar [
Wang et al. , 2004] and Local Ensemble Transform Kalman Filter (LETKF)[
Durazo et al. , 2017]. Photospheric data assimilation, such as the Air Force DataAssimilative Photospheric Flux Transport (ADAPT) model [
Arge et al. , 2010], usesobservations of the magnetic field at the Sun’s surface with physics-based temporalevolution to improve the inner boundary condition to coronal models. The improvedrepresentation of the corona can in turn be used to generate improved inner-boundaryconditions for solar wind models.In this study, we are looking to exploit observations of the solar wind itself to fur-ther improve the inner-boundary conditions for solar wind models (and, as a by prod-uct, provide an observationally constrained validation dataset for the outer boundaryof coronal models). Previous studies applying data assimilation to in-situ observationsof the solar wind, such as
Lang et al. [2017], have focussed on using ensemble-baseddata assimilation methods to improve the representation of the modelled solar wind.The advantage of these methods is that they are relatively easy to implement forcomplex numerical models and recent developments have allowed these methods to beincorporated via Message Passing Interface, MPI [
Browne and Wilson , 2015;
Nerger –2–onfidential manuscript submitted to
Space Weather
Figure 1: A schematic of a DA scheme that updates the model inner boundary (thewhite circle) on the basis of observations from a position behind Earth in its orbit (theyellow star). This enables the updated model conditions (the purple regions) to persistuntil solar rotation brings them to the forecast point at Earth’s location (the black cir-cle). Without this ability, any localised change to the model state is quickly lost from themodel domain by the super-sonic radial solar wind flow. et al. , 2005], in parallel to the numerical models. However, as shown in
Lang et al. [2017], in order to improve the forecast of solar wind conditions in near-Earth space,these methods require observations downstream of the Earth, which are not routinelyavailable. This is due to the localisation of the data assimilation, which means anyimprovements to the model state from the observations are swept out of the modeldomain due to the continual radial outflow of the solar wind. In order to make apersistent change to the model state, the model inner boundary conditions must beupdated. This is shown schematically in Figure 1. Ensemble-based Kalman filtertechniques cannot propagate information back in time (i.e., back towards the Sun),meaning they cannot be used in this way.This study investigates the possibility of using variational data assimilation (DA)methods (e.g.,
Dimet and Talagrand [1986];
Courtier et al. [1994]) for assimilating in-situ observations of the solar wind, specifically to update the inner boundary conditionsof the solar wind model. This is achieved using an adjoint method [
Errico , 1997] to mapinformation from the point of observation back to the model’s inner boundary. Suchchanges to the model state would be persistent and remain within the model domain,long after the observation’s timestep. A simple solar wind model is described in thenext section, followed by a description, in general, of variational data assimilation anda derivation of the particular scheme used in this study. Numerical twin experimentsare presented using synthetic observations to test if the method can reconstruct solarwind speed structures and time-series over one solar rotation (i.e. one CarringtonRotation). The data assimilation method is then applied to real spacecraft data,assimilating Solar-Terrestrial Relations Observatory (STEREO) [
Kaiser et al. , 2008]solar wind observations from a time when the spacecraft were well separated fromEarth (2010-2011), and results verified against Advanced Composition Explorer (ACE)[
Stone et al. , 1998] observations in near-Earth space.Future developments are discussed in the concluding section. Here we note,however, that while the requirement for an adjoint model means that the variationalapproach is not expected to be practical for a full magnetohydrodynamic (MHD) solarwind model like Enlil [
Odstrcil , 2003], it is nevertheless of interest for two reasons.Firstly, it is useful in developing new solar wind data assimilation techniques withsimpler models, as described in Section 5, which may be valuable in forecast situations.Secondly, it enables us to test whether updating the inner-boundary conditions of –3–onfidential manuscript submitted to
Space Weather the solar wind model is indeed effective in improving model forecast capability, asdemonstrated in Section 6.3. It is hoped that in the future, ensemble-based Kalmansmoother methods may be able to provide many of the benefits of variational DAwithout the need for an adjoint, and thus can be readily applied to MHD models ofthe solar wind.
The solar wind is a continuous outflow of plasma and magnetic flux which fillsthe heliosphere (e.g.,
Owens and Forsyth [2013]). The solar wind becomes super-magnetosonic within 10 −
20 solar radii ( r S = 695 , km ). Thus forecasting thenear-Earth solar wind conditions is normally treated as a boundary-value problem,with the near-Sun conditions fixed using empirical relations to the coronal magneticfield [ McGregor et al. , 2011;
Riley et al. , 2015], which itself is determined by extrapo-lation from the observed photospheric magnetic field (e.g.,
Mackay and Yeates [2012];
Linker et al. [1999]). The solar wind is then propagated to Earth, typically by anumerical magnetohydrodynamic (MHD) model such as Enlil [
Odstrcil , 2003], withno further observational constraints. In situ spacecraft provide single-point measure-ments of the solar wind and heliospheric magnetic field which can potentially be usedto constrain the solar wind model.In principle, the DA framework developed here will also be applicable to othersolar wind observations, such as Interplanetary Scintillation (IPS; e.g.,
Breen et al. [2006]) or Heliospheric Imagers (HI; e.g.,
Eyles et al. [2009]). These remote measure-ments of solar wind density structures are subject to line-of-sight integration effects,and estimation of solar wind speed further requires some form of correlation tracking.Thus relative to (single-point) in-situ observations, there is increased uncertainty inboth the measurement and its location, but with the advantage of a more synoptic pic-ture of the solar wind. By explicitly accounting for observational errors, a solar windDA scheme can exploit the positives of both forms of data. Though determining theobservational errors is a significant task which is not addressed in the current study.In the next section we describe a simple solar wind propagation tool which per-mits more rapid development of DA techniques than the complex and computationallyexpensive full-MHD approaches [
Odstrcil , 2003]. This approach was recently used toexplore the forecast potential of large ensembles of solar wind solutions with perturbedinitial conditions
Owens and Riley [2017], as is routinely used for operational Numer-ical Weather Prediction (NWP).
In this study, we use the solar wind propagation model of
Riley and Lionello [2011], which maps the equatorial (i.e., two dimensional) solar wind speed over theheliocentric domain from 30 r S to 215 r S from the Sun: v i +1 ,j ( φ ) = v i,j + ∆ r Ω ROT v i,j (cid:18) v i,j +1 − v i,j C ∆ φ (cid:19) (1)where v i,j is the speed (in km/s ) at radius, r i (the i is the radius coordinate), and atCarrington longitude, φ j (where j is the longitude coordinate). Using the same setupas Owens and Riley [2017], ∆ r = 1 r S is the radial grid resolution (in km), ∆ φ = 2 . ◦ is the latitudinal grid resolution, C = π is a constant representing the conversionfactor from degrees to radians and Ω ROT = π . s − is the solar rotational speed.After this solution is obtained, an additional term, v acci,j , is added to the v i,j torepresent the acceleration of the solar wind in the domain considered, which is givenby: v acci,j = αv ,j ( φ ) (cid:16) − e rirH (cid:17) (2) –4–onfidential manuscript submitted to Space Weather where v ,j is the solar wind speed at the inner-boundary and α = 0 .
15 and r H = 50 r S are constants determined by Riley and Lionello [2011].Given the 2-dimensional nature of the solar wind model, in this study we mustassume the ecliptic plane to be equatorial. In reality, the ecliptic is inclined by 7 . ◦ tothe heliographic equator. We note that a solar wind DA scheme in a full 3-dimensionalsolar wind model could relax this assumption. Data assimilation is the study of combining prior knowledge from a model of asystem with information contained in actual observations of that system in order toobtain an optimal estimate of the truth, including its uncertainty. DA methods can beused to 1) provide better model initial conditions for forecasting (e.g. [
Browne and vanLeeuwen , 2015;
Dee et al. , 2011;
Clayton et al. , 2013]); 2) generate optimal evolutiontrajectories for the system to study important physical processes (e.g. [
Broquet et al. ,2011;
Macbean et al. , 2016]); and 3) improve the model physics by studying model-observation misfits [
Lang et al. , 2016].Variational data assimilation refers to the subset of DA methods used extensivelyin meteorological applications [
Sasaki , 1970] that provide an optimal fit over the wholetime window (the period of time over which the data assimilation is applied). Vari-ational DA typically aims to correct the variables under consideration at the initialtime by making use of all the data available over the entire time window. Sequentialassimilation, on the other hand, provides an optimal fit at the end of the window byconsidering each observation sequentially each time a new observation becomes avail-able. Variational data assimilation methods also tend to find the maximum of theposterior probability distribution (the probability distribution of the state given theobservations, see Appendix A: for more details), whereas sequential methods typi-cally seek the mean of the posterior probability distribution. For linear systems, thismeans that the sequential and variational approaches will lead to the same results atthe end of the assimilation window. For non-linear systems, however, the posteriorprobability distribution may have multiple modes, which may lead to the variationalapproach getting stuck in a local maxima, as opposed to the global maxima of the sys-tem, leading to a poorer analysis. Conversely, multi-modal posterior distributions maylead to the mean of the system occurring in an area of low probability in the posteriordistribution, leading to the sequential assimilation approach being sub-standard. Inthese cases, it is unclear what the ‘optimal’ estimate should be, hence there has beensubstantial efforts to develop Particle Filters [
Ades and van Leeuwen , 2013;
Browneand van Leeuwen , 2015;
Zhu et al. , 2016; van Leeuwen , 2014], which aim to estimatethe full posterior probability distribution as opposed to a single ‘optimal’ estimate.With the development of the 4DVar and adjoint model systems [
Dimet and Tala-grand , 1986;
Lorenc , 1986], variational methodologies have become much more viableand efficient methods than optimal interpolation methods based on finite-differencemethods to calculate the gradient of the cost function (and the Kalman filter methodsthat preceded them), particularly for applications within high-dimensional meteorolog-ical models with large quantities of observations. The 4DVar methodology is typicallyused to estimate the initial condition of a model given observations over a fixed time-window. However, the purpose of using a variational approach for the solar wind isto map information contained within observations back closer to the Sun, where wecan then alter the inner boundary of the model and then re-compute the solar windspeed in the whole domain. Specifically, we wish to estimate the solar wind speedat all Carrington longitudes at the inner boundary (30 r S ) using the observations atgreater heliocentric distances (i.e., beyond 30 r S from the Sun, typically at Earth orbit,approximately 215 r S ). In the remainder of this section, we indicate the general defi- –5–onfidential manuscript submitted to Space Weather nitions of data assimilation, followed by specific definitions when applied to the solarwind model. Following this, we explain the methodology behind the variational dataassimilation proposed in this paper.In DA, the variables within a numerical model, x , are described by a N x -dimensional state vector, which contains the values of the quantities of interest atall gridpoints, N . The variables within a state vector for meteorological applicationscould include the temperature field over Africa, the mean sea level pressure over theUK, wind speed/direction, etc., depending upon the purpose of the NWP model. Thetypical dimension of (number of variables within) the state vector in NWP models isof the order ≈ [ Browne and Wilson , 2015;
Tr´emolet , 2006].For this application, the state vector is defined as a vector that contains the solarwind speed at each Carrington longitude, φ j , in the model domain for a given radiuscoordinate, i , and is written as: v i = ( v i, , v i, , . . . , v i,N ) (3)where i = 0 , . . . ,
185 correspond to radii coordinates from the Sun as r i = 30 r S , . . . , r S and N = 128.For a state x (e.g., solar wind speed), an estimate of the initial conditions beforeany data assimilation is referred to as the prior (or background) state, denoted by x b . The prior state is generated using available prior information about the state. Inmeteorological applications, this typically comes from a previous forecast. The priorstate is assumed to be a random perturbation away from the true state, such that: x t0 = x b + ξ (4)where x t0 is an N x -dimensional discretisation of the true initial state, and ξ representsthe random error in the prior state.For our solar wind model, the prior inner boundary condition, v b0 , is an innerboundary condition that we must provide (i.e. from a previous forecast of the in-ner boundary, such as from a previous coronal solution). The prior inner boundarycondition is assumed to be a random perturbation from the true state, such that v b0 ∼ N ( v t0 , B ), where v t0 represents the true speed at the inner boundary and B represents the prior error covariance matrix (i.e. the covariance matrix of the errorsat the inner boundary).Observations within the time window provide information about the true state ofthe system. In order to merge the information contained by the observations with thestate generated by the numerical model, it is necessary to define a function that mapsfrom the state space to the observation space. This function is called the observationoperator and is defined as: y = H ( x ) + (cid:15) (5)where H : R N x → R N y maps the state into observation space, with y being an obser-vation with the observation error given by (cid:15) .Thus for the solar wind, the k th vector of observations, y k , at radius r k , aredefined as: y k = H k ( v i k ) + (cid:15) k (6)where H k is the observation operator, a function that maps the model solar windspeed vector to the observation space (i.e. what the observation would be for anygiven v ) and (cid:15) k is the random observation error, assumed to be normally distributed, (cid:15) k ∼ N ( , R k ). –6–onfidential manuscript submitted to Space Weather
The numerical model, which approximates the dynamics of the system, is denotedby: x i + = f i ( x i ) + η i (7)where x i is the state vector at point i , the numerical model is represented by f i (whichhere is the solar wind propagation model described in section 3), and η i is an N x -dimensional term representing the model error. The distribution of the model erroris almost always unknown and may contain biases towards particular states (i.e. havenon-zero mean) or be multi-modal, etc.We now adopt the Strong Constraint approach and assume that the numericalmodel is perfect, i.e. contains no model error [ Di Lorenzo et al. , 2007;
Fisher et al. ,2005]. In practice, this typically produces a poorer result than the weak constraintsolution that allows for model error. But the weak constraint problem is much morecomplex as it is impossible to know precisely where the model is incorrect, due to miss-ing physics, the effects of sub-grid processes etc., and hence very difficult to prescribean accurate model error covariance matrix. In addition, including model error can leadto coupled equations that can result in different solutions, depending on which orderthe equations are solved [
Evensen et al. , 1998;
Lang et al. , 2016], leading to additionalwork being done to decouple them [
Bennett , 1992]. Therefore, as this is an introduc-tory study to demonstrate the effectiveness of the variational approach, it is logical tostart with the Strong Constraint approximation. The model evolution equation canthen be rewritten as: v i + = f i ( v i ) (8)= ( f i, ( v i ) , f i, ( v i ) , . . . , f i,N ( v i )) (9)where f i,j ( v i ) = v i,j + ∆ r Ω ROT v i,j (cid:16) v i,j +1 − v i,j C ∆ φ (cid:17) + αv ,j (cid:16) e ri − rH − e rirH (cid:17) if r i (cid:54) = 30 r S v i,j + ∆ r Ω ROT v i,j (cid:16) v i,j +1 − v i,j C ∆ φ (cid:17) + αv ,j (cid:16) − e rirH (cid:17) Otherwise(10)where r i is the radial distance from the Sun at radius coordinate, i = 0 , . . . , r S , adds the acceleration termonto the primary calculation for the speed at the next radial coordinate. However, foreach subsequent radial coordinate, the acceleration term from the previous radial pointmust be removed prior to the addition of the new acceleration to avoid an accumulationof the acceleration terms within the model.Use of the Strong Constraint [ Howes et al. , 2017] allows the observation operatorto map from the speed at the inner boundary, v , to the observation location, suchthat: y k = H k [ f i k − ( f i k − ( . . . f ( v ) . . . ))] + (cid:15) k (11)where k denotes the observation number and the i k are a subset of the radial coordi-nates, i , and are the radial coordinates where the k th observation occurs.This allows a cost function, J , to be written purely in terms of the inner boundarysolar wind speed, such that: J ( v ) = 12 (cid:0) v − v b0 (cid:1) T B − (cid:0) v − v b0 (cid:1) + 12 N y (cid:88) k =1 ( y k − H k [ f i k − ( f i k − ( . . . f ( v ) . . . ))]) T R k − ( y k − H k [ f i k − ( f i k − ( . . . f ( v ) . . . ))])(12)where v is the state vector at the inner boundary (at 30 r S ). –7–onfidential manuscript submitted to Space Weather
The cost function is the sum of the relative contributions of the errors presentwithin the solar wind system. Its derivation is outlined in Appendix A: . The firstterm in the cost function represents the prior errors, (cid:0) v − v b0 (cid:1) , at the inner boundaryweighted by the prior error covariance matrix, B − . The second term represents thesum of the observation errors at each observed radial coordinate, ( y k − H k [ f r k − ( f r k − ( . . . f ( v ) . . . ))]),and is weighted by the observation error covariance matrix R − . Therefore, the opti-mal state, that minimises the errors in the system, is the state that minimises the costfunction. The weighting by the inverse of the error covariance matrices means that ifthere is a high amount of certainty in, for example, the prior state compared to that ofthe observations, then B − will be much larger than R − , meaning that the optimalstate vector will have to move closer to the prior state to minimise the cost function,increasing the dominance of the prior state on the final analysis.Obtaining the v that minimises the cost function is a non-trivial task. Thiscan be done by evaluating the finite differences of J or by evaluating ∇ v J directly[ Bannister , 2007]. However, these methods are impractical for higher dimensionalsystems. In such situations, a more efficient method of calculating the gradient is touse the adjoint method, a method of sensitivity analysis that efficiently computes thegradient of a function [
Errico , 1997]. For the Strong Constraint approach, this involvesusing the method of Lagrange multipliers to generate a set of adjoint equations.Firstly, we rewrite the cost function that we wish to minimise as: J ( v ) = 12 (cid:0) v − v b0 (cid:1) B − (cid:0) v − v b0 (cid:1) T + 12 N y (cid:88) k =1 ( y k − H k ( v i k )) R − ( y k − H k ( v i k )) T (13)subject to the (strong) constraint v r + = f r ( v r ) . (14)To minimise this cost function subject to the constraint, we must minimise theLagrangian [ Arfken et al. , 2011], L ( v r , λ r ), which is defined by: L ( v r , λ r ) = J ( v ) + N r (cid:88) i =0 λ T i + ( v i + − f r ( v i )) . (15)where the λ i ’s are the Lagrange multipliers.By differentiating with respect to the v r and the λ r variables, it is possible toobtain a set of adjoint equations, given by: v r + = f r ( v r ) (16) λ N r +1 = (17) λ r = (cid:26) F r T λ r +1 + H k T R k − ( y k − H k ( v r )) if r = r k is obs. radius F r T λ r +1 Otherwise (18) λ = F T λ + B − (cid:0) v − v b0 (cid:1) (19)where F r is given by the Jacobian of f r , such that: F r = ∂f r, ∂v r ( φ ) ∂f r, ∂v r ( φ ) . . . ∂f r, ∂v r ( φ N ) ∂f r, ∂v r ( φ ) ∂f r, ∂v r ( φ ) . . . ∂f r, ∂v r ( φ N ) ... ... . . . ... ∂f r,Nr ∂v r ( φ Nr ) ∂f r,Nr ∂v r ( φ ) . . . ∂f r,Nr ∂v r ( φ N ) (20) –8–onfidential manuscript submitted to Space Weather
Figure 2: The generation of the prior error covariance matrix from the 576-member en-semble is shown in ( a ). The effects of the localisation with a 15 ◦ localisation length scaleare shown in ( b ).and ∂f r,i ∂v r ( φ j ) = − ∆ r Ω ROT C ∆ φ (cid:16) v r ( φ i +1 ) v r ( φ i ) (cid:17) if ρ r > r S and j = i ∆ r Ω ROT C ∆ φ (cid:16) v r ( φ i ) (cid:17) if ρ r > r S and j = i + 11 − ∆ r Ω ROT C ∆ φ (cid:16) v ( φ i +1 ) v ( φ i ) (cid:17) + α (cid:16) − e ρrρH (cid:17) if ρ r = 30 r S and j = i ∆ r Ω ROT C ∆ φ (cid:16) v r ( φ i ) (cid:17) + α (cid:16) − e ρrρH (cid:17) if ρ r = 30 r S and j = i + 1(21)The Lagrange multiplier at the inner boundary radius gives us the negative ofthe gradient at v , the value we wish to find. This implies that: ∇ v J ( v ) = − λ = − F T λ − B − (cid:0) v − v b0 (cid:1) (22)This gradient allows the use of a plethora of available gradient-based minimisa-tion algorithms, from Newton’s methods to steepest descent methods [ Bazaraa et al. ,2013]. In this study, however, we use the Broyden − Fletcher − Goldfarb − Shanno(BFGS) algorithm, a Quasi-Newtonian method, that is commonly used for its speedand accuracy, even for non-linear problems.
For this data assimilation method, we must provide a set of conditions in order toproceed. These include specifying the prior error covariance matrix and the observationerror covariance matrix. The specification of these conditions is a very important partof the data assimilation process. For example, if the specified prior error uncertainty istoo small, the data assimilation analysis will have too much confidence in the qualityof the prior state and the data assimilated analysis state will have less freedom to moveaway from it.In meteorology, the prior error covariance matrix, B , is built up by operationalcentres studying the misfits between forecasts and reanalysis datasets. In comparison,the space weather forecasting is a relatively young science and we do not have thisback-catalogue of available forecasting and reanalysis datasets to compare against. –9–onfidential manuscript submitted to Space Weather
Therefore, a different interim approach must be taken. (Full specification of B for thesolar wind will be the study of future study and here we make an initial estimate topermit progress at this stage.)We here approximate the prior error covariance matrix using an ensemble ofsolar wind speed states generated in the same way as Owens and Riley [2017]. Theinitial conditions to the solar wind propagation model are near-Sun (30 r S ) solar windspeeds, derived from Carrington rotation solutions of the MAS (Magnetohydrodynam-ics Around a Sphere) coronal model Linker et al. [1999] in which the solar wind speedat 30 r S is determined using empirical relations to the coronal magnetic field configura-tion [ Riley et al. r S at the sub-Earth point to Earth orbit at approximately215 r S . An ensemble of 576 members is further generated using perturbed initial con-ditions. This is achieved by sampling the MAS 30 r S solar wind speed at a range oflatitudes about the sub-Earth point (an example is shown in Figure 4 of Owens andRiley [2017]) and independently propagating each perturbed set of initial conditionsto 215 r S . The resulting ensemble of near-Earth solar wind speeds has been shown toprovide an accurate measure of the uncertainty in the unperturbed forecast [ Owensand Riley , 2017].Here, we use the 576-member ensemble created using MAS solutions to Carring-ton rotation 2100, which spans early August to early September 2010. This intervalcontains no interplanetary coronal mass ejections, but does contain both fast and slowsolar wind. The ensemble is used to approximate the B matrix, similar to that of theEnsemble Kalman Filter (EnKF) and methods used in other studies, such as Pereiraet al. [2006]. The B matrix is approximated by: B ≈ M − M (cid:88) m =1 (cid:20)(cid:16) v ( m ) − v M0 (cid:17) (cid:16) v ( m ) − v M0 (cid:17) T (cid:21) (23)where v M0 = M (cid:80) Mm =1 (cid:104) v ( m )0 (cid:105) and M = 576.A common issue in generating the error covariance matrix using a (finite) ensem-ble is that spurious correlations are introduced, as can be seen in Figure 2. The effectof this is that an observation may have an unrealistic impact on a model state variableat a distant location and degrade the quality of the analysis. Indeed, Hamill et al. [2001] showed that if the error in the covariance estimate provided by the ensemble(the noise) is greater than the true correlation (the signal), the accuracy of an EnKFanalysis would decrease. Additionally, they show that the signal-to-noise ratio is afunction of ensemble size, with larger ensembles representative of the true statistics ofthe system, with less associated noise. Therefore, it is necessary to minimise the effectof these spurious correlations in the estimated prior error covariance matrix, B ens . Todo this, we utilise a technique called localisation, that is used to restrict the covari-ances to more realistic spatial/temporal-scales and hence reduce the effective noise tosignal ratios at longer space/time-scales, increasing the effective size of the ensemblein the process.As a first approximation, we use a standard distance-based Gaussian localisationscheme (see Lang et al. [2017] for a discussion on possible localisation schemes for solarwind data assimilation), and apply it to the ensemble-estimated covariance matrix, B ens , such that the localised prior error covariance matrix, B loc , becomes: B loc = L ◦ B ens (24)where ◦ represents the Schur-product and L is the localisation matrix with entriesgiven by L i,j = exp − ( i − j )∆ φS (25) –10–onfidential manuscript submitted to Space Weather and S is the localisation length scale. As a starting point, we set S = 15 ◦ , whichis the approximate width of the slow wind band, and hence mimics the large-scalespatial variability in solar wind speed (e.g., Owens et al. [2017]). The choice of thisvalue is somewhat arbitrary and further research is required. However, a preliminarysensitivity analysis suggests that at least within the range 10 ◦ to 20 ◦ , the resultspresented below are not significantly affected. A finite S is nevertheless required, as thenoise present on the larger spatial-scales ( ≥≈ ◦ ) can cause the adjoint calculationsto become numerically unstable.The observation error covariance matrix, R , is also an unknown quantity in solarwind modelling. Properly addressing this is a substantial research topic in its own right,far beyond the scope of this initial study. In addition to the measurement uncertainty, y , it also accounts for the uncertainty from the “representativity” error [ Janji´c et al. ,2017], resulting from the approximation of a continuous process in discrete space, theerror resulting from representing an observation in the incorrect location, and the errorresulting from representing a single measurement over the (potentially large) volumeof a grid-cell. Additionally, there is also an implicit component of the model errorwithin H ( x ) which further complicates the specification of this quantity. As a firstapproximation of the observation error (to be further tested in future work), we beginwith the simplest approach, assuming that the observations are all independent ofone another (i.e., a diagonal R matrix). We must also specify an observational errorstandard deviation, which we set at 10% of the mean prior solar wind speed at theobservation radius. This is the same order as the observed variability of near-Earthsolar wind speeds. Again, this is still a somewhat arbitrary value, but it provides astarting point from which we can progress and test such assumptions. Then together,this gives: R kk = (0 . v br k ) = (0 . N (cid:88) i =1 (cid:2) v br ( φ i ) (cid:3) ) . (26) This section describes initial tests of the variational data assimilation methodderived in the previous sections. In order to perform numerical experiments in a con-trolled way, we perform observing system simulation experiments (OSSEs, also knownas identical twin experiments), a typical verification approach for a data assimilationscheme. An OSSE uses the uncertainty present within the model state to generate twomodel states that we refer to as the prior state and the truth state. Specifically, wedraw a prior and truth state from the normal distribution, N (cid:16) v M0 , B (cid:17) .We use our numerical model to propagate this truth state to get a complete truth trajectory . We then take synthetic observations from a fixed point in space fromthis truth trajectory. These ‘observations’ are perturbed by measurement noise tomimic a real data assimilation experiment. The goal is to generate a posterior statethat represents the best estimate of the truth run, using the limited information fromthe observations and uncertain prior information on the boundary conditions. Thisposterior state is then compared to the prior state to evaluate the performance of thedata assimilation scheme used.In our experiments, the true state, v t , is generated as a random sample fromthe normal distribution of the MAS ensemble for CR2100, N (cid:16) v M0 , B loc (cid:17) . This truestate is propagated out to 215 r S by the solar propagation model detailed in Section 3and is shown in Figure 3a. The prior state, v b0 , is also randomly drawn from the same –11–onfidential manuscript submitted to Space Weather
OSSE exp. Init. J ( v ) Final J ( v ) No. ofiter. PriorRMSEs(km/s) Post.RMSEs(km/s)“Optimal”prior, un-shifted 760.879 156.364 458 94.075 26.222Prior shifted62∆ φ φ ≈ ◦ and a uniform prior speed of 500 km/s .Figure 3: The speed of the solar wind in km/s propagated from the 30 r S to the 215 r S using the solar wind propagation model. ( a ) shows the speeds generated from the truthstate, ( b ) shows the speeds generated from the prior state and ( c ) shows the posteriorstate after 50 iterations of the forward and adjoint model have been performed. –12–onfidential manuscript submitted to Space Weather
Figure 4: (a) The differences between the prior and the ‘true’ solar wind speed and (b)the differences between the posterior and the true solar wind speed, in km/s .distribution and is shown mapped by the solar wind propagation model to 215 r S inFigure 3b. In this case, the prior state is reasonably close to the true state, with themost obvious difference being the fast stream around 150-200 ◦ Carrington longitudeis larger in the prior state than in the true state. This is seen as a positive band inthe prior error shown in Figure 4a.Direct observations of the true state were taken at 215 r S every ∆ φ , mimickinga time series of near-Earth spacecraft observations every ∼ (cid:15) , drawnfrom the distribution N ( , R ) is added to mimic the effects of observation error. Theprior state is used as an initial estimate for the ‘optimal state’ that minimises the costfunction (equation (12)). The gradient of the cost function at this point is calculatedby the adjoint equations (equations (17)-(22)). This gradient is then input into theBFGS minimisation algorithm to obtain an estimate for the minimum of the costfunction and the procedure is iterated, using the new estimate as the ‘optimal state’,until the cost function converges to a minimal state, within a tolerance of 10 − (i.e.such that the algorithm is repeated as long as the gradient norm of the cost function,from one iteration to the next, is greater than 10 − ).The posterior state generated by the DA method (i.e., the state that minimisesthe cost function) is shown in Figure 3c. An improvement over the prior state (in Figure3b) can be seen in the form of the high speed stream more closely matching that ofthe true state. Figure 4 shows that the posterior error does not suffer from the largepositive error band of the prior. In addition, there are smaller structures present withinthe true speed that are recreated in the posterior at around Carrington Longitudes50 ◦ , ◦ , ◦ and 350 ◦ . Conversely, the finer structures between 200 ◦ − ◦ cannotbe recreated by the DA due to insufficient observational information to overcome thelarge prior errors in this region.The root mean-square error for the OSSEs are calculated over the whole domain,calculated as: RM SE
OSSE = 1 N r N N r (cid:88) i =1 N (cid:88) j =1 (cid:104)(cid:0) v i,j − v ti,j (cid:1) (cid:105) (27) –13–onfidential manuscript submitted to Space Weather
Figure 5: Solar wind in km/s propagated from the 30 r S to the 215 r S using the solarwind propagation model. ( a ) Speeds generated from the truth state, ( b ) speeds generatedfrom the prior state that is shifted by 62∆ φ and ( c ) shows the posterior state after thedata assimilation has been performed.where v i,j represent the prior or posterior state, dependent upon whether the prioror posterior RMSEs are being generated. Table 1 shows that there is a reduction ofapproximately 72% in RMSE as a result of the variational data assimilation analysis. In this section, we consider the effects of different prior states when using dataassimilation. These prior states will be generated such that they are no longer randomperturbations of the true state (as is commonly the case in all applications, especiallyas we rarely have a true state to compare to in the first place). This will cause thedata assimilation to perform sub-optimally as the assumption regarding the prior statewill no longer be valid. Hence the purpose of this section is to show that whilst thedata assimilation may be compromised by poorly specified priors, it can still providesignificant improvements in the estimation of the true state.In this section, the true state, prior error covariance matrix, observations andobservation error covariance matrix are all specified in the same way as the previoussection. –14–onfidential manuscript submitted to
Space Weather
Figure 6: (a) The differences between the prior and the ‘true’ solar wind speeds and ( s )the differences between the posterior and the true solar wind speeds, in km/s , when theprior state is shifted by 62∆ φ .Figure 7: Solar wind speed in km/s propagated from the 30 r S to the 215 r S using thesolar wind propagation model. ( a ) solar wind speeds generated from the truth state, ( b )the uniform prior state of 500 km/s and ( c ) the posterior state after the data assimilationhas been performed. –15–onfidential manuscript submitted to Space Weather
Figure 8: (a)The differences between the prior and the ‘true’ solar wind speeds and (b)the differences between the posterior and the true solar wind speeds, in km/s , when theprior state is specified as a uniform speed of 500 km/s .The first “poor” prior state we use is the prior state generated in the previoussection, but shifted arbitrarily. We use an angular shift of 62∆ φ ≈ ◦ aroundthe Sun, as this results in a very large difference between the prior and true states.The prior is shown in Figure 5b. The fast wind band is now clearly in the incorrectposition in the prior state. This implies that there are far larger errors in our priorstate, reflected by the much higher prior cost functions and prior RMSE values inTable 1 when compared to the unshifted case in the previous section.After the variational data assimilation method is performed, the RMSEs overthe whole domain have been greatly reduced (by approximately 50%). The posteriorstate obtained from the variational data assimilation scheme is shown in Figure 5c anddetailed in Table 1. We can see that a new fast wind band (between ≈ ◦ − ◦ ) hasbeen included in the posterior state in the correct location, but is slightly narrower thanthe true state’s fast wind band. This results in the reduction of most of the negativeprior error in this area. The reason the full band is not generated at the inner boundarymay be due to the presence of a small fast wind region centred at 200 ◦ in the prior state,which is still present and enhanced by the variational data assimilation scheme. Thissmall region yielded low prior error at the observation location in near-Earth space,therefore the data assimilation scheme will not have removed it. A more noticeableerror in the posterior speed is the remains of the prior’s fast wind band, which the dataassimilation has not been able to fully remove. This results in a large positive windbias in the posterior error in this region (although greatly reduced from the initial priorerror in this region, as can be clearly seen in Figure 6). Nevertheless, the near-Earth’observations’ have clearly been able to update the inner boundary conditions for themodel and result in a persistent change in the model state, something which was notpossible with the ensemble-based Kalman filter approach attempted previously [ Langet al. , 2017].The second ”poor” prior state is a constant solar wind speed of 500 km/s atall points on the inner boundary (see Figure 7), which mimics having no near-Suninformation about the solar wind speed structures. It can be seen from Figure 8 thatthis estimate contains a large positive bias from the true state at all points, exceptthe fast wind band, which is negatively biased. In this case, the data assimilation stillreduces the RMSE over the whole domain (by ≈ –16–onfidential manuscript submitted to Space Weather not been able to fully correct for the large positive biases that were present within theprior error. There are still large positive errors over most regions in the posterior state,albeit reduced compared to the prior state. The posterior state can be seen to recreatethe fast wind band in the correct location with low errors present in that region. Inaddition, some of the slower wind regions present in the true state (at ≈ ◦ , ◦ ,200 ◦ and at 300 ◦ ) have been recreated, albeit still with large errors present.The experiments in this section show while the data assimilation is capable ofreducing the errors in a prior state, specifying an accurate prior state is neverthelessessential to obtain an optimal analysis. We further note that in the future when suchdata assimilation is to be used in conjunction with a more costly numerical model, itmay not be possible to perform enough iterations for the cost function to converge,and specifying an accurate prior state will be even more vital. In this section, the data assimilation scheme described in Section 4 is appliedto real STEREO [
Kaiser et al. , 2008] observations during Carrington Rotation 2100.This Carrington Rotation was chosen because the two STEREO spacecraft have alarge separation on either side of the Earth, shown in Figure 9, with STEREO A 80 . ◦ ahead of the Earth and STEREO B located 72 . ◦ behind. For the purposes of thedata assimilation run with the simple solar wind propagation model, the STEREOspacecraft are assumed to be stationary during the Carrington Rotation and locatedat the same heliographic latitude as Earth. Observations are sampled every 5 hours, asthis corresponds to the approximate amount of time for the Sun to rotate ∆ φ degreesand hence matches the implicit time resolution of our solar wind model. We assimilateonly the solar wind speed, as that is the only parameter contained within the solarwind model.The prior error covariance matrix is produced using an ensemble for CarringtonRotation 2100 via the same methodology as in the OSSE experiments. Similarly, theobservation error covariance matrix, R , is generated using equation (26), as in theOSSE experiments.The prior inner boundary initial condition is drawn from the normal distributionof the MAS ensemble for CR2100, N ( , B loc ) as before (and is not further perturbed inany way). Unlike the OSSE experiments, however, we have no ‘truth’ state to compare,so we test the skill of our data assimilated analysis state against an independentset of observations that are taken from the Advanced Composition Explorer (ACE)spacecraft positioned at L
1, approximately 0 . AU ∼ r S from the Sun on theEarth-Sun line, which are sampled every 5 hours from OMNI-hourly averages, as forthe STEREO A and B spacecraft data.Figure 10 shows the prior (i.e., model with no data assimilation) and posterior(i.e., model with STEREO A and B observations assimilated) states in comparisonto solar wind speed observed at the ACE spacecraft. The band of ‘fast wind’, ap-proximately 600 kms − , in the prior state is slower than is observed. But in general,the prior is in relatively good agreement with observations and thus improving onthis forecast provides a challenge for the DA scheme. After the DA scheme is used,the posterior state still underestimates the peak speed of the fast wind stream. Thegreatest improvements achieved by the DA are at the beginning of the CarringtonRotation (between 315 ◦ − ◦ ), where the posterior is in much better agreementwith the observed ACE data. The smoothness of the prior and posterior states is areflection of simplicity of the numerical model used to represent the solar wind. Therelatively shallow gradient of the fast wind onset in the posterior solution also suggestsnumerical diffusion may be an issue. –17–onfidential manuscript submitted to Space Weather
Figure 9: The average STEREO locations for Carrington Rotation 2100, used for thegeneration of the observation operator, defining the location of the observations. Figuregenerated from https://stereo-ssc.nascom.nasa.gov/where.shtml
Figure 10: Solar wind speed at 213 r S ( L
1) for Carrington Rotation 2100. The observedsolar wind from the ACE spacecraft in near-Earth space is shown in black. The solarwind speeds from the model without any data assimilation (the prior) and after assim-ilation of STEREO A and B observations (the posterior) are shown in blue and green,respectively. The MAS ensemble mean solar wind speed is shown in magenta. Carring-ton longitude is shown decreasing to the right, to mimic a time series observed at a fixedpoint in the heliosphere (e.g., in near-Earth space), assuming the solar wind perfectlyco-rotates with the Sun. –18–onfidential manuscript submitted to
Space Weather
Figure 11: Near-Earth solar wind speed (at L
1) for 12 Carrington rotations spanning Oct2010 to Oct 2011. The black line shows observed solar wind speed in near-Earth spaceby the ACE spacecraft. The blue lines show the prior state, generated using the MAS en-semble, as specified in
Owens and Riley [2017]. The green lines show the posterior state,resulting from assimilation of STEREO A and B observations. The magenta lines showthe mean of the MAS ensemble. –19–onfidential manuscript submitted to
Space Weather
Figure 12: The Root Mean Squared Errors (RMSE) in near-Earth (at L
1) solar windspeed as a function of Carrington rotation. Blue line shows the prior state, generatedusing the MAS ensemble, and the green lines shows the RMSEs of the posterior state,resulting from assimilation of STEREO A and B observations. The magenta line showsthe mean of the MAS ensemble.Carr.Rot STER.A loc.( ◦ fromEarth) STER.B loc.( ◦ fromEarth) Init. J ( v ) Final J ( v ) MASMeanRMSE(km/s) PriorRMSE(km/s) Post.RMSE(km/s) Post./MASMeanred. inRMSE(%) Post./Priorred. inRMSE(%)2100 80.6 -72.8 799.54 381.07 70.60 99.19 74.81 -5.96 24.572101 82.1 -76.0 918.91 570.15 82.69 79.46 71.88 13.07 9.542102 83.2 -79.6 1704.32 768.73 82.41 94.64 65.89 20.05 30.382103 84.1 -83.6 2866.68 496.45 83.50 79.36 99.89 -19.63 -26.872104 85.2 -87.1 1435.50 444.03 93.12 122.96 101.60 -9.09 17.392105 86.0 -90.3 2036.20 515.61 61.76 91.54 67.44 -9.20 26.322106 86.8 -92.9 2429.19 598.49 81.06 106.95 88.67 -9.39 17.092107 87.5 -94.6 1616.64 1030.86 119.06 127.67 125.43 -5.35 1.752108 88.7 -95.3 3086.53 1384.29 78.78 74.09 68.47 13.09 7.592109 90.6 -95.0 1585.78 684.77 92.75 93.14 85.63 7.68 8.062110 93.3 -93.7 2756.18 740.25 112.91 119.32 65.84 41.69 44.832111 96.3 -92.5 4563.84 1461.13 121.141 125.05 79.23 34.60 36.642112 99.1 -92.0 2353.79 1028.48 164.82 154.00 91.65 44.39 40.49Avg. 87.96 -88.11 2165.62 777.25 95.74 105.18 83.57 8.92 18.37Table 2: Table showing prior and posterior cost function values, the RMSEs of the MASMean, Prior and Posterior states for each Carrington Rotation and the percentage reduc-tions in RMSE between the MAS Mean/prior state and the posterior states. The bottomrow shows the average values across all 13 Carrington rotations. –20–onfidential manuscript submitted to Space Weather
In addition, from the first row of Table 2, we see that the overall errors in thesystem, as denoted by the cost function, have been reduced by ≈ RM SE = (cid:118)(cid:117)(cid:117)(cid:116) N ACE N
ACE (cid:88) i =1 (cid:104) ( y i − x i ) (cid:105) (28)where the y i and x i correspond to the ACE observations and state, at their observationtime/location.It can be seen that assimilating the STEREO A and B observations has resultedin a reduction in the near-Earth solar wind speed RMSE of ≈
25% over the wholeCarrington Rotation. While the prior state is a good representation of the true statein this instance, the variational scheme was still able to produce a lasting change inthe model state by correctly updating the inner boundary condition and improvingthe near-Earth solar wind speed. For this particular Carrington rotation, the MASensemble mean (without STEREO data assimilation) provides a very good match tothe observed near-Earth solar wind. In particular, the peak of the fast wind band(between 160 ◦ − ◦ ), is better reproduced by the MAS ensemble mean than by theposterior. The posterior does better capture the structures present in the observedsolar wind speed between 360 ◦ − ◦ and 90 ◦ − ◦ . However, the overall RMSEof the ensemble mean is lower than for the posterior. Thus, in this instance, theimprovements to the prior state provided by the assimilation of STEREO data arenot as great as considering the full ensemble, suggesting the prior error covariancematrix for this Carrington rotation is underestimated, perhaps due to over-aggressivelocalisation.We now expand the analysis to 12 further Carrington rotations spanning theperiod October 2010 to October 2011. For each Carrington rotation, the STEREOspacecraft positions are updated (see 2 nd and 3 rd columns of table 2), as are the B matrices using a new MAS solar wind speed ensemble as described in Section 5.The observation error covariance matrices, R , are also updated for each CarringtonRotation, and are dependent upon the prior states, that are generated as a randomperturbation (distributed by the normal distribution N ( , B )) from the MAS ensemblemean.The results are shown in Figures 11 and 12. The cost function values and RMSEsfor the MAS ensemble mean, posterior and prior are listed in Table 2. For 12 of the13 Carrington Rotations, the use of data assimilation leads to an improvement in theestimate of the near-Earth solar wind speed compared to the no-DA model state. Thereduction in RMSE is in the range 2% − Cane and Richardson , 2003] or at theSTEREO spacecraft during this period.In comparison with the MAS ensemble mean (the best estimate of the truthfrom the MAS ensemble prior to DA), the posterior has greater RMSE for Carring- –21–onfidential manuscript submitted to
Space Weather ton Rotations 2103, 2104, 2105, 2106 and 2107. These RMSE increases are relativelysmall in magnitude compared to the improvements from DA in other CRs, particularlyCarrington rotations 2108 through 2112. The high RMSEs relative to the MAS en-semble mean could be a result of the latitudinal difference in the STEREO and ACEspacecraft (which the 2-dimensional model assumes all lie in the equatorial plane),which could lead to sampling of different solar wind structures. If so, this could beat least partly mitigated by the use of a fully 3-dimensional solar wind model, thoughheliospheric latitudinal localisation may need to be treated in a different manner tothe heliospheric longitudinal localisation considered here. The MAS ensemble meanseems to perform better when there are no transients in the solar wind and when thesolar wind is relatively steady-state. The data assimilated solar wind speed performsincreasingly better than the MAS ensemble mean when there is increased CR to CRvariability, as the solar cycle increases towards the end of the period considered.Averaged across all 13 Carrington rotations considered here, data assimilation ofSTEREO data results in 18 .
4% reduction in near-Earth solar wind RMSE compared tothe prior state and an 8 .
9% reduction in RMSE compared to the MAS ensemble mean.In Figures B.1 through B.4 and Tables B.1 and B.2, we show that at the observationlocations themselves, the gains from DA are greater still. For STEREO A there is a42 .
7% reduction in the prior RMSE and a 35 .
0% reduction in the MAS mean RMSEas a result of the data assimilation. For STEREO B, there is a 38 .
1% reduction in theprior RMSE and a 29 .
8% reduction in the MAS mean RMSE as a result of the dataassimilation.
Whilst this study shows there is a great deal of potential available in the use ofvariational data assimilation when applied to the solar wind, there are many hurdlesthat still need to be overcome. It is also noted that the adjoint method is an extremelypowerful method once generated, though this is a substantially simpler task for thesolar wind propagation tool used here than for a full MHD model of the solar wind,as discussed below.Variational data assimilation methods require the generation of an accurate B matrix, which may be the reason for the lesser improvements noted in CR − B , may not fully represent the errors present in theprior state due to the linear dependence of the ensemble members. This may leadto a prior covariance matrix that is of too low rank for the gradient to be calculatedaccurately by the adjoint method. As an alternative to using a different covariancematrix for each Carrington Rotation, a ‘climatological’ covariance matrix could be gen-erated using MAS ensembles from multiple Carrington Rotations. This will increasethe amount of information contained within the B -matrix and will negate the need forlocalisation. However, this solution negates the advantage of having a flow-dependentcovariance matrix that is specific to each Carrington Rotation. It may be the casethat this ‘flow dependency’ is necessary for effective B matrix generation and that aCarrington-Rotation-specific prior error uncertainty matrix may need to be created foroptimal assimilation. This mirrors the ‘hybrid’ data assimilation schemes currently be-ing investigated in numerical weather prediction [ Goodliff et al. , 2015;
Bonavita et al. ,2016]. In both cases, the accuracy of the B -matrix relies both on the accuracy of thecoronal model for defining the near-Sun solar wind, and the method of sampling thenear-Sun solar wind to produce the ensemble. Further research is required to quan-tify the potential improvements from using Carrington-rotation specific B matrices,‘climatological’ B matrices and the sensitivity to the accuracy of the coronal modelused. The adjoint model is unique to each numerical model and can be extremelydifficult to compute efficiently for high dimensional models. Furthermore, as mentioned –22–onfidential manuscript submitted to Space Weather in Section 4, the variational data assimilation methods rely upon the linearisation ofthe numerical model. The simple numerical model in this study is reasonably linear,and therefore long assimilation windows of one full Carrington Rotation (27 days)are possible. For more complex, higher dimensional models that are more nonlinear,the assimilation window will need to be much shorter to avoid computational issues.For numerical weather prediction, the typical window length is 6 −
12 hours, but it isunclear how long the assimilation window can and should be for complex MHD models,such as Enlil (e.g. [
Odstrcil , 2003;
Odstrcil and Pizzo , 1999;
Odstrcil et al. , 2004]) orEUHFORIA [
Poedts and Pomoell , 2017]. This depends on the spatial resolution of themodel, as higher resolutions typically lead to stronger nonlinearities with shorter errorgrowth timescales, and also how close the model is to the true solution. Numericalweather prediction owes much of its forecast accuracy to the enormous amount ofobservations in each 6-hour window (approximately 10 observations). This also meansthat the initial boundary condition is accurate as well, such that the linearisations workwell for longer. Unfortunately, the sparsity of observations for space weather means itcannot readily adopt this solution. The experiments shown in this study are the first, to the authors’ knowledge,application of a variational data assimilation (DA) method applied to the solar windusing in-situ spacecraft observations. This DA method maps observational informationfrom 215 r S back towards the Sun, to the inner boundary of the solar wind model at30 r S . Twin experiments showed that this variational DA is able to reduce the errorin the dynamical model using 5-hourly observations in the near-Earth space. Byimproving the inner boundary conditions of the solar wind speed, it is possible toretrieve the structure of the solar wind speed produced by the numerical model inthe entire domain between the Sun and Earth (at least within the spacecraft orbitalplane). This allows solar wind structures, such as bands of fast and slow solar wind,to be reconstructed by the data assimilation method.This variational DA scheme requires an estimate of the prior error uncertaintymatrix, B . In meteorological applications, there are many decades of forecasts andreanalysis datasets that enable a more accurate representation of the prior error uncer-tainty matrix. In this study, we have used an ensemble of inner boundary conditionsto estimate the B matrix, over a single Carrington Rotation. It may be more usefulto generate a B matrix with more ensemble members, over multiple Carrington ro-tations, to incorporate more information about possible error structures in the solarwind. This has the benefit of incorporating additional information from the coronalmodel solutions, which generally capture the overall solar wind structure, but sufferfrom relatively small positional errors. Further research is also required in order tobetter determine how to accurately generate the observation error covariance matrix, R , and how to properly incorporate representivity errors into this matrix.By applying the same data assimilation to three different prior states we havedemonstrated the importance of a “good” prior solution from the coronal models. Thusforecast skill gained from assimilation of solar wind observations will further increasewith developments in coronal modelling and initialisation, such as from improved pho-tospheric magnetic field characterisation [ Arge et al. , 2010]. In this study, we focussedon using synthetic and real in-situ spacecraft observations, though in principle thesame data assimilation framework can be applied to remote measurements of the solarwind, such as interplanetary scintillation and white-light heliospheric imager observa-tions. The difficulty, however, lies in accurately representing the observational error,both in the solar wind speed measurement and in the generation of the observationoperator, H (as both techniques involve some degree of line-of-sight integration). Thelatter uncertainty could possibly be approximated within a data assimilation approach –23–onfidential manuscript submitted to Space Weather by weak-but-widespread localisation, but such techniques will form part of a futurestudy.The DA scheme was also used to make hindcasts of near-Earth solar wind speedby assimilating approximately 1 year of in-situ observations from the STEREO space-craft, when they were approximately 80 ◦ ahead and behind the Earth in its orbit withrespect to the Sun. The prior state has been generated from the MAS ensemble, thusthe information contained within the STEREO A and B observations is not carriedon to subsequent Carrington Rotations. The DA analysis can, however, be used toforecast the inner boundary for the subsequent Carrington Rotation. This means thatthere are two possible estimates of the inner boundary for the subsequent Carring-ton Rotation, one using information from within the inner boundary (from the MASmodel) and one estimate using information at heliocentric distances beyond the in-ner boundary (from the assimilation of in situ solar wind data). How we consolidatethese disparate estimates is a question for future research. Nevertheless, even with theapproach used here, the root mean-square error in the near-Earth solar wind speedhindcasts was reduced by 18% by the use of STEREO data assimilation.The interval considered in this study (August 2010 to August 2011) is close tosolar minimum, thus the coronal mass ejection rate is much reduced compared tosolar maximum. This interval, however, still features 9 interplanetary coronal massejections (ICMEs) observed in near-Earth space [ Cane and Richardson , 2003, ; thoughnone with speeds above 600 km/s ]. Given the large angular separation of the STEREOspacecraft from Earth, it is unlikely these same ICMEs would be present in the assim-ilated data. The localisation required for such transient structures will be significantlydifferent (in both space and time) than for the ambient, steady-state solar wind. Forexample, assuming the fast solar wind from an ICME persists at the same Carringtonlongitude for a whole solar rotation will significantly degrade a forecast. Fast ICMEscould result in ”missed” high speed streams in the posterior solution. Conversely,ICMEs seen by the STEREO spacecraft which didn’t encounter Earth, could result in”false” high speed streams in the posterior solution. With very limited spatial sam-pling from in situ spacecraft, such transient events of limited spatial extent will beproblematic for any data assimilation scheme. The solution to dealing with transientstructures will likely be greater observational sampling (both in space and time), suchas with the remotely sensed observations discussed above. But weighting the relativemerits of the two data types, namely precision point measurements and non-localisedsynoptic measurements, will require careful attention in the observation error matrix.In practice, however, the importance of in-situ observations, particularly in terms ofthe magnetic field observations, may mean automated ICME detection algorithms alsoneed to be applied in real time.The major outstanding issue with using a variational data assimilation scheme inan operational setting is that the adjoint method is not scalable to higher-dimensionalor more complex models. The tangent linear and adjoint models are unique to eachnumerical model and can be extremely difficult and require many human-years to cre-ate an efficient scheme (however, once created, the adjoint method is an extremelypowerful and efficient tool, as shown in this paper). This means that without substan-tial investment, it will not be possible to utilise an adjoint model in a full MHD model,such as Enlil or EUHFORIA. This issue indicates that it is perhaps more useful to usethe adjoint-based data assimilation methods with smaller, simpler solar wind/MHDmodels. This approach could be define an optimum set of boundary condition usingsolar wind observations which can then be used to drive more complex models. ForDA within more complex MHD models, in order to map observational informationfrom near-Earth space to the inner-boundary, it is perhaps more informative/usefulto use a hybrid data assimilation method, such as Ensemble-4DVar [
Amezcua et al. ,2017;
Goodliff et al. , 2017], or a smoother-based data assimilation method, such as the –24–onfidential manuscript submitted to
Space Weather
Iterative Ensemble Kalman Smoother [
Bocquet and Sakov , 2014]. These will be testedin future studies.
A: Deriving Strong Constraint Cost Function
In order to combine the observational information within y with the state vector, x , we appeal to Bayes’ Theorem [ Bayes and Price , 1763], which states that: p ( x | y ) = p ( y | x ) p ( x ) p ( y ) (A.1)where p ( x | y ) is the posterior probability distribution function, the probability of thestate, x , occurring given the observation, y and is the distribution we wish to estimatein data assimilation; p ( y | x ) is the likelihood distribution, which is the probability of theobservation occurring for any given state; p ( x ) is the prior distribution, the probabilityof the state occurring and p ( y ) is the probability of the observation occurring and isa constant, normalising factor for any given observation.The posterior distribution gives the full picture of the probability of the statespace, given the observational data and is what we wish to estimate. Typically, how-ever, the state dimension is extremely large, for example, for numerical weather predic-tion the state vector is of the order 10 -dimensional vector. This makes it impossiblefor the full posterior distribution to be computed explicitly. Therefore, assumptionsabout the prior distributions, the likelihood distributions and the numerical modelmust be made to simplify the problem. Furthermore, we must define precisely whatwe mean by the ‘optimal estimate’.In the variational framework, ‘optimal state’ is defined as the state which max-imises the posterior probability distribution (the mode of the posterior distribution)[ Dimet and Talagrand , 1986]. In addition, the prior, observation and model errors areall assumed to be unbiased Gaussian distributions. This means that the probabilitydistributions of the prior error, ξ , the observation error, (cid:15) k , and the model error, η i probability distributions can be fully defined by their mean and respective errorcovariance matrices, such that: p ( ξ ) = (cid:16)(cid:112) π | B | (cid:17) − N x exp (cid:18) − ξ T B − ξ (cid:19) (A.2) p ( η i ) = (cid:16)(cid:112) π | Q i | (cid:17) − N x exp (cid:18) − η i T Q i − η i (cid:19) (A.3) p ( (cid:15) k ) = (cid:16)(cid:112) π | R k | (cid:17) − N y exp (cid:18) − (cid:15) k T R k − (cid:15) k (cid:19) (A.4)where B is the prior error covariance matrix, Q i is the model error covariance matrix foreach point i and R k is the observation error covariance matrix for the k th observation.In this study, as we are performing initial data assimilation experiments in thesolar wind, we shall make the further assumption that the model error is zero, i.e.,the perfect model assumption. This is known as the Strong-Constraint. Whilst it ispossible to define the variational problem with model error [ Lang et al. , 2016;
Scheichlet al. , 2013], the so-called Weak-Constraint approach, it is beyond the scope of thispaper and will not be discussed here.As the model error is assumed equal to zero, x i can be written explicitly interms of the initial condition, x , such that x i = f i − ( f i − ( . . . f ( x ) . . . )). Usingthe probability distribution functions described by equations (A.2)-(A.4), it is possible –25–onfidential manuscript submitted to Space Weather to write the likelihood and prior distributions as: p ( y , . . . , y N y | x ) = N y (cid:89) k =0 (cid:20)(cid:16)(cid:112) π | R k | (cid:17) − N y exp (cid:18) −
12 ( y k − H k ( x )) R k − ( y k − H k ( x )) (cid:19)(cid:21) (A.5) p ( x ) = (cid:16)(cid:112) π | B | (cid:17) − N x exp (cid:18) − (cid:16)(cid:0) x − x b (cid:1) T B − (cid:0) x − x b (cid:1)(cid:17)(cid:19) . (A.6)where H k now implicitly contains the numerical model f i .As p ( y ) is constant, this implies that the posterior probability can be written as: p ( x | y , . . . , y N y ) ∝ exp ( −J ( x )) (A.7)where J ( x ) = 12 (cid:0) x − x b (cid:1) T B − (cid:0) x − x b (cid:1) + 12 N y (cid:88) k =0 (cid:2) ( y k − H k ( x )) R k − ( y k − H k ( x )) (cid:3) (A.8)is the cost function. It can be seen that by minimising J ( x ), we maximise theposterior probability distribution, which is what we wish to find.There are many minimisation algorithms available to do this. The majority of theefficient methods to do this require the gradient/Hessian to be computed, which areoften the most difficult parts to obtain estimates for, especially in higher dimensions. –26–onfidential manuscript submitted to Space Weather
B: Plots/tables of Solar wind speeds at STEREO A and B, over Car-rington Rotations 2101-2112
B.1 STEREO A
Figure B.1: Near-Earth solar wind speed for 12 Carrington rotations spanning Oct 2010to Oct 2011. The black line shows observed solar wind speed at the STEREO A satellitelocation. The blue lines show the prior state, generated using the MAS ensemble, as spec-ified in
Owens and Riley [2017], and the magenta line is the MAS ensemble mean. Thegreen lines show the posterior state. –27–onfidential manuscript submitted to
Space Weather
Figure B.2: The Root Mean Squared Errors (RMSE) in near-Earth solar wind speed as afunction of Carrington rotation. Blue line shows the prior state, generated using the MASensemble, and the green lines shows the RMSEs of the posterior state, resulting fromassimilation of STEREO A and B observations.Carr.Rot. STERA pos. MASRMSE(km/s) PriorRMSE(km/s) PostRMSE(km/s) MAS/PriorRMSERed.(%) Post/PriorRMSERed.(%)2100 80.6 81.12 77.85 50.68 37.52 34.902101 82.1 75.66 77.72 63.25 16.40 18.622102 83.2 67.11 99.11 65.22 2.82 34.192103 84.1 103.26 149.86 52.94 48.73 64.672104 85.2 70.86 89.17 44.88 36.66 49.672105 86 76.61 98.73 51.37 32.95 47.972106 86.8 102.77 130.98 69.08 32.78 47.262107 87.5 95.75 96.69 75.11 21.56 22.322108 88.7 112.13 120.24 81.45 27.36 32.262109 90.6 88.42 95.56 52.15 41.02 45.432110 93.3 120.96 126.39 51.97 57.04 58.882111 96.3 175.45 177.94 81.23 53.70 54.352112 99.1 124.77 120.27 67.46 45.93 43.91Avg. 87.96 99.61 112.35 62.06 34.96 42.65Table B.1: Table showing the RMSEs of the MAS Mean, Prior and Posterior states foreach Carrington Rotation and the percentage reductions in RMSE between the MASMean/prior state and the posterior states at STEREO A’s locations. The bottom rowshows the average values across all 12 Carrington rotations. –28–onfidential manuscript submitted to
Space Weather
B.2 STEREO B
Figure B.3: Near-Earth solar wind speed for 12 Carrington rotations spanning Oct 2010to Oct 2011. The black line shows observed solar wind speed at the STEREO B satellitelocation. The blue lines show the prior state, generated using the MAS ensemble, as spec-ified in
Owens and Riley [2017], and the magenta line is the MAS ensemble mean. Thegreen lines show the posterior state. –29–onfidential manuscript submitted to
Space Weather
Figure B.4: The Root Mean Squared Errors (RMSE) in near-Earth solar wind speed as afunction of Carrington rotation. Blue line shows the prior state, generated using the MASensemble, and the green lines shows the RMSEs of the posterior state, resulting fromassimilation of STEREO A and B observations.CarrRot. STERB pos MASRMSE(km/s) PriorRMSE(km/s) PostRMSE(km/s) MAS/PriorRMSERed.(%) Post/PriorRMSERed.(%)2100 -72.8 63.02 44.34 37.92 39.83 14.482101 -76.0 59.48 65.56 45.23 23.96 31.012102 -79.6 88.71 101.49 68.43 22.86 32.572103 -83.6 74.33 108.31 53.32 28.27 50.772104 -87.1 69.04 96.43 52.06 24.59 46.012105 -90.3 73.83 114.12 79.71 -7.96 30.152106 -92.9 64.47 111.36 39.33 38.99 64.682107 -94.6 102.80 103.28 82.28 19.96 20.332108 -95.3 115.54 137.30 67.63 41.47 50.742109 -95.0 103.43 102.80 69.33 32.97 32.562110 -93.7 125.33 136.90 69.97 44.17 48.892111 -92.5 114.96 122.44 68.35 40.54 44.182112 -92.0 134.54 119.16 84.47 37.22 29.11Avg. -88.11 91.50 104.88 62.93 29.76 38.11Table B.2: Near-Earth solar wind speed for 12 Carrington rotations spanning Oct 2010to Oct 2011. The black line shows observed solar wind speed at the STEREO B satellitelocation. The blue lines show the prior state, generated using the MAS ensemble, as spec-ified in
Owens and Riley [2017], and the magenta line is the MAS ensemble mean. Thegreen lines show the posterior state. –30–onfidential manuscript submitted to
Space Weather
Acknowledgments
References
Ades, M., and P. J. van Leeuwen (2013), An exploration of the equivalent weights parti-cle filter,
Quarterly Journal of the Royal Meteorological Society , (672), 820–840,doi:10.1002/qj.1995.Amezcua, J., M. Goodliff, and P. J. Van Leeuwen (2017), A weak-constraint 4DEnsem-bleVar. Part I: formulation and simple model experiments, Tellus A: Dynamic Me-teorology and Oceanography , (1), 1271,564, doi:10.1080/16000870.2016.1271564.Arfken, G. B., H. J. Weber, and F. E. Harris (2011), Mathematical methods for physi-cists: A comprehensive guide , Academic press.Arge, C. N., C. J. Henney, J. Koller, C. R. Compeau, S. Young, D. MacKenzie,A. Fay, J. W. Harvey, M. Maksimovic, K. Issautier, N. Meyer-Vernet, M. Moncuquet,and F. Pantellini (2010), Air Force Data Assimilative Photospheric Flux Transport(ADAPT) Model, in
AIP Conference Proceedings , vol. 1216, pp. 343–346, AmericanInstitute of Physics, doi:10.1063/1.3395870.Bannister, R. (2007), Elementary 4D-Var,
DARC Technical Report No.2 .Bayes, T., and R. Price (1763), An Essay towards solving a Problem in the Doctrine ofChances. By the late Rev. Mr. Bayes, FRS communicated by Mr. Price, in a letterto John Canton, AMFRS,
Philosophical Transactions (1683-1775) , pp. 370–418.Bazaraa, M. S., H. D. Sherali, and C. M. Shetty (2013),
Nonlinear Programming The-ory and Algorithms , John Wiley & Sons.Bennett, A. F. (1992),
Inverse Methods in Physical Oceanography , Arnold and Caro-line Rose Monograph Series of the American So, Cambridge University Press.Bocquet, M., and P. Sakov (2014), An iterative ensemble Kalman smoother,
QuarterlyJournal of the Royal Meteorological Society , (682), 1521–1535, doi:10.1002/qj.2236.Bonavita, M., E. H´olm, L. Isaksen, and M. Fisher (2016), The evolution of theECMWF hybrid data assimilation system, Quarterly Journal of the Royal Meteo-rological Society , (694), 287–303, doi:10.1002/qj.2652.Breen, A. R., R. A. Fallows, M. M. Bisi, P. Thomasson, C. A. Jordan, G. Wannberg,and R. A. Jones (2006), Extremely long baseline interplanetary scintillation measure-ments of solar wind velocity, Journal of Geophysical Research , (A8), A08,104,doi:10.1029/2005JA011485.Broquet, G., F. Chevallier, P. Rayner, C. Aulagnier, I. Pison, M. Ramonet,M. Schmidt, A. T. Vermeulen, and P. Ciais (2011), A European summertimeCO 2 biogenic flux inversion at mesoscale from continuous in situ mixing ra-tio measurements, Journal of Geophysical Research: Atmospheres , (D23), doi:10.1029/2011JD016202.Browne, P. A., and P. J. van Leeuwen (2015), Twin experiments with the equivalentweights particle filter and HadCM3, Quarterly Journal of the Royal MeteorologicalSociety , (693), 3399–3414, doi:10.1002/qj.2621.Browne, P. A., and S. Wilson (2015), A simple method for integrating a complex modelinto an ensemble data assimilation system using MPI, Environmental Modelling &Software , , 122–128.Bust, G. S., and C. N. Mitchell (2008), History, current state, and future direc-tions of ionospheric imaging, Reviews of Geophysics , (1), RG1003, doi:10.1029/ –31–onfidential manuscript submitted to Space Weather
Journal of Geophysical Research , (A4),1156, doi:10.1029/2002JA009817.Clayton, A. M., A. C. Lorenc, and D. M. Barker (2013), Operational implementationof a hybrid ensemble/4D-Var global data assimilation system at the Met Office, Quarterly Journal of the Royal Meteorological Society , (675), 1445–1461, doi:10.1002/qj.2054.Courtier, P., J. N. Th´epaut, and A. Hollingsworth (1994), A strategy for operationalimplementation of 4D-Var, using an incremental approach, Quarterly Journal of theRoyal Meteorological Society , (519), 1367–1387.Dee, D. P., S. M. Uppala, A. J. Simmons, P. Berrisford, P. Poli, S. Kobayashi, U. An-drae, M. A. Balmaseda, G. Balsamo, P. Bauer, P. Bechtold, A. C. M. Beljaars,L. van de Berg, J. Bidlot, N. Bormann, C. Delsol, R. Dragani, M. Fuentes, A. J.Geer, L. Haimberger, S. B. Healy, H. Hersbach, E. V. H´olm, L. Isaksen, P. K˚allberg,M. K¨ohler, M. Matricardi, A. P. McNally, B. M. Monge-Sanz, J.-J. Morcrette, B.-K.Park, C. Peubey, P. de Rosnay, C. Tavolato, J.-N. Th´epaut, and F. Vitart (2011),The ERA-Interim reanalysis: configuration and performance of the data assimilationsystem, Quarterly Journal of the Royal Meteorological Society , (656), 553–597,doi:10.1002/qj.828.Di Lorenzo, E., A. M. Moore, H. G. Arango, B. D. Cornuelle, A. J. Miller, B. Powell,B. S. Chua, and A. F. Bennett (2007), Weak and strong constraint data assimilationin the inverse Regional Ocean Modeling System (ROMS): Development and appli-cation for a baroclinic coastal upwelling system, Ocean Modelling , (3-4), 160–187,doi:10.1016/J.OCEMOD.2006.08.002.Dimet, F. X. L., and O. Talagrand (1986), Variational algorithms for analysis andassimilation of meteorological observations: theoretical aspects, Tellus A , (2),97–110.Durazo, J. A., E. J. Kostelich, and A. Mahalov (2017), Local ensemble transformKalman filter for ionospheric data assimilation: Observation influence analysisduring a geomagnetic storm event, Journal of Geophysical Research: Space Physics , (9), 9652–9669, doi:10.1002/2017JA024274.Errico, R. M. (1997), What is an adjoint model?, Bulletin of the American Meteoro-logical Society , (11), 2577–2591.Evensen, G., D. P. Dee, and J. Schr¨oter (1998), Parameter estimation in dynamicalmodels, in Ocean Modeling and Parameterization , pp. 373–398, Springer.Eyles, C. J., R. A. Harrison, C. J. Davis, N. R. Waltham, B. M. Shaughnessy, H. C. A.Mapson-Menard, D. Bewsher, S. R. Crothers, J. A. Davies, G. M. Simnett, R. A.Howard, J. D. Moses, J. S. Newmark, D. G. Socker, J.-P. Halain, J.-M. Defise,E. Mazy, and P. Rochus (2009), The Heliospheric Imagers Onboard the STEREOMission,
Solar Physics , (2), 387–445, doi:10.1007/s11207-008-9299-0.Fisher, M., M. Leutbecher, and G. A. Kelly (2005), On the equivalence betweenKalman smoothing and weak-constraint four-dimensional variational data assimi-lation, Quarterly Journal of the Royal Meteorological Society , (613), 3235–3246,doi:10.1256/qj.04.142.Goodliff, M., J. Amezcua, and P. J. van Leeuwen (2015), Comparing hybrid dataassimilation methods on the Lorenz 1963 model with increasing non-linearity, TellusA , (0).Goodliff, M., J. Amezcua, and P. J. Van Leeuwen (2017), A weak-constraint 4DEnsem-bleVar. Part II: experiments with larger models, Tellus A: Dynamic Meteorologyand Oceanography , (1), 1271,565, doi:10.1080/16000870.2016.1271565.Hamill, T. M., J. S. Whitaker, and C. Snyder (2001), Distance-dependent filteringof background error covariance estimates in an ensemble Kalman filter, MonthlyWeather Review , (11), 2776–2790. –32–onfidential manuscript submitted to Space Weather
Howes, K. E., A. M. Fowler, and A. S. Lawless (2017), Accounting for model error instrong-constraint 4D-Var data assimilation,
Quarterly Journal of the Royal Meteo-rological Society , (704), 1227–1240, doi:10.1002/qj.2996.Janji´c, T., N. Bormann, M. Bocquet, J. A. Carton, S. E. Cohn, S. L. Dance, S. N.Losa, N. K. Nichols, R. Potthast, J. A. Waller, and P. Weston (2017), On the repre-sentation error in data assimilation, Quarterly Journal of the Royal MeteorologicalSociety , doi:10.1002/qj.3130.Kaiser, M. L., T. A. Kucera, J. M. Davila, O. C. S. Cyr, M. Guhathakurta, andE. Christian (2008), The STEREO mission: An introduction, in
The STEREO Mis-sion , pp. 5–16, Springer.Kalnay, E. (2003),
Atmospheric modeling, data assimilation and predictability , Cam-bridge university press.Lang, M., P. J. van Leeuwen, and P. A. Browne (2016), A systematic method ofparameterisation estimation using data assimilation,
Tellus A , (0).Lang, M., P. Browne, P. J. van Leeuwen, and M. Owens (2017), Data Assimilation inthe Solar Wind: Challenges and First Results, Space Weather , (11), 1490–1510,doi:10.1002/2017SW001681.Linker, J. A., Z. Miki´c, D. A. Biesecker, R. J. Forsyth, S. E. Gibson, A. J. Lazarus,A. Lecinski, P. Riley, A. Szabo, and B. J. Thompson (1999), Magnetohydrodynamicmodeling of the solar corona during Whole Sun Month, Journal of Geophysical Re-search: Space Physics , (A5), 9809–9830, doi:10.1029/1998JA900159.Lorenc, A. C. (1986), Analysis methods for numerical weather prediction, QuarterlyJournal of the Royal Meteorological Society , (474), 1177–1194.Macbean, N., P. Peylin, F. Chevallier, M. Scholze, and G. Sch¨urmann (2016), Con-sistent assimilation of multiple data streams in a carbon cycle data assimilationsystem, Geosci. Model Dev , , 3569–3588, doi:10.5194/gmd-9-3569-2016.Mackay, D., and A. Yeates (2012), The Sun’s Global Photospheric and Coronal Mag-netic Fields: Observations and Models, Living Reviews in Solar Physics , (1), 6,doi:10.12942/lrsp-2012-6.McGregor, S. L., W. J. Hughes, C. N. Arge, M. J. Owens, and D. Odstrcil (2011), Thedistribution of solar wind speeds during solar minimum: Calibration for numericalsolar wind modeling constraints on the source of the slow solar wind, Journal ofGeophysical Research: Space Physics , (A3), doi:10.1029/2010JA015881.Nerger, L., W. Hiller, and J. Schr¨oter (2005), PDAF-the parallel data assimila-tion framework: experiences with Kalman filtering, in Use of high performancecomputing in meteorology: proceedings of the Eleventh ECMWF Workshop on theUse of High Performance Computing in Meteorology, Reading, UK, 25-29 October2004/Eds.: Walter Zwieflhofer; Geoge Mozdzynski, Singapore [ua]: World , pp.63–83.Odstrcil, D. (2003), Modeling 3-D solar wind structure,
Advances in Space Research , (4), 497–506.Odstrcil, D., and V. J. Pizzo (1999), Distortion of the interplanetary magnetic field bythree-dimensional propagation of coronal mass ejections in a structured solar wind, Journal of Geophysical Research: Space Physics , (A12), 28,225–28,239.Odstrcil, D., P. Riley, and X. P. Zhao (2004), Numerical simulation of the 12 May1997 interplanetary CME event, Journal of Geophysical Research: Space Physics , (A2).Owens, M. J., and R. J. Forsyth (2013), The heliospheric magnetic field, Living Re-views in Solar Physics , (5), 5, doi:10.12942/lrsp-2013-5.Owens, M. J., and P. Riley (2017), Probabilistic Solar Wind Forecasting Using LargeEnsembles of Near-Sun Conditions With a Simple One-Dimensional Upwind Scheme, Space Weather , (11), 1461–1474, doi:10.1002/2017SW001679.Owens, M. J., M. Lockwood, and P. Riley (2017), Global solar wind variations overthe last four centuries, Scientific Reports , , 41,548, doi:10.1038/srep41548. –33–onfidential manuscript submitted to Space Weather
Pereira, M. B., L. Berre, M. B. Pereira, and L. Berre (2006), The Use of an EnsembleApproach to Study the Background Error Covariances in a Global NWP Model,
Monthly Weather Review , (9), 2466–2489, doi:10.1175/MWR3189.1.Poedts, S., and J. Pomoell (2017), EUHFORIA: a solar wind and CME evo-lution model, , , 7396.Riley, P., and R. Lionello (2011), Mapping Solar Wind Streams from the Sun to1 AU: A Comparison of Techniques, Solar Physics , (2), 575–592, doi:10.1007/s11207-011-9766-x.Riley, P., J. A. Linker, R. Lionello, and Z. Mikic (2012), Corotating interactionregions during the recent solar minimum: The power and limitations of globalMHD modeling, Journal of Atmospheric and Solar-Terrestrial Physics , , 1–10,doi:10.1016/J.JASTP.2011.12.013.Riley, P., J. A. Linker, and C. N. Arge (2015), On the role played by magnetic expan-sion factor in the prediction of solar wind speed, Space Weather , (3), 154–169,doi:10.1002/2014SW001144.Sasaki, Y. (1970), Numerical Variational Analysis formulated under the constraints asdetermined by longwave equations and a Low-pass filter, Monthly Weather Review , (12), 884–898, doi:10.1175/1520-0493(1970)098 (cid:104) (cid:105) Large Scale InverseProblems : Computational Methods and Applications in the Earth Sciences. , 216pp., De Gruyter.Stone, E. C., A. M. Frandsen, R. A. Mewaldt, E. R. Christian, D. Margolies, J. F.Ormes, and F. Snow (1998), The advanced composition explorer, in
The Ad-vanced Composition Explorer Mission , vol. 86, pp. 1–22, Springer, doi:10.1023/A:1005082526237.Tr´emolet, Y. (2006), Accounting for an imperfect model in 4D-Var,
Quarterly Journalof the Royal Meteorological Society , (621), 2483–2504.van Leeuwen, P. J. (2014), Particle filters for the geosciences, Advanced Data As-similation for Geosciences: Lecture Notes of the Les Houches School of Physics:Special Issue, June 2012 , p. 291.Wang, C., G. Hajj, X. Pi, I. G. Rosen, and B. Wilson (2004), Development of the GlobalAssimilative Ionospheric Model,
Radio Science , (1), doi:10.1029/2002RS002854.Zhu, M., P. J. van Leeuwen, and J. Amezcua (2016), Implicit Equal-Weights ParticleFilter, Quarterly Journal of the Royal Meteorological Society ..