A multiple flux-tube solar wind model

We present a new model, MULTI-VP, that computes the three-dimensional structure of the solar wind which includes the chromosphere, the transition region, and the corona and low heliosphere. MULTI-VP calculates a large ensemble of wind profiles flowing along open magnetic field-lines which sample the whole three-dimensional atmosphere or, alternatively, on a given region of interest. The radial domain starts from the photosphere and extends, typically, to about 30 $R_{sun}$ . The elementary uni-dimensional wind solutions are based on a mature numerical scheme which was adapted in order to accept any flux-tube geometry. We discuss here the first results obtained with this model. We use Potential Field Source-Surface (PFSS) extrapolations of magnetograms from the Wilcox Solar Observatory (WSO) to determine the structure of the background magnetic field. Our results support the hypothesis that the geometry of the magnetic flux-tubes in the lower corona controls the distribution of slow and fast wind flows. The inverse correlation between density and speed faraway from the Sun is a global effect resulting from small readjustments of the flux-tube cross-sections in the high corona (necessary to achieve global pressure balance and a uniform open flux distribution). In comparison to current global MHD models, MULTI-VP performs much faster and does not suffer from spurious cross-field diffusion effects. We show that MULTI-VP has the capability to predict correctly the dynamical and thermal properties of the background solar wind (wind speed, density, temperature, magnetic field amplitude and other derived quantities) and to approach real-time operation requirements.


INTRODUCTION
The solar wind separates into fast and slow wind streams whose large-scale distribution evolves markedly during the solar activity cycle. Slow solar wind flows are usually found in the vicinity of streamer / coronal hole boundaries (S/CH), while the fastest wind flows stream out of coronal holes. The distribution of fast and slow wind streams hence follows the cyclic variations of the magnetic field structure of the solar corona (Mc-Comas et al. 2008;Smith 2011;Richardson & Kasper 2008). In fact, slow wind flows are essentially confined to the equatorial regions of the Sun during solar minima, with the fast wind flows covering all other latitudes. During solar maxima, slow wind flows also occur at higher latitudes following the appearance of nonequatorial streamers and pseudo-streamers. During the rise and decay phase of the cycle, fast wind flows also make incursions into the equatorial regions (McComas et al. 2003). The terminal wind speeds seem to be determined to a great extent by the geometrical properties of the magnetic flux-tubes through which the solar wind flows (Wang & Sheeley 1990;Pinto et al. 2016). In fact, the properties of the surface motions (assumed as energy sources for the heating and acceleration of the wind) are more uniform across the solar disk than the variations in amplitude of the wind flows that stream out, suggesting that it is the coronal environment that causes the segregation between fast and slow solar wind flows. Other than their asymptotic speed, slow and fast wind streams are differ consistently in terms of density, mass flux, heavy ion composition and thermal structure.
Understanding solar wind acceleration and predicting the terminal solar wind speed (together with other properties) has been the subject of intense research over the last decades. The solar wind speed is commonly associated with simple parameters describing the variations of the cross-sections of the flux-tubes as a function of height, namely the expansion factor where A 0 and A 1 are the cross-section of a given elemental flux-tube respectively at the surface of the Illustration of the operation of the model. The first row corresponds to CR 2056 (April -May 2007, close to solar minimum) and the second one to CR 2121 (Mars 2012, close to solar maximum). The first column shows the input WSO magnetogram rendered in gray-scale over the surface of the Sun and a sample of the magnetic field lines obtained via PFSS extrapolation used to initiate the model. The transparent yellow surface indicates the coronal hole boundaries (the closed-field regions are excluded from the domain). The second column shows a close-up of the wind speeds in the low corona, represented in colour-scale (from dark blue at 250 km/s to dark red at 650 km/s). The third column shows the same information at global scale (truncated at a radius of 15 R and with one octant removed).
Sun (r = r 0 ) and at some point higher in the corona (r = r 1 > r 0 ) above which the flux-tubes expand radially outwards (and not super-radially). The total expansion ratio is equal to 1 for a flux-tube expanding radially, while a very strongly expanding flux-tube has f >> 1. The continued exploitation of potential field extrapolations with source-surface (PFSS) of magnetogram data lead to associating r 1 with the radius of the source-surface, commonly placed at a fixed height of r SS = 2.5 R , and to evaluation the expansion factor at this height: f ≡ f SS (Wang & Sheeley 1990). This source-surface radius is the one that produces the best matches between the geometry of the extrapolated magnetic fields and the shapes of the coronal structures observed in white-light during solar eclipses, especially the size of the streamers and coronal hole boundaries (Wang et al. 2010;). However, matching quantities such as the open magnetic flux seems to require defining r SS as a function of the solar activity (Lee et al. 2011;Arden et al. 2014), or more generally as a function of the properties of the global coronal magnetic field (Réville et al. 2015). Suzuki (2006) suggested that the terminal wind speed would be better predicted by a combination of the expansion factor and the magnetic field amplitude at the foot-point of any given flux-tube or, equivalently, to the open magnetic flux (see also Fujiki et al. 2015). Other authors also invoke empirically derived parameters such as the angular distance from the foot-point of a given magnetic flux-tube to the nearest streamer / coronal-hole (S/CH) boundary (parameter θ b ; Arge et al. 2003Arge et al. , 2004McGregor et al. 2011b). Recent studies by Li et al. (2011);Lionello et al. (2014a);Pinto et al. (2016); Peleikis et al. (2016) indicate that fieldline curvature and inclination also have an impact on the wind speed. These results altogether motivate the formulation and adoption of semi-empirical strategies for predicting the state of the solar wind. The so-called Wang-Sheeley-Arge (WSA) model is the most success- Figure 2. Maps of the magnetic field at the solar surface, close to the source-surface and at 21.5 R obtained after extrapolating the WSO magnetograms via PFSS and applying the corrections discussed in the text. The orange dots on the first column show the foot-points of the sample of magnetic field-flux tubes used to compute the solar wind solutions. These flux-tubes cover the whole latitude -longitude map uniformly beyond the source-surface, with and angular resolution of 5 • . The amplitude of the open magnetic field is very variable at the source-surface, but becomes nearly uniform faraway from the Sun. Radial profiles of the magnetic field amplitude for a small subset of flux-tubes, at two different Carrington rotations (2055 and 2070). The black continuous lines represent the amplitudes after the correction to the high coronal field was applied to the background PFSS field (which makes the open magnetic field amplitude uniform above 12 R ). For comparison, the blue dotted lines on the second panel represent the same flux-tubes before the correction was performed.
ful and the most widely used relation of this kind. It provides predictions of the solar wind speed close to the Earth's environment at any time, and it does so very fast. But it requires empirical calibration, and cannot be guaranteed to work properly outside its narrow range of validity (typically restricted to the ecliptic plane near 1 AU).
An alternative approach consists of modelling the solar wind acceleration and propagation based solely on physical principles. Comprehensive modelling of the solar wind at a global scale is, however, a daunting challenge. Fully three-dimensional MHD wind models take into account the whole magnetic topology of the corona at the expense of important simplifications to the solar wind thermodynamics and of a large CPU time consumption (Gressl et al. 2013;Lee et al. 2009;Yang et al. 2012, and many others). Some uni-dimensional solar wind models consider much more sophisticated heat transport models, but over-simplify the geometry of the magnetised corona (Pinto et al. 2009; Woolsey & Cran- Radial profiles of the expansion factor for the same subset of flux-tubes and for the same Carrington rotations as in Fig. 3. Continuous black lines and dotted blue lines have the same meaning. mer 2014).
We propose a new solar wind model which assumes a new approach, in between those of the traditional MHD global-scale models and of the more specialized uni-dimensional models. The model consists of computing many 1D wind solutions which sample the whole solar atmosphere (or any smaller sub-domain of interest). The individual 1D solutions are based on a previous uni-dimensional wind model (Pinto et al. 2009;Grappin et al. 2010), modified in order to take the full magnetic flux-tube geometry (expansion, inclination and amplitude of the field) and different heating functions (see Sect. 2). The background magnetic field geometry is currently obtained via potential field source-surface extrapolation (PFSS) from publicly available magnetogram data, but the model is ready to use any other data source (real data, coronal reconstructions or modeled data).
We expect this approach to bring many benefits over the current solar wind models. Our method lets us map the wind speed, density, temperature and magnetic field amplitude across the whole atmosphere, as well as derived quantities such as the characteristic phase Wind speed, plasma temperature and density as a function of height (r − R ) across the whole domain for a random subset of flux-tubes with asymptotic speeds falling close to 270, 360, 500 and 650 km/s. The lines are coloured according to the asymptotic wind speed values using the colour-scheme used in Figs. 1 (second and third columns) and 6 (first column). speeds, wind ram pressure and mass fluxes. The model is computationally light (much lighter than full-fledged 3D MHD models) and has a perfect parallel scaling. Another advantage of our method is that, as there is no cross-field diffusion in the model, the heliospherical current sheets remain thin and do not experience the enhanced (and spurious) resistive broadening typical of global MHD simulations (at the scale of at least a few grid cells). The main limitations of our approach are that the realism of the background magnetic field depends on the reconstruction method used, cross stream interactions are neglected, and so are small-scale crossfield diffusive processes (the effects and strategies for addressing these caveats are discussed in Sects. 2 and 3).

METHODS
The MULTI-VP model consists of a large ensemble of contiguous uni-dimensional wind solutions used to derive the three-dimensional structure of the solar wind. The numerical model is therefore a multiplexed version of a mature one-dimension numerical model representing the heating and acceleration of the solar wind (called, hereafter, the baseline model). MULTI-VP operates, typically, under the following work flow (see Fig. 1): 1. Choice of a magnetogram data source. The source magnetograms can be full carrington maps (e.g from Wilcox Solar Observatory -WSO) or adaptative/forecast magnetograms (for example, from the Helioseismic and Magnetic Imager -HMI -, or from the Air Force Data Assimilative Photospheric Flux Transport model -ADAPT -, at much higher temporal cadences).
2. Choice of a coronal field reconstruction method. We reconstruct the coronal magnetic field and sample out an ensemble of open magnetic fluxtubes extending from the surface of the Sun up to about 30 R and covering all latitudes and longitudes of interest.
3. Computation of field-aligned wind profiles for each one of the sampled flux tubes. The wind model takes full account of the magnetic field amplitude, areal expansion and inclination in respect to the vertical direction along the flux tubes. The model includes a simplified chromosphere, the transition region (TR) and the corona (more details are given below).
4. Assemble the flow profiles and the wind speed, temperature, density and magnetic field amplitude at all the positions desired. We routinely produce maps at 21.5 R which can be used to initiate heliospheric propagation models (such as ENLIL and EUHFORIA).
This framework was designed to be fully modular, such that any of the points above can easily be replaced by other data sources and models depending on the scientific application, and as newer methods become available. In this manuscript, we will use Potential Field Source-Surface extrapolations from WSO synoptic maps covering several Carrington rotations both at solar minimum and at solar maximum (CR 2055 -2079 and CR 2130 -2149). We set a constant source-surface radius R SS = 2.5 R for the PFSS extrapolations and follow the general method and polar-field correction of Wang & Sheeley (1992). We trace an ensemble of open magnetic field lines starting from the source-surface down to the solar surface with a standard angular resolution of 5 • . Each field-line is, at first, assigned a purely radial expansion above the source-surface. This leads to an interplanetary magnetic field amplitude which is very variable within each magnetic sector. The Ulysses mission showed that the radial field component is, however, uniformly distributed with latitude much unlike the amplitudes measured by space probes in the interplanetary field (Balogh et al. 1995). We correct for this by adding an additional flux-tube expansion profile which smoothly (and asymptotically) transforms the very nonuniform source-surface field at r = R SS into a uniform field at about r ≥ 12 R . The correction conserves the total open magnetic flux. Its effects on the properties of the wind flow are thoroughly discussed in Sect. 3. The baseline model used to compute the solar wind profiles follows closely the strategy described in Pinto et al. (2009) and Grappin et al. (2010), albeit with a number of modifications. The numerical code solves the system of equations describing the heating and acceleration of a wind stream along a given magnetic flux-tube where ρ is the mass density, u is the wind speed, and T is the plasma temperature. The wind profiles are computed on a grid of points aligned with the magnetic field (with curvilinear coordinate s), α is the angle between the magnetic field and the vertical direction (cf. Li et al. 2011;Lionello et al. 2014b), and r represents the radial coordinate (distance to the center of the Sun). The divergence operator is defined as where B (s) ∝ 1/A (s), A (s) being the fluxtube's cross sectional area and B (s) the magnetic field amplitude. The ratio of specific heats is γ = 5/3. The terms F h , F c denote the mechanical heating flux and the Spitzer-Härm conductive heat flux, which are both field-aligned. The radiative loss rate is Λ (T ). For simplicity, we define here the mechanical heating flux F h as a function which parametrizes the effects of the coronal heating processes (rather than being the result of small-scale turbulent dissipation; cf. e.g. Sokolov et al. 2013). We assigned it a phenomenological form inspired by those discussed by Withbroe (1988), McKenzie et al. (1995) and Habbal et al. (1995), but depending here on the basal magnetic field amplitude |B 0 |, on the flux-tube expansion ratio f , and on the curvilinear coordinate s The heating coefficient F B0 is proportional |B 0 |, and H f represents an arbitrary damping length which is anticorrelated with the expansion ratio f SS in the low corona (which is expected to lead to higher temperature peaks occurring at lower coronal altitudes and higher mass fluxes in the slower wind flows; see discussions by Cranmer et al. 2007;Wang et al. 2009). We ran a series of trial runs to calibrate the free parameters in Eq. (6), and retained the forms The actual heating processes are still under debate, with physical mechanisms driven or triggered by MHD waves that propagate along open flux-tubes in the corona being favoured by the community. Some of these invoke turbulent dissipation, resulting from the interaction of the upward propagating wave fronts with their reflected counterparts that feeds a turbulent cascade channeling energy from the larger injection scales to the smaller dissipative scales (cf. Tu & Marsch 1995;Matthaeus et al. 1999;Verdini & Velli 2007;Verdini et al. 2012;Cranmer et al. 2007). At the small-scale end of the turbulence spectrum, kinetic effects are likely to be responsible for the actual plasma heating, and specially for the selective heating of different species (Axford & McKenzie 1992;Kohl et al. 1998;Maneva et al. 2015). Another proposed scenario is that low frequency Alfvén waves, which are non-compressible modes, non-linearly convert energy to compressible modes, which dissipate much faster in the corona than their mother waves (Suzuki & Inutsuka 2005;Matsumoto & Suzuki 2014). In all the cases, the injected energy flux density and heat dissipation height are related to the geometry of each flux-tube (e.g steepness of the stratification, gradients of wave phase speeds) as well as to the amplitude of the magnetic field (e.g Poynting flux associated with transverse waves). The formulation for F h in Eq. (6) provides a simple description of these phenomena that we can apply generally (for all flux tubes simulated), but that can be easily substituted in future work by a more precise phenomenology.
The radiative loss function is given by where χ(T ) equals unity for T > 0.02 MK and varies linearly from 0 to 1 for 0.01 < T < 0.02 MK. T M equals 0.2 MK. We employ a nonuniform grid of 640 points between the solar surface (where ∆r = 10 −4 R ) and 31.5 R (where ∆r = 0.3 R ). Time integration is done with a Runge-Kutta scheme of order 3, while an implicit finite-difference scheme of order 6 is used for the spatial dimension, except when computing temperature gradients in the conductive term, for which an explicit scheme of order 2 is applied. Numerical filtering is employed to increase the stability of the schemes (Lele 1992). The top and bottom boundary conditions are transparent to perturbations propagating away from the domain, and are set in terms of the characteristic form of the hyperbolic part of equations (2)-(4), as in Grappin et al. (1997). The non-hyperbolic terms in eqs. (3) and (4) are negligible at both numerical boundaries, except for F h , which is prescribed a priori. Placing the lower boundary at r ≈ 1 R lets us avoid imposing strong constraints on the mass and heat fluxes at the source region of the solar wind. The plasma temperature is practically uniform and constant in the lowest layers of our domain, remaining close to 6000 K and making both conduction and optically-thin radiative cooling vanish there. Density variations are also negligible there, unlike at the top of chromosphere (see Fig. 5  and Grappin et al. (2010) (cf. also background wind solution of Verdini et al. 2012). To speedup the relaxation process, we start off from typical slow and fast wind solutions calculated in advance rather than performing full ab initio computations (i.e, the initial state already has a stable chromosphere, transition region and transsonic wind flow). As the core computing task consists of an ensemble of independent calculations, the parallel speedup scales linearly with the number of processes. We currently compute full synoptic maps (at a 5 • angular resolution) in about 6 hr using 320 parallel computing cores. In other words, the model can run almost in real-time for source magnetic field maps with a cadence larger or equal to 6 hr, or of four maps per day (which is the cadence of the standard PFSS maps currently available via SolarSoft). Other geometrical setups are possible, using for example different angular resolutions or covering smaller fractions of the synoptic map. The total computing load will, in all cases, depend mostly on the number of uni-dimensional samples used and on the typical wind speeds on the regions of the solar atmosphere considered (due to the CFL condition on the integration time-step; wind solutions on regions dominated by slow wind converge faster than those on regions dominated by fast wind).

Magnetic field
The three-dimensional geometry of the coronal magnetic field is given directly by PFSS extrapolations and is represented in Fig. 1 for Carrington rotations 2056 (at the end of solar cycle 23, during the solar minimum of 2008) and 2121 (∼ 1.5 years before the peak of cycle 24). The former is characterized by equatorial streamers covering all longitudes, large coronal holes rooted at the poles, and a well-defined heliospheric current sheet (HCS) which remains close to the ecliptic plane (albeit with a noticeable warp). The latter shows, in contrast, the presence of high-latitude streamers and of large coronal holes rooted close to the equator and a double-lobed HCS. The figure also shows that the resulting slow and fast wind distribution relate to the large-scale geometry of the magnetic field, with slow wind confined to low latitudes during solar minimum and with slow and fast wind streams mixed up in latitude during solar maximum. Figure 2 shows full Carrington maps of the magnetic field obtained after extrapolating the WSO magnetograms via PFSS and applying the corrections described in Sect. 2 at different heights. Each row represents a different Carrington rotation (repectively CRs 2055, 2068, 2079 and 2136). The first three cover a period of time close to solar minimum (years 2007 and 2008). The fourth column corresponds to the time of reversal of the global magnetic polarity at about the sunspot number maximum. In this particular example, the HCS splits up into two separate structures with an almost cylindrical shape (Wang et al. 2014, as discussed by). The last row corresponds roughly to the beginning of the decaying phase of cycle 24. The first column shows the magnetic field at the surface of the Sun (WSO data), the second column shows the magnetic field amplitude just below the source surface, and the last column at about 21.5 R . The orange dots show the positions of the foot-points of all the open flux-tubes considered and mark the coronal holes (the orange dots are not plotted on the 2nd and 3rd rows as they would sample the whole plane uniformly there). As expected, the mix of polarities at the surface arrange into two (or into a few) magnetic sectors higher up in the corona, with one well defined polarity inversion line (PIL) marking the position of the heliospheric current sheet (HCS). The open magnetic field amplitude becomes progressively uniform within each sector with height, but only in the upper part of the corona, well above the source-surface (cf. second and third columns of Fig. 2). This is in good agreement with the measures of the radial magnetic fields in the interplanetary medium by space probes such as ULYSSES. However, the PFSS coronal field reconstruction does not on its own generate such uniform open magnetic fields, and a straight-forward radial extension of the field would propagate the large scale non-uniformities in the magnetic flux up to all heliospheric heights. This justifies the correction applied to the magnetic field geometry discussed in the beginning of Sect. 2, which consists of smooth adjustments to the flux-tube expansion ratios in the higher corona (between 2.5 and 12 R ) such that the open magnetic flux becomes uniform (while keeping the total open flux invariant). The HCS remains everywhere thin, in opposition to global MHD models which require and enhancement of cross-field diffusion on scales of the order of the transverse grid size, which are well above those characteristic of cross-field turbulence and reconnection processes (see e.g Lazarian & Vishniac 1999). Figures 3  and 4 show the radial profiles of the unsigned magnetic field amplitude and of the expansion ratio for a small sample of open flux-tubes at two different Carrington rotations. The blue dotted-lines show the same profiles without the correction to the magnetic flux-tubes in the high corona. These adjustments to the flux-tube expansion ratios in the high corona have a moderate effect on the terminal wind speed distribution, a negligible effect on the wind temperature, but a strong effect on the wind density, and in particular on the correlation between density and terminal wind speed (see Sects. 3.2 and 3.3).

Solar wind properties
Multiple density, temperature and wind speed profiles we obtained are shown in Fig. 5 as a function of atmospheric height (h = r − R ) across the whole computational domain. The log-log scale lets us visualize the full atmospheric stratification, especially in its lower layers. The curves represented in the plot correspond to a sample of about 300 individual wind streams selected randomly from all the CR's in Figs. 2 and 6 that have asymptotic speeds falling close to 270, 360, 500 and 650 km/s. The curves are coloured according to the values of their asymptotic speeds using the same colourscale as the two last columns in Fig. 1 and the first column in Fig. 6 (that is, from dark blue at 250 km/s to dark red at 650 km/s). The slower wind flows converge more slowly to their asymptotic state and exhibit a higher degree of variability than their faster counterparts. The slowest of these streams (dark blue lines) can go through one or more zones of deceleration, unlike the fastest wind flows (dark red lines). These can occur either below the source-surface on flux-tubes with sharp variations of their cross-sections (over-expansion, reconvergence) and inclination (cf. Pinto et al. 2016), and also above due to the smoother variations of their expansion factors in the high corona. Slow wind flows are usually denser than fast flows across all coronal heights, but a few of them can be under-dense between h = 0.01 and ≈ 2 R . Plasma temperature remains invariable in the bottom half of the chromosphere (to machine precision), and density shows only very small variations between all the computed solutions. However, variations in basal magnetic field amplitude and flux-tube expansion alter the heating amplitudes and scale-heights (see Eq. 6), and lead to significant differences in the density of the various wind streams in the higher part of the chromosphere and in the low corona. Typically, slow wind flows are associated with stronger over-expansions in the low corona that lead to reduced H f , to shallower than average density fall-offs with height across the transition region ("evaporation"), and to higher mass fluxes. This effect can be offset or reinforced by variations of |B 0 |, which span a higher range for slow than for fast wind flows. The balance between these variations in heating rate, density stratification, thermal conduction and radiative cooling determines the position of the TR and the temperature profile across the corona. Figure 6 shows a sequence of maps of the wind speed, plasma temperature and density at 21.5 R for several Carrington rotations (at the same as in Fig. 2). The vast majority of the wind streams are close to their assymptotic state, the exceptions being those on the lower tail of the wind speed distribution, such as the 200 − 250 km/s flows appearing on the 4th row of the figure (cf. Sanchez-Diaz et al. 2016). Slow wind flows tend to concentrate on the close vicinity of the HCS, but also occur in the regions surrounding pseudo-streamers. One such example is the feature on the western part of the maps on the second and third lines of the figure. The fastest wind flows always occur within large coronal hole, which are rooted at the polar region during minimum, and sometimes at equatorial latitudes during maxima. The last two rows of the figure show two examples of this situation at the the peak and decay phase of cycle 24 (cf. Wang et al. 2014). The plasma temperature is generally well correlated with wind speed, while the plasma density anti-correlates (see Sect. 3.3 for a quantitative analysis). The velocity and density maps shown in Figure 6 for solar minimum show a close resemblance to those in (Yang et al. 2012), but with the band of slow wind being much thinner (its thickness is closer to that in the WSA model maps).
The combination of the magnetic field amplitude of the corona with the obtained wind speeds and densities lets us deduce the distribution of the characteristic magneto-hydrodynamical phase speeds. Figure 7 represents the Alfvén Mach number M A , which is the ratio between the wind speed and the Alfvén speed (measured in the direction parallel to the magnetic field, equal to the radial direction above the source-surface), in the meridional plane parallel to the plane-of-sky in the middle of Carrington rotations 2055 and 2133 (or equivalently, the meridional planes which cross longitudes 90 • and 270 • in the Carrington map). The white shades indicate the positions for which M A = 1, delineating the Alfvén surface. During solar minimum (top panel), the Alfvén surface assumes typically a prolate double-lobed ellipsoidal shape, with the Alfvén radius r A reaching values close 10 − 12 R over the poles, 5 − 6 R at mid-latitudes, and showing strong incursions at low latitudes. Thin outward "spikes" reaching heights above 15 R can sometimes appear in the close vicinity of the HCS, where the magnetic field amplitude and the local Alfvén speed quickly approach zero. The Alfvén surface becomes closer to spherical during solar maxima (bottom panel of the figure), albeit with some irregularities, with an average Alfvén radius in the range 5 − 7 R . These variations are comparable to those found by Pinto et al. (2011) using idealized solar dynamo and wind models, but with a smaller contrast between the values of r A at minimum and maximum. Overall, the values of r A we obtain are overall lower than those reported by DeForest et al. (2014) from the analysis of inbound propagation of density perturbations in coronograph data (see also Tenerani et al. 2016;Sanchez-Diaz et al. 2017). The effects of setting a fixed source-surface radius R SS on the open magnetic flux amplitude discussed by Réville et al. (2015) probably contribute to these disparities (see also Lee et al. 2011;Arden et al. 2014). Figure 8 compares the plasma density and wind speed along each individual flux-tube at different heights with their asymptotic values for two different Carrington rotations, in order to show how these quantities evolve during the propagation (see McGregor et al. (2011a) for similar scatter plots at higher heliospheric heights). The dashed lines mark the positions for which no net evolution occurs during the propagation between the two heights considered. Points placed above this line indicate that the plotted quantity has decreased during the trajectory, while points below the line indicate an increase. The indices ss, 5, 15 and ip in the axis labels correspond respectively to the height of the source-surface (2.5 R ), 5, 15 and 30 R . The figure shows that all of the fast wind streams are already close to their asymptotic state at 15 R (to within 8% of the asymptotic speed). The slow streams are accelerated more progressively, with a fraction of them falling within an envelope 40% below the value of their asymptotic speed at the same height. Note that the plots show the absolute difference to the asymptotic speeds, which can be of the same order for both fast and slow wind flows, unlike the relative differences (see also Fig. 5). The general trend in our simulations is that the slow wind is accelerated more progressively than the fast wind, in accordance to, e.g., Habbal et al. (1995) and Sheeley et al. (1997). Interestingly, a non-negligible fraction of the wind streams are decelerated in the intervals 2.5 − 21.5 R and 5 − 21.5 R (but not above 15 R ). The decelerations are more frequent and more significant in the slow wind regime, for which the most extreme cases drop from about 400 to 200 km/s between 2.5 and 15 R . The wind speed diagrams show a more complex structure during phases of high solar activity, sometimes indicating the presence of two populations of wind flows in the low corona which converge higher up into a smoother distribution (cf. CR 2138 in the figure). The evolution of the plasma density (normalized for the radial cross-field expansion) shows a rather regular trend, with a progressive reversal of the slope of the diagram between r = R SS and 15 R . Both the initially overdense and underdense streams converge progressively and monotonically to their asymptotic values below 15 R .

Correlations between speed, expansion, density and temperature
The terminal wind speed we obtained correlate with the flux-tube expansion factors, but the degree and sign of the correlation depends on the reference height at which the latter are evaluated. In particular, the wellknow speed-expansion factor anti-correlation is only clearly observed in our solutions if the expansion factor is evaluated at (or close to) the source-surface. Figure 9 shows the dependence of the terminal wind speed (V wind ) on the flux-tube expansion factor measured at two different heights (f SS at r = 2.5 R represented as black diamonds, and f tot at r = 31.5 R represented as red crosses) for CR 2070, which illustrates this disparity particularly well (in general, f tot is simply not well correlated with V wind , with the corresponding scatterplot showing slopes which vary but remain close to flat on average). The expansion that the flux-tubes suffer in the high corona contributes directly to this disparity (f ss would be equal to f tot otherwise), even though it remains smaller than a factor 2 for the vast majority of the flux-tubes considered. The effect on the high-coronal expansion on the actual wind speeds is small. This result is in agreement with the ideas of Wang & Sheeley (1990), who suggested that f SS (expansion ratio measured at the height of the source-surface rather than the total expansion factor) is a good predictor for the terminal wind speed, in spite of the significant spread. However, this particular height does not have a well-defined phys- ical meaning, other than the fact that we are basing our computations on PFSS coronal field extrapolations. It cannot be guaranteed that the anti-correlation between V wind and f SS remains valid for more general types of coronal field geometry. Figure 10 shows scatter-plots of the plasma density and of the wind speed at 3 and at 31 R for CR 2070. The black dots correspond to the final wind solutions we obtained, and the blue dots to the solutions un-corrected for the expansion in the high corona. In the final solutions, the n-V relation shows a high scatter at 3 R and evolves into a well-defined anti-correlation, which is in agreement to the in-situ measurements of HELIOS and ULYSSES. The un-corrected solutions (which have a coronal magnetic field strictly aligned with the vertical direction above the source-surface) show a very different behaviour during the propagation of the wind flow into the high corona. The n-V distributions start off with similar configuration at low coronal heights, but then evolve towards a state which can reach a flat or even positive slope. In this incorrect scenario, the slowest wind streams account for the largest disparities, remaining under-dense all the way through to the upper limit of the domain. The correction applied to the magnetic field turns the non-uniform magnetic field amplitude at r ss into uniform in the heliosphere and, while doing that, emulates to a certain extent the effects of cross-stream interactions on the wind density (see discussion on Sect. 4.2). Figure 11 show scatter-plots of the plasma density (top panel) and of the temperature (bottom panel) against the wind speed at 31 R for the whole set of flux-tubes considered in the present study. The terminal wind speed V wind is very clearly anti-correlated with density for all the data-set analyzed. The plasma temperature is positively correlated with the terminal wind speed, in agreement with available in-situ data, although with a higher scatter (see Elliott et al. 2012). We used the FORWARD tool-set (Gibson et al. 2016) to deduce the white-light emission from our wind model and to build synthetic images of the corona. Figure 12 shows synthetic white-light polarised brightness (WLpB) images of the corona obtained from our simulations for Carrington rotations 2079 and 2136 sided by SoHO/LASCO-C2 images at the corresponding dates. The synthetic images we filtered with a Normalizing Radial Graded Filter (NRGF; Morgan et al. 2006) to enhance the contrast of the coronal features and ease the qualitative comparison. We found that the positions and widths of the main features are very well matched by our simulations for configurations typical both of solar minimum and solar maximum. The main differences between the synthetic and the real coronagraph images relate to the low angular resolution of the magnetograms we have used (5 • × 5 • ), meaning that we cannot capture the finer structure of the streamers and pseudo-streamers, and the absence of transient events (the magnetic and wind models are stationary). The lack of angular resolution is visible on the coronal features near the equator (both east and west) in the first set of images (for CR 2079). The CME visible in the LASCO-C2 image on the bottom row is of course absent in the corresponding synthetic image. We note furthermore that some of the streamers are not strictly aligned with the vertical direction in the C2 images, which probably corresponds to a temporary deflection due to CME activity (see e.g. Rouillard et al. 2012). The corresponding features in the synthetic images are perfectly vertically-aligned. Figures 13 and 14 show Carrington maps of synthetic WLpB built using west limb cuts at two different heights (r = 7 and 13 R ) with the corresponding real maps constructed using STEREO-B/COR2 data made available by the Naval Research Laboratory. Once again we observe that the main features of the WL maps are very well reproduced (positions, slopes and widths) in our simulations, except for the signatures of coronal transients (CME) which appear as vertical traces in the COR2 maps.

Comparison with coronograph observations
4. DISCUSSION 4.1. Strategy, strengths and caveats of the model MULTI-VP adopts a new approach that complements past and present efforts both on modelling the solar wind at global scales using full 3D MHD (Yang et al. 2012(Yang et al. , 2016van der Holst et al. 2014;Lee et al. 2009;Gressl et al. 2013;Oran et al. 2013, among many others) and on modelling the heating and transport processes occurring at smaller scales on the wind flow (e.g Verdini & Velli 2007;Woolsey & Cranmer 2014;Lionello et al. 2014b;Cranmer et al. 2007;Maneva et al. 2015;Pinto et al. 2009;Grappin et al. 2011). MULTI-VP computes detailed solutions of the background solar wind on a arbitrarily large bundle of open flux-tubes extending from the bottom of the chromosphere up to the high corona (typically up to ∼ 30 R ). The model is able to sample large regions of the solar atmosphere (up to a full spherical domain) with more detailed thermodynamics and with significantly smaller computational requirements than the current full MHD global models. MULTI-VP is furthermore unaffected by numercal resistive effects such as the spurious broadening of the HCS. We currently compute the state of the whole corona in about 6 hrs with moderate angular resolution (5 • × 5 • ) and with a moderate number of allocated computing cores. But the total execution can be significantly reduced, as the model is nearly perfectly scalable, and real-time operation can be envisaged. The downsides of the MULTI-VP strategy are that it relies on coronal field reconstruction methods (or any other more or less realistic magnetic field model), it neglects cross-stream effects on the wind, and is only well defined for stationary flows. The underlying numerical model is in fact fully time-dependent, but the setup used is however not well adapted to the study of large perturbations to the background flow, particularly in the direction transverse to the magnetic field. In this manuscript, the magnetic field geometry is obtained externally by means of PFSS extrapolations, which imply the coronal magnetic field to be stationary. We trace out open field-lines, and build an ensemble of magnetic flux-tubes fully taking into account the geometry of each one of them. This translates into prescribing the magnetic field amplitude and polarity, areal expansion rate, field-line inclination and radial height as a function of the distance to the foot-point of each flux-tube. The geometrical description is as general as possible, allowing for any type of field-line bend, kink, over-expansion or re-convergence profile. This is of utmost importance, as both the location and spatial extent of field-line bends and regions with over/underexpansion play a crucial role on the properties of the wind flow which propagate along them (see discussion by Pinto et al. 2016). The model was tested against extreme scenarios which included sequences of exponential expansion and re-convergence, and field-line switchbacks.

Other model parameters
The PFSS extrapolations produce coronal magnetic field of varying complexity up to the height of the sourcesurface (placed at 2.5 R ). If one assumes radial and spherical expansion of the open magnetic field lines from the source surface outwards, the resulting magnetic fields at 30 R will be highly non uniform. This is contrary to in-situ measurements of the interplanetary magnetic field made by the Ulysses mission which showed that the latitudinal distribution of the radial field component is uniform (Balogh et al. 1995). Furthermore, the plasma density distributions we calculated using radially expanding magnetic fields above the source-surface are incorrect. The correlation between wind speed and density is, in particular, very different from the expected one, with a tendency for the simulated solar wind speed to be correlated with density, in contrast with the wellknown anti-correlation measured at 1 AU. After careful analysis, we realized that both issues are related. Flux-tubes with higher than average magnetic field amplitude tend to be over-dense, while flux-tubes with lower than average field strength tend to be under-dense. Also, flux-tubes with higher than average magnetic field amplitude at the source-surface are bound to suffer some additional expansion in the high corona, while flux-tubes with lower than average field strength must undergo the inverse process. It is not possible otherwise to generate an approximately uniform magnetic field (apart from polarity inversion at sector boundaries), as measured consistently in the interplanetary space. We then proposed a correction to the radial variation of the magnetic field which consists of switching from spherical expansion to a smooth flux-tube expansion in the high corona (from a little below the source-surface and up to 12 R ). The individual expansion factors are defined such that the total unsigned open magnetic flux is conserved. In practice, the required additional expansions factors between r = 2.5 R and 12 R are small in respect to the ones for the lower corona (they remain smaller than 2, while the expansion factors between the surface and the source-surface can be as high as several hundreds, see Fig. 3 and 4). This simple correction to the fluxtube profiles in the high corona was enough to produce the correct anti-correlation between the wind speed and density (see Fig. 10).
In more physical terms, the uniformisation of the magnetic field with height in the high corona and heliosphere is most likely a consequence of magnetic pressure balance, with neighbouring flux-tube adjusting their crosssections in order to eliminate transverse pressure gradients. As the magnetic pressure still dominates the transverse pressure budget at these heights, eliminating the pressure gradients implies that the magnetic field becomes close to uniform, with the density of the channeled flows accommodating for the subtle variations in cross-section with height. As described in Sect. 3, these variations in the expansion rate of the high coronal part of the flux-tubes produces a significant effect on the density, but a more moderate effect on the terminal wind speeds obtained. The response of the wind flow speed in respect to variations in tube cross-section A (s) is, to first approximation, proportional to its logarithmic gradient ∂ s log A (s) and to the ratio v/ 1 − M 2 , where v is the wind speed and M is the sonic Mach number (Wang 1994;Pinto et al. 2016). Most of the additional expansion occurs well above the sonic point, such that v/ 1 − M 2 1, and the log-gradient of A (s) remains small.
Having fixed the magnetic field geometry, the only free parameters left in the model relate to the coronal heating functions. The heating functions are empirical parametrizations of the actual coronal heating processes, as the small-scale dissipation mechanisms are still under debate and cannot be accounted for self-consistently in the model. We chose to use a simple formulation based on commonly used phenomenological forms for the heating function with a few extra parameters. More specifically, we used a heating formulation similar to that of Withbroe (1988) with two main modifications (see eq. 6). We made the coefficient proportional to the amplitude of the basal magnetic field B 0 , to account for the energy input to the solar wind resulting from horizontal surface motions at the surface of the Sun. The corresponding energy flux density (the Poynting flux) is proportional to B 0 v ⊥ √ ρ at the surface. We also made the dissipation scale-height depend on the flux-tube expansion ratio, following ideas from Wang et al. (2009). This varying dissipation scale-height increased the contrast between the slow and fast wind speed, making the solutions match in-situ measurements more closely. We have retained the same final form (see Eq. 6) for all calculations to keep our analysis as general as possible.
We did not make any other adjustement to the heating phenomenology, amplitude or scale-height depending on latitude, moment of the cycle or terminal wind speeds obtained. Discriminating between different heating phenomenologies currently under debate is beyond the scope of this manuscript and will be subject of future work.

SUMMARY AND PERSPECTIVES
We present and discuss the design, implementation and testing of a new solar wind model, called MULTI-VP. The model calculates the dynamical and thermal properties of the solar wind from 1 R up to about 30 R , and can cover the totality or a fraction of a spherical domain representing the three-dimensional open-field corona. The model is initiated using an externally prescribed magnetic field geometry. In the current study, we used magnetograms from the Wilcox Solar Observatory and performed standard PFSS extrapolations to derive the structure of the coronal magnetic field. We defined simple phenomenological forms for the heating flux which result in correct angular distributions and amplitudes of slow and fast wind flows. The model provides estimates of the wind speed, density and tempera-ture (as well as any derived quantity such as MHD phase speeds) at a very moderate computational cost when compared to global 3D MHD simulations, and without resorting to semi-empirical hypothesis. The model was designed to be as flexible and modular as possible, such that other sources of magnetogram or coronal field data, different heating scenarii (theoretical or data driven), and other types of domain can be easily setup.
We produced of a series of Carrington maps of the solar wind speed, density, temperature and magnetic field amplitude, which display a correct distribution of slow to fast wind flows at different moments of the solar cycle. The model solutions were used to derive synthetic white-light images of the corona, which compare very well with LASCO-C2 and STEREO COR2 data. The model reproduces correctly the well-know correlations between wind speed, density and temperature measured in-situ by spacecraft in the interplanetary medium. We found that the inverse correlation between the density and wind speed is established in the high corona and that it is a consequence of the small adjustments that neighboring open flux-tubes undergo in order to maintain pressure balance in the transverse direction. This effect is intrinsically related to the uniformisation of the magnetic flux amplitude which takes place between the corona and the interplanetary medium.
Future work will study the inclusion of alternative magnetogram data sources and coronal field reconstruction methods and of more sophisticated heating scenarios (theoretical or data-driven). We plan on performing detailed comparisons between the results of our model with other well-established methods on the near future. We aim at increasing the integration with other solar and heliospheric data and models, and at progressively approaching the requirements of real-time forecasting.
R. F. P. and A. P. R. acknowledge funding by the FP7 project #606692 (HELCATS). The numerical simulations of this study were performed using HPC resources from CALMIP (Grant 2016-P1504). We are grateful to Y.-M. Wang for enlightening discussions and comments on the manuscript, and to R. Grappin for his commitment to the development of the baseline 1D wind model VP.