STUDY ON ARCTIC MELT POND FRACTION RETRIEVAL ALGORITHM USING MODIS DATA

During spring and summer, melt ponds appear on the sea ice surface in the Arctic and play an important role in sea ice-albedo feedback effect. The melt pond fraction (MPF) can be retrieved using multi-band linear equations, but the calculation is complicated by the ill-conditioned reflectance matrix. In this paper, we calculated the condition numbers which represent the degree of the illconditioned reflectance matrix in the results of the MPF from a MODIS-based unmixing algorithm. The condition number is introduced here as a criterion for the sensitivity of the solution in the system to the error in the input value. By combining 3 bands among 5 visible and near-infrared bands of MODIS data, the results show that the three-band combination with the lowest sensitivity to the error of input is B245. To improve the algorithm, we introduce pre-processing to remove open water from the four surface types and then remove one reflectance equation from the original set. The best two-band combination algorithm is B15. Compared with the discrimination results from Landsat5-TM, the RMS is 0.14. This algorithm is applied in pan-Arctic scale, the MPF results are larger than that from University of Hamburg, especially in the Pacific sector.


INTRODUCTION
The albedo of melt ponds is between open water and sea ice. Melt ponds occupy a large fraction of the Arctic sea ice surface during spring and summer. The fraction and distribution of melt ponds play an important role in the sea ice-albedo feedback. Accurate melt pond fraction (MPF) data is critical for calculating solar radiation absorption in sea ice and ocean underneath.
Optical imagery is one of the primary sources to acquire MPF information. Observations have shown that the melt pond fraction can vary by more than 60% throughout the melt season and by up to 40% depending on years and locations (Polashenski et al., 2012;Landy et al., 2014;. The Arctic melt ponds are dominated by small-sized ones (normally <200m 2 ) (Perovich et al., 2002;Tschudi et al., 2001;Lu et al., 2010). If the sizes of several melt ponds are larger than the pixel area of high-resolution optical imagery, a basic principle of retrieving MPF is used to discriminate pixels of snow/bare ice from open water and melt ponds, based on their different reflective properties (Markus et al., 2002;Rösel et al., 2011). For the medium-resolution sensors with 500m or 1km resolution, such as MODerate-resolution Imaging Spectroradiometer (MODIS), unmixing is certainly required. Although the calculation of the unmixing formula has problems (for example, when the equations are approximately linearly correlated to each other, ill-conditioned reflectance matrix will occur and the solutions are sensitive to errors in the input data), its advantage of broad coverage and daily repeat-frequency gives people the reason to use it for the research of Arctic variability and to improve the algorithm constantly.
The wildly used MODIS MPF product (Rösel et al., 2012, hereafter called R12) is based on a spectral unmixing procedure which first proposed by Tschudi et al. (2008), hereafter called T08. In the R12 algorithm, they applied a neural network to determine the coefficients of algorithm. Zege et al. (2015) developed a MPF algorithm using the analytical solution of optical thickness layer to describe the albedo of white ice to calculate the BRDF of white ice and melt ponds based on MEdium-Resolution Imaging Spectrometer (MERIS) data. This algorithm is more physically grounded than that of R12. However, their MPF results only cover part of the Arctic .
In this paper we calculated the condition numbers which represent the degree of the ill-conditioned reflectance matrix happened in the results of the MPF based on T08 using MODIS data. Some improvements were made to reduce the sensitivity of the linear functions. And the results were compared with Landsat data as well as the MPF product from Hamburg University (Rösel et al., 2012).

DATA
MODIS surface reflectance data (MOD09 and MOD09A1) is used as basic dataset. MOD09 is L2 data and stored by strips. MOD09 provides MODIS surface reflectance for bands 1 and 2 (at 250m), bands 1-7 (at 500 m), bands 1-16 (at 1 km resolution) and geographic information. MOD09 TrueColor image ( Figure  1) which covered Beaufort Sea on June 13, 2004 was used following the original case of T08. MOD09A1 is L3 grid data. Each MOD09A1 pixel includes the best possible L2G observation during an 8-day period as selected on the basis of high observation coverage, low view angle. MOD09A1 contains bands 1-7 at 500-meter resolution with quality control. It provides the broad coverage of the whole Arctic data and has removed the strong effects of cloud, shadow, and aerosol compared with MOD09GA data, which is a grid track data.
The MPF product from University Hamburg (R12) is used to compare with these results in the year of 2007. Landsat5-TM data archived by the United States Geological Survey (USGS) serves as the high-resolution data for validating MODIS products.

METHOD
Th spatial resolution of MODIS data is 500m, much greater than the common melt pond size. Therefore, it is impossible to resolve individual pond directly in almost all cases. T08 algorithm assumes every surface type in a given pixel reflects independently, and the total energy they reflect at a given wavelength determines the pixel's corresponding reflectance. For each pixel, (1) where ai = the coverage of each surface type ri = each surface type's reflectance R = the "in-situ" reflectance from MOD09 data i = indicates surface types k = MODIS bands In the original algorithm design, there are 4 types: snow, ice, melt pond, and water. This means 4 independent linear equations are needed to solve for coverage. The sum of all the types' coverage should be 1, The other 3 equations are acquired from bands 1, 2, 3 and the reflectance coefficients ri are from in-situ observations. According to T08, the reflectance coefficient matrix is as following: Where columns represent melt pond, ice, snow and water (from left to right), and rows denote reflectance coefficients in bands 1, 2, and 3, as well as the coefficients in the equation (2). The first 3 elements in R are reflectance values from MOD09 and the last one is 1.00, or the sum of all the types' coverage. The different observations report gave similar values (Grenfell and Gary A., 1977;Tucker et al., 1999;Tschudi et al., 2008) for the coefficient matrix. Finally, Lower-Upper (LU) decomposition is applied to solve the equation set.
In the T08 algorithm, the reflectance coefficients of four types in channels 1 and 3 are close, so the system of linear equations is approximately linearly correlated which causes the illconditioned problem to the solution of the reflectance matrix. More specifically, the solutions to the linear equations are sensitive to biases in the input data and result in large numeric errors and computational instability. To detect this problem, the condition number is introduced here as a criterion for the sensitivity of the solution in the system to the error in the input value.
The condition number of the matrix is equal to the product of the norm of the matrix and the norm of the inverse matrix. The smaller the condition number is, the lower the sensitivity of the matrix. By combining 3 bands among 5 visible and nearinfrared bands of MODIS data, their respective condition numbers were obtained. And the new band combination with the lowest sensitivity (lowest degree of linear correlation) will be selected.
To improve the algorithm, we introduce pre-processing and also remove one reflectance equation from the original set. And the same condition number process was taken afterwards by combining 2 bands among 5 MODIS bands.

Condition numbers for the three-band combination algorithms and the MPF retrieval results
To check the condition numbers of different combinations of three bands, the reflectivity of each surface type for different MODIS bands need to be given first. In   Norm∞  B123  284  140  209  B124  198  95  136  B125  129  81  169  B134  1355  785  1398  B135  161  109  193  B145  114  77  134  B234  555  283  388  B235  84  53  112  B245  74  46  97  B345  395  265  469  Table 2 there, which is physically impossible. Accordingly, the MPF algorithm using the above three combinations was applied to MOD09 data. The results show that B245 and B235 give more reasonable MPF than the original T08 combination of B312 (Figure 2b-d).

Condition numbers for the two-band combination algorithms and the MPF retrieval results
Even for the combination of B245, the conditional number is still between 46 and 97, so there still exists a linear correlation problem. To further reduce the sensitivity of the matrix, the next step is to reduce one equation for the experiment. This means we have to exclude one surface type to make the equation solvable.
Here we remove open water by ice-water discrimination based on the reflectance of band 4, which has the largest difference between the reflectance of water and that of other types.

Comparison and validation
According to the classification method of Markus et al. (2002Markus et al. ( , 2003 from Landsat-7ETM+, we discriminate different surface types based on bands 1, 2, and 3 reflectance of Landsat5 TM. r and 3 r stand for the reflectance of Landsat Band 1, Band 2 and Band 3, respectively. Figure 4a shows the combination TrueColor image of B123 from Landsat data. We projected Landsat data onto the MODIS data grid and statistically calculated the proportion of the areas covered by melt ponds by discriminating bare ice, snow, melt pond and water, and the MPF results of classification are shown in Figure 4b. There is an area with relatively high MPF in low left corner of the image. The comparisons were made between the retrieval MPF of five experiments in section 4.2 and the result from Landsat classification ( Figure 5). All the MPF results obtained by the two-band algorithms from MODIS data are smaller in high MPF regions, while relatively larger in low MPF regions than those obtained by surface type discrimination method from Landsat data.   (Figure 6). The main difference between the two algorithms is that our results have a relatively higher MPF in the regions covered by first-year ice, especially in the Pacific sector. Further comparison with aerial photographs and other high-resolution optical data is needed. Figure 5. MPF retrieval from two-band algorithms (left panels) and the difference compared with that of Landsat discrimination (right panels) Figure 6 Comparison between MPF from B15 algorithm and that from University Hamburg

CONCLUSION AND DISCUSSION
Melt ponds have an important effect on the decrease of Arctic sea ice in summer. Due to their low albedo, they accelerate the positive feedback process of sea ice-albedo. In this paper, we analyse a regional MPF algorithm proposed by T08 algorithm which is based on solving a linear equation set consisting of multi-channel reflectance and area faction. We developed the MPF retrieval algorithm with some experiments to reduce the sensitivity of the linear functions based on MODIS data. Comparison is made with Landsat solutions and the MPF product from University Hamburg. Several conclusions and discussions are summarized as the following.
The solution of linear unmixing equation is very sensitive to input data errors due to its ill-conditioned reflectance matrix. This issue is caused by the strong linear correlation between reflectance vectors of Band 1 and Band 3 in T08's algorithm.
The strong correlation, along with limitations of the earth surface in-situ observation hinders the improvement path of adding new band equations. We propose to improve T08 algorithm by introducing the condition number as well as a preprocessing module (remove open water from surface types) and removing an equation from the original linear system.
The results show that B245 has the smallest condition number and lead to more reasonable estimates of the MPF than the original T08 combination of B123. All the two-band results have lower MPF values in high MPF regions while have larger MPF in low MPF regions than those from three-band algorithms. Among them, B15 has the smallest condition number and also has the smallest bias from the MPF from Landsat discrimination.
The new MODIS algorithm outperforms the original version in reducing unrealistic solutions. It can also generate reasonable MPF distribution over broad area. Further comparison with MPF from University Hamburg shows that the B15 algorithm can generate reasonable outputs across the Arctic basin with a relatively larger value especially in Pacific sector.
The comparisons of results with Landsat-5 TM classification suggest that we still need to enhance the new algorithm's skill in treating some fine-resolution surface features; meanwhile, there may be some biases in the derived MPF values.
Investigations that combine in-situ and airborne observations will help to improve the algorithm. We are also working on using airborne images to validate the algorithm.