Reconstructing coronal hole areas with EUHFORIA and adapted WSA model: optimising the model parameters

The adopted WSA model embedded in EUHFORIA (EUropean Heliospheric FORecasting Information Asset) is compared to EUV observations. According to the standard paradigm coronal holes are sources of open flux thus we use remote sensing EUV observations and \textsc{catch} (Collection of Analysis Tools for Coronal Holes) to extract CH areas and compare them to the open flux areas modelled by EUHFORIA. From the adopted WSA model we employ only the Potential Field Source Surface (PFSS) model for the inner corona and the Schatten Current Sheet (SCS) model for the outer (PFSS+SCS). The height, $R_{\rm ss}$, of the outer boundary of the PFSS, known as the source surface, and the height, $R_{\rm i}$, of the inner boundary of the SCS are important parameters affecting the modelled CH areas. We investigate the impact the two model parameters can have in the modelled results. We vary $R_{\rm ss}$ within the interval [1.4, 3.2]$R_{\rm \odot}$ with a step of 0.1$R_{\rm \odot}$, and $R_{\rm i}$ within the interval [1.3, 2.8]$R_{\rm \odot}$ with the same step, and the condition that $R_{\rm i}<R_{\rm ss}$. This way we have a set of 184 initial parameters to the model and we assess the model results for all these possible height pairs. We conclude that the default heights used so far fail in modelling accurately CH areas and lower heights need to be considered.


Introduction
The ambient solar wind (SW) and the frozen-in open magnetic field sculpt the medium through which coronal mass ejections (CMEs) and solar energetic particles (SEPs) propagate. Consequently, modeling the SW with high accuracy is a significant component toward more reliable space weather forecasts. Fast SW and open magnetic field are primarily associated with coronal holes (CHs; Levine et al., 1977;Schwenn, 2006aSchwenn, , 2006bWang et al., 1996), which in X-ray and extreme-ultraviolet (EUV) imaging observations can be identified as dark structures in the solar corona. Their dark appearance at these wavelengths is a result of having significantly lower density and temperature compared to their surroundings. Although a plethora of both empirical and intricate magnetohydrodynamic (MHD) models exists, due to its simplicity and comparability to computationally expensive MHD ones Riley et al., 2006), the semiempirical Wang-Sheeley-Arge (WSA) model (Arge & Pizzo, 2000) is the model most commonly employed for estimating the SW flow close to the Sun. Its magnetic field model comprises the Potential Field Source Surface (PFSS; Altschuler & Newkirk, 1969;Schatten et al., 1969) and the Schatten Current Sheet (SCS; Schatten, 1971) models to extract the SW speed and the magnetic field in the corona and then by using a 1-D kinematic scheme attempts to predict the SW speed at Earth. Being an empirical model, the WSA consists of a parametrization including several free parameters, which throughout the decades have been fitted and refitted to conform in situ observations of the SW plasma at 1 AU. One such parameter is the height of the outer boundary of the PFSS model, known as the source surface. In the traditional WSA model it serves as the inner boundary to the SCS model (Arge & Pizzo, 2000;Wang & Sheeley, 1990, 1992. The source surface and its distance from the Sun, R ss , divide the corona into two zones; the inner and outer. Field lines piercing through the source surface are considered as open, while those forming closed loops below it are identified as belonging to closed magnetic structures. Open field lines are stretched out and away in the heliosphere by the SW forming the interplanetary magnetic field (IMF) (e.g., Wang & Sheeley, 2003). The PFSS does not account for dynamic effects that are especially important at the top of streamers (e.g., the high closed loops; Arge et al., 2002). Due to the association of open field to CHs, the areas on the source surface that are defined by concentration of open field lines map down to (one or more) CHs. The location/height of R ss determines strongly the total area of open field regions rooted in the low corona. Lowering R ss will result in more field lines being open, and subsequently in larger areas of CHs. On the contrary, increasing it will reduce the number of open field lines and thus shrinking CH areas. This is visualized in Figure 1, where lowering the source surface expands the areas of open flux. Qualitatively, some field lines that were forming closed loops for the case of a source surface at height 3.0 solar radii are crossing through the source surface placed at height 2.5 solar radii and therefore are not considered closed anymore. Decreasing the source surface height even lower to 1.7 solar radii will result in increasing the number of closed loops that are becoming open and, therefore, in further growing of the open flux areas.
Since fast SW and open magnetic field are related to CHs, for the WSA model to successfully reconstruct the in situ measured SW properties and the open magnetic flux, it should model the size and geometry of CHs as accurately as possible. The dependence of the latter on the source surface height makes the choice of R ss a definitive element in the model's success.
Selecting the optimal source surface height has been a long-standing debate. Earlier studies (e.g., Altschuler & Newkirk, 1969) indicated that the source surface is best placed at 2.5R ⊙ (R ⊙ = 1 solar radius). This value has been adopted by many as the canonical height; however, a range of values from 1.6R ⊙ to 3.25R ⊙ have been used (e.g., Arden et al., 2014). Criteria for the choice of these values have been (i) the sector structure of the Heliospheric Current Sheet (HCS; Hoeksema et al., 1982Hoeksema et al., , 1983Hoeksema & Scherrer, 1986), (ii) in situ and remote observations of the HCS by Ulysses (Phillips et al., 1994;Wang & Sheeley, 1995) and during total solar eclipses (Schatten et al., 1969), (iii) the IMF polarity (Hoeksema et al., 1983), (iv) the CH boundaries, (v) the IMF strength (e.g., Lee et al., 2011), and (vi) other coronal signatures that can be sources of open field (Levine, 1982). While plenty of studies were restricted to and/or concluded the use of a singular value for R ss (Altschuler & Newkirk, 1969;Arge & Pizzo, 2000;Hoeksema et al., 1983;Linker et al., 2017;Luhmann et al., 2002Luhmann et al., , 2009Luhmann et al., , 2013McGregor et al., 2008;Riley et al., 2006Riley et al., , 2015Wang, 2009), others have been suggesting a solar cycle-dependent source surface height (e.g., Arden et al., 2014;Lee et al., 2011;Levine, 1982), a concept known as the "breathing" source surface (termed by Arden et al., 2014). In addition to that concept, the idea of a nonspherically shaped source surface has also been proposed (e.g., Levine et al., 1982;Riley et al., 2006;Schulz et al., 1978;Schulz, 1997Schulz, , 2008. Aiming to eliminate a known discontinuity in the WSA model, namely, kinks in the field lines occurring at the transition from the low-coronal PFSS to the high-coronal SCS domain, caused by the PFSS solution forcing the magnetic field to be radial at the source surface, McGregor et al. (2008) repositioned the inner boundary of the SCS model at a radius R i lower than the source surface (R i < R ss ). This concept is illustrated 10.1029/2019JA027173 in Figure 1, where the red dotted curve encircling the Sun represents the new inner boundary of the SCS model and the gray dashed curve is the source surface placed at height R ss . As discussed in McGregor et al. (2008), these field line kinks are produced by a discontinuity in the latitudinal and longitudinal components of the field at the source surface. This implies the presence of surface current flows there. McGregor et al. (2008) introduced the lower SCS model inner boundary to minimize the currents flowing, and they demonstrated that not only they achieved that to some extent but, in addition, these changes improved SW speed predictions at Earth. Consequently, placing the inner boundary of the SCS modeling domain lower than the source surface boundary lead to an overall advancement of the WSA model. In establishing the optimal pair of [R i , R ss ] heights for their new adaptation of the WSA model, McGregor et al. (2008) aimed for conservation of the absolute total magnetic flux between the one from the original WSA model and that from their updated version. The purpose of that was to facilitate the comparison between the two models, original and improved WSA. The pair they deduced in their study is [2.3, 2.6]R ⊙ for R i and R ss , respectively. EUHFORIA (EUropean Heliospheric FORecasting Information Asset) is a novel heliospheric wind and CME evolution model (Pomoell & Poedts, 2018) McGregor et al. (2008). In a comprehensive study of modeling the background SW with EUHFORIA using the adopted WSA model and default heights (Hinterreiter et al., 2019), it became evident that the capability of the model to predict a high-speed stream at Earth is often quite low, with predicted high-speed streams either shifted in time, having lower than expected amplitude in the SW speed profiles or being overall absent. Similar results have been reported in other physics-based reconstructions of the ambient SW (Gressl et al., 2014;Jian et al., 2011;Lee et al., 2009;Owens et al., 2008). The novelty of our study is the assessment of the PFSS+SCS capability to reconstruct CH areas with accuracy, knowing the relation between CHS and high-speed streams.
In the model the source surface and the SCS inner boundary positions are adjustable parameters. Considering the freedom in the choice of the source surface height, already brought forward by other studies discussed above, we focused first in qualitatively assessing the default values of [R i , R ss ] and subsequently in investigating the possibility of other height choices. To do that, we selected a sample of 15 CHs located around the central meridian of the solar disk as viewed from Earth, but at different latitudes (see section 2.2). Using catch (Collection of Analysis Tools for Coronal Holes; Heinemann et al., 2019), we extracted the boundaries of these CHs based on their appearance in EUV filtergrams at 193 Å (section 2.3). We reconstructed their areas with EUHFORIA, by varying the paired values for the heights of the source surface and the inner boundary of the SCS model, in order to find optimal values that produce a better match between modeled and observed CH areas (section 2.4). A quantitative analysis of the results was performed by defining two parameters, the coverage (cov) and the connected pixels area of open flux associated to the CH, of which we consider the percentage that lies within the EUV defined CH boundaries (CA in ). By coverage, in the context of this paper, it is meant the percentage of the CH area observed in EUV that is successfully modeled as open field area. And by connected pixels we consider the open field pixels that lie both inside and outside the EUV-defined CH boundary, and which form a continuous open field region that has the same polarity as the CH in study. Both parameters are analytically described in section 3.
From this study we conclude that the PFSS+SCS part of the adopted WSA model for the default heights of [2.3, 2.6]R ⊙ does not properly reconstruct CH areas. Lowering the heights [R i , R ss ] improved the model results; that is, open field lines are rooted to areas in the lower corona that better represented the CHs observed in EUV images. However, it also led in the presence of open field lines that mapped down to the corona at areas lying outside the CH boundaries defined via observations. This is an undesirable modeling effect, and we attempted to identify optimal height pair candidates by selecting those who lead to a modeling result that has improved CH area reconstruction without expanding the CH area beyond the expected boundaries.

Observational Data
To compute the SW and magnetic field in the corona, the WSA model requires a magnetogram as input. The model performs a global corona simulation and outputs open-closed field regions for the entire solar surface. However, we focus our analysis to the results for the front side of the Sun, from Earth's perspective, where the selected CHs are located. We also do not analyze results toward the limb, with one exception, that of the CH on 3 January 2017, which is linked to a polar CH. Even though a variety of magnetograms are available, we used exclusively magnetograms provided by the Global Oscillation Network Group (GONG) and developed with the Air Force Data Assimilative Photospheric flux Transport (ADAPT) model. The aim was a coherent output and to avoid any effects on the modeled results due to quantitative differences among magnetograms , as a result of being obtained at different observatories, with different instruments, and constructed following various distinct techniques. These are synchronic magnetograms, meaning they are developed such as to resemble with high accuracy the magnetic field of the solar surface on a particular time. A detailed description of the ADAPT model can be found in Arge et al. (2010) and Hickmann et al. (2015), while the GONG-ADAPT database from which magnetograms were acquired is available online (ftp://gong. nso.edu/adapt/maps/gong/).
For extracting the CH boundaries we used EUV observations during 2012-2017 from the Atmospheric Imaging Assembly (AIA; Lemen et al., 2012) on board the Solar Dynamic Observatory (SDO; Pesnell et al., 2012). For this study we utilized the full disk coronal images (1,024 × 1,024) at a wavelength of 193 Å, which are made available by the Joint Science Operation Center (JSOC). Although it is argued that CHs appear differently when observed at different wavelengths, that is, varying boundaries, the 193 Å filter is most commonly used. The response function of this filter is most sensitive in the temperature range of quiet sun and CH regions resulting in a high contrast. Using these data, we extract the CH area and its boundary in a manner described in section 2.3.

Selected CHs
The 15 selected CHs are shown in chronological order in Figure 2. As mentioned earlier, all CHs have longitudinal position of center of mass (CoM) located along or within 10 • of the central meridian of the Sun. The CoM longitude is calculated from the central meridian viewed from Earth in Stonyhurst coordinate system, where the central meridian has longitude equal to 0 • . Although negative values of the CoM longitude are possible, the CH sample we worked with, nondeliberately, includes CHs with positive CoM longitude. The latitudinal position of their CoM varies from low to high latitudes. Limiting the sample to CHs located only along the central meridional zone provides the opportunity to focus our assessment on the possible impact of the CH latitudinal position in modeling their areas. In addition, the sample covers the extended maximum, declining, and early minimum phases of Solar Cycle 24, from 2012 to 2017, enabling the investigation of potential solar activity effects in reconstructing CHs and thus in the choice of the source surface height, an idea that has been suggested before (Arden et al., 2014;Lee et al., 2011).
Out of the 15 CHs, 11 are latitudinally and 3 longitudinally elongated, while one is rather small and bears a more round appearance. Regarding their CoM latitude, 5 and 3 CHs are low latitude ones positioned in the Northern and Southern Hemispheres, respectively, while 5 and 2 CHs are midlatitude ones in Northern and Southern Hemispheres accordingly. Due to known weaknesses on magnetic field and EUV observations in polar regions (i.e., line-of-site effects), no polar CHs were selected for this study; however, the CH on the 3 January 2017 is clearly connected to a polar CH and is an exception to our sample. In this particular case we do consider the polar part that is visible in the EUV image, acknowledging the impact this will have on the result. The CH on 20 December 2016 appears to be associated to a polar one as well. A faint channeling between the CH and the polar one is visible in the EUV images. Nevertheless, we assess it individually, as its boundaries can be defined and are separated from those of the polar CH. The sample also consists of CHs exhibiting different level of patchiness that is defined by the presence within the CH of large areas with a dipole magnetic field configuration (i.e., 3 January 2017) or due to the CH consisting of a cluster of small dark regions (i.e., 27 March 2013).

Extracting CH Boundaries With CATCH
The CH boundaries were extracted using catch (Heinemann et al., 2019), an intensity threshold algorithm that is modulated by the intensity gradient along the CHs boundary. The SSW-IDL-written GUI enables the user to optimize the threshold used for the CH extraction by considering the average intensity profile at the boundary to find the highest gradient. catch operates under the basic assumption that the optimal boundary lies where the intensity gradient is strongest. Additionally, catch is the first algorithm that estimates the uncertainties in the extraction to quantify the boundary. Generally, catch can be used to detect and extract CHs from multiple spacecrafts (SDO, SOHO, and STEREO), but for this study we used only the EUV 193 Å filtergrams from AIA/SDO. For better quantifying the uncertainties in the PFSS+SCS computed open magnetic field when comparing to the coronal remote sensing observations, in addition to the optimal CH boundary (as defined by catch), we also define lower and upper bounds for each CH. The lower and upper bounds present a significantly overestimated and underestimated threshold for the extraction. For further information on the algorithm, we refer to Heinemann et al. (2019). Line-of-sight effects can introduce errors in the boundary extraction of the CH with CATCH; however, these errors are of the order of 2 to 8 arcseconds. EUHFORIA modeling domain is divided in a mesh grid that has resolution of 2 • per pixel. Thus, the errors introduced by CATCH in the boundary location is below the EUHFORIA output resolution.

Reconstructing CH Areas With EUHFORIA
The reconstruction of open flux areas was done using EUHFORIA. For each of the selected CHs we run the coronal model of EUHFORIA for a large sample of [R i , R ss ] heights. To define the height pairs, we started with R i taking the value of 1.3R ⊙ up to 2.8R ⊙ with a step of 0.1R ⊙ , and for each R i we varied R ss from 1.4R ⊙ to 3.2R ⊙ with the same step, and the condition that R i < R ss . ADAPT magnetograms, at date-time same as that of the EUV images analyzed, were used as input to EUHFORIA, which then produced maps of the open and closed magnetic regions for the specified input heights R i and R ss .
To create these open-closed flux maps, the solar surface is divided in a regular mesh consisting of pixels covering 2 • × 2 • per pixel. For each pixel the procedure assigns a field line and traces it from the solar surface upwards toward the source surface. If the tracked field line is curved below the source surface and can be traced back down to the solar surface, it is assigned as a closed one. Per contra, if the traced field line pierces through the source surface, it is accredited as an open one. Examples of open-closed field areas are given in Figure 3, where dark blue and red colors represent areas of open flux, with positive and negative polarity, respectively, while light blue and orange areas are correspondingly closed flux areas. The light green contours are the optimal CH boundaries derived from EUV observations using catch, as described in the previous section. The results presented in this figure are the output of the model run for the default pair of heights, [2.3, 2.6]R ⊙ . To better visualize the modeled CH area and the overplotted EUV-based boundaries, especially in the cases of very small CHs, we only plot in Figure 3 the front side of the Sun (Earth view). In addition, for CHs that do not extend to polar regions, high latitudes of Northern and Southern Hemispheres are excluded from the images.

CH Reconstruction Using the Default Model Boundary Heights
The examples in Figure 3 display maps based on the default pairs of heights [2.3, 2.6]R ⊙ used as default values in EUHFORIA. These values are the ones concluded in McGregor et al. (2008). It is clear from this figure that, for the majority of the CHs, the model runs based on the default heights [2.3, 2.6]R ⊙ fail in reconstructing well the area and geometry as compared to the EUV extraction. Not only are the CHs modeled thinner and overall smaller but also for some cases (i.e., 9 June 2012, 23 May 2013, and 18 January 2014) the CHs appear to be shifted. For the CH of 24 June 2014 the model is unsuccessful in modeling open flux both within the expected boundaries and in the surrounding area. This CH seems to be invisible for the PFSS+SCS models in the adopted WSA model. A CH that appears to be well modeled by EUHFORIA is the one on 3 January 2017, which is a southern polar CH with a large extension that reaches to low Northern Hemisphere latitudes. It is worth pointing out, however, that the CH boundaries extracted with catch are only for the part of the CH that is visible from Earth's line of sight. This explains why also in the EUHFORIA output open/closed field map we do not focus on longitudes beyond ±90 • .
In order to quantify how successful the model is in reconstructing CH areas, we define the coverage parameter, cov, which is given by where N open is the number of open flux pixels contained within the optimal catch boundaries (green contour) and N total is the total number of pixels enclosed by that same boundary (open and closed flux pixels).  Thus, the coverage is the percentage of the expected open field area, enclosed by the green outline in Figure 3, which the model has successfully predicted as open field area. Low coverage indicates not only that the area size of the CH is not correctly modeled but also that the CH could possibly be modeled shifted in space. High coverage suggests that both area size and location are possibly well captured by the model reconstruction of the CH. The coverage for the model results using the default source surface and SCS inner boundary heights, [2.3, 2.6]R ⊙ , is given in the right panel of Figure 4 as green colored circles. The x axis is the nominal number of the CHs when they are accounted in chronological order (Column 1 of table 1). It can be seen that for 13 out of the 15 CHs we studied, the model gives a coverage below 60%, with the majority of CHs ranking below 30% coverage. This is a clear indication that the model under the default setup does not properly model the CH areas.
We also investigated the possibility of comparing the modeled CH areas to the ones defined by the smaller and larger area boundaries obtained from the EUV images using catch. An example of how these boundaries differ from each other is shown in the left panel of Figure 4. In yellow we show the smaller area boundaries and in magenta the larger area boundaries. For most CHs the differences between the areas  This concludes that there is no systematic improvement if smaller area boundaries were to be obtained from the EUV observations when a smaller threshold is considered. In a similar manner, larger area boundaries do not worsen the result for all CHs.

CH Characteristics and Their Effect
The open-closed flux maps given in Figure 3 and the coverage, cov, shown in Figure 4 indicate that the model performs better for some CHs comparing to others when the model runs are setup using the default pair of heights. In order to exclude the possibility of the model showing preference in better modeling CHs bearing particular features, we assess how the coverage relates to CH characteristics. The elongation, latitudinal position of CoM, patchiness, area size, and mean intensity are apparent features of a CH that are of interest to this study.
The elongation of a CH can have an effect to the modeled results due to the way the selected magnetograms are constructed. Even though dynamical processes are applied to them, synchronic ADAPT magnetograms best represent the magnetic field on the Sun along the central meridian as viewed from Earth. This suggests that CHs, which are latitudinally elongated and lie within the central meridional zone, can be possibly better modeled when selecting a magnetogram from the date the CH was located there. Thus, a longitudinally elongated CH might not be as accurately modeled at its full length. So when studying a CH on a particular moment, it can also have an impact to the modeled output based on that perspective.
CH mean intensity and area are parameters extracted using catch, as described in the section 2.3. The mean intensity is automatically computed with catch based on the threshold used for extracting the CH boundaries. It is expected that active regions in the vicinity of a CH will have an impact on its configuration, which is also imprinted on magnetograms, and thus has the potential to affect the modeled CH areas. The sample's average mean intensity is 31.3 DN, and only 6 out of the 15 CHs studied have mean intensity above that value. The size of a CH is expected to affect the width/duration of a high-speed stream as well as its speed and open flux and thus is a parameter we investigate (Hofmeister et al., 2018;Nolte et al., 1976;Rotter et al., 2012;Vršnak et al., 2007). In terms of the area the CH sample consists of both small and large CHs, with the biggest on 3 January 2017 extending from polar to solar equatorial latitudes. The sample mean is 8.7 × 10 10 km 2 , and four CHs have area above this value.
In Figure 5, the left panel presents the coverage, cov, with respect to the CoM latitude. From the level of scatter in the points, it is clear that the latitudinal position of the CH does not affect the performance of the modeled result. The same applies for the other two parameters shown in the middle panel (mean intensity) and the right most panel (area). We know that in estimating the mean intensity of CHs, we have to be aware of multiple factors that influence the uncertainty of the absolute values of the measured intensity. Any instrument (in this case AIA) suffers inevitable depletions (e.g., decreasing sensitivity, aging filters, and others), distortion of optical system, and scattered light to name just a few. But in the context of this study, especially regarding the CH detection, we do not believe that these uncertainties introduce significant errors, as each CH was manually tuned to account for such instrumental and observational changes (see for more details Heinemann et al., 2019). For the case of large CHs we have only one point so the relation between coverage and area remains inconclusive. We also checked whether the visual characteristics of elongation and patchiness of a CH have any effect. No such conclusion could be made; however, the number of CHs is too low for the result to be compelling. We also accounted for the possibility of a coupling of different parameters, for example, area or intensity and CoM latitude, which could have an effect in the modeled results. But even for this scenario no trend is apparent with respect to the coverage. One conclusion that could be made is that all CHs located in the Southern Hemisphere of the Sun have low mean intensity, while those in the northern have high. This can possibly be a solar cycle related effect, since we studied CHs from one SC only (SC24), and indeed, from 2013 to 2015 Southern Hemisphere clearly dominated both in the number and in the area size of sunspots (Li et al., 2019; see also sunspot number information at http://www. sidc.be/silso/ssngraphics). Our sample size is rather small for an unambiguous conclusion to be made, so a more detailed analysis of a larger CH sample is necessary.

Assessing for Systematic Shifting and Solar Cycle Trends
Maps of open-closed flux indicated the possibility of shifting of the modeled CH areas with respect to the expected location. To assess this, we investigated the likelihood of a systematic shifting. All directions (eastward/westward and northward/southward) were investigated. In addition, we considered both a 2 • and a 4 • shifts in each direction. From this analysis no systematic effect could be identified. The role of CH characteristics, such as mean intensity, elongation, size, CoM latitude, and hemispheric position, in shifting effects was evaluated independently. Regardless of whether the sample is assessed as a whole or divided in groups based on the CH characteristics, the conclusion for possible shifting remained negative for all cases.
In earlier studies (Arden et al., 2014;Lee et al., 2011;Levine, 1982) the concept of a solar cycle-dependent source surface height has been suggested, that is, varying within the course of the solar cycle, as well as from one solar cycle to another. We investigated this prospect at first for the whole CH sample and then by separating CHs to groups based on their apparent characteristics mentioned also in the shifting investigation. We found no solar cycle trends when assessing the result for the default heights, but also for our full set of [R i , R ss ] heights discussed below. We note that the number of considered CHs is rather small and, therefore, no conclusive results can be drawn about the systematic shift. In addition, the period over which the sample extends does not cover a complete solar cycle, so to understand solar cycle effects on the modeling output, it is important for such a study to consider a sample that spans over the entire solar cycle.

Finding the Optimum Paired Values for [R i , R ss ]
After running the model for all 184 pairs of [R i , R ss ] heights, it was made evident that lowering the source surface height significantly improved the modeled CHs. One very clear example is that of 24 June 2014, shown in Figure 6. This particular CH was invisible in the open-closed flux maps created with the boundaries placed at the default pair of heights. As can be seen from the right panel in Figure 6, the model nicely maps the CH area when a lower pair of heights is considered. We specifically compare the default heights output to this one because the source surface height of 1.8R ⊙ has been suggested by other studies as a better choice for particular periods and has been used in earlier studies for comparison purposes (e.g., Lee et al., 2011).  the modeled result; that is, the modeled open flux areas grow to cover more of the total expected area outlined by the EUV defined CH boundaries and thus better capture it. This is reflected in the improvement of the coverage, cov, parameter plotted on the y axis of Figure 7. For some cases the coverage even rises from 0% to almost 100%. Overall, a lower source surface height has the same effect; however, for a fixed low value R i a decrease in the source surface height does not necessarily have a high impact. The large CH that consisted of the southern polar one and its long extension to equatorial latitudes was very well modeled by EUHFORIA run with the default values. Lowering the values of the pair of heights below the default ones did result in improvement of the coverage; however, as it was already high the improvement was not significant (within 10%). The CH on 27 March 2013 is a very patchy CH that is not well modeled by EUHFORIA's adaptation of the WSA. The coverage is below 60% regardless of the selected pair of heights considered for running EUHFORIA. Another patchy CH, but significantly less patchy than that on 27 March 2013, is the CH on 20 December 2016. This CH is a bit better modeled when lower values of the pair of heights are considered, but the coverage remained below 70%. Reconstruction of these strongly (27 March 2013) and moderately (20 December 2016) patchy CHs suggests that the model might have difficulty in reconstructing patchy CHs. It will be interesting for this to be investigated further by considering a larger sample of patchy CHs. Also, it will be interesting to assess possible impact in modeling high-speed streams originating from these CHs.
Even though there is a clear indication that a pair of lower heights will result in a significantly improved coverage for all CHs, not a single pair could be specified as the ideal one. In addition, lowering the two heights results in opening more flux to the heliosphere, not only within the CH areas but also outside; namely, open field areas appear to be overestimated. This is illustrated in Figure 6 where the modeled CH area extends beyond the EUV-defined boundaries, green outline, both to the north and to the south. This total open field area, appearing in blue color on the map, we call the connected pixels area. The portion of this connected area that lies outside the EUV-defined CH boundary we call connected pixels outside, while the portion that lies within the CH boundary we call connected pixels inside. This feature is present in all CH cases and rises the question of how low one can place the two heights, [R i , R ss ], in order to best model the CH area, but without overestimating the open flux outside the EUV defined CH boundaries. To answer that question, we consider the connected pixels area. To identify these areas, we apply the region growth method by selecting seed pixels within the CH area that are assigned as open field. We calculated the total number of connected open flux pixels, N total , and divided them to those that lie within the boundaries, N in , and those that lie outside, N out . From these we define the percentage of connected area within the boundary, CA in , which is given by In Figure 8 we plot this parameter against the coverage. We can see that for the majority of the CHs both parameters at first increase but then the system reaches a point where the percentage of connected area that lies within the boundaries starts decreasing. This implies that although open flux pixels grow, the majority lies outside the EUV defined boundaries. If flux was opening only inside the boundaries, then the curves would be approximately linear. This does not happen though and a saturation point is reached from which further decreasing of the two heights opens more flux outside the boundaries than inside. Following the standard paradigm that the primary source of open flux are CH areas, we can conclude that too low heights lead to nonphysical results due to cutting actually closed loops. This is an important finding since it limits how low the two heights can be. From the same figure one can notice that for some CHs there is a double saturation limit.
These saturation limits identified in Figure 8 can be used in order to constrain the possible pairs of boundary heights to those that provide an improved result both in coverage, cov, and CA in . We define two criteria with the first being that the coverage parameter needs to be above 60% and the second is that CA in should lay no less than 5% from the saturation limit. Although, for most CHs the coverage is above 70%, for four of them (i.e., 27 March 2013, 12 December 2013, 20 December 2016, and 28 December 2016, the maximum coverage is well below that value, which is the reason for this 60% limit. Especially for the CH on 27 March 2013 the coverage parameter is below 55%; thus, a special limit of minimum 50% coverage is applied. For 8 out of the 15 CHs the value of CA in is below 50%. To this, we include the CH on 9 June 2012 for which the high values of CA in are only for when the coverage is below 10%, which indicates that the number of open flux pixels is extremely low. A CA in < 50% is the result of having more open flux lying outside the CH area defined by the EUV boundaries for all the model runs performed. For some of the CH, a secondary saturation limit can be identified less than 10% from the first one, which, however, results in better coverage without affecting greatly CA in . For those cases we do consider the second saturation limit for the CA in in investigating the optimal [R i , R ss ] pair. After analyzing the results for each CH, we collected the pairs of heights that justify the criteria posed. The table is presented in Figure 9. The color map indicate the number of CHs for which the pair fitted the criteria. Although all the height pairs fulfill the criteria for at least one CH, some appear as good candidates more frequently. As can be seen in the

Conclusion and Discussion
As already discussed in other works (e.g., Arden et al., 2014;Linker et al., 2017;Wallace et al.,  In this work we presented a collective study assessing the capability of the PFSS+SCS part of the adopted WSA model by EUHFORIA forecasting tool in modeling CH areas. Although values between 1.6 to 3.25 solar radii are considered allowed for the source surface height, a commonly used value is 2.5R ⊙ (e.g., Arge & Pizzo, 2000;Riley et al., 2015;Wallace et al., 2019). On the other hand, McGregor et al. (2008) concluded that for the modified WSA model that they presented, and which is also adopted in EUHFORIA, the optimal heights of the source surface and the SCS model inner boundary are [2.3, 2.6]R ⊙ , which we considered as the default pair of heights. Other studies suggested an even lower or higher source surface height (Lee et al., 2011;Phillips et al., 1994) or even a solar cycle varying source surface (Arden et al., 2014;Levine, 1982). We considered 184 pairs of heights aiming to investigate whether the modeled result can be improved. We selected 15 CHs having CoM latitude within the central meridional zone of the Sun, and a variety of morphology defined by their elongation, area size, intensity, and level of compactness. Using EUHFORIA, we reconstructed the CH area based on model runs for the 184 pairs. The modeled areas were compared to those derived with catch using EUV imaging observations.  Levine (1982), and Lee et al. (2011) investigated the possibility of a solar cycle varying source surface height and concluded that a lower source surface height during solar minimum periods is required. We investigated this possibility, but we identified no solar cycle trends in the success of the model in reconstructing CH areas; however, capturing the whole cycle could improve the results. In addition, we investigated a possible latitude dependence; however, the latitudinal position of the CH did not have any effect on the modeled results. These are interesting results, which, however, need to be considered with care due to rather small sample size. A more detailed analysis, considering a larger number of CHs, spanning over more than one solar cycle and over all solar cycle phases, and having more variety in their latitudinal positions, can potentially lead to a different outcome.
An interesting remark is that for the CH on 24 July 2014 a channeling to a polar CH was present regardless of the selected pair of heights. This was also appearing for other CHs in the sample when lower heights were selected for the source surface and the inner boundary of the SCS model. The idea of CHs being linked to polar ones via narrow passages in the photosphere has been discussed in, for example, Antiochos et al. (2007Antiochos et al. ( , 2011. It will be interesting to investigate this aspect on a larger sample of CHs. At this point it is worth pointing out that, although comparing model reconstructed open field areas to EUV extracted CH areas is a concept already incorporated in other studies (e.g., Lowder et al., 2014;Linker et al., 2017), the two do not necessarily capture the same phenomena. We assume that a CH, and subsequently an open flux area, can be well identified in EUV images as a dark region. However, the areas surrounding CHs might also contain open field, even though due to strong presence of closed field there they appear bright in EUV images and are thus not recognized as part of the CHs (Linker et al., 2017). As it has been suggested, open flux could also originate from outside CHs (Fisk & Zurbuchen, 2006), such as foot points of active regions, which do not appear dark in emissions (Linker et al., 2017;Luhmann et al., 2008, and references therein). Therefore, due to these obscuration effects, the EUV images provide an open flux area that can be an underestimate of the actual one. Moreover, magnetograms are momentary representations of the magnetic field configuration of the photosphere and thus do not capture the short-scale dynamic processes that can contribute to the total open field areas. This can lead to underestimated modeled open flux areas. The effect of a nontime-evolving input boundary has also been suggested in Riley et al. (2019), where they are comparing MHD and PFSS reconstruction of the structure of the coronal magnetic field using the same input magnetogram. Although Riley et al. (2019) reduced the source surface height to 2.0R ⊙ for the PFSS model, still, the predicted open flux areas were smaller than the MHD prediction, which has been considered in the past as more definitive. The impact of time-dependent phenomena that are obscured in EUV, of the rapid evolution of the photospheric field, of the steady state assumption of the PFSS model and of the synoptic maps considered should not be neglected (Luhmann et al., 2002(Luhmann et al., , 2008(Luhmann et al., , 2013Wallace et al., 2019).
It is necessary to mention that the outcome of this study, and other similar studies before it, is focused on how the two free parameters involved in the PFSS+SCS, and subsequently in the WSA model, can affect the modeled results. Saying that we aim to bring forward that the missing open flux in the model outputs might likewise require fine tuning of the other free parameters involved in the model, an idea also discussed in Wallace et al. (2019). EUHFORIA currently implements an expansion of spherical harmonics method to analytically solve the Laplace equation of WSA model. However, that method can give rise to two effects in the modeling output (Tran, 2009;Tóth et al., 2011). The first one is a ring-like pattern present in low latitudes, that can be visible in the open-closed field maps in Figure 3. This pattern is due to strong isolated pixels associated to active regions. Therefore, if an active region is close to a CH, it can generate the ring-like pattern that can affect the shape of the modeled open field area. The second effect is a ring-like patterns near the poles similar to Gibbs phenomenon in the Fourier analysis. The choice of the number of spherical harmonics used in the model is important, since using higher spherical harmonics can reduce the first effect but enhances the second. Thus, it is necessary to work in improving the model to reduce the ring patterns. Future work should also focus on improving the statistics and also on different aspects of the model outcomes, that is, CH areas, open flux at Earth, velocity of high-speed streams, and current sheet position. 859729 (SWAMI). EUHFORIA is developed as a joint effort between the University of Helsinki and KU Leuven. The validation of solar wind and CME modeling with EUHFORIA is being performed within the BRAIN-be project CCSOM (Constraining CMEs and Shocks by Observations and Modelling throughout the inner heliosphere; http://www.sidc.be/ ccsom/). The results presented here have been achieved under the framework of the Finnish Centre of Excellence in Research of Sustainable Space (Academy of Finland Grant 312390), which we gratefully acknowledge. This work benefited from open access to GONG magnetograms developed with the ADAPT model, and the EUV filtegrams obtained by AIA on board SDO and provided by JSOC. Last but not least, we thank the anonymous reviewers for their valuable comments.