11institutetext: Max-Planck Institut für Sonnensystemforschung, Justus-von-
Liebig-Weg 3, 37077 Göttingen, Germany
11email<EMAIL_ADDRESS>22institutetext: School of Space Research, Kyung Hee
University, Yongin, 446-701 Gyeonggi, Korea 33institutetext: School of Space
and Environment, Beihang University, Beijing 100191, China 44institutetext:
Key Laboratory of Solar Activity, National Astronomical Observatories, Chinese
Academy of Sciences, Beijing 100012, China

# The relationship between bipolar magnetic regions and their sunspots

K. L. Yeo 11 S. K. Solanki 1122 N. A. Krivova 11 J. Jiang 3344

###### Abstract

Context. The relationship between bipolar magnetic regions (BMRs) and their
sunspots is an important property of the solar magnetic field, but it is not
well constrained. One consequence is that it is a challenge for surface flux
transport models (SFTMs) based on sunspot observations to determine the
details of BMR emergence, which they require as input, from such data.

Aims. We aimed to establish the relationship between the amount of magnetic
flux in newly emerged BMRs and the area of the enclosed sunspots, and examine
the results of its application to an established SFTM.

Methods. Earlier attempts to constrain BMR magnetic flux were hindered by the
fact that there is no extensive and reliable record of the magnetic and
physical properties of newly emerged BMRs currently available. We made use of
the empirical model of the relationship between the disc-integrated facular
and network magnetic flux and the total surface coverage by sunspots reported
in a recent study. The structure of the model is such that it enabled us to
establish, from these disc-integrated quantities, an empirical relationship
between the magnetic flux and sunspot area of individual newly emerged BMRs,
circumventing the lack of any proper BMR database.

Results. Applying the constraint on BMR magnetic flux derived here to an
established SFTM retained its key features, in particular its ability to
replicate various independent datasets and the correlation between the model
output polar field at the end of each cycle and the observed strength of the
following cycle. The SFTM output indicates that facular and network magnetic
flux rises with increasing sunspot magnetic flux at a slowing rate such that
it appears to gradually saturate. This is analogous to what earlier studies
comparing disc-integrated quantities sensitive to the amount of faculae and
network present to sunspot indices had reported. The activity dependence of
the ratio of facular and network flux to sunspot flux is consistent with the
findings of recent studies: although the Sun is faculae-dominated (such that
its brightness is mostly positively correlated with activity), it is only
marginally so as facular and network brightening and sunspot darkening appear
to be closely balanced.

###### Key Words.:

Sun: activity — Sun: magnetic fields

##  1 Introduction

The global solar dynamo produces the cyclic emergence of bipolar magnetic
regions (BMRs) on the solar surface which underlies the 11-year activity cycle
(Charbonneau, 2020). Emergent BMRs range from the larger sunspot-bearing
magnetic bipoles, which go on to form active regions, to the smaller spot-free
and relatively short-lived (timescale of hours) ephemeral regions (ERs)
(Harvey, 1993, 2000; Hagenaar, 2001; Hagenaar et al., 2003; Thornton &
Parnell, 2011). How the larger BMRs relate to the sunspots they enclose and
how ER emergence relates to the number of sunspots present are of interest as
they connect directly to how solar magnetism and activity relate to sunspot
prevalence, with implications for the study of a wide range of solar
phenomena.

The large-scale photospheric magnetic field evolves as BMRs emerge, their
magnetic flux is transported and dispersed by systematic and turbulent
motions, and magnetic elements of opposite polarity approach and cancel each
other (van Driel-Gesztelyi & Green, 2015). This process can be described via
the surface flux transport equation (Leighton, 1964; DeVore et al., 1984) or,
at a simpler level, by a closed system of ordinary differential equations
(Solanki et al., 2000, 2002). Models have been developed to reproduce the past
and forecast the future evolution of the large-scale photospheric magnetic
field by solving the transport equation for a given set of initial conditions
and details of BMR emergence, aptly termed surface flux transport models
(SFTMs) (Jiang et al., 2014). Some studies have simulated the evolution over
the past few centuries by prescribing to SFTMs the details of BMR emergence
over this period, as determined from sunspot number or sunspot area records
(e.g. Wang et al., 2005; Cameron et al., 2010; Jiang et al., 2011a, b).
Sunspot number and sunspot area records are employed as they are the only
direct observations of solar surface magnetism to go this far back in time.
Otherwise, SFTMs usually make use of solar magnetograms (e.g. Whitbread et
al., 2018; Virtanen et al., 2018; Yeates, 2020) or Ca II K spectroheliograms
(e.g. Virtanen et al., 2019) instead. While magnetogram and Ca II K
spectroheliogram observations do not extend as far back in time as sunspot
observations, the details of BMR emergence can be established more directly
and accurately from such data.

Since the large-scale photospheric magnetic field is the observable part of
the global dynamo and determines the structure of the coronal magnetic field,
SFTMs represent a powerful tool in the study of these solar phenomena (Cameron
et al., 2012; Mackay & Yeates, 2012; Jiang et al., 2013). Models of solar
irradiance variability are often employed to generate the historical solar
forcing input required by climate simulations (Haigh, 2007) and to investigate
how the brightness variability of the Sun compares to that of other cool stars
(Shapiro et al., 2014, 2015, 2016; Nèmec et al., 2020a, b). Through their
incorporation into models of solar irradiance variability, discussed next,
SFTMs are also relevant to the study of the Earth’s climate and of the solar-
stellar connection.

Solar irradiance variability on timescales greater than a day is believed to
be mainly driven by the variation in the prevalence and distribution of dark
sunspots and bright faculae and network on the solar disc as the photospheric
magnetic field evolves (Solanki et al., 2013; Yeo et al., 2017). Accordingly,
solar irradiance models describe the variation in this quantity due to the
intensity deficit from sunspots and the excess from faculae and network, and
this information is determined from observations of solar magnetism (Domingo
et al., 2009; Yeo et al., 2014a).

The variation in solar irradiance over the past few centuries, a period of
particular interest to climate studies, is often modelled from sunspot number
and sunspot area records (e.g. Lean, 2000; Krivova et al., 2007, 2010; Dasi-
Espuig et al., 2014, 2016; Coddington et al., 2016; Wu et al., 2018).
Modelling solar irradiance variability with sunspot observations means
inferring not just the effect of sunspots, but also the effect of faculae and
network on solar irradiance from information about sunspots alone. This can be
aided by SFTMs aimed at reproducing the evolution of the large-scale
photospheric magnetic field from the same sunspot observations. For example,
Coddington et al. (2016) incorporated into their solar irradiance
reconstruction the decadal variability in the SFTM reported by Wang et al.
(2005). In another instance, the solar irradiance reconstruction by Dasi-
Espuig et al. (2014, 2016) is based on synthetic full-disc magnetograms
generated from sunspot area and sunspot number records using the SFTM of
Cameron et al. (2010) and Jiang et al. (2011a, b) (hereafter CJSS).

Various studies found the facular area to increase with increasing sunspot
area at a diminishing rate such that it appears to gradually saturate (Foukal,
1993, 1996, 1998; Chapman et al., 1997; Shapiro et al., 2014). Solanki &
Fligge (1998) and Yeo et al. (2020) (hereafter YSK), compared various disc-
integrated quantities that are sensitive to the amount of faculae and network
present to sunspot area and/or sunspot number and returned analogous results.
As first pointed out by Foukal (1993), these results suggest that the balance
between the competing effects of bright faculae and network and of dark
sunspots on solar irradiance might shift towards the latter as we go from low
to high activity levels. This is believed to extend to cool stars in general,
such that as we go from less to more active stars, we transit from a faculae-
dominated regime where facular and network brightening dominates starspot
darkening to a spot-dominated regime where the converse is true. This explains
why less active stars appear to become brighter with increasing activity,
while more active stars exhibit the opposite trend (Lockwood et al., 2007;
Hall et al., 2009; Shapiro et al., 2014; Radick et al., 2018). The balance
between bright faculae and network and dark sunspots on the solar disc, and
how that might change with the activity level and the inclination of the Sun
to the observer, is of interest for the relevance to the understanding of the
brightness variability of the Sun and other cool stars. In recent studies
Nèmec et al. (2020a, b) extended the model of solar irradiance variability by
Dasi-Espuig et al. (2014, 2016), which incorporates the SFTM by CJSS, to probe
the inclination-dependence of this balance.

In order to reproduce the evolution of the large-scale photospheric magnetic
field over a given period of time, SFTMs require as input the physical
properties of each BMR that had emerged over this time interval, including the
timing, location, and amount of magnetic flux at emergence. However, deriving
this information from observations of sunspots is less than straightforward.
As ERs do not contain or emerge synchronously with sunspots and are difficult
to observe due to their relatively small sizes and short lifetimes, the
relationship between ER emergence and sunspot prevalence is not exactly known.
As for sunspot-bearing BMRs, it is less than straightforward to infer the
details of their emergence from sunspot number records as these are
essentially sunspot and sunspot group counts, and therefore do not contain
information on individual sunspot groups. While sunspot area records do
tabulate the location and area of individual sunspot groups, how the amount of
magnetic flux in a given newly emerged sunspot-bearing BMR, $\Phi_{\rm BMR}$,
depends on the area of the enclosed sunspots, $A_{\rm BMR,S}$, is not well
constrained. By the term ‘newly emerged’ we refer to the point where the
sunspot-bearing BMR has fully emerged but not yet started to decay. This paper
focuses on the issue of the relationship between $\Phi_{\rm BMR}$ and $A_{\rm
BMR,S}$. For the sake of brevity, hereafter the term BMR denotes sunspot-
bearing BMRs alone, excluding ERs.

The difficulty in establishing BMR emergence from sunspot observations leaves
SFTMs to make varying assumptions. The Wang et al. (2005) model, based on the
Hoyt & Schatten (1998) sunspot group number record, accommodated the fact that
there is no information on individual sunspot groups in such data by assigning
the same amount of magnetic flux to each of the emergent BMRs in a given solar
cycle. In the CJSS model, detailed in Sect. 3.1, the $A_{\rm BMR,S}$ of each
BMR is taken directly from sunspot area records. We let $\Delta\beta$ denote
the separation between the two poles of a BMR. The model combines an empirical
relationship between $\Delta\beta$ and $A_{\rm BMR,S}$ and the assumption that
$\Phi_{\rm BMR}$ is linearly proportional to $\left(\Delta\beta\right)^{2}$ to
determine $\Phi_{\rm BMR}$ from $A_{\rm BMR,S}$.

Few studies have attempted to constrain the relationship between $\Phi_{\rm
BMR}$ and $A_{\rm BMR,S}$. From an examination of the faculae in active
regions and newly emerged BMRs, Schrijver & Harvey (1994) postulated that
$\Phi_{\rm BMR}$ scales linearly with the total area extended by the faculae
and sunspots in each BMR, $A_{\rm BMR}$. This can be combined with the
empirical relationship between facular area and sunspot area (Foukal, 1993,
1996, 1998; Chapman et al., 1997; Shapiro et al., 2014) to roughly relate
$\Phi_{\rm BMR}$ to $A_{\rm BMR,S}$. In a more recent study, Lemerle et al.
(2015) compared the database of newly emerged BMRs compiled by Wang & Sheeley
(1989) to the USAF-NOAA sunspot area record and found $\Phi_{\rm BMR}$ to
conform to a power-law relationship with $A_{\rm BMR,S}$ with a power-law
exponent of 0.71. The analyses of Schrijver & Harvey (1994) and Lemerle et al.
(2015) (hereafter SH and LCC), were restricted by certain limitations of the
observations available to them, discussed in Sect. 2.3. It would be
straightforward to derive the relationship between $\Phi_{\rm BMR}$ and
$A_{\rm BMR,S}$ from an extensive and reliable record of the magnetic and
physical properties of newly emerged BMRs. There is, however, no such record
currently available, as also discussed in Sect. 2.3.

In this study we aim to establish the relationship between the amount of
magnetic flux in each newly emerged BMR, $\Phi_{\rm BMR}$, and the area of the
enclosed sunspots, $A_{\rm BMR,S}$. We make use of the empirical model
relating the disc-integrated facular and network magnetic flux,
$\Sigma\Phi_{\rm FN}$, to the total surface coverage by sunspots,
$\Sigma{}A_{\rm S}$ presented by YSK (Yeo et al., 2020), which builds on the
earlier work by Preminger & Walton (2005, 2006a, 2006b, 2007), abbreviated
here as PW. As we demonstrate, the set-up of the YSK model allows us to
establish from it an empirical relationship between $\Phi_{\rm BMR}$ and
$A_{\rm BMR,S}$ while circumventing the lack of a proper database of newly
emerged BMRs. In the following, we describe the YSK model (Sect. 2.1), derive
the empirical relationship between $\Phi_{\rm BMR}$ and $A_{\rm BMR,S}$ (Sect.
2.2), and compare it to the SH and LCC results (Sect. 2.3). Then, we apply
each of these constraints on $\Phi_{\rm BMR}$ to the SFTM by CJSS (Sect. 3.1)
and examine the effect on the model output (Sect. 3.2). Making use of the
surface flux transport simulation based on the empirical relationship between
$\Phi_{\rm BMR}$ and $A_{\rm BMR,S}$ derived here, we will look at the balance
between the amount of faculae and network and the amount of sunspots on the
solar disc, and how that changes with solar activity level (Sect. 4). Finally,
we provide a summary of the study in Sect. 5.

##  2 The relationship between the amount of magnetic flux in newly emerged
bipolar magnetic regions and the area of the enclosed sunspots

###  2.1 Yeo et al. (2020) model

In this section we derive an empirical relationship between the amount of
magnetic flux in each newly emerged bipolar magnetic region or BMR, $\Phi_{\rm
BMR}$, and the area of the enclosed sunspots, $A_{\rm BMR,S}$. For this
purpose we employ the empirical model of the relationship between the disc-
integrated facular and network magnetic flux, $\Sigma\Phi_{\rm FN}$, and the
total surface coverage by sunspots, $\Sigma{}A_{\rm S}$, recently reported by
YSK (Yeo et al., 2020). The YSK model is given by

$\Sigma\Phi_{\rm FN}=\left(\Sigma{}A_{\rm S}\right)^{h_{1}}\otimes{}H+h_{2},$
(1)

where $h_{1}$ and $h_{2}$ are fit parameters. The convolution of
$\left(\Sigma{}A_{\rm S}\right)^{h_{1}}$ with $H$, a finite impulse response
(FIR) filter, describes the variation in $\Sigma\Phi_{\rm FN}$ associated with
sunspot-bearing BMRs. The $h_{2}$ term represents network magnetic flux that
did not originate from sunspot-bearing BMRs; it is assumed in the model to be
constant. The FIR filter is given by

$H\left(t\right)=h_{3}\max\left[\exp\left(-\frac{|t|}{h_{4}}\right)\cos\left(\frac{2\pi{}t}{t_{\sun}}\right),0\right],$
(2)

where $-0.25t_{\sun}\leq{}t\leq 14.25t_{\sun}$, $t_{\sun}=26.24\ {\rm days}$
denotes the synodic rotation period of the Sun, and $h_{3}$ and $h_{4}$ are
fit parameters.

Yeo et al. (2020) extended the composite time series of daily $\Sigma{}A_{\rm
S}$ by Balmaceda et al. (2009) (Fig. 1a) and of daily $\Sigma\Phi_{\rm FN}$ by
Yeo et al. (2014b) (Fig. 1b) to June 2017 and constrained the fit parameters
with the result. Yeo et al. (2014b) had derived daily $\Sigma\Phi_{\rm FN}$
back to 1974 by taking daily full-disc magnetograms from the Kitt Peak Vacuum
Telescope (Livingston et al., 1976; Jones et al., 1992), SoHO/MDI (Scherrer et
al., 1995), and SDO/HMI (Scherrer et al., 2012), accounting for the systematic
differences between the observations from the various telescopes due to
instrumental factors and isolating the faculae and network present in each
magnetogram. YSK could not extend the Balmaceda et al. (2009) time series
beyond June 2017 as it is based on USAF-NOAA sunspot areas for the period of
1986 on, and these measurements are available only up to this time. For
$\Sigma{}A_{\rm S}$ in millionths of the solar hemisphere ($\mu$Hem) and
$\Sigma\Phi_{\rm FN}$ in units of Weber (Wb), $h_{1}=0.843$, $h_{2}=1.82\times
10^{14}{\rm\ Wb}$, $h_{3}=1.99\times{}10^{11}{\rm\ Wb}$, and $h_{4}=33.5{\rm\
days}$. The robustness of the model is visibly evident in Fig. 1b; the model
reconstruction of $\Sigma\Phi_{\rm FN}$ from $\Sigma{}A_{\rm S}$ (red) is very
similar to the Yeo et al. (2014b) time series (black), and the 365-day running
mean of the difference between the two (blue) is close to null and does not
indicate any clear systematic departures. The Pearson’s correlation
coefficient, $R$, between the two time series is 0.950, meaning the model
reproduces about $90\%$ of the observed variability. We refer the reader to
Yeo et al. (2020) for a detailed discussion of the uncertainties in the model
and how it compares to earlier models.

Figure 1: The observed $\Sigma{}A_{\rm S}$ and $\Sigma\Phi_{\rm FN}$, the YSK
model (Equation 1) reconstruction of $\Sigma\Phi_{\rm FN}$ from
$\Sigma{}A_{\rm S}$ and the $H$ term in this model (Equation 2). a) Composite
time series of the total surface coverage by sunspots, $\Sigma{}A_{\rm S}$, by
Balmaceda et al. (2009). b) Composite time series of the disc-integrated
facular and network magnetic flux, $\Sigma\Phi_{\rm FN}$, by Yeo et al.
(2014b) (black), the reconstruction of this quantity from $\Sigma{}A_{\rm S}$
with the YSK model (red), and the 365-day running mean of the difference
between measurement and model (blue). The Balmaceda et al. (2009) time series
goes back to 1874 but the range of the plot is confined to 1974 to 2017, the
period where it overlaps with the Yeo et al. (2014b) time series. The gaps in
the Yeo et al. (2014b) time series correspond to periods where there no
suitable magnetogram observations with which to determine $\Sigma\Phi_{\rm
FN}$ were available. c) The $H$ term, which describes the response of
$\Sigma\Phi_{\rm FN}$ to the appearance of a sunspot per unit sunspot area.
The dotted line follows the exponential envelope of $H$ and $t_{\sun}$ denotes
the synodic rotation period of the Sun (see Sect. 2.1 for details).

The FIR filter $H$ (Equation 2) is based on the earlier findings of PW
(Preminger & Walton, 2005, 2006a, 2006b, 2007). Sunspots form some time after
the BMRs they belong to have started emerging on the solar surface; they are
short-lived compared to active regions and their decay products, such that the
latter persists long after their sunspots have decayed (van Driel-Gesztelyi &
Green, 2015). Consequently, the response of any quantity that is sensitive to
the presence of BMRs, active regions, and their decay products, including
$\Sigma\Phi_{\rm FN}$, to the appearance of a sunspot similarly starts some
time before it appears and persists for a period after it dissipates. PW
deconvolved $\Sigma{}A_{\rm S}$ from the average unsigned magnetic flux
density over the solar disc, solar irradiance, and various chromospheric and
coronal indices, deriving FIR filters that elucidate the response of each
quantity to the appearance of a sunspot and how this response varies with time
due to the phenomenon just described. Preminger & Walton (2007) captured the
key features of these empirical FIR filters with a model FIR filter (Equation
2 in their paper), which YSK modified slightly to more accurately reflect the
form of the empirical FIR filters, yielding $H$. Depicted in Fig. 1c, $H$
describes a series of sinusoidal lobes, $0.5t_{\sun}$ wide and $t_{\sun}$
apart, of diminishing amplitude, as delineated by the exponential envelope of
time constant $h_{4}$ (dotted line). The first lobe has an amplitude of
$h_{3}$ and is centred on $t=0$. The sequence of lobes is truncated at
$t=14.25t_{\sun}$. The FIR filter charts the amount of facular and network
magnetic flux associated with a sunspot-bearing BMR, per unit sunspot area,
from $0.25t_{\sun}$ before it fully emerges to $14.25t_{\sun}$ after. The
rising edge of the first lobe corresponds to the BMR emerging, the peak of the
first lobe to full emergence, and the falling edge and the subsequent lobes to
its decay and modulation with solar rotation.

If the amount of facular and network magnetic flux associated with a sunspot-
bearing BMR scales linearly with the area of the enclosed sunspots, then the
variation in $\Sigma\Phi_{\rm FN}$ due to sunspot-bearing BMRs would be given
by $\Sigma{}A_{\rm S}\otimes{}H$. However, as noted in the introduction,
various studies have found the facular area to increase with increasing
sunspot area at a diminishing rate, such that it is at a lower level than if
the two quantities are linearly related (Foukal, 1993, 1996, 1998; Chapman et
al., 1997; Shapiro et al., 2014). This implies that the larger an active
region, the higher (lower) the proportion of the total area occupied by
sunspots (faculae). Based on this observation, YSK argued that it is unlikely
that the amount of facular and network magnetic flux associated with a
sunspot-bearing BMR conforms to a linear relationship with the area of the
enclosed sunspots. They introduced the fit parameter $h_{1}$, such that the
variation in $\Sigma\Phi_{\rm FN}$ is given by $\left(\Sigma{}A_{\rm
S}\right)^{h_{1}}\otimes{}H$ instead. This and the modifications introduced to
the model FIR filter reported by PW produced a closer agreement between
modelled and measured $\Sigma\Phi_{\rm FN}$, detailed in Yeo et al. (2020).

After the time of the YSK study, the $\Sigma{}A_{\rm S}$ composite by
Balmaceda et al. (2009) was revised by Mandal et al. (2020). Both the
Balmaceda et al. (2009) and Mandal et al. (2020) time series are based on, up
to 1976, Royal Greenwich Observatory (RGO) sunspot area measurements. They are
therefore exactly identical to one another up to this time, and differ only in
the observations used to extend the RGO time series to more recent times.
Mandal et al. (2020) noted that these two time series differ from one another
by about 6% between 1977 and 1985, and that there are otherwise no noticeable
systematic differences between the two. The revision of the YSK model with the
Mandal et al. (2020) time series is outside the scope of the current study. In
any case, we do not expect the result of applying the Mandal et al. (2020)
time series to the YSK model to differ significantly from what was obtained
with the Balmaceda et al. (2009) time series as input. The disparity between
these two time series between 1977 and 1985 is relatively small, especially in
comparison to the random error in sunspot area measurements (20% to 50%, Sofia
et al., 1982). In addition, less than 10% of the data points in the daily
$\Sigma\Phi_{\rm FN}$ composite by Yeo et al. (2014b) lie in this period.

###  2.2 Derivation

The total magnetic flux in each newly emerged BMR is the sum of the magnetic
flux in its faculae and in its sunspots, or

$\Phi_{\rm BMR}=\Phi_{\rm BMR,F}+\Phi_{\rm BMR,S}.$ (3)

Here we establish how $\Phi_{\rm BMR,F}$ and $\Phi_{\rm BMR,S}$ relate to
$A_{\rm BMR,S}$ in turn. To derive the relationship between $\Phi_{\rm BMR,F}$
and $A_{\rm BMR,S}$, let us isolate from the YSK model the contribution to
$\Sigma\Phi_{\rm FN}$ by the faculae in newly emerged BMRs, $\Sigma\Phi_{\rm
BMR,F}$. As noted in Sect. 2.1, the peak of the first lobe of the FIR filter,
which is at $t=0$, corresponds to the response of $\Sigma\Phi_{\rm FN}$ to a
sunspot-bearing BMR when it has fully emerged and has not started to decay.
Taking Equation 1, we isolate $\Sigma\Phi_{\rm BMR,F}$ by replacing $H$ with
$H\left(t=0\right)=h_{3}$ and factoring out the basal level, $h_{2}$, arriving
at

$\Sigma\Phi_{\rm BMR,F}=h_{3}\left(\Sigma{}A_{\rm S}\right)^{h_{1}}.$ (4)

The PW and YSK approach, by modelling the sunspot-related variation in a given
disc-integrated quantity as the convolution of $\Sigma{}A_{\rm S}$ with a FIR
filter and with the various quantities at daily cadence, treats the appearance
of sunspots as one-day events. In other words, it implicitly assumes that each
sunspot exists for just one day, such that the total surface coverage by
sunspots, $\Sigma{}A_{\rm S}$, is treated in the model as the total surface
coverage by sunspots in newly emerged BMRs, $\Sigma{}A_{\rm BMR,S}$. Replacing
$\Sigma{}A_{\rm S}$ in Equation 4 with $\Sigma{}A_{\rm BMR,S}$, we get

$\Sigma\Phi_{\rm BMR,F}=h_{3}\left(\Sigma{}A_{\rm BMR,S}\right)^{h_{1}}.$ (5)

Assuming this relationship applies to individual BMRs, then

$\displaystyle\begin{split}\Phi_{\rm BMR,F}&=h_{3}A_{\rm BMR,S}^{h_{1}}\\\
&=1.99\times{}10^{11}A_{\rm BMR,S}^{0.843}{\rm\ Wb}.\end{split}$ (6)

Cho et al. (2015) examined how the amount of magnetic flux in sunspots and
pores varies with their size. The authors noted that for fully formed
sunspots, the amount of magnetic flux appears to scale linearly with its area,
confirming the early result by Sheeley (1966). Assuming the mean magnetic flux
density over sunspots is 1350 G or $4.11\times{}10^{11}\ {\rm Wb}\ \mu{\rm
Hem^{-1}}$ (Solanki & Schmidt, 1993; Lites et al., 1993; Keppens & Martinez
Pillet, 1996), this means

$\Phi_{\rm BMR,S}=4.11\times{}10^{11}A_{\rm BMR,S}\ {\rm Wb}.$ (7)

Taking Equations 6 and 7 into Equation 3, we get

$\displaystyle\begin{split}\Phi_{\rm BMR}&=4.11\times{}10^{11}A_{\rm
BMR,S}+1.99\times{}10^{11}{}A_{\rm BMR,S}^{0.843}\ {\rm Wb}\\\
&=4.11\times{}10^{11}\left(A_{\rm BMR,S}+0.485A_{\rm BMR,S}^{0.843}\right)\
{\rm Wb}\\\ &=4.11\times{}10^{11}\tilde{A}_{\rm BMR}\ {\rm Wb},\end{split}$
(8)

which relates the amount of magnetic flux in newly emerged BMRs to the area of
the enclosed sunspots. The $\tilde{A}_{\rm BMR}$ term, representing $A_{\rm
BMR,S}+0.485A_{\rm BMR,S}^{0.843}$, denotes for a given BMR the area of the
sunspot that would contain the same amount of magnetic flux (compare Equations
7 and 8). It is worth noting here that the YSK model (Sect. 2.1) and the
analysis based on this model that we have just presented pertains only to
sunspot-bearing BMRs and does not apply to ephemeral regions.

In Fig. 2 we chart $\Phi_{\rm BMR,S}$, $\Phi_{\rm BMR}$ (dashed and solid
black lines, top panel), and the ratio of the two, $\frac{\Phi_{\rm
BMR,S}}{\Phi_{\rm BMR}}$ (black line, bottom panel) as a function of $A_{\rm
BMR,S}$. As is visibly evident, the empirical relationship between $\Phi_{\rm
BMR}$ and $A_{\rm BMR,S}$ presented in Equation 8 indicates that most of the
magnetic flux in newly emerged BMRs is in the sunspots they enclose. The
quantity $\frac{\Phi_{\rm BMR,S}}{\Phi_{\rm BMR}}$ increases gradually from
about 0.7 at $A_{\rm BMR,S}=1\ \mu{}{\rm Hem^{-1}}$ to 0.9 at $6000\ \mu{}{\rm
Hem^{-1}}$ (left axis, Fig. 2b). Schrijver (1987) established an empirical
relationship between the proportion of the magnetic flux in active regions
contained in their sunspots and the proportion of the area of active regions
occupied by their sunspots. Assuming it applies to newly emerged BMRs, then

$\frac{\Phi_{\rm BMR,S}}{\Phi_{\rm BMR}}=5.1\frac{A_{\rm BMR,S}}{A_{\rm
BMR}},$ (9)

meaning the trend in $\frac{\Phi_{\rm BMR,S}}{\Phi_{\rm BMR}}$ with $A_{\rm
BMR,S}$ (left axis, Fig. 2b) is accompanied by $\frac{A_{\rm BMR,S}}{A_{\rm
BMR}}$ rising in proportion from about 0.14 to 0.18 (right axis). In other
words, the smallest newly emerged BMRs are about 70$\%$ sunspots by magnetic
flux and 14$\%$ sunspots by area, and this increases steadily to about 90$\%$
and 18$\%$, respectively for the largest newly emerged BMRs.

Figure 2: The relationship between BMR magnetic flux and sunspot area. a)
Amount of magnetic flux in newly emerged BMRs, $\Phi_{\rm BMR}$, as a function
of the area of the enclosed sunspots, $A_{\rm BMR,S}$, as indicated by the
empirical relationship between $\Phi_{\rm BMR}$ and $A_{\rm BMR,S}$ presented
in this study (black solid line, Equation 8) and the constraints on $\Phi_{\rm
BMR}$ reported by SH (blue, Equation 10) and LCC (red, Equation 13). For the
estimate of $\Phi_{\rm BMR}$ based on the SH result, the lower and upper
bounds, corresponding to $k=1$ and $k=2$, are given. The dashed line follows
the amount of magnetic flux in newly emerged BMRs contained in the embedded
sunspots, $\Phi_{\rm BMR,S}$ (Equation 7). b) Corresponding plot of
$\frac{\Phi_{\rm BMR,S}}{\Phi_{\rm BMR}}$. To an approximation,
$\frac{\Phi_{\rm BMR,S}}{\Phi_{\rm BMR}}$ is directly proportional to the
fraction of the area of each BMR occupied by their sunspots, $\frac{A_{\rm
BMR,S}}{A_{\rm BMR}}$, indicated on the right axis (Equation 9) (see Sects.
2.2 and 2.3 for details).

We note the following. Firstly, active regions and newly emerged BMRs are not
equivalent to one another, meaning the application of the findings of
Schrijver (1987) to the latter (Equation 9) is not strictly correct. Secondly,
some of the smallest BMRs might contain pores instead of fully developed
sunspots (e.g. Švanda et al., 2016), such that the amount of sunspot magnetic
flux they enclose might be smaller than indicated by Equation 7. Nonetheless,
the assertion here that the sunspots in a given newly emerged BMR can enclose
most of its magnetic flux while occupying only a small proportion of its area
is reasonable and within expectation given that sunspots and pores are
magnetically denser than faculae.

###  2.3 Comparison to Schrijver & Harvey (1994) and Lemerle et al. (2015)

In this section we compare the empirical relationship between $\Phi_{\rm BMR}$
and $A_{\rm BMR,S}$, derived here (Equation 8) to the constraints on
$\Phi_{\rm BMR}$ reported by SH (Schrijver & Harvey, 1994) and by LCC (Lemerle
et al., 2015). Taking Kitt Peak magnetograms, SH isolated the faculae in
active regions and newly emerged BMRs and found the facular magnetic flux in
individual regions to scale linearly with the area extended by these features.
Based on this result and assuming an average BMR magnetic flux density of 100
G or $1.52\times{}10^{10}\ {\rm Wb}\ \mu{\rm Hem^{-1}}$, SH estimated that

$\Phi_{\rm BMR}=1.52\times{}10^{10}kA_{\rm BMR}\ {\rm Wb},$ (10)

where $k$ is a free parameter that the authors did not constrain. SH
interpreted the quantity $k-1$ for a given newly emerged BMR as the sum of
$\Phi_{\rm BMR,S}$, and the magnetic flux that dissipated before the BMR had
fully formed as a proportion of $\Phi_{\rm BMR}$. Since this ratio is unlikely
to exceed unity, at least not by a wide margin, $k$ most probably lies between
one and two. Instead of establishing the relationship between $\Phi_{\rm BMR}$
and $A_{\rm BMR}$ directly, SH had estimated it from the relationship between
facular magnetic flux and area. This was necessary because, as noted in
precursor studies to the SH study conducted by the same authors (Schrijver,
1987; Harvey & Zwaan, 1993), sunspot magnetic flux cannot be taken directly,
at least not reliably, from Kitt Peak magnetograms.

Ideally, to relate $\Phi_{\rm BMR}$ to $A_{\rm BMR,S}$ using the relationship
between $\Phi_{\rm BMR}$ and $A_{\rm BMR}$ reported by SH (Equation 10), it
should be combined with the relationship between $A_{\rm BMR}$ and $A_{\rm
BMR,S}$. Unfortunately, the relationship between $A_{\rm BMR}$ and $A_{\rm
BMR,S}$ is not known. Instead, let us assume the empirical relationship
between the total surface coverage by faculae and by sunspots reported by
Chapman et al. (1997) (hereafter CCD) can be taken as the relationship between
the area of the faculae and of the sunspots in individual newly emerged BMRs,
such that

$A_{\rm BMR}=414+22A_{\rm BMR,S}-0.0036A_{\rm BMR,S}^{2}.$ (11)

Inserting Equation 11 into Equation 10, we obtain

$\Phi_{\rm BMR}=1.52\times{}10^{10}k\left[414+22A_{\rm BMR,S}-0.0036A_{\rm
BMR,S}^{2}\right]\ {\rm Wb}.$ (12)

The estimate of $\Phi_{\rm BMR}$ from this expression at $k=1$ and $k=2$
represents the lower and upper bound of this quantity based on the combination
of the SH and CCD findings.

Lemerle et al. (2015) compared the $\Phi_{\rm BMR}$ measurements in the
database of newly emerged BMRs compiled by Wang & Sheeley (1989) to USAF-NOAA
sunspot areas. They could not match each newly emerged BMR in the Wang &
Sheeley (1989) database to the corresponding sunspot group in the USAF-NOAA
record, which excluded a direct comparison and compelled them to compare the
distribution of values in the two records instead. From this, LCC estimated
that

$\Phi_{\rm BMR}=1.12\times{}10^{12}A_{\rm BMR,S}^{0.714}\ {\rm Wb}.$ (13)

Taking Equations 12 and 13, we recalculate $\Phi_{\rm BMR}$ and
$\frac{\Phi_{\rm BMR,S}}{\Phi_{\rm BMR}}$, drawn in Fig. 2 (blue and red
lines) along the values from the empirical relationship between $\Phi_{\rm
BMR}$ and $A_{\rm BMR,S}$ derived in this study (black solid lines).

Figure 2a shows, at lower levels of $A_{\rm BMR,S}$, that the estimates of
$\Phi_{\rm BMR}$ based on the SH (blue) and LCC results (red) are higher than
the current estimate (black solid line). However, $\Phi_{\rm BMR}$ increases
with increasing $A_{\rm BMR,S}$ at a slower rate in the former, such that at
higher levels of $A_{\rm BMR,S}$ the estimates of $\Phi_{\rm BMR}$ based on
the SH and LCC results drop below the current estimate and also below
$\Phi_{\rm BMR,S}$ (dashed line). This translates into the analogous behaviour
of $\frac{\Phi_{\rm BMR,S}}{\Phi_{\rm BMR}}$ illustrated in Fig. 2b. At low
$A_{\rm BMR,S}$ the estimates of $\frac{\Phi_{\rm BMR,S}}{\Phi_{\rm BMR}}$
based on the SH (blue) and LCC results (red) are lower than the current
estimate (black), but they rise more rapidly with increasing $A_{\rm BMR,S}$
such that they eventually exceed not just the current estimate, but also
unity. The amount of magnetic flux in a given BMR has to be greater than the
amount carried by the sunspots it encloses, meaning $\Phi_{\rm BMR}$ going
below $\Phi_{\rm BMR,S}$ or, equivalently, $\frac{\Phi_{\rm BMR,S}}{\Phi_{\rm
BMR}}$ going above unity is unphysical. The observation here that the
estimates of $\Phi_{\rm BMR}$ based on the SH and LCC results drop below
$\Phi_{\rm BMR,S}$ towards higher $A_{\rm BMR,S}$ is unlikely to be from any
uncertainty in the latter. The margin is so wide that even if the mean
magnetic flux density over sunspots is only one-quarter the value 1350 G we
adopt here, $\Phi_{\rm BMR,S}$ will still be greater than the estimates of
$\Phi_{\rm BMR}$ based on the SH and LCC results for the largest BMRs. The
constraint on $\Phi_{\rm BMR}$ reported by SH, at least when combined with the
findings of CCD as we do here (Equation 11), and that by LCC appear to grossly
underestimate the amount of magnetic flux in larger BMRs.

The discrepancies between the results reported in this study and earlier by SH
and LCC that we have just highlighted could have arisen from the following. SH
examined active regions and newly emerged BMRs together even though the two
are not equivalent. Furthermore, SH had to constrain $\Phi_{\rm BMR}$
indirectly as sunspot magnetic flux cannot be determined directly and reliably
from the magnetogram observations they utilised. How the SH results compare to
the LCC and current results could have also been affected by the fact that we
had to combine the SH result with the findings of CCD in order to make this
comparison, and the latter does not strictly apply to newly emerged BMRs.

As for the LCC study, the authors could not compare the Wang & Sheeley (1989)
and USAF-NOAA records directly because they could not match each newly emerged
BMR in the former to the corresponding sunspot group in the latter. The Wang &
Sheeley (1989) database includes spot-free BMRs, which are not represented in
the USAF-NOAA record and should have been excluded. Wang & Sheeley (1989) had
also determined $\Phi_{\rm BMR}$ approximately from photographic magnetograms,
introducing further uncertainty to the LCC analysis. SH and LCC were hampered
by limitations in the data available to them. As stated in the introduction,
it would be straightforward to derive the relationship between $\Phi_{\rm
BMR}$ and $A_{\rm BMR,S}$ from an extensive and reliable database of the
magnetic and physical properties of newly emerged BMRs, but no such record
currently exists. The Wang & Sheeley (1989) database is clearly not ideal for
this purpose. The BMR databases compiled more recently by Yeates et al. (2007)
and by Yeates (2020), although of excellent quality, are not suitable either
as they tabulate the physical properties of each BMR at the time it crosses
the central meridian, not when it newly emerges.

The YSK model and the analysis based on it presented in this study have their
own assumptions and limitations, discussed in Sects. 2.1 and 2.2. Nonetheless,
it is worth noting that the YSK model allowed us to circumvent the lack of any
proper database of newly emerged BMRs, an issue that severely limited the SH
and LCC analyses, and constrains the relationship between $\Phi_{\rm BMR}$ and
$A_{\rm BMR,S}$ from $\Sigma\Phi_{\rm FN}$ and $\Sigma{}A_{\rm S}$ (i.e. disc-
integrated observations) instead.

##  3 Application to a surface flux transport model

In this section we look at the result of applying the empirical relationship
between $\Phi_{\rm BMR}$ and $A_{\rm BMR,S}$ presented in this study (Equation
8) and the constraints on $\Phi_{\rm BMR}$ reported by SH (Equation 10) and by
LCC (Equation 13) to the SFTM by CJSS. We generated three surface flux
transport simulation runs, denoted CJSS-SH, CJSS-LCC, and CJSS-YEA (Table 1).
The CJSS-SH run is generated with the CJSS model as it is and set up following
the ‘reference case’ presented by Cameron et al. (2010). As we explain below,
the CJSS model (incidentally and not by design) already complies with the
constraint on $\Phi_{\rm BMR}$ reported by SH. The CJSS-LCC and CJSS-YEA runs
are generated with the CJSS model modified to comply with the constraints on
$\Phi_{\rm BMR}$ reported by LCC and derived here instead. In the following,
we describe the CJSS model and the simulation set-up (Sect. 3.1) before
discussing the results (Sect. 3.2).

Table 1: Surface flux transport simulation runs generated in the current
study.

Simulation run | $\tilde{B}_{\rm max}\ [{\rm G}\ \mu{}{\rm Hem^{-1}}]$ | $f$ | Constraint on $\Phi_{\rm BMR}$ | Reference
---|---|---|---|---
CJSS-SH | 0.0156 | $A_{\rm BMR}$ | Equation 10 | Schrijver & Harvey (1994)
CJSS-LCC | 2.19 | $A_{\rm BMR,S}^{0.71}$ | Equation 13 | Lemerle et al. (2015)
CJSS-YEA | 0.223 | $\tilde{A}_{\rm BMR}$ | Equation 8 | Current study

111The value of $\tilde{B}_{\rm max}$, the definition of $f$ (Equation 18),
and the constraint on $\Phi_{\rm BMR}$ applied to the run by the definition of
$f$ are listed. The total area of each newly emerged BMR, $A_{\rm BMR}$, and
the area of the sunspot that would contain the same amount of magnetic flux,
$\tilde{A}_{\rm BMR}$, are related to the area of the enclosed sunspots,
$A_{\rm BMR,S}$, via Equations 11 and 8, respectively.

###  3.1 Simulation set-up

The governing equation of the CJSS model is

$\displaystyle\begin{split}\frac{\partial{}B}{\partial{}t}=&-\Omega\left(\lambda\right)\frac{\partial{}B}{\partial\phi}-\frac{1}{R_{\odot}\cos\lambda}\frac{\partial}{\partial\lambda}\left[\upsilon\left(\lambda\right)B\cos\lambda\right]\\\
&+\eta_{\rm
S}\left[\frac{1}{R^{2}_{\odot}\cos\lambda}\frac{\partial}{\partial\lambda}\left(\cos\lambda\frac{\partial{}B}{\partial\lambda}\right)+\frac{1}{R^{2}_{\odot}\cos^{2}\lambda}\frac{\partial^{2}B}{\partial\phi^{2}}\right]\\\
&+D\left(\eta_{\rm R}\right)+S\left(\lambda,\phi,t\right),\end{split}$ (14)

where $R_{\sun}$, $\lambda$, and $\phi$ denote the radius of the Sun,
latitude, and longitude, respectively. This describes the passive transport of
the radial component of the photospheric magnetic field, $B$, under the
influence of differential rotation, $\Omega$; meridional flow, $\upsilon$; and
diffusion, which is characterised by the surface and radial diffusivity,
$\eta_{\rm S}$ and $\eta_{\rm R}$ (DeVore et al., 1984; Wang et al., 1989;
Mackay et al., 2000; Schrijver et al., 2002; Baumann et al., 2004). Radial
diffusion is introduced through a linear operator, $D$ (Baumann et al., 2006),
and BMR emergence through a source term, $S$. As noted in the previous
paragraph, we set up the simulation runs following the reference case
presented by Cameron et al. (2010). We refer the reader to this earlier paper
for details on the definition of the initial condition and the $\Omega$,
$\upsilon$, $\eta_{\rm S}$, and $\eta_{\rm R}$ parameters. For the CJSS-SH
run, we set up the source term following Cameron et al. (2010). For the CJSS-
LCC and CJSS-YEA runs, we modified the source term in such a way as to conform
the model to the constraints on $\Phi_{\rm BMR}$ reported by LCC (Equation 13)
and derived here (Equation 8). We describe the definition of the source term
and these modifications next.

The source term (i.e. details of BMR emergence) is inferred from the composite
of the RGO and USAF-NOAA sunspot area records by Hathaway & Choudhary (2008).
The location and area of each sunspot group when the latter is at its maximum
is taken to denote the location and $A_{\rm BMR,S}$ of a newly emerged BMR.
Following van Ballegooijen et al. (1998), each BMR is described as two
circular magnetic patches of opposite polarity. Let
$B_{+}\left(\lambda,\phi\right)$ represent the magnetic flux density of the
positive polarity patch at $\left(\lambda,\phi\right)$,
$\left(\lambda_{+},\phi_{+}\right)$ the centre of the patch, and
$\beta_{+}\left(\lambda,\phi\right)$ the heliocentric angle between
$\left(\lambda,\phi\right)$ and $\left(\lambda_{+},\phi_{+}\right)$, and let
$B_{-}\left(\lambda,\phi\right)$, $\left(\lambda_{-},\phi_{-}\right)$ and
$\beta_{-}\left(\lambda,\phi\right)$ represent the same for the negative
polarity patch. The magnetic flux density profile is given by

$B_{\pm}\left(\lambda,\phi\right)=B_{\rm
max}\left(\frac{0.4\Delta\beta}{\delta}\right)^{2}\exp\left(-\frac{2\left[1-\cos\beta_{\pm}\left(\lambda,\phi\right)\right]}{\delta^{2}}\right),$
(15)

where
$\Delta\beta=\beta_{+}\left(\lambda_{-},\phi_{-}\right)=\beta_{-}\left(\lambda_{+},\phi_{+}\right)$
denotes the heliocentric angle between the two patches and $\delta=4^{\circ}$
the size of each patch. The constant factor $B_{\rm max}$ modulates the
absolute scale of the model output and is fixed at 374 G by matching the
output average unsigned magnetic flux density to the observed values from the
Mount Wilson Observatory (MWO) and Wilcox Solar Observatory (WSO) (Arge et
al., 2002), illustrated in Fig. 3a. We note that BMR magnetic flux density,
$B_{\pm}$, and therefore the amount of magnetic flux in each BMR, $\Phi_{\rm
BMR}$, are assumed to scale linearly with the square of the angular
separation, $\left(\Delta\beta\right)^{2}$. The angular separation,
$\Delta\beta$, is given by

$\Delta\beta=6.46\times 10^{-2}\sqrt{A_{\rm BMR}},$ (16)

where $\Delta\beta$ and $A_{\rm BMR}$ are in units of degrees and $\mu$Hem,
respectively, and $A_{\rm BMR}$ is calculated from $A_{\rm BMR,S}$ using the
CCD relationship (Equation 11). Cameron et al. (2010) established this
empirical relationship with the $\Delta\beta$ and sunspot umbral area
measurements from the Kodaikanal Solar Observatory (KSO) and MWO. The tilt
angle of each BMR is calculated from the observed cycle-average tilt angle,
also from KSO and MWO, through the procedure detailed in Cameron et al.
(2010). The source term, and therefore the simulation runs, extend from 1913
to 1986 (solar cycle 15 to 21), the period over which sunspot group location,
area, and tilt angle measurements are available. The simulation runs and the
analysis based on them (see the following sections) are not sensitive to which
sunspot area composite we take into the CJSS model as input. The Hathaway &
Choudhary (2008), Balmaceda et al. (2009), and Mandal et al. (2020) composites
are essentially different efforts to extend the RGO record, which ends in
1976, to more recent times with sunspot area measurements from other
observatories. The various composites are therefore exactly identical to one
another up to this time (i.e. over most of the simulation period of 1913 to
1986).

Figure 3: The surface flux transport simulation runs. a) Average unsigned
photospheric magnetic flux density in the CJSS-SH (blue), CJSS-LCC (red), and
CJSS-YEA (black) simulation runs (Table 1). The CJSS-SH and CJSS-LCC time
series are largely hidden by the CJSS-YEA time series due to the close
similarity. The circles and squares denote the observed values from MWO and
WSO (Arge et al., 2002). b) Open flux corresponding to each run, derived by
extrapolating the model output photospheric magnetic field into the
heliosphere, and the reconstruction of this quantity from the $aa$ index by
Lockwood et al. (2009) (green solid line) and by Lockwood et al. (2014) (green
dashed line). In a) and b) the CJSS-LCC and CJSS-YEA time series are dashed
and dotted, respectively, to aid visibility. c) Mean magnetic field over the
north (solid lines) and the south polar caps (dashed lines) in each run. The
green crosses indicate the observed timing of the polar field reversals
reported by Makarov et al. (2004) (see Sects. 3.1 and 3.2 for details).

The CJSS model relies on the combination of the assumption that $B_{\pm}$
scales linearly with $\left(\Delta\beta\right)^{2}$ (Equation 15) and the
empirical relationship between $\Delta\beta$ and $A_{\rm BMR}$ (Equation 16)
to define $B_{\pm}$. Substituting Equation 16 into Equation 15 we obtain

$B_{\pm}\left(\lambda,\phi\right)=\tilde{B}_{\rm max}A_{\rm
BMR}\exp\left(-\frac{2\left[1-\cos\beta_{\pm}\left(\lambda,\phi\right)\right]}{\delta^{2}}\right),$
(17)

where $\tilde{B}_{\rm max}=B_{\rm max}0.4^{2}\left(6.46\times
10^{-2}\right)^{2}\delta^{-2}=0.0156\ {\rm G}\ \mu{}{\rm Hem^{-1}}$ is a
constant factor that serves the same purpose as $B_{\rm max}$. The consequence
of this approach is that $B_{\pm}$, and therefore also $\Phi_{\rm BMR}$,
scales linearly with $A_{\rm BMR}$ in the model. Incidentally, this is similar
to the constraint on $\Phi_{\rm BMR}$ reported by SH, where this quantity also
scales linearly with $A_{\rm BMR}$ (Equation 10).

Let us generalise Equation 17 by rewriting it as

$B_{\pm}\left(\lambda,\phi\right)=\tilde{B}_{\rm
max}{}f\exp\left(-\frac{2\left[1-\cos\beta_{\pm}\left(\lambda,\phi\right)\right]}{\delta^{2}}\right),$
(18)

where $f=A_{\rm BMR}$ determines the $A_{\rm BMR,S}$-dependence of $B_{\pm}$
and $\Phi_{\rm BMR}$. In the CJSS-LCC run we modify the CJSS model to comply
with the constraint on $\Phi_{\rm BMR}$ reported by LCC (Equation 13) by
redefining $f$ as $A_{\rm BMR,S}^{0.71}$. Similarly, in the CJSS-YEA run we
apply the constraint on $\Phi_{\rm BMR}$ derived in this study (Equation 8) to
the CJSS model by redefining $f$ as $\tilde{A}_{\rm BMR}$. By this
redefinition of $f$, we replace the assumption in the CJSS model that
$B_{\pm}$ scales linearly with $\left(\Delta\beta\right)^{2}$ with direct
observational constraints on the relationship between $\Phi_{\rm BMR}$ and
$A_{\rm BMR,S}$. As in the CJSS-SH run, $\tilde{B}_{\rm max}$ is set in the
CJSS-LCC and CJSS-YEA runs at the value that brings the modelled average
unsigned magnetic flux density to the absolute scale of the observed values
(Fig. 3a). The value of $\tilde{B}_{\rm max}$ and the definition of $f$ in
each simulation run is summarised in Table 1.

###  3.2 Simulation results

####  3.2.1 Total photospheric magnetic flux

In terms of the total photospheric magnetic flux, represented in Fig. 3a by
the average unsigned magnetic flux density, the three simulation runs (black,
blue, and red curves) are very similar. The agreement in absolute scale is
moot since $\tilde{B}_{\rm max}$ is set in each run at the level that brings
the model output to the absolute scale of the observed values (circles and
boxes). What is notable is that the time variability is also very similar even
with the different $A_{\rm BMR,S}$-dependence of $\Phi_{\rm BMR}$ imposed for
each run.

####  3.2.2 Solar open flux and polar field reversals

Cameron et al. (2010) and Jiang et al. (2011b) had verified their model by
demonstrating its ability to reproduce the evolution of the solar open flux
and the timing of the polar field reversals. The authors derived the open flux
corresponding to the model output photospheric magnetic field by extrapolating
it into the heliosphere using the Current Sheet Source Surface (CSSS) model by
Zhao & Hoeksema (1995b, a). This was matched to the independent reconstruction
of the open flux from the $aa$ geomagnetic activity index by Lockwood et al.
(2009). They compared the timing of the polar field reversals in the model
output to the timings reported by Makarov et al. (2004) based on an
examination of the polar filaments in H$\alpha$ synoptic charts.

Following CJSS, we applied the CSSS model to each simulation run to calculate
the corresponding solar open flux. In Fig. 3b we compare the results (black,
blue, and red) to the Lockwood et al. (2009) reconstruction (green solid line)
and the update by Lockwood et al. (2014) (green dashed line). The various runs
do not differ from one another more than they do from the Lockwood et al.
(2009, 2014) time series.

In Fig. 3c we chart the evolution of the polar field in the simulation runs as
given by the mean magnetic flux density over the north ($\lambda>75^{\circ}$,
solid lines) and the south polar caps ($\lambda<-75^{\circ}$, dashed lines).
Polar field reversal is indicated by the mean polar field switching polarity.
In the same plot we indicate the observed timing of the polar field reversals
reported by Makarov et al. (2004) (green crosses). Analogous to the solar open
flux, the various runs do not differ more from one another than they do from
the Makarov et al. (2004) observations. Applying the constraints on $\Phi_{\rm
BMR}$ presented by LCC (Equation 13) and derived in this study (Equation 8) to
the CJSS model retained its ability to reproduce the evolution of the solar
open flux and the timing of the polar field reversals.

####  3.2.3 Correlation between the polar field and solar cycle strength

A key finding of the CJSS studies is that the amount of polar flux at the end
of each solar cycle in the model appears to be highly correlated to the
observed strength of the succeeding cycle. This is consistent with the
existence of a Babcock–Leighton-type solar dynamo (Babcock, 1961; Leighton,
1964, 1969), where the polar field, which peaks around the end of each cycle,
is the main source of the toroidal field in the convection zone (see
discussion in Cameron & Schüssler, 2015). We examine if this feature of the
CJSS model is sensitive to the $A_{\rm BMR,S}$-dependence of $\Phi_{\rm BMR}$
assumed. Taking each simulation run, we compare the peak of the mean unsigned
magnetic flux density over both polar caps at the end of each cycle to the
maximum total sunspot area over the following cycle, representing modelled
polar flux and measured cycle strength (Fig. 4). The Pearson’s correlation
coefficient, $R$, between the two quantities is 0.933 for the CJSS-SH run,
0.925 for the CJSS-LCC run, and 0.937 for the CJSS-YEA run. The different
$A_{\rm BMR,S}$-dependence of $\Phi_{\rm BMR}$ imposed on the various runs
made no appreciable difference to the correlation between the model output
polar field at the end of each cycle and the observed strength of the
succeeding cycle.

Figure 4: For each simulation run, the peak of the mean unsigned magnetic flux
density over the polar caps at the end of solar cycle $n$ vs the maximum total
surface coverage by sunspots over cycle $n+1$, where $n$ ranges from 15 to 21.
We excluded solar cycle 14 as the end of this cycle occurred at the start of
the simulation, and is therefore sensitive to the arbitrary initial condition.
The crosses and lines denote the individual cycle values and the corresponding
linear fit. The Pearson’s correlation coefficient, $R$, is indicated (see
Sect. 3.2.3 for details).

Whether in terms of the time variation of the total photospheric magnetic flux
(Fig. 3a) and the solar open flux (Fig. 3b), the timing of the polar field
reversals (Fig. 3c), or the correlation between the polar field at the end of
each cycle and the observed strength of the succeeding cycle (Fig. 4), the
various simulation runs are very similar. It appears that the CJSS model is,
at least to the extent tested here, not particularly sensitive to the $A_{\rm
BMR,S}$-dependence of $\Phi_{\rm BMR}$ assumed. This is likely, at least in
part, due to the fact that the magnetic flux in each newly emerged BMR can
persist for weeks or months, such that the state of a simulation run at a
given time depends not only on the BMRs injected then but also on all the BMRs
that were introduced in the weeks to months prior. The amount of magnetic flux
in individual BMRs depends strongly on the $A_{\rm BMR,S}$-dependence of
$\Phi_{\rm BMR}$ assumed. However, the aggregate magnetic flux introduced over
an extended period by all the BMRs that had emerged then would not be as
sensitive to the $A_{\rm BMR,S}$-dependence of $\Phi_{\rm BMR}$ assumed due to
the summation. Nonetheless, it is noteworthy that the application of the
empirical relationship between $\Phi_{\rm BMR}$ and $A_{\rm BMR,S}$ presented
here to what is an established SFTM retained its key features, including its
ability to replicate various independent datasets.

##  4 Balance between bright faculae and network and dark sunspots at various
solar activity levels

As discussed in the Introduction, the balance between bright faculae and
network and dark sunspots on the solar disc and how that changes with solar
activity level is of interest for the relevance to the understanding of the
brightness variability of the Sun and other cool stars. Making use of the
CJSS-YEA run (Sect. 3), let us examine how the amount of magnetic flux on the
solar disc in faculae and network, $\Sigma\Phi_{\rm FN}$, and in sunspots,
$\Sigma\Phi_{\rm S}$, relate to one another and how the ratio
$\frac{\Sigma\Phi_{\rm FN}}{\Sigma\Phi_{\rm S}}$ changes with solar activity
level. The quantity $\Sigma\Phi_{\rm S}$ is calculated from the Balmaceda et
al. (2009) $\Sigma{}A_{\rm S}$ time series (Fig. 1a) using Equation 7. We
obtained the total magnetic flux on the solar disc from the CJSS-YEA run and
subtracted $\Sigma\Phi_{\rm S}$ from it to yield $\Sigma\Phi_{\rm FN}$. In the
following discussion we examine the annual mean of each quantity.

In Fig. 5a we chart $\Sigma\Phi_{\rm FN}$ estimated from the CJSS-YEA run as a
function of $\Sigma\Phi_{\rm S}$ (black crosses). As a check, we compare this
to the YSK model reconstruction of $\Sigma\Phi_{\rm FN}$ from $\Sigma\Phi_{\rm
S}$ (red crosses). Before we discuss the trend in $\Sigma\Phi_{\rm FN}$ with
$\Sigma\Phi_{\rm S}$, let us address the following. Firstly, both the CJSS and
YSK models describe the evolution of the photospheric magnetic field due to
sunspot-bearing BMRs alone, ignoring ERs. Secondly, the two estimates of
$\Sigma\Phi_{\rm FN}$ differ in absolute scale and we have to rescale the
values from the YSK model by a factor of 1.46 to bring them to the absolute
scale of the values from the CJSS-YEA run. Thirdly, plotted against
$\Sigma\Phi_{\rm S}$, the values from the CJSS-YEA run are markedly more
scattered than the values from the YSK model.

Figure 5: The disc-integrated facular and network magnetic flux,
$\Sigma\Phi_{\rm FN}$, as indicated by the CJSS-YEA run (black crosses) and by
the YSK model (red crosses), as a function of the disc-integrated sunspot
magnetic flux, $\Sigma\Phi_{\rm S}$. The first solar cycle in the CJSS-YEA run
(cycle 15, 1913 to 1923) has been excluded as the simulation is sensitive to
the arbitrary initial condition here. The annual mean over the period of 1924
to 1986, where there is data from both models, is shown. The black line and
the grey region represents the quadratic polynomial fit to the CJSS-YEA run
values and the corresponding 95$\%$ confidence interval, while the red line
represents the fit to the YSK model values. See Sect. 4 for details.

Even though the CJSS model neglects ERs, considering the ability of the CJSS-
YEA run to reproduce various independent datasets (Sect. 3.2), the effect of
this limitation on the corresponding estimate of $\Sigma\Phi_{\rm FN}$ is
likely minor. A similar argument applies to the YSK model estimate of
$\Sigma\Phi_{\rm FN}$. The absolute scale of the CJSS model is set to that of
the average unsigned magnetic flux density values reported by Arge et al.
(2002), which is based on MWO and WSO magnetograms (circles and squares, Fig.
4a). The absolute scale of the YSK model is set to that of the
$\Sigma\Phi_{\rm FN}$ time series reported by Yeo et al. (2014b), which is
based on Kitt Peak, MDI, and HMI magnetograms (black, Fig. 1b). We attribute
the different absolute scale of the CJSS and YSK model outputs to systematic
uncertainties in how Arge et al. (2002) and Yeo et al. (2014b) derived their
respective datasets and in the absolute magnetogram signal scale of the
various solar telescopes (Riley et al., 2014). In the YSK model, the variation
in $\Sigma\Phi_{\rm FN}$ due to sunspot-bearing BMRs is given by the
convolution of $\Sigma{}A_{\rm S}^{h_{1}}$ with a FIR filter, $H$ (Equations 1
and 2). The implicit assumption is that all BMRs exhibit the same time
evolution as delineated by $H$. This ignores the fact that even BMRs
containing the same amount of magnetic flux can evolve rather differently due
to factors such as latitude, tilt angle, and the magnetic field in their
surroundings. These factors are taken into account in the CJSS model,
producing the greater scatter in the corresponding $\Sigma\Phi_{\rm FN}$
versus $\Sigma\Phi_{\rm S}$ plot.

To more clearly show the trend in $\Sigma\Phi_{\rm FN}$ with $\Sigma\Phi_{\rm
S}$, we fit a quadratic polynomial in $\Sigma\Phi_{\rm S}$ to each of the two
estimates of $\Sigma\Phi_{\rm FN}$. The fit to the CJSS-YEA run values (black
line, Fig. 5) is given by

$\displaystyle\Sigma\Phi_{\rm FN}=$ $\displaystyle\left(0.52\pm
0.11\right)+\left(3.33\pm 0.53\right)\Sigma\Phi_{\rm S}$ (19)
$\displaystyle-\left(0.73\pm 0.46\right)\Sigma\Phi_{\rm S}^{2}$

and the fit to the YSK model values (red line) by

$\displaystyle\Sigma\Phi_{\rm FN}=$ $\displaystyle\left(0.334\pm
0.019\right)+\left(3.367\pm 0.086\right)\Sigma\Phi_{\rm S}$ (20)
$\displaystyle-\left(0.425\pm 0.074\right)\Sigma\Phi_{\rm S}^{2}.$

This revealed that $\Sigma\Phi_{\rm FN}$ increases with increasing
$\Sigma\Phi_{\rm S}$ at a diminishing rate such that it appears to gradually
saturate, analogous to what earlier studies noted comparing facular area and
various disc-integrated quantities that are sensitive to the amount of faculae
and network present to sunspot area and number (Foukal, 1993, 1996, 1998;
Chapman et al., 1997; Solanki & Fligge, 1998; Shapiro et al., 2014; Yeo et
al., 2020), discussed in the Introduction. The fit to the CJSS-YEA run values
(black line) and the fit to the YSK model values (red line) do differ
slightly, though the latter lies within the 95$\%$ confidence interval of the
former (shaded region). This means that, taking the scatter in the CJSS-YEA
run values (black crosses) into account, that there is little separating the
two estimates of $\Sigma\Phi_{\rm FN}$ in terms of the $\Sigma\Phi_{\rm
S}$-dependence. The CJSS-YEA run and the YSK model are not entirely
independent of one another as we had constrained the former using the
empirical relationship between $\Phi_{\rm BMR}$ and $A_{\rm BMR,S}$ derived
from the latter. Nonetheless, it is significant that two models that adopt
vastly different approaches to modelling the evolution of the photospheric
magnetic field returned similar results.

Next, we examine what the apparent trend in $\Sigma\Phi_{\rm FN}$ with
$\Sigma\Phi_{\rm S}$ implies for the balance between the amount of faculae and
network and the amount of sunspots on the solar disc and how that changes with
solar activity level. Sunspot indices are often employed as proxies of solar
activity, while stellar activity is often characterised by the quantity
$\log{}R^{\prime}_{\rm HK}$ defined in Noyes et al. (1984). Taking this into
consideration, let us examine $\frac{\Sigma\Phi_{\rm FN}}{\Sigma\Phi_{\rm
S}}$, as estimated from the CJSS-YEA run, as a function of $\Sigma{}A_{\rm S}$
and of $\log{}R^{\prime}_{\rm HK}$ (crosses, Fig. 6). We make use of the model
reconstruction of $\log{}R^{\prime}_{\rm HK}$ from $\Sigma{}A_{\rm S}$
presented by YSK.

Figure 6: The variation in the balance between bright faculae and network and
dark sunspots with solar activity. a) The ratio $\frac{\Sigma\Phi_{\rm
FN}}{\Sigma\Phi_{\rm S}}$, where $\Sigma\Phi_{\rm FN}$ is obtained from the
CJSS-YEA run, as a function of $\Sigma{}A_{\rm S}$ (crosses). As in Fig. 5,
the annual mean over the period of 1924 to 1986 is depicted. The solid line
follows the power-law fit (Equation 22) and the dotted line gives the position
of the horizontal asymptote as given by the constant term in this fit. b) The
same, but for $\frac{\Sigma\Phi_{\rm FN}}{\Sigma\Phi_{\rm S}}$ as a function
of $\log{}R^{\prime}_{\rm HK}$. The vertical dashed line indicates
$\log{}R^{\prime}_{\rm HK}=-4.9$ and the horizontal dashed lines gives the
value of $\frac{\Sigma\Phi_{\rm FN}}{\Sigma\Phi_{\rm S}}$ at
$\log{}R^{\prime}_{\rm HK}=-4.9$ and $\log{}R^{\prime}_{\rm HK}=-4.7$, as
given by the power-law fit (Equation 23) (see Sect. 4 for details).

The gradual saturation in $\Sigma\Phi_{\rm FN}$ as we go from low to high
$\Sigma\Phi_{\rm S}$ (Fig. 5) translates into $\frac{\Sigma\Phi_{\rm
FN}}{\Sigma\Phi_{\rm S}}$ decreasing with increasing solar activity level at a
diminishing rate such that it appears to approach a horizontal asymptote (Fig.
6). Let us model $\frac{\Sigma\Phi_{\rm FN}}{\Sigma\Phi_{\rm S}}$ as
conforming to a power-law relationship with activity level such that

$\frac{\Sigma\Phi_{\rm FN}}{\Sigma\Phi_{\rm
S}}=g_{1}\left[\alpha+g_{2}\right]^{g_{3}}+g_{4},$ (21)

where $\alpha$ represents activity as given by either $\Sigma{}A_{\rm S}$ or
$\log{}R^{\prime}_{\rm HK}$ and $g_{1}$ to $g_{4}$ are fit parameters. The
$g_{1}\left[\alpha+g_{2}\right]^{g_{3}}$ term describes the variation in
$\frac{\Sigma\Phi_{\rm FN}}{\Sigma\Phi_{\rm S}}$ with $\alpha$, while $g_{4}$
denotes the position of the horizontal asymptote. We fit Equation 21 to the
$\frac{\Sigma\Phi_{\rm FN}}{\Sigma\Phi_{\rm S}}$ versus $\Sigma{}A_{\rm S}$
plot (Fig. 6a) and the $\frac{\Sigma\Phi_{\rm FN}}{\Sigma\Phi_{\rm S}}$ versus
$\log{}R^{\prime}_{\rm HK}$ plot (Fig. 6b) in turn. For the former, $g_{2}$ is
fixed at null. At convergence we arrive at

$\frac{\Sigma\Phi_{\rm FN}}{\Sigma\Phi_{\rm S}}=\left(690\pm
110\right)\Sigma{}A_{\rm S}^{\left(-0.876\pm 0.043\right)}+\left(2.63\pm
0.29\right)$ (22)

and

$\displaystyle\frac{\Sigma\Phi_{\rm FN}}{\Sigma\Phi_{\rm S}}$ $\displaystyle=$
(23) $\displaystyle\left(0.038\pm 0.060\right)\left[\log{}R^{\prime}_{\rm
HK}+\left(4.9405\pm 0.0046\right)\right]^{\left(-1.35\pm 0.45\right)}$
$\displaystyle+\left(2.39\pm 0.72\right),$

represented in Fig. 6 by the solid lines. In both instances the $g_{4}$ term
or the position of the horizontal asymptote (dotted lines) is above unity,
implying that the amount of magnetic flux on the solar disc in sunspots never
exceeds that in faculae and network. This is consistent with Parnell et al.
(2009) and Thornton & Parnell (2011), who came to a similar conclusion by
examining the distribution of magnetic features with the amount of flux in
each feature.

The observation here that most of the magnetic flux on the solar disc is in
faculae and network does not contradict our earlier assertion that most of the
magnetic flux in newly emerged BMRs is in the enclosed sunspots (Sect. 2.2).
The magnetic flux in BMRs persists long after the comparatively short-lived
sunspots they enclose have dissipated, and manifests as faculae and network
(van Driel-Gesztelyi & Green, 2015). In combination with the fact that the
photospheric magnetic field at a given time is comprised of not just the BMRs
that had emerged then, but also what remains from all the BMRs that had
emerged in the weeks to months prior, the two assertions we are making do not
exclude one another.

In the context of the brightness variability of cool stars, the transition
between the faculae-dominated and spot-dominated regimes is estimated to be in
the range of $\log{}R^{\prime}_{\rm HK}$ of $-4.9$ to $-4.7$ (Lockwood et al.,
2007; Hall et al., 2009; Shapiro et al., 2014; Radick et al., 2018; Reinhold
et al., 2019). Looking at Fig. 6b, the $\log{}R^{\prime}_{\rm HK}$ of the Sun
(crosses) is in the range of $-4.93$ to $-4.84$, around the lower end of the
range of activity where this transition occurs. The extrapolation of the
power-law fit (solid line) indicates that in the range of
$\log{}R^{\prime}_{\rm HK}$ of $-4.9$ to $-4.7$, $\frac{\Sigma\Phi_{\rm
FN}}{\Sigma\Phi_{\rm S}}$ varies slowly with activity, such that it stays
within a relatively narrow range of about 3 to 5 (dashed lines). Taken
together, these two observations imply that the $\frac{\Sigma\Phi_{\rm
FN}}{\Sigma\Phi_{\rm S}}$ level of the Sun might not be far above the level
where facular and network brightening and sunspot darkening would cancel one
another out. This concurs with the findings of Shapiro et al. (2016) and
Radick et al. (2018), who noted that while the Sun is faculae-dominated, it is
only marginally so as facular and network brightening and sunspot darkening
appear to be closely balanced.

##  5 Summary

The relationship between sunspot-bearing bipolar magnetic regions (BMRs) and
the sunspots they enclose, connected to how solar magnetism and activity
relate to sunspot prevalence, is important. It is particularly relevant to
surface flux transport models (SFTMs) that derive the details of BMR emergence
they require as input from sunspot observations. However, attempts to
constrain how the amount of magnetic flux in each newly emerged BMR,
$\Phi_{\rm BMR}$, relates to its area (Schrijver & Harvey, 1994) or the area
of the enclosed sunspots, $A_{\rm BMR,S}$ (Lemerle et al., 2015), were limited
by the fact that there is no extensive and reliable database of the magnetic
and physical properties of newly emerged BMRs currently available. In the
absence of proper observational constraints, existing SFTMs adopt various
simplifications to determine the details of BMR emergence from sunspot
observations. In this study we aimed to establish the relationship between
$\Phi_{\rm BMR}$ and $A_{\rm BMR,S}$.

To this end, we made use of the empirical model of the relationship between
the disc-integrated facular and network magnetic flux, $\Sigma\Phi_{\rm FN}$,
and the total surface coverage by sunspots, $\Sigma{}A_{\rm S}$, presented by
Yeo et al. (2020). The structure of the Yeo et al. (2020) model enabled us to
circumvent the lack of any proper BMR database, and instead derive an
empirical relationship between $\Phi_{\rm BMR}$ and $A_{\rm BMR,S}$ (Equation
8) from $\Sigma{}A_{\rm FN}$ and $\Sigma{}A_{\rm S}$ (i.e. disc-integrated
quantities). The empirical relationship between $\Phi_{\rm BMR}$ and $A_{\rm
BMR,S}$ indicates that for newly emerged BMRs, the enclosed sunspots occupy a
small proportion of their area (14$\%$ to 18$\%$), but contain most of their
magnetic flux (70$\%$ to 90$\%$) (Fig. 2b). We found evidence that the earlier
results by Schrijver & Harvey (1994) (Equation 10) and Lemerle et al. (2015)
(Equation 13) might underestimate the $\Phi_{\rm BMR}$ of larger BMRs (Fig.
2a).

To verify the empirical relationship between $\Phi_{\rm BMR}$ and $A_{\rm
BMR,S}$ presented here, we examined the result of applying this and the
earlier results by Schrijver & Harvey (1994) and Lemerle et al. (2015) to a
SFTM. We made use of the SFTM reported by Cameron et al. (2010) and Jiang et
al. (2011a, b). This model, incidentally and not by design, already complies
with the Schrijver & Harvey (1994) result. We generated three simulation runs,
one with the model as it is and the other two with the model modified to
comply with the Lemerle et al. (2015) results and the current results. We
compared the evolution of the solar open flux to the independent
reconstruction of this quantity from the $aa$ geomagnetic activity index by
Lockwood et al. (2009, 2014) (Fig. 3b) and the timing of the polar field
reversals to the observed timing reported by Makarov et al. (2004) (Fig. 3c).
The various runs reproduce the Lockwood et al. (2009, 2014) and Makarov et al.
(2004) results relatively well and almost equally. An important feature of the
Cameron et al. (2010) and Jiang et al. (2011a, b) model is that the polar
field at the end of each cycle is highly correlated to the observed strength
of the following cycle. This is significant as it is consistent with the
existence of a Babcock–Leighton-type solar dynamo. In terms of this
correlation, the various runs are again very similar (Fig. 4). Surface flux
transport simulations do not appear to be, to the extent tested here,
particularly sensitive to the $A_{\rm BMR,S}$-dependence of $\Phi_{\rm BMR}$
assumed. Nonetheless, it is worth noting that applying the current result to
what is an established SFTM retained its key features.

Making use of the surface flux transport simulation based on the relationship
between $\Phi_{\rm BMR}$ and $A_{\rm BMR,S}$ presented here, we examined how
the amount of magnetic flux on the solar disc in faculae and network,
$\Sigma\Phi_{\rm FN}$, and in sunspots, $\Sigma\Phi_{\rm S}$, relate to one
another and the relationship between the ratio $\frac{\Sigma\Phi_{\rm
FN}}{\Sigma\Phi_{\rm S}}$ and the solar activity level. We found
$\Sigma\Phi_{\rm FN}$ to increase with increasing $\Sigma\Phi_{\rm S}$ at a
diminishing rate such that it appears to gradually saturate (Fig. 5). This
complies with the analogous results from the earlier studies that compared
facular area and various disc-integrated quantities sensitive to the amount of
faculae and network present to sunspot area and number (Foukal, 1993, 1996,
1998; Chapman et al., 1997; Solanki & Fligge, 1998; Shapiro et al., 2014; Yeo
et al., 2020). Taking $\Sigma{}A_{\rm S}$ and $\log{}R^{\prime}_{\rm HK}$ as
proxies of solar activity, we found $\frac{\Sigma\Phi_{\rm
FN}}{\Sigma\Phi_{\rm S}}$ to decrease with increasing activity at a
diminishing rate such that it appears to approach a horizontal asymptote (Fig.
6). We modelled $\frac{\Sigma\Phi_{\rm FN}}{\Sigma\Phi_{\rm S}}$ as conforming
to a power-law relationship with activity (Equations 22 and 23). The power-law
model indicates that, at least within the range of activity examined, the
proportion of the magnetic flux on the solar disc that is in sunspots remains
smaller than that in faculae and network. This does not contradict our claim
that most of the magnetic flux in newly emerged BMRs is in the enclosed
sunspots (see Sect. 4). The extrapolation of the power-law model to higher
activity levels indicates that the $\frac{\Sigma\Phi_{\rm FN}}{\Sigma\Phi_{\rm
S}}$ level of the Sun might be close to the level where facular and network
brightening and sunspot darkening would cancel one another out, consistent
with Shapiro et al. (2016) and Radick et al. (2018).

###### Acknowledgements.

The authors are grateful to R. H. Cameron of the Max-Planck Institut für
Sonnensystemforschung for the useful discussions and for granting access to
the surface flux transport model featured in this work. KLY and NAK received
funding from the German Federal Ministry of Education and Research (project
01LG1909C), KLY and SKS from the European Research Council through the Horizon
2020 research and innovation program of the European Union (grant 695075) and
JJ from the National Science Foundation of China (grant 11873023 and 11522325)
and the Fundamental Research Funds for the Central Universities of China.

## References

  * Arge et al. (2002) Arge, C. N., Hildner, E., Pizzo, V. J., & Harvey, J. W. 2002, J. Geophys. Res., 107, 1319
  * Babcock (1961) Babcock, H. W. 1961, ApJ, 133, 572
  * Balmaceda et al. (2009) Balmaceda, L. A., Solanki, S. K., Krivova, N. A., & Foster, S. 2009, J. Geophys. Res., 114, 7104
  * Baumann et al. (2006) Baumann, I., Schmitt, D., & Schüssler, M. 2006, A&A, 446, 307
  * Baumann et al. (2004) Baumann, I., Schmitt, D., Schüssler, M., & Solanki, S. K. 2004, A&A, 426, 1075
  * Cameron & Schüssler (2015) Cameron, R. & Schüssler, M. 2015, Science, 347, 1333
  * Cameron et al. (2010) Cameron, R. H., Jiang, J., Schmitt, D., & Schüssler, M. 2010, ApJ, 719, 264
  * Cameron et al. (2012) Cameron, R. H., Schmitt, D., Jiang, J., & Isık, E. 2012, A&A, 542, A127
  * Chapman et al. (1997) Chapman, G. A., Cookson, A. M., & Dobias, J. J. 1997, ApJ, 482, 541
  * Charbonneau (2020) Charbonneau, P. 2020, Living Rev. Sol. Phys., 17, 4
  * Cho et al. (2015) Cho, I. H., Cho, K. S., Bong, S. C., et al. 2015, ApJ, 811, 49
  * Coddington et al. (2016) Coddington, O., Lean, J. L., Pilewskie, P., Snow, M., & Lindholm, D. 2016, Bull. Am. Meteor. Soc., 97, 1265
  * Dasi-Espuig et al. (2014) Dasi-Espuig, M., Jiang, J., Krivova, N. A., & Solanki, S. K. 2014, A&A, 570, A23
  * Dasi-Espuig et al. (2016) Dasi-Espuig, M., Jiang, J., Krivova, N. A., et al. 2016, A&A, 590, A63
  * DeVore et al. (1984) DeVore, C. R., Boris, J. P., & Sheeley, N. R., J. 1984, Sol. Phys., 92, 1
  * Domingo et al. (2009) Domingo, V., Ermolli, I., Fox, P., et al. 2009, Space Sci. Rev., 145, 337
  * Foukal (1993) Foukal, P. 1993, Sol. Phys., 148, 219
  * Foukal (1996) Foukal, P. 1996, Geochim. Res. Lett., 23, 2169
  * Foukal (1998) Foukal, P. 1998, ApJ, 500, 958
  * Hagenaar (2001) Hagenaar, H. J. 2001, ApJ, 555, 448
  * Hagenaar et al. (2003) Hagenaar, H. J., Schrijver, C. J., & Title, A. M. 2003, ApJ, 584, 1107
  * Haigh (2007) Haigh, J. D. 2007, Living Rev. Sol. Phys., 4, 2
  * Hall et al. (2009) Hall, J. C., Henry, G. W., Lockwood, G. W., Skiff, B. A., & Saar, S. H. 2009, AJ, 138, 312
  * Harvey (1993) Harvey, K. L. 1993, PhD thesis, Utrecht University
  * Harvey (2000) Harvey, K. L. 2000, in Encyclopedia of Astronomy and Astrophysics, ed. P. Murdin (Nature Publishing Group)
  * Harvey & Zwaan (1993) Harvey, K. L. & Zwaan, C. 1993, Sol. Phys., 148, 85
  * Hathaway & Choudhary (2008) Hathaway, D. H. & Choudhary, D. P. 2008, Sol. Phys., 250, 269
  * Hoyt & Schatten (1998) Hoyt, D. V. & Schatten, K. H. 1998, Sol. Phys., 181, 491
  * Jiang et al. (2013) Jiang, J., Cameron, R. H., Schmitt, D., & Isık, E. 2013, A&A, 553, A128
  * Jiang et al. (2011a) Jiang, J., Cameron, R. H., Schmitt, D., & Schüssler, M. 2011a, A&A, 528, A82
  * Jiang et al. (2011b) Jiang, J., Cameron, R. H., Schmitt, D., & Schüssler, M. 2011b, A&A, 528, A83
  * Jiang et al. (2014) Jiang, J., Hathaway, D. H., Cameron, R. H., et al. 2014, Space Sci. Rev., 186, 491
  * Jones et al. (1992) Jones, H. P., Duvall, Thomas L., J., Harvey, J. W., et al. 1992, Sol. Phys., 139, 211
  * Keppens & Martinez Pillet (1996) Keppens, R. & Martinez Pillet, V. 1996, A&A, 316, 229
  * Krivova et al. (2007) Krivova, N. A., Balmaceda, L., & Solanki, S. K. 2007, A&A, 467, 335
  * Krivova et al. (2010) Krivova, N. A., Vieira, L. E. A., & Solanki, S. K. 2010, J. Geophys. Res., 115, A12112
  * Lean (2000) Lean, J. 2000, Geochim. Res. Lett., 27, 2425
  * Leighton (1964) Leighton, R. B. 1964, ApJ, 140, 1547
  * Leighton (1969) Leighton, R. B. 1969, ApJ, 156, 1
  * Lemerle et al. (2015) Lemerle, A., Charbonneau, P., & Carignan-Dugas, A. 2015, ApJ, 810, 78
  * Lites et al. (1993) Lites, B. W., Elmore, D. F., Seagraves, P., & Skumanich, A. P. 1993, ApJ, 418, 928
  * Livingston et al. (1976) Livingston, W. C., Harvey, J., Pierce, A. K., et al. 1976, Appl. Opt., 15, 33
  * Lockwood et al. (2007) Lockwood, G. W., Skiff, B. A., Henry, G. W., et al. 2007, ApJS, 171, 260
  * Lockwood et al. (2014) Lockwood, M., Nevanlinna, H., Barnard, L., et al. 2014, Ann. Geophys., 32, 383
  * Lockwood et al. (2009) Lockwood, M., Owens, M., & Rouillard, A. P. 2009, J. Geophys. Res., 114, A11104
  * Mackay et al. (2000) Mackay, D. H., Gaizauskas, V., & van Ballegooijen, A. A. 2000, ApJ, 544, 1122
  * Mackay & Yeates (2012) Mackay, D. H. & Yeates, A. R. 2012, Living Rev. Sol. Phys., 9, 6
  * Makarov et al. (2004) Makarov, V. I., Tlatov, A. G., & Callebaut, D. K. 2004, Sol. Phys., 224, 49
  * Mandal et al. (2020) Mandal, S., Krivova, N. A., Solanki, S. K., Sinha, N., & Banerjee, D. 2020, A&A, 640, A78
  * Nèmec et al. (2020a) Nèmec, N. E., Isık, E., Shapiro, A. I., et al. 2020a, A&A, 638, A56
  * Nèmec et al. (2020b) Nèmec, N. E., Shapiro, A. I., Krivova, N. A., et al. 2020b, A&A, 636, A43
  * Noyes et al. (1984) Noyes, R. W., Hartmann, L. W., Baliunas, S. L., Duncan, D. K., & Vaughan, A. H. 1984, ApJ, 279, 763
  * Parnell et al. (2009) Parnell, C. E., DeForest, C. E., Hagenaar, H. J., et al. 2009, ApJ, 698, 75
  * Preminger & Walton (2005) Preminger, D. G. & Walton, S. R. 2005, Geochim. Res. Lett., 32, L14109
  * Preminger & Walton (2006a) Preminger, D. G. & Walton, S. R. 2006a, Geochim. Res. Lett., 33, L23108
  * Preminger & Walton (2006b) Preminger, D. G. & Walton, S. R. 2006b, Sol. Phys., 235, 387
  * Preminger & Walton (2007) Preminger, D. G. & Walton, S. R. 2007, Sol. Phys., 240, 17
  * Radick et al. (2018) Radick, R. R., Lockwood, G. W., Henry, G. W., Hall, J. C., & Pevtsov, A. A. 2018, ApJ, 855, 75
  * Reinhold et al. (2019) Reinhold, T., Bell, K. J., Kuszlewicz, J., Hekker, S., & Shapiro, A. I. 2019, A&A, 621, A21
  * Riley et al. (2014) Riley, P., Ben-Nun, M., Linker, J. A., et al. 2014, Sol. Phys., 289, 769
  * Scherrer et al. (1995) Scherrer, P. H., Bogart, R. S., Bush, R. I., et al. 1995, Sol. Phys., 162, 129
  * Scherrer et al. (2012) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207
  * Schrijver (1987) Schrijver, C. J. 1987, A&A, 180, 241
  * Schrijver et al. (2002) Schrijver, C. J., De Rosa, M. L., & Title, A. M. 2002, ApJ, 577, 1006
  * Schrijver & Harvey (1994) Schrijver, C. J. & Harvey, K. L. 1994, Sol. Phys., 150, 1
  * Shapiro et al. (2014) Shapiro, A. I., Solanki, S. K., Krivova, N. A., et al. 2014, A&A, 569, A38
  * Shapiro et al. (2015) Shapiro, A. I., Solanki, S. K., Krivova, N. A., Tagirov, R. V., & Schmutz, W. K. 2015, A&A, 581, A116
  * Shapiro et al. (2016) Shapiro, A. I., Solanki, S. K., Krivova, N. A., Yeo, K. L., & Schmutz, W. K. 2016, A&A, 589, A46
  * Sheeley (1966) Sheeley, N. R., J. 1966, ApJ, 144, 723
  * Sofia et al. (1982) Sofia, S., Oster, L., & Schatten, K. 1982, Sol. Phys., 80, 87
  * Solanki & Fligge (1998) Solanki, S. K. & Fligge, M. 1998, Geochim. Res. Lett., 25, 341
  * Solanki et al. (2013) Solanki, S. K., Krivova, N. A., & Haigh, J. D. 2013, ARA&A, 51, 311
  * Solanki & Schmidt (1993) Solanki, S. K. & Schmidt, H. U. 1993, A&A, 267, 287
  * Solanki et al. (2000) Solanki, S. K., Schüssler, M., & Fligge, M. 2000, Nature, 408, 445
  * Solanki et al. (2002) Solanki, S. K., Schüssler, M., & Fligge, M. 2002, A&A, 383, 706
  * Thornton & Parnell (2011) Thornton, L. M. & Parnell, C. E. 2011, Sol. Phys., 269, 13
  * van Ballegooijen et al. (1998) van Ballegooijen, A. A., Cartledge, N. P., & Priest, E. R. 1998, ApJ, 501, 866
  * van Driel-Gesztelyi & Green (2015) van Driel-Gesztelyi, L. & Green, L. M. 2015, Living Rev. Sol. Phys., 12, 1
  * Virtanen et al. (2019) Virtanen, I. O. I., Virtanen, I. I., Pevtsov, A. A., et al. 2019, A&A, 627, A11
  * Virtanen et al. (2018) Virtanen, I. O. I., Virtanen, I. I., Pevtsov, A. A., & Mursula, K. 2018, A&A, 616, A134
  * Švanda et al. (2016) Švanda, M., Brun, A. S., Roudier, T., & Jouve, L. 2016, A&A, 586, A123
  * Wang et al. (2005) Wang, Y.-M., Lean, J. L., & Sheeley, Jr., N. R. 2005, ApJ, 625, 522
  * Wang et al. (1989) Wang, Y. M., Nash, A. G., & Sheeley, N. R., J. 1989, Science, 245, 712
  * Wang & Sheeley (1989) Wang, Y. M. & Sheeley, N. R., J. 1989, Sol. Phys., 124, 81
  * Whitbread et al. (2018) Whitbread, T., Yeates, A. R., & Muñoz-Jaramillo, A. 2018, ApJ, 863, 116
  * Wu et al. (2018) Wu, C. J., Krivova, N. A., Solanki, S. K., & Usoskin, I. G. 2018, A&A, 620, A120
  * Yeates (2020) Yeates, A. R. 2020, Sol. Phys., 295, 119
  * Yeates et al. (2007) Yeates, A. R., Mackay, D. H., & van Ballegooijen, A. A. 2007, Sol. Phys., 245, 87
  * Yeo et al. (2014a) Yeo, K. L., Krivova, N. A., & Solanki, S. K. 2014a, Space Sci. Rev., 186, 137
  * Yeo et al. (2014b) Yeo, K. L., Krivova, N. A., Solanki, S. K., & Glassmeier, K. H. 2014b, A&A, 570, A85
  * Yeo et al. (2020) Yeo, K. L., Solanki, S. K., & Krivova, N. A. 2020, A&A, 639
  * Yeo et al. (2017) Yeo, K. L., Solanki, S. K., Norris, C. M., et al. 2017, Phys. Rev. Lett., 119
  * Zhao & Hoeksema (1995a) Zhao, X. & Hoeksema, J. T. 1995a, Space Sci. Rev., 72, 189
  * Zhao & Hoeksema (1995b) Zhao, X. & Hoeksema, J. T. 1995b, J. Geophys. Res., 100, 19