Multifractal analysis of glaciers in the Lombardy region of the Italian Alps

Glaciers are continuously monitored to detect their spatial extension and time evolution since they are the best witnesses of climate changes. There is a particular interest for Italian glaciers in the Alps as there is evidence that they are melting at a faster rate than those located in other regions of the globe. The determination of the perimeters of glaciers represents an effective method to evaluate the area covered by them. The availability of data for the perimeters encompassing several years suggests the opportunity of correlating the morphological variations in time with the properties of their shrinkage. In this work, we investigate the multifractal properties of the perimeters of the Lombardy glaciers in the Italian Alps. We characterize the area and perimeter distributions of the population of glaciers and we show that the distribution of perimeters exhibits a marked peak, not present in the distribution of areas. We analyze the area-perimeter relation, which is characterized by a power-law behavior that indicates a fractal structure of the perimeters with fractal dimension 1.2, independently from the size of the glaciers. We investigate the multifractal spectra of perimeters and we show that their features are strongly correlated with the area of the glaciers. Finally, we study the time evolution of the area and perimeter of glaciers and we detect the existence of a large class of glaciers whose perimeters increase while their areas decrease. We show that this behavior has a well definite counterpart in their multifractal spectra.


Introduction
The study of glaciers and the monitoring of their evolution in time is a problem of great interest, as it is directly related to the effects of global warming [1][2][3][4], although the melting rates are not trivially explained by climate variability. The cryosphere holds about 69% of the world's freshwater, providing an invaluable resource for civil and industrial use, and feeding large river systems in high-mountain and densely populated regions of the world where other sources of water are scarce [5]. For example, Central Chile is semi-arid with mostly dry summer, and rain in winter, and snow and glaciers' melt contributes runoff in the summer up to 70-80% [6]. The economy of the Hindu-Kush-Karakoram-Himalaya region relies upon agriculture, and likely more than 50% of the water in the Indus river originating from the Karakoram comes from snow and glacier melt [7]. For example, for a period of 18 days in summer 2011, Senese et al [8] quantified a total water volume of 1.54 km 3 derived from ice melting, corresponding to 11% of the reservoir capacity of the Tarbela Dam, a very large dam on the Indus River that plays a key role for irrigation, flood control and the generation of hydroelectric power for Pakistan. On Italian Alps, and especially in the province of Sondrio in the Lombardy region where twentyseven hydropower dams are located, about one-half of the total water for hydropower plants is provided by snowfall and ice melt [9]. Observations of glacier coverage and evolution are then essential to understand the role of the cryosphere in influencing the regional hydrology and water resources in these regions.
Recent studies have demonstrated that the shrinkage of glaciers is particularly severe upon the Alps [1][2][3][4]10]. Most alpine glaciers have experienced a more or less continuous mass loss since the end of the Little Ice Age (ca 1850) [4,11] and over recent years this phenomenon has been intensifying [12]. In addition, Italian glaciers lost 30% of their area, from 526.88 km 2 in the 1960s to 369.90 km 2 in the past decade [10,13]. The spatial extension of glaciers is monitored regularly by determining the perimeters of glaciers from highresolution orthophotos and storing the results in inventories. So far the study of the morphology of glaciers has mainly involved their orientation and shape anisotropy. It could be particularly interesting to characterize the morphologies of glaciers and try to find some relations with the evolution of their shrinkage.
During the last 40 years, several methods have been developed for the statistical characterization of the geometrical properties of natural systems. In particular, the evidence of fractality and multifractality in many natural systems, has stimulated the effort to get information about physical properties through the use of mathematical modelization [14].
The concepts of fractal geometry have been applied to investigate the geometry of a wide variety of systems, such as coastal profiles [15], medical data [16], colloidal aggregates, and critical phenomena [17][18][19][20], only to cite some examples. The limitations of the fractal analysis became apparent when it was realized that systems as the ones mentioned above were not completely described by a unique fractal dimension, but instead required a continuous spectrum of exponents. This led to the development of multifractal analysis [21][22][23][24][25][26], which has been widely applied in many branches of physics, such as astrophysics, medical physics, fluid dynamics [27,28], and natural and geological systems [29,30]. However, notwithstanding the ample utilization of multifractal analysis in several scientific fields, so far its application to the investigation of glacial landscapes has been limited to the investigation of the elevation profiles of glaciers, in comparison to those of fluvial landscapes [31].
Besides the application to natural phenomena, research in this field has also advanced for the mathematical aspects and has led to the availability of more and more refined tools for studying fractal and multifractal sets [14,[32][33][34].
In this work, we present an application of multifractal methods to the study of the morphological properties of Italian glaciers in the Lombardy region in the period from 2003 to 2012. We perform a statistical characterization of the population of glaciers, determine the multifractal spectrum of glaciers, and correlate its features with the area and its time evolution. We will present a concise and partial description of the methods of fractal and multifractal analysis for which we follow the approach used in similar applicative works [15,27,31] and adopt the mathematical formalism used in reference [35]. Then, we describe the procedure used to collect and analyze the data and we finally discuss the results obtained.

Fractal analysis
Fractal structures are characterized by diffuse irregularities that cannot be properly described by Euclidean geometry and require to introduce non-integer fractal dimension, different from the Euclidean one D. One representative example of a natural fractal is that of coastal profiles introduced by Mandelbrot [36], but the same line of reasoning applies more generally to an open set of area A bounded by a connected perimeter of length L. In order to measure the length of a jagged coast, one can utilize the box-counting method, where the coastline is first covered by square boxes of size δ and the coastal length is estimated from the relation L TOT = N(δ)δ, where N(δ) is the number of boxes of size δ needed to cover the coastline. If the coastline was a one-dimensional line, its length would be independent of the size of the boxes used to cover it, because N(δ) would increase inversely with their size. At variance, the coastal length of natural islands depends on the size as L(δ) ∝ δ D−D f , where D f is the box-counting (Minkowski-Bouligand) fractal dimension. In the limit of small δ, the number of boxes needed to cover the coastline scales as: This power-law scaling reflects the self-similarity of the coastline, which does not exhibit a characteristic length scale. The presence of a fractal coastline also affects the size of the area A bounded by it. In the case of convex polygons A ∝ L D . Conversely, it can be shown that when an area A of a surface in a bi-dimensional space is bounded by a connected fractal perimeter of length L, the two are related by a power-law relation [35]: where 1 < D f < 2.

Multifractal analysis
The fractal analysis portrays some of the geometrical features of systems that exhibit scale invariance. However, most of the natural systems, if observed in deep detail, can be better described only through a more general approach. In particular, many systems present irregularities whose statistical properties cannot be described by means of a single fractal dimension. For these systems, the scaling exponents are sensitive on both the dilation factor and the position. The multifractal formalism is effective in describing such structures, as it is able to account for the local properties and not only for the global ones [27,28,32,35]. Let us consider again the problem of the measurement of the coastal length. By adopting a box-counting method like the one described in the previous section, all the square boxes that intercept the coastal profile are counted with the same weight, without making any difference between squares containing long, short, or multiple portions of the curve. In other words, the box-counting approach used to unveil the fractal structure only considers the possibility of a true/false occupancy of a box. With a more accurate description, a measure could be assigned to each square, able to account for other properties than the simple presence of a curve's section inside it.
To clarify this point, we will rely on a practical example. Let us consider a digital representation of the profile of Italy ( figure 1(a)). The profile is represented by a finite set of black pixels arranged on a square lattice, in the customary bitmap representation of digital images. By dividing the lattice into square boxes of size δ, we can count the number of boxes N(δ) that contain at least one black pixel and determine the fractal dimension of the profile from equation (1). By adopting this simple box-counting method, any information about the level of occupancy of different squares would be lost. A better way to characterize this signal is to count the number of black pixels n i (δ) in each square of size δ and divide it by the overall number n(δ) of pixels available in the box. This allows to define the integrated measure of the probability of occupancy of each box μ i (δ) = n i (δ)/n(δ). In the limit of small δ, this probability diverges as where α i is the Lipschitz-Hölder exponent that characterizes the local divergence of the probability of occupancy of the box i. The divergence is apparent in the sequences (b)-(e) and (f)-(i) provided in figure 1, where the measure is evaluated by progressively reducing the size δ of the squares used to cover the profile and is represented both with false colors ranging from black to yellow and with a histogram in sequence (f)-(i). From the figure, one can appreciate that by decreasing the size of the boxes used to cover the profile, localized spikes appear on it. In the limit of small δ, the number N(α) of boxes with Lipschitz-Hölder exponent comprised between α and α + dα diverges as [24,26]: where f(α) is usually named multifractal singularity spectrum [14,32,35]. From the comparison of equations (1) and (4) we see that f(α) may be interpreted as the fractal dimension of the subset of boxes with the same Lipschitz-Hölder exponent α. The multifractal spectrum f(α) is a continuous function defined for α > 0, characterized by a negative second derivative, and maximum f max (α) that coincides with the fractal dimension of the set that supports the values of the measure.
Another way to account for the different measures μ i of the boxes used to cover a multifractal object is to start from the moments of order q of the measures [14,24,26,35]. This formalism allows to define two other multifractal spectra, τ (q)-the sequence of mass exponents-and D(q)-the generalized dimension function-both containing the same information of f(α), but somehow easier to be interpreted. By defining the weighted number of boxes N(q, δ) = i μ q i as the sum of the moments of order q of the probability distribution μ i , in the limit of small δ, it diverges as: We can appreciate that, for the moment's order q = 0, we recover the box-counting divergence of equation (1), as μ q i = 1 and N(q = 0, δ) = N(δ) is the number of boxes needed to cover the multifractal object. As a consequence, τ (q = 0) corresponds to the fractal dimension D f of the set. Remarkably, f(α) and τ (q) are simply related through a Legendre transformation: Therefore, they both contain exactly the same information about the multifractal set. Finally, starting from τ (q), the generalized dimension D(q) can be defined as [22,26]: D(q) is a generalization of the fractal dimension to the case of a multifractal set and consequently, it is often used in substitution to the previous functions. For a non-fractal object, D(q) is constant and simply coincides with the Euclidean topological dimension, while for a monofractal object D(q) = D f at all the values of q. At variance, for a multifractal object, different q values provide information about distinct regions of the system, and D(q = 0) = D f is the fractal dimension of the system. For high q values (q 1) the contributions from boxes with high values of μ i are more favored. On the other side, for q −1 contributions from boxes with a lower measure μ i are more relevant. In fact, for q > 1 one has that μ q i > μ q j when μ i > μ j , while when q < −1 the opposite occurs. In this paper, we will use D(q) for the characterization of the multifractal properties of glaciers because it allows a more straightforward interpretation of results.

Methods
A glacier is a dynamic and fragile ice body, always moving forward under its own weight, and it forms through the natural metamorphosis of snow.
In this paper, we have applied the multifractal analysis to the glaciers of Lombardy, a northern Italian region, monitored in the years 1983, 2003, 2007, and 2012. In figure 2, as an example, the glacier named Pisgana Ovest is sketched from the perimetral coordinates contained in the four inventories. The progressive shrinkage in years can be well appreciated. In order to investigate the dependence of the fractal and multifractal properties from the size of glaciers, we took advantage of the availability of inventories where the area of glaciers spans about five orders of magnitude. This allowed us to analyze an extensive and statistically significant population of glaciers. Incidentally, it has to be remarked that an ice mass can be considered to be a glacier when its minimum mapped area is 0.01 km 2 , as recommended by Paul et al [37]. Even if bodies smaller than 0.1 km 2 are classified as a glacier, they cannot feature the dynamics and behavior typical of bigger glaciers.
For mapping 2003, 2007, and 2012 glacier boundaries we analyzed color orthophotos, which are purchasable products distributed by Compagnia Generale Riprese Aeree Spa (BLOM-CGR SPA, flight code IT2000). These products are derived by processing high-resolution color aerial photos acquired in the summertime during days of clear sky conditions or featuring a very low cloud coverage (i.e. minor than 30%). Moreover, the photos were acquired when the snow coverage in mountain areas is minimum thus reducing the uncertainty in detecting and mapping glacier limits (see also [13,[38][39][40][41]). The digital orthophotos, which feature a planimetric resolution of 1 pixel (pixel size: 0.5 m × 0.5 m), were managed through a GIS (geographic information system) where they represented the base layer to detect and map glacier limits. To assess the reliability and the accuracy of the obtained polygons representing glacier boundaries, we performed several crosschecks by comparing our data with other ones derived from different sources (e.g. local glacier database, technical maps, and regional cartography, data from dedicated field investigations, etc). As regards the 1983 glacier boundaries, they were obtained by analyzing couples of aerial photos with an optical stereoscopic system. Through this instrument, the scientist can obtain a 3D view of the study area thus being able to detect and map the glacier limits. The glacier boundaries were reported in the GIS project as polygons using as raster base the technical regional map (named CTR) of the Lombardy region, at the scale 1 : 10, 000, derived from an aerial survey performed in 1983. The availability of a map surveyed in the same period of the aerial photos used to detect and map glaciers permitted us to assess the accuracy and reliability of our work. The final planimetric accuracy of the 1983 glacier database was evaluated to be equal to ±5 m. Due to both the lower accuracy of 1983 data (with respect to 2003, 2007, and 2012 data) and the different method applied to map glacier boundaries, this database was not used for the multifractal analysis. Beyond perimeter coordinates, inventories contain much information about each glacier such as their name, identification number, area, perimeter exposition, presence of supraglacial debris. Should a glacier have fragmented due to shrinkage, the various bodies are, since then, numbered as the original glacier identified in 1983, followed by a numerical extension. This convention allows keeping track of the history of each glacier in time.
In order to calculate the multifractal properties of the glaciers, we have customized FracLac, a plugin of the OpenSource program ImageJ, to enable it to read the inventories, extract the perimeter coordinates, make multifractal analysis and produce the three functions defined in equations (4), (5) and (7).

Area and perimeter distributions
The Lombardy region during years 2003-2012 comprised about 300 glaciers with a wide distribution of areas, the surface of the smallest glacier being of the order of 500 m 2 (Pizzo del Ferro Inferiore II-2007) and the largest 2 × 10 7 m 2 (Adamello-2003). Statistical characterization of the distributions of area and perimeter of the glaciers can be achieved by plotting their probability density functions (pdf). We consider a property X of a population, where X in the case of glaciers can be either the area or the perimeter of a glacier. The probability that X is comprised between x and x + dx is given by p(x)dx, where p(x) is the pdf. The pdf is defined by the relation p(x) = dN N dx , where N is the total number of elements of the population and dN is the number of elements with property X comprised between x and x + dx. The log-log plots of the pdfs for the population of glaciers ( figure 3) show that the distributions of area and perimeters are dominated by a large number of small glaciers. As it commonly occurs in the case of natural phenomena encompassing a wide range of length scales, the pdfs for the area and perimeter follow asymptotically a power-law behavior p(x) ∝ x −γ at large x, indicating the absence of a typical value of x in the range where the power-law behavior holds. An accurate quantitative characterization of the pdfs is prevented by the fact that their determination requires to group the population into bins corresponding to ranges of X, and the fitted parameters for the pdfs are significantly affected by the chosen binning representation. As discussed in detail by Newman [42], an accurate characterization of the probability distributions can be achieved by using a cumulative probability function (cpf) instead of a pdf. In the case of the area of glaciers, the cpf P(A) can be defined as the probability that the area of a glacier of the population is larger than a chosen area A. Therefore, p(A) and P(A) are directly connected to each other by the relations: and The advantage of representing the statistical distribution of a population by using a cpf rather than a pdf is that the first is defined for every element of the population and, as a consequence, a binning procedure is not needed and the number of degrees of freedom is not reduced. Log-log plots of the cpfs of the area and perimeter of glaciers are shown in figure 4.
From their asymptotic behavior at large x one can appreciate that they also obey a power law relation p(x) ∝ x −(γ−1) . The cpfs can be characterized accurately by fitting the data in figure 3 with the empirical function: where x 0 and x m are the smallest value and the median value of the population, respectively. By fitting the data with a statistical weighting where all the points have the same relative indetermination, one can appreciate that the model function defined in equation (10) is in excellent agreement with the experimental data, both for the cpfs of area and perimeter. From the values of the parameters obtained by fitting the cpfs, detailed in table 1, one can appreciate that the median of the distribution of area decreases in time, providing a statistical signature of the shrinkage of the area of the population, while that for perimeters is fairly constant. All the cpfs for the area are compatible with a unique value of the power-law exponent γ a ≈ 1.9, independently from the year, and the same occurs for those of the perimeter, but with a different power-law exponent γ p ≈ 2.4. The high exponents indicate that the populations are extremely unbalanced towards the smaller sizes. Noticeably, there is no evidence of substantial changes in the shapes of the curves in years.  By combining equations (9) and (10) we obtain an analytic expression suitable to describe the pdf: which, in combination with the value of the parameters detailed in table 1, provides an accurate empirical model to describe quantitatively the pdfs (solid line in figure 3). Although the cdfs of area and perimeters are both defined by equation (10), from figure 3 one can appreciate that the shapes of the pdfs of area and perimeter differ significantly from each other. More in detail, the pdfs of area decrease monotonically as the area increases, while the pdfs of perimeter exhibit a peak at values of the order of 400-500 m. This marked difference between the pdfs is determined by the fact that the cdf defined by equation (10) and used to fit the experimental results exhibits a transition for γ = 2. In fact, for γ > 2 the curve exhibits an inflection point, which determines a peak in the associated pdf, while for γ < 2 the inflection point is not present. The two dependencies lead to very different statistical behaviors of the area and perimeter distributions, with the meaningful consequence that the absolute maximum peak of the pdf of the area corresponds to the smallest glaciers, while the one for the pdf of perimeter occurs for glaciers whose perimeter is a factor 4-5 larger than that of the smallest glaciers.

Fractal analysis
A first understanding of the fractal properties of the population of glaciers can be achieved by investigating how the areas of glaciers scale with their perimeters. According to equation (2), in the presence of a fractal perimeter the area should scale as a power-law of the perimeter, with a non-integer exponent directly related to the fractal dimension of the population. We have plotted the glaciers' area as a function of perimeter (figure 5) and noticeably the data of the different years are fully compatible with a power-law behavior. The fitting of the data shown in figure 5 with equation (2) allows to establish that the average fractal dimension of the populations  of glaciers is D f = 1.21 ± 0.01, a value compatible with the fractal dimensions of all the years considered. To achieve accurate fitting, data have been weighted by assuming that they all had the same relative statistical error. This procedure allows to avoid that fitting parameters are biased by the glaciers with the largest areas. The analysis of the perimeter-area relation provides a first confirmation that the perimeters of glaciers have fractal geometry. The scale invariance associated to the fractal geometry of the perimeters would imply that a small glacier is entirely equivalent to a larger one but, as we will see, this is not the case.

Multifractal analysis
A deeper understanding of the morphological properties of glaciers can be attained through the use of a multifractal approach. In figure 6 we show the three multifractal spectra obtained for a typical glacier of small area (i.e. 850 m 2 ). Much equivalent information on a glacier can be extracted from the three multifractal spectra. The fractal dimension can be alternatively determined from the values τ (q = 0), D(q = 0) or f max (α). As already pointed out, a simple fractal would be described by a function D(q) = D F for any q. On the contrary, the large difference between the asymptotic values D −∞ and D ∞ of D(q) for q −1 and q 1, respectively, is an indication of a multifractal behavior, as it means that the system presents an ample range of fractal behaviors. The same information can be extracted from either of the two alternative spectra, namely the function f(α) and the sequence of mass exponents τ (q) defined in equation (5) [24,26,35]. Indeed the two slopes of τ (q) for q 1 and q 1 coincide with D −∞ and D ∞ , respectively. Analogously, the two values α min and α max where f(α) = 0 are respectively equal to D ∞ and D −∞ .
The main purpose of our study is to find whether there is any correlation between multifractal properties and the glaciers' size. To investigate this point, all the glaciers belonging to each inventory have been analyzed. The generalized dimension functions are plotted superimposed in figure 7 for each of the three inventories. Different colors correspond to different glaciers' sizes, as detailed in the color-bar. Remarkably, the curves show similar trends depending on glaciers' size. For large glaciers, D −∞ and D ∞ are much closer than in the case of small glaciers, which exhibit a more varied behavior. The little difference between the two asymptotic values of the curves referring to glaciers larger than approximately 10 5 m 2 , suggests that large glaciers are those that could be better approximated by simple fractals, described by flat curves. At variance, the values of D −∞ increase when glaciers become smaller, while the values of D ∞ span a smaller range and do not show an evident trend. These behaviors can be due to different dynamic processes typical of glaciers smaller than 10 5 m 2 compared to larger glaciers. Indeed, the smaller the glacier, the faster it will respond to changes in climate. Moreover, the smaller the glacier, the thinner the glacier thickness will be, and the ice thickness is known to have a great influence on glacial movement and ablation. Therefore, depending on thickness, area, and external climate, the dynamic process involved in how glaciers change and whether glaciers' changes in melting tend to be stable or irregular will be different.
The theoretical investigation of closed loops generated by synthetic rough surfaces has shown that the fractal properties of the loops are intimately related to the roughness of the surface [43]. In particular, the investigation of the area of the regions bounded by the closed loops generated by self-affine surfaces showed that their cumulative distribution obeys a power-law, whose exponent is linearly related to the scaling exponent for the roughness of the surfaces [44]. Similarly, the investigation of synthetic multifractal rough surfaces has shown that the loops exhibit multifractal properties determined by those of the surfaces [45]. The perimeters of glaciers do not represent contour loops obtained by cutting a rough surface, because the elevation of the points changes along the perimeter. However, the theoretical study mentioned above suggests that the multifractal properties of the perimeters of glaciers should be in principle closely related to the multifractal properties of the surface of the Earth in the regions populated by glaciers, and paves the way for future in-depth analysis of the relation between the multifractal properties of the perimeters and those of the orography, which is beyond the aim of this work.
The analysis of the features of the generalized dimension of the perimeters of glaciers allows to provide further insights about the origin of their multifractal structure. In figure 8  As anticipated at the beginning of the paper, almost all glaciers are shrinking and their area is decreasing in years. Surprisingly, however, the same does not always occur for perimeters. For a not negligible amount of glaciers, the perimeter increases while the area shrinks. These glaciers are almost 20% of the whole in the time range 1983-2012 but, when the year interval is restricted to 2007-2012, they become 40% of the total. Due to the relation between area and perimeter described in equation (2), we wondered whether this particular behavior somehow affects the generalized fractal dimension. We have therefore compared the variation of D (−20) in time with that of the fractal dimension D F = D(0). In figure 9(a)    of the corresponding color marks the centroid of the distribution of points, as can be better appreciated in the expanded plots shown in figures 9(c) and (d)). The two groups of glaciers for which areas increase in time comprise only a low number of small glaciers. More interesting is to compare the different behavior of the two more densely populated categories: that of the glaciers for which both area and perimeter decrease in time, and that of the glaciers whose perimeters grow while the areas shrink. For both classes, we see that the centroid lies in the semi-plane of positive ΔD (−20). This means that, on average, ΔD(−20) increases while the area diminishes, as already evidenced in figure 7.
A different behavior is observed for the fractal dimension. When both area and perimeter decrease in time, D(0) does not change much (green crosses in figure 9), while when perimeters increase in spite of area shrinkage, the fractal dimension D(0) increases in time (red crosses in figure 9). A possible interpretation is that a growth of perimeters accompanied by a shrinkage of areas gives rise to a significant increase of curves' curdling, reflected by the behavior of D(0). Remarkably, for these glaciers, the generalized fractal dimension simply shifts towards larger values without appreciably changing its shape. This is an indication that for glaciers belonging to this category the degree of multifractality does not change much in time. At variance, when both area and perimeter decrease, multifractality increases, as the value D(−20) rises while D (20) and the fractal dimension D(0) do not change significantly in time. Furthermore, recalling that for q −1, D(q) is affected by contributions from boxes with a relatively low value of the measure μ i , we can state that, when both area and perimeter decrease, there is a significant increase of boxes characterized by a low degree of occupancy without any significant variation of other properties.

Conclusions
We have characterized the perimeters of Lombardy glaciers of the years 2003, 2007, and 2012 using fractal and multifractal formalisms.
We have studied the area and perimeter distributions of the population of glaciers and we have shown that the perimeters' distributions exhibit a marked peak, not present in the distributions of areas. We have shown that, independently from their size, glaciers are well described as fractals with an average fractal dimension 1.2, as evidenced by the power-law area-perimeter relation. The investigation of the multifractal spectra of perimeters evidenced the existence of a correlation between multifractal properties and glaciers' size. Small glaciers, typically smaller than 10 5 m 2 , exhibit a strong multifractal morphology, characterized by a wide range of generalized dimensions D(q) that, on the contrary, is narrower in the case of larger glaciers that are more similar to monofractals. The multifractal spectra of small and large glaciers are significantly different at negative values of the moment q, and this leads us to conclude that the perimeters of small glaciers are more affected by curdling than those of large ones. This can be due to different response time, thickness, and dynamic processes, featured by smaller glaciers compared to the larger ones. We finally identified four glaciers' classes, determined by the four combinations of growth and shrinkage of area and perimeter that can occur within an observation period. Glaciers whose perimeter increases, while their area shrinks, are a significant amount of the total and have a peculiar behavior, with an increase of fractal dimension D(0) but no appreciable change in the degree of multifractality, defined by ΔD = D(−20) − D (20). In conclusion, our study has laid the grounds for the utilization of multifractal analysis as a promising tool for the investigation of glaciers, since their evolution over time is linked in a non-trivial way to the properties of their multifractal spectra. Currently, the study has been limited to the glaciers of Lombardy in the Italian Alps, and it could be that some of the features of the analysis reflect the orography of this region. It would be interesting to compare the multifractal spectra of glaciers from different regions of the world, to assess whether the features of the perimeter and area distributions and the correlation between the multifractal spectra and the size of the glaciers are a common feature, or instead are related to local details of the landscape and the climate of the region.

Data availability
The data that support the findings of this study are openly available at https://doi.org/10.5281/zenodo. 4301867.

Author contributions
GAD and AV conceived the study, MC and AV designed the study, analyzed results, and wrote the paper, AR processed the multifractal spectra, GAD, AS, and DM provided the glacier boundaries inventories, AR, AS, GAD, MC, and AV discussed results. All the authors contributed to the paper. be uploaded on the official regional WEBGIS portal (www.cartografia.regione.lombardia.it). AV and MC gratefully acknowledge discussions with D Brogioli.