Detection of Snow Cover from Historical and Recent AVHHR Data—A Thematic TIMELINE Processor

Global snow cover forms the largest and most transient part of the cryosphere in terms of area. On the local and regional scale, small changes can have drastic effects such as floods and droughts, and on the global scale is the planetary albedo. Daily imagery of snow cover forms the basis of long-term observation and analysis, and only optical sensors offer the necessary spatial and temporal resolution to address decadal developments and the impact of climate change on snow availability. The MODIS sensors have been providing this daily information since 2000; before that, only the Advanced Very High-Resolution Radiometer (AVHRR) from the National Oceanographic and Atmospheric Administration (NOAA) was suitable. In the TIMELINE project of the German Aerospace Center, the historic AVHRR archive in HRPT (High Resolution Picture Transmission) format is processed for the European area and, among other processors, one output is the thematic product ‘snow cover’ that will be made available in 1 km resolution since 1981. The snow detection is based on the Normalized Difference Snow Index (NDSI), which enables a direct comparison with the MODIS snow product. In addition to the NDSI, ERA5 re-analysis data on the skin temperature and other level 2 TIMELINE products are included in the generation of the binary snow mask. The AVHRR orbit segments are projected from the swath projection into LAEA Europe, aggregated into daily coverages, and from this, the 10-day and monthly snow covers are finally calculated. In this publication, the snow cover algorithm is presented, as well as the results of the first validations and possible applications of the final product.


Introduction
The importance of the product "Area covered by snow" is emphasized by the fact that it is considered one of 54 Essential Variables (EV) by the Global Climate Observing System (GCOS) [1]. As part of the cryosphere, snow forms its largest component in terms of area, but it is also the most short-lived part of the cryosphere. The analysis of the development of the snow cover between different years, as well as within different years, enables statements on the development and the recognition of trends. Due to the short lifespan of the snow cover and its spatial heterogeneity, both a high repetition rate and an appropriate geometric resolution are required. Therefore, medium-resolution optical remote sensing systems are most suitable for providing this information [2]. The National Snow and Ice Data Center (NSIDC) has provided daily data on global snow cover since 2000 from the Moderateresolution Imaging Spectroradiometer (MODIS) on the Terra satellite. Since 2002, these data have been made available by MODIS on the Aqua satellite [3]. Although MODIS already covers a large period of time, other sensors must be used for periods that are further in the past, such as climatologically relevant periods.
In 2013, the Earth Observation Center (EOC) of the German Aerospace Center (DLR), initiated the TIMELINE project which stands for "TIMe Series Processing of Medium Resolution Earth Observation Data assessing Long-Term Dynamics In our Natural Environment". The objective is to generate a long and homogenized National Oceanic and Atmospheric Administration (NOAA) and Meteorological Operation Satellite (MetOP) Advanced Very High Resolution Radiometer (AVHRR) time series over Europe and North Africa using the historical AVHRR archive dating back to 1982 [4]. The AVHRR sensor has a long history in determining snow cover through multispectral classification [5] but also in detecting sea ice [6] and estimating the snow-water equivalent [7]. Above all, the snow-covered area is often a by-product of processors for cloud detection, such as the SPARC algorithm [8] or the APOLLO algorithm [9] developed at DLR. However, the studies with AVHRR data on snow covered areas mostly take place on a regional scale, such as on the Tibetan Plateau [10], in Turkey [11], or in the Alps [12,13]. In complex terrain, attempts are often made to obtain information at a subpixel level [10,14]. Since AVHRR data are often used to extend the MODIS time series, there are also different approaches to harmonize the two datasets, such as in Central Asia [10,15].
In the TIMELINE project, the processors are divided into three levels. The snow cover processor presented in this article is made available as a Level 2 and a Level 3 product. The initial data set is the Level 1b product, which is atmospherically corrected in the satellite swath projection. In addition, the processor requires preceding products such as the cloud mask [16] and the water mask [17]. Both the use of the AVHRR sensor to determine snow cover and the analysis of the time series of snow cover are established methods. The novelty of the TIMELINE project and the snow product is that there is no medium-resolution time series for Europe that covers this long period of almost 40 years and is based on data recorded several times a day (high temporal resolution). This and the associated opportunities are unique. While the MODIS snow product delivers different results for the Terra and Aqua platforms [3], our unique characteristic is that the harmonization of 12 platforms (NOAA-07 to NOAA-19, NOAA-13 never became operational) leads to platform-independent results [4]. Such a time series is of immense importance for research into climate change, as it almost doubles the previous standard for daily recorded medium-resolution snow cover (the MODIS time series for the last 21 years). In addition, the derivation of the snow product, which is based on the MODIS NDSI-NDVI algorithm [3,18], allows a direct comparison and an expansion of the two products. Our research question is the following: Can the TIMELINE AVHRR snow product be used to complete the MODIS snow cover time series available since 2000 up to the 1980s? In the following we will show the theoretical basis and the methodology of the Level 2 and Level 3 products, as well as the first results of the validation with the MODIS snow product.

Normalized Difference Snow Index (NDSI)
A very common method for detecting snow surfaces from the reflective part of the electromagnetic spectrum is to use snow's unique spectral properties. Figure 1 shows the typical reflection properties of snow compared to vegetation (conifers in this example). In addition, the spectral response functions of the first three bands of the sensor AVHRR/3 (on NOAA-17) are shown. In the reflection spectrum of snow, there is almost complete reflection in the visible spectral range (VIS: 0.38-0.78 µm). The reflection decreases rapidly in the near infrared (NIR: 0.78-1.4 µm), and in the short-wavelength infrared (SWIR: 1.4-3.0 µm), the radiation is almost completely absorbed. This extremely different reflection between VIS and SWIR is used in many methods for determining snow surfaces. Most methods for distinguishing between clouds and snow use a decision tree with fixed thresholds [8,9], which have to be adjusted individually depending on the region under consideration, such as the Alps [12,13]. The snow products from MODIS [3] are no longer based on fixed threshold values; instead, a normalized index is used that still maps the reflection gradient between VIS and SWIR well, even with lower reflections. Equation (1) shows the calculation of the Normalized Difference Snow Index (NDSI) [18] based on the AVHRR/3 spectral bands: where is the atmospherically corrected reflectance in AVHRR Band 1 and the reflectance in Band 3a. The third generation of the AVHRR sensor has two options for using band 3: With variant 3a, the spectral range between 1.58 and 1.64 µm is recorded, but with variant 3b, however, the mid-wavelength infrared (MWIR) is between 3.55 and 3.93 µm. The two bands cannot be operated synchronously, so that there is either reflectance data or thermal data. This causes a problem in the NDSI calculation, since the reflectance data of band 3a is required (Equation (2)). This problem is solved by determining the reflective component of band 3b. The blackbody radiation of the brightness temperature determined from band 5 is used for this from the wavelength range of band 3b. This radiation is then subtracted from the recorded radiation of band 3b and normalized using the solar irradiance [8]: In the reflection spectrum of snow, there is almost complete reflection in the visible spectral range (VIS: 0.38-0.78 µm). The reflection decreases rapidly in the near infrared (NIR: 0.78-1.4 µm), and in the short-wavelength infrared (SWIR: 1.4-3.0 µm), the radiation is almost completely absorbed. This extremely different reflection between VIS and SWIR is used in many methods for determining snow surfaces. Most methods for distinguishing between clouds and snow use a decision tree with fixed thresholds [8,9], which have to be adjusted individually depending on the region under consideration, such as the Alps [12,13]. The snow products from MODIS [3] are no longer based on fixed threshold values; instead, a normalized index is used that still maps the reflection gradient between VIS and SWIR well, even with lower reflections. Equation (1) shows the calculation of the Normalized Difference Snow Index (NDSI) [18] based on the AVHRR/3 spectral bands: where R 1 is the atmospherically corrected reflectance in AVHRR Band 1 and R 3a the reflectance in Band 3a. The third generation of the AVHRR sensor has two options for using band 3: With variant 3a, the spectral range between 1.58 and 1.64 µm is recorded, but with variant 3b, however, the mid-wavelength infrared (MWIR) is between 3.55 and 3.93 µm.
The two bands cannot be operated synchronously, so that there is either reflectance data or thermal data. This causes a problem in the NDSI calculation, since the reflectance data of band 3a is required (Equation (2)). This problem is solved by determining the reflective component of band 3b. The blackbody radiation of the brightness temperature determined from band 5 is used for this from the wavelength range of band 3b. This radiation is then subtracted from the recorded radiation of band 3b and normalized using the solar irradiance [8]: where L 3b is the recorded radiance of band 3b (units: W/ m 2 str µ 0 ), B(λ 3 , T 5 ) is the blackbody radiance at wavelength λ 3 and temperature T 5 (i.e., brightness temperature of band 5), B sun is the solar radiance corrected for the Earth-sun distance, and µ 0 is the cosine of the sun zenith angle. The sensor-specific values of the wavenumbers and the solar radiance for the calculation of B(λ 3 , T 5 ) are taken from Trishchenko (2006) [19]. From the NDSI, the percentage of area covered with snow can also be determined. The NDSI threshold of 0.4 corresponds to a fractional snow cover of 50%, and all NDSI values above this threshold indicate more than 50% snow coverage [3,20,21]. Version 5 of the MODIS snow product used this threshold value to make a binary classification between snow (NDSI ≥ 0.4) and non-snow (NDSI < 0.4). Since version 6 (the latest version is 6.1), there is no longer a binary distinction, but all positive NDSI values are displayed in the "NDSI_snow_cover" layer, scaled between 0 and 100. This is due to the fact that a strict threshold value of 0.4 can lead to an underestimation of the detected snow surfaces in (especially coniferous) forests, which has been found by Klein et al. (1998) [22] by combining a snow reflectance model with a canopy reflectance model (GeoSAIL). Therefore, in addition to the NDSI, the Normalized Difference Vegetation Index (NDVI) is now also used to identify these snow areas in forests. The NDVI is calculated using Equation (3): where R 1 is the atmospherically corrected reflectance in AVHRR Band 1 and R 2 the reflectance in Band 2. Based on NDSI and NDVI, depending on the NDVI value, a linear (ndsi linear ) and an exponential function (ndsi exponential ) is used to create an area for the "snow within forest" class. The following functions (4a) to (4c) are based on Klein et al. (1998) [22]: Based on the decision from Equation (4c), the area that forms the "snow within forest" class is shown in Figure 2. This additional area is also used for processing the AVHRR data.

Snow Cover-TIMELINE Level 2 Processor
Based on the atmospherically corrected reflectance of the AVHRR Level 1b product, the Level 2 cloud mask [16], the Level 2 water mask [17], a static water mask, and ERA5 reanalysis data on the skin temperature [23] are required as input data. The reanalysis

Snow Cover-TIMELINE Level 2 Processor
Based on the atmospherically corrected reflectance of the AVHRR Level 1b product, the Level 2 cloud mask [16], the Level 2 water mask [17], a static water mask, and ERA5 reanalysis data on the skin temperature [23] are required as input data. The reanalysis data are needed because distinguishing between clouds and snow using AVHRR is often challenging [24]. Therefore, the hourly skin temperature data are interpolated to the exact mean acquisition time of the scene, and the data set is converted into the swath projection of the orbit. The skin temperature is also linearly interpolated spatially in order to then create a difference raster between the brightness temperature of band 4 and the skin temperature. The cloud mask of APOLLO_NG [16] no longer specifies exact classes like APOLLO [9] but follows a probabilistic approach and provides probabilities for the occurrence of clouds.
This probabilistic cloud mask, the water masks (L2 product and static water mask), the calculated temperature difference, and the NDSI/NDVI classification explained in the previous chapter form the basis of the Level 2 snow processor. The steps are shown in Figure 3. The first part of the processing is the masking of obvious cloud and water areas. This is performed with the Level 2 products mentioned, the static ocean/sea water mask, the NDVI, and the temperature difference. Areas with a solar zenith angle greater than 85° are classified as "night". Then, according to Equation (4c), a provisional snow determination takes place. All land areas outside are classified as "snow free". In the second part of the snow determination, additional probability thresholds are used, and snow pixels are reclassified if necessary. First, a temperature-height screen is applied: If the brightness temperature of band 4 is above 281 K and the terrain altitude is below 1300 m, the pixel is assigned to the "snow-free" class. The next screen concerns high SWIR reflectance: Pixels with more than 45% reflectance are determined as "snow free". The last screen concerns the reflection in the visible range: If this is below 7% in band 1 or 2, the pixel is also declared as "snow-free". For all remaining "snow-free" pixels, the cloud mask is used again, and all pixels with a cloud probability of over 70% are assigned to the "cloud" class. The cloud probability is an output layer of the Level 2 product cloud mask, which is computed with the probabilistic cloud processor APOLLO_NG [16]. The threshold value for the cloud probability of 70% was determined empirically, since clouds are best recognized here, and snow surfaces in this value range are usually not misclassified as clouds.
This division into classes forms the output layer "snow_mask" of the Level 2 product. In addition, the layer "raw_ndsi" is provided, which contains the raw NDSI values-the basis for snow determination. In addition, a "snow_quality_flag" layer is supplied, which hierarchically maps the quality of the classification in big-endian 8-bit coding ( Table 1).
The threshold values correspond to those of the MODIS snow product [3]. The daily ag- The first part of the processing is the masking of obvious cloud and water areas. This is performed with the Level 2 products mentioned, the static ocean/sea water mask, the NDVI, and the temperature difference. Areas with a solar zenith angle greater than 85 • are classified as "night". Then, according to Equation (4c), a provisional snow determination takes place. All land areas outside are classified as "snow free". In the second part of the snow determination, additional probability thresholds are used, and snow pixels are re-classified if necessary. First, a temperature-height screen is applied: If the brightness temperature of band 4 is above 281 K and the terrain altitude is below 1300 m, the pixel is assigned to the "snow-free" class. The next screen concerns high SWIR reflectance: Pixels with more than 45% reflectance are determined as "snow free". The last screen concerns the reflection in the visible range: If this is below 7% in band 1 or 2, the pixel is also declared as "snow-free". For all remaining "snow-free" pixels, the cloud mask is used again, and all pixels with a cloud probability of over 70% are assigned to the "cloud" class. The cloud probability is an output layer of the Level 2 product cloud mask, which is computed with the probabilistic cloud processor APOLLO_NG [16]. The threshold value for the cloud probability of 70% was determined empirically, since clouds are best recognized here, and snow surfaces in this value range are usually not misclassified as clouds.
This division into classes forms the output layer "snow_mask" of the Level 2 product. In addition, the layer "raw_ndsi" is provided, which contains the raw NDSI values-the basis for snow determination. In addition, a "snow_quality_flag" layer is supplied, which hierarchically maps the quality of the classification in big-endian 8-bit coding ( Table 1). The threshold values correspond to those of the MODIS snow product [3]. The daily aggregation of the Level 2 products in the Level 3 processor is based on this layer. As can be seen from Table 1, the "snow_quality_flag" layer usually contains the same information as the "snow_mask" layer ( Figure 3), but these are hierarchically structured and some differ: For example, if the brightness temperature of band 4 exceeds 281 K, the pixel is only reclassified to "snow free" if the elevation is below 1300 m-however, quality bit 5 is set in any case. If bit 4 is set, this does not cause a reclassification, but has a major impact on the detection of snow, since large sun zenith angles lead to poorer snow detection. Wang & Zender (2010) [25] have already determined quality losses from a sun zenith angle of more than 55 • . However, since this threshold value would mean a large loss in surface area in winter, a sun zenith angle of over 70 • was chosen for the MODIS snow product [3], from which point the quality flag is set. This is due to the unreliable albedo from 70 • , which has been proven by radiative transfer and BRDF modelling [26][27][28][29]. The same is true for large satellite zenith angles, since at AVHRR/3, the pixel size increases from 1.08 km in nadir to 6.15 km at the edge (across-track) [30], and with it, the length of the path through the atmosphere. Since the cloud processor, for example, also has problems in these extreme peripheral areas of the scene (satellite zenith angle of~68 • ), the flag "unfavorable illumination and/or observation geometry" was set for areas when the satellite zenith angle exceeds the empirically determined threshold of 65 • .

Snow Cover-TIMELINE Level 3 Processor
In contrast to the Level 2 products, the Level 3 products are available in geographical projection. ETRS89-extended/LAEA Europe (EPSG: 3035) was chosen as the target coordinate system, whereby the extent of the TIMELINE area (Europe, Middle East, and North Africa) was divided into four tiles (Figure 4). ellite zenith angle exceeds the empirically determined threshold of 65°.

Snow Cover-TIMELINE Level 3 Processor
In contrast to the Level 2 products, the Level 3 products are available in geographical projection. ETRS89-extended/LAEA Europe (EPSG: 3035) was chosen as the target coordinate system, whereby the extent of the TIMELINE area (Europe, Middle East, and North Africa) was divided into four tiles (Figure 4).  The basis of all Level 3 products is the daily composite of all scenes available on that day. These are first projected into the respective tile. An output pixel size of 1000 m is selected, and "nearest neighbor" is selected as the resampling method. Only the "snow_mask" and "snow_quality_flag" layers are projected, with a no-data value of 0 or 255 being set. This gives the "snow_quality_flag" the highest value at this point, which is critical for further processing, since the processor is looking for the lowest value to fill the day composite. Two three-dimensional arrays are then created and filled with the two projected layers. The position with the minimum value is now searched for in the "snow_quality_flag" stack, and the final "snow_mask" is filled with the corresponding values of the "snow_mask" stack.
In addition to the daily composites, a 10-day product and a monthly product are also created. For this purpose, stacks are formed from 10 consecutive days or the whole month. From this, the minimum snow cover, the maximum snow cover, and the snow cover percentage of each pixel (without cloud cover) for the period under consideration are determined.

Validation with MODIS
The daily MODIS snow products [3] from Terra (MOD10A1, obtained from http://nsidc. org/data/MOD10A1, accessed on 19 December 2021) and Aqua (MYD10A1, obtained from http://nsidc.org/data/MYD10A1, accessed on 19 December 2021) in the latest version 6.1 are used to validate the snow mask. The entire TIMELINE area is covered by about 27 MODIS tiles (h16v01 to h21v05). These are downloaded, mosaicked, and converted into the target projection and pixel size for the period to be compared. A confusion matrix (Table 2) of the classes "snow" and "snow-free" is now created for all pixels that are cloud-free both in the AVHRR data and in the MODIS data. The following indicators of classification accuracy can be determined from the confusion matrix: the true positive rate (TPR), the true negative rate (TNR), the positive predictive value (PPV), the negative predictive value (NPV), and the accuracy (ACC) [31]. The following Equations (5a)-(5e) show their determination from the confusion matrix:

Results
The Level 2 processor is applicable to all AVHRR Level 1B data of the TIMELINE Project, which is available at the German Aerospace Center (DLR) in the Local Area Coverage (LAC) with a 1 km spatial resolution. In the LAC mode, the recording time is limited to 10 min, but usually, the area is sufficient for the meridional expansion of the TIMELINE area (Figure 4). Although a channel in the SWIR (band 3a) is required for the NDSI calculation (Equation (1)), it can be calculated according to Equation (2) from the band 3b. We can therefore use the entire AVHRR archive that ranges from October 1981 to December 2018-this is a time series of 37 years. Since only a small part of the entire data set could be used for processor development, we chose 2009 as the test year. The winter of 2009 (and spring transition) showed all imaginable snow conditions in the examination area. For this period, MODIS data of the Terra and Aqua platforms are also available, which are used for validation.

Timeline Snow Cover-Level 2 and Level 3 Products
For testing, we focused on the months of January to April 2009. For the snow cover classifications, only day recordings can be used. However, since the number of Level 2 products created is very large (241 in January), Figure 5 shows an example of the processing steps from a scene on 8 January 2009. are used for validation.

Timeline Snow Cover-Level 2 and Level 3 Products
For testing, we focused on the months of January to April 2009. For the snow cover classifications, only day recordings can be used. However, since the number of Level 2 products created is very large (241 in January), Figure 5 shows an example of the processing steps from a scene on 8 January 2009.  Figure 5A shows snow-covered areas in the color cyan due to the RGB channel combination making clouds appear gray or white. In the upper section of Figure 5B (beginning of the recording), one can see a sharp transition due to the switching from band 3b (night  Figure 5A shows snow-covered areas in the color cyan due to the RGB channel combination making clouds appear gray or white. In the upper section of Figure 5B (beginning of the recording), one can see a sharp transition due to the switching from band 3b (night mode) to band 3a (day mode). The NDSI image of Figure 5B also highlights the problem associated with this index for snow detection: In addition to snow, this index also takes a high value for water surfaces and some clouds [24]. The final classification result shows a satisfactory snow mask for this scene in Figure 5C. Clouds over seas were eliminated by the static water mask.
As examples for daily Level 3 products, 15 January 2009 (mid-winter) and 22 February 2009 (maximum snow coverage) were selected and are shown in Figures 6 and 7. mode) to band 3a (day mode). The NDSI image of Figure 5B also highlights the problem associated with this index for snow detection: In addition to snow, this index also takes a high value for water surfaces and some clouds [24]. The final classification result shows a satisfactory snow mask for this scene in Figure 5C. Clouds over seas were eliminated by the static water mask.
As examples for daily Level 3 products, 15 January 2009 (mid-winter) and 22 February 2009 (maximum snow coverage) were selected and are shown in Figures 6 and 7.   mode) to band 3a (day mode). The NDSI image of Figure 5B also highlights the problem associated with this index for snow detection: In addition to snow, this index also takes a high value for water surfaces and some clouds [24]. The final classification result shows a satisfactory snow mask for this scene in Figure 5C. Clouds over seas were eliminated by the static water mask.
As examples for daily Level 3 products, 15 January 2009 (mid-winter) and 22 February 2009 (maximum snow coverage) were selected and are shown in Figures 6 and 7.    Figures 6 and 7 show linear artifacts in the overlap area of two level 2 products, which are due to insufficient cloud recognition on the image edge. The problem occurs with absolute satellite zenith angles larger than 65 • and is labelled as a "bad satellite-viewinggeometry" in the quality layer (which also applies to a sun-zenith angle greater than 70 • ). It is also striking that wide parts of Fenno Scandia and northern Russia were classified as "snow-free", which is unlikely at this season. In addition to the poor observation geometry, this is also due to a low reflection in the visible spectral range, whereby these areas are reclassified to "snow-free" despite a high NDSI. This leads to the accuracy analysis, which was carried out for these 4 months with MODIS snow data.

Validation Result
For validation, the daily MODIS snow products of Terra (MOD10A1) and Aqua (MYD10A1) had to be prepared first. For this purpose, the data of 27 MODIS tiles of both platforms were downloaded for the first 120 days of 2009. Then, the information of the layer "NDSI_snow_cover" was translated to the classes of the AVHRR snow product (Figure 3). The data of Terra and Aqua were then united to close possible data gaps through clouds [32]. Finally, the 27 tiles were mosaicked and transferred to the TIMELINE projection ( Figure 4). The pixel size of 500 m was also scaled to 1 km, using a nearest neighbor resample algorithm. All preprocessing was performed with GDAL [33].
A confusion matrix (according to Table 2) was created for each day, which contains all pixels that have either the "snow" or "snow-free" classes in both the MODIS dataset and the AVHRR dataset. The "snow" class is considered a "true" value, and the "snow-free" class is considered a "false" value. Figure 8 shows the boxplots of the classification indicators of true positive rate (TPR), true negative rate (TNR), positive predictive value (PPV), negative predictive value (NPV), and accuracy (ACC) for all 120 days. Figures 6 and 7 show linear artifacts in the overlap area of two level 2 products, which are due to insufficient cloud recognition on the image edge. The problem occurs with absolute satellite zenith angles larger than 65° and is labelled as a "bad satellite-viewinggeometry" in the quality layer (which also applies to a sun-zenith angle greater than 70°). It is also striking that wide parts of Fenno Scandia and northern Russia were classified as "snow-free", which is unlikely at this season. In addition to the poor observation geometry, this is also due to a low reflection in the visible spectral range, whereby these areas are reclassified to "snow-free" despite a high NDSI. This leads to the accuracy analysis, which was carried out for these 4 months with MODIS snow data.

Validation Result
For validation, the daily MODIS snow products of Terra (MOD10A1) and Aqua (MYD10A1) had to be prepared first. For this purpose, the data of 27 MODIS tiles of both platforms were downloaded for the first 120 days of 2009. Then, the information of the layer "NDSI_snow_cover" was translated to the classes of the AVHRR snow product (Figure 3). The data of Terra and Aqua were then united to close possible data gaps through clouds [32]. Finally, the 27 tiles were mosaicked and transferred to the TIMELINE projection ( Figure 4). The pixel size of 500 m was also scaled to 1 km, using a nearest neighbor resample algorithm. All preprocessing was performed with GDAL [33].
A confusion matrix (according to Table 2) was created for each day, which contains all pixels that have either the "snow" or "snow-free" classes in both the MODIS dataset and the AVHRR dataset. The "snow" class is considered a "true" value, and the "snowfree" class is considered a "false" value. Figure 8 shows the boxplots of the classification indicators of true positive rate (TPR), true negative rate (TNR), positive predictive value (PPV), negative predictive value (NPV), and accuracy (ACC) for all 120 days.  Figure 8 shows that the true positive rate (TPR) and negative predictive value (NPR) indicators are always close to 1 (100%). This is due to the fact that the proportion of snowfree areas in the entire TIMELINE area is always very high, even in winter-especially in northern Africa and the Middle East. Only the indicators true negative rate (TNR), positive predictive value (PPV) and accuracy (ACC) are suitable for a differentiated consideration of the results. Figure 9 therefore shows their seasonal development over the months of January to April.  Figure 8 shows that the true positive rate (TPR) and negative predictive value (NPR) indicators are always close to 1 (100%). This is due to the fact that the proportion of snowfree areas in the entire TIMELINE area is always very high, even in winter-especially in northern Africa and the Middle East. Only the indicators true negative rate (TNR), positive predictive value (PPV) and accuracy (ACC) are suitable for a differentiated consideration of the results. Figure 9 therefore shows their seasonal development over the months of January to April. While no seasonality was found for the positive predictive value (PPV) indicator, it exists for true negative rate (TNR) and accuracy (ACC). There is an increase in both over time. On the one hand, this is due to the general decrease in snow cover, but on the other hand, it is also due to its better detection. So, it is worth having a closer look at the difference grids for the lowest and highest values of TNR, PPV, and ACC, respectively. The date 13 January 2009 stands out as a negative example, as this is when TNR, PPV, and ACC reach their lowest values. Figure 10 shows the misclassifications.  While no seasonality was found for the positive predictive value (PPV) indicator, it exists for true negative rate (TNR) and accuracy (ACC). There is an increase in both over time. On the one hand, this is due to the general decrease in snow cover, but on the other hand, it is also due to its better detection. So, it is worth having a closer look at the difference grids for the lowest and highest values of TNR, PPV, and ACC, respectively. The date 13 January 2009 stands out as a negative example, as this is when TNR, PPV, and ACC reach their lowest values. Figure 10 shows the misclassifications. While no seasonality was found for the positive predictive value (PPV) indicator, it exists for true negative rate (TNR) and accuracy (ACC). There is an increase in both over time. On the one hand, this is due to the general decrease in snow cover, but on the other hand, it is also due to its better detection. So, it is worth having a closer look at the difference grids for the lowest and highest values of TNR, PPV, and ACC, respectively. The date 13 January 2009 stands out as a negative example, as this is when TNR, PPV, and ACC reach their lowest values. Figure 10 shows the misclassifications.    Figure 10 shows a large underestimation of snow cover by AVHRR in large parts of Eastern Europe (Ukraine, Belarus, Russia), as well as in the Alps. Due to the earlier point in time in the year (low sun altitude), this is again mainly due to the poor observation and illumination geometry. The later the time in the year, the better the illumination conditions and thus the results. Thus, on 28 April 2009, TNR and ACC show their highest (best) value ( Figure 11). PPV also has its maximum in spring: on 11 April 2009. The snow cover is only underestimated in small areas of Scandinavia, the remaining areas are mapped well ( Figure 12).
The examination of the difference grid shows that the improvement in the indicators is related to a better recognition of the snow surface and is independent of an increase in the "snow-free" area. PPV also has its maximum in spring: on 11 April 2009. The snow cover is only underestimated in small areas of Scandinavia, the remaining areas are mapped well ( Figure 12).

Discussion
Numerous methods have been developed to detect snow from AVHRR data. Snow cover is usually a by-product of cloud masking, as both classes have similar spectral properties [9,[34][35][36]. The determination of the classes "snow" and "snow-free" is carried out either via a hierarchical classification [9], via probabilities (or scores) [8,13,37], or from The examination of the difference grid shows that the improvement in the indicators is related to a better recognition of the snow surface and is independent of an increase in the "snow-free" area.

Discussion
Numerous methods have been developed to detect snow from AVHRR data. Snow cover is usually a by-product of cloud masking, as both classes have similar spectral properties [9,[34][35][36]. The determination of the classes "snow" and "snow-free" is carried out either via a hierarchical classification [9], via probabilities (or scores) [8,13,37], or from multichannel thresholds combined with indices [38,39]. Although good results were achieved with hierarchical approaches [40], it turned out early on that they are not suitable for the present AVHRR data set. At least for the test period considered, we often received unusually low snow reflections in the VIS spectral range (far below the 80% given in Figure 1) for areas with a high sun zenith angle. This may be due to too much impact from the harmonization factors applied to the data set [4,41]. Since the SPARC algorithm [8,13], like APOLLO [9], also works with fixed reflection threshold values of the reflective bands, these were discarded, and an index-based algorithm was chosen with regard to comparability with MODIS [24]. The NDSI method makes snow detection independent of absolute reflection and BRDF effects over snow, which are known to have a strong influence on AVHRR [42]. Numerous evaluated threshold values [3,8,9,24] were used in the post-classification to eliminate particularly low reflections in the VIS, excessive brightness temperatures and too high reflections in the SWIR, which makes the result comparable with other methods for snow cover detection from AVHRR data.
The results of the validation against the MODIS snow product show that the classification quality increases with increasing sun elevation (decreasing sun zenith angle). Sun zenith angles of more than 70 • are generally not well suited for determining snow from optical remote sensing data [25]. Due to the mentioned underestimation of the reflection in the VIS of the Level 1b products due to atmospheric correction and harmonization, snowcovered areas are reclassified as "snow-free". In addition, there are already inaccuracies in the Level 2 products used, so that cloud detection in the area of the scene edges (satellite zenith angle > 65 • ) can fail or snow areas can be recognized as clouds. The larger area at the edge of the scene and the longer path through the atmosphere also lead to a loss in radiation, particularly in the short-wave spectrum. There, the VIS-reflections can fall below the threshold of 7% and be reassigned to the "snow-free" class. This can be seen in Figure 12 in the area of Sweden and Norway, where both a large sun zenith angle and peripheral areas of a scene (high satellite zenith angle) lead to incorrect classification.
This leads to the question of how we can solve the identified problems? A drastic possibility would be that all "inaccurate" Level 2 scene excerpts-which do not meet the quality standards mentioned in Table 1-are excluded from the Level 3 product generation. However, this would mean an immense loss of information, since large areas of the snow-covered northern hemisphere would not be mapped. Eliminating the VIS threshold would improve snow detection, but this decision should be made using a larger data set (several years, platforms, and sensor generations). The Level 3 processor is currently being implemented in the TIMELINE processing environment [4], and then the entire data set will be processed with it. This will give us more information on how good the results are for the different platforms and sensor versions. High-resolution Landsat data will then be used for validation, which will enable subpixel analysis [43,44].
When the entire data series has been processed, a daily cloud-free AVHRR product covering a unique time series of almost 40 years is created from the daily data analogous to DLR's Global SnowPack [45]. This enables the observation of hydrological processes at the catchment area level [15,46,47] and, for example, the investigation of snow changes in mountains [11,12], which was previously not possible for this long period of time. The data can also be used with other environmental data for hydrological modeling in order to investigate the occurrence of extreme situations in climatologically relevant periods [48,49].

Conclusions
In this manuscript, the first results of the thematic snow cover processor of the TIME-LINE project were presented. In TIMELINE, a data set from almost 40 years of AVHRR data is reprocessed according to up-to-date processing methods. One of the thematic products is snow cover, the largest and most ephemeral part of the cryosphere. The results for the test period presented here (the first 4 months of 2009) show good agreement between AVHRR snow cover and the daily MODIS product used as a reference. In general, however, AVHRR underestimates snow under poor lighting conditions (high sun zenith angle). The processor is now applied to the entire dataset to re-validate the results and adjust them if necessary. The time series created in this way will offer a variety of options for investigating the change in snow cover for a wide variety of questions. Data Availability Statement: The TIMELINE Level 2 data used will be made available via the EOC GeoService portal (https://geoservice.dlr.de/web/) in the near future, as will the snow cover data generated (Level 2 and Level 3). The ECMWF Reanalysis v5 (ERA5) data are available from https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5 (accessed on 1 June 2021).