Multifractal analysis of the perimeters of glaciers in the Svalbard Archipelago

The investigation of the evolution of glaciers largely relies on the characterisation of extensive quantities like their mass, area, and perimeter. In this work we use fractal and multifractal analysis to investigate the non-extensive structural properties of the perimeters of glaciers in the Svalbard Archipelago. We show that the perimeters of the glaciers exhibit a fractal structure with a fractal dimension D f ≃ 1 . 25, independently from the area of the glaciers. The investigation of the multifractal properties of the perimeters shows that small glaciers exhibit a more pronounced multifractal structure, as witnessed by the larger range of generalized dimensions D q needed to characterise them. The range ∆ D q of generalised dimensions required to characterise the multifractal perimeter of a glacier exhibits a power-law dependence with exponent − 1.2 from the area, and represents a non-extensive parameter able to grab effectively the dependence of the multifractal structure of the perimeters on the size of glaciers. The comparison with similar results obtained in a previous study performed on glaciers in the Lombardy region of the Italian Alps confirms the robustness of the analysis performed, which does not appear to be affected by the morphology of the substrate or by climate conditions.


Introduction
Glaciers have a significant local and global impact on the environment, as they affect many hydrological, physical, chemical and biological processes [1,2].In high latitude areas, the importance of glaciers is also determined by the fact that meltwater influences terrestrial and marine ecosystems and fjord circulation [3].Moreover, glacier retreat leads to significant topographic changes, such as increased fjord length and glacier forefield area [4,5].In high altitude regions, glaciers are the main source for water availability [6] playing a crucial role in regional water cycles in a variety of river systems around the world [7,8].Consequently, they influence also human activities, such as agriculture and energy production through the amount of water available for irrigation or hydroelectric dams [9].On a global scale, glacier melting is responsible, at least in part, for sea level rise [10].For instance, the collective volume of water stored in Arctic glaciers has the potential to rise global sea level by ≃0.3 m if melted completely [11].Although this is a much smaller volume than the sea-level equivalent of the ice sheets on Greenland and Antarctica, the contribution from Arctic glaciers (−213 ± 29 Gt y −1 ) currently accounts for about one third of the eustatic sea level change [12].Studies on the evolution of glaciers show that the rate at which they have lost mass during the early 21st century has increased compared to earlier periods: the current rate of loss of mass is four times the one of the period 1851-1900, three times the one of 1901-1950, and twice the one of 1951-2000 [13].In addition, the global climate observing system includes the area, elevation change, and mass change of glaciers among the essential climate variables (ECVs), i.e. variables that reliably characterize the Earth's climate [14].
For all these reasons, it is important to investigate the dynamical evolution of glaciers through the characterization of relevant parameters, which is usually accomplished using the three ECVs mentioned above, but also using other quantities such as volume and perimeter length.Almost all of the parameters currently used for the characterisation of glaciers are of extensive nature, as they scale with the size of the system.In this work we investigate the effectiveness of non-extensive parameters in achieving a quantitative characterisation of the geometrical structure of the glaciers in relation to their size.To attain this, we apply the fractal and multifractal formalism to the perimeters of glaciers.The fractal formalism has been widely used to investigate the geometric properties of natural systems in a variety of scientific fields, including biology and medicine [15].In the field of Earth Sciences, methods based on the fractal and multifractal formalism have been largely used for the investigation of landforms [16], the study of river networks [17,18], the measurement of coastal perimeters [19], and even recently glaciology [20].
The application of the concepts of fractal and multifractal analysis to the perimeters of a sample of Alpine glaciers in Italy [21] has recently shown that they exhibit a fractal structure with fractal dimension D f ≃ 1.2 and that their multifractal spectra differ significantly as a function of their size.In particular, the range of multifractal exponents needed to describe the structure of the perimeters of small glaciers (i.e.those with area < 10 5 m 2 ) is systematically larger than the one needed for wider glaciers and represents an effective non-extensive statistical property for the characterisation of their morphology.Since glaciers are formed through the accumulation of snow metamorphosed into ice over a solid bedrock substrate, a legitimate question is whether the fractal and multifractal properties reported for the Italian Alps (IA) glaciers could be significantly affected by the orographic properties of the region or can largely be considered a property of the glaciers themselves.To shed light on this fundamental question, we performed a fractal and multifractal analysis of the glaciers in the Svalbard Archipelago (SA), which is characterised by a different orography and climate with respect to the IA.
In the following we will first recall some well established concepts and methods of fractal and multifractal analyses.We will then proceed to the application of these concepts to the investigation of the perimeters of the glaciers in the SA.

Concepts of fractal analysis
Benoit Mandelbrot, considered to be one of the fathers of fractal geometry, was reluctant in defining the concept of fractal unambiguously, preferring to rely on examples [22].Other authors have tried to find a satisfactory definition, but without success.Among them, Kennet Falconer was of the idea that the definition of fractal should be regarded in the same way of that of life, that cannot be expressed with a hard-and-fast definition, but relies on a list of features that a living being should possess [23].Although a simple definition of a fractal system is difficult, a fractal should exhibit some form of self-similarity, that is, its details should resemble the whole part.These observations suggest that the properties of a fractal system can be described starting from the concept of self-similarity.

Box-counting dimension
Similarity relations for a set are usually not known in advance.Let us imagine, for example, that we want to measure the length of a river, the perimeter of a glacier, or the coastline of an island [19].These sets have an irregular shape, and every time one chooses a scale on a map, the details on smaller scales are lost.A further problem is represented by the fact that the length of sets like these turns out to be ill-defined.Indeed, if one would attempt a direct measure of one of these sets, this measure would strongly depend on the length of the yardstick used as a unit length to perform the measurement.The complexity of these systems, and the way their measurement varies with the length of the stick used, are well characterized by the box-counting dimension D B .Let us imagine creating a grid composed of equal squares of side ε and superimposing it to the set.We define a measure on each grid based on the number of squares N ε that contains part of the set to be measured.We perform this same procedure as ε varies for multiple grids.The investigation of the relation between N ε and ε for natural systems like those cited above leads to a power-law relation N ε ∝ ε −DB .We can then define the box-counting dimension as the power-law exponent: where C is a constant.The similarity in this case holds in the statistical sense and depends on the measure that is used.The advantage of the box-counting dimension is that it can be measured directly on physical systems of various nature, without having any particular prior information about the system under investigation.

Hausdorff-Besicovitch dimension
By introducing the box-counting dimension we are able to assign a fractal dimension to a large variety of sets, but this does not ensure that N ε (ε) follows a power-law.To solve this problem, we use the basic idea of the box-counting process and cover the set with a finite (or infinite countable) number of sets U i of size ε.
We endow these sets with a measure h(ε) = γ(d)ε d where γ(d) depends on the choice made on the measure, and on the covering set, and d is a positive real number called the dimension of the measure.The set that has been covered will have measure This measure behaves differently as ε tends to zero, depending on the value of d, and has a critical dimension D HB called Hausdorff-Besicovitch dimension, such that for values of d larger and smaller than D HB the measure M d either tends to 0 or to +∞: The Hausdorff-Besicovitch dimension is generally more complex to derive than the box-counting dimension, but it relies on a more rigorous definition, which can be applied to any set.We also note that equations ( 2) and (3) imply a definite asymptotic behaviour for N(ε): According to this relation, to find D HB we can, instead of using the definition directly, study the asymptotic behavior of N(ε): With the introduction of this new dimension, we can rely on the definition of fractal provided by Mandelbrot, which states that a fractal is by definition a set having a Hausdorff-Besicovitch dimension strictly greater than its topological dimension.This definition, while not definitive for the reasons addressed above, highlights the importance of the Hausdorff-Besicovitch dimension in fractal geometry.Usually for the fractal dimensions of a set the following inequality is valid where D T is the topological dimension and D F the fractal dimension referred to.For example, a fractal line (such as a Koch curve [24]) will have dimension D F between 1 and 2, while a fractal surface between 2 and 3.

Concepts of multifractal analysis
It is not always sufficient or possible to adequately characterize a fractal set through a single number, because sometimes more fractal dimensions are needed, if not entire functions.These functions that generalize the concept of fractal dimension are called multifractal spectra, and they describe the extent to which the fractal dimensions of the set are present on the set itself.

Copper minerals
An intuitive practical example provided by Mandelbrot [25] showing the need of more than one fractal dimension, is the distribution of copper on the globe.Copper ores can be classified according to the percentage of the element copper (Cu) present in them.Minerals containing larger amounts of copper are called high-grade copper, while those containing little compared to the others are low-grade.The high-grade ores are the rarest and are found only in a few regions of the world with a certain distribution.Suppose we look at one of the regions where high-grade minerals are found.In this region we can assume that the minerals are distributed unevenly and that they follow the same statistical distribution present on a larger scale.By dividing this region into subregions and looking at the ones in which high-grade copper is concentrated, we can carry out the same reasoning and continue iteratively.In this way the high-grade copper will be distributed on a self-similar fractal set in the narrow sense, which will possess its own fractal dimension.We can perform a similar reasoning for a lower grade of copper ore by finding that it too will be distributed on another fractal set with its own dimension.In this case the fractal dimension will be larger, since the fractal set will occupy a larger amount of space, due to the fact that copper of a lower quality is less rare.As a result, the characterization of the distribution of copper ores in relation to their quality requires the introduction of a spectrum of fractal dimensions.

Spectrum τ (q)
A multifractal spectrum that can be introduced naturally by generalizing the Hausdorff-Besicovitch dimension is the spectrum of the mass exponents τ (q).The HB dimension does not provide information about the internal structure of the fractal.In analogy with equation ( 2), we find a covering of the set S on which we define a measure.Each set S i of the covering has a probability of occupancy µ i of the ith cell, normalised so that ∑ i µ i = 1.We introduce an overall measure for the set S, similar to that in equation ( 2), but in which we introduce moments of order q acting on the distribution µ i : where The exponent d that allows to obtain a finite measure M d (q) in the limit of vanishing ε is spectrum of the mass exponents τ (q), which is a function of the order q of the moment of the probability µ i .The advantage of introducing moments of order q of µ i in equation (7), lies in the fact that as q varies, the contribution of the individual probabilities µ i to the total measure M d (q) is weighted in a different way.For example, a large value of q would yield µ q i ≫ µ q j in the case where µ i > µ j , i.e. the trend of N(q, ε) will be influenced mostly by the largest values of µ i .In the case of negative values of q with large modulus, the reasoning is reversed.As in the case of the Hausdorff-Besicovitch dimension we can write the power-law relation and find the critical dimension explicitly via Choosing q = 0 we find the Hausdorff-Besicovitch dimension, in fact N(q = 0, ε) = N(ε) and thus τ (0) = D HB .In the case q = 1 we get N(q = 1, ε) = ∑ N i µ i = 1 leading to the conclusion that τ (1) = 0.

Spectrum D q
We introduce the spectrum of generalised dimensions D q , the one of interest for this work, which can be defined from the spectrum of mass exponents τ (q): Let us briefly recall some of its properties to explain the reasons to adopt this spectrum to investigate the multifractal properties of the perimeters of glaciers.A first important property of D q is that D 0 = D HB .For a nonfractal set, the fractal dimension must be equal to the Euclidean dimension, and the spectrum D q solves a similar consistency problem for multifractal sets: a set with uniform internal structure, which is not a multifractal, will have a constant spectrum equal to the fractal dimension of the set D q = D F ∀q.In this case, a set S of fractal dimension D F , will be covered by N = ε −DF boxes of side ε.If the associated multifractal measure is uniform, since it has no complex internal structure, the probability of each covering set will be µ = ε DF .Thus the normalization of the probability According to these considerations, the following relation holds: and then using the definition provided in equation ( 10) Note that if instead of a generic fractal with uniform structure, we had considered a nonfractal set with Euclidean dimension D E , the reasoning would have remained unchanged by replacing D F with D E .
We can then state that the spectrum D q for non-multifractal sets has a constant value D F .Defining the aperture of the spectrum D q as ∆D q := D(−∞) − D(+∞), we can use this range of generalised dimensions as a parameter suitable to quantify the extent to which a studied set exhibits a multifractal structure [21].

The study area
Svalbard is a high Arctic archipelago located north of mainland Norway between Greenland and Novaya Zemlya (75 • − 82 • N, see figure 1).It is covered for about 60% by glaciers, including valley and cirque glaciers, ice fields and the large ice caps on Nordaustlandet [26], corresponding to about 10% of the glacier area in the Arctic, outside the Greenland ice sheet [3].The total volume of Svalbard glaciers is ca.6200 km 3 [27], corresponding to a potential sea-level equivalent of 1.5 cm [3].About 15% of all glaciers on Svalbard by number and as much as 60% by area [4] are tidewater glaciers (i.e.terminating into fjord or ocean water), while the remain glaciers terminate on land.
The largest island, Spitsbergen (37 500 km 2 ), is dominated by steep mountains surrounded by low flat terrain around the coastline and is covered for about 60% by glaciers.The two largest ice masses in Svalbard, Austfonna (8000 km 2 ) and Vestfonna (2450 km 2 ) are located on the island of Nordaustlandet (14 500 km 2 ).The two smaller islands, Barentsøya (1300 km 2 ) and Edgeøya (5100 km 2 ), are dominated by plateau-type terrain [28], and are covered by 2800 km 2 of mainly smaller ice caps at lower elevations.
The SA is characterized by a relatively warm and variable climate as compared to other regions at the same latitude.The winter sea ice cover of the Arctic Ocean to the North limits moisture supply, while in the region South of Svalbard cyclones gain strength as storms move northward [29].These geographical and meteorological conditions make the climate of Svalbard not only extremely variable spatially and temporally, but also sensitive to deviations in both the heat transport from the south and the sea ice conditions to the north [30].
The retreat of glaciers in the SA is influenced by several key factors.Rain-on-snow (ROS) events have become more frequent due to a warming climate, impacting the glaciers' mass balance and leading to varying interpretations in ROS climatology across datasets [31].Specifically, the southern and western coastal areas of Spitsbergen experience the highest ROS frequency during early winter.Additionally, glacier mass balance assessments using laser altimetry ICESat-2 data indicate an accelerated glacier mass loss in Svalbard, with the Northwestern Spitsbergen experiencing the most significant mass loss rates [32].This increased mass loss is attributed to rising temperatures and an increase in glacier surge events.An additional factor affecting the retreat is determined by the peculiar geometry of several glaciers in the SA, which are frequently bound, at least partially, by a coastline.New high-resolution calving front dataset for marine-terminating glaciers in Svalbard revealed widespread calving front retreats over the past four decades, with complex patterns of glacier surging events and seasonal calving cycles identified [33].The melting of glaciers generates glacial groundwater, which releases methane from deep beneath the glaciers and permafrost, contributing to atmospheric greenhouse gas emissions [34].Seasonal variations show that in summer, deep methane-rich groundwater is diluted and partly oxidized by shallow oxygenated groundwater, reducing potential methane emissions.However, in winter, reduced subsurface methane oxidation due to freezing could lead to increased methane emissions, with future changes expected in the ratios of groundwater sources due to warming climate impacts .These findings highlight the significant impact of climate variability, increased temperature, and glacier surging events on the accelerated retreat of Svalbard's glaciers.

The dataset
The dataset contains the perimeters of glaciers of the SA in shapefile format.This dataset is part of the most detailed Svalbard land covering map dataset, called S100 [35], which corresponds to the map series Svalbard 1:100 000.These glacier perimeters were manually digitalized on aerial photos mostly acquired in 2008-2012 period (78% of the total) with a planimetric accuracy of 2-5 meters.Only few areas are based on other aerial imagery from 1961 and 1970 (2.5%) and 1990 (17%), located in the southern Spitsbergen (i.e.Heer Land, Nathorst Land, Wedel Jarlsberg Land and Torell Land), the biggest island of the archipelago.Due to the relevant dynamics of glaciers, especially in the terminus, these portions were successively modified and updated by the operators (over the areas based on 20th century images) using Landsat 7 ETM+ images.
The choice of this dataset was suggested by the fact that the orography and climate of the region are significantly different from those of the IA, where the first study on the multifractal properties of glaciers was conducted [21].Additionally, among the existing Svalbard glacier inventories [26] this is the only recent based on aerial images at a high spatial resolution, and not on satellite images at tens meter resolution [36].These features of the dataset, along with the fact that the glaciers were manually digitalized, are similar to the ones of the previously studied Italian glacier inventory [21].
The inventory is divided into three shapefiles containing 3020 small glaciers with an area lower than 5 km 2 (they are numerically 97% of the total, but cover only 3% of the total area), 78 medium glaciers with area between 5 km 2 and 50 km 2 (numerically 2.5% of the total, and cover 3% of the total area), and 19 large glaciers with area higher than 50 km 2 (numerically 0.5% of the total, but cover 96% of the total area), for a total of 3117 glaciers.The smallest glacier (ID: 35 084) has an area of 1.2 ×10 −4 km 2 , while the largest (ID: 31 552) has an area of approximately 8197 km 2 .By comparison, the largest glacier in the IA dataset (Adamello Glacier) had an area of about 20 km 2 (data from 2003).This shows a further difference between the two groups of glaciers: a limited number of glaciers in the SA are significantly larger than the IA ones.
We proceeded to visualize the dataset (figure 1) noting that the glaciers are not all simple polygons, but may have holes (called nunataks and corresponding to bedrock outcrops or rock windows inside a glacier area).Therefore, the perimeter of some glaciers, generally the largest ones, consists of an outer perimeter and one or more inner perimeters.The nunataks are generally increasing in size and are becoming very frequent due to the ongoing glacier retreat and thinning [37].The rock has a lower albedo (or reflectivity) than the glacier ice.Consequently, mountain slope areas external to the glacier surface, with no snow cover during summer, absorb more solar radiation and then warm faster than the glacier surface during the day.The surface energy balance of lateral slopes results higher and influences the nearby glacier surface through increased outgoing longwave radiation flux and air temperature.The result is a promoted snow/ice melting in the areas surrounding the nunataks causing an enlargement of the nunataks [38].As a consequence of this process, the glacier could undergo a fragmentation into two or more ice bodies [7].As these nunataks were considered to be significant for the structure of the glacier itself, they were included in the final multifractal analysis.
To obtain a dataset of glaciers of areas comparable to those of the IA, we decided to remove from the analysis the five smallest glaciers in the entire SA dataset and the 19 largest.The removal of the smallest glaciers does not pose any particular problem, as other slightly larger glaciers are present abundantly in the dataset.By comparison, the smallest glacier in the SA dataset has an area of 6.68×10 −4 km 2 , compared with 4.6 ×10 −4 km 2 for the smallest IA glacier (i.e.Pizzo del Ferro Inferiore II, data from 2007).Instead, the removal of the larger glaciers results in the loss of glaciers that significantly differentiate the Svalbard dataset from the IA dataset.These removed glaciers are as much as two orders of magnitude larger in size than the largest Italian glaciers, and this would have resulted in a larger range of areas in which correlation with the D q spectrum could have been studied.The largest SA glacier included in the analysis has area of 36 km 2 versus about 20 km 2 for Adamello Glacier, the largest glacier in IA dataset.However, the perimeters of these large glaciers removed from the SA dataset are almost entirely dominated by the coastline of SA.As such, they represent level lines, and are not a meaningful term of comparison with the glaciers of the IA dataset.Comparing analyses between these two datasets covering a comparable range of areas, therefore, remains a good way to test whether the correlation between area and multifractal D q spectrum reported for alpine glaciers [21] are region-dependent, or a property of glaciers in general.

Validation of the dataset
When dealing with the production of glacier inventories through a manual digitalization on aerial photos, inaccuracies may occur due to operator's misinterpretation.These depend upon the image resolution and the meteorological and environmental conditions at the time of acquisition (i.e.cloud-and snow-cover, presence of shadows and debris), hampering the discrimination between pixels corresponding to glacier ice and those corresponding to rock.
Since image resolution influences the accuracy of glacier mapping, the final planimetric precision value was assessed considering the uncertainty due to the sources (aerial photos).Following Minora et al [39], the linear resolution error should be half the resolution of the image pixel, i.e. in our case ranging from 1 to 2.5 m [40].This error may be too low for debris-covered glacier areas because the glacier limits are more difficult to distinguish when ice is covered by debris.Therefore, we set the error for debris pixels to be three times that of clean ice.
As our study focuses on the effectiveness of the application of multifractal analysis on glacier perimeters, we carried out a further investigation on possible artefacts determined by the manual mapping of glacier outlines.Indeed, the vertices of the polygons constituting the shapefiles are plotted one at a time by the operator, who reproduces as accurately as possible the perimeter of the glaciers from raster aerial images.This process could lead to different perimeter accuracy, depending on the glacier and on the operator who performed the task.In the undesirable case of the presence of an accuracy depending systematically on the glacier sizes, the multifractal analysis and any correlations related to glacier area would be affected, and this would spoil the reliability of the analysis.For this reason, we performed a further analysis investigating systematically and quantitatively the relation between the accuracy used for mapping outlines of a glacier and its area.We define the perimeter accuracy as the number of perimeter vertices of a glacier in a shapefile, divided by the length of the perimeter.The perimeter accuracy plotted as a function of the area does not exhibit any significant correlation (figure 2).This result provides a validation of the dataset, which can be used reliably to determine the fractal and multifractal spectra of glaciers.It may be interesting to note that the results in figure 2 are clustered into three main groups.This can be due probably to three different operators, each with a different average precision, or to three sources of data.This latter seems the most probable cause since the SA aerial photos were acquired in 2008-2012 period (78% of the total), 1990 (17%) and 1961-1970 (2.5%).

Rasterisation
To perform multifractal analyses, it is necessary to use an outline of the perimeter of glaciers in a raster format; this is because FracLac, the software adopted for the multifractal analysis, requires raster images as an input.The process of conversion from vector to raster image is called rasterisation, and can be performed using the dedicated algorithms present in the GDAL library [41] of QGIS, an open source GIS (Geographic Information System) application [42] made available by the Open Source Geospatial Foundation [43].Rasterisation was performed for each glacier in the dataset adopting a resolution of the raster images such that each pixel corresponds to a square of side 0.5 m in the real world.This choice is dictated by the need of maintaining the same precision used in analyses of alpine glaciers.

FracLac
FracLac is an ImageJ plugin that allows fractal, multifractal, and gap analysis [44].We chose to make use of this plugin for multifractal analyses because it enables a reliable optimisation of the parameters, and it has been used effectively to perform the analysis of the glaciers of the IA which motivates this study.

Algorithms of multifractal spectra
There are a variety of algorithms to implement fractal and multifractal analysis, each with its own peculiarities, advantages and disadvantages.Fractal dimensions and multifractal spectra derived from different algorithms do not necessarily coincide, an aspect that requires some precautions when comparing results obtained with different methodologies.Unlike a measurement, for example of length, where one can compare results by means of error estimates on individual measurements, comparing two fractal dimensions of different sets derived from different procedures is not straightforward.Although the literature regarding possible algorithms for deriving dimensions and spectra is abundant, the one regarding comparisons of the results obtained with different methods is less conspicuous.Therefore, special attention was dedicated on reproducing the analyses through a methodology as similar as possible to the one adopted for alpine glaciers.The box-counting dimension lends itself well to an algorithmic implementation, and following the same basic idea it is easy to produce the τ (q) spectrum.The image to be analyzed is covered with grids of different side ε, as in the case of box-counting, while the probability of occupancy P i associated with each grid square is where m i is the number of black pixels (i.e. according to the conventions belonging to the perimeter being analysed) in the square region, and M is the number of total black pixels in the perimeter.According to this definition P i obeys the normalisation condition ∑ i P i = 1.For each fixed q and for each grid size ε one can compute where N is the number of square regions containing at least one black pixel.As a result, one can obtain the multifractal spectrum τ (q) by adopting a procedure very similar to the box-counting one: Similarly to what occurs in box-counting, the value of τ (q) is obtained as the exponent of the power-law dependence between N(q, ε) and ε at each fixed q.
FracLac, in its implementation, follows the procedure laid out in [45], which relies on the determination of the multifractal spectrum f(α) [24,46], where α is the Lipschitz-Hölder Exponent that characterises the dimension of one of the interwoven fractals sets that form the multifractal system and f(α) is the fractal dimension of the set of exponent α.As before, grids with different gauges are used and the P i measurement is calculated for each grid square at fixed q.Then, the measure defined as is also calculated for each square, and the quantity is calculated for each gauge.The value of α is derived as the exponent of the power-law relation between A ε and ε.Repeating the process for each q leads to the determination of the function α(q).For values of f(α(q)) the procedure is similar.One starts from the quantity and the values f(α(q)) represent the exponents of the power-law relation between F ε and ε.Once that the multifractal spectrum f(α) is known, the spectrum τ (q) can be obtained through the Legendre transformation τ (q) = f[α(q)] − qα(q), and the spectrum D q through equation ( 10).

Scan multifractal
FracLac was used to produce multifractal spectra data from raster images of the glacier outlines.Due to the batch process capabilities of FracLac, no special codes had to be implemented to automate the process.The results were used to produce plots of the multifractal spectra and to perform some tests using the Python programming language [47].Before performing the analyses with FracLac, a tweaking of the optimal parameters was needed.In fact, the analyses depend on the parameters set to accomplish them, and with this idea in mind we tried to replicate as much as possible the ones used for the analyses performed on the alpine glaciers.

Analysis and discussion
This section presents and discusses the results obtained from the analyses performed on glaciers in the SA, which are compared with the ones on the multifractal analysis of the perimeters of Italian glaciers [21].In addition to the multifractal analysis, the set of glaciers has been characterised by determining the probability distributions of areas and perimeters and the relationship between area and perimeter.

Distributions of area and perimeter
The basic geometric properties of the dataset were characterised statistically by determining the cumulative probability distribution (Cdf) and the probability density (Pdf) for the perimeters and areas of the glaciers.We followed the approach proposed by Newman [48], which consists in first deriving the cumulative distribution, which can be defined for each point in the dataset without binning the data.Denoting with P(x) the cumulative distribution, where we are using the convention that it accounts for the probability of a value larger or equal to x, and denoting with ρ(x) the probability density, the following relations hold between the two quantities [21]: The cumulative distribution was modeled and fitted with the empirical function [21]  where x 0 and x m denote the smallest value and the median of the dataset, respectively.This function was chosen because in the case of the IA data it reproduced the Cdfs accurately.The probability density was derived combining equations ( 20) and ( 21): Observing figures 3(a) and (b), it can be appreciated that for large values of the area and perimeter the cumulative distributions exhibit a power-law dependence.The γ values found by the fit for area and perimeter distributions are 2.01 and 2.61 , respectively.For IA glaciers, on the other hand, the following exponents γ were obtained: An interesting feature, reported also for the glaciers in the IA, is that the perimeter distribution (figure 3(c)) exhibits a marked peak, which is not apparent in the area distribution (figure 3(d)).The peak in the perimeter distribution shows the presence of a preferential length of the perimeter of glaciers in the dataset, which is not accompanied by the presence of a characteristic area.

Area-perimeter relationship
In the case of a fractal set, a definite relation holds between the area of a plane region and its perimeter.
A plane region with area A bounded by a fractal contour of length L and fractal dimension D F follows power-law relation: where D E is the Euclidean dimension [24].
We can then derive the average D F value for the entire glaciers dataset from the power-law exponent of the Area-Perimeter relation (figure 4).The average fractal dimension of the perimeter of the glaciers in the SA dataset is D F ≃ 1.26, close but larger than the value 1.21 found for the IA.We can then state that Svalbard glaciers tend to have more irregular perimeters.

Multifractal analysis
To investigate the relation between the multifractal spectra of a glacier and its area we processed the D q spectrum for all the SA glaciers (figure 5(a)).The spectra obtained for glaciers in the SA exhibit features similar to those in the IA dataset (figures 5(b)-(d)).In both cases, the spectra D q superimpose at values of q > 0, indicating the lack of effectiveness of this region of the spectra in discriminating between glaciers of different area.In the study on the glaciers of the IA dataset a negative correlation was reported between the area of the glaciers and the value of the D q spectra at negative values of q.This correlation is clearly visible from the plots in figure 5 for q ∈ [−20, −5], both for the spectra of the IA and SA glaciers.This result provides a first confirmation of the effectiveness of the D q spectrum in discriminating the different structure of the perimeters of glaciers of different area.
To achieve a better quantitative characterisation of this correlation, in figure 6 we plot the relation between the value of the D q spectrum at the three values q = −20; 0; +20 as a function of the area of the SA glaciers.From the figure, one can appreciate the presence of a negative correlation between the D(−20) value and the area.In addition, the absence of correlations between the area and the D(0) and D(+20) values can be observed, a result compatible with the one reported for alpine glaciers.Recalling that the value D(0) is equal to the fractal dimension of the set, from figure 6 it can be seen that the average value of the fractal dimension is around 1.2 − 1.3, a value fully compatible with the one found through the area-perimeter relation in the previous paragraph (figure 4).Thus, from the point of view of fractal analysis all the perimeters are equivalent.More in detail, the fractal dimension characterizes the scale invariant properties of a perimeter, that is, a small part of the perimeter when enlarged statistically exhibits the same corrugation as a larger part of it.Therefore, when looking at a small part of the perimeter, in the absence of a spatial reference it is impossible to infer what might be its spatial extent, because in the presence of scale invariance the perimeter does not exhibit a characteristic length scale.On average, all glaciers in a population exhibit a fractal structure of the perimeter with the same fractal dimension, regardless of their spatial extent.This result shows that, according to the fractal analysis, all the glaciers in the dataset have a similar scale invariant structure, which does not depend on their size.
The aperture ∆D q := D(−20) − D(20) of the D q spectrum allows to achieve a quantitative characterisation of the correlation between the multifractal spectrum and the area (figure 7).The relationship follows a power-law dependence ∆D q ∝ A β , with a characteristic exponent β = −1.2close, but significantly different, from the exponent β = −1.0 reported for the IA glaciers.This observation allows us to conclude that small glaciers exhibit in general a more marked multifractal structure than larger glaciers, as witnessed by the larger range of D q values needed for their description.In other words, the perimeters of  small glaciers are statistically more variably corrugated than those of large glaciers, which are closer to purely fractal perimeters.This implies that the perimeters of small glaciers exhibit a wider variety of conformations than the larger ones.The transition appears to be more marked for the glaciers in the SA, due to the more pronounced steepness of ∆D q (A) as a function of the area.The fact that the values of the D q spectra for negative q values are affected largely more than those at positive q provides an indication of the structural differences between the perimeters of small and large glaciers.In fact, D q at negative q is affected more by the smallest µ i weights, i.e. those grid squares that contain few black pixels in comparison to others.

Conclusions
In this study, we compared the fractal and multifractal properties of glaciers in the SA with those of Alpine Italian glaciers in the Lombardy region.The characterization of the area and perimeter distributions of the dataset shows that the distribution of perimeters exhibits a marked peak, similarly to the one found for the IA glaciers.The area and perimeter of glaciers in the SA exhibit a power-law relationship, leading to an average fractal dimension of 1.26, similar but slightly larger than the one found for the IA data of 1.21.This allows us to conclude that the perimeters of the Svalbard glaciers are more corrugated than the ones of the Alpine Italian glaciers.This result can be explained by the fact that Svalbard is currently among the fastest warming regions on Earth [30] leading to a general glacier shrinkage.As reported in a previous study of the evolution of the multifractal properties of IA glaciers [21], retreating glaciers feature a growing perimeter, which becomes increasingly irregular and then more corrugated.In addition, it was found [36] that the retreat of a sample of 400 glaciers in south and west Spitsbergen was greater after 1990, while glacier area change rates were higher in the decades before 1990, corresponding to more lateral wastage in the early period (affecting the shape of glacier outlines).The authors suggest that the spatio-temporal analysis of change signals can be achieved more effectively using a typical lenght scale rather than the area change.This falls in line with previous investigations [49] and may enhance the use and incorporation of glacier inventory changes, for example, into numerical models [50] and/or temperature reconstructions [51].Within this context, our glaciological application of fractal and multifractal analyses can represent a further attempt to understand the change signals of glaciers in the ongoing climate change and water impact.
The multifractal analysis of SA glacier perimeters showed the presence of a negative correlation between the aperture ∆D q of the D q spectrum and the area of the glaciers, compatible with a power-law relation ∆D q ∝ A β with exponent β = −1.2.The robustness of the correlation analysis method based on the aperture of the D q spectrum allows us to conclude that for the Svalbard dataset, glaciers with small area (< 5 km 2 ) have perimeters with a more marked multifractal character than the one of larger glaciers.This result is in agreement with the one found for glaciers in the IA [21].Since the two datasets were purposely chosen from regions with different orographic characteristics of the bedrock and climate conditions, we conclude that the more marked multifractal character of small glaciers is a feature that does not depend on the orography of the region and, as such.represents a characteristic property of the glaciers themselves.Finally, this multifractal analysis confirms that the range of generalized dimensions D q needed for the description of glacier perimeters represents a meaningful non-extensive parameter for the statistical characterization of the structure of glaciers, negatively correlated to their area.
To further confirm the conclusions of this work, it would be desirable to carry out similar studies using data from other regions of the world.Moreover, studies on the temporal characterization of the multifractal properties of glaciers could be useful in formulating physical models that explain the dynamic evolution of glaciers.More in detail, a physical model that predicts an increase or decrease in the area of a glacier would also have to explain why on average there is an increase or decrease in a non-extensive property like the multifractal structure of perimeters.

Figure 1 .
Figure 1.Map of the entire dataset of glaciers in The Svalbard Islands: large glaciers in green, medium glaciers in red, and small glaciers in blue.

Figure 2 .
Figure 2. Dependence of perimeter accuracy on glacier area: no significant correlation can be appreciated.The data are clustered into three main groups, which correspond likely to three different sources of data.

Figure 3 .
Figure 3. Distributions of perimeter and area of glaciers in the Svalbard Islands.(a) cumulative distribution of perimeters; (b) cumulative distribution of areas; (c) distribution of perimeters; (d) distribution of area.The solid lines represent the best fit with equations (21) and (22).

Figure 4 .
Figure 4. Relation between area and perimeter of glaciers in the Svalbard Islands.The power-law exponent is DE DF .

Figure 6 .
Figure 6.Plots of D(−20), D(20) and D(0) as a function of glacier area.A negative correlation is particularly noticeable for the values of D(−20), while the fractal dimension of glaciers D(0) does not exhibit any significant correlation.

Figure 7 .
Figure 7. Power-law relation between ∆Dq and the area of glaciers.The power law exponent is 1.2.