© T. Dudok de Wit, Published by EDP Sciences 2022

1 Introduction

One of the major challenges in solar irradiance studies is the production of long records that are accurate and stable in time. Climate modellers, for example, rely on the solar UV input into the atmosphere to quantify ozone production in the tropical stratosphere, which is one of the levers by which solar variability affects climate (Gray et al., 2010). Without accurate expressions of instrument errors, decadal climate change signals cannot be confidently assessed (Ohring et al., 2005; Fox et al., 2011; Wielicki et al., 2013; Kopp, 2014). Solar physicists need to understand long-term solar variability to properly constrain past solar variations (Coddington et al., 2016), while satellite operators need an accurate UV flux to estimate satellite drag (Vourlidas & Bruinsma, 2018). Of particular interest for all these users is the spectrally-resolved solar irradiance, or solar spectral irradiance (SSI), because different wavelengths impact the Earth’s environment in distinct ways.

The need for radiometric accuracy and stability is challenged by the difficulty in making SSI observations from space, where instruments suffer continuous degradation of instrumental sensitivity and have a limited lifetime (BenMoussa et al., 2013). Different strategies have been developed to mitigate these adverse effects (see Woods et al., 2018 and references therein). A common strategy consists of duty cycling between identical sensors that are not exposed to the Sun for the same amount of time. Corrections can then be applied by assuming the degradation to be proportional to accumulated exposure time or to the fluence, and also to mission time (Kopp et al., 2005; Harder et al., 2009; Mauceri et al., 2020). A second approach consists in performing in-flight calibration to track instrument degradation, either by using stable stars (e.g. Snow et al., 2010), internal lamps (e.g. Meftah et al., 2016), independent observations at a different wavelength (e.g. Wehrli et al., 2013) or of an integral quantity, such as the total solar irradiance (TSI) to help constrain the SSI (e.g. Mauceri et al., 2018). Finally, sounding rockets with a copy of the instrument can be used (e.g. Didkovsky et al., 2012).

In spite of these solutions, the routine production of accurate SSI data remains a complex and demanding task. The consequences of our incomplete understanding of instrument degradation are well illustrated by the anomalous solar cycle amplitude that has been reported for the Spectral Irradiance Monitor (SIM) onboard SORCE (Harder et al., 2009) and the debate it has fuelled (Lean & DeLand, 2012; Deland & Cebula, 2012; Ball et al., 2016b). In this context, any additional constrains on SSI observations is a precious ally. This is our prime motivation for proposing a different approach for detecting and quantifying undocumented trends.

1.1 Stability versus solar variability

A simple metric for determining the significance of a trend consists in comparing its magnitude to the level of solar variability. The two most conspicuous time scales of interest are the solar or Schwabe cycle (~11 years) and the modulation due to solar rotation (~27 days). Figure 1 compares the relative amplitude of these two modulations, as estimated from observations made by the Spectral Irradiance Monitor (SIM, Harder et al., 2010) and by the SOLar-STellar Irradiance Comparison Experiment (SOLSTICE, Snow et al., 2010). In both cases, the modulation amplitudes increase at short wavelengths (Lean, 1991). While solar rotational effects are too fast for instrumental trends to be a concern, solar cycle modulations are harder to disentangle from trends. As a rule of thumb, the proper detection and quantification of a decadal trend from observations that span one solar cycle require the magnitude of the trend to be comparable to or larger than that of the solar cycle. Consider, for example, the 2% cycle amplitude at 200 nm in Figure 1. Trends at that wavelength can be detected if their magnitude typically exceeds 2% over an 11 year period, or 0.2% per year if one assumes a linear trend. Unfortunately, this value is also comparable to the stated instrument stability, which illustrates how delicate the production of irradiance data can be.

Fig. 1 Relative amplitude of the solar cycle variation and the solar rotational modulation, as estimated from SORCE/SIM and SORCE/SOLSTICE measurements between May 2003 and April 2019. SOLSTICE observations are used below 300 nm, and SIM observations above. To be less affected by non-physical trends, the amplitudes are estimated by Fourier analysis and not by taking the usual difference between extrema.

1.2 What trends?

So far, we have used the word trend in a loose sense and indeed, there is no crisp definition for it. The primary origin of instrumental trends is the gradual decrease of instrument sensitivity over time. Trends, however, are not necessarily monotonic since there may be phases of sensitivity increase (e.g. Ball, et al., 2016b). Discontinuities can occur, for example, following a solar proton event or after a sudden change in the properties of the instrument.

A considerable amount of literature has been devoted to detecting such trends. Trend detection is a major issue in the homogenisation of climate records (Chandler & Scott, 2011). Applications can also be found in the detection of ozone trends (Staehelin et al., 2001; Laine et al., 2014; Ball et al., 2017), ionospheric properties (Bremer, 2005; Laštovička, 2017), solar UV radiation (Lean, 2010; Morrill et al., 2014; Woods et al., 2018) and in TSI or total solar irradiance (Fröhlich, 2009; Kopp et al., 2012).

While the solutions that have been proposed are context-dependent, most of these studies have some features in common:

1. In the absence of a rigorous definition, a trend denotes a drift or prevailing tendency in the instrument response occurring on time scales that are comparable to or larger than the duration of the observations. In the case of SSI observations, trends typically occur on time scales of years.

2. The change in the instrument response does not have to be monotonic or regular, yet regularity is frequently assumed. The detection of abrupt changes, better known as changepoint detection (Reeves et al., 2007), is a different problem that goes beyond the scope of our study.

3. Trends are also affected by measurement uncertainties. Leroy et al. (2008) provide estimates for the minimum duration to detect a trend using a real observing system.

The majority of trend detection techniques are based on multiple regression by modelling variations in the observable interest with one or several solar proxies. The difference (or residual error) between the model fit and the observations should then capture the trend if there is any. Implicit to this assumption is the good stability of the model inputs, i.e. the absence of non-solar trends in the proxies. A more delicate issue, which is rarely addressed in the literature, is whether the proxy model is actually capable of revealing the trend. For example, a linear model may fail to capture the nonlinear relationship between proxy and irradiance or the presence of hysteresis. For these reasons, it is important to find a model that makes the least possible commitment to the connection between proxies and irradiance. In what follows, we shall introduce a non-parametric model that goes in that direction. While this model cannot prove the existence of a trend in an unambiguous way, it helps provide stronger evidence against or in favour of one. Our objective is to reveal trends, and not to seek to model the irradiance.

1.3 Trend detection in the SSI

There is a long history of trend correction in the SSI, and even more so in the TSI. The usual strategy involves one or several of the aforementioned methods, namely the use of redundant sensors, in-flight calibration and multiple regression with solar proxies. All these methods involve some modelling and therefore come with assumptions. For example, Mauceri et al. (2020) consider piecewise linear fits to compare instrument channels that have different duty cycles. The main challenge is to include the uncertainty that stems from these assumptions in the correction. While this is relatively natural in a Bayesian setting (d’Agostini, 2003), it is much more difficult to find in the literature realistic uncertainty estimates associated, for example, with regression models of solar proxies. A model that matches the SSI remarkably well may have large error bars because of a lack of uniqueness. In this context, what is missing is a more transparent approach, in which the methods are made publicly available, and the results come with meaningful uncertainties, allowing others to replicate the outcome and test its significance. The method we propose belongs to the category of semi-parametric regressions since it relies on a series of proxies to reconstruct the observations. The method provides a systematic way to verify trends; it does not require explicit modelling of the SSI, which is more challenging. The concept was initially proposed by Morrill et al. (2014). Woods et al. (2018) later used it to produce the Multiple Same-Irradiance-Level (MuSIL) technique for detecting trends. MuSIL has since been used by Meftah et al. (2020) to correct irradiance observations from the SOLAR instrument. Here, we formalize and expand the concept by bringing two improvements.

In Section 2 we describe the methodology, followed by its practical implementation in Section 3. Two applications are then provided, one to the MgII index in Section 4 and one to SSI observations from SORCE/SIM and SORCE/SOLSTICE in Section 5. Conclusions follow in Section 6.

2 Methodology

Consider a time series of solar irradiance s(t) whose trend we want to investigate. In the following s_o(t) denotes the observed irradiance that is possibly affected by an unknown trend, and s_t(t) the true but unknown irradiance. The ratio between the two(1)will henceforth be called correction and defines the relative amount by which the observations have drifted away from the true value. This correction is our main quantity of interest. We assume that the true and observed irradiance are equal when the instrument starts operating at time t₀ so that f(t₀) = 1.

To estimate the correction we make the four assumptions, whose validity will be discussed later:

1. We have an independent and stable measurement of a reference quantity r(t) that can be considered a proxy for the true irradiance s_t(t). That is, we assume the existence of a deterministic link between the two that can be described by a time-independent function s_t(t) = F(r(t)). This is equivalent to saying that the two are in phase or antiphase. The function F(r(t)) does not have to be linear.

2. The functional relationship F must be monotonic so that each value of r(t) corresponds to a unique value of s_t(t) and vice-versa. This monotonicity requirement can be partly alleviated. However, the time-invariance of F (and thus the requirement for both quantities to vary in phase) is essential.

3. Both the reference and the irradiance need to be observed for more than half a solar cycle to guarantee that some levels of solar activity are observed at least twice and are several months or years apart. Our trend detection cannot be applied to a rising (or decaying) phase only of the solar cycle.

4. The long-term stability of the reference proxy must be better than the magnitude of the trend to be determined. This assumption cannot be tested rigorously since we have no ground truth for the proxies either. Even the sunspot number, whose drifts are constrained by its positivity, is known to be affected by long-term changes in its amplitude (Clette et al., 2016). Such amplitude changes may mimic a trend.

Assumptions 2. to 4. can be fulfilled by a proper choice of the proxy and the time span. The monotonicity requirement is satisfied by most proxies except when these are too indirectly related to the irradiance. The brightness of sunspots in the visible band, for example, is known to vary in a non-monotonic way with sunspot area (e.g. Solanki et al., 2013). However, when the relationship between proxy and irradiance becomes that nonlinear, then the relevance of the proxy becomes questionable.

The most delicate assumption is the first one, which requires the reference and the irradiance to be in phase. Preminger & Walton (2007) and Dudok de Wit et al. (2018) have shown that this is not the case for time scales below approximately three solar rotations. On longer time-scales most proxies can be considered to be in phase. For that reason, in what follows, we ignore time scales shorter than 81 days and apply a lowpass filter to all quantities.

To illustrate how our method works, we consider in Figure 2 a reference r(t) and two dates t₁ and t₂ for which it reaches the same value before and after the solar minimum, i.e. r(t₁) = r(t₂). Since the reference and the observed irradiance are assumed to be in-phase, we should also have s_o(t₁) = s_o(t₂). Figure 2 shows that this is not the case, and the difference Δs_o (t₁, t₂) = s_o (t₂) − s_o (t₁) may then be the signature of an underlying trend.

Fig. 2 Illustration of the method, showing the reference r(t) and the observed solar irradiance so(t).

Morrill et al. (2014) and Woods et al. (2018) used that idea to detect trends in the SSI during the latest solar minimum. Their approach is relevant for estimating variations between specific dates t₁ and t₂. To obtain a more global picture, however, multiple pairs of dates need to be compared, and so we are left with a large ensemble of values for Δs_o (t_i, t_j). While all these values formally constrain the trend, in the same way, the shape of a function can be guessed from multiple estimates of its derivative, they do not give us immediate access to the magnitude of the trend and its uncertainty. Furthermore, the trend can be guessed at specific dates only. To overcome these limitations, we carry this approach one step further and introduce a parametric expression of the trend. We estimate the model parameters that give us access to the trend and its uncertainty at any time. If for, some reason, the trend is poorly constrained at some given time, then this will be manifested by a large uncertainty of the correction.

Let us consider a pair of dates (t₁, t₂) such that r(t₁) = r(t₂), and therefore s_t (t₁) = s_t(t₂). From equation (1) we have(2)

This ratio and the value of f(t) can be determined at specific times only. To estimate f(t) at any time we need a parametric model. For a simple linear trend the model would be(3)where T is the total duration of the record. This model can be made more flexible by adding more terms. Polynomial expansions are popular but Fourier series are numerically better behaved(4)

To ease notation, let us define a reduced time τ = (t − t₀)/T. In this new time frame observations start at τ = 0 and end at τ = 1. Our model then becomes(5)

The set of Fourier modes constitute an orthonormal basis on [0, 1], which offers considerable numerical advantages over polynomial expansions. The main free parameter is the number of terms N in the Fourier series; N determines the level of details that be described by the correction since the range of time-scales goes from T/N to T. We shall see later that N is primarily constrained by the amount of observations and by the duration T.

The procedure now consists in identifying a large number of pairs of dates {τ_i, τ_j} for which the references are identical. Next, we use equation (4) to define(6)which can be cast as a linear system(7) (8)

In matrix notation(9)where μ_ij(k) = s_o (τ_i)sin(πkτ_j) − s_o(τ_j)sin(πkτ_i). The indices (i, j) run over all pairs of observations such that r(τ_i) = r(τ_j). Under normal conditions, the number of dates M is considerably larger than the number N + 1 of unknown parameters, which allows us to solve the linear and overdetermined set of equations in equation (7) in a total least-squares sense to obtain {a₀, a₁, …, a_N}. Inserting these in equation (4) gives us the desired correction f(t). Because our model is parametric the correction can be estimated at any time, even during periods when no pairs of dates are available.

3 Practical implementation

We now provide some practical guidelines on implementing the method.

3.1 Uncertainty of the distribution

The robustness of the correction is a crucial aspect of the solution. There are two major contributions to the model uncertainty: noise or errors in the observations (and in the corrections applied by the instrument teams), which generates a natural scatter in the model parameters. The other is the mismatch between our model and the smoothed observations. Usually, only the first one is mentioned in the literature. When the data come with uncertainties, then the uncertainty on the outcome may be estimated by error propagation. The second contribution, which reflects the assumptions we make, is often overlooked. Unfortunately, it can easily exceed the first contribution.

We estimate both contributions to the model uncertainty by doing block bootstrapping (Hastie et al., 2009). Bootstrapping is a resampling method that consists of drawing a large number of pairs of dates M times randomly from the total number of M pairs available with replacement. For each draw, the correction is estimated; the dispersion of these corrections then provides an estimate of their uncertainty. In block bootstrapping, we correct for the presence of serial correlation by considering groups of consecutive dates rather than individual dates. In all the cases we encountered, the dispersion was found to follow a Gaussian distribution. For that reason, the uncertainty will be quantified in what follows by the standard deviation.

3.2 Model order

The three main tuneable parameters of the method are the cutoff frequency for removing solar transients, the number of amplitude levels to determine pairs of dates, and the number N of Fourier modes or model order. As already mentioned, transients occurring on time scales shorter than approximately 3 months can be ignored, so we set the cutoff at 81 days.

The number of amplitude levels directly impacts the number M of pairs of dates. Its choice is a compromise between a proper statistical coverage of the variability and computational expense. To a lesser degree, it also depends on the relative error because inputs with large noise amplitudes require fewer levels. We found 100–300 levels to give satisfactory results with SSI data. Larger numbers do not bring any improvement and rather lead to overfitting. In our implementation of the method, these levels are regularly spaced between the extrema of s_o (t).

The only free parameter that requires careful tuning is the model order N. Small values of N will cause important features of the correction to go unnoticed. Conversely, large values will lead to overfitting and generate unstable solutions. Eventually, the model order is limited by the shortest time-scale, which cannot be less than 81 days, i.e. N ≤ T/81. A classical method for selecting the best order is cross-validation (Hastie et al., 2009): one part of the values is used to estimate the model, which is then tested on the other part. The reconstruction error usually drops with increasing N until it reaches a minimum before it increases again because of overfitting. The Bayesian Information Criterion (BIC) is often used in this context as it seeks a compromise between minimum error and model complexity (Ljung, 1997).

Another strategy relies on the condition number (Press et al., 2002) of the regression matrix in equation (9), which quantifies the degree of collinearity between the columns of the matrix. Ideally, this number should be equal to one while values larger than one are indicative of collinearity and, therefore, error amplification. Examples will be provided in the next section.

After performing many tests, we obtained the best results by using the BIC (i.e., cross-validation) under the condition that the condition number remains below a predefined threshold. The reason for this is that with excessively large numbers M of pairs (for example, when the spacing of the amplitude levels is too small), the randomly selected values will not be independent anymore, which will lead to an underestimation of the reconstruction error. Although this effect can be corrected, in the following, we prefer to rely on the condition number because of its simplicity.

3.3 The procedure in practice

The procedure for estimating the corrections can be now be summarised as follows:

1. Remove short transients that may be caused, for example, by flares, solar rotation, sunspot deficit (in the visible range) or anomalous values. We consider a running median filter that is typically 27 days (one solar rotation) wide and is centred on the day of interest. Median filters perform better than classical running averages in the presence of sunspot darkenings in the visible band or flares in the UV band.

2. Since different bands of the solar spectrum do not evolve in phase, we lowpass filter them with an 81-day cutoff.

3. Identify pairs of dates at which the reference r(t) has the same value. We bin the filtered reference r(t) into N_lev equispaced levels, with N_lev = 100–300. For each level, we determine the times at which r crosses that particular level. With n dates, we obtain sets {τ_i, τ_j, s_o (τ_j), s_o (τ_i)}, which are stored. Occurrences that are less than 81 days apart are ignored.

4. Once the set is complete, we estimate the correction f(t) by total least squares and its confidence intervals by block bootstrapping (Hastie et al., 2009).

4 Application to the MgII index

Let us first illustrate our method with the MgII index, which is the core-to-wing ratio of the MgII K line at 280 nm. This index is widely used as a UV proxy for upper atmospheric modelling and has been measured almost continuously and on a daily basis since November 7, 1978. Here we consider two composites that have been built by stitching together observations from different sets of instruments: the composite made by the Laboratory for Atmospheric and Space Physics (hereafter named LASP composite) and the one by the International University of Bremen (Bremen composite). The two composites rely on different sets of instruments (Snow et al., 2019). Both agree until early 2011 and then start diverging. The exact reason for this is not known, although stray light in the SORCE/SOLSTICE instrument, which is the main contributor to the LASP composite after 2003, might be one cause. For recent years the Bremen composite relies instead on the MgII index from METOP/GOME-2A (starting in 2007) and from METOP/GOME-2B (starting in 2012). The time interval of interest runs from November 7, 1978, until July 14, 2013, when the LASP composite was discontinued.

In this first example, we compare two identical observables, and so formally, our method is not needed because we can just as well take the ratio of the two composites. However, this example provides us with an opportunity to better understand how the method works. Without loss of generality we take the Bremen composite as the reference proxy r(t) and the LASP composite as the observable s_o (t).

Both composites are dimensionless quantities but are scaled differently because of the way the core-to-wing ratio is estimated. To ease their visualisation, we rescale the LASP composite to the Bremen composite by means of an affine transformation. Figure 3a displays the two composites (after performing a lowpass filtering with an 81-day cutoff) and reveals a small but significant discrepancy after approximately 2010. We now decompose both records into 100 equispaced levels and estimate the correction for orders N = 1 − 50. The lowest order correction captures variations with a time-scale of 34 years, whereas with N = 50 the time-scale is 8 months only.

Fig. 3 Illustration of the MgII index correction with: a) rescaled indices from LASP and Bremen, b) the BIC for model orders up to 50, c) the condition number for the same orders, d) comparison of the correction (from a 12th order model) and the ratio of the two indices and e) comparison of the index from Bremen with the corrected index from LASP, using the correction of d). All time series are lowpass filtered at 81 days.

Figures 3b and 3c show respectively the Bayesian Information Criterion (BIC) and the condition number. The best model is the one that has the lowest BIC or, alternatively, the one for which the BIC stops decreasing significantly. Here the best order is N = 12. There is no equivalent criterion for the condition number since this depends on the desired numerical accuracy. In our context, a reasonable limit is a condition number between 100 and 500. From these two criteria, we conclude that the model should be at least of order 12. Adding more terms does not improve the fit. Let us, therefore, keep N = 12.

The estimated correction f(t) is shown in Figure 3d with its ±1σ confidence interval. The correction remains close to 1, as expected. After 2009, the correction is larger than one, and indeed the LASP composite has been found to overestimate the true MgII index (Snow et al., 2014). This is an example wherein confidence intervals are important as they enable us to assess the significance of the trend. For example, we can show that the small local increase in 2003–2005 (which is when the LASP composite started using SORCE/SOLSTICE data) is significantly different from one at a level of 5%.

For this particular example, we can estimate the correction directly, simply by taking the ratio between both composites. This would not be possible in the general case when the records represent different physical quantities. Figure 3d shows that this ratio is in excellent agreement with the estimated correction, given the confidence interval of the latter. The ratio shows small-scale variability that cannot be captured by our model because it has not been designed for that. Note that the confidence intervals tend to be narrower near solar maximum, where the correction is better constrained by a large number of pairs of dates with the same level of solar activity.

Once we have the reference, we can reconstruct the corrected (“true”) composite from the observed one by using the definition s_t(t) = s_o(t)/f(t). Without surprise, the corrected and reference composites fully overlap in Figure 3e. Making such a correction may be problematic for quantities that drop to zero at minimum, such as the sunspot number. There are two workarounds. One is to add an artificial offset and keep this in mind when using the correction. A better solution is to define the correction to be additive rather than multiplicative, i.e. s_t(t) = s_o(t) + f(t). This leads to a non-dimensionless estimator, which differs from the one of equation (5). The concept, however, remains unchanged.

5 Application to SORCE SSI observations

In this second and more realistic example, we consider sixteen years of daily SSI observations made between May 2003 and April 2019 by the SORCE/SOLSTICE and SORCE/SIM instruments. Their spectral ranges are respectively 115.5–309.5 nm and 240.0–2412 nm. Here we restrict ourselves to the UV and visible range up to 700 nm. Our time interval is constrained by the launch of SORCE in 2003 and the time span of one of the proxies. The SORCE/SOLSTICE (level 3, version 17) and SORCE/SIM (level 3, version 27) data are available at http://lasp.colorado.edu/home/sorce/data/. Because of battery degradation, after 2013, SORCE could be operated only during the sunlit portion of the orbit. This led to larger thermal amplitudes that affected the performance of its instruments. The mission ended on 25 February 2020.

A few years after SORCE was launched, both SORCE/SOLSTICE and SORCE/SIM teams reported a solar cycle modulation that was unusually strong as compared to earlier missions and SSI models. This discrepancy has given rise to a continuing debate (Harder et al., 2009; Lean, 2010; Lockwood, 2011; Deland & Cebula, 2012; Ermolli et al., 2013; Morrill et al., 2014; Marchenko et al., 2016). Indeed, if truly solar, such variations would significantly affect the radiative forcing of climate (Haigh et al., 2010; Wen et al., 2013; Shapiro et al., 2013; Bolduc et al., 2015; Ball et al., 2016a). Although successive corrections have progressively reduced these discrepancies, the existence of an undocumented instrumental trend remains unsettled.

The detection of such trends is challenged by the complex degradation correction process of the instruments and the absence of ground truth. Comparisons with independent observations were made Marchenko et al. (2016) using AURA/OMI data after July 2004 in specific bands between 264 and 504 nm, by Lockwood (2011) using UARS/SUSIM between 115 and 411 nm until August 2005, and by Mauceri et al. (2020) for TSIS/SIM after 2018. These studies confirmed an unusual variability for observations made by SORCE but also suffered from a limited overlap in time.

5.1 Solar proxy selection

In the absence of simultaneous SSI measurements with the same spectral and temporal coverage, we resort to a set of solar proxies that have been selected for their high stability and their ability to mimic solar spectral variability. Let us stress that none of these proxies, except for one, are direct SSI measurements. In addition, they track different physical properties of the Sun. Therefore, even though they are highly correlated with SSI, we do not expect them to reproduce all of its features (Dudok de Wit et al., 2009). In the following, our working hypothesis will be:

The solar proxies we consider in this study are:

• The sunspot number (SSN, version 2.0) from the Royal Observatory of Belgium (Clette et al., 2016). Although it is only loosely connected to the SSI, the sunspot number remains a popular proxy for solar activity. The data are available at https://wwwbis.sidc.be/silso/datafiles.

• The daily sunspot area (DSA) on the solar disk after (Mandal et al., 2020). This is a proxy for magnetic flux emergence, which is a driver of SSI variability (Preminger & Walton, 2007; Dudok de Wit et al., 2018). These data are available up to 2019 at http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/640/A78.

• The solar radio flux at 10.7 cm (or F10.7) from Penticton Observatory (Tapping, 2013). This flux is widely used as a proxy for Extreme UV (10–120 nm). The F10.7 index is available at https://lasp.colorado.edu/lisird/data/penticton_radio_flux/.

• The solar radio flux a 30.0 cm (or F30) from Nobeyama Radio Observatory. The flux at that wavelength receives a stronger contribution from thermal Bremsstrahlung and has been shown to be a better proxy than F10.7 for the Far-UV (120–200 nm) and Medium-UV (200–300 nm) (Dudok de Wit and Bruinsma, 2017). The F30 index is available at https://spaceweather.cls.fr/services/radioflux/.

• The Mg index from Section 4. This is a proxy for the Far and Medium UV. We consider the Bremen composite, which is the only one that covers the full-time interval. The proxy can be retrieved from https://www.iup.uni-bremen.de/UVSAT/Datasets/mgii.

• The total solar irradiance (TSI), which is not a proxy but the integral of the SSI, is measured by an independent instrument. The TSI variability is dominated by contributions from the visible band (400–700 nm). Here we consider an update of the composite by Dudok de Wit et al. (2017), which is available from https://spot.colorado.edu/koppg/TSI/. Although this composite is public, it is not (yet) the official one. Note that the close connection between TSI and SSI has already been used by several authors to correct for instrumental effects in the SSI (e.g. Mauceri et al., 2018).

While there exist other proxies, the ones we have selected stand out by their high stability. Each of them captures a different aspect of solar variability and thus should lead to a different correction. In their MuSIL method, Woods et al. (2018) combine several of these proxies together into one single “super proxy”. We recommend instead testing each individual proxy so that each correction can be documented and evaluated in light of its physical properties. With an open science approach in mind, we want these corrections to remain fully traceable.

When determining the best order from the BIC and the condition number, we find a broad minimum with 50% of the values located between 5 and 13 for SORCE/SOLSTICE and 5 and 12 for SORCE/SIM. Differences mainly arise because of the temporal evolution of the correction. Most of the corrections for SORCE/SIM can be approximated by a linear slope so that a low order model suffices. For SORCE/SOLSTICE the corrections often have a more complex shape and therefore require higher orders. To ease comparisons, in what follows, the order of all corrections will be 9. With this, the shortest accessible time-scale in the correction is 2 years.

5.2 Observed trends

Figure 4 illustrates the SSI from SORCE/SOLSTICE at a specific wavelength of 142.5 nm by showing the original and corrected spectral irradiance. Also shown are the corrections from each of the seven above-mentioned proxies. Note the scatter of the corrections, which is not surprising since each proxy quantifies a different aspect of solar variability. This scatter is also reflected by the uncertainty of the reconstructions, which is not shown. These results already highlight the importance of considering several proxies.

Fig. 4 Illustration of the SSI at 142.5 nm from SORCE/SOLSTICE and its correction, with a) the observed spectral irradiance before and after applying the correction and b) the corrections estimated from each of the proxies, using 9th order models, and their median value. The latter has been used for correcting the SSI in plot a). The uncertainties of the corrections (not shown) are comparable to or smaller than the dispersion.

Figure 4 reveals a feature that is common to most other wavelengths, with a tendency of the F10.7, F30 and MgII indices to give comparable corrections while the daily sunspot area (DSA) and the sunspot number (SSN) often lead to different corrections. This is not surprising as the latter two are of a different nature and always return to the same level at the solar minimum.

Another noteworthy result is the presence of a common initial decay for all proxies. Indeed, all corrections initially drop below 1. Such a common correction strongly supports the existence of an undocumented drop in the spectral irradiance at 142.5 nm. Figure 5 illustrates a similar analysis for SORCE/SIM at 244.56 nm. In contrast to the previous example, we have a remarkable agreement between all proxies, which provides strong evidence for an uncorrected drop of 5%.

Fig. 5 Same plot as Figure 4, but for SORCE/SIM at 244.56 nm.

Corrections at other wavelengths actually reveal a great diversity of patterns. While most share an initial drop, only a few are monotonic. In an attempt to classify these patterns, we performed a k-means classification (Hastie et al., 2009) and found that there are essentially four main patterns. The ones found for SORCE/SOLSTICE are illustrated in Figure 6. All of them show an initial drop, except for pattern no. 2. The dominant pattern for each wavelength is shown in Figure 6b.

Fig. 6 Illustration of the correction applied to SORCE/SOLSTICE using the MgII index and a 3rd order model, showing: a) the four main patterns that can be identified in the corrections (with arbitrary amplitudes), and in b) the magnitude of the correction after 10 years. The colour of the correction corresponds to the dominant pattern, using the same colour code as in plot a).

Several details corroborate the instrumental nature of the observed trend. For example, the discontinuity at 180 nm coincides with the transition between two channels of the SORCE/SOLSTICE instrument. Above 280 nm, the sensitivity of the instrument gradually decreases, probably because of an overcorrection of the degradation. In addition, Figure 6b shows that the patterns are not randomly distributed in wavelength but occur in spectral bands. Note the diversity of patterns, a few of which can be reduced to a simple exponential decay, and their marked wavelength dependence.

Let us now compare the magnitude of the correction to the solar cycle variability, see Figure 7. For each wavelength, we consider the median of all corrections and also display the range between the smallest and largest correction. We prefer the median to the mean because we are looking for a central tendency in the corrections, not their average. The downside of the median is the occasional presence of discontinuities in its slope. Such a discontinuity was visible in Figure 4 in 2019. The two instruments exhibit different wavelength dependence in the spectral range in which they overlap; this was to be expected as they are independent.

Fig. 7 Median correction estimated from SORCE/SOLSTICE and SORCE/SIM in dimensionless units (a) and after normalisation by the solar cycle amplitude (b). The shaded area represents the range between smallest and largest correction. We use a 9th order model with all six proxies and the correction shown is the value averaged over the full time interval.

The strongest corrections occur in the UV bands, with values between −5% and +5% after 15 years. In the visible band, they are considerably smaller and remain within ±0.5%. Although each wavelength has a different correction, there is a tendency for most corrections to drift during the first few years of the mission and then to stabilize. When averaged over one solar cycle, these corrections correspond to stability of approximately 0.5% per year for SORCE/SOLSTICE and 0.05% per year for SORCE/SIM. These numbers are in good agreement with the stability specifications, which are 0.5% per year for SORCE/SOLSTICE and 0.03% per year for SORCE/SIM. Small as they may be, these corrections actually are considerable when we compare them to the solar cycle amplitude, see Figure 7b. This figure shows that the largest corrections can easily be several times larger than the solar cycle modulation.

There are two lessons to be learnt from these comparisons. First, our corrections are consistent with other estimates and with instrument specifications. The good agreement obtained between corrections from different proxies (or, conversely, the small dispersion in Fig. 7b) strongly supports the idea that the degradation of instrumental sensitivity may not be fully accounted for by the current calibration.

Secondly, while the corrections may seem small in absolute terms they are actually comparable to or larger than the solar cycle magnitude. Therefore, great care should be taken in interpreting variations reported by SORCE/SOLSTICE and SORCE/SIM. This is particularly true for variations observed between 2003 and 2009 when the correction was largest and partly mimicked the shape of the decaying solar cycle.

To illustrate the impact of these corrections, we compare in Figure 8 three different versions of the SSI: 1) the original data from SORCE/SOLSTICE or SORCE/SIM, 2) the original data with our correction, and 3) the original data with the MuSIL correction (Woods et al., 2018). The latter two reconstructions agree on the magnitude of the trend but still show significant differences, especially for SORCE/SOLSTICE at 249.5 nm. Notice, however, that they are not strictly comparable because MuSIL relies on an earlier release of SORCE data. Also shown on this plot is the confidence interval that results from bootstrapping. The scatter between reconstructions that are based on different proxies provides another measure of the error, which is comparable to the former. However, as stated before, this scatter is misleading because, depending on the spectral band, not all proxies are supposed to offer a proper representation of solar variability.

Fig. 8 Comparison of different SSI reconstructions. Plot (a) compares our reconstruction from SORCE/SOLSTICE and SORCE/SIM at the wavelength nearest to 249.5 nm. For easier comparison, the record from SORCE/SIM has been shifted downwards by 0.01 mW/nm/m2. Plots (b–f) compare reconstructions based on our method and on MuSIL only. In each plot we use a 9th order model and consider the median correction of all six proxies. The reconstruction, according to MuSIL (available from https://lasp.colorado.edu/lisird/), is based on an earlier release of SORCE data and therefore is not strictly comparable to our reconstruction. Short dropouts observed at 500.1 nm are primarily due to sunspot darkening. The shaded area around our reconstruction represents a ±2σ confidence interval that has been estimated by bootstrapping.

If we apply our correction to the MuSIL dataset and compare it to the reconstruction from original data, we find a good agreement (not shown), given their uncertainty. This important result means that our method is robust with respect to degradation corrections applied by the instrument teams as long as these can be approximated by our low order model.

As a final sanity check, let us compare the reconstruction at 249.5 nm from two independent measurements from SORCE/SOLSTICE and SORCE/SIM, see Figure 8a. While the original records show substantially different trends (not shown), their reconstructions agree remarkably well apart for an offset. The main discrepancy between the two reconstructions is a weak yearly residual modulation in SORCE/SIM data, which is difficult to correct because it would require a high-order model. This example further illustrates the capacity of our method to recover a realistic solar cycle variability regardless of the trend that affects the original data. In some sense, this is not surprising because our method attempts to recover the same variability as in proxies.

5.3 Impact on radiative output

Let us finally consider how the reconstruction impacts the solar radiative output by considering the peak-to-peak variation of the SSI during the solar cycle 23. SORCE was launched two years after solar maximum, so as a substitute for the latter, we consider the first 365 days of the mission, starting on May 15, 2003. For solar minimum, we consider the year 2008. Finally, we compute the difference between the two periods and integrate over specific spectral bands. Any imbalance between different spectral bands may significantly impact stratospheric ozone levels since UV radiation below 242 nm produces ozone by photolysis, while below 320 nm it leads to a catalytic loss.

A similar approach had been used by Haigh et al. (2010) and Ermolli et al. (2013), who revealed an unusually large cycle modulation of the SSI from SORCE, with a negative variation (i.e. anticorrelation) in the visible range between 400 and 700 nm. Our results are based on more recent releases of SORCE data, which therefore are not strictly comparable to earlier plots. What matters, however, is how much the corrected SSI affects the conclusions made by Haigh et al. (2010) and Ermolli et al. (2013).

Figure 9 summarizes the solar cycle variability after and before trend correction, together with results obtained from two widely used SSI models: the Naval Research Laboratory’s SSI model, NRLSSI2 (Coddington et al., 2016) (version 02r01) and the Spectral And Total Irradiance REconstruction for the Satellite era, SATIRE-S (Yeo et al., 2014). Note a much closer agreement, both in sign and in amplitude, between the corrected SSI and the two SSI models. Although this again does not prove that the trend correction is correct, it provides additional evidence in favour of it. The estimated error on our reconstruction (not shown on the plot) is comparable to the difference between our corrected SORCE data and the two models. From this, we conclude that there is no significant difference between the solar cycle variability of the SSI derived the SORCE data – after applying our correction – and the SATIRE-S and NRLSSI2 models.

Fig. 9 Difference between the average irradiance near solar maximum (May 2003–May 2004) and at solar minimum (July 2008–July 2009) for the observed SSI, the SSI after trend subtraction, and two SSI models. Four different spectral bands are shown. The values in the Far-UV band are amplified to facilitate comparison.

6 Conclusions

We provide a mathematical framework for extracting undocumented trends from solar irradiance data by comparing the observed variability with several solar proxies such as the F10.7 index, the MgII index and the sunspot number. Compared to previous approaches by Morrill et al. (2014) and Woods et al. (2018), we provide a parametric expression of the correction, which 1) allows the correction to be estimated at any time and 2) provides confidence intervals, which are vital for testing the significance of the results.

The main assumptions are the in-phase variation of the solar irradiance with the proxies and the stability of the latter. The shape of the trend is approximated by a Fourier series, and we provide recipes for tuning the free parameters of the method. An important restriction is the necessity to observe at least half a solar cycle so that some levels of solar activity occur twice.

The application of the method to 16 years of solar spectral irradiance observations from SORCE/SOLSTICE and SORCE/SIM reveals significant trend corrections that can be as large as several times the solar cycle amplitude. Comparisons with the solar cycle variability of the NRLSSI2 and SATIRE-S models reveals a good agreement after the correction, especially in the Middle-UV band, which is important for modelling catalytic loss of stratospheric ozone. Among the various proxies that were used for this analysis, the MgII index (Bremen composite) and the F10.7 and F30 radio indices generally lead to the smallest uncertainty in the trend correction.

No comparison with solar proxies can unambiguously prove the instrumental origin of the trends we observe in SSI data. However, different sanity checks give us confidence in the method. In particular, the reconstructions are resilient (given their uncertainty) to degradation corrections applied by the instrument teams. Likewise, reconstructions of the SSI observed at the same wavelength by two different instruments also give comparable results.

The main limitation of the method is the need to specify the order of the model, which impacts the shape of the correction. Higher-order models with more detailed trend evolutions can be inferred when the spectral irradiance correlates well with the proxies. Conversely, in the visible range, where relative solar cycle variability is low, and the signal-to-noise ratio is poor, the model is less well constrained by the observations. In these cases, only a linear trend can be inferred, while in the UV, we observed a variety of patterns that do not necessarily decay monotonically. Gaussian processes (Rasmussen & Williams, 2006) may provide an interesting alternative to our parametric Fourier model, although they come with a much higher computational cost.

Finally, our method provides confidence intervals that reflect both the uncertainties in the data and the errors that result from the modelling of the trend. The latter errors are often neglected while the robustness of the model is vital for allowing the reconstruction to be tested for significance. In an open science approach, we make the Matlab routine and some test data available at https://github.com/tddwit/solar_trend_detection.

The observation of the same trend in SORCE data regardless of the type of proxy lends strong support to the presence of undocumented instrumental effects. Let us stress again, however, that our method should be regarded as an independent and complementary tool for detecting uncorrected variations. It is then up to the instrumental team to decide whether to actually correct these variations, based on the additional evidence from the instrument. Our method performs best in the UV range because there is a physical motivation for using proxies such as the MgII index and also because the relative variability of the SSI is large as compared to the intrinsic error of the method. For that same reason, the correction of the SSI in the visible and near-infrared ranges requires much greater care.

Acknowledgments

This work was supported by the French Space Agency (CNES). I am grateful to Marty Snow (LASP, now at the South African Space Agency) for making comments on an earlier version of this manuscript. I also wish to thank all the instrument teams for making their high-quality data products available: the SORCE team (LASP), Penticton Observatory, Nobeyama Radio Observatory, NOAA Space Weather Prediction Center, Solar Influences Data Centre (Brussels), Greg Kopp (LASP, TSI composite) and Mark Weber (Bremen, MgII composite). The Nobeyama Radio Polarimeters are operated by Nobeyama Radio Observatory, a branch of the National Astronomical Observatory of Japan. This research has made use of NASA’s Astrophysics Data System. The editor thanks Odele Coddington and Matthew DeLand for their assistance in evaluating this paper.

References

All Figures

Fig. 1 Relative amplitude of the solar cycle variation and the solar rotational modulation, as estimated from SORCE/SIM and SORCE/SOLSTICE measurements between May 2003 and April 2019. SOLSTICE observations are used below 300 nm, and SIM observations above. To be less affected by non-physical trends, the amplitudes are estimated by Fourier analysis and not by taking the usual difference between extrema.

In the text

	Fig. 2 Illustration of the method, showing the reference r(t) and the observed solar irradiance so(t).
In the text

Fig. 3 Illustration of the MgII index correction with: a) rescaled indices from LASP and Bremen, b) the BIC for model orders up to 50, c) the condition number for the same orders, d) comparison of the correction (from a 12th order model) and the ratio of the two indices and e) comparison of the index from Bremen with the corrected index from LASP, using the correction of d). All time series are lowpass filtered at 81 days.

In the text

Fig. 4 Illustration of the SSI at 142.5 nm from SORCE/SOLSTICE and its correction, with a) the observed spectral irradiance before and after applying the correction and b) the corrections estimated from each of the proxies, using 9th order models, and their median value. The latter has been used for correcting the SSI in plot a). The uncertainties of the corrections (not shown) are comparable to or smaller than the dispersion.

In the text

	Fig. 5 Same plot as Figure 4, but for SORCE/SIM at 244.56 nm.
In the text

	Fig. 6 Illustration of the correction applied to SORCE/SOLSTICE using the MgII index and a 3rd order model, showing: a) the four main patterns that can be identified in the corrections (with arbitrary amplitudes), and in b) the magnitude of the correction after 10 years. The colour of the correction corresponds to the dominant pattern, using the same colour code as in plot a).
In the text

	Fig. 7 Median correction estimated from SORCE/SOLSTICE and SORCE/SIM in dimensionless units (a) and after normalisation by the solar cycle amplitude (b). The shaded area represents the range between smallest and largest correction. We use a 9th order model with all six proxies and the correction shown is the value averaged over the full time interval.
In the text

Fig. 8 Comparison of different SSI reconstructions. Plot (a) compares our reconstruction from SORCE/SOLSTICE and SORCE/SIM at the wavelength nearest to 249.5 nm. For easier comparison, the record from SORCE/SIM has been shifted downwards by 0.01 mW/nm/m2. Plots (b–f) compare reconstructions based on our method and on MuSIL only. In each plot we use a 9th order model and consider the median correction of all six proxies. The reconstruction, according to MuSIL (available from https://lasp.colorado.edu/lisird/), is based on an earlier release of SORCE data and therefore is not strictly comparable to our reconstruction. Short dropouts observed at 500.1 nm are primarily due to sunspot darkening. The shaded area around our reconstruction represents a ±2σ confidence interval that has been estimated by bootstrapping.

In the text

	Fig. 9 Difference between the average irradiance near solar maximum (May 2003–May 2004) and at solar minimum (July 2008–July 2009) for the observed SSI, the SSI after trend subtraction, and two SSI models. Four different spectral bands are shown. The values in the Far-UV band are amplified to facilitate comparison.
In the text