In this paper, the approximate solution of the solar corona magnetic field is calculated by using the PFSS model with GONG (NSO/GONG: Data Access) measurements of the photospheric magnetic field as the boundary condition. GONG has six stations around the world, and they are ground-based stations that use the helioseismological principle to study the solar interior and can satisfy the near-continuous observation of solar oscillations. The GONG obtains the magnetic field of the entire photospheric surface based on the full-disk magnetogram of the six stations on the photosphere. The magnetogram files used in our study were retrieved from https://gong2.nso.edu/archive/patch.pl?menutype=zeroPoint.

The PFSS model is usually used to solve for the coronal magnetic field. The PFSS model assumes that there is no current in the corona and only uses the radial component of the magnetic field (Altschuler and Newkirk, 1969). We start from the momentum equation of the ideal MHD theory and assume that the corona is in a quasi-static equilibrium and has low-plasma beta. Low-plasma beta means that the magnetic pressure dominates over plasma pressure. In the corona region, we also assume that the magnetic pressure dominates non-magnetic forces (such as gravity for example). By eliminating now in the momentum equation all the terms that become 0 based on these three assumptions, one ends up with only one term: J ⃗ × B ⃗ = 0 . This implies that the corona is free of Lorentz forces. Now, we can substitute in this equation J ⃗ from Ampere’s Law ( J ⃗ = ∇ ⃗ × B ⃗ ), and we end up with: ∇ ⃗ × B ⃗ × B ⃗ = 0 . This equation has two solutions: the first one is that ∇ ⃗ × B ⃗ = 0 , which is known as the potential field approximation that becomes the PFSS solution. From this solution, we now substitute ∇ ⃗ × B ⃗ = 0 in Ampere’s Law and get J ⃗ = 0 . Assuming a scalar potential function ψ with B ⃗ = − ∇ ψ , then solving for ψ can work out the magnetic field of the corona. With the divergence condition ( ∇ ⃗ ⋅ B ⃗ = 0 ), the Laplace equation of ψ is ∇ 2 ψ = 0 . Assuming that the boundary conditions are ∂ ψ ∂ r | r = 1 R s u n = − B r θ , φ ,   ψ | r = R s s = 0 (1) where B is the magnetic field, J ⃗ is the current density, θ ∈ 0 , π , φ ∈ 0,2 π , and the source surface radius R s s is the outer boundary for the potential model. The source surface is the spherical shell where the closed magnetic lines of force within the source surface exist and open magnetic lines of force outside the source surface. The source surface is defined as the surface where the potential becomes 0, and thereby the (open) magnetic field lines emerge orthogonal from this surface. Then, the function ψ for the coronal magnetic field can be calculated.

The simulations of the PFSS approximation are performed employing the POT3D code, which solves the Laplace equation using the finite-difference numerical scheme. The finite-difference approach can match better with the data resolution than with the harmonic approach and can achieve better high resolution for localization (Caplan et al., 2021; Tóth, et al., 2011). Therefore, we use the finite-difference approach to solve the PFSS model. The reader is referred to Caplan et al. (2021) for a detailed understanding of the mathematical representation of PFSS and the solution in the finite-difference format.

We set the grid resolution of r × θ × φ in POT3D to 120 × 180 × 360 , where θ and φ have the same resolution as the input photospheric magnetic field. The three components of the coronal magnetic field B r , B t , and B p , as well as the corresponding coordinates r , θ , φ , are calculated for CR2069 and CR2217 at the source surface radius R s s = 2 R s , R s s = 2.5 R s , and R s s = 3 R s , respectively. The magnetic field components of B x , B y , and B z can be obtained through the transformation of spherical coordinates to Cartesian coordinates in the heliocentric Carrington coordinate system.

The POT3D code adopts the FORTRAN programming language with MPI parallel programming on high-performance clusters for efficiency (Caplan et al., 2021). We run the code in the Chinese supercomputing Tianjin Tianhe No.1A cluster environment Linux system (Wang and Yuan, 2021) (more details about the computer center are available at https://www.nscc–tj.cn/). The POT3D code was installed and configured following the accompanying README instructions.

The input of this code is an approximate map of the photospheric magnetic field in HDF5 and free parameter variables in the DAT format. The output includes three components of the magnetic field with B r , B t and B p in HDF5, the position information of the grid point, the time required for the operation, and the running log files. Then, we convert the three components of the magnetic field into TXT format for subsequent post-processing work.

By obtaining the magnetic field at the source surface, we can calculate f s , V r , and other parameters with the WSA empirical model (Arge et al., 2003). The DCHB model combines f s in the WS model with θ b (Arge et al., 2003). First, the magnetic lines of force are traced from the surface of the Sun to determine the boundary of the coronal hole. If the magnetic lines eventually return to the surface of the Sun, the area is called the closed area; if the magnetic lines reach the source surface, the area is called the open area. When all grid points are traced, the boundary between the open and closed areas is marked as the boundary of the coronal hole. Then, the magnetic lines are traced downward from a certain altitude to the surface of the Sun to determine the location of the footpoint of the magnetic lines. Based on the location of the footpoint of the magnetic line, the minimum angular distance from the boundary of the coronal hole is calculated as θ b (Yang et al., 2018). Through the following equation, f s = R s 2 R s s 2 B r R s , θ s , φ s B r R s s , θ s s , φ s s (2) substituting B r R s , θ s , φ s and B r R s s , θ s s , φ s s , f s can be obtained. The solar wind velocity can be calculated from V r = V s + V f 1 + f s a 1 a 2 − a 3 ⁡ exp ⁡ ⁡ − θ b a 4 a 5 a 6 (3) where a 1 to a 6 are six free parameters. According to the previous research (Li et al., 2019; Yang et al., 2018), these parameters can be set as a 1 = 2.0 9.0 , a 2 = 1.0 , a 3 = 0.8 , a 4 = 0.8 , a 5 = 1.05 , V s = 240.0   k m / s , and V f = 675.0   k m / s for CR2069. Moreover, we wrote our own code to trace the field lines.

With the PFSS model and the WSA empirical model, we can calculate the distribution of the solar wind velocity on the source surface of the corona. Before we compare the simulated V r with the measured velocity at 1AU, first, the measured V r at the first Lagrange point (L1) has to be extrapolated back to the corona. The heliocentric Carrington coordinate system (HECAR) used in the PFSS model has the origin at the center of the Sun, the Z-axis perpendicular to the solar equatorial plane, the X-axis pointing to the direction of Carrington longitude 0° on the solar equatorial plane, and the Y-axis defined by the right-hand rule (Thompson, 2006). The solar wind propagates radially into the heliosphere. Due to the rotation of the Sun, the trace consecutive solar wind parcels originating at the same solar source region left in space corresponds to a Parker spiral (Parker, 1958). The radial component ( B r ) of the observed interplanetary magnetic field (IMF) data is extrapolated back to the source surface along the Parker spiral, provided that the solar wind speed is constant and the magnetic flux is conserved. The relationship between the heliocentric longitude ϕ S S of the HECAR coordinate system back-projected to the source surface, and the actual longitude ϕ L 1 at point L1 is ϕ S S = ϕ L 1 + ω R L 1 − R S S v (4) where R L 1 and R S S are the heliocentric distances between point L1 and the source surface, respectively; ω is the angular velocity of solar rotation; and v is the measured solar wind speed at L1.

Statistical parameters such as the mean square error (MSE), root mean square error (RMSE), and correlation coefficient (CC) are selected for quantitative evaluation of the deduced solar wind parameters from the model.

The MSE reveals the error between the simulated data and the observed data (Allen, 1971), and the formula is M S E = 1 N ∑ i = 1 N y i s − y i o 2 (5) where N is the number of simulated values, y i s is the i-th simulated value, and y i o is the i-th observed value.

The RMSE reflects the error between simulated and observed values (Chai and Draxler, 2014). Its formula is expressed as R M S E = 1 N ∑ i = 1 N y i s − y i o 2 (6) where N is the number of simulated values, y i s is the i-th simulated value, and y i o is the i-th observed value.

The CC represents the correlation between observed and simulated data (Asuero et al., 2006). The value of CC ranges from 1 to − 1 . The closer is CC to 0, the weaker the correlation is. The closer CC is to 1, the more positive correlation between observed and simulated data. If CC is closer to − 1 , it indicates that the observation data are negatively correlated with the simulation data. The formula for CC is given as: C C Y s , Y o = C o v Y s , Y o V a r Y s V a r Y o (7) where Y s and Y o are the ensembles of simulated and observed data, respectively; C o v Y s , Y o is the covariance of simulation and observation; and V a r Y s and V a r Y o are the variances of simulated and observed results, respectively.

By using the POT3D code, we can obtain the extrapolated coronal magnetic field components B x , B y , and B z . The calculated B x is shown in Figure 1. Figure 1 exhibits the coronal magnetic field structure of CR2069 and CR2217 for the cases where the source surface was placed at heights R s s = 2 R s , R s s = 2.5 R s , and R s s = 3 R s . The black lines represent the magnetic field line, and the arrow indicates the direction of the magnetic field. This suggests that the magnetic field lines connecting the surface of the Sun to the source surface come mainly from the polar regions. In addition, there are also several magnetic field lines at the solar surface with footpoints located at low latitudes. Through adjusting R s s , the magnetic field extrapolation results reveal inconspicuous variations in coronal streamers, such that the figures are almost the same for each CR. This is, in general, still consistent with a reasonable model result for solar minima (Badman et al., 2020).

There are eight free parameters in the WSA model such as a 1 − a 6 , V f , and V s . These eight parameters will vary with different model initiations, data sources, or the period under study. According to previous studies (Li et al., 2019; Yang et al., 2018), these parameters can be set as a 1 = 2.0 9.0 , a 2 = 1.0 , a 3 = 0.8 , a 4 = 0.8 , a 5 = 1.05 , V s = 240.0   k m / s , and V f = 675.0   k m / s . The magnetic field structure at R s s as well as f s , θ b , and V r distributions are then obtained for CR2069 (Figure 2) and CR2217 (Figure 3).

In the fourth row of Figure 2, we can see that the magnetic field in the South Pole region extends north at longitude 0 ° ∼ 230 ° , and the magnetic field in the North Pole region extends south at longitude 280 ° ∼ 330 ° . This is because in 2008, the extension of polar coronal holes to low latitudes and independent coronal holes of low latitudes were more common (Wang et al., 2009). In the fourth row of Figure 3, the magnetic field in the North Pole region extends southward at longitude 250 ° ∼ 280 ° and in the South Pole region extends northward at longitude 300 ° ∼ 330 ° . This is similar to the phenomenon that coronal holes in the polar region extended to low latitudes in 2008. By comparing the polar magnetic field in Figures 2, 3, it can be seen that the magnetic field is evidently reversed. This is due to the solar magnetic polarity cycle changes, resulting in the reversal. In Figures 2, 3, low-latitude regions are covered by low-speed solar wind, while high-latitude regions are covered by high-speed solar wind. This corresponds to the structure of the magnetic line of force corresponding to Figure 1. In Figure 2M, we can see that a low-speed solar wind structure is formed at the longitude of 240 ° ∼ 300 ° and near the latitude of 20 ° ∼ 30 ° , which may be related to the pseudo-streamer. Although the foot of the pseudo-streamer is connected with coronal holes of the same polarity, its foot is also close to the boundary of the coronal hole, and the value of the magnetic flow tube θ_b is small, so it can form low-speed solar wind.

Figures 2, 3 yield the following results:

(1) The magnetic field at the source surface will decrease with increasing R s s . Through the conservation of magnetic flux (e.g., conservation of the spherical integral of B r r 2 ), increasing R s s will cause the attenuation of magnetic field intensity. Meanwhile, the number of open magnetic field lines that reach the source surface will reduce with increasing R s s . The corresponding decrease in the open flux on the source surface further aggravates the decrease in the magnetic field intensity. This conclusion is consistent with that in several previous research studies (Arden et al., 2014; Asvestari et al., 2019; Asvestari et al., 2020; Badman et al., 2020; Lee et al., 2011; Panasenco et al., 2020).

According to a previous empirical model (Arge et al., 2003), the solar wind velocity V r can be calculated as Eq. 3. The a 1 reflects the effect of f s on V r , and a 2 − a 6 reveal the influence of θ b on V r . In the WSA model, the optimal values of a 1 − a 6 , V f , and V s will change over time (Riley et al., 2015).

When the source surface radius R s s = 2.5 R s for the CR2069, we set a 1 = 2.0 / 9.0 , a 2 = 1.0 , and a 3 = 0.8 . Then, we look for the best values of the parameters through adjusting a 4 , a 5 , V f , and V s . According to previous studies (Li et al., 2019; Yang et al., 2018), we set a 4 = 0.8 , a 5 = 1.05 , V s = 240.0   k m / s , and V f = 675.0   k m / s . Afterward, the effect of each parameter on the simulation results is analyzed by adjusting one variable each time, that is, a 4 ± 0.2 , a 5 ± 0.2 , V f ± 50.0   k m / s , and V s ± 50.0   k m / s . Then, the comparison between the observed and simulated V r at the source surface is shown inFigures 4, 5.

Figures 4, 5 reveal the following results: (1) a decrease in a 4 leads to an overall increase in V r at middle and low latitudes, and vice versa. However, the change in a 4 does not cause significant changes in V r over high latitudes. (2) A decrease in a 5 decreases V s at middle and low latitudes for faster solar wind and increases it for slower solar winds and vice versa for a 5 . However, the change in a 5 leads to an insignificant change in V r at high latitudes. (3) An increase in V s will result in an overall increase in V r , and a decrease in V s will result in the opposite. (4) The increase in V f will result in increment of the high-speed flow velocity, while a decrease in V f does the opposite. However, the variation in V f has little effect on the low-speed flow.

On the other hand, the variations in MSE, RMSE, and CC are monitored through adjusting different parameters (shown in Table 1). The simulation reveals the results as follows: (1). the decrease in a 4 will cause a reduction in MSE and RMSE, and an increase in CC, vice versa. (2). Decreasing a 5 will make MSE and RMSE smaller and CC larger, while increasing a 5 will do the opposite. (3). An increase in V s will reduce MSE and RMSE, but CC does not change with V s . (4). Increasing V f makes MSE and RMSE smaller, but CC does not change. Therefore, for further optimization, we will decrease a 4 and a 5 , and increase V s and V f based on the initial setting of a 4 = 0.8 , a 5 = 1.05 , V s = 240.0   k m / s , and V f = 675.0   k m / s .

In this section, for CR2069 the parameters for the source surface radius R s s = 2.5 R s are tuned and refined in order to find the local optimal solutions. The V s and V f have a significant impact on the results, which is that V s determines the extreme value of low–speed solar wind and V f determines the extreme value of high-speed solar wind. Thus, it can be adjusted through visual comparison (Li et al., 2019). Figure 4A shows that when a 4 = 0.8 , a 5 = 1.05 , V s = 240.0   k m / s , and V f = 675.0   k m / s , there is a significant difference between the simulated low velocity and the observed low velocity. Therefore, the low-speed flow is adjusted visually to V s = 330   k m / s to be consistent with the observation data and then V f = 1000   k m / s to make the high-speed flow match the observation data, as shown in Figure 6.

To simplify the calculation, only parameters a 4 and a 5 will be tuned, with no change after the adjustment of V s = 330   k m / s and V f = 1000   k m / s . The range of a 4 and a 5 is from 0 to 1.3 in the step of 0.05, with a total of 729 combinations. Then, the two-dimensional (2D) distribution of RMSE and CC with different a 4 and a 5 is given in Figure 7, in which we can find a rough parameter space with smaller RMSE and larger CC. In addition, the details of the good results are listed in Table 2. As shown in Figure 7, a 4 = 1.2 and a 5 = 0.2 correspond to RMSE between 80 and 100 with C C = 0.42 . When a 4 = 0.4 and a 5 = 0.6 , it corresponds to between 100 and 120, with C C = 0.38 . For decreasing a 4 or a 5 individually, the RMSE decreases and CC increases. However, when a 4 and a 5 change together, the RMSE and CC respond irregularly. Table 2 presents the best results among the 729 simulations, with smaller RMSE and larger CC. Considering the evaluation results, a 4 = 0.4 and a 5 = 0.1 and a 4 = 0.25 and a 5 = 0.1 are selected for the graph, as shown in Figure 8. The RMSE and CC values are significantly optimized through the tuning. From the comparison of Figures 8, 4A, it can be seen that the variation trend of solar wind speed in Figure 8 better agrees with the observation than that in Figure 4. Before and after tuning, RMSE decreased by approximately 135 and CC increased by approximately 0.09. There are two notable structures of solar wind high-speed flow in Figure 8, which originate from low-latitude high-speed flow near longitudes 140 ° and 240 ° . This is consistent with the research of Li et al. (2019). Our simulation results can reproduce these high-speed flows with lower peak velocities, but the duration of high-speed flows is consistent with the observed data, and the arrival time of high-speed flows is similar. The maximum velocity of the high-speed flow is evidently underestimated. On one hand, this could be because the tuning parameters have a certain limitation; on the other hand, it may be due to the quality of the observed photospheric magnetic field (Li et al., 2019), and the other limitation is the simplicity of the PFSS model itself.

From the previous studies, it is clear that the height of the source surface in the coronal magnetic field model affects the values of f s and θ b (Arden et al., 2014; Hoeksema et al., 1983; Lee et al., 2011; Sun and Hoeksema. (2009)). Therefore, we set the source surface height as R s s = 2 R s and R s s = 3 R s to study the specific effects of source surface height variation on physical parameters (Arden et al., 2014; Kruse et al., 2021). The optimization steps are the same as in the R s s = 2 . 5 R s case, through tuning a 4 , a 5 , V s , and V f . When R s s = 2 R s , first set V s = 330   k m / s to make the low-velocity flow consistent with the observed data, and then, we set V f = 700   k m / s to ensure that the high-speed flow and the observed data match. Based on the fixed V s and V f , the adjustment range of a 4 and a 5 is from 0 to 1.3, with the step length of 0.05 and a total combination of 729 simulations. The 2D distribution plots of the RMSE and CC of a 4 and a 5 (Figures 9A, B) and the distribution table (Table 3) are analyzed. As mentioned in Section 3.4, the decrease in a 4 and a 5 will give better simulation results. Thus, we select a 4 = 0.05 and a 5 = 0.05 to produce Figure 9C.

When R s s = 3 R s , we can make a similar adjustment of V s = 330   k m / s and V f = 1000   k m / s and then obtain the simulated results given by Figure 10. The 2D distributions of RMSE and CC for a 4 and a 5 are shown in Figures 10A, B, respectively, and the details of the distribution are exhibited in Table 4. In order to get smaller MSE and RMSE with larger CC, we apply a 4 = 0.25 and a 5 = 0.15 to obtain Figure 10C.

It can be seen that the surface heights of the three sources all reproduce the two high-speed flows well, but the peak values and the duration of high-speed flows are not the same. When R s s = 2 R s , the peak value of the high-speed flow is smaller relative to R s s = 2 . 5 R s , but the change trend is relatively good. When R s s = 2 . 5 R s and R s s = 3 R s , the peak value of high-speed flow is relatively large, but it is worse than the trend of actual observation data. The reason for underestimating the peak of high-speed flow is the same as in Section 3.4.

Then, the comparison for the MSE, RMSE, and CC at R s s = 2 R s , R s s = 2 . 5 R s , and R s s = 3 R s shows the following results: (1). when R s s = 3 R s , the desirable value of RMSE is 78.034. When R s s = 2.5 R s , the RMSE is relatively small, with a best value of 74.80. When R s s = 2 R s , the value equals 77.90. Overall, the RMSE exhibits an unremarkable difference at different R s s . (2). When R s s = 3 R s , CC equals 0.34. When R s s = 2.5 R s , CC is approximately 0.42. When R s s = 2 R s , CC is approximately 0.55. It is notable that the CC is increasing with decreased R s s . Thus, for the PFSS–WSA model, the reduction in R s s can give better simulation results.

Arge and Pizzo. (2000) used the WSA-2000 model to forecast the results for 3 years before and after 1996. The CC between the predicted and observed values is 0.4, and the average relative error is 15%. Owens et al. (2005) used the observed magnetic map of NWO and the PFSS + SCS model to forecast the results from 1995 to 2002. The results show that the RMSE of different years ranges from 75   k m / s to 115   k m / s . Gressl et al. (2014) simulated the background solar wind in 2007 using three models: magnetohydrodynamic algorithm outside a sphere/magnetohydrodynamic algorithm outside a sphere (MAS/MAS), MAS/ENLIL, and Wang–Sheeley–Arge/ENLIL (WSA/ENLIL), in which MAS/MAS and MAS/ENLIL significantly overestimated the density of low-velocity flow. Yang et al. (2018) used the MHD model to improve the CC of each parameter. Li et al. (2019) used the MHD model to establish an automated method for systematic quantitative evaluation of simulation results. Compared with these models, the model in this paper adjusts very few parameters and can reproduce the structure of high-speed flow and reflect the variation of source surface height. In addition, when R s s = 2 R s , the CC can reach 0.42, and the correlation is improved. We believe that the solar source surface should drop appropriately with respect to R s s = 2 . 5 R s during the low-solar activity phase of the solar cycle 23. This is consistent with the findings of Lee et al. (2011).