Madograms help to quantify mountain frontal zones—An approach towards comparative spatial analysis of complex landforms

This article explores the value of madograms to characterize the geomorphology of frontal belts of mountain ranges. Madograms are a geostatistical, variogram‐related concept, allowing for the analysis of altitude variability as a function of distance. After simple idealized digital elevation models were tested to reveal characteristic madogram shapes for various landform configurations, madograms were calculated for three frontal zones of different origin and morphology, located in the Sudetes range (Central Europe). Isopleth madograms and Ward's clustering method were applied to investigate similarities and differences within and between frontal zones. The results proved consistent with the accumulated, largely qualitative knowledge about geomorphology of the frontal zones and confirm the usefulness of the method, which allows for more insightful quantification of the morphology of mountain frontal zones, recognition of their internal spatial structure, and discrimination between fronts of various morphology and origin. Various recommendations for future users are also offered.

problematic delimitation of frontal zones, especially in low-activity areas, where the fronts have suffered from severe degradation (slope flattening) and erosional dissection. Likewise, the choice of parameters employed in the analysis remains arbitrary, whereas methods and techniques of comparative analysis should respect spatial autocorrelation of the input data. The latter limits applicability of classic statistical methods such as tests of significance, to evaluate differences between morphometric variables at both intra-and inter-front scale. One way to overcome these limitations is to explore geostatistical tools and methods for mountain front analysis.
Geostatistical methods, including semivariogram-based ones, have been long and widely applied to assess spatial variability of environmental processes and phenomena, referred to as "environmental variables" (Hengl, 2009).
In geomorphology, the application of variogram-based analysis to landforms or their assemblages seems much less frequent. Mulla (1988) adopted it to compare areas of different topographic patterns and found it useful for identifying the correlation scale and orientation of relatively small, non-periodic landscape features. The morphological differences between cockpit karst landscapes and other karstic areas in Jamaica were reflected in the different values of theoretical variogram models (Lyew-Ayee, Viles, & Tucker, 2007). The latter were used to compare morphometric differences between landform classes distinguished on the basis of fuzzy k-means classification (Hengl, 2009), although an inverse approach, in which variograms were incorporated into the classification procedure, can also be encountered (Azañón, Delgado, & Gómez, 2004). The measures of spatial continuity, whether variograms or their modifications, were also used to assess surface roughness/texture, most frequently for fineor medium-scale landforms such as scree slopes or gravel river beds (Huang, Atkinson, & Wang, 2018;Trevisani & Cavalli, 2016;Trevisani, Cavalli, & Marchi, 2009;Trevisani & Rocca, 2015).
The above review shows that geostatistical tools have hardly been used in tectonic geomorphology and, more broadly, to provide insights into the characteristics of larger and complex landforms. Hence, the main objective of the research presented below is to develop and test an analytical procedure based on madograms that would enable new insights into the topography of mountain frontal zones. The analysis is based on the assessment of altitude variability as a function of distance and direction. Our working hypothesis states that madogram-based analysis of altitude variability allows for more insightful quantification of mountain frontal zones, recognition of their internal spatial structure, and discrimination between fronts of various morphology and origin. The key questions to be answered include: (a) which specific features of the method should be particularly considered in the context of mountain frontal zones; and (b) how to process and present madograms so that quantification and comparative analysis can be effectively carried out. This study builds on our pilot study, where madogram analysis was part of a bigger story and found a promising avenue of research, worth further exploration (Jancewicz, Różycka, Szymanowski, & Migoń, 2020;Różycka et al., 2021). In those studies only one model setup was used, theoretical background was not discussed, and the analysis was carried out for one mountain area. Likewise, the visualization of results was limited to simple madogram curves.
As this article is primarily methodological in scope, the following structure is adopted. Section 2 includes a short terminological discussion about mountain frontal zones and problems of their delimitation, crucial for further procedures. Section 3 gives a brief presentation of madograms, and Section 4 outlines the relevant characteristics of the study areas, located in the Sudetes (Central Europe). Section 5 is divided into three subsections. First, madograms derived from simulated DEMs, used as training sets of elevation data, are generated and discussed.
Second, we have an analysis of two frontal zones delimiting the Sowie Mountains block, where different parameters of the method are tested. The focus on two frontal zones helped to assess the sensitivity of the method to identify differences between fronts of likely different origin (fault-generated vs. fault-line escarpments). Third, the fault-controlled NE frontal zone of the Sowie Mountains is analyzed jointly with the non-tectonic, denudational escarpment of the Stołowe Mountains to see whether the method is capable of making fundamental discriminations between fronts which are vastly dissimilar. In Section 6, we discuss the added value, advantages, and limitations of the method. Section 7 concludes and outlines possible avenues for further exploration.

| FRONTAL ZONE S AND THEIR DELIMITATI ON
Transitions from lowlands to uplands or mountains may be gradual or abrupt. The latter occur in the form of an escarpment or mountain (range) front, which provides a link between low-relief foreland and much higher, steeper, and usually more rugged topography. Although both these terms are in common use, attempts to define them more formally are few. Young (2004, p. 337) observed that "escarpment, or scarp, has been applied traditionally to a steep, often single slope, of considerable length, that dominates a section of landscape." Further on, he noticed that escarpments may have various origins, with some being denudational, typically formed as a result of differential rock resistance, and others produced by tectonic uplift along faults or thrusts. For the latter, the term "mountain front" is widely used and considered synonymous with "fault-generated mountain front" (thus being essentially its shortened variant), although Różycka et al. (2021) suggested that the adoption of such a meaning is unnecessarily restrictive and does not consider uncertainties regarding the origin of the mountain/foreland (M/F) junction. Moreover, the hinterland is not necessarily a rugged dissected terrain (mountains); it could also be a high-elevation plateau. Another observation is that escarpment slopes may remain relatively undissected, except for large and widely spaced canyons, whereas faulted mountain fronts are typically considerably dissected. Realizing these terminological ambiguities and trying to keep it simple, in this article we use the generic term "frontal zone" to include all marginal zones of the highlands, immediately above the M/F junction, whatever their origin and morphological details.
Defining the location and width of the frontal zone is strongly scale-dependent. At small scale, it may be approximated by a line, which can be straight or sinuous. However, this highly generalized approach is only helpful to examine the large-scale ground plan, whereas any other morphological details remain unknown. Therefore, it is more useful to consider the frontal zone as an area located immediately upslope from the M/F junction. Its morphology may be relatively simple or complex, reflecting the following key morphological attributes. These are, among others, height, slope inclination, sinuosity, distance to the water divide, density of erosional dissection, and the presence of minor landforms such as triangular facets along fault-generated mountain fronts (Wallace, 1978), or sapping cirques (amphitheaters) within some denudational escarpments in sedimentary rocks (Laity & Malin, 1985).
Delimitation of the frontal zone may be approached in various ways. Defining the M/F junction, or the base of the frontal zone (outer boundary), is usually more straightforward, especially for faulted margins, and involves tracing a distinct break of slope, where flattish terrain gives way to much more inclined surfaces. Wide valley mouths may cause problems and force some arbitrary decisions as to how large foreland embayments into the mountains should be considered. The choice is much more problematic in respect of the inner boundary, where a clear change of morphology (upper slope break, reflected by an evident increase in profile curvature) may not occur. In fact, except for denudational escarpments with caprock-built rock wall in the uppermost section, such evident changes are at best of very local occurrence.
However, even if they can be identified, the resultant line will likely be very complex and hence, the entire frontal zone will show variable width, causing problems in comparative analysis. Therefore, a different approach is followed in this study, where frontal zones are understood as 2,000 m-wide belts extending from the M/F junction upslope. Thus, in GIS analysis they are equal to buffer zones, which include various morphological elements and hence, may show significant morphological diversity. However, the adopted width of 2,000 m is not a universal recommendation, but was selected for the purpose of this study, after expert-based geomorphological assessment of all frontal zones considered in this analysis. This buffer ensured that the entire frontal zone was included, although in specific situations (narrow frontal zones, cuesta-type escarpments) this also resulted in crossing a water divide. In other geomorphic contexts a different width of the frontal zone might be a better choice, depending on the height of the marginal escarpment.

| ME THODS
Spatial variation occurs when the quantity is location-dependent, that is its values vary between observation sites.
It is a typical feature of the natural environment, just as altitude variability is a typical feature of any mountain area, even though it is manifested by significant differences in range, spatial scales, and directions. Quantitative assessment of spatial variability (or continuity), frequently based on variograms, is one of the basic tasks of geostatistics. A variogram is defined as the variance (2 ) of differences between field values Z at two locations i separated by a distance h (Cressie, 1991): The classical estimator of the variogram for a stationary field, as proposed by Matheron (1962), is: where N (h) is the number of pairs of observations separated by distance h.
The Matheron estimator is unbiased, but it is not robust-it behaves poorly if there are outliers in the data due to magnifications of the differences through squaring. Therefore, in accordance with the recommendations from various publications (e.g., Gneiting, Ševčíková, & Percival, 2012), we decided to use one of the variogram-related concepts-the madogram (variogram of order 1; Olea, 1991). This decision also allowed us to conduct analysis in units of altitude, in close interpretation reference to the concept of relative height, even if considered as a function of distance and direction. For this purpose, the squared term of the variogram, which allows for the calculation of the average absolute altitude difference ΔZ (h) for all N pairs of points separated by the distance h, computed in the horizontal plane. In this study, however, the analysis of spatial variability of altitude was based on empirical madograms, for which separation distance bins (spatial lags) h ± are used rather than exact distances: where 2 is the range of each lag. Finally, the vector of mean absolute altitude differences for every spatial lag was computed and plotted.
Madograms can be calculated as omnidirectional, typically used if the analyzed variable is isotropic, that is it does not show directional tendencies in space. If anisotropy is found, directional madograms are more useful.
The direction is then indicated and a certain angular tolerance 2 is defined such that all pairs of points whose directions of the separation vectors lie within the angular tolerance of the indicated direction ( ± ) are included in the calculations for each spatial lag ( Figure 1).
The empirical madogram is the basic tool used in this work, but not the only one. Additionally, we used two methods for further processing of the madograms. The first led to the transformation of empirical madogram vectors into a spatial form, the so-called isopleth madograms. The second allows us to search for similarities and differences in the shapes of madograms by using the hierarchical, agglomerative method of clustering-the Ward's (Everitt, Landau, Leese, & Stahl, 2011). The notion of the "shape" of the madogram refers to the course of the line that results from joining individual madogram values in successive spatial lags (also termed "madogram curve" below). However, as our study is generally methodological in nature and consists of several successive steps, a more in-depth explanation of the specific methods applied, including calculations, analysis, and further processing of the madograms, will be presented in the following sections of the text.
Morphometric properties of frontal zones considered in this study were analyzed on the basis of a 5 m × 5 m DEM, resampled from a 1 m-resolution LiDAR model in order to remove artefacts. Resampling was done using the cubic convolution method. Calculations and spatial analyses were performed using ArcGIS 10.7 (including Geostatistical Analyst, Spatial Analyst, and 3D Analyst extensions) and Statistica 13.1 software.

| S TUDY ARE A
The area selected for this study is located in the middle part of the Sudetes range (Poland and Czechia, Central Europe) (Figure 2), which represents an example of polygenetic topography, shaped by concurrent differential vertical crustal movements and rock-controlled denudation (Migoń, 2011;Szymanowski et al., 2019). Along with considerable lithological diversity of the Sudetes, these circumstances dictated the origin of numerous topographic escarpments of sub-regional extent (>10 km long). Among them are fault-controlled mountain fronts (Badura et al., 2007;Krzyszkowski et al., 1995;Różycka & Migoń, 2017), fault-line escarpments along inactive faults , and denudational escarpments, which owe their origin to the presence of resistant caprock overlying softer rock formations beneath (Duszyński & Migoń, 2015;Migoń & Zwiernik, 2006). The proximity of three frontal zones of apparently different origin and morphology, but similar length and height, provides a good opportunity to attempt the methodological study, since spatial scale becomes a factor that is not predicted to influence the results in a significant way. Otherwise, testing the working hypothesis that diverse morphologies and, indirectly, diverse origins may be reflected in madograms seems feasible and worthwhile.
The NE front of the Sowie Mountains ( Figure 2) is part of a much longer marginal escarpment of the Sudetes, which is more than 150 km long in total (Badura et al., 2007) and reflects late Cenozoic uplift along the Sudetic Marginal Fault (SMF). Geomorphic evidence for ongoing uplift during the Quaternary includes terrace staircases, truncated alluvial terraces, very low value of front sinuosity index (Krzyszkowski et al., 1995;Krzyszkowski & Olejnik, 1998;Krzyszkowski & Pijet, 1993), and displaced late Quaternary sediments elsewhere along the fault (Štěpančíková et al., 2010). The footwall is built exclusively of gneiss, although localized remnants of Carboniferous

F I G U R E 2 Study area
The SW frontal zone of the Sowie Mountains ( Figure 2) is superficially similar to the NE front in terms of sinuosity, height, and certain geomorphometric signatures, but evidence for neotectonic faulting is lacking. Its course, although it follows an old fault zone, is also rock-controlled in that the mountainous area is built of gneiss, whereas the foreland is underlain by sedimentary rocks of Carboniferous age, juxtaposed against each other along the fault. Therefore, Różycka et al. (2021) tend to interpret the SW front as a fault-line escarpment, whereas certain morphological similarities between the two fronts are due to identical lithology of the mountainous zone.
The third frontal zone considered in this study, the NE escarpment of the Stołowe Mountains (Figure 2), is very different. It represents a topographic margin of a sedimentary tableland (SE part) and a prominent cuesta (NW part). The escarpment is capped by moderately resistant Upper Cretaceous sandstones, which are exposed in 10-30 m-high cliff lines along a considerable length of the scarp (Duszyński & Migoń, 2015). The rim of the escarpment is very clear along most of its length. The middle slope truncates older Cretaceous strata and the lower slope is underlain by Permian deposits. Thus, in contrast to the fronts of the Sowie Mountains, the frontal zone of the Stołowe Mountains is built exclusively of sedimentary rocks. The escarpment is moderately sinuous in plan, with large amphitheaters (sapping cirques; Pulinowa, 1989) in the eastern part. However, fluvial incision is very limited and no large valleys cross the frontal zone.

| Experiment design
Madograms, which reflect altitude variability as a function of distance and direction, may take different, fairly complicated shapes for natural land surface configurations . Therefore, their interpretation is not intuitive since the distance relationships are defined in respect to distances between points in the dataset, without taking into account the exact locations of these points inside the area subject to analysis. To facilitate better understanding of madogram shapes and to see how various land shapes may be reflected, we used simulated, idealized DEMs which represent specific morphological configurations to be expected along a frontal zone. In addition, this exercise helped to define the parameters of the method, to be applied in the next stage. In generating simulated DEMs, the following parameters were selected: 1 m × 1 m spatial resolution, 2,000 m × 2,000 m dimensions, 400 m elevation difference in respect to the frontal zone base (equal to mean slope of 11.3°) ( Figure 3).
These parameter values are not very much different from the real-world characteristics of mountain fronts in the Sowie Mountains, considered in our pilot study (Jancewicz et al., 2020;Różycka et al., 2021) and subject to further analysis in this article.
The following model land surface configurations were simulated: and at the lowest tolerance angle and bandwidth, for which the mean difference of altitudes could be derived for each lag. Hence, tolerance = 10° ( = 5°) and bandwidth = 0.5 lag range ( = 50 m) were part of the setup.

| Experiment results
For a planar slope, the directional madogram for = 0° is essentially identical to the longitudinal profile of the terrain. The inclination of the madogram curve equals the slope inclination calculated from the local altitude variability. The madogram for = 90° retains a steady value slightly above 0 (~4.7 m), except for the first two lags. This properties to those derived for planar slopes. It is important to note that madograms are insensitive to curvature, since they are solely based on altitude variability between two points, whilst the shape of the slope between these two points is irrelevant. Hence, directional madograms for concave and convex slopes with identical curvature radii are practically identical (Figure 3s). However, if the curvature changes along the slope and a concave segment grades into a convex one, the shape of the madogram deviates from a straight line and shows higher mean altitude F I G U R E 3 Madograms for simulated DEMs: (a-r) simulated DEMs and their directional (0°, 45°, 90°) and omnidirectional empirical madograms; and (s) 0° madograms of all model land surface configurations differences for the respective lags, particularly in the initial lags ( Figure 3n). This effect is observed in both omnidirectional madograms and all directional ones, except = 90°, but is particularly evident for = 0°. The variability is even higher for initial lags in the = 0° madogram (perpendicular to the frontal zone) if a front-parallel ridge occurs somewhere inside the frontal zone rather than at its inner boundary (Figure 3q). In this case, the highest values typify the lag that coincides with the distance between the axis of the ridge and the frontal zone base, whereas for more distant lags the mean altitude differences gradually diminish (Figure 3r).
In the remaining cases, our idealized terrains comprise additional elements trending perpendicular to the frontal zone or, if a valley network is added (Figure 3k), running both perpendicular and parallel to the front base. They show as secondary ridges (Figure 3g) or main and/or secondary valleys (Figures 3i,k,o). Their presence results in the flattening of madogram curves due to a decrease in altitude variability. If such minor landforms follow a direction broadly perpendicular to the front base, then madograms for = 90° take characteristic shapes. They include an initial rise towards the maximum, which reflects the axial part of this secondary landform (irrespective of whether it is a valley or a ridge), followed by a steady decrease in more distant lags (Figures 3h,j,l,p).
To summarize, an empirical madogram, which is a vector of mean absolute altitude differences, appears an effective measure of variability as a function of distance and direction. Changing values of a madogram indicate relief differentiation, which in turn results from the concurrent action of various geomorphic processes, both exogenous and endogenous. However, relating spatial patterns of processes to the shape of madograms is not straightforward and may even be counterintuitive. For example, a systematic increase in elevation from the base into the frontal zone shows as a relatively steep madogram curve, whereas the presence of secondary valleys and ridges may bring about an initial rise of the curve (local effects), followed by flattening of the curve due to diminishing mean altitude differences. Thus, a higher degree of erosional dissection of an escarpment, which points to high efficacy of surface processes, results in lower values of a madogram. By contrast, if uplift along a boundary fault is poorly counterbalanced by erosion (limited dissection), the respective madogram is steep and attains higher values. Madograms inform about the range of altitude variability, but do not directly reveal hillslope shape or the occurrence of particular landforms. General relief features analyzed in the context of altitude variability are better portrayed by directional madograms, set at the direction of maximum variability. For a frontal zone, = 0° works best. However, an analysis may be supplemented by examining madograms running parallel to the front, which will inform about the degree of dissection. Omnidirectional madograms, as well as those following directions other than 0° and 90°, are less useful due to interpretation ambiguities.

| Zonal and sectoral madograms
The key prerequisite to apply a madogram-based analysis is to delimit the area for which the variability is to be quantified. In our case, it is the extent of the frontal zone. Following Section 2, the frontal zone is equal to the 2,000 m-wide buffer zone that extends from the front base line into the highland terrain. However, in the southernmost part of the area the mountain block is so narrow that sectors of the two opposite fronts overlap inside the mountain range. To avoid situations where a sector includes fragments of both frontal zones, the inner boundaries of the frontal zones were set at the main water divide of the Sowie Mountain and in consequence, the width dimensions of such sectors are less than 2,000 m.
In further analysis two approaches are implemented. In the first, each frontal zone is considered as one spatial unit, characterized by one madogram derived for the entire zone. This allows for the preliminary characterization of the frontal zones and recognition of major similarities and differences. However, the emerging picture is averaged and does not allow one to trace intra-zonal variability (Figure 4). Therefore, in the second, more detailed approach each frontal zone is sliced into sectors of equal length, for which separate madograms are calculated. In addition, some decisions regarding parameters of the method have to be made at this early stage to proceed with calculations. In our real case, we adopted an approach analogous to the one applied in theoretical considerations (Section 5.1). The frontal zones were divided into 20 lags, with 100 m lag range each. The decision should reflect parallel considerations of lag range and DEM resolution. We propose that the lag range should be at least several times bigger than the model resolution, to ensure that a sufficient number of pairs of points exist to estimate the mean difference. For a 5 m × 5 m spatial resolution and 100 m lag size, the proportion is 1:20. Consequently (see Section 5.1), the tolerance angle and the bandwidth were set as 10° ( = 5°) and 0.5 lag range ( = 50 m). The shapes of madograms for the entire NE and SW frontal zones of the Sowie Mountains reflect certain differences between them. The NE frontal zone shows larger altitude differences in general, although they are noted for more distant lags. This suggests that altitude maxima are located closer to the inner boundary of the zone. For the SW zone, higher madogram values in proximal lags indicate higher local variability. In turn, the lower position of the averaged madogram for the entire SW zone in respect to the one derived for the NE zone in the distal lags (>1,100-1,200 m) may signify smaller distances to the high-altitude planar surfaces and ridge axial zones along the former (Figure 4).
To allow for an intra-zonal analysis or a more detailed comparative analysis of different frontal zones, they should be divided into smaller units, further on called "sectors." For each sector, an individual madogram can be derived. The decision about the length of sectors is arbitrary, but specific features of the study area should be considered. The length of 2,500 m adopted here (see Różycka et al., 2021) takes into account that the mean distance between valley mouths of individual drainage basins is less than 1 km. The 2,500 m length increases the probability that each sector would include both valley forms and undissected front faces and hence, individual results would be comparable. Since in a few sectors the inner boundary coincides with the water divide rather than observing the 2,000 m width of the buffer zone, calculations had to be based on a lower number of pairs of points, and were not possible for the most distant lags, especially along the SW frontal zone (e.g., sectors So_SW21, So_SW22) ( Figure 4).
The results reveal both differences between two frontal zones, as well as within them (Figure 4). The SW frontal zone shows considerably higher altitude variability than the NE zone, in the respective lags. In the NE zone

| Isopleth madograms
In principle, the madogram quantifies properties of spatial continuity and variability. However, a madogrambased analysis of mountain frontal zones, which are characterized by considerable difference of horizontal dimensions (length vs. width), implies discretization of the approach and division of the entire zone into sectors. Sector-wise analysis of madogram shapes may be supplemented by an analysis of continuous madogram graphs, the so-called isopleth madogram, which may be likened to the map of altitude variability as a function of distance. The graph shown in Figure 5 utilizes a simple interpolation solution. Having madogram values for n lags and madograms for m sectors, an m × n grid is built, with coordinates of node points defined as follows: where (5) (x, y) = (L × i − 0.5L, 2 × j − ) (6) i = 1, …, m; j = 1, …, n; L = sector length; 2 = lag size The next step consists of spatial interpolation between node points, so that isopleths of mean absolute altitude differences are created. Depending on the purpose of the study, different interpolation methods may be used, such as exact, smoothing, local, or global. In this study we derived isopleth madograms for the frontal zones using the following parameters: 25 m isopleth interval, kriging method, 2,500 m sector, length and 20 lags, observing the location of the inner boundary along a water divide. Exaggeration by five times on the y axis was applied ( Figure 5). This method of visualization allows one to trace changes in mean absolute altitude differences in different lags along the frontal zone base, helping to identify sectors with decreasing or increasing variability. It also makes it possible to compare variabilities in adjacent sectors, as well as in sectors located on the opposite sides of the range. For example, sectors So_SW15 and So_SW18, presented earlier (Figure 4), are not only very different from each other, but also from the surrounding sectors ( Figure 5).
An isopleth madogram can be calculated depending on the purpose of the analysis. If a more generalized picture is recommended, the relevant input parameters would be longer sectors and fewer lags. By contrast, if the examination needs to be more detailed, the number of lags can be increased and the length of the sectors can be set at a lower value. In addition, exact rather than smoothing spatial interpolation would be preferred. The issue of selecting sector parameters is more elaborately commented on below.

| Problem of sector delimitation
For the analysis of frontal zones in the Sowie Mountains sector, dimensions were set at 2,500 m long and 2,000 m wide (unless the water divide was reached) (Jancewicz et al., 2020; Różycka et al., 2021). This was an arbitrary F I G U R E 5 Isopleth madograms for 2,500 m-long sectors along the NE (a) and SW (b) frontal zones of the Sowie Mountains decision, which however respected the morphological features of the frontal zones, such as valley spacing and width of secondary ridges (Section 5.2.1). Thus, the question is valid whether a different choice of method parameters would affect the results. To address this issue, the following additional setups were considered: 1. Sector length variants 1,500, 2,000, and 3,000 m, retaining width of 2,000 m, respecting the presence of water divides.
2. Sector length variants 2,000 and 3,000 m, retaining width of 2,000 m, disregarding the presence of water divides.
3. Sector width variants 1,000, 1,500, and 2,500 m, retaining length of 2,500 m, respecting the position of water divides.
4. Sector length 2,500 m and width 2,000 m, using moving window with 1,500 m shift.
While testing options 1 to 3, sectors are considered separately, whereas in option 4 they partly overlap.
The first test was focused on a comparison of madogram shapes for the variants listed above for the NW and SE extremities of the SW frontal zone, that is So_SW13 and So_SW22 sectors. The decision to select the SW front was dictated by its considerable narrowing in the SE part, which was expected to result in significant differences in madograms calculated with or without taking into account the presence of the water divide as the inner boundary. The division into sectors of predefined length was always executed starting from the NW end of the frontal zone. This means that in each variant, the first sector considered was always of predefined length, according to the variant, whereas the length of the last one could differ from the setup value, depending on how the frontal zone was sliced into full-length sectors. In addition, the water divide in the first sector (So_SW13) is located more than 2,000 m from the frontal base, so that the inner boundary is equal to the extent of the 2,000 m-wide buffer zone. By contrast, the relatively close distance between the water divide and the front base in the SE part (<2,000 m) results in additional differences in the width of the last sector (So_SW22) in each variant that respects the presence of the divide. The 100 m lag range was retained in calculations, so that the number of lags (=length of madogram vector) varied from 10 to 20, depending on the variant (Figure 6).
The first (NW) sector is represented by linear madograms showing a regular, systematic increase of mean absolute altitude differences with distance. Differences between madograms for different variants are small (a few F I G U R E 6 Madograms for (a) the first (NW extremity) and (b) the last (SE extremity) sector of the SW frontal zone of the Sowie Mountains. Variants of sector setups (L = length, W = width): 1, L1,500 m/W2,000 m; 2, L2,000 m/W2,000 m; 3, L2,000 m/W2,000 m unless water divide is reached earlier; 4, L2,500 m/W2,000 m; 5, L2,500 m/W2,000 m unless water divide is reached earlier; 6, L3,000 m/W2,000 m; 7, L2,500 m/W1,000 m; 8, L2,500 m/W1,500 m; 9, L2,500 m/W2,000 m, with overlapping windows meters range), except one, for the 1,500 m sector length. However, the differences in relation to the remaining ones do not exceed +20-25 m for the more distant lags. In the 1,500 m length variant, the first sector terminates at a valley axis and the ridge overlooking the valley on the SE falls into the next sector. This results in an increase of madogram values in the more distal lags (Figure 6a).
In the last (SE) sector, most madograms show similar shapes, but three of them are distinctly different.
Ignorance of the water divide significantly alters madogram shapes, irrespective of sector length (2,000, 2,500 m). This is because the sector becomes wider and includes less elevated terrain beyond the water-dividing ridge, resulting in lowering madogram values in each lag. By contrast, an increase in distal lags is observed in the variant with length set at 3,000 m. Here, the last sector is very short, c. 1,300 m. This decrease in length results in an increase in madogram values, as in the NW sector. For the remaining variants, madogram shapes are similar, although they diverge more than in the first sector (Figure 6b).
The change in length implies rearrangement of the entire spatial patterns of sectors, so it is obvious that madogram shapes will change. The question may be asked whether this impacts the spatial picture of altitude differences, which in turn may imply different interpretations. Therefore, isopleth madograms for different sector dimensions can be examined. Now we turn to the NE frontal zone, since it is less diverse morphologically and hence even minor changes induced by sector pattern change may be easier to detect ( Figure 7).
The general pattern of isopleths in different variants is similar, although following expectations, the longer are the sectors, the more generalized is the emerging picture and the differences between sectors become blurred. This can be exemplified by clear lowering of madogram values in distal lags for sectors located 12-13 km from the NW end. Lower values are clear in the 1,500 and 2,000 m length variants, forming "embayments." They are less evident for the 2,500 m variant, whereas for the 3,000 m length this "embayment" merges with another low-value area, initially detected at 4-5 km from the NE end.
It is worth emphasizing the visual, and hence interpretative, value of the picture obtained from the application of the moving window, with partial overlaps (Figure 7e). In this variant, the madogram map returns well the structure recognized in small sector-length variants, simultaneously smoothing the isopleths, enhancing the visual quality of the diagram. In quantitative analysis, though, one should bear in mind that the same pairs of points were used in calculations for overlapping areas inside sectors.
The above analysis was based on a visual assessment of madogram maps. It was also attempted to quantify the variability of madogram values for different variants in separate lags. Therefore, for each of 20 lags, selected statistics were computed (mean, range, standard deviation) and juxtaposed on graphs (Figure 8).
In line with the idea of a madogram, the values of each statistic tend to increase towards more distal lags.
However, measures of variability (range, std dev) for different variants are more diversified after passing 8-10 lags (i.e., 800-900 to 1,000-1,100 m; Figure 8). This would not be surprising if the 1,500 m sector length variant is typified by the highest variability. In reality, the 2,000 m length variant proves the most sensitive to terrain attributes. Hence, it can be asserted that this variant is better fitted to the real topographic features than the originally selected 2,500 m variant. This approach can inform the choice of an optimal selection of parameters, but this would first require calculations of madograms for all sectors in several length variants, which limits its applicability.
The second key parameter is sector width, whose selection is important because experimental madograms may tend to fluctuate increasingly as the separation distance increases. Thus, Chilès and Delfiner (2012) recommended only calculating variograms for separation distances lower than half the domain size (here, zone or sector width). However, if this proposal is followed, the analysis would be limited to the part of the frontal zone in closest proximity to the front base. Quantification of the frontal zone would then remain incomplete, and specific terrain features relevant for comparative analysis within and between zones may not have been revealed. Therefore, leaving the discussion about the actual width of the frontal zone aside, the influence of sector width on the values of madograms needs emphasis. At a given lag size, reduction in width results in a Sector So_NE8 was selected to illustrate these effects, as here the change of sector width has particularly clear consequences due to morphological characteristics of the terrain (Figure 10). A madogram involving 10 lags only shows a steady rise, reflecting the steepness of the mountain front towards the first row of local elevations and its limited dissection. However, a front-parallel valley occurs some 1.3 km inside the mountain massif, influencing the madogram for 1,500 m width and causing a decrease of its values in distal lags. The same effect is evident for 2,000 m width, whereas it is less pronounced for the widest sector (2,500 m), being counterbalanced by a further increase in altitude towards the inner boundary of the buffer zone. This example shows the sensitivity of madograms to the presence of front-parallel geomorphic elements, both convex and concave. In terms of geomorphological interpretation, such wavy madograms may inform about the degree of dissection of the frontal zone and, depending on the valley depth, the lows in the intermediate/distal lags will be more or less evident.

| Considering morphological differences between frontal zones
Although the frontal zones analyzed so far are likely to have different origins, they nevertheless share many similarities, arising from tectonic control (the presence of a fault line), whether direct or indirect .
Therefore, to provide more informed evaluation of the madogram-based analysis and its potential, we now add another frontal zone (NE front of the Stołowe Mountains), which is morphologically and genetically different, representing a denudational escarpment, with the steepest slope segment near the crest. However, its height F I G U R E 8 Selected statistics of madogram values for the NE frontal zone of the Sowie Mountains, for different variants of sector length range is similar, around 300 m, facilitating interpretation. Following the emerging recommendation to slice the NE frontal zone of the Sowie Mountains into sectors 2,000 m long (Section 5.2.3), identical slicing is applied to the Stołowe Mountains (Figure 11b). Madogram values have been converted to isopleth madograms, according to the procedure presented in Section 5.2.2 (Figures 11c,d).
To identify intra-and inter-frontal similarities and differences, the analytical procedure was expanded to include sector clustering based on vectors of madogram values. Clustering was carried out for the entire population of 29 sectors, from both frontal zones. As a first step, madograms were standardized, since the clustering aimed at recognition of sectors represented by madograms of similar shape, whereas the altitude range was considered less important. Two unsupervised clustering methods were tested, that is k-means and Ward's agglomerative, F I G U R E 9 Isopleth madograms for the NE frontal zone of the Sowie Mountains and variously set width of the zone, from 1,000 to 2,500 m. NW end of the frontal zone is on the left side of the graphs hierarchical clustering. Realizing limitations of the k-means method, which include the necessity of an a priori decision regarding the number of groups and restricted sensitivity to madogram shapes, Ward's method was selected for further work. The method is also known as Ward's minimum variance method, because it aims to minimize the error sum of squares of any two clusters that may be formed at each stage of grouping (Ward, 1963). The resultant dendrogram (tree diagram) allows one to trace how individual madograms join to form clusters ( Figure 12). The main finding is the distinctive association of two frontal zones with two branches of the tree, each forming two groups. One group is formed essentially by sectors from the Sowie Mountains only, the northwesternmost sector of the Stołowe Mountains (St_NE1) included in one of these branches being an exception. The other group consists of sectors from the other frontal zone, save So_NE6 and So_NE7. Division into four main branches was applied at the merging distance equal to 3.5 ( Figure 12a) and reflects fundamental differences in the morphology of  and rather shallow (Migoń & Duszyński, 2018;Pulinowa, 1989). By contrast, in the NW part, sedimentary beds are dipping southwestward and the frontal zone is the front of a cuesta. The longitudinal profile of the front slope does not change much and remains generally concave, but the rim of the front marks the position of the water divide. Beyond this rim the terrain is sloping SW and altitudes decrease. Locally, considerable dissection of the backslope occurs (Jancewicz, Migoń, & Kasprzak, 2019). Thus, madogram signatures for plateau margins and cuestas, both including part of the terrain beyond the plateau/cuesta edge, are effectively different.

| D ISCUSS I ON
This section addresses the usefulness of a madogram-based approach in the specific context of landform analysis, in which information about elevation constitutes the input data. In our approach, empirical madograms were used as the primary dataset for a clustering procedure that aimed at identification of front sectors characterized by similar morphometric properties. In a similar way, parameters of theoretical variogram models were incorporated into landform classifications by Azañón et al. (2004), although in their approach calculations were made within the moving window instead of predefined spatial units, like in this study. Both approaches differ from another one present in the literature, within which variograms are used for quantification of differences between previously distinguished landform classes (Hengl, 2009).
The procedure followed in our analysis helped us to identify various critical steps and these will be discussed below, considering both limitations of the method and opportunities offered by this approach. Although we have focused on one specific type of landform, that is frontal zones of mountains and highlands, more general evaluations and recommendations for future use can also be offered.

| General observations
In common with many other geomorphometric and geostatistical methods, the application of the madogram method requires various arbitrary decisions regarding procedural setup, which, inevitably, influence partial and final results. One group concerns the parameters of the method itself, including direction of madogram, tolerance angle, bandwidth, number and range of spatial lags. Another one includes decisions about the reference area and spatial units, for which madograms are calculated. In each case, these decisions have to be site-and objectspecific, adjusted to the characteristics of relevant landform(s) and reflecting the purpose of the study. In the study reported here, it involved decisions on how to delimit mountain frontal zones, set their width, and slice them into sectors of a certain length. Therefore, the proposal presented in this article cannot be considered universal in terms of numerical setup, although the principles of the analysis apply to mountain fronts and escarpments of any dimension.
The informational value of madograms for both the entire study area and smaller sectors (madogram bundles) can be enhanced through their further processing. In our case, this increased the possibilities to detect intra-frontal similarities/differences and allowed for comparative analysis of different frontal zones. Two methods have been found useful. First, spatial interpolation of individual madogram values results in isopleth madograms ("maps" of madograms), which allow one to interpret discretized madograms while respecting spatial continuity of the analyzed property (in this case, elevation differences). The selection of interpolation method depends on the purpose of the study and the expected level of detail. Second, clustering of vectors of madogram values facilitates spatial analysis of similarities and differences between spatial units. In our case, these were different sectors of frontal zones. Various clustering methods are possible, and the choice depends on the properties to be highlighted, expected to have relevant informational value. In our study, the shape of the madogram was considered important, as it provides information about the morphological diversity of the area, and therefore Ward's agglomerative method was used and is recommended. The k-means clustering approach emphasizes the mean elevation difference, which is less useful in interpretation. In addition, this property can easily be quantified from DEMs, without the need to calculate madograms.

| Advantages and added value of madogram-based approach
The following general advantages of the method are recognized: 1. It is based on one input parameter, which in our case is elevation (more precisely, elevation difference).
Therefore, it is relatively simple and easily applicable.
2. It may be used to analyze parameters which are characterized by spatial autocorrelation and hence, expected differences between sub-populations cannot be tested for statistical significance using classic statistical methods.
Further advantages are more specific to the object of study (mountain frontal zones), but they show that the method is capable of providing added value to landform analysis. First, the method quantifies the morphological complexity of frontal zones using a spatially continuous raster dataset, although a similar analysis based on discreetly arranged point data (e.g., LiDAR point cloud, geodetic/GNSS measurements, data from vectorization of topographic maps) can also be performed. Before, a quantitative approach to the geomorphology of frontal zones was based mostly on simple measures calculated from attributes, which in vector models of geographical data are individual points (such as extreme altitudes within a drainage basin) or lines (such as the length, width, or perimeter of a basin). Second, morphological differences between three frontal zones were highlighted by madograms, despite their similar height. Third, the method proved useful in detecting intra-frontal variability, otherwise rather difficult to show using zonal statistics, which does not allow one to quantify differences as a function of distance and direction. Positive evaluation of the method is supported by the comparison of the results of our study, especially the clustering part, with accumulated geomorphological knowledge about frontal zones subject to analysis (Badura et al., 2007;Duszyński & Migoń, 2015;Krzyszkowski et al., 1995;Krzyszkowski & Olejnik, 1998;Krzyszkowski & Pijet, 1993;Migoń & Duszyński, 2018;Migoń & Zwiernik, 2006;Różycka et al., 2021). In general, the results of clustering are consistent with this knowledge and this statement can be resolved as follows: 1. The method applied to frontal zones which constitute the NE margins of the Sowie Mountains and Stołowe Mountains clearly shows that they are morphologically different, confirming visual, qualitative recognition that straight slope profiles and considerable dissection dominate the NE front of the Sowie Mountains, whereas the NE marginal escarpment of the Stołowe Mountains is generally concave and much less dissected.
2. Two major sub-zones can be distinguished in the Stołowe Mountains. These correspond to the plateau margin (SE part) and cuesta (NW part), respectively.
3. The NE faulted margin of the Sowie Mountains is morphologically fairly homogeneous, which suggests that it represents one segment within a much longer morphostructure of regional extent. 4. Given that these two frontal zones are of different origin (faulting vs. denudation and escarpment retreat) and evolve along different pathways, the method proved sensitive to morphological dissimilarities arising from these geomorphic histories.
Moreover, the method highlighted differences between the opposite fronts of the Sowie Mountains, although madograms for the two fronts were not as strikingly different as in the Sowie/Stołowe Mountains comparison.
Generally, they showed less variability within the fault-generated NE front, where geomorphic evolution was controlled by one major driving factor (uplift) and more variability within the fault-line SW escarpment.

| Limitations
Since the method is based on elevation differences between two points, it ignores terrain characteristics, which occur in between these two points, whether the sloping surface is convex, concave, or more complex.
Consequently, despite diverging slope morphologies, the resultant madograms could be similar and this is confirmed by our experimental madograms for simulated DEMs (Section 5.1; Figure 3). However, real situations are more complex than simple convexity or concavity, and it seems that the sensitivity of the method increases with the complexity of the terrain. We also observed that the shapes of madograms do not directly translate into certain terrain characteristics. For example, distinct wavy shapes are produced (a) if the mountain block is incised by valleys striking parallel to the front, forming a rectangular drainage network with the main valley mouth at the front line, but also (b) if a water divide runs close to the front base, so that drainage beyond it is away from the frontal zone. Thus, the meaning of the first inflection point of the madogram cannot be extracted from its shape alone. However, if such a water divide is set as an inner boundary, then the problem is obviously overcome.
Depending on the morphological complexity of the frontal zone, the choice of sector width influences madogram shape and, consequently, regional interpretation. This is illustrated by the exercise reported in Section 5.2.3.
Finally, although the presented approach is relatively simple, with only one parameter (elevation) involved, the interpretation of madograms may not be immediately intuitive given that mean absolute differences are calculated for pairs of points irrespective of their actual location within the frontal zone, whether close to the baseline or not.

| CON CLUS IONS
The exercise reported in this article was driven by the assumption that madograms can be good representations of elevation variability that occurs within mountain/highland frontal zones (mountain fronts, escarpments). As such, they might allow one to quantify altitude relationships, to subdivide frontal zones on an objective basis, and to develop comparative analysis. Frontal zones are not trivial landforms, but are of global occurrence and are extremely diverse in terms of morphology, origin, evolutionary history, etc. Therefore, the approach proposed here can be applied to analogous features elsewhere, although the parameters of the method have always to be tuned to regional/local conditions such as the height and width of the frontal zone.
We maintain that our working hypothesis was verified positively and madograms, including isopleth madograms, proved useful in better, more objective characterization of frontal zones selected as test areas. They also helped to recognize the internal spatial structure of fronts, highlighting similar and dissimilar sectors. Finally, they clearly confirmed that frontal zones, which limit the Sowie Mountains and Stołowe Mountains on their NE sides, are morphologically different and these differences can be quantified. Despite various limitations, the method can be recommended for other frontal areas, although interpretation cannot be disassociated from real terrain conditions, and for further testing in different geomorphological settings. Avenues to explore include testing the behavior and usefulness of madograms for characterization of landforms of completely different shape and spatial distribution (e.g., glacial cirques, doline fields, volcanic cones), calculating madograms of DEM derivatives (e.g., slope, curvature, topographic wetness index), and testing madograms representing directions different from those experiencing the highest altitude variability. For example, madograms representing a direction parallel to the extension of the frontal zone may help to quantify the degree of erosional dissection.

ACK N OWLED G M ENTS
This article is a contribution to Research Project No. 2015/19/B/ST10/02821, supported by the National Science Centre. LiDAR data for Poland are provided by the Head Office of Land Surveying and Cartography (GUGiK).
Digital elevation data for Czechia are based on the Digital Terrain Model of the Czech Republic of the 5th generation (DMR 5G), processed at an earlier stage of the project under license 96357/16 issued by the Land Survey Office of the Czech Republic (ČÚZK). We are grateful to three reviewers of the journal whose comments helped us to clarify the approach and presentation of results.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data may be shared at reasonable request.