Spectral Aging Model Applied to Meteosat First Generation Visible Band

Abstract

The Meteosat satellites have been operational since the early eighties, creating so far a continuous time period of observations of more than 30 years. In order to use this data for climate data records, a consistent calibration is necessary between the consecutive instruments. Studies have shown that the Meteosat First Generation (MFG) satellites (1982–2006) suffer from in-flight degradation which is spectral of nature and is not corrected by the official calibration of EUMETSAT. Continuing on previous published work by the same authors, this paper applies the spectral aging model to a set of clear-sky and cloudy targets, and derives the model parameters for all six MFG satellites (Meteosat-2 to -7). Several problems have been encountered, both due to the instrument and due to geophysical occurrences, and these are discussed and illustrated here in detail. The paper shows how the spectral aging model is an improvement compared to the EUMETSAT calibration method with a stability of 1%–2% for Meteosat-4 to -7, which increases up to 6% for ocean sites using the full MFG time period.

1. Introduction

The continuous and global coverage provided by satellites in space, is an important asset for climate research and monitoring. Due to the geostationary orbit around the Earth and the high temporal data frequency, climate data records from the Meteosat satellites, exploited by the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT), are extremely useful. With an image every half hour, the Meteosat Visible and Infrared Imager (MVIRI) on the Meteosat First Generation (MFG) satellites covered the full Meteosat field-of-view (FOV) at its nominal position around 0° longitude using 3 spectral bands: the visible (VIS), water vapor (WV), and the infrared (IR) band (see

Table 1. Characteristics of the three MVIRI channels.

Property	VIS	WV	IR
λmin–λmax	0.4–1.1 μm	5.7–7.1 μm	10.5–12.5 μm
Temporal frequency	25 min scanning + 5 min retracing scan mirror
Number of detectors	2	1	1
Pixel resolution	2.5 × 2.5 km	5 × 5 km	5 × 5 km
Image size	5000 × 5000 pixels	2500 × 2500 pixels	2500 × 2500 pixels

for more information on these channels).

To create climate data records, the measurements from consecutive instruments are often combined to reach a minimum time length of observations which is necessary to do climate research (25–30 years). In such data records, consistency between the different instruments is then of vital importance. The Global Climate Observing System (GCOS) has put requirements on the stability and accuracy of long-term Fundamental and Thematic Climate Data Records (FCDRs and TCDRs) derived from e.g., space observations. For the WV and IR channels of the MVIRI instruments, different vicarious calibrations have been derived in order to reach this stability [1–3]. For the VIS channel, stable calibrations have been developed using known Earth targets like desert, ocean or convective clouds [4–9] or through satellite intercalibrations [10–12].

Soon after the Meteosat instruments became operational, it was clear that the VIS channel of the imagers suffered from in-flight degradation [6,13] due to a build-up of contamination on the mirrors [14]. The SEVIRI Solar Channel Calibration (SSCC) [9] was adopted as the official EUMETSAT calibration for the MVIRI VIS channel (Even though the SSCC method was originally developed for the Spinning Enhanced Visible and Infrared Imager (SEVIRI) onboard the Meteosat Second Generation (MSG) satellites [15], it proved to work well for the MFG instruments too.). This vicarious calibration technique corrects the decrease in signal due to degradation, by allowing the calibration coefficient to increase linearly in time. However, validation showed that, first of all, there is a spectral component to the contamination-induced degradation of the VIS channel. This means that, in time, the spectral response of the instrument does not decrease uniformly over the spectrum of the channel, but that there is a wavelength dependency in the process. On top of that, the signal of the imagers does not keep decreasing at a constant rate, but saturates after a given amount of time. As neither effects were removed using the time dependent SSCC calibration, Decoster et al. [16] created a method which is able to correct the VIS Level 1.5 images of Meteosat-7 for these two effects. In this method, the authors modeled the temporal change of the pre-launch spectral response curve ϕ(λ) of the VIS channel of the satellite. Due to the effect of photodeposition [18], the spectral response of the imager decreases more strongly in the shortest VIS wavelengths than in the longer ones. For that reason, Decoster et al. [16] allowed the modeled spectral response curve ϕ(λ, t) to decrease differently at time t for different wavelengths λ in order to correct for this spectral degradation. An exponentially time dependency was added to correct for the saturation of the decreasing signal. The mathematical formula of this spectral aging model is given by

$$\varphi (\lambda ,t)=\varphi (\lambda ,0)({\text{e}}^{-\alpha t}+\beta (1-{\text{e}}^{-\alpha t}))(1+\gamma t(\lambda -{\lambda }_0))$$

(1)

where α, β and γ are the instrument dependent model parameters, ϕ(λ, 0) is the spectral response curve of the instrument as it was measured before launch, and λ₀ is the central wavelength of this pre-launch characterised spectral response. The first part of the model, e^−αt + β (1 − e^−αt), represents the wavelength independent degradation process where α is the decay rate of grey degradation and β the sensitivity of a fully degraded mirror. In this grey degradation, e^−αt represents the sensitivity of the clear non-contaminated part of the mirror, and β (1 − e^−αt) the sensitivity of the degraded (contaminated) part of the mirror. The wavelength dependency of the degradation is modeled through 1 + γt (λ − λ₀), which shows the spectral part of the ageing as a simple linear function in time and wavelength, with γ the rate of spectral degradation. For a positive γ value, this part of the model becomes less than 1 for wavelengths smaller than the central wavelength λ₀, and larger than 1 for wavelengths larger than λ₀. In Decoster et al. [16], the three model parameters α, β and γ were computed for Meteosat-7 using a set of clear-sky and cloudy targets which show a low level of variability during the time period that the satellite was operational (06/1998–07/2006) and which represent different parts of the VIS spectrum. These targets were corrected for degradation effects by converting their filtered degrading reflectance time series into unfiltered non-degrading reflectance time series. This conversion was done using a linear regression between simulations of these filtered and unfiltered reflectance observations, where the filtering of the simulations was done using the modeled time dependent spectral response curve of Equation (1). The regression coefficients that were calculated using the simulated filtered and unfiltered reflectances, were then applied to the filtered target observations. As the results and validation for Meteosat-7 looked promising, the same method is now applied to derive the model parameters for all six MFG satellites (Meteosat-2 up to -7 (Meteosat-1 failed after two years in orbit due to a design fault. As its data was never transcribed into the archive, it is not used here.)).

Section 2 briefly repeats how the stable clear-sky and cloudy targets used in this study are processed. This is in general similar to what was done in Decoster et al. [16], with a few differences which are explained here. The original degrading time series are shown in Section 3 for all six instruments. Many problems, both intrinsic to the instrument and extrinsic, like volcanic eruptions, can be observed here. These problems are addressed and corrections are made where possible. Next, the model parameters are derived, and the resulting degradation corrected time series are shown for each MFG satellite separately in Sections 4 up to 9. The best model parameters are finally summarized in Section 10, proving the importance of the spectral aging model over the SSCC calibration. In this section the full dataset is also discussed in terms of long term stability. In the future, the model and model parameters presented in this paper will be used to generate FCDRs and TCDRs for the full MFG dataset, attempting to reach the GCOS requirements for stability and accuracy as good as possible.

2. Data Selection and Processing

Between 1981 and 1997, six MFG satellites (Meteosat-2 up to -7) have been launched. The launch date and used time period of observations at a certain longitude are shown in the third and fourth columns of

Table 2. For each of the MFG satellites the following information is given: the longitude at which the instrument was operational, the launch date, the data period at that longitude that is used in this paper, the gain level, the calibration coefficient (in W·m⁻² sr⁻¹/ DC), the time averaged offset value (in DC), and the solar irradiance (in W·m⁻²).

Satellite	Longitude	Launch Date	Data Period Used	Gain Level	Calibration Coefficient	Offset	Solar Irradiance
Meteosat-2	0°	19/06/1981	12/02/1982–11/05/1987	0	0.6519	3.730	499.9
			12/05/1987–09/08/1988	1	0.5454	3.685	499.9
Meteosat-3	0°	15/06/1988	11/08/1988–27/06/1989	1	0.6277	3.712	602.2
			13/01/1990–09/12/1990	0	0.7571	4.001	602.2
(ADC)	50°W		01/08/1991–22/01/1993	0	0.7571	4.001	602.2
(XADC)	75°W		21/02/1993–22/05/1995	0	0.7571	4.001	602.2
Meteosat-4	0°	06/03/1989	19/06/1989–03/02/1994	4	0.7320	4.661	599.5
Meteosat-5	0°	02/03/1991	20/01/1994–03/02/1997	5	0.8142	4.460	690.6
(IODC)	63°E		01/07/1998–05/04/2007	5	0.8142	4.460	690.6
Meteosat-6	0°	20/11/1993	29/01/1997–13/06/1998	5	0.8376	5.542	691.4
Meteosat-7	0°	02/09/1997	03/06/1998–11/07/2006	6	0.9184	4.840	690.8
(IODC)	57°E		01/11/2006–20/09/2012	6	0.9184	4.840	690.8

.

Only one image a day is selected in this work, i.e., the one at 1200 UTC where all pixels in the Meteosat FOV are in daylight. As two detectors are used instead of one in the VIS channel (see Table 1), these images are reduced from the original size of 5000 × 5000 pixels to 2500 × 2500 pixels to average out the effect of the slightly different spectral responses of the two detectors. The original values are then converted from digital counts (DC) into reflectance ρ. In the conversion from DC to radiance L, the calibration coefficients are kept fixed at the values given on the EUMETSAT webpage ([17]), and the offset values are the mean of the different values for each instrument, as given by EUMETSAT. Both calibration and offset values are given in Table 2 for all MFG instruments, in columns 6 and 7 respectively. When converting the incoming radiation into electrical currents, either 6 (Meteosat-2 and -3) or 8 bits (Meteosat-4–7) were used, where the 6 bit values were then multiplied by 4 to get the same range of digital counts between 0 and 255 as for the more recent MVIRI instruments. To optimise this digital range, for each satellite, an electronic gain level between 0 and 15 was selected by EUMETSAT. This way, saturation could be avoided by lowering the gain, or the gain could be increased when the highest measured value at local noon was less than 212 (The gain steps have an approximate ratio of 1.2, so an increase in gain level can only be done when the maximum daily value becomes less than 212). The gain levels used for each satellite are also added to Table 2, and can also be found on the EUMETSAT webpage [17]. As changes in the gain level were made during the period Meteosat-2 and -3 were operational, the calibration and offset values had to be adapted, and so separate entries are made in the table for these time periods. In the conversion from radiance L to reflectance ρ, the Filtered Solar Irradiance (FSI) values were also calculated by EUMETSAT (see column 8 of Table 2), and the changing distance between the Sun and the Earth is taken into account. Not all pixels in the Meteosat FOV are useful. The sun glint region is removed over ocean to avoid saturation, and apart from that, some sites are too variable in time for this study. As the intention is to see how the observed Meteosat signal decreases in time due to degradation, it is necessary to reduce the noise on the time series as much as possible to be able to see the trend, and so only a selection of pixels is used.

2.1. Selection of Targets

This selection process is similar to the one explained in Decoster et al. [16], and is based on the Meteosat-7 data alone (06/1998–07/2006) as this turned out to lead to the most stable targets. The clear-sky targets are selected from clear-sky images, created using the method of Ipe et al. [19] with a temporal frequency of 10 days. To reduce the spatial variability, the clear-sky pixels are replaced by the mean values in boxes of 25 × 25 pixels around it. The pixels with the lowest variability in time are the local minima in boxes of 51 × 51 pixels in the image that shows for each pixel the ratio of the variance to the mean value, calculated from linear fits through the pixel time series. Figure 1a shows the positions of the 298 clear-sky targets, while Figure 1b indicates the scene type class into which each target is put: ocean, dark vegetation, bright vegetation, dark desert, bright desert, and snow.

Contrary to Decoster et al. [16], each local minimum is selected for its own scene type only, instead of using all land scene types together. Because of this, it can happen that two targets are located close to each other in the Meteosat FOV, but have different scene types. Like this, targets from more different regions in the FOV are selected as they cannot be eliminated due to slightly less variable sites from other scene types in their neighborhood. Apart from the clear-sky targets, convective cloudy pixels are selected in the InterTropical Convergence Zone (ITCZ) of the original images. For this, all pixels are replaced by the mean values in boxes of 7 × 7 pixels around it. The convective cloudy targets are selected as the local maxima in boxes of 151 × 151 pixels.

2.2. Reflectance Ratio r Time Series

For each satellite, the targets are converted from reflectance ρ to reflectance ratio r. This process first consists of unfiltering the reflectances ρ using simulated spectral radiances L(λ) for different scene types. These spectral radiances are first converted into filtered and unfiltered radiance values, L and L_u respectively

$$L={\int }_{\text{VIS}}L(\lambda )\varphi (\lambda ,t)\text{d}\lambda$$

(2)

$$L_{\text{u}}={\int }_{0.25-5\mu \text{m}}L(\lambda )\text{d}\lambda$$

(3)

where the filtered radiance is integrated over the spectrum of the VIS spectral response curve, and the unfiltered radiance is integrated over the full solar reflectance spectrum. Next, they are converted into filtered and unfiltered reflectances ρ and ρ_u using the same FSI as for the observations, but a fixed Sun-Earth distance d =1 as the simulations were made for d =1AU. A linear relation exists between these two reflectances

$${\rho }_{\text{u}}=a+b\rho$$

(4)

which is used to unfilter the observed reflectances of the same scene types. This unfiltering process is necessary for the spectral aging correction because of Equation (2), through which the modeled spectral response curve ϕ(λ, t) is applied. In order to decrease the variability of the time series as much as possible, the unfiltered reflectance values are divided by an anisotropy factor R and albedo A to correct for the difference in Sun-Earth-satellite geometry and albedo for the different surface types

$$r=\frac{{\rho }_{\text{u}}}{R({\theta }_0,\theta ,\psi )\hspace{0.17em}\hspace{0.17em}A({\theta }_0)}$$

(5)

where θ₀ is the solar zenith angle, θ the viewing zenith angle and ψ the relative azimuth angle. For the clear-sky sites, R and A are provided by the Angular Dependency Models (ADMs) that have been empirically estimated by Loeb et al. [20] using data of the Clouds and Earth’s Radiant Energy System (CERES) onboard the Tropical Rainfall Measuring Mission (TRMM) [21]. The validation and accuracy measurement of these models was done by Loeb et al. [22]. The availability of these empirically derived anisotropy factors and albedo values is the reason why the integration of Equation (3) is done over a broader spectral range than Equation (2). As the CERES short wave band is broader than the VIS MVIRI channel, the cloudy R and A contain the effects of the deep ice and water absorption bands between 1 and 2 μm. For this reason the anisotropy and albedo models for deep convective clouds are not taken from Loeb et al. [20]. Instead they are modelled with the radiative transfer program libRadtran [23] using the Key et al. [24] ice crystal parameterization, assuming solid column shaped particles.

For each of the 298 clear-sky targets and the 6 cloudy sites per image, reflectance ratio time series are created. This is done for the full time period of each satellite at the nominal position of 0° longitude, where the exact time period was given in Table 2. Because the correction for the surface reflectance anisotropy is fixed in time, there are is still a seasonal cycle in the data. To decrease the variability of the time series as much as possible, these residual seasonal effects are corrected by subtracting the difference between the reflectance ratio value and the monthly mean annual cycle from the time series (i.e., deseasonalization), for which a minimum of two years of data is needed. As these time series are degrading in time, the monthly mean annual cycle is calculated with respect to a linear first order fit. Next, the clear-sky time series are averaged out according to 5 scene types: ocean, dark vegetation, bright vegetation, dark desert and bright desert (As there are barely any all year round snowy pixels in the Meteosat FOV, this scene type is not used here). The second column of

Table 3. For each surface type, the number of sites used in each time series is given together with the weight values used in the cost function of Equation (6).

Surface Type	Number of Sites	Weight wi
clouds	60	0.6562
ocean	55	0.1611
dark vegetation	57	0.0252
bright vegetation	102	0.0554
dark desert	47	0.0268
bright desert	37	0.0753

shows the number of sites used to create each clear-sky time series.

The six convective cloudy pixels are averaged out every 10 days to get the same temporal frequency as the clear-sky time series. This then finally leads to 6 time series with different spectral characteristics.

In the most perfect situation, the r values should be equal to 1 for all 6 time series. In practice this is, for different reasons, not always the case. First of all, the Meteosat pre-launch characterized spectral response curve on which the spectral aging model is based, is not perfect, with possibly larger errors for the older instruments when the characterization was done less accurately than nowadays [25]. Secondly, the unfiltering correction relies on simulations which may not exactly represent the observed surface type. The final reason is that the CERES TRMM ADMs used to convert the unfiltered clear-sky reflectance values ρ_u into reflectance ratio r, are not perfect for this use here. As they are global tropical models, they are adequate as an average over the tropical region (latitude between 35°S and N), but might slightly misrepresent targets in the Meteosat FOV.

2.3. Selecting Model Parameters

The 6 time series are now used for each satellite independently to fit the model parameters α, β and γ in Equation (1). These are found by minimising the mean square error of the 6 time series with respect to the mean reflectance ratio u_i

$$C=\sum _{i=1}^6{w}_i\left(\frac{1}{N}\sum _{j=1}^N{(r_{ij}-u_i)}^2\right)$$

(6)

using the Powell technique [26], where index i runs over the 6 scene types, index j over all N points in the time series, w_i are the weights given to the 6 scene types (see Table 3) based on their presence in the FOV, and r_ij is the reflectance ratio for time series i and time j. The Powell routine requires for each parameter an initial value and an initial step size. Through a bi-directional search for each of the parameters, starting at the initial value in the direction normal to the parameter axis using the initial step size, the algorithm works toward the parameters which minimize the cost function (Equation (6)). Another difference from Decoster et al. [16], is that in this minimization process, instead of finding the minimum of the cost function for the model parameters (α, β, γ), in practice α is replaced by the slope s, which is defined as

$$s=-\alpha (\beta -1).$$

(7)

This means that the Powell minimization looks for the best model parameters (s, β, γ), and that α is derived afterwards from s and β. The reason for this modification is that, for short time series, no saturation has been reached yet and so it is difficult to know if either α and β are both big or both small. The slope, however, will always stay the same, and becomes a more reliable, and numerically stable and independent, variable in that case. The full minimization routine is summarized in the following roadmap:

(1)

Simulate the spectral radiance L(λ) for different scene types, cloudiness types and geometries

(2)

Set the model parameters (s, β, γ) to an initial value

(3)

Calculate L and L_u using Equations (2) and (3) with the given values for s, β and γ

(4)

Convert these simulated radiances into reflectances

(5)

Do the unfiltering through Equation (4), fitting the a and b values for these simulated reflectance values

(6)

Use this fit to convert the observed reflectances ρ into unfiltered reflectances ρ_u (Equation (4))

(7)

Transform ρ_u to reflectance ratio r using Equation (5)

(8)

Calculate the cost function (Equation (6))

(9)

If the variance is not yet the lowest possible, the Powell routine returns a new set of (s, β, γ) parameters and goes back to step 3.

To estimate the accuracy of the four model parameters (s, α, β, and γ), the standard deviation on each value is computed by running the same minimization technique on 30 different random subsets of 100 time series (instead of 299). As there is only one convective cloud time series, no subset can be taken here. Instead, three different time series are created by taking once the 6 highest cloudy pixels to average out per day (as usual), once the 5 highest, and once the 4 highest. These three time series are randomly added to the 30 subsets in such a way that each subset contains 1 convective cloud and 99 clear-sky time series.

3. Original Degrading Time Series

For each of the MFG satellites, six reflectance ratio time series are created as shown in Figure 2. The label in the figure indicates for which satellite the data is shown. The time series are clearly all decreasing in time due to the degradation. Of the six scene types used here, the ocean is the one that reflects most radiation in the shortest visible wavelengths. This means that the spectral effect in the degradation should be best visible by comparing the slope of the ocean time series with the slope of the other time series. For Meteosat-5 and -7 it can be seen in Figure 2 that the ocean time series indeed decreases more than for example the vegetation time series. For the other satellites, this difference is less clear due to several reasons. The time periods of observation of Meteosat-3 and -6 at 0° longitude were very short, making it difficult to do the deseasonalization. As said, at least two years of data is necessary, which was not the case for Meteosat-6 (explaining the high noise level in its time series), and was barely the case for Meteosat-3, which was split up in two periods of about one year with a break of several months in between (see Table 2). For Meteosat-2 and -4, the time series, and especially the ocean time series, were influenced by the eruptions of two volcanoes: El Chichón in Mexico (28 March–4 April 1982) and Mount Pinatubo in the Philippines (June 1991).

In the following, some of these problems, which involve several satellites, are discussed. This is, first of all, the correction done on the ocean time series to eliminate the majority of the effect of the volcanic eruptions. A second problem is caused by the 6-bit digitization that was used for Meteosat-2 and -3, and a last problem is the saturation that occurs for the same two instruments.

3.1. Aerosol Correction

Figure 3 shows aerosol optical depth (AOD) measurements of both tropospheric and stratospheric aerosols, averaged out over the 55 clear-sky ocean targets for the time period of August 1981 until June 2006. The AOD values are retrieved from the Global Aerosol Climatology Project (GACP) [27,28]of the National Aeronautics and Space Administration (NASA) Goddard Institute of Space Studies (GISS).

This dataset was created using channel-1 and -2 AVHRR data, supplemented by data from other satellites, field observations, and chemical-transport modeling. As the time period covers the full MFG dataset and the measurements were made globally over ocean, these 1° × 1° monthly mean AOD values are used to correct the MFG ocean dataset in the following way.

First of all, the AOD dataset is smoothened out in time to eliminate extreme noise in the original data. Every AOD value is replaced with the median of three values: the data point itself, the value of the same site the month before and the value of the month after. This monthly mean AOD dataset is then coupled for each of the clear-sky ocean targets, and for each of the six MVIRI instruments, to monthly mean reflectance ratio time series of these targets. As for AOD values smaller than 1 the reflectance over ocean increases linearly with the presence of aerosols [29], the same relation is valid between AOD and reflectance ratio r. Following this, the intercept at AOD =0 of the linear regression of reflectance ratio and AOD should then represent the reflectance ratio that ocean target would have when no aerosols were present. To find this linear relation for each target, the data of all satellites is combined. This way, the range in AOD values is increased, as only for Meteosat-2 and -4 high volcanic AOD values were measured and for the other satellites the amount of aerosols measured was low. The linear relation between r and AOD is calculated through the least-squares fitting of

$$r=a_i+b_{i}t+c_i{t}^2+\text{AOD}\frac{\partial r}{\partial (\text{AOD})}.$$

(8)

In this equation, there are ten parameters that need to be fitted: the satellite dependent coefficients a_i, b_i and c_i of the second degree polynomial used to approximate the degradation for each satellite (i), and the satellite independent slope of the linear relation between r and AOD, expressed by ∂r/∂(AOD). The most important parameter here is the slope ∂r/∂(AOD). This parameter is found by fitting Equation (8) on the observed reflectance ratio r and AOD values for each site, but using the data of all six satellites together. On average, over the different sites, it is equal to 0.61 ± 0.19. Knowing the value of this parameter for each target, allows to subtract the product of ∂r/∂(AOD) and AOD from each reflectance ratio ocean time series, which is equal to the intercept of the linear regression. This is done for the full ocean time series of each of the 6 satellites. In Section 3.4 it is shown that the large variability on ∂r/∂(AOD) makes it really difficult to remove all aerosol effects in the time series.

As there is no land AOD data for this time period, the 4 land time series cannot be corrected. The following short study has proven that the lack of correction over land does not pose problems in finding the model parameters. From the GACP AOD dataset, it can be shown that the root mean square AOD value is on average equal to 0.051. With the help of simulations, showing how reflectance changes with increasing AOD values for the different scene types used, the root mean square reflectance ratio can be calculated, due to a change in AOD. This results in the following values: RMS(r_ocean) = 4.8 × 10⁻², RMS(r_{dark vegetation}) = 4.0 × 10⁻³, RMS(r_{bright vegetation}) = 4.1 × 10⁻⁴, RMS(r_{dark desert}) = 7.3 × 10⁻⁴ and RMS(r_{bright desert}) = 5.7 × 10⁻⁶, which shows that only for ocean there is a significant contribution of aerosols.

3.2. 6-bit Digitisation

Meteosat-2 and -3 were part of the pre-operational phase. As said before, the discretization that was used for these satellites, when translating the electrical current into digital counts, was different for the VIS channel than for the IR channel: the output of the VIS channel was converted using 6 bits (values 0–63), while for the IR channel, 8 bits were used (values 0–255). To make the data of the two channels more comparable, the VIS data was multiplied by 4 [4]. This leads to the same range of values as for the other channel, but with steps of 4 DC. Figure 4 is made using the 55 ocean targets for the Meteosat-2 and -3 time period. The light blue crosses show the daily minimum value of the targets (instead of the mean), while the dark blue full curve represents the ocean reflectance ratio time series of Figure 2. The crosses nicely show the discrete levels, with the majority being a multiple of 4, and some intermediate values which are introduced in the rectification process when the Level 1.0 data is converted into Level 1.5. As there is an offset of about 4, the lowest non-zero ocean values for Meteosat-2 are around 8, where the volcanic eruption is visible in the time span where the values jump from 8 to 12. The gain level change of 0 to 1 for Meteosat-2 is indicated by the vertical line, separating Meteosat-2(a) from (b). It can be seen that a number of pixels jumps back up from 8 to 12 after this gain change. For Meteosat-3, the values start at 12 and 16 (Meteosat-3(a)), and fall back to 8 and 12 after the gain level went from 1 to 0 (Meteosat-3(b)). This digitization problem is worst for ocean data as ocean reflects only little light, and so almost all dark ocean values arrive at a digital count value of 8 or 12. For the other scene types, the relative difference is smaller as the signal is higher. Converting this ocean data from digital count to reflectance ratio leads to jumps in the reflectance ratio values. There is a very good correspondence between the 8 to 12 jumps in DC during the volcanic eruption of El Chichón, and the jumps in the reflectance ratio curve for that same time period. This leads to more variance on the reflectance ratio time series and probably also to a bigger effect of the El Chichón eruption in 1982 and the following years. The jumps in DC during the gain changes are not visible in the reflectance ratio time series, as this is corrected for by the different calibration coefficients used.

There is no way to solve this discretization problem yet at this point. One way to have an idea on the true value of a certain ocean pixel, can be to look at its diurnal cycle, and see how, and if, the values change. There is also strong doubt about the offset value, which might be too high due to the conversion of 6 to 8 bits, but this discussion will come back later on in Section 10, when all the results are shown.

3.3. Saturation

A third issue that is discussed here, is the saturation that took place after the gain level change of 0 to 1 for Meteosat-2 and before the gain level change of 1 back to 0 for Meteosat-3 (see Table 2 for the exact dates). Figure 5 shows the 60 convective cloudy targets at top of the image and the 37 bright desert targets at the bottom expressed in digital counts for both Meteosat-2 and -3. The upper part of the figure shows in dark blue the 6 values that were selected in each image as a convective cloudy site, while in light blue the daily mean is given of these 6 values on top of it. The bright desert time series in the figure (red curve), is made up of the daily mean values of all 37 time series, with a temporal frequency of 10 days. The jumps show the date when the gain level was changed. Normally these jumps disappear when converting the DC into radiance as different calibration and offset values need to be used before and after the gain change. Due to the saturation, however, the calibration might overcorrect for the convective clouds data as the relative jump is smaller for the saturated pixels than for the non-saturated ones.

This relative jump was calculated for both Meteosat-2 and -3. The bright desert time series in Figure 5 are used as a reference for the saturated convective cloudy time series as the SSCC calibration was based on stable bright desert scenes [15]. For the convective clouds, linear fits are made through the daily averaged time series (light blue). For the desert time series only the tops are used for the fit. The reason for this is that the desert signal is much lower than the signal for the convective clouds, and so it needs to be measured more accurately. As the tops are very stable in shape, they are more reliable to use than the full time series in the small time periods. For Meteosat-2, the relative differences between the time series before and after the gain are

$$\text{Clouds}:\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{230-200}{200}=0.150$$

(9)

$$\text{Desert}:\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{86.7-73.5}{73.5}=0.180$$

(10)

resulting in a difference between both jumps of 0.03, indicating a 3% loss of signal due to saturation, and for Meteosat-3

$$\text{Clouds}:\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{205-240}{240}=-0.146$$

(11)

$$\text{Desert}:\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{79-95.6}{79}=-0.174$$

(12)

with a difference in jump of 0.028, leading to a loss of signal due to saturation of almost 3%. Due to the fact that the convective cloudy targets are averaged out every 10 days, these effects might be less visible in the reflectance ratio time series.

3.4. Aerosol Corrected Time Series

With these issues in mind, the original time series are shown in Figure 6 after seasonal correction, where the ocean time series is corrected using the GACP AOD dataset. The reason why the peaks in the ocean time series of Meteosat-2 and -4 have not disappeared, is because it is extremely difficult to find the right value for the reflectance ratio to AOD slope ∂r/∂(AOD) as there is a lot of variance in the GACP data, especially for the lowest AOD values. Small changes to the ∂r/∂(AOD) value result in big differences in the aerosol corrected ocean time series. Apart from that, the reason that the peaks for Meteosat-2 are still higher than for Meteosat-4, is most likely due to the fact that the Meteosat-2 data was digitized using 6 bits while 8 bits were used for Meteosat-4. Even though the aerosol correction is not perfect, it does decrease the time period in which the volcanic eruptions affect the ocean time series, which improves the minimization process in finding the model parameters. For all six satellites, an overall decrease of the ocean time series is visible when comparing to Figure 2, without large temporal changes. To have an idea on the amount of degradation per year for each satellite, linear fits are made through the time series. The ratio of slope (yr⁻¹) to intercept then gives the percentage of decrease per year. These values are shown in

Table 4. Table giving the amount in % per year that the slopes of fits through the time series change for each MFG satellite before and after aging correction, together with the standard error on these values.

Surface Type	Meteosat-2		Meteosat-3
	Before	After	Before	After
convective clouds	−0.9411 ± 0.0489	−0.0678 ± 0.0477	−0.8336 ± 0.3232	0.0075 ± 0.3264
ocean	−0.3510 ± 0.1219	0.1316 ± 0.1222	−0.7546 ± 0.3728	−0.1708 ± 0.3769
dark vegetation	−0.7939 ± 0.0659	−0.0262 ± 0.0509	−2.1408 ± 0.2955	−1.4830 ± 0.2999
bright vegetation	−0.5011 ± 0.0426	0.2081 ± 0.0436	−1.7529 ± 0.2232	−1.1043 ± 0.2293
dark desert	−0.3485 ± 0.0474	0.4394 ± 0.0520	−1.6432 ± 0.2326	−0.9330 ± 0.2402
bright desert	−0.5794 ± 0.0437	0.2440 ± 0.0464	−1.2364 ± 0.1558	−0.4679 ± 0.1625
Surface Type	Meteosat-4		Meteosat-5
	Before	After	Before	After
convective clouds	−2.1375 ± 0.0652	−0.1832 ± 0.0649	−0.8667 ± 0.1033	0.1179 ± 0.1057
ocean	−1.8306 ± 0.1168	0.0407 ± 0.1153	−1.4873 ± 0.1062	−0.6385 ± 0.1413
dark vegetation	−1.1990 ± 0.0880	0.3420 ± 0.0958	−1.1213 ± 0.1158	−0.4815 ± 0.1169
bright vegetation	−1.3130 ± 0.0733	0.2756 ± 0.0812	−1.0542 ± 0.0922	−0.3321 ± 0.0926
dark desert	−1.5212 ± 0.0881	0.2430 ± 0.0939	−0.6019 ± 0.0894	0.2502 ± 0.0903
bright desert	−1.7410 ± 0.0822	0.1453 ± 0.0868	−0.7289 ± 0.0708	0.2000 ± 0.0717
Surface Type	Meteosat-6		Meteosat-7
	Before	After	Before	After
convective clouds	−0.9775 ± 0.6855	1.0503 ± 0.6779	−1.8838 ± 0.0293	−0.1016 ± 0.0261
ocean	−2.2195 ± 0.7879	−0.3318 ± 0.7890	−2.0222 ± 0.0297	−0.3107 ± 0.0268
dark vegetation	−1.5746 ± 0.6235	−0.3872 ± 0.6344	−1.4012 ± 0.0221	−0.1232 ± 0.0201
bright vegetation	−1.5382 ± 0.4842	−0.1884 ± 0.5002	−1.4073 ± 0.0215	−0.0281 ± 0.0201
dark desert	−1.3293 ± 0.5439	0.2368 ± 0.5553	−1.5427 ± 0.0226	0.01731 ± 0.0236
bright desert	−1.0094 ± 0.2533	0.7208 ± 0.2620	−1.6924 ± 0.0195	−0.0026 ± 0.0209

for each satellite, where the standard errors on these percentages are calculated through error propagation of the slope and intercept standard errors. These values are discussed in Section 10.

In the following sections the time series will be corrected for degradation effects. The results are shown for each satellite individually, comparing the time series before and after the aging correction, and giving the best values of the model parameters.

4. Meteosat-2

In June 1981, Meteosat-2 was launched as the second satellite in the Meteosat pre-operational phase. The data used at 0° longitude runs from February 1982 to August 1988. As explained, a gain change was made, resulting in different calibration coefficients before and after (Another gain change was performed earlier on 20 October 1981 going from level 1 to level 0, but this data was not used here). Even though the VIS channel makes use of 2 detectors, for the satellites in the pre-operational phase, one of the two was turned off every second image for the images at the hour (i.e., at 09:00 UTC, 10:00 UTC, etc.) because of transmission bandwidth limitations [30]. The 6 original reflectance ratio time series for Meteosat-2 are shown in Figure 7 before (a) and after (b) the seasonal correction, and (c) after the ocean time series are corrected for aerosols.

There seems to be some yearly increased noise in the second half of the year in the clear-sky time series of Figure 7b,c. Comparing with Figure 7a, this falls together with the strong dips of the residual seasonal cycle. The reason for this added variance comes from the fact that the seasonal correction is computed monthly and, in this case, the variation per month is big and does not correct all of the 3–4 images per month used. The effect of El Chichón is most clearly visible in the reflectance ratio time series of the ocean data in Figure 7c, but also for the 4 land time series a temporal increase in r can be seen in the time period of April 1982–June 1984. As aerosols both reflect and absorb solar radiation, the lack of decrease in the time series confirms the results from Kondratyev and Galindo [31] that only a negligible amount of absorption occurred for the type of aerosols that were present after the El Chichón eruption. No effect of the eruption is visible in the convective cloud time series, as the reflecting tops have such high albedos that neither reflecting or absorbing effects of the aerosols are visible. From the figures it is also clear that the aerosols spread quickly over the full equatorial range of the Earth (about a month), but that it took years before the effect of the aerosols disappeared from the data.

The 6 time series for the six different scene types are corrected for aging by minimizing the cost function of Equation (6). Here, the parts of the time series which are most affected by the volcanic eruption are not used as they would lead to overcorrections. The part that is avoided for the five clear-sky time series is the time period between 4 June 1982 (This is a few months after the eruption took place, because this seems to be the moment around which the aerosols were spread far enough around the globe to start affecting all ocean data in the clear-sky images.) and 2 October 1984, which is a balance between taking away enough data and still have enough to fit the parameters to. The same period is avoided when doing the seasonal correction. Due to the high level of variance on the ocean time series, it is not possible to find a positive γ parameter using this minimization technique. Figure 7d shows the six aging corrected time series, where the model parameters used are

$$\text{slope}=-0.017\pm 0.008\hspace{0.17em}{\text{yr}}^{-1}$$

(13)

$$\alpha =0.00044\pm 0.00034\hspace{0.17em}{\text{day}}^{-1}$$

(14)

$$\beta =0.90\pm 0.26$$

(15)

$$\gamma =0.00\pm 0.00\hspace{0.17em}\mu {\text{m}}^{-1}{\text{day}}^{-1}$$

(16)

The percentages of change per year of the linear fits through the time series before and after ageing correction are shown in Table 4, together with the standard deviations on these values. In the calculation of these slopes, not the full time series are used. To avoid having the volcanic eruption affecting the values of the slopes too much, the time periods where the effect of El Chichón is biggest is not used, i.e., 4 June 1982–2 October 1984.

5. Meteosat-3

The third Meteosat satellite—and at the same time the last of the pre-operational programme—was launched in June 1988. The satellite replaced Meteosat-2 as the main operational satellite in August of the same year and stayed at 0° longitude for about 10 months. After a period of 7 months, it started to collect data again for a second period of a bit more than a year. For this second period the gain level was changed back from 1 to 0. Figure 8 shows the original time series for Meteosat-3 before (a) and after (b) seasonal correction, and (c) after aerosol correction for the ocean time series. More variation is visible in the time series of the second period. This could be due to the deseasonalisation, which is done using both parts together but is in total barely enough to calculate the mean annual cycle. Apart from that, there is also the 6-bit digitisation problem which introduces jumps in the reflectance ratio time series, especially in the second part (as could also be seen from Figure 4).

The ageing model is applied to the time series, with the result shown in Figure 8d where the parameters used are

$$\text{slope}=-0.009\pm 0.027\hspace{0.17em}{\text{yr}}^{-1}$$

(17)

$$\alpha =0.00010\pm 0.00029\hspace{0.17em}{\text{day}}^{-1}$$

(18)

$$\beta =0.75\hspace{0.17em}(\text{fixed})$$

(19)

$$\gamma =0.0000\pm 0.0004\hspace{0.17em}\mu {\text{m}}^{-1}{\text{day}}^{-1}$$

(20)

It was not possible to find significant γ and β values, so the value for γ is equal to the smallest positive number, and, as the time period is really short, the β parameter was fixed at 0.75, the value found in Decoster et al. [16] for Meteosat-7. The residual slopes of the 6 time series, expressed in percentage per year, are shown in Table 4 before and after ageing correction.

5.1. Atlantic Ocean Data Coverage (ADC/XADC)

In August 1991, Meteosat-3 was moved over the Atlantic Ocean for the Atlantic Ocean Data Coverage (ADC) at 50° W, in loan to the National Oceanic and Atmospheric Administration (NOAA) to take over the operations from the Geostationary Operational Environmental Satellite GOES-6 which had failed in 1989. After one year and a half, Meteosat-3 was relocated to the GOES-EAST nominal position of 75° W, to run the Extended-ADC (XADC) from February 1993 until November 1995. The FOVs for the ADC and XADC are shown in Figure 9a,b respectively.

In an attempt to improve the model parameters of Meteosat- 3 by increasing the length of the Meteosat-3 time series, reflectance ratio time series are created for the ADC and XADC datasets for the same 6 scene types. As there is barely any overlap between the FOV above the Atlantic Ocean and the FOV at the nominal position of 0° longitude, new sites are selected. This is done in the same way as explained above, but looking for overlapping sites which are present in both FOVs. The selected 134 sites are indicated by the white boxes in Figure 9a,b. Figure 10 shows the resulting time series after aerosol correction for the nominal position, ADC and XADC data. The time period for the ADC is too short to perform the seasonal correction. The correction factors of the XADC period were used here, but they are not good enough to correct the ADC data for residual seasonal effects. On top of that, the time period of the ADC was also right in the aftermath of the Pinatubo eruption of June 1991, of which the results and effects are further explained in the following section. It was not possible to take the effects of Pinatubo out of the data, due to a combination of the fact that the seasonal correction was not good enough, the Meteosat-3 data was processed using 6 bits, and the fit of reflectance ratio r with respect to AOD was too unstable (At the nominal position, it was possible to use all 6 satellites to do this fit, and even then the method was very sensitive to small AOD changes. In this case, only four years of data are available to fit the r to AOD relation on.). The latter is also visible in the XADC period, where the AOD correction introduced extra variability in the ocean time series. The fact that the convective cloudy time series do not couple over the three time periods of Meteosat-3 data, hints to a possible other gain change when Meteosat-3 was moved over the Atlantic Ocean. An attempt was made to correct these time series for spectral aging, but due to the fact that the ADC time period is not useable, and the ocean time series of the XADC is not useful either, this does not lead to an improvement of the Meteosat-3 model parameters.

6. Meteosat-4

The Meteosat Operational Programme (MOP) started with the launch of Meteosat-4 in March 1989. The operational period used here runs from June 1989 until February 1994. The calibration coefficient and offset value used for the Meteosat-4 time period were shown in Table 2. For Meteosat-4 and the rest of the MFG satellites, no more gain changes were performed while the satellites were operational. The original reflectance ratio time series for Meteosat-4 are shown in Figure 11a before and (b) after seasonal correction, and in (c) after aerosol correction for the ocean time series. The effect of Pinatubo is visible as a steep increase in Figure 11b at the moment of eruption, and a gradual decrease afterwards as aerosols leave the stratosphere. The reflectance ratio increase due to the eruption is not entirely compensated for after the aerosol correction, but has decreased more than for Meteosat-2. The two desert time series also seem to be influenced by the presence of volcanic aerosols: a dip can be seen at the same time as the ocean time series starts to increase. This indicates that the aerosols that were present in the atmosphere after the Pinatubo eruption, had stronger absorbing properties than for El Chichón, reducing the amount of reflectance over the brightest land time series. For the convective clouds, the effect is negligible as the reflectance of convective clouds is still almost twice that of desert land sites. The added variance that was visible in the time series of Meteosat-2 due to the deseasonalization, can be seen here too. When comparing with Figure 11a, it is clear that it coincides with the steep slopes in the seasonal cycles.

In Figure 11d, the time series have been corrected for spectral aging. Similarly as for El Chichón in the Meteosat-2 data, the period where the effect of Pinatubo is the largest was not used in the seasonal, nor in the aging correction of the Meteosat-4 time series. The ocean data was not considered between 21 June 1991 and 8 July 1993, while the land time series were not used between 21 June 1991 and 26 December 1991. This leads to the following model parameters:

$$\text{slope}=-0.026\pm 0.002\hspace{0.17em}{\text{yr}}^{-1}$$

(21)

$$\alpha =0.000276\pm 0.000075\hspace{0.17em}{\text{day}}^{-1}$$

(22)

$$\beta =0.743\pm 0.145$$

(23)

$$\gamma =0.000049\pm 0.000037\hspace{0.17em}\mu {\text{m}}^{-1}{\text{day}}^{-1}$$

(24)

Table 4 shows how the fits through both the original aerosol corrected and aging corrected data evolve in time. These values were calculated avoiding the same time periods for the ocean and land time series where the effect of Pinatubo is visible in the data.

7. Meteosat-5

As the second MOP satellite, EUMETSAT launched Meteosat-5 in March 1991. The satellite became officially operational in May 1991, but until January 1994 Meteosat-4 and -5 alternated at the nominal position as operational satellite. There is not enough Meteosat-5 data available to use in this period to create the clear-sky images, so the dataset used runs only from January 1994 until February 1997. Following the suggestion of Govaerts [32], the spectral response curve of Meteosat-7 was used in the unfiltering process of the Meteosat-5 data, for which the calibration coefficient and offset value had been adjusted on the EUMETSAT webpage, and are used here in this work (see Table 2). In Figure 12a,b, the 6 original time series are shown for Meteosat-5, before and after seasonal correction respectively. There seems to be again some yearly increased noise around the months September–October in the time series of the (b) part of the figure. Comparing with Figure 12a, this falls together again with the deeper dips of the seasonal effects. Figure 12c shows the ocean time series which have been corrected for aerosols.

After correcting the time series with the aging model, the minimization process leads to the parameters

$$\text{slope}=-0.0060\pm 0.0073\hspace{0.17em}{\text{yr}}^{-1}$$

(25)

$$\alpha =0.000066\pm 0.000081\hspace{0.17em}{\text{day}}^{-1}$$

(26)

$$\beta =0.75\hspace{0.17em}(\text{fixed})$$

(27)

$$\gamma =0.00015\pm 0.00011\hspace{0.17em}\mu {\text{m}}^{-1}{\text{day}}^{-1}$$

(28)

for which the corrected time series are shown in Figure 12d. The β value was kept fixed again, as the Powell method did not lead to a plausible value.

7.1. Indian Ocean Data Coverage (IODC)

Meteosat-5, -6 and -7 were all three relocated over the Indian Ocean after their operational period at 0° longitude, to support the Indian Ocean Experiment (INDOEX). Meteosat-5 was moved to 63° E, and was the operational satellite for the IODC from June 1998 until December 2006. The data over the Indian Ocean is treated in the same way as the data at the nominal position, in order to increase the time period and improve the model parameters. Different targets need to be found for the IODC dataset, as the IODC FOV barely overlaps with the FOV at 0° longitude. Figure 13 shows the Meteosat FOV over the Indian Ocean at 63° E with the selected set of 338 clear-sky targets. The time series for the six different scene types are shown in Figure 14a before seasonal correction, (b) after seasonal correction, and (c) after the ocean time series are corrected for aerosols, for both the data at the nominal position and the IODC. The time series for the IODC data period do not connect with the first set of time series because different sites are used. Either way, it is clear that the IODC time series keep degrading, in a similar way as the time series at the nominal position. The end of the IODC time series seems to show an added variance, which might be due to an increased inclination angle of the orbit at the end of life when fuel became scarce.

The full set of time series (0° + 63°) are now corrected for aging, where this time the cost function that is being minimized contains 2 × 6 time series so that the variance of both data coverages needs to be as low as possible. Figure 14d shows the aging corrected time series, using the parameters

$$\text{slope}=-0.0110\pm 0.0007\hspace{0.17em}{\text{yr}}^{-1}$$

(29)

$$\alpha =0.000121\pm 0.000009\hspace{0.17em}{\text{day}}^{-1}$$

(30)

$$\beta =0.75\hspace{0.17em}(\text{fixed})$$

(31)

$$\gamma =0.000055\pm 0.000020\hspace{0.17em}\mu {\text{m}}^{-1}{\text{day}}^{-1}$$

(32)

These parameters are really close to the ones found using the observations at the nominal position alone, and so they are used to correct the full Meteosat-5 data set. The percentages of the slopes per year at 0° longitude are given in Table 4, while the IODC slopes are shown in

Table 5. The yearly percentage of decrease of the Meteosat-5 and Meteosat-7 IODC time series, before and after ageing correction.

	Meteosat-5 IODC		Meteosat-7 IODC
Surface Type	Before Ageing Correction	After Ageing Correction	Before Ageing Correction	After Ageing Correction
convective clouds	−0.9833 ± 0.0267	−0.2384 ± 0.0279	−0.917 ± 0.047	−0.087 ± 0.0476
ocean	−0.5651 ± 0.0453	0.2848 ± 0.0451	−0.612 ± 0.051	0.263 ± 0.053
dark vegetation	−0.6889 ± 0.0257	−0.1921 ± 0.0268	−0.340 ± 0.050	0.205 ± 0.052
bright vegetation	−0.6563 ± 0.0226	−0.1042 ± 0.0236	−0.481 ± 0.041	0.102 ± 0.042
dark desert	−0.7204 ± 0.0240	−0.0856 ± 0.0255	−0.546 ± 0.046	0.148 ± 0.049
bright desert	−0.9257 ± 0.0193	−0.2355 ± 0.0206	−0.698 ± 0.029	0.069 ± 0.031

.

8. Meteosat-6

Meteosat-6 was launched in November 1993 as the third and last satellite in the Meteosat Operational Programme, and remained the back-up for Meteosat-5 until January 1997 due to a high noise level in the WV channel which was observed soon after launch. Only one year and a half of data is available to work with (January 1997–June 1998). For Meteosat-6 too, the suggestion of Govaerts [32] was followed to replace the spectral response curve at launch by the one of Meteosat-7 in the spectral model, and use the adjusted calibration coefficient and offset values given on the EUMETSAT website (see Table 2 for the exact values). Figure 15a shows the original time series before seasonal correction, (b) after seasonal correction, and (c) after aerosol correction. As there was not enough data available to calculate the mean annual cycle, the seasonal correction factors of Meteosat-7 were used, but by comparing figures (a) and (b), it is clear that this does not remove the Meteosat-6 seasonal effects.

When correcting the time series, the minimization technique is not able to find any significant model parameters. Even when β is kept fixed to 0.75, the Powell method returns zero values for the three other parameters. As long as the seasonal correction cannot be performed, it is not possible to find good parameters using the minimization technique. As mentioned earlier, Meteosat-6 was also operational over the Indian Ocean, but unfortunately this was only as back-up, and not enough data is available to do the clear-sky process and extend this way the Meteosat-6 time series. In the discussion at the end of the paper, it is explained how Meteosat-6 can still be corrected making use of the time series of Meteosat-5 and -7.

9. Meteosat-7

The last of the first generation satellites was launched in September 1997 and was part of the Meteosat Transition Programme (MTP). Meteosat-7 was operational at 0° longitude from June 1998 until July 2006. The 6 original time series are shown in Figure 16a before seasonal correction, in Figure 16b after seasonal correction and Figure 16c after aerosol correction.The same periodic variance is visible in Figure 16b,c as for Meteosat-2, -4 and -5, due to the method of deseasonalization. The percentage to which the fits through the time series of Figure 16c decrease per year, is shown in Table 4. The spectral effect is again clearly visible from these values, where the ocean time series is decreasing more strongly.

After correcting these aerosol corrected time series with the aging model from Equation (1), they become the ones shown in Figure 16d. The model parameters that come out of minimizing the cost function, are

$$\text{slope}=-0.0295\pm 0.059\hspace{0.17em}{\text{yr}}^{-1}$$

(33)

$$\alpha =0.000327\pm 0.000009\hspace{0.17em}{\text{day}}^{-1}$$

(34)

$$\beta =0.7529\pm 0.0053$$

(35)

$$\gamma =0.000125\pm 0.000014\hspace{0.17em}\mu {\text{m}}^{-1}{\text{day}}^{-1}$$

(36)

9.1. Indian Ocean Data Coverage (IODC)

In December 2006, Meteosat-7 took over the operations of Meteosat-5 in the IODC at 57° E, and is expected to continue at least until 2016. By extending the time series at the nominal position with the IODC time series, the model parameters can be improved, having now, like Meteosat-5, about 13 years of data of the same satellite. New sites are needed over the Indian Ocean at 57° E, as there is almost no overlap with the 0° longitude FOV. The 314 clear-sky selected targets are shown in Figure 17. Figure 18 shows the original time series for the targets at the nominal position and the IODC targets (a) before, and (b) after the seasonal correction. The time series of the same scene type do not connect for the two different time periods due to the fact that different targets were used. The full data set nicely shows the exponential behavior, where the degradation saturates in the IODC part of the time series. No aerosol correction is done because the GACP dataset only runs until 2006. As there are no clear effects of volcanic eruptions in the data, i.e., no prominent peaks or dips in the time series, it should not be a problem to minimize the time series without the aerosol correction. The model parameters can still be found in the same way as before.

The time series are corrected for aging, where again 2 × 6 time series are used in the minimisation process. The model parameters that come out of the Powell method are

$$\text{slope}=-0.031918\hspace{0.17em}{\text{yr}}^{-1}$$

(37)

$$\alpha =0.000374\hspace{0.17em}{\text{day}}^{-1}$$

(38)

$$\beta =0.766187$$

(39)

$$\gamma =0.000074\hspace{0.17em}\mu {\text{m}}^{-1}{\text{day}}^{-1}$$

(40)

These parameters are comparable to the ones found for the 0° time series of Meteosat-7, and are here too used to correct the full Meteosat-7 dataset. Figure 18c shows the aging corrected time series. The values of the slopes in percentage per year of the aerosol corrected time series at the nominal position are given in Table 4, while the ones for the IODC (not corrected for aerosols) are given in the Table 5.

10. Discussion

Figure 19a shows the aging corrected time series for all 6 MFG satellites. Compared to Meteosat-4, -5 and -7, the Meteosat-6 time series have consistently lower values. The reason for this is the fact that no significant model parameters were found in the minimization process and so no aging correction could be done, while the used calibration was extrapolated at launch, which is three years prior to the start of the time series, using the SSCC drift. Furthermore, spectral degradation is clearly present, as the clear ocean and cloudy time series are decreasing even more than the ones over land. Even though it is not possible to model the aging based on the 1.5 years of data, the aging parameters can be derived in such a way that the Meteosat-6 time series agree as closely as possible with the previous (i.e., Meteosat-5) and following (i.e., Meteosat-7) ones. The best consistency is obtained for Meteosat-6 using the parameters

$$\text{slope}=-0.023\hspace{0.17em}{\text{yr}}^{-1}$$

(41)

$$\alpha =0.000250\hspace{0.17em}{\text{day}}^{-1}$$

(42)

$$\beta =0.750000$$

(43)

$$\gamma =0.000100\hspace{0.17em}\mu {\text{m}}^{-1}{\text{day}}^{-1}$$

(44)

Figure 19b shows the new aging corrected time series, where, this time, the Meteosat-6 data connects reasonably with the others. Small shifts are still present in the bright desert time series, which were used for the calibration. This is due to the fact that the calibration coefficients of the SSCC method are used, but the two models correct the data for aging differently: linear grey modeling in the SSCC method, and exponential spectral modeling in the spectral aging model. These small shifts are corrected for by normalizing the bright desert time series, based on the one of Meteosat-7, and multiplying each time series of the same satellite by the same value. Figure 19c shows the normalized time series, where the normalization coefficients are given in

Table 6. The normalisation coefficients for each satellite with respect to the Meteosat-7 bright desert time series.

Satellite	Normalisation Coefficient
Meteosat-2	0.9889
Meteosat-3	1.0096
Meteosat-4	1.0195
Meteosat-5	1.0303
Meteosat-6	0.9903
Meteosat-7	1.0000

. This normalization lies on average between 1% and 3%, and, as expected, the highest correction is done for Meteosat-5, for which the time between launch and operations was the longest. The calibration coefficients at launch, after correction, are given in the last column of

Table 7. For each satellite, the optimal values for the 4 model parameters, the square root of the minimized costfunction (stddev), and the normalized calibration coefficients are given.

Satellite	Slope (yr−1)	α (day−110−3)	β	γ (μm−1day−1 10−3)	stddev	Normalised Calibration
Meteosat-2	−0.017 ± 0.008	0.44 ± 0.34	0.90 ± 0.26	0.000 ± 0.000	0.023 ± 0.012	0.645/0.539
Meteosat-3	−0.009 ± 0.027	0.10 ± 0.29	0.75 (fixed)	0.000 ± 0.396	0.024 ± 0.047	0.634/0.764
Meteosat-4	−0.026 ± 0.002	0.28 ± 0.08	0.74 ± 0.15	0.049 ± 0.037	0.020 ± 0.014	0.746
Meteosat-5	−0.011 ± 0.001	0.12 ± 0.01	0.75 (fixed)	0.055 ± 0.020	0.017 ± 0.003	0.839
Meteosat-6	−0.023 ± 0.018	0.25 ± 0.05	0.75 (fixed)	0.100 ± 0.025	0.030 ± 0.001	0.829
Meteosat-7	−0.032 ± 0.002	0.37 ± 0.06	0.77 ± 0.02	0.074 ± 0.011	0.017 ± 0.011	0.918

.

The final model parameters used for each of the satellites in Figure 19b,c are also given in this table. Having a look at the first column, it appears that the slopes are quite variable from one satellite to the other, ranging from about −0.9% up to −3.2% yr⁻¹, where the two oldest instruments seem to have somewhat smaller slopes than the two newest ones. When it was possible to derive them, β values of about 0.75 were obtained, which represent the asymptotic sensitivity in the aging model. The higher β value for Meteosat-2 (0.90) exhibits a higher uncertainty, which is probably due to the 6-bit digitization. For the spectral aging (column 5), relatively similar values are derived for Meteosat-4, -5 and -7, where γ ranges between 0.000049 and 0.000074 μm day⁻¹. In terms of residual standard deviation, the most stable satellite data records after aging correction are the ones of Meteosat-5 and -7 (about 1.7% at 1σ). Due to the lack of good deseasonalization, a higher variance is observed for Meteosat-6 (3%). As expected, higher standard deviations are observed for Meteosat-2 and -3 due to the 6-bit digitization.

From Figure 19c, it is clear that the major trends in the time series have been reduced through the aging correction. This is confirmed by comparing the values of the slopes before and after the correction in Table 4, where the spectral aging effect has been removed.

Apart from Meteosat-2 and -3, all 6 time series are now also at about the same level in the figure. For Meteosat-3 there seems to be a difference in level for the convective clouds and the ocean when comparing with the more recent satellites, while for Meteosat-2 this is mainly the case for the vegetation and the ocean data. Part of this is due to the 6 to 8-bit conversion of these two satellites. Even the lowest electrical currents, observed by these two instruments, arrive at a digital count value of 1 (nothing smaller is possible). On average, these values are more likely equal to about 0.5. In the 6 to 8-bit conversion, all values are multiplied by 4, so that in the 8-bit version, the offset is equal to 4, but might be closer to 2 in reality Govaerts [33]. This effect is largest for the darkest surface types which reflect the least visible radiation, like ocean and vegetation. By decreasing the offset, the darkest time series shift upwards. This has been tested, and indeed, a better agreement is observed, but it does not solve everything. Another important issue is the pre-launch characterization of the spectral response curve of these instruments. If this characterization was not done accurately enough, it could explain the lower ocean and vegetation time series of Meteosat-2 and -3. However, the investigation is difficult due to the loss of information in the signal quantification of 6 bits. Decoster et al. [25] did a study on the Meteosat-7 visible spectral response curve at launch and found evidence of a similar, though probably much smaller, problem with the pre-launch characterization.

For the other satellites, only the dark desert time series still shows a decrease for Meteosat-6, which is due to the fact that it is difficult to find the right parameters for this dataset. The ocean time series are on the whole, for the full MFG dataset, much lower than the other scene types, which is due to the fact that the data was corrected for the presence of aerosols, while the CERES ADMs that were used to divide the unfiltered reflectance with to create the reflectance ratio, contain an average amount of aerosols. The overall standard deviation is computed for each scene type to check the long-term stability of the time series. The values are shown in

Table 8. The standard deviation of each full MFG time series for the Meteosat-4 up to 7 time series in the first column, and the Meteosat-2 up to 7 time series in the second column.

Surface Type	Meteosat-4–7 (17 yrs)	Meteosat-2–7 (24 yrs)
convective clouds	0.0123	0.0239
ocean	0.0167	0.0611
dark vegetation	0.0140	0.0437
bright vegetation	0.0120	0.0266
dark desert	0.0142	0.0230
bright desert	0.0098	0.0099

. In the first column, only the data of Meteosat-4 up to 7 is used, as these have the least differences from one satellite to another. The stability ranges here between about 0.9% and 1.6%. Adding Meteosat-2 and -3, leads to values between 0.9% and 6%, with the largest ones for the dark scenes (ocean and dark vegetation), which is due to the 6-bit digitization of these two satellites.

The calibration coefficients used in this work for each of the 6 MFG satellites, were calculated using the SSCC method, which is based on bright desert targets, and validated with clear ocean targets. To compare the spectral ageing model with the linear SSCC calibration increase, the grey SSCC correction is applied to the same targets, where the ocean time series is also corrected for aerosols. Figure 19d shows the SSCC corrected time series. The difference for the Meteosat-7 time series is most clearly visible. First of all, there is a slight bending in the time series corrected using the SSCC method, which is removed by the exponential decay of the spectral ageing model [16]. Apart from that, the ocean time series decrease more strongly in the SSCC corrected version for Meteosat-6 and -7. Overall, the ocean time series of Meteosat-4 to -7 are also more connected to the same level in Figure 19c, which supports the need of a spectral model.

11. Conclusions and Future Prospects

The spectral aging model that was developed by Decoster et al. [16], has been applied to 299 clear-sky and cloudy target time series. Using these selected scenes, the model parameters have been derived for all 6 MFG satellites. For Meteosat-2 and -3 no significant γ parameter was found due to the problems encountered (volcanic eruption, 6-to-8bit conversion, saturation), but for the other parameters, values were derived with acceptable standard errors. To reduce the variability in the ocean time series, the effect of the volcanic eruptions of El Chichón and Pinatubo was reduced by correcting these time series for aerosols, based on the GACP dataset. It was also proven, that the lack of AOD data over land, does not pose any problems to the derivation of the model parameters. The data over the Indian Ocean was able to improve the model parameters for Meteosat-5 and -7, while for Meteosat-3 it was shown that the data over the Atlantic Ocean was not good enough to increase the time series at the nominal position. Apart from Meteosat-6, it was possible to derive model parameters for all MFG satellites using the Powell minimization method. For Meteosat-6, the parameters were found by comparing the dataset to the ones of Meteosat-5 and -7. The resulting long-term stability in the corrected time series is equal to about 1%–2% for Meteosat-4 up to -7, but increases up to 6% when adding the Meteosat-2 and -3 time series. The reason for this increase is a combination of a signal digitization on only 6 bits (converted to 8 bits through multiplication), which decreases the time series of the darkest scene types due to a possible problem with the calibration offset, and a pre-launch spectral response characterization with a significantly larger uncertainty. As long as the calibration offset problem and the inaccurate spectral response characterization are not corrected for, it is not recommended to use the Meteosat-2 and -3 time series in FCDR and TCDRs for the VIS channel of MVIRI.

In the near future, the spectral aging model will be used to derive the TOA radiation GERB-like database from the MVIRI VIS images, in support of the CM SAF. This dataset will be based on an empirical unfiltering between the GERB SW images of Meteosat-8 and MVIRI VIS images of Meteosat-7 during the overlap period of two years (2004–2006). Based on these empirically derived GERB-like Meteosat-7 images, the spectral aging model will be applied through theoretical unfiltering regressions to derive the full MVIRI GERB-like dataset. In collaboration with EUMETSAT, the MFG FCDR could also be derived using this model. The model would then be applied through linear regressions based on a reference spectral response curve like, for example, the Meteosat-7 pre-launch characterized spectral response curve, or the Meteosat-8 High Resolution Visible (HRV) spectral response curve [25]. As, with respect to the original Level 1.5 images, the aging corrections will be smaller, in this case it might be better to express the FCDR in reflectance ρ instead of DC. This way, the discretization effect will disappear, and the aging correction will be done in floating point. Finally, the spectral aging model can also be provided in its raw format, with the right model parameters for each satellite. Even though this is more complicated for the user, the model can be used this way to create LUTs that vary in time for the retrieval of, for example, aerosol or cloud properties where the full corrected images are not necessary.