Simulation of Solar Wind Turbulence near Corotating Interaction Regions: Superposed Epoch Analysis of Simulations and Observations

The effect of the turbulence that is associated with solar wind corotating interaction regions (CIRs) on transport of galactic cosmic rays remains an outstanding problem in space science. Observations show that the intensities of the plasma and magnetic fluctuations are enhanced within a CIR. The velocity shear layer between the slow and fast wind embedded in a CIR is thought to be responsible for this enhancement in turbulent energy. We perform physics-based magnetohydrodynamic simulations of the plasma background and turbulent fluctuations in the solar wind dominated by CIRs for radial distances between 0.3 and 5 au. A simple but effective approach is used to incorporate the inner boundary conditions for the solar wind and magnetic field for the periods 2007–2008 and 2017–2018. Legendre coefficients at the source surface obtained from the Wilcox Solar Observatory library are utilized for dynamic reconstructions of the current sheet and the fast and slow streams at the inner boundary. The dynamic inner boundary enables our simulations to generate CIRs that are reasonably comparable with observations near Earth. While the magnetic field structure is reasonably well reproduced, the enhancements in the turbulent energy at the stream interfaces are smaller than observed. A superposed epoch analysis is performed over several CIRs from the simulation and compared to the superposed epoch analysis of the observed CIRs. The results for the turbulent energy and correlation length are used to estimate the diffusion tensor of galactic cosmic rays. The derived diffusion coefficients could be used for more realistic modeling of cosmic rays in a dynamically evolving inner heliosphere.


Introduction
It is well known that the solar wind carries fluctuations in particle phase-space density and magnetic field that are distributed over a broad range of spatial and temporal scales and are commonly referred to as "turbulence." In situ observations by the Mariner (Jokipii & Coleman 1968;Belcher & Davis 1971) and Helios (Bavassano et al. 1982;Roberts et al. 1987) missions provided initial glimpses of the turbulent environment close to the Sun (down to 0.3 au), while the Ulysses mission provided a global view of the solar wind turbulence in high-latitude regions of the heliosphere (Smith et al. 1995;Horbury et al. 1996;Bavassano et al. 2001), and the most recently launched Parker Solar Probe continues to deliver valuable information about the state and spatiotemporal evolution of the turbulence in the outermost regions of the solar corona (Adhikari et al. 2020;Chen et al. 2021;Telloni et al. 2021). Fluctuations are generated on the Sun by convective flows in the photosphere (the so-called pre-existing turbulence) and subsequently in the solar wind, often referred to as the in situ component, through interactions between fast and slow solar wind streams and compression effects caused by shock structures that are also often associated with stream interactions . The ionization of interstellar material and the subsequent isotropization of the resulting pickup ion (PUI) is also an important source of Alfvénic fluctuations more than 5 au from the Sun (Williams et al. 1995). Generally, it is difficult to distinguish between pre-existing and locally generated fluctuations (Horbury & Schmidt 1999;Adhikari et al. 2015), creating a certain degree of ambiguity for the models. Turbulence plays a significant role in various aspects of plasma behavior in space, including solar wind acceleration, solar energetic particle transport, plasma heating, and galactic cosmic-ray (GCR) modulation (Markovskii & Vasquez 2010;Zank et al. 2012;Bruno & Carbone 2013;Engelbrecht & Wolmarans 2020;van den Berg et al. 2020;Chen et al. 2021). It is also thought that turbulent dissipation is responsible for the nonadiabatic radial temperature profile in the solar wind Smith et al. 2001).
The solar wind is comprised of plasma streams with different radial velocities, but it is primarily bimodal (i.e., either "fast" or "slow"). Fast streams originate from open interplanetary magnetic field (IMF) lines in high-latitude regions on the Sun, while slow streams originate from solar equatorial regions. The turbulent properties, such as the fluctuation energy and the correlation length, are also different between the fast and slow solar wind regimes (Bruno & Carbone 2013). Owing to solar rotation, these regimes interact and form more complex structures in terms of the topology of the IMF and plasma features. The most common of these are corotating interaction regions (CIRs). A CIR is generated in the solar wind when a fast stream overtakes the foregoing slow wind (Crooker et al. 1999;Gosling & Pizzo 1999). A shear layer between the fast and slow wind within a CIR is categorized as a tangential discontinuity (TD) and is called a stream interface (SI). Neugebauer (1985) and Neugebauer et al. (1986) reported that velocity fluctuations become aligned with magnetic field fluctuations across TDs, suggesting macroscopic plasma instabilities of the interface to be the cause. Due to the nature of the SIs, both Rayleigh-Taylor and Kelvin-Helmholtz (KH) instabilities could be expected to operate when two streams with different plasma properties slip past each other (Korzhov et al. 1984;Neugebauer et al. 1986;Schwenn & Marsch 1990;Odstrcil 1994). However, the only in situ observation of the KH instability in the solar wind has been reported in Solar Orbiter data by Kieokaew et al. (2021). The observations showed rotation in the direction of the fluctuations of solar wind velocity and magnetic field within a vortical structure, and an enhancement in the fluctuation energy was observed in the power spectrum.
The work presented in this paper can be described as an effort to construct an ab initio model of GCR propagation in turbulent solar wind flows that are dominated by CIRs. The model is based on a 3D magnetohydrodynamic (MHD) simulation of the solar wind plasma and bulk turbulence properties. The simulations were conducted by using partially assumed and partially data-driven physical geometry and boundary conditions to generate CIRs inside the domain, spanning heliocentric distances between 0.3 and 5 au. The entire physical domain can be covered by a single model, because the inner boundary is well outside the Alfvén critical point (∼10 R s ). Of the three sources of turbulence mentioned above, the present model emphasizes the interactions of fast and slow streams, but excludes the PUI source, which is too weak this close to the Sun. The investigation covers two periods near the past two solar minima (2007-2008 and 2017-2018, respectively). These are the same periods that were studied in the observational work of Ghanbari et al. (2019), which serve as a benchmark for the simulations.
Numerous models have been developed over the years to describe the (radial) evolution of the energy range turbulence inside the heliosphere. The simplest WKB model (Barnes 1979) assumes that the turbulence consists of a superposition of linear Alfvénic modes. This model predicts that for an r −2 plasma density radial profile (r being the heliospheric radial distance), the fluctuation energy decreases as r −3 (Belcher & Burchsted 1974;Roberts et al. 1990). The WKB model predictions were generally inconsistent with observations, because the model ignores the processes of wave reflection and dissipation to a turbulent cascade (Hollweg 1974). The next generation of models, based on incompressible MHD (Marsch & Tu 1989;Zhou & Matthaeus 1990a, 1990bMatthaeus et al. 1994;Zank et al. 1996;Matthaeus et al. 1999;Smith et al. 2001Smith et al. , 2006, was much more successful. Observations show that the solar wind generally resembles the behavior of an incompressible magnetofluid (Bruno & Carbone 2013). The higher-frequency fluctuation spectra in the solar wind were found to be power laws that permitted selfsimilar descriptions of the turbulent dissipation process in the energy equations (Kolmogorov 1941;Kraichnan 1965). Nonetheless, it was also recognized that the solar wind contains weak density fluctuations and that the incompressible description was not always appropriate, especially in the low-plasma beta regime. This prompted Zank & Matthaeus (1992) to develop a two-component model, where the leading-order fluctuations are 2D incompressible structures that are distributed in planes normal to the magnetic field, while the secondary component consists of Alfvén waves (the so-called slab component) and (compressible) magnetosonic waves.
Observations have revealed that the solar wind turbulence is often bimodal, consisting of a dominant (∼80%-90%) 2D component, which varies in directions normal to the magnetic field, and a minority (∼10%-20%) slab component, which varies along the field (Bieber et al. 1996). Regardless of which approach is used, the fluid (MHD) description of turbulence has been generally successful in explaining the bulk properties of low-frequency fluctuations in the solar wind (e.g., Goldstein 1995;Tu & Marsch 1995;Matthaeus 2000;Bruno & Carbone 2013;Adhikari et al. 2015;Shiota et al. 2017;Zank et al. 2017;Usmanov et al. 2018).
The present model implements the turbulence evolution equations of Zank et al. (2012), with some modifications, which are described in Section 2. This model contains separate energy equations for the background solar wind and the fluctuating component, using scale separation and Reynolds averaging; the fluctuations are not assumed to be small. According to convention, the fluctuations are described in terms of the Elsässer variables z ± = δu ± δv a (Elsasser 1950), which are sensitive to the direction of propagation of the fluctuations and the alignment between the fluctuating velocity and the Alfvén speed fields δu and d d pr where δB is the fluctuating magnetic field component and ρ is the background flow density), respectively. The starting point is the momentum conservation equation for the fluctuations (Marsch & Mangeney 1987;Usmanov et al. 2011): Here, u and v a are the flow velocity and the Alfvén velocity of the background flow, respectively. Because the focus is on the super-Alfvénic regimes in the solar wind, where v a = u, the terms proportional to the Alfvén speed of the background flow (v a ) are neglected. The resulting momentum conservation equation is simplified to The left-hand side of Equation (2) describes the effects of advection, compression, and shear, while the right-hand side contains nonlinear terms coupling z + and z − together. Note that Equation (2) contains no local (ad hoc) sources of turbulence. As stated above, PUI driving is weak for heliocentric distances under 5 au. Therefore, it is excluded from our model. This model does incorporate pre-existing turbulence amplification by shear (via a mixing of large-scale and short-scale flow variables), but not the production of new fluctuations via ad hoc source terms. This is consistent with the approach of Usmanov et al. (2011) andUsmanov et al. (2016), who also excluded such ad hoc source terms from their turbulence transport model (they did include a PUI driving term because their simulation covered large heliospheric distances). A significant number of 2D and 3D simulations combining the solar wind background and turbulent fluctuation transport have previously been published (e.g., Oughton et al. 2006;Breech et al. 2008;Oughton et al. 2011;Usmanov et al. 2011;Guo & Florinski 2016;Wiengarten et al. 2016;Shiota et al. 2017;Usmanov et al. 2018). Breech et al. (2008) presented the derivation and numerical solution of the three-variable system of turbulence equations for the energy density, correlation length, and normalized cross helicity in the meridional plane (r, θ). They took into account the effect of PUI protons and included driving terms for stream shear interactions and the mixing of the waves. A bimodal solar wind, with slow solar wind below the latitude of 15°and fast solar wind above the latitude of 35°, with a transition region in between, was considered for the background solar wind. The turbulent energy was increased at these latitudes, due to the large gradient of the solar wind speed. The model-derived turbulent quantities were in reasonable agreement with the observations. Wiengarten et al. (2016) developed their two-component turbulence model based on the work by Oughton et al. (2006Oughton et al. ( , 2011. Their model is a 3D treatment of the evolution equations for the high-frequency parallel propagating wave-like and low-frequency perpendicularly cascading quasi-2D components of turbulent fluctuations. They took into account the effect of the mixing of terms involving Alfvén velocity, stream shears, and PUI turbulence driving. They compared the results of their simulation with the results of Breech et al. (2008), Usmanov et al. (2011), andOughton et al. (2011), for the case of a spherically symmetric background solar wind. Their model was also applied to energetic particle transport by modeling the diffusion mean free paths (mfps) and drift length scales. Guo & Florinski (2016) studied the modulation of GCRs by turbulent flows near CIRs. They used a titled bimodal boundary condition and used a three-variable turbulence model, similar to Breech et al. (2008) and Usmanov et al. (2011). In these simulations, the energy difference was constant and equal to one-third of the total turbulent energy. In the results of Guo & Florinski (2016), the turbulent energy did not have an enhancement at the SI within the CIRs, which was inconsistent with some of the observations (Borovsky & Denton 2010;Ghanbari et al. 2019). The effect of dissipation was overestimated by using a large value for the de Kármán-Taylor constant. Shiota et al. (2017) presented the results for multiple simulations with different configurations of the background solar wind, IMF inhomogeneity, and turbulence sources and their effects on the distribution of the turbulence during a solar minimum case. They took into account the effects of turbulence transport perpendicular to the local IMF direction, which has an essential role in determining the latitudinal distribution of the turbulence. In addition to the physical source terms, they included a phenomenological source term on the right-hand side of Equation (1), as a function of radial distance. For the case of a nonaxisymmetric solar wind speed producing a typical CIR pattern, they reported an enhancement of the turbulent energy in the compression regions associated with those CIRs. They proposed that the enhancement of the turbulent energy near the CIRs was produced by the mixing term containing the Alfvèn speed. Specifically, a larger amplitude and gradient of the Alfvèn speed with respect to the flow speed causes the cross helicity increase (decrease), which, in turn, leads to a decrease (increase) of the turbulent energy. They found that regions of velocity convergence also contribute to the enhancement of the turbulent energy. Usmanov et al. (2018) used separate equations for protons and electrons in the region of space from the base of the solar corona (1R e ) up to 5 au, with a tilted dipole-like inner boundary condition. Unlike most of the previous work, they used a finite Alfvén speed and took into account the effects of the Reynolds stress tensor on the solar wind momentum balance. Their results for the solar wind plasma were in good agreement with the Ulysses data, except for the electron temperature.
The model presented here offers several advantages over previous work. The dynamic inner boundary condition is an improvement that enables us to generate CIRs that are reasonably comparable with observations. The turbulence system of equations is solved for six turbulent quantities in three dimensions, which allows us to investigate the variations of each quantity within a CIR. The model does not use a parametric source term to describe in situ turbulence production by shear, relying instead on the enhancement of pre-existing turbulence, according to the transport equation. The turbulence model is described in detail in Section 2, which also specifies the geometry of the simulation and the boundary conditions. The key results of this work are reported in Section 3, while their context and significance are presented in Section 4.

Framework
The underlying numerical framework is based on a 3D geodesic mesh, consisting of a 2D hexagonal tessellation extruded radially in a concentric fashion, with a variable stepping. A detailed description of the framework can be found in Florinski et al. (2013). The background component of the model solves the standard system of MHD conservation laws for the flow density, momentum and energy, and magnetic flux. Control of the magnetic field divergence is accomplished through the use of the generalized Lagrange multiplier method of Dedner et al. (2002). The system of MHD equations is solved with a finite-volume second-order (space and time) method on hexagonal prisms (Florinski et al. 2013). Reconstruction is employed to achieve second-order accuracy, using 1D and 2D versions of the weighted essentially nonoscillatory limiter, in both radial and tangential directions. Fluxes are computed using the HLLC Riemann solver (Li 2005). This solver consists of four states, including the left and right unperturbed states and two intermediate states that are separated by a TD.
The simulation domain was restricted to supersonic solar wind between 0.3 and 5 au, which greatly simplifies the treatment of both inner and outer boundary conditions. Two sets of simulations were conducted for this paper. The first set includes the results of a run using a level 5 geodesic grid with 10,242 hexagonal cells on the surface of a sphere, with 256 spherical layers in the radial direction (low resolution), for comparison of the simulation results with the observations at 1 au on Earth's orbit. The second set includes a level 6 geodesic grid with 40,962 hexagonal cells on the surface of a sphere, with 512 spherical layers in the radial direction (high resolution), to investigate the results in the meridional plane as well as a spherical slice at 1 au. The results for both sets are presented in Section 3.
The Parker (1958) magnetic field solution was imposed at the inner boundary. Two cases of monopole and split monopole fields were tested in this simulation. In the case of a monopole field, a passive variable was defined to track the magnetic field polarity, by being advected by the plasma flow to solve Equations (5)-(9). In the case of a split monopole field, a sign change was applied to the components of the magnetic field, using the same passive variable to define a neutral current sheet in the simulation. The results of both cases are compared for their plasma properties in Section 3. Fast and slow streams were introduced at the inflow boundary in a time-dependent manner, leading to the formation of CIRs, as discussed in 2.3.

Turbulence Transport Equations
The turbulence model is based on the six-equation formalism of Zank et al. (2012), which consists of equations for the densities of energy-like quantities and the correlation lengths. The energy variables are the energy sum E t , the cross helicity H c , and the energy difference E d , defined as In the above, l ± are the correlation lengths of the fluctuations that are associated with the two Elsässer variables, and l d is the correlation length corresponding to the energy difference variable.
There is no feedback from the turbulence on the background flow. The transport equations themselves are derived by multiplying Equation (2) by z ± , evaluated at a lagged spatial position, and performing ensemble averaging. The nonlinear terms from Equation (2) turn into triple correlations in the energy equations, which are based on some prescribed shortscale behavior. We use the turbulence transport equations in a quasi-conservative form, where the transported quantities have the units of energy density and the fluxes describe advection with the solar wind; the remaining terms describing the coupling with the large-scale flows are treated as source terms. The forms of the turbulence transport equations that are implemented in the model are slightly different from those of Zank et al. (2012). We use only "extensive" variables, which correspond to the energy densities and integrals defined in Equations (3) and (4), times the plasma density, ρE t , ρH c , ρE d , ρL ± , and ρL d . The evolution equations are given as . An important coefficient in these equations is the so-called de Kármán-Taylor constant α, which is generally thought to be of order one (e.g., de Kármán & Howarth 1938;Breech et al. 2008;Chhiber et al. 2021). However, in this model α = 0.05, because larger values resulted in correlation lengths that were an order of magnitude larger than expected, based on the observations, in agreement with the same conclusion reached by Usmanov et al. (2018). The nonlinear terms are all modeled as in Zank et al. (2012), except for the energy difference, which is based on Zank et al. (2017). This is because the dissipation term from the former paper could generate unphysical solutions, as was demonstrated in Dosch et al. (2013).

Boundary Conditions
A common method of generating CIRs in MHD simulations is to define the trace of the neutral line for the magnetic field at the inner boundary of the simulation with a great circle tilted with respect to the rotation axis (Pogorelov et al. 2007). The slow solar wind conditions are then applied at latitudes less than a given angle from the current sheet (in the tilted frame), while the fast-wind conditions are imposed at higher latitudes in both hemispheres. As the inner boundary rotates, a CIR structure develops, which eventually becomes fully periodic in space and time. Observed CIRs, however, lack this periodicity, because the magnetic field topology is rarely that simple. To obtain a more realistic CIR topology that could be compared with observations, then, a simple but effective boundary condition model is developed. A standard synoptic magnetic field map is not directly usable, because it is not periodic in longitude, owing to the time evolution of the field. However, the Wilcox Solar Observatory 1 also provides coefficients of the Legendre expansion for the potential field source surface (PFSS) model for all recorded Carrington rotations (CRs) from 1976 until the present time, up to the ninth-degree harmonic. By constructing surface magnetic field maps from the Legendre expansion at (R = 2.5R s ), the current sheet is identified as the location where the polarity of the magnetic field switches sign. Some maps have isolated "islands," which represent regions of a stronger magnetic field surrounded by weaker-field regions. This enables the simulation to model a fast-wind region that is entirely surrounded by a slow-wind region. This feature is not available in previous methods for distinguishing the slow and fast wind by using a high-order polynomial trace of the heliospheric current sheet (HCS; e.g., Guo & Florinski 2016;Usmanov et al. 2018).
The maps are normalized in magnitude, and then the regions with a normalized factor of less than 0.09 are set to be slow wind, with the remaining regions being set to have the properties of the fast wind. Since this factor is a threshold for determining the slow-fast transition line, we call it the slowfast threshold factor. Figures 1(a) and (b) show two reconstructed periodic maps for CRs 2051 and 2052. In these figures, the y-axis represents the latitude on the Sun (or the inner boundary of our simulation). The gray shaded area is the slow wind, while the remaining areas are considered as the fast solar wind. As expected, the north pole has a negative field polarity (blue) and the south pole has a positive polarity (red), due to the negative polarity of the Sun during the 2007-2008 solar minimum. The black line in the middle of the shaded area is the projection of the current sheet onto the inner boundary. An "island" is also visible in CR 2052, as shown in Figure 1(b). In this context, an "island" means a region that is specified as fast stream being surrounded by regions that are defined as slow streams. Linear interpolation between two consecutive CRs is performed in order to obtain the field at an arbitrary moment in time. The central meridian correction is also applied at each time step. The boundary values for the solar wind plasma background and turbulence are given in Table 1. The magnitude of the magnetic field is taken to be 2.4 × 10 −5 G at 1 au in the ecliptic slow wind. The components of the magnetic field vector are obtained from the Parker model, where the radial component of the magnetic field is obtained from Gauss's law and the azimuthal component is obtained from Faraday's law. The solar wind velocity vector is assumed to be purely radial at the inner boundary of the simulation. Note that the energy difference shown in Table 1 is negative, as expected from the dominance of the magnetic energy over the kinetic energy, which is characteristic of the solar wind turbulence.

Results
1 minute in situ spacecraft data obtained from the Omniweb interface 2 were used for comparison with the simulation results. One of our goals was to evaluate the effects of the HCS on the turbulence evolution within the CIRs. To isolate the physics specific to the HCS, we first ran the simulation with a monopole magnetic field, in which the field lines were directed outward everywhere. In the split monopole case, the magnetic field components were reversed (i.e., became inward) wherever the polarity factor was negative. This factor was determined at the inner boundary, using the observed polarity of the IMF based on the reconstructed maps of B r . The magnitude and the r and f components of the IMF at 1 au near Earth's orbit in the simulation are shown in the top and bottom panels of Figure 2, respectively. The monopole case is shown with the solid lines and the split monopole case is shown with the dashed lines. The magnitudes of the field (|B|) in these two cases differ from one another in the sharp dips that are present in the split monopole configuration. This effect is caused by the fact that in HCS crossings, all the components of the magnetic field pass through zero. These sharp dips are absent from the monopole configuration, as expected, and instead the height of the peak in |B| is larger. The turbulent properties in the simulation did not change appreciably between the two configurations for the magnetic field. Because the monopole case produced slightly higher turbulent energies inside the CIRs, we have used that model for the rest of this paper.
The time profiles of both the plasma and the turbulent properties at Earth's orbit obtained from the low-resolution runs are presented in Figures 3, 4, 5, and 6 with the simulation results shown in blue and the observations shown in red.   shows the magnetic field polarity factor. The positive values of this quantity represent the outward IMF and vice versa. Since this quantity is part of the large-scale structure of the solar wind, it is expected to be well reproduced (Hoeksema & Scherrer 1986). Some of the HCS crossings are within a day of the observed HCS crossings. The away-to-toward (+ to -) crossings occurred earlier in the simulation than in the observations, at days 15, 26, 43, 53, 127, 135, 138, 164, and 284. The azimuthal expansion of the fast streams by radial distance, especially at the trailing edges of the CIRs, could be amplified by the factors excluded from this simulation. In reality, the transition region between the fast and slow streams can be affected by transient events on the way from the Sun to the Earth. This could possibly cause a later occurrence of the away-to-toward HCS crossings in the observation. The solar wind speed is shown in the second panel from the top. The observed CIRs are marked with the passage of their SIs by the vertical brown dotted lines in this panel and the panel below. The large-scale quasiperiodic structures (CIRs with their respective HCS crossings) are generally reproduced to within a few hours to a day of their appearance in the observations. The widths of the fast streams are different in different CIRs in both the simulations and the observations. A better agreement between the two can be achieved by setting a dynamically variable (with time) slow-fast threshold factor at the inner boundary of the simulation, to set the width of the slow-stream band. Some events have clearly not been reproduced in the simulation, such as those at days 170, 180, 200, 222, 227, 249, 257, and 276. In these CIRs, the fast streams had smaller speeds than those that are reproduced. These slower CIRs could have been affected by transient events in the solar wind that were not included in the simulation.
The third panel of Figure 3 illustrates the IMF magnitude. Due to the highly fluctuating nature of the IMF, only the IMF peaks associated with the CIRs are reproduced, and the fluctuations are not present in the simulation results. There is a peak in the IMF magnitude at each CIR, which is a typical signature of an SI. Finally, the proton density is shown in the bottom panel. This quantity is also reproduced reasonably near the CIRs with a high-density slow wind just before a peak, followed by a fast wind with lower density after the peak. The magnitudes and trends of the simulated quantities are close to the observed values. Overall, the magnitudes and heights of the peaks of the density are comparable between the simulations and the observations. Figure 4 shows the results for the first 300 days of the year 2017 (near the positive cycle minimum). The top panel shows the magnetic field polarity change over this time interval. There are more missed crossings in the IMF polarity in this period than in the previous period. This is mostly due to the variations in the shape of the current sheet and the abundance of magnetic islands during CRs 2188 and 2189. As a result, rapid IMF polarity switches are observed on Earth, as seen in the top panel of Figure 4, between days 80 and 140. These fluctuations are not reconstructed in the model output. This has also affected the simulation results of the solar wind speed plotted in the second panel of Figure 4. The magnetic field magnitude is shown in the third panel. One can see that the observations have sharp features that are absent from the model results. The peaks in the IMF are evident at the locations of the SIs in the simulated CIRs. The fourth panel shows the proton number density. The smooth increase in the density from the fast to slow streams, the jump right before the SI, and the subsequent drop in this quantity at the SI of the simulated CIRs are evident. In summary, the simulation reproduces more CIRs during the first period than during the second period. Figure 5 shows the observed and simulated turbulent quantities, along with the diffusion coefficients and mfps, at the Earth's location for the same time interval as in Figure 3. The observed turbulent quantities were computed over 6 hr intervals, with 2 hr steps. The observed turbulent energy E t , shown in the top panel in red, is the sum of the variances of the solar wind speed and magnetic field fluctuations obtained for each 6 hr interval. The values in the slow and fast streams are determined by the respective boundary conditions in each region. While they are qualitatively in agreement, the turbulent energy in the fast streams is slightly greater in the model compared with the data. The simulation could not properly resolve the sharp increases in E t at the SIs, which was probably a consequence of insufficient grid resolution. In this period, the fast-stream turbulent energy is about 1800 km 2 s −2 , while the slow-stream turbulent energy is around 200 km 2 s −2 . The heights of the peaks of E t at the SIs are variable for the different CIRs, but on average they reach up to 2500 km 2 s −2 at the highest point.
The correlation length (l t = L t /E t ; second panel from the top) was inferred from the bendover length of the measured power spectrum of the magnetic field fluctuations for each interval, using the technique of Ghanbari et al. (2019). The observations and simulation results for this quantity have the same order of magnitude. This quantity fluctuates between 0.008 and 0.013 au in the slow and fast streams, respectively. The l t computed from the data changes very rapidly in the observations. Unlike the energy, there are no sharp peaks in this quantity at the SIs in the simulation results, as in the E t . Some small sudden increases are seen at the leading and trailing edges of the CIRs that are small in magnitude.
The third panel from the top in Figure 5 shows the normalized cross helicity, defined as the covariance of the wind velocity and magnetic field fluctuations divided by the turbulent energy, σ c = H c /E t . It is reasonably well reproduced in the simulation. This quantity varies between 1 and −1, where the former implies the dominance of the parallel propagating waves, and vice versa. Since this quantity depends on the field polarity, it is also an alias for the latter. Some significant differences between the simulation results and the observations in this panel occur between days 250 and 260, as well as days 270 and 290. In these two time intervals, the growths of the forward waves are gradual in the observations, whereas the simulation results show a sudden change from parallel to antiparallel waves at days 253 and 280.
The parallel (κ ∥ ) and perpendicular (κ ⊥ ) diffusion coefficients are computed according to quasilinear theory (Jokipii 1966) and the nonlinear guiding center theory (Matthaeus et al. 2003;Zank et al. 2004) formulations, respectively. In brief, the parallel mfp is inversely proportional to the power spectral density of the magnetic fluctuations resonant with the particle, while the perpendicular mfp, in the high-rigidity limit, is proportional to the cubic root of the parallel mfp and the total magnetic variance of the fluctuations raised to the power of 2/ 3. The turbulent energy and correlation/bendover length obtained from our simulation were used to compute these diffusion quantities. The perpendicular diffusion coefficient is shown on the left y-axis, in black. It is typically in the range 1-12 × 10 20 cm 2 s −1 , and is higher in fast streams and lower in slow streams. This quantity has a small sharp drop right at the SI, which is caused by the small sudden decrease in the correlation length. It also increases toward the trailing edges of the CIRs, and has a peak when the fast stream comes after the slow stream. The parallel diffusion coefficient is shown on the right y-axis, in magenta. This quantity takes the values of 2-9 × 10 22 cm 2 s −1 , with lower values in fast streams and higher values in slow streams. This behavior is opposite to that of the perpendicular diffusion coefficient. Instead of a small sharp dip, it has a sharp peak at the SIs. The larger peak in κ ∥ compared to the smaller dip of κ ⊥ at the SI could be due to the fact that k μ l t 2 3 , while k µ - In the bottom panel, the parallel and perpendicular mfps are plotted. The perpendicular mfp is shown on the left y-axis, in black, and it fluctuates between 0.002 and 0.009 au. The parallel mfp is shown on the right y-axis, in magenta. This quantity varies between 0.1 and 0.6 au. Since these quantities are linearly dependent on the corresponding diffusion coefficients, they follow the same trend and have the same features as the diffusion coefficients. Figure 6 shows the same quantities as in Figure 5, but for the first 300 days of 2017. The first panel shows the turbulent energy in both the observations (blue) and simulations (red). The simulated turbulent energy in the fast streams is about 1800 km 2 s −2 , which is slightly higher than the fast-stream energy observed during this period. The heights of the peaks in the turbulent energy obtained from the observations reach higher than 4000 km 2 s −2 , but are not reproduced in the simulation. The second panel from the top illustrates the correlation length, which fluctuated between 0.007 and 0.013 au. The highest observed correlation length in this period is about 0.03 au, which is never reached in the simulation. Since this quantity is computed indirectly from L + and L − , it is possible to obtain better estimates for it by using a slightly different boundary condition for the L +/− . For the same reason as mentioned in the discussion following Figure 4, the normalized cross helicity shows a very perturbed cycle during this period. [70][71][72][73][74][75][76][77][78][79][80][150][151][152][153][154][155][156][157][158][159][160][161][162][163][164][165][180][181][182][183][184][185][186][187][188][189][190] are clearly very perturbed for this quantity, based on the observations. These perturbations are not reproduced in the simulation.
The parallel (right; magenta) and perpendicular (left; black) diffusion coefficients for this period are shown in the fourth panel from the top. κ ⊥ varies from 3 × 10 20 cm 2 s −1 in slow streams to 10 × 10 20 cm 2 s −1 in fast streams. Following the same trend as the fourth panel from the top in Figure 5, it increases within the fast stream, until the trailing edge of the CIR, where it has a sharp peak, before dropping to the slowstream values. κ ∥ varies between 2 and 8 × 10 22 cm 2 s −1 in the fast and slow streams, respectively. The sharp peaks of this quantity are found at the same times as the dips in κ ⊥ , right at the leading edges of the CIRs. The bottom panel shows the parallel and perpendicular mfps. The trend for this quantity is the same as for the diffusion coefficients. The high and low values of λ ⊥ are 0.002 and 0.007 au, while those for λ ∥ are 0.2 and 0.7 au, respectively. Figure 7 shows the 3D features of the simulated structures in the solar wind background plasma at day 95 of year 2007. In this figure, the top panels present the meridional cuts, while the bottom panels show spherical slices at r = 1 au. The left panels depict the solar wind speed (V sw ), the middle panels show the magnitude of the magnetic field (|B|), and the right panels show the proton density (ρ).
Panel (a) shows two slow streams (blue) originating from the equatorial regions; the one on the left is trailed by a fast stream (red). Their interaction produced a CIR. V sw is about 700 km s −1 in the fastest streams and 250 km s −1 in the slowest stream in this snapshot. There is no smooth transition between the slow and fast streams in the latitudinal direction at the inner boundary of the simulation. Larger gradients in velocity can contribute to the enhancement of the turbulent energy at the SIs within the CIRs. Panel (d) shows a spherical slice of V sw at 1 au. A compression region and a rarefaction region are marked with the black and green lines, respectively. A weak broadening of the slow-fast stream transition region is generated due to the numerical effects at 1 au, which can be seen in panel (d).
The contoured regions have higher (compressed) and lower (rarefied) magnetic field magnitudes compared to the neighboring regions, as shown in panel (e). Since this snapshot is obtained from the split monopole field configuration, one can see a narrow line in |B| at the trace of the current sheet, where |B| falls to ∼3 × 10 −6 G. This dip was absent from the monopole field case. The meridional view of the magnetic field in panel (b) shows the complex topology of the magnetic field near the equatorial regions, in comparison to the polar regions, where |B| decreases smoothly with radial distance.
Panel (c) shows the plasma density, which is about 100 cm −3 near the inner boundary and 0.1 cm −3 near the outer boundary at high latitudes. The meridional view of the density in panel (c) also shows how the high-density slow wind interacts with the low-density fast wind at low latitudes, while the polar regions have a smooth decrease in density outward in the radial direction. Panel (f) shows a spherical slice of the density profile at 1 au. One can see in this panel that the peak of the density is within the slow stream, close to the SI, in the region shown with the black contour in panel (d). A rarefaction and, consequently, a lower-density region is visible in panel (f), corresponding to the green contour in panel (d). These variations in the plasma features are consistent with the characteristic signatures of CIRs that are seen in the observations (Jian et al. 2011). Figure 8 shows the turbulent quantities in the meridional plane (top three panels) and on the surface of a sphere with a diameter of 2 au (bottom three panels). The turbulent energy is shown on the right, the correlation length in the middle, and the normalized cross helicity on the left. The normalized cross helicity σ c = H t /E t is shown in panels (a) and (d). Since the latter depends on the polarity of the magnetic field, this variable also traces the HCS that is clearly visible as the line separating the positive σ c (red and yellow) regions from the negative σ c (blue and green) regions. Note that due to the variable boundary conditions in the simulation, the shape of the current sheet also changes dramatically with time.
The meridional variations of the correlation length l t are shown in panel (b). This quantity is approximately 0.01 au at 1 au and 0.06 au at the outer boundary of the simulation in the polar regions. It increases with radial distance, which is in agreement with the observations. From panel (e), one can see that the correlation length increases right at the SI and wherever there is a shear in the flow. This quantity has a value of 0.02 au in the fast streams and 0.008 au in the slow solar wind regions.
The magnitude of the turbulent energy is presented in meridional (panel (c)) and spherical (panel (f)) views in Figure 8. According to the boundary conditions, the turbulent energy is only slightly higher in the fast wind than in the slow wind, at 0.3 au. The energy decreases with increasing radial distance, but the effect is much stronger in the slow wind, which has a higher rate of turbulent dissipation than the fast wind, so by 1 au the two differ by about a factor of 10. The turbulent energy is large in the compression region that is associated with the CIR, as shown with the area of the black contour in panel (f). It is lower in the rarefaction region, which is shown with the green contour in the same panel. Some ripples are generated at the interface between the slow and fast regions, as shown with the red contour, resembling the features of KH instability. However, because the length scales in the model are larger than the characteristic length scale of KH instabilities, these ripples are most likely artificial effects that have been generated through the visualization process. No such ripples are generated within the CIR shown with the black contour. The ripples are more prominent in panel (e), especially in the same regions that are shown with the red contour in panel (f). The highest-fluctuation energy is reached on the fast-stream side of the SI (3000 km 2 s −2 at 1 au).
We next performed superposed epoch (SPE) analysis of the turbulent properties of both the spacecraft data and our simulations. In this method, multiple time records of the quantities of interest are aligned with the zero epoch (a reference point representing a key event) and stacked on top of one another, to produce an average. In this case, the reference point is the SI passage. The SPE analysis of 55 CIRs in the simulation results for the first period (2007)(2008) were compared with the SPE analysis of 90 CIRs that were observed during the same period. The results, for an 8 day interval around the SI, are shown in Figure 9.
The SPE analysis results for the turbulent energy are shown in the top left panel. The mean trend is well reproduced in both the slow-and fast-stream sides of the SI, as well as the peak at the SI. Although, as was seen in the time profile results, the heights of the peaks for this quantity are significantly lower than the observed values, the two are much closer to one another in the SPE analysis. This, of course, is a consequence of the averaging that is performed by the SPE analysis, which tends to eliminate small-scale features. The top right panel shows the SPE analysis of the correlation length computed from the observations and the model. The length is shorter in the slow wind, because it is more evolved, and the turbulent cascade has developed a longer inertial range (Borovsky & Denton 2010). The simulated SPE values are about 10% smaller on both sides of the SI, compared to the SPE analysis of the observations.
The bottom left panel of Figure 9 compares the parallel diffusion coefficients as derived from the simulation results (solid lines) and the data (dashed lines) and subjected to SPE averaging for protons with energies of 100 MeV, 200 MeV, 500 MeV, and 1 GeV. For all energies, κ ∥ has its maximum nearly a day before the SI. This maximum corresponds to the minimum of the correlation length, as shown in the top right panel. The peaks are due to the inverse dependence of κ ∥ on the correlation length. The bottom right panel illustrates the perpendicular diffusion coefficient computed from the SPE analysis of the simulations (solid) and the observations (dashed) for the same energies as used for κ ∥ . Although the trend for this quantity is followed correctly, its magnitude is 10%-20% smaller in the results obtained from the simulations as compared to the observations, on the slow-stream side of the SI. This difference increases to 50% at the SI, then decreases to less than 10% on the fast-stream side of the SI. There is a shallow dip about a day before the SI, which occurs at the same time as the minimum of the correlation length. This feature is also a consequence of the direct dependence of κ ⊥ on the correlation length.
The same analysis was performed on 75 observed CIRs during the second period (2017-2018), compared with the SPE analysis of 45 modeled CIRs, as shown in Figure 10. The top left panel shows the turbulent energy. The mean energy in the modeled CIRs (blue) is slightly lower on the slow-wind side of the SI. The mean turbulent energy is lower in the SPE analysis of the observed CIRs on the fast-stream side of the SI. The height of the peak in the turbulent energy is almost the same for the SPE analyses of the simulations and observations, half a day after the SI. The top right panel shows the SPE analyses of the correlation length for simulated (blue) and observed (red) CIRs. The same trend as for the first period is also seen for the second period. The mean correlation length is approximately 30% higher in the SPE analysis of observations on the slowstream side of the SIs. This quantity rises almost half a day earlier in the SPE analysis of the simulated CIRs than in the observed CIRs, and unlike the observations it has its peak right at the SI. The SPE analysis of the observed CIRs reaches its peak at around two days after the SI in the fast stream. The mean correlation length is roughly 20% higher in the SPE analysis of observations. The parallel diffusion coefficients obtained from our model (solid lines) and from the observations (dashed lines) are illustrated for protons with same energies as Figure 9 in the bottom left panel of Figure 10. In all energies, κ || has its maximum nearly a day before the SI, which coincides with the minimum of the correlation length. The modeled values of this quantity have slightly different behaviors in the second period, compared to the first period, which are due to the marginally different trends of the turbulent energies and correlation lengths in these two periods. The bottom right panel compares the perpendicular diffusion coefficients for protons with the same energies. The mean values of κ ⊥ obtained from both the observations and the model decrease before the SI, until half a day before the SI. In all energies, the minimum is half a day before the SI. The values then start rising, with the maximum being reached two days after the SI in the fast stream. The observed values are 50% higher than the values derived from the model results, before the SI. The separation between these two decreases to around 20%, for 1 GeV protons, and 10%, in the case of 100 MeV protons.

Discussion and Conclusions
In this paper, we have presented and demonstrated the capabilities of a new computer model of turbulence evolution in the solar wind, during the periods of low solar activity that are dominated by a recurrent CIR structure. It has been shown that the method of determining the magnetic field polarity, the slow-fast stream boundaries, and the corresponding background solar wind and turbulence parameters, based on source surface magnetic field maps, can be effective in modeling the large-scale structure of the solar wind near the 2007-2008 negative solar minimum. On the other hand, the presence of solar transient events at the beginning of the most recent positive solar minimum (2017), have caused significant discrepancies between the simulation results and the data.
The ultimate purpose of this model is to enable physicsbased modeling of quiet-time cosmic-ray transport in the heliosphere in three dimensions, driven by observational inputs. While the only data incorporated in the present version involve the shape of the solar magnetic equator, they produce a rich evolving structure of magnetic sectors and stream boundaries that are characteristic of actual CIRs. GCR transport critically depends on the turbulent properties of the medium. By 1 au, the turbulent quantities have a bimodal structure, where the fluctuations are on average weaker and smaller in size in the more evolved slow wind, compared with the fast wind, which has a much lower rate of turbulent decay with heliocentric distance. The interaction between the two regions produces the dominant time-dependent pattern of weak/smallstrong/large fluctuations. As a rule, the particle scattering rate is higher in regions with stronger fluctuations, but an increasing correlation length has the opposite effect, because more fluctuations are now out of the resonance interaction range with GCRs, with energies between a few hundred MeV and a few GeV. As a result, the changes in the parallel diffusion coefficient across the CIR structure are somewhat damped. The perpendicular mfp has a significantly stronger response, varying by as much as a factor of 4 or more across the SI. The most important feature is the minimum that typically precedes the SI passage by a few hours. There is a strong correlation between this phenomenon and the relatively large decrease in GCR intensity between the slow and the fast sides of the SI (Ghanbari et al. 2019). The simulations presented here feature very significant drops in κ ⊥ near the SIs, and the authors believe that GCR transport models based on these results could confirm the cause and effect relationship between the two trends.
Another effect worthy of note is the increase in the amplitude of the cosmic-ray modulation by the CIR structure that occurred in the second half of 2007 and persisted through early 2008 (Leske et al. 2011). This is evident in the neutron monitor and the ACE Cosmic Ray Isotope Spectrometer (CRIS) measurements reflecting GCR intensities, but also in the ACESolar Isotope Spectrometer (SIS)data that record the lower-energy anomalous cosmic rays. Inspecting Figure 5 shows that there is a corresponding change in the diffusion pattern that started around day 180. Where, at earlier times, we observe a succession of two magnetic sectors of roughly equal duration per rotation, in the latter period, the two sectors are very close together, separated by a long stretch of quiet slow solar wind. As a result, instead of a 13 day low-amplitude periodicity, cosmic-ray intensities have a 27 day highamplitude periodicity during the latter period. It should be possible to reproduce this phenomenon with cosmic-ray transport models that are based on the diffusion input from the present model.
Any model that is based on a series of synoptic maps requires time interpolation between consecutive maps, to obtain the boundary condition at an arbitrary time. This could potentially create a time shift between the observed feature and its simulated counterpart. The reconstructed maps are most accurate at the central meridian and least accurate on the far side of the Sun, which is also a source of error in the model. In particular, streams that are generated near the central meridian are reproduced with the highest fidelity, but the accuracy decreases by moving away (both eastward and westward) from that point. In some cases, this lower accuracy results in a significant loss of information at the inner boundary. In some cases, a CIR is not produced at all, while it is generated in others, but with a time shift relative to the observational data. Transient and short-timescale events, such as coronal mass ejections, are not present in this model.
Although the main purpose of this work was not to reproduce the solar wind during the given periods, a statistical analysis of the simulation results has shown that 60% of the HCS crossings were simulated within a few hours of the times at which they were observed, while the remaining 40% were within a day of the respective observed crossing times during the first period. Nearly 60% of the HCS crossings of the second period were missed by the simulations. Here, 45% of the CIRs were simulated within a day of their observation, about 20% were observed less than 2 days from their observations, and about 35% were absent from the simulation results. The CIRs that were not reproduced had a smaller difference in the speeds of the slow and the fast streams than those that were reproduced (ΔV ∼ 150 km s −1 ). The model could be improved by using a time-variable slow-fast threshold factor for the normalized magnetic field maps. In the current implementation, this threshold is constant and equal to 0.09. Because the solar magnetic field on the Sun is changing with time, latitude, and longitude, a time-varying slow-fast threshold might result in more accurate predictions of HCS and CIR arrivals.
A more accurate boundary condition for the turbulent quantities could also improve the results of our simulations. Specifically, the correlation lengths for the forward and backward Elsässer energies are poorly known at the heliocentric distance corresponding to the inner boundary in our model. The boundary conditions for the two integral quantities that are used in the current simulations were obtained experimentally, by testing a range of values to find the set giving the best estimates of the turbulent energy and correlation length (those closest to the observed values of these quantities) near 1 au. An obvious improvement would be to use boundary conditions with more observational inputs. This would result in more realistic estimates of the correlation length and, consequently, the diffusion coefficients within CIRs.
We have neglected the effects of the finite Alfvén speed in the turbulent transport component of our model. As discussed in Shiota et al. (2017), the mixing terms involving this variable could possibly amplify the increase in the turbulent energy that is seen within the simulated CIRs. The increase in the turbulent energy in our simulations is about 15%, and is right behind the SI in the fast wind within the CIRs. While this is consistent with the results of the SPE analysis in Ghanbari et al. (2019; see also Borovsky & Denton 2010;Shiota et al. 2017), the enhancements that are measured in individual events are much larger. Test simulations that were performed for short time periods showed that the heights of the peaks increase with the increasing spatial resolution of the numerical mesh. However, using a geodesic mesh with division higher than seven turned out to be impractical (numerically expensive) for simulations covering periods of the order of 1 yr in duration. Even so, this work has shown an improvement in comparison with Guo & Florinski (2016), where no enhancement of the turbulent energy was seen at the SIs. The same applies to other turbulent quantities, such as the correlation length. It would be fair to say that the presented results reproduce the average trends of the turbulent quantities observed in the solar wind at 1 au.
One of the important quantities governing the evolution of turbulence with radial distance is the de Kármán-Taylor constant that enters into Equations (5)-(7) (de Kármán & Howarth 1938;Matthaeus et al. 1996;Zank et al. 2012;Usmanov et al. 2018). While the value of this constant is close to one in some work (Hossain et al. 1995;Matthaeus et al. 1996;Usmanov et al. 2011;Dosch et al. 2013), Usmanov et al. (2018 adopted a much smaller value of 0.128. We have found that making this constant even smaller (0.05) gives the best agreement with the 1 au observed quantities. More studies involving observations at different radial distances would be required to narrow down the possible range of this parameter. Finally, our model does not include local sources of turbulence, i.e., all turbulence is injected at the inner boundary. Turbulence is modified in regions with a high degree of flow velocity shear (such as near SIs), as a result of energy exchange with the plasma background, but this is not the same as injecting new fluctuations from some sort of instability. In particular, we note that our model does not have sufficient resolution to properly simulate the growth of the KH instability at fast-slow SIs. A quantitative investigation of the occurrence of the KH instability and the generation of the vortices in the shear layer would require an extremely high spatiotemporal Figure 10. The same as Figure 9, but for the second period (2017-2018).