Fractal-Based Retrieval and Potential Driving Factors of Lake Ice Fractures of Chagan Lake, Northeast China Using Landsat Remote Sensing Images

: A thorough understanding of the freshwater ice process received considerable critical at-tention due to increasing winter recreations and ice engineering. The development of the lake ice process of Chagan Lake was monitored using MODIS and Landsat images over eight consecutive snow seasons from October 2013 to April 2021. We derived the lake ice phenology from an eight-day time series of lake water skin temperature (LWST) provided by MODIS, including freeze-up date, break-up date, and ice cover duration. We discovered a large-scale fracture extending from northwest to southeast that repeatedly appeared on Landsat images since 1986. A novel fractal-based auto-extraction is proposed to extract the length and angle of these fractures. We also carried out a field campaign and an ice ridge was found at the southernmost part of what we observed from the images. Moreover, we explained the fracturing development by thermal changes, wind in lake, and underlying flow. Results show that the lake ice fracture is nearly perpendicular to the dominant wind direction in the cold season, which indicates the crucial role of wind on lake ice fracturing.


Introduction
The freeze-thaw cycle of lake ice is an essential process of the cryosphere and hydrosphere, and serves as a sensitive indicator to climate changes [1,2]. About 65% of lakes throughout the world experience the seasonal freeze-thaw change during the winter, approximately 2.4 million km 2 [3]. The thaw-freeze cycle of lake ice is mainly controlled by thermal exchange, the displacement of water and ice, and water level variations, leading to lake ice fracturing. The existence of ice cracks could modify the optical transmission and then influence the algal communities and biological productivity beneath the surface ice [4]. During different stages of the ice process, the expansion and contracts of lake ice could exert a huge thrust on slope protection because of the bank, which could endanger the human infrastructure along the rivers and lakes. A thorough understanding of the lake ice process has significant implications for the geophysical, engineering, and biological utilization of freshwater ice during the cold season [5,6].
Previous work has mainly been attentive to the physical and thermal properties of ice cracks and has ignored the morphological features of lake ice cracking. The physical properties of snow and ice on freshwaters and the related mechanical strength tests have

Study Area
The Chagan lake, located in the southwest of Songnen Plain, Northeast China, is the largest natural lake in the Jilin Province and one of the ten largest lake areas in China, which have important implications for the fishery, agriculture uses, and winter recreations [31,32]. The geographic location of the study area is shown in Figure 1. It has a coverage area of 372 km 2 and a volume of 5.98 × 10 8 m 3 [33]. The water depth has an average value of 2.5 meters and a maximum of 4.5 meters [34]. The nineteenth edition of the Chagan Lake Ice and Snow Fishing and Hunting Cultural Tourism Festival was successfully held in 2020, which have created great social and economic benefits.
The Chagan Lake belongs to a mid-temperate, semi-arid, continental monsoon climate with well-marked seasons [32,33]. The region has an average temperature of 5.5 ℃, total precipitation of 450 mm, and an annual sunshine duration of 2880 hours [33]. The average runoff and evaporations are 451 and 1206 mm, respectively [24]. In winter, the Chagan Lake begins to freeze in late November and stops melting in the middle of April. The ice cover exists for as long as 150 days with a completely frozen duration of 130 days [25]. The ice thickness ranges from 0.8 meters to 1.1 meters [35]. The primary water source of Chagan Lake mainly consists of precipitation, groundwater, and irrigation discharge from three adjacent districts, including Da' an, Qian'an, and Qianguo [36]. The inflow and outflow of ground water in Chagan is 0.90 x 10 8 m 3 and 0.24 x 10 10 8 m 3 in summer, respectively, and the flow direction of groundwater is displayed in Figure 1. We carried out a winter investigation around Chagan lake for ice regime in December 2020, and a fall investigation for water depth on September 17, 2021. An artificial dike has been built on the lake's eastern shore, where a lakeside road with several viewing platforms and country lanes enabled us to investigate along the east coast of the lake. The fishing field in Figure 1 is the main place for winter recreation. The surface of lake ice along this side was flatter and smoother. The western shore has not been developed by human beings where wetlands and marshes are prevalent in summer. The surface of lake ice along this side was uneven; ice ridges and stacked ice are widely distributed. Figure 2 displays the recent spatial distribution of water depth measured using hand-held sonar detector SpeedTech SM-5 on September 17, 2021. The water depth measured ranged from 1.9 to 3.9 meters with an average value of 3.3 meters.

Landsat
Landsat 8 Satellite was launched on February 11, 2013, with Operational Land Image (OLI) and Thermal Infrared Sensor (TIRS). The satellite products have been widely used in freshwater monitoring and information extraction with eleven bands. The spatial resolution of bands 1-7 and 9-11 is 30 meters, and the spatial resolution of band 8 is 15 meters; the temporal resolution is 16 days. We downloaded the Landsat 8 OLI images from USGS. Global visualization viewer web (http://glovis.usgs.gov/), with the cloud cover less than 20% (Table 1). Forty-two images with fractures covered the cold seasons from 2013 to 2020. A cold season ranges from October to April of the following year, which characterizes the annual cycle of lake growth and decay. The row and column of satellite images are 119 and 29. The pre-processing steps include radiation correction, atmospheric correction, geometric correction, and mosaic. MODIS LST products averaged daily LWST during eight days and offered reliable LWST with a spatial resolution of 1 km. The data processing followed the strategies proposed by Song et al. (2016) and developed by Du et al. (2020), including MRT pre-processing, mosaicking based on quality control filter, band math, and time-series generation [38,39]. The mean average LWST is defined as the mean values of all the pixels in the Chagan lake for a given day. The time series of average daily LWST is calculating the average values of day-time and night-time LWST, and diurnal temperature differences that are calculated by day-time LWST minus night-time LWST.

Climate Records
Daily air temperature (AT), ground surface temperature (GST), wind direction, and wind speed of three stations around the Chagan Lake are provided by China Meteorological Data Service Centre for Climate Data (http://data.cma.cn/) for eight recent cold seasons since October 2013. Table 2 lists the details of three climate stations used in this paper, including Da'an, Qian'an, and Qianguo. MODIS-derived LWST were validated by the AT and GST values of three meteorological stations, marked as the red circle in Figure 1. The wind speed and wind direction are measured at an altitude of 10 meters above the surface. Table 2. The list of climate stations around Chagan Lake and the spatial distribution is shown in Figure 1.

Methods
With the development of remote sensing, the geospatial resolution has become more refined, and the texture information of remote sensing images has been greatly enriched. To start with, we extracted the lake ice phenology, which provides a temporal reference to check whether the textures on images are caused by snow and ice. Then, we proposed a fractal-based method to extract the length and angle for the large-scale fractures. Finally, we calculated the fractal dimension D to characterize the spatial distribution of lake ice fracture.

The Extraction of Lake Ice Phenology
Lake ice can be directly derived from a diverse set of physical properties, such as albedo, brightness temperature, backscatter coefficient, lake surface water temperature, and so on [25,28]. Herein, we derived the lake ice phenology from the time series of MODIS LWST. A variety of definitions of the lake ice phenology have been developed [28,40]. We define the freeze-up date as the day when the LWST fall below 0℃ in early winter, and the break-up date as the day when the LWST rise above 0°C in the spring of the following year; and ice cover duration describes the time duration between freeze-up date and break-up date [38,39].

The Fractal-Based Auto-Extraction
We implemented the following steps to find the lake ice fracture on satellite images, and calculated the longest length and angle, as shown in Figure 3. The blue band with a wavelength ranging from 0.450 to 0.515 μm is prepared and transformed into binary images for further analysis after pre-processing. Firstly, we extracted the linear structure using the Canny operator and generated the initial boundary. The linear structure is considered as the boundary of homogenous regions with different texture features. Secondly, we carried out the edge detection and excluded the outer boundaries. The removing part is determined by opening, filling, and eroding sequentially. Thirdly, the contour of the maximum connected domain of inter part was extracted. Finally, we located the south and north end of this fracture to calculate the angle and chose the shortest path for calculating the crack length. The angle referenced the 16-degree angle of wind direction and is defined as the angle between lake ice fracture and north direction in a clockwise direction.

The Calculation of Fractal Dimension
The fractal dimension of strictly self-similar objects can be derived mathematically and is given by Mandelbrot [13,14]: (1) where Nr represents an object of N parts scaled down by a ratio of r. The D derived from Equation (1) is the shape's similarity dimension [13][14][15]. We used the differential box counting method to calculate D [15] . The image with a size of M×M has been scaled down to a new image of s×s，then: the maximum value and minimum value calculate the gray level of each grid: Nr is the sum of contributions of each grid, Substituting Equation (2-4) into Equation (1) gives,

The Changes of MODIS-Derived LWST
We compared the MODIS-derived LWST with GST and AT from the three climate stations, and the scatter plot is shown in Figure 4. MODIS products provide eight-day records by averaging LWST within specific eight days. If the date of field falls within these eight days, we matched the LWST with in-situ climate records. Three hundred and fortysix in-situ records were selected for comparing the difference of MODIS LWST, GST, and AT from October 2013 to April 2021. The coefficient of determination (R 2 ) between GST and MODIS LWST had high values of 0.92 with the of RMSE (4.049). Meanwhile, a higher R 2 of 0.94 was found between MODIS LWST and AT with the RMSE of 3.432. These results reveal that the MODIS-derived LWST has good agreement with both GST and AT. Therefore, the MODIS LWST is reliable for describing the spatial-temporal changes of lake water skin temperature. We present the annual and seasonal changes of AT, GST, and mean LWST in Figure  5. AT showed a slightly decreasing trend in view of yearly changes, and the GST and LWST displayed a flat trend. Compared with GST, the LWST was closer to AT. What is interesting is that AT, GST, and LWST are almost identical in the cold season 2019. The seasonal changes in Figure 4(b) suggested that the difference between AT and GST is minor, common in the relatively shallow natural lakes. The air is warmer than water in summer, and the water is warmer than the air in winter, suggesting the temperature index is a good indicator of thermal changes. Generally, the LWST and AT showed the same trend except for the cold season of 2015. The errors source of MODIS-derived LWST could be explained by the time difference of satellite overpass and field measurement, the scale difference between pointed-based filed measurement and pixel-based satellite images, emissivity estimation errors, the spatial gap caused by cloud, and so on [33,41,42].   Figure 6 displays the comparison of the auto-extraction and visual interpretation. The angle based on auto extraction had mean values of 335°28′34″ ± 0°13′34″, and the length had mean values of 21141.57 ± 68.36 km. We compared the angle and length of lake ice cracks using auto-image extraction and visual interpterion, as listed in Table 4. The R 2 of angle and distance between auto-extraction and visual interpterion was 0.96 and 0.98; the MAE was 3.99% and 2.48%. These results show that the auto-extraction method agrees with visual interpterion well, and the auto-extraction results are reliable.

The Results Based on Landsat Images
The fractal analysis has a wide range of applications in diverse fields, such as geoscience, materials, meteorology, and so on [15]. We calculated D values of 42 Landsat images for these eight cold seasons. Figure 7 shows the boxplot of D values and the Landsat images with the highest D values within the given cold season. Generally, the fractal dimension D at the lake-based scales varied slightly within the range of 1.8569 and 1.8756 during these eight cold seasons. The cold season from 2014 to 2018 showed a relatively wide distribution, ranging from 1.8596 to 1.8759, with a mean value of 1.8663. The cold season from 2018 to 2019 exhibited a concentrated distribution, ranging from 1.8586 to 1.8605 with mean values of 1.8595. The Landsat image of this cold season also has good image quality due to scarce snow and clouds and offers excellent material for fracture analysis. The fractal demission D expressed an overall statistical indicator for the whole image. The monofractal dimension index is not sufficient to explain the complex nature phenomena beneath the lake ice surface. We checked the quick images of three sensors of Landsat since 1986, including TM, ETM+, and OLI. A similar linear structure had been found on the Landsat quick images during 18 of those 35 cold seasons. Figure 8 illustrates the images with the best quality in each cold season and corresponding dates. One scene with the best image quality is presented herein for a given cold season. Liu et al. (2018) also identified similar lake fractures and analyzed the impact factors of the mechanical processes on three lakes in northeast China. These findings pointed out that image crack is a common phenomenon that has not attracted the interest of humans. The infiltrated lake snow ice played an unavoidable role in the development of lake ice fracture. Substantial snowfall is prevalent in temperate areas like northeast China. When the weight of snowfall exceeds the capacity of surface lake ice, it generates the lake ice fracture [43]. To sum up, the lake ice cracks may be caused by snow cover, self-similarity, and image quality from the view of satellite images. Previous work showed that the large-scale fractures usually appear in the blue ice regions, where snow cover is blown away by the wind and showed a rarely distributed [5]. No matter if snow-covered or not, the fractal analysis provided an effective solution to analyze the texture features of lake ice surface on satellite images.
Self-similarity and fractal dimension, a key concept of fractal theories, means that the fractal dimension is supposed to be stable when applied to remote sensing images. The D dimension showed a minor change in Figure 6, which suggests that the remote sensing data used are not strictly fractal. Previous work concluded that natural phenomena have a statistical self-similarity over a limited range of scales [15]. Besides, we also used fractalbased edge detection to calculate the length and angle of image fractures. Self-similarity is not a prerequisite when applying fractal theories to extract textural information [44,45].

The field investigation
We carried out a field investigation to justify our findings at the end of December 2020. We located two spots around the southernmost edge of the lake ice fracture to identify the lake ice fracture on December 20, 2020. There existed an ice ridge around Spot A (124°18′35.25″ E, 45°9′56.71″ N), as shown in Figure 9. As seen from Figure 9(b) and Figure  7(c), unshaped oval ice cubes were widely distributed in this area. There existed a transparent frozen water path with a lighter tone and smoother surface than the neighborhood around Spot B (124°17′59.2476″E, 45°9′59.68″N), as shown in Figure 9(c). The action of squeezing pressure made the crushed ice clutter together and formed a hill-like pile of ice. More field photos can be found in Figure A1. The rugged surface suggested that the lake ice deforms along the shore due to thermal and mechanical stress. Therefore, we concluded that the image fracture generated from ice ridges or stacked ice is one of the reasonable explanations for what we observed from Landsat.

The wind field of Chagan Lake
We examined the recent eight cold seasons' wind speed and direction data collected by three meteorological stations around the lake area. Figure   The wind rose in Figure 11 displays spatial distribution and frequency of wind direction and velocity. Chagan Lake is normally characterized by strong easterly and southerly in the warm season. The dominant wind direction differed among different months. The wind domain from August to October was SSW, and that from November to March was SW, and that of April was SSW. Moreover, SW, SSW, and WSW had the three highest frequencies in the winter. The lake fracture had an average angle of 335°28′34″, the WSW and SW direction had the angle of 247°30′ and 225°. The angle between fracture and WSW is 87°58′34″, nearly perpendicular. We treated the fractures as straight lines instead of curves when calculating the angle. A closer look at the cracks on the ice surface reveals that the cracks were not straight lines, but along several major cracks that extended, diverged, crossed with other cracks, and then extended and bifurcated. The direction of lake fracture is nearly perpendicular to the wind domain direction in the cold season, which is responsible for the direction of lake ice fraction. The prevailing direction also forced the lake ice to drift from the southwest to the middle, revealing that the southwest side is unsuitable for human recreation or fishing. Interestingly, the direction of lake ice fractures is also perpendicular to the flow direction of groundwater. Both inflow and outflow get slow in winter and make the water level come down, leading to perpendicular cracks.

Discussion
A comprehensive understanding of the lake ice process is increasingly significant due to prevailing human recreation and fishing. The evolution of the lake ice is controlled by a balanced exchange of mass, heat, and momentum between lake and atmosphere, which largely depends on the local climatology [5]. Although the earliest lake ice phenology could trace back to the 1800s, the long-term record is available for limited lake areas in cold regions [46]. Temperature has been considered the most significant impact factor of the frozen-thaw cycle of lake surface [22,25]. Ice formation is governed by intense radiative and turbulent heat losses between the warmer lake surface and the colder atmosphere. Herein, we considered three types of temperature (AT, GST, and LWST). All of them underwent minor changes during the past eight cold seasons, showed the same trend with lake ice phenology. In other words, the heat changes rather than the air temperature index control the lake ice process.
During the freezing process of Chagan Lake, the water body is subjected to different degrees of expansion forces, resulting in dry cracks, wet cracks, and wide cracks in sequence [5]. The surface ice contracts as the surface temperature decrease, resulting in dry cracks. The water comes up, fills the air gap, and freezes under the wind and waves' synoptic effect. The dry crack just happens on the upper layer of ice without penetrating. Then, narrow wet fractures accumulate rapidly and supply continuous materials for a steady ice sheet. After that, wind waves with long and narrow ice lead appear and last for several days until it complete freeze. Meanwhile, a thin ice bridge usually formed, and the ice thickness is supposed to be less than 20 cm [5], which was consistent with our results from satellite images and further confirmed by the fieldwork.
As the temperature increases in the break-up process, the surface ice expands, generating small displacement and occasional ice drift. A steady ice sheet could be broken by the wind. Taking Chagan Lake as an example, the wind blows in the southeastern direction synchronously and creates intense pressure to the windward side of the lake. Thus, the Chagan Lake begins to melt along the southeast side and expand across the lake area. Herein, we used the climate records to describe wind changes, and the difference of wind field between lake areas and surrounding land is ignored. The pointed-based climate record had limited abilities to characterize lake-based or basin-based climate changes. Wind, air temperature, and humidity together influence the turbulent transfer over the air-lake interaction. The semi-elliptical ice egg found in field investigation also suggested that the lake ice deforms along with the lake-land interaction under wind and waves. Notably, the formation of recurrent ice fracture closely associates with the shoreline geometry and topography [5], which requires a detailed bathymetric measurement in future work.
The role of the flow field under the ice cover has been ignored in previous work, especially before freeze-up. Liu et al. (2020) constructed a hydrodynamic model for Lake Chagan based on the Delft3D platform [32]. The hydrodynamic field of Lake Chagan in October was driven by the wind direction and the digital elevation model (DEM) of the lake bottom and calibrated by the water depth and water temperature from field measurement. From the planar distribution map of flow velocity ( Figure A2), the lake water showed a clockwise flow along the western shore and a counterclockwise flow from along the eastern shore, and the junction line showed a northwest-southeast direction. Under the influence of lake hydrodynamics, the dynamics interfacial reaction of lake ice is intense when the lake ice freezes, which provides a beneficial condition for lake ice fractures. We concluded that lake hydrodynamics before the freeze-up also provide reasonable explanation the lake ice fracture of Chagan Lake shows a northwest-southeast trend caused by lake hydrodynamics before the freeze-up process.

Conclusions
Remote sensing provided a promising tool for exploring the large-scale lake ice fracturing over different lakes, which may have promising applications in ice engineering. Chagan lake froze on November 13 in winter and melted on April 1 in spring, leading to the ice cover duration of 140.44 ± 12 days during these eight years. The large-scale lake ice fracture appeared repeatedly in the Chagan Lake, 18 times of 35 cold seasons since 1986. We proposed a fractal-based method to calculate the angle and length, and the large-scale fractures had an average angle of 335°28′34″ ± 0°13′34″ and the length of 21141.57 ± 68.36 km. Fractal geometries supply a helpful tool for exploring the spatial distribution in remotely sensed images, especially for the linear structure.
Moreover, we carried out a field investigation to validate the existence of the lake ice feature in December 2020. The lake ice surface along the southwestern side of the Chagan lake differed greatly from the other side due to complex underwater topography and lakeside landscape. An ice ridge has been found on the southernmost edge of fractures on satellite images and explained what we have observed. We also analyzed the influence of temperature, snow, and wind on fracture generation, and found out the lake ice fracture is nearly perpendicular to the domain wind direction in winter. The flow field of October showed a northwest-southeast direction, which implies that lake hydrodynamics before the freeze-up process also played a crucial role in fracturing.
The work herein reveals that remote sensing provides a macro perspective to monitoring the large-scale lake ice fracturing using fractal methods. We plan to carry out a bathymetric measurement along the southeast shore to find the relationship between surface morphology and underwater topography. Further exploration explores more texture measures and the related computing methods for multifractal modelling and consider the scale effects of multi-source remote sensing data with the multi-band information and spectral features on physically-based ice processes.

Acknowledgments:
The authors would like to thank the Earthdata Search for providing MODIS, and Landsat data. Thanks are given to Professor Zhijun Li for the guidance for field validation. The anonymous reviewers to improve the quality of this manuscript are greatly appreciated.

Conflicts of Interest:
The authors declare no conflict of interest.