Generating a Long-Term Spatiotemporally Continuous Melt Pond Fraction Dataset for Arctic Sea Ice Using an Artificial Neural Network and a Statistical-Based Temporal Filter

Abstract

The melt pond fraction (MPF) is an important geophysical parameter of climate and the surface energy budget, and many MPF datasets have been generated from satellite observations. However, the reliability of these datasets suffers from short temporal spans and data gaps. To improve the temporal span and spatiotemporal continuity, we generated a long-term spatiotemporally continuous MPF dataset for Arctic sea ice, which is called the Northeast Normal University-melt pond fraction (NENU-MPF), from Moderate Resolution Imaging Spectroradiometer (MODIS) data. First, the non-linear relationship between the MODIS reflectance/geometries and the MPF was constructed using a genetic algorithm optimized back-propagation neural network (GA-BPNN) model. Then, the data gaps were filled and smoothed using a statistical-based temporal filter. The results show that the GA-BPNN model can provide accurate estimations of the MPF (R² = 0.76, root mean square error (RMSE) = 0.05) and that the data gaps can be efficiently filled by the statistical-based temporal filter (RMSE = 0.047; bias = −0.022). The newly generated NENU-MPF dataset is consistent with the validation data and with published MPF datasets. Moreover, it has a longer temporal span and is much more spatiotemporally continuous; thus, it improves our knowledge of the long-term dynamics of the MPF over Arctic sea ice surfaces.

1. Introduction

The melt pond fraction (MPF) of Arctic sea ice is an important geophysical parameter for the climate and surface energy budget of the cryosphere [1,2,3,4], and it has been widely used in sea ice evolution and general circulation models (GCMs) [5,6]. During the summer in the Arctic, melt ponds appear when the snow and/or sea ice on the top of sea ice surface melt [7,8]. These ponds can achieve coverage of up to 50–60% over the entire Arctic sea ice zone [9,10,11,12,13], and up to 75% coverage of first-year ice surfaces [6,11,14].

The formation and disappearance of melt ponds are closely related to changes in the sea ice albedo [10,15,16], and the expansion of melt ponds during summer can significantly decrease the surface albedo of Arctic sea ice [17]. This decrease in the surface albedo can lead to more absorption of solar radiation, which can further amplify the climate change in the Arctic via the sea-ice feedback mechanism [18]. It has also been reported that the MPF is considered to be a good predictor for the Arctic sea ice extent in September and can be used to improve forecasts of seasonal sea ice changes and climate models [19,20,21,22,23].

Therefore, it is crucial to monitor the spatiotemporal dynamics and evolution of the MPF over the Arctic sea ice surface. In recent decades, the MPF of Arctic sea ice has typically been measured through ship-based visual observations, aerial photography, and satellite remote sensing [24,25,26,27,28,29,30]. However, owing to the limited spatial coverage of surface measurements, the spatiotemporal variations in the MPF cannot be represented well by ground measurements [24,31,32,33,34], such as the data collected by floating buoys [31], fixed ice stations [33], and cruise ships [34,35].

Therefore, high-spatial-resolution images obtained from tethered balloons [25], unmanned aerial vehicles (UAVs) [36], helicopters [28,37], planes [38], and satellite platforms [26,39,40,41,42,43] are more commonly used to derive the MPF of Arctic sea ice. The various classification methods were developed to estimate the MPF of Arctic sea ice based on various types of high-spatial-resolution images, e.g., aerial photographs [37,44,45], Landsat [39,46], Sentinel-2 [43], QuickBird [42], and WorldView-2 [37,41] data.

However, owing to their limited spatial coverage and relatively long revisit cycles, these high-spatial-resolution data cannot well represent the spatiotemporal variations in the MPF over the entire Arctic sea ice area [38,40,45,47,48,49]. Therefore, moderate-spatial-resolution satellite data with much more frequent observations, such as Moderate Resolution Imaging Spectroradiometer (MODIS) and Medium Spectral Resolution Imaging Spectrometer (MERIS) data, are usually considered to be more suitable for generating long-term MPF datasets for Arctic sea ice [50,51,52,53].

Recently, many MPF datasets have been generated from data obtained by multiple remote sensing platforms and sensors. Tschudi et al. [50] proposed a linear spectral unmixing method to derive the fractional coverages of different surface types (e.g., white ice, snow-covered ice, melt ponds, and open ocean water) from MODIS data, and their results revealed that the estimated MPF in the Beaufort/Chukchi Sea region during the summer of 2004 is in good agreement with the measurements acquired by aerosonde unpiloted aerial vehicles and a surface campaign.

To reduce the computation cost of the linear spectral unmixing method, an artificial neural network method was used to speed up the procedure for estimating the MPF [54], and a MPF dataset over the Arctic sea ice surface from 2000 to 2011 was derived by researchers at the University of Hamburg (hereinafter referred to as the UH-MPF dataset) [55]. Although the linear spectral unmixing method can be used to estimate the MPF, it results in large estimation uncertainty due to the prior fixed spectral reflectance of the sea ice components, especially when the MPF and melt pond depth change significantly [56].

To overcome this issue, Zege et al. [56] proposed a melt pond detector (MPD) algorithm to retrieve the MPF of Arctic sea ice from MERIS data. An asymptotic radiative transfer theory (AART) model [4,57,58] was developed to simulate the bidirectional reflectance factor (BRF), black-sky, and white-sky albedo over the surfaces of snow, white ice, and melt ponds. Then, the MPF and albedo over the Arctic region were accurately retrieved from the MERIS data using an iterative procedure.

Finally, a daily MPF dataset for Arctic sea ice from 2002 to 2011 was generated from MERIS Level 1B data by researchers at the University of Bremen (hereinafter referred to as the UB-MPF dataset) [51,52,59]. As MERIS data were not obtained after 2011, the researchers at the University of Bremen also used the data acquired by the Ocean and Land Color Instrument (OLCI) onboard Sentinel-3 to generate an MPF dataset for Arctic sea ice (hereinafter referred to as the UB-OLCI dataset) from 2017 to the present based on a procedure similar to the MPD [60].

In addition, the MPF can also be estimated by establishing the relationship between the satellite observations and MPF data collected from multiple sources. Ding et al. [22] proposed an ensemble-based deep neural network method for estimating the MPF and released a long-term MPF dataset for Arctic sea ice from 2000 to 2019 (as it was derived by researchers at Beijing Normal University; it is hereinafter referred to as the BNU-MPF dataset).

Although many MPF datasets for Arctic sea ice have been generated from a variety of satellite observations, the temporal span and spatiotemporal continuity of these datasets still need to be improved. Among these datasets, the UB-MPF dataset was considered to have a higher estimation accuracy because it was derived using a physical-based model [51,56]. However, owing to its relatively short temporal span (2002–2011), the UB-MPF dataset is not suitable for analyzing the long-term trends and evolutions of the MPF. In addition, most of these MPF datasets are spatiotemporally discontinuous due to cloud obscuration, and the analysis results based on the datasets with numerous data gaps are also questionable.

To extend the temporal span and improve spatiotemporal continuity of the MPF datasets, we generated a long-term spatiotemporally continuous MPF dataset for Arctic sea ice from 2000 to 2020, which is called the Northeast Normal University-melt pond fraction (hereinafter referred to as the NENU-MPF) dataset, using an artificial neural network and a statistical-based temporal filter. In the following sections of this article, we present the data and methodologies used in this study in Section 2. We describe the validation and compare the results with satellite-based and ship-based measurements as well as the published MPF datasets in Section 3. We provide a short summary of this study in Section 4.

2. Materials and Methods

2.1. Data

2.1.1. MODIS Daily Surface Reflectance Data

In this study, we used the MODIS Terra daily surface reflectance product (MOD09GA, Version 6) to generate an MPF dataset for Arctic sea ice from 2000 to 2020. The MOD09GA data provide daily atmospheric corrected surface reflectance data for MODIS bands 1–7 (the central wavelengths of bands 1–7 are 648, 859, 466, 554, 1244, 1631, and 2119 nm, respectively), ancillary incident/viewing geometries (e.g., the solar zenith angle (SZA), view zenith angle (VZA), solar azimuth angle (SAA), view azimuth angle (VAA)), and quality control (QC) information with spatial resolutions of 500 m and 1 km. The MOD09GA data over the Arctic region from 2000 to 2020 were downloaded from the Earthdata website (https://www.earthdata.nasa.gov, accessed on 1 October 2021).

Then, the pixels covered with clouds and open ocean water were screened out based the QC information and sea ice concentration data, and only the cloud-free surface reflectance data over the Arctic sea ice zone were retained. Finally, the surface reflectance data were further reprojected from the sinusoidal projection to the polar stereographic projection and were aggregated to a spatial resolution of 12.5 km. In the subsequent training and prediction processes, only the cloud-free surface reflectance data over the Arctic sea ice zone from 8 May to 24 September of each year were employed.

2.1.2. MPF Data

We used the currently available MPF datasets as the training data for the neural network model and as intercomparison data (

Table 1. MPF datasets for Arctic sea ice derived from satellite observations.

Name	Sensor	Temporal Span	Temporal Resolution (d)	Spatial Resolution (km)	References
UH-MPF	MODIS	2000–2011	8	12.5	[54,55]
UB-MPF	MERIS	2002–2011	1	12.5	[51,52,56]
UB-OLCI	OLCI	2017 to the present	1	12.5	[60]
BNU-MPF	MODIS	2000–2019	8	12.5	[22,61]

). The UB-MPF dataset for 2003 was used as training data for the artificial neural network; and the UB-MPF datasets for other years (from 2004 to 2011), the UB-OLCI dataset, the UH-MPF dataset, and the BNU-MPF dataset were used for comparison with our results. In this study, all of the MPF datasets were reprojected to the polar stereographic projection with a spatial resolution of 12.5 km for comparison. In the comparison, the NENU-MPF dataset was directly compared with the UB-MPF and UB-OLCI datasets; however, it was aggregated to an 8-day composite temporal resolution for comparison with the UH-MPF and BNU-MPF datasets first.

2.1.3. Validation Data

We used two satellite-based and three ship-based MPF measurements to validate our results, and the spatial distributions of these measurements are shown in Figure 1. The Arctic sea ice melt pond statistics and maps (ASIMPSM) dataset was obtained from high-spatial-resolution (1 m) satellite images using the maximum likelihood supervised classification method [62]. The ASIMPSM datasets were collected in the Beaufort Sea, East Siberian Sea, Canadian Arctic, and Fram Strait during 1999–2001, with scene sizes of 3 to 10 km.

The measurements of Earth data for environmental analysis (MEDEA) dataset was obtained from high-spatial-resolution (1 m) satellite images based on the supervised classification method [63]. The MEDEA dataset was collected at 69–86°N during 1999–2014 for sea ice model validation, with scene sizes of 5 to 25 km. In this study, the original images of the ASIMPSM and MEDEA were aggregated to 12.5 km for comparison with our results.

The IceWatch dataset was compiled from archival visual sea ice observations made by ships in the Northern Hemisphere. The Arctic Shipboard Sea Ice Standardization Tool (ASSIST) was used to document the sea ice conditions (e.g., the ice type, sea ice concentration, snow depth, MPF, and melt pond depth). In this study, the observations collected at 69–89°N in July to August 2011, July 2017, and September 2020 were employed for validation.

The Trans-Arctic survey of the Arctic Ocean in transition (TransArc) dataset was derived from hourly visual observations of sea ice conditions (e.g., sea ice type, sea ice thickness, sea ice concentration, and MPF) acquired by the Polarstern icebreaker from August to September 2011 [34]. The Healy Oden Trans-Arctic Expedition (HOTRAX) dataset was compiled from measurements of ice station, ship, and helicopter aerial images acquired every two hours [9]. In this study, the HOTRAX data acquired at 69.55–89.98°N from August to September 2005 was used for validation.

2.2. Methods

In this study, a genetic algorithm optimized back-propagation neural network (GA-BPNN) and a statistical-based temporal filter were used to generate the NENU-MPF dataset. The main workflow for generating the NENU-MPF dataset is shown in Figure 2.

2.2.1. GA-BPNN Model

In this framework, the GA-BPNN model was used to determine the non-linear relationship between the MODIS directional reflectance/geometries and the MPF derived from the MERIS data (UB-MPF dataset). The daily UB-MPF dataset from May to September 2003 and the corresponding MODIS surface reflectance data were used in the GA-BPNN model, and 70%, 15%, and 15% of the pixels were used as training, validating, and testing data, respectively. The input variables were the surface directional reflectance of the seven MODIS bands and four incident/viewing geometries (SZA, VZA, SAA, and VAA).

The objective output variable was the MPF of the Arctic sea ice. It has been widely reported that the snow, melt ponds, and ocean water are non-Lambertian, and the directional surface reflectance of Arctic sea ice varied significantly with the incident/viewing geometries [4,64]. To capture the reflectance anisotropic properties, the incident/viewing geometries were used as input variables for establishing the relationship between surface reflectance and shortwave surface albedo of Arctic sea ice using an ensemble back-propagation neural network [65].

In this study, a similar approach was carried out for predicting the MPF of Arctic sea ice, and the incident/viewing geometries were also used as input variables. The framework and structure of the GA-BPNN model are shown in Figure 3. As the relationship between the directional surface reflectance/geometries and MPF varied during the different stages of sea ice melting, five GA-BPNN models were separately constructed for each month (from May to September) in this study.

To obtain the optimal structure of the artificial neural network, the following experiments were conducted to determine the structure of the GA-BPNN model. First, the number of hidden layers was assumed to be one, and the influence of the total number of neurons on the prediction accuracy was analyzed. To minimize the influence of the initial values of the weights, we ran the BP neural network model 100 times using random weights and biases assigned for each epoch. Then, the maximum number of hidden layers was assumed to be six, and the influence of the neural network’s structure on the prediction accuracy was examined.

The results of the experiment show that the minimum value of the median root mean square error (RMSE) occurred when the total number of neurons was 22 (the neural network model for May; Figure 4). When the total number of neurons was 22, 20 hidden layer structures could be obtained: (22), (11, 11), (6, 6, 10), (6, 10, 6), (10, 6, 6), (5, 5, 5, 7), (5, 5, 7, 5), (5, 7, 5, 5), (7, 5, 5, 5), (4, 4, 4, 4, 6), (4, 4, 4, 6, 4), (4, 4, 6, 4, 4), (4, 6, 4, 4, 4), (6, 4, 4, 4, 4), (7, 3, 3, 3, 3, 3), (3, 7, 3, 3, 3, 3), (3, 3, 7, 3, 3, 3), (3, 3, 3, 7, 3, 3), (3, 3, 3, 3, 7, 3), and (3, 3, 3, 3, 3, 7).

The results revealed that the minimum value of the median RMSE occurred when the hidden layer structure was (10, 6, 6) (the neural network model for May, Figure 5). Finally, the optimal number of total neurons for the neural networks for May to September were determined to be 22, 26, 25, 23, and 24, respectively, and the optimal structures of the neural networks for May to September were determined to be (10, 6, 6), (13, 13), (25), (8, 5, 5, 5), and (12, 12), respectively.

Moreover, a genetic algorithm [66] was also used to optimize the values of the weights and biases in the BP neural network model and to minimize the prediction errors due to the random assignment of weights and biases. The results show that the performances of the BPNN models were significantly improved by the genetic algorithm, and the determination coefficients (R²) of the GA-BPNN models (0.46, 0.73, 0.68, 0.62, and 0.47 for May to September, respectively) were higher than those of the BPNN models (0.42, 0.72, 0.64, 0.6, and 0.44 for May to September, respectively) (Figure 6).

2.2.2. Statistical-Based Temporal Filter

To generate a spatiotemporally continuous MPF dataset for Arctic sea ice, a statistical-based temporal filter was utilized to fill in the data gaps in the GA-BPNN predictions. In our previous study, a statistical-based temporal filter was used to generate spatiotemporally continuous datasets of the land surface albedo [67,68]. The daily land surface albedo can be reconstructed based on a weighted combination of the albedo estimated from remote sensing data, the albedo estimated based on the relationships between neighboring days, and the albedo climatology (multi-year mean albedo).

However, because the MPF of Arctic is highly spatiotemporally dynamic and has large inter-annual variability, the MPF climatology was not used in this study, and the missing data were mainly filled based on the relationships between neighboring days. It is clear that the MPFs of neighboring days are closely correlated, and this can be explicitly expressed as a linear relationship [67]:

$$MP{F}_k=a_{∆k}MP{F}_{k+∆k}+b_{∆k}+{\xi }_{∆k},$$

(1)

where $MP{F}_k$ is the MPF on k^th day, $MP{F}_{k+∆k}$ is the MPF at (k + ∆k)^th day, ∆k is the temporal interval between the neighboring days (−4 ≤ ∆k ≤ 4 for this study), $a_{∆k}$and $b_{∆k}$ are the regression coefficients of the linear regression, and ${\xi }_{∆k}$ is the regression residuals. These parameters can be derived based on the priori statistics of the MPF climatology [67]:

$$a_{∆k}={\rho }_{∆k}\frac{{\sigma }_{∆k}}{{\sigma }_{k+∆k}},$$

(2)

$$b_{∆k}={\mu }_k-a_{∆k}{\mu }_{k+∆k},$$

(3)

$${\xi }_{∆k}=\sqrt{\left(1-{\rho }_{∆k}^2\right)}{\sigma }_k,$$

(4)

where ${\mu }_k$(k = 1, 2, …, 140, k denotes the sequential days from 8 May to 24 September) and${\mu }_{k+∆k}$ are the mean MPFs of the multi-year data on k^th and (k + ∆k)^th days, ${\sigma }_k$ and ${\sigma }_{k+∆k}$ are the standard deviations of the MPF on k^th and (k + ∆k)^th days, ${\rho }_{∆k}$($∆k$= −4, −3, …, 3, 4) is correlation coefficient for the MPFs of the neighboring days with an interval of ∆k, and ${\xi }_{∆k}$ is the RMSE between the MPFs of k^th and (k + ∆k)^th days. The above-mentioned priori statistics of the MPF climatology can be derived from statistical analysis of the UB-MPF dataset.

Finally, the MPF of k^th day can be reconstructed through the weighted combination of MPFs estimated using the GA-BPNN model and the MPFs predicted by the statistical relationships between the neighboring days:

$$\widehat{MP{F}_k}=\frac{{{\displaystyle \sum }}_{∆k=-4}^{∆k=+4}\frac{a_{∆k}MP{F}_{k+∆k}+b_{∆k}}{{\xi }_{∆k}^2+a_{∆k}^2{\epsilon }_{k+∆k}^2}+\frac{MP{F}_k}{{\epsilon }_k^2}}{{{\displaystyle \sum }}_{∆k=-4}^{∆k=+4}\frac{1}{{\xi }_{∆k}^2+a_{∆k}^2{\epsilon }_{k+∆k}^2}+\frac{1}{{\epsilon }_k^2}},$$

(5)

where $\widehat{MP{F}_k}$is the reconstructed MPF on k^th day, ${\xi }_{∆k}$is the RMSE between the MPFs on k^th and (k + ∆k)^th days, and ${\epsilon }_k$ and ${\epsilon }_{k+∆k}$ are the standard derivations of the MPFs estimated using the GA-BPNN model for k^th and (k + ∆k)^th days.

3. Results and Discussion

3.1. Accuracy of the GA-BPNN Model

The sophisticated non-linear relationship between the directional surface reflectance (seven bands of MODIS data), incident/viewing geometries (SZA, VZA, SAA, and VAA), and MPF can be efficiently constructed using the GA-BPNN model. The validation results show that the GA-BPNN models have good performance in terms of predicting the MPF of Arctic sea ice, with an R² of 0.76 and an RMSE of 0.05 (Figure 7). The predicted RMSEs of the GA-BPNN models for May to September are 0.032, 0.054, 0.059, 0.048, and 0.049, respectively (Figure 8). The GA-BPNN predictions are in good agreement with the validation dataset (i.e., the true MPF) in June, July, and August, and most of the estimates plot close to the 1-to-1 line.

The validation results exhibit relatively lower R² values for May and September, which is mainly because the MPFs in May and September are much lower than those in June, July, and August, and larger uncertainties in the MPF can be obtained when the MPF is lower.

3.2. Accuracy of Statistical-Based Temporal Filter

The MPF datasets for Arctic sea ice before and after gap filling are shown in Figure 9 and Figure 10, respectively. The original dataset generated using the GA-BPNN model is not spatiotemporally continuous owing to cloud obscuration. However, the data gaps can be efficiently filled by the statistical-based temporal filter, and the generated NENU-MPF dataset represents the spatiotemporal dynamics of MPF of the Arctic sea ice well.

From May to early June, the MPFs were relatively low, and the melt ponds mainly appeared in the northern coastal areas of Asia and North America. Then, from mid-June to mid-July, the MPFs increased significantly in the first-year ice zone (Beaufort Sea, East Siberian Sea, and Laptev Sea) as the air temperature increased. Finally, from late July to September, the MPFs decreased significantly and reverted to the level during the pre-melting season.

An example of MPF data processed using the statistical-based temporal filter is shown in Figure 11. As can be seen from Figure 11, numerous data gaps existed in the original MPF dataset derived using the GA-BPNN model; however, a spatiotemporally continuous dataset without missing data was obtained by applying the statistical-based temporal filter.

Moreover, the temporal variations in the NENU-MPF dataset are similar to the reference UB-MPF dataset, which indicates that the MPF data missing due to cloud obscuration can be well reconstructed using the statistical-based temporal filter. Figure 12 presents a comparison of the gap-filled MPF with the reference UB-MPF dataset. The comparison results show that the gap-filled MPF derived using the statistical-based temporal filter is consistent with the UB-MPF dataset, with an R² of 0.43, an RMSE of 0.047, and a bias of −0.022. Its accuracy is lower than that of the GA-BPNN model. However, with the advantage of generating spatiotemporally continuous dataset, the reconstruction accuracy of the statistical-based temporal filter is acceptable at most times.

3.3. Validation Results

The validation results of the NENU-MPF dataset with the ASIMPSM dataset are shown in Figure 13. The NENU-MPF dataset is roughly consistent with the MPF derived from the images collected in the Beaufort Sea, East Siberian Sea, Canadian Arctic, and Fram Strait during 2000–2001, with an RMSE of 0.08 and an R² of 0.58. The MEDEA, IceWatch, TransArc, and HOTRAX datasets were used to validate the NENU-MPF dataset, and the UB-MPF dataset was also used as a reference for comparison. The validation results show that the NENU-MPF dataset is roughly consistent with the MEDEA (R² = 0.41, RMSE = 0.102, and bias = −0.04), IceWatch (R² = 0.60, RMSE = 0.105, and bias = 0.05), TransArc (R² = 0.027, RMSE = 0.21, and bias = 0.174), and HOTRAX (R² = 0.20, RMSE = 0.098, and bias = 0.007) datasets (Figure 14).

The results of the validation against the UB-MPF dataset indicate that it has an accuracy comparable with that of the NENU-MPF dataset and is roughly consistent with the MEDEA (R² = 0.39, RMSE = 0.119, and bias = −0.07), IceWatch (R² = 0.26, RMSE = 0.126, and bias = 0.002), TransArc (R² = 0.007, RMSE = 0.219, and bias = 0.177), and HOTRAX (R² = 0.004, RMSE = 0.149, and bias = 0.106) datasets. The satellite-based measurements (ASIMPSM and MEDEA datasets) are more consistent with our results than with the ship-based measurements (IceWatch, TransArc, and HOTRAX datasets). We also found that the TransArc dataset is not consistent with the NENU-MPF and UB-MPF datasets, which is mainly due to the discrepancy between the spatial resolutions of the remote sensing datasets (12.5 km) and the footprints of the ship-based measurements (50 m to 10 km, which mainly depends on the horizontal visibility and height of the observation position).

3.4. Comparison Results

The results of the comparisons of the NENU-MPF dataset with the UB-MPF, BNU-MPF, UB-OLCI, and UH-MPF datasets are shown in Figure 15. The results show that the NENU-MPF dataset is consistent with the UB-MPF dataset, with an R² of 0.70, an RMSE of 0.052, and a bias of −0.02, which is reasonable because the UB-MPF dataset was used as the training dataset for the GA-BPNN model. The NENU-MPF dataset is also consistent with the BNU-MPF dataset, with an R² of 0.40, an RMSE of 0.06, and a bias of −0.011.

While the consistencies between the NENU-MPF and UB-OLCI datasets (R² = 0.27, RMSE = 0.114, and bias = −0.064) and the NENU-MPF and UH-MPF datasets (R² = 0.27, RMSE = 0.09, and bias = −0.025) are relatively poor. The discrepancies between the NENU-MPF and UH-MPF datasets are mainly due to the different algorithms used to estimate the MPF. Moreover, the spectral unmixing method used in the UH-MPF dataset assumes that the spectral reflectance of the sea ice components is invariant [54], which may result in a large error when the depth and area of the melt ponds increase.

The annual mean MPF values derived using the UB-MPF, UH-MPF, NENU-MPF, and BNU-MPF datasets from 2002 to 2020 are shown in Figure 16. As there are large data gaps in the UH-MPF, BNU-MPF, and UB-MPF during 2000–2001, the annual mean MPF values of the Arctic sea ice during 2000–2001 are not shown in Figure 16.

To eliminate the missing data, the original UB-MPF dataset was filtered using a running mean of 8 days before the annual mean MPF was calculated. To minimize the influence of the missing data, only the pixels with valid MPF values were used. In addition, the UB-OLCI dataset was not compared with the other datasets due to the limited spatial coverage of this dataset. The results show that the NENU-MPF, UB-MPF, BNU-MPF, and UH-MPF exhibit MPF trends of −0.0015/yr, −0.00005/yr, –0.0014/yr, and 0.0018/yr, respectively.

The NENU-MPF dataset has a much longer temporal span than the UB-MPF and UH-MPF datasets, and its decreasing trend is consistent with that of the BNU-MPF dataset from 2002 to 2020. It is noticeable that the values of the UH-MPF dataset are much higher than those of the other datasets, and it exhibits a contrast increasing MPF trend, which is not consistent with the MPF trends derived from the UB-MPF, NENU-MPF, and BNU-MPF datasets. The discrepancies between these datasets can be explained by the differences in the estimation methods used, as well as the degree of spatiotemporal continuity (or data gaps) of these datasets.

The multi-year mean temporal evolutions of the MPF derived from the UH-MPF, UB-MPF, NENU-MPF, and BNU-MPF datasets from May to September during 2003–2011 are shown in Figure 17. For all of the MPF datasets, the MPF initially increases and then decreases with time. The MPF of the Arctic sea ice is lowest in early May, rapidly increases in June, and reaches the maximum MPF value in late June to mid-July. Then, the MPF decreases in August and September, and it falls to the lowest value in September. These datasets exhibit similar temporal trends in the MPF, and the UH-MPF and NENU-MPF datasets present the highest values in the peak season, followed by the BNU-MPF and UB-MPF.

In addition, the UH-MPF dataset exhibits an approximate one-month lag in the decline of the MPF, which is not consistent with the other three datasets. The temporal variations in the MPF of the Arctic sea ice derived from the NENU-MPF datasets from 2000 to 2020 are shown in Figure 18. The regional mean MPF varies from approximate 6% to 20% exhibits a large standard derivation from late June to early August. The newly generated NENU-MPF dataset accurately characterizes the temporal evolution of the MPF and captures the phenomena of abnormal high MPF values in 2007.

3.5. Discussion

The temporal span and spatiotemporal continuity are crucial characteristics of MPF datasets for Arctic sea ice, which influence the reliability of the long-term trends derived from satellite observations. As most of the geophysical parameters of the Earth vary with climatic and environmental factors, the trends of the geophysical parameters can be significantly influenced by the general circulation and atmospheric oscillations [69], as well as the data gaps caused by cloud obscuration [67], and thus datasets with longer temporal span and more spatiotemporal continuity are required to more accurately determine the MPF trend.

The comparison results indicate that the newly generated NENU-MPF dataset has a much longer temporal span (from 2000 to 2020) than that of UB-MPF (from 2002 to 2011) and UH-MPF (from 2000 to 2011) datasets and is more spatiotemporally continuous compared to the other datasets. Thus, it has potential advantages in term of representing the dynamics and evolution of the MPF of Arctic sea ice more realistically.

For example, the long-term trend in MPF of the Arctic sea ice derived from the NENU-MPF dataset can be considered more accurate and reliable, and the NENU-MPF dataset can be used as a reference for calibrating and developing the melt ponds evolution models [70]. In addition, there are also many potential applications for the newly generated MPF dataset of Arctic sea ice [50], such as improving the treatment of sea ice albedo in the global climate models, qualifying the uncertainty of sea ice concentration dataset, analyzing the interannual variability in melt processes of Arctic sea ice, and estimating the contribution of melt ponds to the sea ice albedo feedback mechanism.

Although differences in the MPF values of the different datasets were identified, a significant decreasing MPF trend can be found in most of these datasets (Figure 15). The large bias and opposite trend of the UH-MPF dataset are mainly due to the assumption in the spectral unmixing method that the spectral reflectance of the sea ice components is invariant [54], which resulted in large uncertainties in estimating the MPF when the sea ice surface changed significantly during the melting season.

In addition, the data gaps due to cloud obscuration also increased the uncertainty of the MPF trends derived from these datasets. The significant decreasing trend of the MPF in the Arctic can be explained by the decline in the sea ice extent, and the fact that the melted area of first-year ice (where the melt ponds mainly appear) has increased dramatically due to global warming, leading to a decrease in the MPF of the Arctic sea ice in recent decades [70].

In this study, a GA-BPNN model and statistical temporal filter were applied to improve the temporal span and spatiotemporal continuity of the MPF dataset. Compared to previous studies that did not consider the effect of the reflectance anisotropic of the melt ponds, in this study, the bidirectional reflectance distribution (BRDF) characteristics of the melt ponds were fully considered using the directional reflectance/geometries and GA-BPNN model. Our knowledge of the reflectance anisotropic properties of the melt ponds was derived from the physical-based UB-MPF dataset. In addition, the missing data due to cloud obscuration were efficiently filled using a statistical-based temporal filter, and the spatiotemporal continuity of the MPF dataset was significantly improved.

The shortcomings of this study and the issues that need to be addressed in the future are listed as follows:

(1) The spatial resolution and validation accuracy of the proposed NENU-MPF dataset need to be improved in the future. Owing to the inconsistency of the footprints of the MPF dataset and the validation data, our results are only roughly consistent with the validation dataset. To address this issue, a labeled dataset for melt ponds in the Arctic derived from high-spatial-resolution satellite imageries is needed [71]. There is a potential to generate a more accurate MPF dataset using deep learning and a training dataset, and the validation of the MPF datasets would also benefit from a labeled dataset of melt ponds. In addition, developing new validation and uncertainty qualification methods [72,73] over Arctic sea ice region are also required in future studies.

(2) Although data gaps can be efficiently filled using the statistical-based temporal filter, the gap filling methods with higher accuracy are still required for generating spatiotemporally continuous MPF datasets. In this study, only the relationship between the neighboring days was used in the gap filling of the MPF dataset, and this assumption may be invalid due to snow fall/melt, and sea ice drifts, which would result in a relatively poor gap filling accuracy under certain circumstances. In the future, the temporal correlation, spatial autocorrelation, and the relationship with the passive remote-sensing data can be incorporate to generate a seamless MPF dataset [74]. In addition, the deep-learning and inpainting methods can also be used for developing new gap filling methods, which have shown promising results in recent studies [75,76].

(3) Satellite observations acquired by multiple platforms and sensors can be jointly used to extend the temporal span of the MPF dataset. The method proposed in this study can be adapted to other remote sensing data, e.g., Advanced Very High Resolution Radiometer (AVHRR) and Visible Infrared Imaging Radiometer Suite (VIIRS) data. The combined use of AVHRR, MODIS, and VIIRS data has the potential to extend the temporal span to more than 40 years (from 1981 to the present).

4. Conclusions

In this study, we generated a long-term spatiotemporally continuous MPF dataset for Arctic sea ice from MODIS data from 2000 to 2020. First, a non-linear relationship between the MODIS directional reflectance/geometries and the UB-MPF dataset was successfully constructed using a GA-BPNN model. Then, the missing data due to cloud obscuration were efficiently filled using a statistical-based temporal filter. Finally, the generated NENU-MPF dataset was validated using satellite-based and ship-based measurements and was compared with other published datasets, such as the UB-MPF, UB-OLCI, UH-MPF, and BNU-MPF datasets.

The MPF values of Arctic sea ice were accurately predicted using the GA-BPNN model (R² = 0.76, RMSE = 0.05), and the missing data were filled in using the statistical-based temporal filter. The validation and comparison results showed that the NENU-MPF dataset was consistent with the validation measurements as well as the UB-MPF (R² = 0.70, RMSE = 0.052) and BNU-MPF (R² = 0.40, RMSE = 0.06) datasets. The newly generated NENU-MPF dataset was demonstrated to be more spatiotemporally continuous with a longer temporal span (21 years).

The results of this study improve our knowledge of the long-term variations in Arctic sea ice over recent decades. The newly generated NENU-MPF dataset demonstrated great advantages in providing accurate and reliable trends and interannual variability in the MPF of Arctic sea ice during the recent two decades and can also be used as a fundamental dataset in studies of global climate change and the surface energy budget.