Modelling solar irradiance from ground-based photometric observations

Total solar irradiance (TSI) has been monitored from space since 1978. The measurements show a prominent variability in phase with the solar cycle, as well as fluctuations on timescales shorter than a few days. However, the measurements were done by multiple and usually relatively short-lived missions making the possible long-term trend in the TSI highly uncertain. While the variability in the UV irradiance is clearly in-phase with the solar cycle, the phase of the variability in the visible range has been debated. In this paper, we aim at getting an insight into the long-term trend of TSI since 1996 and the phase of the solar irradiance variations in the visible part of the spectrum. We use independent ground-based full-disc photometric observations in Ca~II~K and continuum from the Rome and San Fernando observatories to compute the TSI since 1996. We follow the empirical San Fernando approach based on the photometric sum index. We find a weak declining trend in the TSI of -7.8$^{+4.9}_{-0.8}\times10^{-3}$ Wm$^{-2}$y$^{-1}$ between the 1996 and 2008 activity minima, while between 2008 and 2019 the reconstructed TSI shows no trend to a marginally decreasing (but statistically insignificant) trend of -0.1$^{+0.25}_{-0.02}\times10^{-3}$ Wm$^{-2}$y$^{-1}$. The reference TSI series used for the reconstruction does not significantly affect the determined trend. The variation in the blue continuum (409.2 nm) is rather flat, while the variation in the red continuum (607.1 nm) is marginally in anti-phase, although this result is extremely sensitive to the accurate assessment of the quiet Sun level in the images. These results provide further insights into the long-term variation of the total solar irradiance. The amplitude of the variations in the visible is below the uncertainties of the processing, which prevents an assessment of the phase of the variations.


Introduction
The solar irradiance (SI) is an important parameter in studies of the solar variability and Earth's climate. The spectrally-resolved SI, called the spectral solar irradiance (SSI), is defined as the solar radiative energy flux per unit area and wavelength as measured at the top of the Earth's atmosphere at a mean distance of one astronomical unit. The total solar irradiance (TSI) is the SSI integrated over the whole spectrum. Current knowledge of SI comes from space-based measurements and SI models.
Measurements collected since 1978 by a series of instruments show that the TSI varies at all discernible timescales from minutes to decades. Most prominent are a clear change by ∼0.1% in phase with the solar cycle (Kopp, 2016) and fluctuations by up to 0.2-0.3% on timescales shorter than a few days (Domingo et al., 2009).
On timescales longer than the solar cycle, TSI changes are rather uncertain due to the uncertainties in instrumental calibrations and the limited duration of the individual space experiments. In particular, there are currently four composite records of TSI measurements, presented by Willson (1997); Dewitte et al. (2004); Frhlich (2006); and Dudok de , showing partly conflicting secular trends over the last 4 solar cycle minima.
The TSI is the integral over the entire solar radiation spectrum that ranges from X-rays to radio waves. The bulk of solar radiant energy is in the visible (Vis, 400-700 nm), followed by the infrared (IR, 700 nm-1 mm) and ultraviolet (UV, 100-400 nm) bands. Measurements of the SSI also started in the 1970s, but are less continuous and more uncertain than those of TSI. Until just less than 20 years ago, they mainly covered the UV band (see Ermolli et al., 2013, for a review). Nevertheless, all present data clearly show that the form and amplitude of SSI variations are strongly wavelength dependent. Although the various experiments produced different estimates of the magnitude of UV changes, all available data agree on that they are in phase with the solar cycle (e.g. Floyd et al., 2003;Rottman et al., 2006). In contrast, SSI measurements at the IR and longer wavelengths show variations that are in anti-phase with the solar cycle. In regard to Vis, conflicting results from two experiments make SSI variations at this spectral range still uncertain. In particular, measurements from the Spectral Irradiance Monitor (SIM) on-board the SOlar Radiation and Climate Experiment (SORCE) satellite show Vis changes that are in anti-phase relation with the solar cycle (Harder et al., 2009), while those from the Variability of solar IRadiance and Gravity Oscillations (VIRGO) experiment on-board the SOlar and Heliospheric Observatory (SOHO) display an in-phase variation (Wehrli et al., 2013). However, it has been suggested that the SORCE/SIM data before 2010, which are the ones showing the anti-phase relation, are unreliable due to potential instrumental issues (Haberreiter et al., 2017;Mauceri et al., 2018Mauceri et al., , 2020. Models that ascribe variations in SI at timescales greater than a day to solar surface magnetism have proved to be successful in reproducing existing measurements of TSI and of SSI at given spectral ranges and timescales (e.g. Domingo et al., 2009;Shapiro et al., 2017;Yeo et al., 2014Yeo et al., , 2017b. There are various SI models presented in the literature. They can be grouped in the so-called proxy (e.g. Hudson et al., 1982;Chapman et al., 2013;Tebabal et al., 2015;Georgieva et al., 2015;Yeo et al., 2017a;Lean, 2018;Choudhary et al., 2020), semi-empirical (e.g. Fligge et al., 2000;Ermolli, 2001;Krivova et al., 2003;Crouch et al., 2008;Shapiro et al., 2010;Fontenla et al., 2011;Bolduc et al., 2012;Yeo et al., 2014;Wu et al., 2018), and physical (Yeo et al., 2017b) models. The proxy models use linear combinations of indices of solar magnetic features and regression to actual TSI measurements, while semi-empirical and physical models reconstruct the SI variations by employing results from spatially-resolved full-disc observations and from radiative transfer calculations performed on either atmosphere models or outcomes from magnetohydrodynamic simulations of the solar atmosphere, respectively.
Results from the various SI models show considerable disagreement in the long-term evolution of TSI with estimated changes since the Maunder minimum ranging from 0.7 Wm −2 to 6 Wm −2 (Wang et al., 2005;Krivova et al., 2007Krivova et al., , 2010Tapping et al., 2007;Steinhilber et al., 2009;Vieira et al., 2011;Shapiro et al., 2011;Judge et al., 2012Judge et al., , 2020Egorova et al., 2018;Wu et al., 2018;Lockwood & Ball, 2020). The various SI models also estimate different SSI variations on timescales longer than a few solar rotations (see Ermolli et al., 2013, for a review of this topic). In particular, for the SSI in the Vis range, most models suggest an in-phase variation with the solar cycle (e.g. Krivova et al., 2006;Ball et al., 2011;Unruh et al., 2012;Yeo et al., 2014Yeo et al., , 2017aMauceri et al., 2019), although there are also models indicating an anti-phase change (e.g. Preminger et al., 2011;Fontenla et al., 2011;Fontenla & Landi, 2018). In the UV, all SI models return variability in phase with the solar cycle.
Long and reliable time series of TSI and SSI are urgently needed to allow accurate quantification of solar contribution to the evolution of Earth's climate. To improve our understanding of the solar-cycle and long-term variability of SI, here we analyse the TSI and SSI variations in two spectral intervals in the Vis over the last two solar cycles derived from a careful analysis of two independent series of ground-based solar observations. In particular, we exploit the two longest time-series of modern full-disc photometric observations of the solar atmosphere, performed at the Rome and San Fernando observatories, and process these data with the method developed by Chatzistergos et al. (2018bChatzistergos et al. ( , 2019bChatzistergos et al. ( , 2020b to reconstruct the TSI and SSI variations over the period 1996-2020 by applying the proxy model by Chapman et al. (2013).
The paper is structured as follows. In Sect. 2 we describe the various full-disc solar observations and irradiance series used in our study, as well as the methods employed to process the images and to reconstruct SI variations. In Sect. 3 we analyse the solar cycle variability in different spectral bands covered by the full-disc observations and their phase relations. We present our results for the long-term variation of TSI over the last three solar activity minima (Sect. 4). We discuss the results in Sect. 5, before summarising them and drawing our conclusions in Sect. 6.

Photometric sums and irradiance modelling
We follow the approach by Chapman et al. (2012Chapman et al. ( , 2013, who used an empirical model based on the so-called photometric sums derived from Ca II K and continuum images. The assumption is that the irradiance variability is due to the modulation by dark and bright features, such as sunspots and faculae (in white light). Faculae are much better seen in chromospheric Ca II K observations, where they have a higher contrast than in photospheric observations and are called plage. We will refer to both faculae and plage as faculae, independently of the type of their observations. Sunspots have a higher contrast in the VIS and IR and therefore images in the red continuum in the 600-700 nm range are used to describe their effect. Faculae are much better represented by Ca II K observations (see e.g. Skumanich et al., 1984).
The photometric sum is defined as ). Here C i is the contrast value of pixel i defined as C i = (I i − I QS i )/I QS i , with I i being the intensity value of this pixel and I QS i the intensity of the quiet Sun (QS) regions. I CLV i is the average QS centre-to-limb variation (CLV) normalised such that at the disc centre it has the value of unity. Thus C i is dimensionless, and the photometric sum of quiet sun regions is equal to 0. By normalising I CLV i to unity we render the photometric sums derived from Ca II K and continuum observations analysed here directly comparable. We note, that the photometric sums as defined here depend on the pixel scale of the images. The images and their processing are described in the next section.
The TSI is calculated as a linear combination of the photometric sums in the red, Σ r , or blue, Σ b , continuum, and the one in the Ca II K line, Σ k (Chapman et al., 2013): The three parameters, a, b and c of the model are determined through a linear regression of the photometric sums to the actual TSI measurements described in Sect. 2.3

Full-disc observations
We analysed full-disc solar observations acquired with the Precision Solar Photometric Telescope at the INAF Osservatorio Astronomico di Roma (Rome/PSPT, hereafter) and with the Cartesian Full-Disk Telescope 2 at the San Fernando observatory (SFO/CFDT2, hereafter). Table 1 lists key characteristics of the analysed full-disc observations. In particular, we considered Rome/PSPT observations taken in the Ca II K line at 393.37 nm (0.25 nm bandwidth), blue continuum at 409.20 nm (0.25 nm bandwidth; blue hereafter), as well as in the red continuum at  [-0.5,0.5] for the continuum and Ca II K images, respectively. In the contrast images, regions outside the ranges listed above are saturated and appear black or white.
607.1 nm (0.5 nm bandwidth; red hereafter). The SFO/CFDT2 filtergrams were acquired in the Ca II K line at 393.37 nm (0.9 nm bandwidth). Regular monitoring of the solar atmosphere in Ca II K and the blue continuum with Rome/PSPT 1 started in May 1996. Observations in the red continuum started in September 1997 (Ermolli et al., 1998(Ermolli et al., , 2007. Two significant instrumental changes occurred in August 1997 and September 2001, involving the replacement of the CCD device and the Ca II K filter, and the telescope relocation from the Monte Mario to Monte Porzio Catone site, respectively. These changes introduced slight dis-1 Available at https://www.oa-roma.inaf.it/pspt-in-rome/ continuities in the series of the Rome/PSPT observations, which seem not to significantly affect the high photometric accuracy of the Rome/PSPT data (Ermolli et al., 2011;Chatzistergos et al., 2019b). The observations are stored as 1024×1024 pixel 2 images with 2 /pixel size. They are available as single exposure images (RPS, hereafter), as well as the sum of 25 exposures (RP, hereafter). The RPS data have better spatial resolution, while the RP images have an improved photometric accuracy. In this study we considered the RP data obtained from 16 May 1996 to 12 June 2020, and the RPS images acquired from 21 August 2000 to 12 June 2020. Even though RP images are sums of 25 exposures, individual exposures were not stored before the 21st of August of 2000. Notes. Columns are: name of the observatory, abbreviation used in this study, number of summed frames, central wavelength, spectral width, and period of observations, average pixel scale of the images, total number of available images, and the bibliography entry. The SFO/CFDT2 observations 2 analysed in this study cover the period 1996-2013. Regular solar monitoring started at the San Fernando Observatory with the Cartesian Full Disk Telescope 1 (SFO/CFDT1) in 1986. This telescope produces 512×512 pixel 2 images with ∼5.1 /pixel size. A second telescope, SFO/CFDT2, was introduced in 1992 ( Chapman et al., 2013) providing higher spatial resolution data. The images acquired with SFO/CFDT2 have dimensions 1024×1024 pixel 2 with ∼2.6"/pixel size. The images are the sum of two exposures taken roughly 7.5 minutes apart (Chapman et al., 2013). Here, we restricted our analysis to the SFO/CFDT2 data (SF, hereafter), because of their higher resolution compared to SFO/CFDT1, as well as potential instrumental issues with SFO/CFDT1 images (Chatzistergos et al., 2020b). We also restricted our analysis to SFO/CFDT2 data taken since 1996 due to the availability of the Rome/PSPT blue continuum observations.
All observations were calibrated by the observing teams for the instrumental response. We have processed the images with the methods described by Chatzistergos et al. (2016Chatzistergos et al. ( , 2018aChatzistergos et al. ( ,b, 2019aChatzistergos et al. ( ,b,c, 2020b,a, hereafter referred to as "our" method). This allowed us to consistently compensate the RP, RPS, and SF data for the intensity CLV of the QS regions and to account for any artificial intensity patterns affecting the data up to 0.99r, where r is the solar disc radius. By this compensation we obtained contrast images. A slightly reduced disc area up to 0.98r was used for the computation of the photometric sums in order to remove the effect of any possible artefacts very close to the limb. Figure 1 shows examples of the intensity images analysed in our study, as well as the corresponding contrast images obtained from our processing. Also shown (bottom panels) are profiles along the horizontal and vertical cross-sections of the contrast images (marked by the dashed green lines in the contrast images). In particular, these profiles traverse facular and sunspots regions on the contrast image. The contrast of faculae in the RP continuum is clearly lower than in RP Ca II K observations. In the RP blue, the contrast is higher than in the RP red. In both RP continuum bands, the contrast in sunspots is considerably higher than in faculae. Sunspot contrast is also higher compared to that of faculae in SF Ca II K, though the difference is less pronounced than in the RP continuum data. In the RP and RPS Ca II K data, 2 Available at http://www.csun.edu/ SanFernandoObservatory/sfosolar.html the facular contrast is higher than the sunspot one. The facular contrast in RP and RPS data is higher than in SF Ca II K data. These relations between contrast values in the various data are consistent with the narrower bandwidth used for the RP and RPS observations compared to that employed for the SF data.

TSI series
In this study, we considered two original instrumental TSI measurement series, three TSI composites, and three model TSI reconstructions. The original instrumental records are the daily TSI measurements from the SORCE/TIM 3 (version 18, February 2020; Kopp & Lawrence, 2005) and SOHO/VIRGO 4 (version 6.5, May 2018; Frhlich et al., 1997) covering the time intervals from 2003 to 2020 and from 1996 to 2018, respectively. The considered TSI composites are the ACRIM 5 (Active Cavity Radiometer Irradiance Monitor, which is the instrument taken as the reference by Willson, 1997), PMOD 4 (named after Physikalisch-Meteorologisches Observatorium Davos; version 42.65; Frhlich, 2006), and RMIB 6 (named after Royal Meteorological Institute of Belgium, in french called IRMB; Dewitte et al., 2004;Dewitte & Nevens, 2016), covering the periods 1978-2013, 1978-2018, and 1984-2020, respectively. These TSI composites are created by daisy-chaining, i.e. using the daily measurements of one instrument as the reference to calibrate the daily records of the second one, which then acts as the reference for the next one, etc. Finally, to compare the results of our study with others, we also considered the daily TSI reconstructions with the SATIRE-S 7 (Spectral And Total Irradiance REconstruction; Yeo et al., 2014)  (bottom) observations. Daily values are shown in black, while 180-day running averages are shown in red. We note that the photometric sum indices were adjusted, following the method by Nesme-Ribes et al. (1996), to take into account the bias of the QS level (see Section 3). The dashed orange line shows the 180-day running average of the photometric sums prior to this adjustment, while the dash-dotted ciel line (only for the continuum data) was obtained following the approach by Preminger et al. (2002). The dashed blue horizontal line denotes the 0-level of photometric sum.
All these models describe the irradiance variability through the varying contributions of sunspot darkening and facular brightening. SATIRE-S is a semi-empirical model. It employs the areas and positions on the solar disc of the various features at a given time as well as their time-invariant brightness contrasts as input. The brightness contrasts as a function of wavelength and position on the solar disc are obtained with the radiative transfer simulations from semi-empirical atmospheric models (Unruh et al., 1999). The disc area coverage by the surface features is derived from magnetograms and continuum filtergrams acquired from ground-and space-borne telescopes since 1974 (see e.g. Yeo et al., 2014, for more details). EMPIRE and NRLTSI are empirical models that linearly combine proxies of sunspot dark-ening and facular brightening to match the model outcome to the observed TSI changes. Whereas NRLTSI uses the ordinary least squares regression, EMPIRE uses the orthogonal distance regression accounting for the error in both variables. Both models use the Mg II index and sunspot areas as input over the period considered in this study.
3. Variation of the photometric sums with the solar cycle Figure 2 shows the Σ k , Σ b , and Σ r series obtained from the processing of the RP Ca II K, blue, and red images. All of them vary with the solar cycle, although with different amplitudes. In particular, the amplitude of the variation of Σ b and Σ r is considerably lower than that of Σ k , with the maximum absolute value of ∼ 1500, which is an order of magnitude lower than that of Σ k . The orange line shows the originally derived Σ k , Σ b , and Σ r . As was already suggested by Fig. 1, the contribution of facular regions in Ca II K images is considerably higher than that of sunspots, leading to an in-phase variation of Σ k with the solar cycle, i.e. the Sun in the Ca II K line appears brighter during periods of activity maxima. In the red and blue continuum images, where the contrast of faculae is generally lower than that of sunspots, the balance is rather shaky and Σ b and Σ r might vary in anti-phase with the solar cycle, i.e. the Sun in the specific continuum intervals might appear marginally darker during periods of activity maxima. This is in agreement with the result of Preminger et al. (2011). However, as discussed below and as also argued by Peck & Rast (2015), this result turns to be extremely sensitive to the image processing, in particular, to the accuracy of the determination of the CLV.
In the following we assess the potential uncertainties in the derived photometric sum indices. We use two different approaches to identify the QS level on the contrast images and to evaluate its variation with time. This allows us to test whether the processing had introduced any bias in the QS level. Here, we focus on the continuum data since such a test on Ca II K observations has already been performed by Chatzistergos (2017) and Chatzistergos et al. (2018b) using synthetic images.
Both methods make the underlying assumption that the distribution of the image contrast values is Gaussian, with active regions contributing to the wings of the distribution, but they take different approaches to identify the mean level of the QS regions. The first method was proposed by Nesme-Ribes et al. (1996) and was also used by Chatzistergos et al. (2019bChatzistergos et al. ( , 2020b). An iterative process is used to estimate the level of the QS. Initially, the mean contrast,C, and the standard deviation, σ, of the contrast values within the disc are computed. Then the mean and the standard deviation are re-calculated iteratively, each time leaving out regions with values outside the intervalC ± kσ, where k is a constant. The value of k is taken between 0.1 and 3.0, and using a bisector approach, we search for the lowest value of the mean contrast among all k values. This lowest mean contrast value is then adopted as the one representing the QS.
The second approach is an adaptation of the method used by Preminger et al. (2002) to analyse the potential bias in the SF data. For this, we compute the histogram of contrast values  . At the bottom of each panel, the corresponding difference between the PMOD TSI and the reconstructed TSI (i.e. PMOD minus our reconstruction) is shown (black circles for daily values and orange line for the 180-day running mean). TSI series are reconstructed using the following data: RP Ca II K and RP blue (a); RP Ca II K and RP red (b); RPS Ca II K and RPS red (c); and SF Ca II K and RP red (d) observations. For all reconstructions, the PMOD TSI composite was used as the reference when deriving the regression coefficients.
for each RP red and blue image, where the contrast values are grouped into bins of 10 −5 width. A consistent image processing would return profiles centred at a value constant over the solar cycle. Following Preminger et al. (2002), we fit two Gaussians to the histogram of contrast values for each processed image, one for the central part, representing the QS regions, and one for the wing representing the bright regions. To confine the central part, we iteratively reduce the considered interval and fit a Gaussian function, searching for the contrast range for which the lowest χ 2 of the Gaussian fit is achieved. The location of the peak of the central Gaussian is taken to represent the mean level of the QS regions. The residual between the distribution and the central Gaussian in the bright wing is then fit with another Gaussian function. Preminger et al. (2002) performed the fit on 7 annuli equally spaced in µ rather than the whole disc. They also did not iterate. Despite this modification, we will still refer to this approach as the one by Preminger et al. (2002).
As an example, Fig. 3 shows the histograms of the RP red and blue images presented in Fig. 1 along with the evolution of the QS level computed with the two approaches mentioned above for all continuum images analysed in this study. Both approaches indicate a weak although non-negligible variation of the QS level in anti-phase with the solar cycle. The QS level derived with the method of Preminger et al. (2002) varies by more than Fig. 5. Differences between the various reference TSI datasets (see the legend) and the corresponding TSI reconstructions from RP Ca II K and blue observations. The residuals have been offset to improve visibility. The dotted red lines denote differences of 0, while the dashed lines are linear fits to the residuals. when derived with the method of Nesme-Ribes et al. (1996). The mean and the standard deviation of the QS level obtained with the method of Preminger et al. (2002) are −1 ± 4.9 × 10 −4 and −2.5±3.8×10 −4 for the RP red and blue images, respectively. The mean and the standard deviation of the QS level determined with the method of Nesme- Ribes et al. (1996) are −0.7 ± 0.3 × 10 −4 and −1.2 ± 6.5 × 10 −4 for the RP red and blue images, respectively. These findings suggest that the QS level in the RP red and blue images calculated with both approaches is affected by the solar activity. In particular, the QS level is underestimated during activity maxima. Around activity minima, both methods return relatively accurate results. The variation is stronger in the RP blue than in the RP red images. We found no dependence of the bias on the heliocentric angle. However, when applied on narrow annuli, the results are less accurate for both methods. The reason is the significantly reduced statistics compared to the case when the entire disc is considered. It is worth noting that the results obtained with the method by Nesme-Ribes et al. (1996) are qualitatively similar to those obtained by Chatzistergos (2017) from the analysis of Ca II K data, for which we found a bias of ∼ 1±1×10 −4 with a minute anti-phase relation to the solar cycle.
In Chatzistergos (2017) we showed that the accuracy of recovering the pixel-by-pixel contrast values in synthetic Ca II K images was on average 10 −3 (the RMS error reached up to 3 × 10 −3 ). We also showed that our approach worked better than other published methods. This pixel-by-pixel error is one order of magnitude higher than the QS bias reported here with both methods. These processing errors can affect the determination of the QS level. To estimate the sensitivity of the applied methods on the pixel-by-pixel errors of our processing, we performed a Monte Carlo simulation. In particular, we used all RP contrast images and imposed a random offset on contrast values in each pixel. We then recomputed the QS level with the two considered methods. We have done this for 100 different sets of uniformly distributed random contrast offsets within the range -3-3×10 −3 . The distribution of the differences between the QS level computed from the offset and from the original images gives us an estimate of the uncertainty of the QS bias estimation due to processing errors. These distributions for annual bins are shown in Fig. 3. We find the uncertainty due to the pixel-by-pixel errors to typically lie below 10 −4 and 3 × 10 −5 for the methods by Preminger et al. (2002) and Nesme-Ribes et al. (1996), respectively. This suggests that the estimation of the QS bias is not strongly affected by the pixel-by-pixel processing errors of our method. The uncertainty with both methods is comparable, although the effect is slightly stronger when using the method by Preminger et al. (2002). The derived QS bias for RP blue data is higher over the years 1996 and 1997, thus rendering our results over these two years rather uncertain. We note, however, that this uncertainty estimate of the QS bias is only valid for images processed with our method and it is expected to be higher when using other published processing methods.
The results presented in Fig. 3 suggest that to accurately determine the photometric sums Σ b and Σ r , we need to account for the bias in the level of the QS definition. Therefore, we repeated the computation of the Σ series after the appropriate adjustment of the QS level in each image, using both the methods by Preminger et al. (2002) and Nesme-Ribes et al. (1996). These adjusted photometric sum series are shown in Figure 2. The RMS difference between the Σ series adjusted with the method by Nesme-Ribes et al. (1996) and the original Σ series is 115 for both the red and the blue data. Using the method by Preminger et al. (2002) the RMS differences are 294 and 384 for the red and blue observations, respectively. For comparison, the corresponding RMS difference for the Ca II K data when using the method by Nesme-Ribes et al. (1996) is 146. The adjustment to account for the bias of the Σ k series has barely any effect on the series since the adjustment is two orders of magnitude lower than the Σ values. However, for the blue and red observations, the correction is comparable to the Σ values.
After the correction, the temporal profile of Σ b and Σ r series changes both in terms of the amplitude and the shape (see Fig. 2). In particular, with the QS level corrected with the method by Nesme-Ribes et al. (1996, red curves) the profile of Σ b becomes essentially flat, with only a tiny decrease over the maximum of solar cycle 24, whereas the variation of the Σ r remains in anti-phase with the solar cycle. After the correction of the QS level with the method by Preminger et al. (2002) the variation of Σ b turns to be in-phase with the solar cycle, while the profile of Σ r becomes essentially flat. However, while the adjusted Σ r and Σ b are, on average, close to zero at minima when the method by Nesme-Ribes et al. (1996) is used, the corresponding values for Σ b with the corrections following Preminger et al. (2002) are roughly around 100. This points to a potential over-correction with the method by Preminger et al. (2002). We thus conclude, that the uncertainty in the definition of the QS bias, independently of the method used, as well as the potentially remaining dependence of the QS bias on solar activity render the phase of the variation of Σ r and Σ b significantly uncertain. We consider the method by Nesme-Ribes et al. (1996) to return more accurate results than the method by Preminger et al. (2002) as applied here. In the following we will use the Σ series adjusted for the QS bias with the method by Nesme-Ribes et al. (1996) being the intermediate case found here.
We note that Preminger et al. (2002) also reported a variation of Σ r derived from SF red continuum images similar to what we obtain here. However, they argued that the offset was too small to have any effect on the solar-cycle phase of the Σ r variability. Specifically, they found a mean bias in the QS level during an activity maximum of only ∼ 5×10 −4 , although this offset doubled towards the limb. They argued that this bias introduced a mean offset in Σ r of 120 and -13 for 1989 and 1996, corresponding to activity maximum and minimum periods, respectively. These values are consistent in amplitude with the ones derived here.
Puiu (2019) analysed the accuracy of Σ evaluation by processing a small sub-sample of RPS observations with methods different than those employed in this study (see Sect. 5 for a brief description of their method). In particular, they studied the effect of various parameters of the processing on their results and especially on the photometric sums. Their final adopted parameters returned contrast images favouring the variation of Σ r in anti-phase with the solar cycle. However, their analysis showed that generally the image processing can result in a reversal of the phase of the relation between Σ r and the solar cycle.

Long-term trends in reconstructed TSI
Using the computed photometric sums and Eq. 1, we now reconstruct the TSI. The free parameters of the multiple linear regression model are determined by fitting various combinations of the photometric sums, Σ, to the PMOD TSI composite. We consider the Σ series derived from the following 5 pairs of observations: RP Ca II K and RP red, RP Ca II K and RP blue, RPS Ca II K and RPS red, as well as SF Ca II K and RP red and SF Ca II K and RP blue. Figure 4 shows the reconstructed TSI along with the PMOD TSI composite. Table 2 summarises the parameters of the fit, along with the RMS difference, and linear correlation coefficient, R, for each reconstruction. We found that the best TSI reconstruction is achieved with RP Ca II K and RP blue observations, with a RMS difference and R to the PMOD TSI series of 0.15 Wm −2 and 0.96, respectively. The TSI reconstructions from RP red observations exhibit a slightly higher scatter than those derived from RP blue data. The quality of the TSI reconstructions is similar to the previous combinations when RPS or SF Ca II K data are used. The PMOD TSI shows a somewhat stronger declining trend than the reconstruction. This is discussed in more detail below.
The empirical TSI reconstructions, based on the regression of solar activity indices to the TSI measurements, are obviously subject to the uncertainties in those measurements. Therefore, we have done various TSI reconstructions by varying the TSI reference series used for the regression. Table 2 lists the parameters of the fit, along with the RMS difference, and R for each reconstructions, whereas Fig. 5 shows the residual between each considered TSI reference record and the corresponding reconstruction from RP Ca II K and RP blue observations. We generally found a good agreement between our TSI reconstructions and the various published TSI series (see Sect. 2.3). Most reconstructions show the same characteristics. In particular, as in the case of PMOD, a slightly declining trend is seen in the residuals with time. Exceptions are the reconstructions that use the SORCE/TIM data as reference, which display a slightly increasing trend of 1.3 × 10 −2 Wm −2 y −1 . This might partly be because the SORCE/TIM data do not cover the maximum of cycle 23. However, we note that Dewitte & Nevens (2016) and Woods et al. (2018) reported an increasing drift of 3.43 × 10 −2 Wm −2 y −1 in the TSI values from SORCE/TIM over the same period. We also found that our TSI reconstruction with the ACRIM composite as the reference displays the strongest trend of residuals to the measurements. This is most likely due to issues with the ACRIM composite series. We also computed the R and RMS differences between our reconstructed TSI and the various reference series after detrending our series. For detrending we used the residual trend listed in Table 2. We find the RMS differences to generally decrease slightly in all cases, however qualitatively we report the same findings.
For a more quantitative comparison, we have performed a linear fit to the reconstructed TSI over each two subsequent solar activity minima (that is 1996 vs. 2008 and 2008 vs. 2019). For this, we used 12-month intervals centered at the minimum of each of the 3 cycles. The dates of the minima were taken as 01 August 1996, 01 December 2008, and 01 October 2019 (Hathaway, 2015). We note that at the time of writing this manuscript there is no official starting date for solar cycle 25 and hence the value used here might change somewhat, although probably not much (see e.g. Alterman et al., 2020). Therefore, we use the last year for which the data are currently available for the last minimum. Figure 6 shows the reconstructions for the various TSI measurement records employed as the reference, along with the linear fits to the TSI values during the activity minima. For ease of comparison, we offset the reconstructed series so to match their mean values over a 12-month interval centred at the minimum of 2008 to that of the PMOD TSI series. While the various reconstructions differ slightly in the magnitude of the solar cycle variation, they agree remarkably well during activity minima, showing essentially the same minimum-to-minimum changes. Thus, the choice of the reference TSI dataset is not decisive for the long-term trend in the reconstructed record. We also note that the long-term trends of the TSI series derived from the PMOD and RMIB TSI composite data employed as the references are remarkably close to each other, with differences being less than 2.1 × 10 −4 Wm −2 y −1 . However, this is mainly due to the similarity of the PMOD and RMIB TSI composites over the selected periods.
The TSI reconstruction based on RP Ca II K and RP blue images that uses the PMOD data as the reference shows a long-term trend of -7.8 +4.9 −0.8 × 10 −3 and -0.1 +0.25 −0.02 × 10 −3 Wm −2 y −1 between the minima 1996/2008 and 2008/2019, respectively, where the uncertainty ranges account for the spread in results derived using all other published TSI data as the reference. The uncertainty of the fit is roughly 10 −3 Wm −2 y −1 and is lower than the uncertainty due to using different TSI reference series. The declining trend over the 1996/2008 minima is statistically significant based on a Student's t-test, but this is not the case for the declining trend over the 2008/2019 minima. However, using different input data returns slightly different values for the long-term trend. For example, for the change between 1996 and 2008 we found −11.9 +1.2 −0.06 ×10 −3 Wm −2 y −1 when using SF Ca II K and RP blue observations, which is slightly higher than obtained from RP Ca II K and blue data. Between the 2008 and 2019 minima, the trend marginally changes sign and is 2.6 +0.3 −0.1 ×10 −3 Wm −2 y −1 with RPS Ca II K and red observations, while it becomes 1.6 +0.7 −0.1 ×10 −3 Wm −2 y −1 with RP Ca II K and red observations. Overall, our reconstructions suggest a weak decline in TSI from 1996 to 2008, and essentially no trend from 2008 to 2019.
Finally, for completeness, we have also compared our reconstructions to the widely used models, namely SATIRE-S, NRLTSI and EMPIRE (see Sect. 2.3 for the description of the datasets and Table 2 and Figs. 5 and 6 for the results). Thereby, for an unbiased comparison, we also used the corresponding series as the reference. The lowest RMS difference and highest R are obtained in comparison to SATIRE-S and EMPIRE reconstructions. The slope of the residual is somewhat closer to zero when the two empirical reconstructions, EMPIRE and NRLTSI, are used as reference. These models are similar in concept to ours, i.e. they are also based on linear regressions of the solar activity indices to the measured TSI, except that they use discintegrated proxies. The slope of the residual between our reconstruction and SATIRE-S is closer to that obtained for most recon-structions using observational records as the reference (except SORCE/TIM). All differences are, however, rather small.

Discussion
In Sect. 3 we analysed the sensitivity of the Σ series to the computation of the QS level in the images when processed with our method. We now consider also the effect of the QS computation on the Σ series and the TSI reconstruction when using different processing methods. In particular, we consider the methods presented by Brandt & Steinegger (1998), Walton et al. (1998), Worden et al. (1998), andPuiu (2019). All of these studies used full-disc Ca II K observations to reconstruct TSI.
Briefly, Brandt & Steinegger (1998) and Walton et al. (1998) derive the QS CLV by considering image pixels grouped into concentric annuli, with the value of the QS within each annulus determined with a threshold in the cumulative histogram of intensity values and the median of intensity values, respectively. Since both of those methods assume that the CLV is perfectly radially symmetric, they are unable to account for some image artefacts that often affect the full-disc observations, such as the linear gradient of intensity values caused by the atmospheric extinction in the RP and RPS data. To account for this, we have added an extra step to perform a 2D linear fit to the contrast images resulting from the processing by Brandt & Steinegger (1998) and Walton et al. (1998). This fit is then subtracted from the contrast images. The processing by Worden et al. (1998) incorporates 5th degree polynomial fits along lines in the vertical and horizontal directions, as well as in segments with the orientation of 45 • and 135 • . Finally, Puiu (2019) calculates the background as Walton et al. (1998), but regions within and outside 0.7r are processed separately. Additionally, active regions are excluded from the QS computation with an iterative procedure that follows the one proposed by Chatzistergos et al. (2018b) and used in this study. As the last steps, Puiu (2019) applies a linear vertical fit and a running window median filter to account for any residual intensity gradient in RPS data.
In Chatzistergos (2017) and Chatzistergos et al. (2018b) we studied the accuracy of the methods proposed by Brandt & Steinegger (1998), Walton et al. (1998), andWorden et al. (1998) to remove the limb-darkening on synthetic Ca II K data. We showed that these methods performed worse than the ones employed in this study, by introducing artefacts that are subject to the solar cycle due to an inaccurate exclusion of active regions when defining the QS CLV. In particular, when using the methods by Brandt & Steinegger (1998), Walton et al. (1998), andWorden et al. (1998), mean offsets for the QS level were up to 2, 1, and 2 orders of magnitude higher than the one obtained with the method employed here, respectively. Furthermore, Puiu (2019) compared the accuracy of their processing to ours and also found the latter to be more accurate.
It is important to note that these error estimates were limited to Ca II K images, for which a significant source of errors is inaccurate accounting for active regions. On continuum images, where the disc coverage by active regions is significantly lower than in Ca II K images, other methods might fare better. Therefore, to check this, we used all these four methods to pro-cess the RPS Ca II K and red images, produce photometric sum series and reconstruct the TSI. We restricted this test to the RPS data for consistency with the analysis by Puiu (2019). Figure 7 (left panels) shows the Σ k series derived from RPS Ca II K observations processed with our and three of the earlier methods. We do not show the results with the method by Brandt & Steinegger (1998), because they are very similar to those by Walton et al. (1998). The Σ k series obtained from our method exhibits the highest variability over the solar cycle. This suggests a potential suppression of the contrast of the bright regions by the other processing techniques, which is in agreement with the finding that these methods do not accurately exclude active regions when computing the QS CLV (Chatzistergos et al., 2018b). Figure 7 (right panels) shows the Σ r series derived from RPS red images by using our and three of the tested methods. All series exhibit a variation in anti-phase with the solar cycle. However, the series obtained with the methods by Brandt & Steinegger (1998) and Walton et al. (1998) show high scatter and include many extreme values that are likely due to image artefacts unaccounted for by the image processing. For this reason, we ignored 77 images for which the processing by the method of Walton et al. (1998) produced Σ r values >2000. We note that these results were obtained by applying the 2D linear fit, which might not accurately remove the linear gradient in the data. However, the results without the additional 2D fit were considerably worse, lacking any clear cycle variation. We also note that the 5th degree polynomial used for the fitting along linear segments in the method by Worden et al. (1998) fails to reproduce the CLV in RPS red observations. This is in agreement with the findings by Puiu (2019). We found that the method by Puiu (2019) results in a distinct anti-phase relation between the Σ r series and the solar cycle. The Σ r series computed with the method by Puiu (2019) shows the best agreement among the various tested methods to that derived with our method, exhibiting a high correlation coefficient of 0.88 (see Table 3). This suggests that the method by Puiu (2019) performs better than the methods by Brandt & Steinegger (1998), Walton et al. (1998), andWorden et al. (1998). The spread of Σ r values during minima is greater when using the method by Puiu (2019) than with our method. Table 3 summarizes the fit parameters (Eq. 1) when using the Σ series derived from all tested methods. Also listed are the RMS difference, R, the slope of the residuals to the reference series, as well as the long-term trend between the last solar cycle minima in 2008 and 2019. Figure 8 shows the TSI variations reconstructed from RPS Ca II K and red data processed with the various tested methods. The TSI reconstruction from RPS data processed with our method is plotted in Fig. 4c.
TSI reconstructed from the images processed with our method scores better than when applying other processing methods tested here, followed by the method by Puiu (2019) with slightly worse RMS differences to the PMOD series. All other processing techniques give worse match to the TSI reference series, with RMS differences increasing to 0.35-0.4 Wm −2 instead of 0.17 Wm −2 obtained with our method. Interestingly, the longterm change in TSI between the minima in 2008 and 2019 returned by different methods shows opposite signs. Specifically, the trend in TSI determined with the processing following Puiu (2019) and Worden et al. (1998) is negative (decreasing TSI), and

Notes.
Columns are: Ca II K and continuum observations employed to compute Σ k , Σ b , and Σ r , TSI reference series used for the regression, best fit parameters of Eq. 1, RMS difference, linear correlation coefficient R, and slope of residual between the reconstructed TSI series and the reference one, and long-term trend as defined with a linear fit during activity minimum periods.
These tests suggest that for reliable irradiance reconstructions, accurate processing of the images is crucial. Uncertainties of the processing methods might considerably affect the determined long-term trends in the reconstructed TSI.

Summary
Solar irradiance has been measured from space for over four decades now. Despite significant progress in our understanding of the irradiance variability, its long-term trend is still under debate. To address this question, we have reconstructed the TSI over the period 1996-2020 covering the last two solar cycles.
We have used an empirical model which computes the irradiance changes through a linear regression of solar activity indices. We have followed the approach by Chapman et al. (2013) who used the so-called solar photometric sums as input prox- Fig. 6. TSI reconstructed with various reference series using the following images: RP Ca II K and RP blue (a); RP Ca II K and RP red (b); RPS CA II K and RPS red (c); and SF Ca II K and RP blue (d). Colours are the same as in Fig. 5. The solid line is 81-day running mean values for the series with PMOD as the reference, while the grey shaded surface shows the range covered by the reconstructions with all other TSI reference series. The dashed and dotted lines are linear fits to the reconstructed and reference TSI values, respectively, during subsequent activity minima. All series were offset to match their mean value within a 12-month interval centered at the minimum of 2008 to that of the PMOD TSI composite. The value of the PMOD TSI composite over this period is marked with the horizontal dash-dotted pink line. The various dashed and dotted lines extend only to the intervals for which their respective data cover the complete two surrounding minima. Notes. Columns are: processing technique applied on the RPS Ca II K and RPS red continuum data to compensate the CLV, linear correlation coefficient R 0 between the Σ k and Σ r series computed with our processing and those from the various tested methods, the best fit parameters of Eq. 1, RMS difference, linear correlation coefficient R, and slope of residual between the reconstructed TSI series and the PMOD TSI composite, and long-term trend as defined with a linear fit during the activity minimum period of 2008-2019.
ies. The photometric sums were produced from full-disc Ca II K, red, and blue continuum observations acquired with the Precision Solar Photometric Telescope at the INAF Observatory of Rome and Ca II K observations taken with the Cartesian Full-Disk Telescope 2 at the San Fernando Observatory.
We find a good agreement between our reconstructed TSI series and the direct TSI measurements by SOHO/VIRGO and SORCE/TIM, as well as the RMIB and PMOD TSI composites. The agreement is worse with the ACRIM TSI composite. We also find a good agreement with the empirical NRLTSI and EMPIRE models, as well as with the semi-empirical SATIRE-S model, whereby the agreement with SATIRE-S and EMPIRE is somewhat better.
We find a weak declining trend of -7.8 +4.9 −0.8 × 10 −3 Wm −2 y −1 between the 1996 and 2008 minima. Between the minima in 2008 and 2019, the trend is -0.1 +0.25 −0.02 × 10 −3 Wm −2 y −1 , i.e. the TSI might have decreased slightly although this change is within the error-bars. Furthermore, measurements from SORCE/SIM (Harder et al., 2009) and from SOHO/VIRGO (Wehrli et al., 2013) suggested conflicting solar-cycle trends in the irradiance variations in the visible part of the spectrum. Based on the analysis of photometric sum series from SF red and blue observations, Preminger et al. (2002) argued that the irradiance variations in the visible are in anti-phase with the solar cycle, in agreement with the SORCE/SIM measurements. Here we used RP observations in the blue and red continuum to produce series of photometric sums and study the phase of their variations. By carefully analysing the accuracy of the image processing applied to the data, we found a small bias in the QS level of RP red and blue continuum observations of the order of 5×10 −4 in anti-phase with the solar cycle. We argue that this bias renders the photometric sum indices rather uncertain, such that the correlation to solar activity can change from being negative to marginally positive after the correction. The amplitude of the variations in the visible is very low and is below the uncertainties, which prevents an assessment of the phase of the variations. Results obtained with alternative techniques are even more uncertain.
This approach to irradiance reconstruction is unfortunately limited to a more recent period only, due to the lack of narrowband continuum observations before the 1980's. Further modifications to the method, e.g. accounting for sunspots using alternative data, such as white light observations (e.g. Willis et al., 2016;Pal et al., 2020), could allow irradiance reconstructions from properly calibrated historical Ca II K observations, e.g. those from Meudon (Malherbe & Dalmasse, 2019) or Kodaikanal (Chatzistergos et al., 2019c) observatories, to extend the reconstruction back to 1892 (Chatzistergos et al., 2020b).  Walton et al. (1998, panel c);and Worden et al. (1998, panel d). The PMOD TSI composite was taken as reference for all those reconstructions. Find more details in Section 5. At the bottom of each panel, the corresponding difference between the PMOD TSI and the reconstructed TSI is shown (black circles for daily values and orange line for the 180-day running mean).