Spatial organization of drumlins

Ice‐sheets flowing over soft sediments produce undulations in the bed, typically of metres in relief, of which drumlins are the most abundant and widely investigated. Consensus regarding their mechanism of formation has yet to be achieved. In this paper we examine the spatial organization of drumlins in order to provide an improved description of the phenomenon and to guide hypotheses of their formation. We review the literature highlighting contradictory findings regarding drumlin spatial organization and then use this to motivate our study based on a large sample (42 488) of drumlins from Canada, Britain and Norway. Are there typical arrangements in drumlin positioning and are they organized in a regular spatial manner (patterned) or are they distributed randomly? We recognize that drumlin fields are inherently patchy and therefore apply inhomogeneous spatial statistics in order to study their distribution. This shows that whilst drumlins are occasionally randomly placed, their main state is non‐ random. They exhibit a strong and statistically significant signal of regularity across lengths scales of 100 to 1200 m. We conclude that patterning is a near ubiquitous property of drumlins. This finding of regularity demonstrates spatial self‐organization in the bedforming process with drumlins as an emergent manifestation of sub‐glacial sediment mobility. Kilometre‐scale interactions between drumlins must occur as they evolve, or interactions may arise as a consequence of growth or migration. Hypotheses or models are required that can explain the regular spacing of drumlins. We highlight three suggestions for such self‐organization: instability in the coupling of ice flow–sediment flux–bed shape; local feedback between sediment mobility and relief; and coarsening by growth or migration. © 2017 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd.


Introduction
Somewhat counter-intuitively, ice-sheets abhor flat beds when flowing over soft sedimentary substrates. They produce undulated surfaces, metres to a few tens of metres in relief and with length-scales of hundreds of metres to kilometres, that have been categorized into a number of easily recognized states with well-known geomorphological labels; those oriented transverse to flow (ribbed or Rogen moraine; Dunlop and Clark, 2006, to mega-scale transverse ribs; Greenwood and Kleman, 2010) or parallel to flow and varying in length, from drumlins (Clark et al., 2009;Eyles et al., 2016) to mega-scale glacial lineations (Clark, 1993;Spagnolo et al., 2014). They can be thought of as a family of landforms called sub-glacial bedforms (Aario, 1977;Rose, 1987;Ely et al., 2016) with an implied genetically related formation process, or each type (drumlin, ribbed moraine, etc.) or sub-type (e.g. ribbed moraine, Rogen moraine, Blattnick moraine, minor ribbed moraine, etc.; Hättestrand, 1997) can be regarded as fundamentally different landforms that were formed by different processes.
These sub-glacial landforms are estimated to cover around half the area previously occupied by palaeo-ice sheets (Clark et al., 2009), with drumlins by far the most abundant type. Drumlins have recently been discovered beneath the Antarctic Ice Sheet (King et al., 2009) and emerging from an Icelandic ice cap (Johnson et al., 2010;Jónsson et al., 2014).
Consensus regarding the formation of drumlins has yet to emerge in spite of two centuries of investigation (Hall, 1815). Rather, there are many hypotheses, ideas and conceptual models, seeking to explain their internal sedimentary structure and the shape, size and pattern of the landforms (overviews by Menzies, 1979;Patterson and Hooke, 1995;Shaw, 2002;Clark, 2010;Stokes et al., 2011;Eyles et al., 2016). Progress has matured with some model-based development of the physics of ice-sediment-water interactions at the base of ice sheets (Hindmarsh, 1998;Fowler, 2000;Iverson, 2000;Hooke and Medford, 2013), but these have yet to be developed to yield many predictions for testing, although see Dunlop et al. (2008), Fowler et al. (2013) and McCracken et al. (2016) for such model-data comparisons. As with many seemingly intractable scientific problems, researchers from various disciplines (geomorphologists, glaciologists, sedimentologists, and theoreticians) have approached the task from different angles, leading to useful observations, theory and insights and some disagreement regarding how these connect. The so-called 'drumlin problem' (Menzies and Rose, 1987) remains difficult, and we do not solve it here. Our approach in this paper is to make progress by investigating the spatial organization of drumlin positioning within their flow-sets; are there typical arrangements and are drumlins organized in a regular spatial manner or are they distributed randomly? The former is indicative of an overall organization within drumlin flow-sets, suggesting they are a patterned phenomenon, of which nature provides numerous examples (Ball, 1999). Conversely, the latter is consistent with a 'sitespecific' view whereby each drumlin forms in isolation from its neighbours. Data and findings answering these questions should guide the development of theory or could be used to test numerical process-models that are both starting to make quantitative predictions and to yield three-dimensional model simulation movies of bedform organization, evolution and migration (Chapwanya et al., 2011;Barchyn et al., 2016). Do we need to explain the development of the individual bump (the drumlin) or the wavelength (spacing) and relief of the undulating surface; the drumlin field? Does the overall bumpiness arise in a single phase of development (triggered simultaneously, like measles) or arise from incremental additions of new bumps that 'know of' the existing ones? Rippled surfaces created on a beach in one tide would be an example of the former and Barchan dunes growing and migrating as individuals across a desert are the latter. Unlike these non-glacial examples, we have yet to observe a flat surface beneath an ice sheet grow into drumlins, thus allowing many questions to remain, and highlighting why investigating spatial properties of the final forms might be useful. We review the literature to highlight observations and analyses regarding drumlin spatial organization and then use these to motivate our study based on a large sample of drumlins.

Existing observations and analyses on spatial organization
The literature on drumlins abounds with comments on the relationship of a drumlin to its neighbours (see reviews by Menzies, 1979;Patterson and Hooke, 1995;Clark, 2010;Eyles et al., 2016). An en echelon arrangement (Figure 1a) has often been mentioned (e.g. Armstrong and Tipper, 1948;Clapperton, 1989) although a survey of such findings (Patterson and Hooke, 1995) notes that it is not common and no field of drumlins has been found where such positioning dominates. Much more Figure 1. Types of drumlin spatial organization. (a) En echelon arrangement with the next lateral drumlin positioned slightly downstream (like geese in flight, with each bird positioned to avoid turbulence from an upwind flock member). (b) Transverse banding in drumlin densities. (c) Packed versus widely-spaced distribution. (d) Random, clustered or regularly spaced positioning within a flow-set, the latter having a preferred spacing expressed as a single narrow peak in the frequency distribution. In spatial statistics an homogeneous distribution of features is taken to be that they are uniformly distributed across an area (i.e. no large gaps). As argued in the text and apparent here and in Figure 2. we show that drumlins are more appropriately considered as being inhomogeneously distributed in that they often occur in patches with large gaps between them. C. D. CLARK ET AL.
common is the observation of transverse banding of drumlins ( Figure 1b) such that in the downstream direction there are alternating bands of high and low densities. Hill (1973) for example quantified such variations for c. 3900 drumlins in Ireland. Elsewhere, the size and shape of drumlins within each band have been reported to be similar (Shaw and Kvill, 1984).
The spacing between neighbouring drumlins is a useful measure that quantifies an aspect of their overall spatial organization; a preference for a specific spacing implying regularity ( Figure 1d). Typically, drumlins are usually found to be widely-spaced with large gaps between them, but in places they have also been reported as tightly packed (i.e. centre-tocentre spacing controlled by drumlin size; Rose and Letzer, 1977) as depicted in Figure 1c. Reed et al. (1962) reported 204 measurements of drumlin-to-drumlin spacing in three drumlin fields in the United States, finding their spacing to be rather varied (and multimodal) but with some tendency towards periodicity (i.e. a single preferred spacing). With a larger sample size of drumlins (733) in Ireland, Vernon (1966) found no preferred spacing and that they varied throughout the fields investigated. Drumlin spacing in Poland (Baranowski, 1969) and Canada (Clark and Wilson, 1994) was again found to be varied, but in Ireland and with a larger sample size of 3900, Hill (1973) reported a preferred spacing of around 300 m. Consistent with the earlier mentioned papers, and according to the review by Menzies (1979), for drumlins sampled in the United States, Canada, Ireland and Poland the frequency distributions of spacing were sometimes found unimodal and normally-distributed, or multimodal and with wide variations in average spacing leading to the conclusion of no overall preferred spacing for drumlins. It is possible however, that if any actual preferred spacing of drumlins exists that some combination of varied mapping approaches, or non-representative sample sizes might be masking the answer. Or perhaps, genuine differences in preferred spacing exist from place to place whose reasons (substrate geology or glaciological controls?) could be usefully sought.
The arrangement of drumlins within flow-sets may provide clues for theories of drumlin formation; inferences of 'process from pattern'. Are they positioned randomly or as regular arrays or as a series of clusters ( Figure 1d)? Random positioning appears to have been the most commonly presumed arrangement, in part promoted by the nearest neighbour statistical analyses of Smalley and Unwin (1968), and which may have steered theory development towards drumlins developing from (randomly distributed) sediment inhomogeneities (e.g. patches of stiff till or gravel) or from bedrock bumps. Francek (1991) using a much larger sample (> 4000) than the 1968 investigation also reported drumlins to be distributed mostly randomly when analysed by statistical analyses. Contrary to these statistically-based findings, many have remarked on the apparent regularity expressed by drumlins (e.g. Werth, 1909, cited in Baranowski, 1969Baranowski, 1977;Carl, 1978;Smalley and Warburton, 1994;Fowler, 2000;Clark, 2010) thereby arguing for some form of regularizing process inherent in drumlin formation. Examples of drumlins that appear to exist with a repetitive, systematic arrangement are not difficult to find (e.g. Figure 2) and from such observations many might simply infer that drumlins are positioned regularly. However, as Table I highlights published investigations involving statistical analyses of drumlin distribution shows that there is no consensus regarding whether drumlins are placed regularly, randomly, or clustered. Different areas and statistical techniques have produced all three results.

Data set
In order to address the uncertainties in drumlin spatial organization highlighted earlier we used a geographic information system (GIS) database comprising 42 488 drumlins spread across 72 different drumlin fields, three palaeo-ice sheets and countries, and across a wide diversity of geological substrates. Most of these (36 222) were mapped across Britain from a high resolution (5 m) digital elevation model (Hughes et al., 2010). Additional mapping was conducted from Landsat ETM+ imagery (15 m resolution) in Alta, Norway (1483) and Ungava Bay, Canada (5903). An advantage of this database over those used in previous studies is that it has a larger sample size, covers three palaeo-ice sheets, and uses drumlins from highly geologically diverse locations such that drumlins in the sample cover a range of geological substrates and with widely varying drift thicknesses, although we do not analyse these here [see Greenwood and Clark (2010) for such an analysis in Ireland].
The dataset (Figure 3) was divided into flow-sets (i.e. fields interpreted as forming under single phases of ice flow) based upon their morphology, parallel conformity, spacing and orientation (Clark, 1999).

Transverse banding
We used our database to address Hill's (1973) query as to whether the transverse banding he reported for some Irish drumlins is a ubiquitous characteristic. A simple procedure manually linking (by drawing lines in the GIS) laterally adjacent drumlins was adopted. The rules being that an adjacent drumlin belongs to a transverse band if it is positioned laterally, is of similar size and orientation and is not further away than around four times the local drumlin width. This latter rule avoids linking drumlins that are too distant. A quantitative rule-based GIS algorithm could no doubt be devised to accomplish this task with less potential for bias, but we judge our manual method appropriate to a reconnaissance level exploration of the extent to which transverse banding exists. Of arguably greater concern than any such bias is that randomly positioned objects will exhibit alignments by chance and we seek to know if there is any spatial organization above this  level. Figure 4 illustrates this issue by comparing a real drumlin field to a simulated one whereby ellipses of equivalent size and orientation have been randomly positioned across the same domain.

Nearest neighbour distances
To assess any tendency for drumlins to exist at a preferred spacing (see Figure 1d) nearest neighbour distances were calculated as Euclidean distance (within software ArcGIS) between the centroids defined for each drumlin. Centroids were calculated automatically in ArcGIS as the geometric centre of the drumlins, i.e. the arithmetic mean position of all the points in the drumlin's polygon shape. Across-flow spacing of drumlins was also calculated using the mean ice-flow direction (from 10 neighbouring drumlins) and measuring the shortest orthogonal spacing between the centroids of adjacent drumlins. Measurements were used to derive frequency histograms of drumlin spacing.

Spatial statistical methods for assessing regularity
In order to investigate the nature of spatial organization in the distribution of drumlins (random, clustered, regular) statistical methods are required, but with spatial phenomena this is not trivial (Kulldorff, 2006). Tests for randomness, such as the Clark-Evans statistic used in the first published investigation of drumlin distribution (Smalley and Unwin, 1968) and the more advanced Ripley's K and L functions (Ripley, 1988), compare measured distances between drumlins to expectations of complete spatial randomness (CSR) producing a metric indicative of dispersion (regularity), randomness or clustering. The K-function [K(r)] is: whereby λ is the mean density of points, e ij is the Euclidean distance between the ith and jth points in the dataset, r is the radius of the search area and n is the number of points within the dataset. With CSR, the expected value is: For data analysis the related L-function is more appropriate to use as it is more stable and is defined as: These functions have been applied in numerous fields such as ecology and cell biology (Dixon, 2012). Unfortunately, these homogenous tests are sensitive to the area being tested, and do not consider changes in point density across a domain. This is an acute problem for drumlins because we have no a priori knowledge of the precise spatial domain. Should this be defined by a convenient box or by the limits of the furthest drumlins? For the latter there is rarely a clear edge, just a decrease in density. It is noteworthy that this issue of sensitivity to the chosen domain led to the later retraction (Unwin, 1997) of the original and widely assimilated conclusion of drumlins having a random distribution (Smalley and Unwin, 1968).
A further problem for any analyses is that areas exist within drumlin fields where drumlins are absent (e.g. Figure 2b). Should these gaps be regarded as 'no observations' and thus be included or as 'no data' because drumlinization could not occur here and so should not be part of the domain? Such absences often exist due to a lack of sediment cover, although preservation and visibility issues caused by lakes, bogs and urban areas also occur.
We undertook extensive experimentation using Ripley's K [K(r)] and L [L(r)] functions for analysing spatial point patterns. We used these functions (in R statistical software) to compare measured distances between drumlins to expectations of CSR, producing a metric indicative of dispersion (regularity), randomness or clustering. To further this approach it is common to use simulation envelopes to compare the observed result Figure 4. Comparison of transverse banding evident in a real drumlin field (a; Leitrim Ireland) with a simulated one (b) using ellipses of the same size and orientation but positioned randomly. A few red lines are drawn in to illustrate banding, many more are visible. Note that in (b) some transverse banding and en echelon positioning is evidentpurely by chancebut is less prevalent and it is much rarer to find similarly-sized drumlins comprising a band. [Colour figure can be viewed at wileyonlinelibrary.com] SPATIAL ORGANISATION OF DRUMLINS to Monte-Carlo simulations of CSR derived from a Poisson (or other) process (Baddeley et al., 2014). Defining the bounding box is an important consideration, as it forms the spatial domain over which simulations of CSR are implemented and . The problem of patchiness (density variations) in drumlin fields for analysis of spatial organization. The preference for drumlins to occur in patches of higher density amongst more sparsely distributed forms yields a clustering in the spatial arrangement, and with some gaps often arising from lack of sedimentary cover. (a) Typical example from northern England (flow-set 66 from Hughes et al., 2014) with large gaps and patches of varying density. (b) To illustrate a point we show several patches of graph paper (each with obvious regularity). Ripley's L function applied to both the drumlins (a) and the graph paper (b) yield clustering as the signal (c), and yet for the graph paper (b) we know that regularity actually exists. The inhomogeneous pairwise correlation function however, takes account of the inherent clustering and shows (d) that for both the drumlins and the graph paper there is a strong preference (attraction) for points to be spaced a certain distance apart. Clustering has not masked the regularity that exists. [Colour figure can be viewed at wileyonlinelibrary.com] C. D. CLARK ET AL. thus influences the null hypothesis to be tested. We used a Ripley-Rasson bounding box, which is a convex hull enlarged by a factor dependent upon the number of points within the study area (Ripley and Rasson, 1977). Our analyses across the drumlin database led to confusing results, reported in Ely (2015). By varying the domain size over which the function is performed, we found it possible to yield different answers as to the spatial organization of drumlins, with clustering as the most common conclusion (see Maclachlan and Eyles, 2013). Typically, with the domain encompassing the whole flow-set, the results tend towards clustering, but if a reduced domain is positioned over a prominent patch of higher density drumlins then regularity is often revealed such as in Figure 5, but the sample size for such an assessment is low.
It is this seemingly contradictory finding that led the realization that drumlin fields are actually more appropriately characterized as patches of high density within a background of lower density and with many large gaps (Figure 6a). Such inhomogeneity turns out to be critical with regard to how the distributions are analysed. By way of explanation, if Ripley's functions are performed on points on a regularly spaced grid (with obvious regularity) but with areas rubbed out ( Figure 6b) the dominant result will be clustered, even though we can clearly see that there is strong regularity within the clusters. We suspect that the widely varied results reported in the literature (see Table I) is a consequence of drumlin fields typically being patchy and the methods used being overly sensitive to this patchiness. This issue of inhomogeneity has been widely recognized in ecology, where variables such as soil availability and nutrition introduce inherent clustering of tree positioning (Perry et al., 2006) that is independent of any biophysical spacing controls between individual trees. If we are not to be misled, a technique is therefore required which accounts for the inherent patchiness within a drumlin flow-set; a method of spatial analysis for inhomogeneous point patterns rather than homogenous patterns is required. This is important because arbitrarily defining clusters can be problematic (Couteron et al., 2003) and if one experiments using the Ripley method within individual clusters (e.g. Figure 5) then the sample sizes typically become too low to gain meaningful results. Furthermore, Ripley's L function is a cumulative function, accounting for all previous measurements in subsequent iterations of the search radius. Thus, inferences at specific scales are difficult to make. Therefore, a different method is required for testing whether spatial positioning of drumlins is regular, random or clustered because they have inhomogeneous spatial distributions.
Inhomogeneous pairwise correlation function (g inhom ) Ripley's K-function and the L-function are cumulative in the sense that K(r) and L(r) count all points (drumlins) with distance at most r to a randomly selected point of the pattern. If one wants to infer at specific distances r, one should use the pair correlation g(r) function instead (e.g. Perry et al., 2006). It is derived from Ripley's K-function in the following manner: The pair correlation function tells us how likely it is to encounter another point of the pattern at a distance r from a 'typical' point, as compared to the corresponding likelihood under complete spatial randomness. Thus, g(r) = 1 means that, seen from a typical point of the pattern, the world at distance r looks completely random (i.e. like a Poisson point process), implying no interactions between points. If, however, g(r) < 1, then points at distance r are less likely to occur than they would under CSR; this is termed 'repulsion' between points at distance r, or that points have been inhibited from occurring at this scale. Conversely, g(r) > 1 means 'attraction' between points at distance r; they are more likely to occur than under CSR. Regular point patterns are characterized by repulsion at small r, followed by a peak of attraction at the typical neighbour distance, while clustered point patterns show attraction already for small values of r (Illian et al., 2008, p. 240).
The original pair correlation function is as sensitive to patchiness and choice of domain boundaries as is the K-function. However, we will see that this problem can be overcome by using the inhomogeneous pair correlation function, g inhom (Baddeley et al., 2000; Figure 6d). The g inhom function is thus appropriate for deciphering the spatial distribution of drumlins. In order to account for inhomogeneities in a dataset, the function takes into account the intensity (density) of the observed points. This is advantageous over arbitrarily defining smaller regions of perceived homogeneity and then performing a homogenous test (Couteron et al., 2003). Instead, g inhom weights the point counts according to the local intensity of points within the observed pattern. It is equivalent to comparing the observed pattern with inhomogeneous CSR of same (inhomogeneous) intensity.
In order to assess the spatial arrangement of drumlins over a long range, we computed g inhom (r) for a range of r-values. We then simulated CSR through Monte-Carlo simulations of a Poisson process producing the same number of observations and point density, to derive a significance envelope (Baddeley et al., 2014). Deviation from this envelope indicates that the observed patterns are different from CSR. Recently, a technique for assigning significance to this deviation has been developed (Myllymäki et al., 2016), and is implemented here.
The earlier approach was implemented in the open-source statistical software R. We employed a Ripley-Rasson bounding box (Ripley and Rasson, 1977) and the same edge corrections as used for the Ripley L tests performed earlier. The libraries 'spatstat' (Baddeley et al., 2015) and 'spptest' (Myllymäki et al., 2016) were used to perform the tests and generate Monte-Carlo significance envelopes. A typical example output from the procedure and labelling is shown in Figure 7.
To account for the clustering in inhomogeneous point patterns an estimation of the point intensity (density) is required. In spatstat, this is done by a kernel density estimator yielding a pixel image. The intensity value in a given pixel is obtained as the weighted number of points within a given search radius. This distance is also called the 'kernel band width'. As a parameter of the estimation procedure, the search radius strongly influences the resulting intensity estimate (O'Sullivan and Unwin, 2010, p. 70), and thus the estimate of g inhom . Therefore, sensitivity analysis was conducted upon this parameter in order to determine an appropriate value for analysis. Flow-set 29 in southern Scotland (Hughes et al., 2010(Hughes et al., , 2014 was chosen for sensitivity analysis due to the high number of drumlins (n = 1473). Figure 8 reports the sensitivity of g inhom (r) to perturbations of the default value (as a percentage) of the search radius. A selection of results (80%, 40%, 30%, 5%, Figure 8) are presented but were explored at 10% increments between 100 and 10% and then also at 5%. For the default value (100%) the intensity of points across the flow-set is smooth increasing away from the edges of the flow-set due to edge corrections and appears similar to the 80% case shown in Figure 8. Small-scale dispersion is apparent, followed by a strong peak in attraction and final decay to CSR in a manner similar to Figure 7. This is the typical signal across most values in the sensitivity analysis (100-30%). With further reduction in the search radius (lower percentage) clusters within the data become apparent (column two of Figure 8). Reduction of the search radii below 30% of the default value increases the size of the CSR envelope at smaller values (< 500 m). By 5% of the default value the peak no longer occurs above CSR, as the intensity estimate forces simulations of CSR to occur close to the drumlin pattern, almost mirroring it. However, at 5% the intensity map has identified almost as many clusters as there are drumlins such that often only a single drumlin or pair of drumlins is within a single intensity cluster and this is clearly below the length-scale at which we are investigating spatial organization. For all other flow-sets, a similar pattern of an initial smoothing across the flow-set with an increasing clustering at 30% is apparent. Thus, 30% of the default bandwidth parameter was chosen for all flow-sets.

Transverse banding
Banding is found to be widespread, with 48% of the drumlins analysed for banding (the British examples, n = 36 222) found to belong to a transverse band (Figure 9). This extends the findings of Hill (1973) showing that such banding is a common feature of the drumlins we have examined and shows that there is a strong element of spatial organization in drumlin positioning. The mean downstream spacing of bands (n = 1885) is 2026 m and the frequency distribution approximates that of a lognormal distribution with a mode of around 800 m. Although we do not demonstrate it here, our visual examination of drumlins elsewhere in the world (using Google Earth) suggests that transverse banding is commonplace. We note that our definition of what constitutes a band is subjective as was our visual method of defining them, so the results quantifying the spacing should be taken as indicative only, pending further investigation. Later we report quantitative and more robust investigations on spatial organization.

Nearest neighbour spacing
The nearest neighbour spacing of drumlins ( Figure 10; British drumlins) reveals a strongly preferred concentrationa unimodal peakaround 200 to 450 m. This provides a strong hint that regularity exists, but is not proof because relative spatial positions are not analysed and the preferred spacing could be a consequence of drumlin dimensions. A second analysis, designed to assess for the effect of drumlin size was performed; is the preferred spacing simply a consequence of drumlins being packed into a limited area (i.e. Figure 1c)? The analogy here is how ball-size in a container full of balls controls their spacing (i.e. spatial inhibition). Nearest neighbour analysis conducted Figure 9. Of all British drumlins in our database 48% were found to belong to transverse bands suggesting that this style of spatial organization is a common feature of drumlin flow-sets. Drumlins classified by our method as belonging to a band are depicted as black.
SPATIAL ORGANISATION OF DRUMLINS orthogonal to ice flow direction revealed drumlin spacing (centroid to centroid; mean of 386 m) to exceed drumlin width (mean of 226 m). This rules out drumlin dimensions being the primary control on spacing. Visual inspection of the arrangement of typical drumlin fields (e.g. Figure 2) leads to the same conclusion in that they are usually found to be widelyspaced rather than packed together and touching. Spatial statistical methods for assessing degree of regularity When applied to a typical flow-set the plot of g inhom (r) exhibits a clear signal of attraction, well above that which occurs due to random fluctuations and with a prominent peak at 330 m ( Figure 11). This demonstrates that there is a strong and statistically significant attraction of points within this nearest neighbour distance, between the length scales of 160 to 1250 m. This combination of short range repulsion followed by a strong peak in attraction is typical of regular patterns. Moreover, of the 42 488 drumlins analysed, 74% of the flow-sets, including those in Canada and Norway, displayed regularity as summarized in Figure 12, representing 93% of the drumlins. By accounting for inherent patchiness and gaps in drumlin fields, this is the first robust demonstration that regularity exists and that it is the predominant signal of drumlin spatial organization. Furthermore, those drumlins that were found to be not regularly distributed tended to occur in flow-sets that contained small numbers of drumlins (n < 300). The attraction peak of g inhom (r) occurs at a mean value of~600 m, with the curve often continuing to 2000 m, indicating that drumlins prefer to be positioned a nonrandom distance apart across this length scale. Such regularity demonstrates drumlins to be a patterned phenomenon.

Discussion
The finding of widespread banding and a uni-modal peak in the nearest neighbour analyses suggest that drumlins are not randomly placed, but these simple approaches are not a robust demonstration of this. However, by correctly accounting for the patchiness and gaps in drumlin fields by using an appropriate statistical method for analysing such inhomogeneous phenomena we have robustly demonstrated that whilst drumlins can be randomly distributed (7% of drumlins analysed; 26% of flowsets), their most common state (93% of drumlins analysed; 74% of flow-sets) is non-random. They exhibit a strong and statistically significant signal of regularity across lengths scales of 100 to 1200 m indicating that inter-drumlin interactions occur over this scale range. Regularity also occurs between 1200 to 10 000 m, but some of this might arise from regularity in the spacing of clusters in addition to longer range inter-drumlin interactions. In our database the finding of regularity holds across different ice sheets, countries, and across a variety of substrate geologies. Because of this, and the size of the database, we suggest regularity to be a near ubiquitous property of drumlins. Given that it took ourselves and many others so long to realize that the gaps in drumlin fields (due for example to lack of sediment and bedrock outcropping) needed accounting for we now wish we had paid more attention to the words of Menzies in his review (1979, p. 338); 'It must be further noted that in most studies of drumlin spacing, little or no account is taken of other influencing structures such as roches moutonees or bedrock protrusions [i.e. gaps]. It is the writer's opinion that such influencing features are equally important in such morphometric studies.' Although the spatial statistical analysis used in this paper demonstrates strong regularity the method says nothing about Figure 10. The spacing of drumlins. Frequency distribution of nearest neighbour distances (avoiding double-counting) between drumlin summits in the UK database (n = 36 222). The uni-modal peak centred on 325 m is interpreted as indicating a clearly defined preference for the spacing of drumlins, suggestivebut not provingthe existence of patterned spatial organization. Figure 11. The inhomogeneous pairwise correlation function (PCF) for flow-set 29 in the Southern Uplands of Scotland comprising 1473 drumlins. The function g inhom (r) is plotted against r, with increasing radii of analysis (in metres) away from each drumlin. The black line plots the function for the drumlins which should be compared against the red line and grey envelope representing complete spatial randomness (CSR). The peak at 330 m and extending between 160 and 1250 m indicates a statistically significant demonstration of regularity (i.e. above grey envelope) in drumlin positioning across this length scale; 'attraction' of points to this distance has occurred, typified by a regular spatial arrangement. Values plotting within the grey envelope are not significantly different from a random spatial arrangement. [Colour figure can be viewed at wileyonlinelibrary.com] C. D. CLARK ET AL.
the spatial structure of the inherent organization (i.e. how the pattern is organized). This is because the analysis works by assessing the departure from randomness using circlesincreasingly long radiifrom each drumlin. So we know that regularity exists but not the form it takes. Connecting this finding with the observations made earlier (albeit less quantitative or secure) about the transverse banding, we interpret that these are likely a large influence on the measured regularity, but also that lateral spacing contributes. The idealized form of spatial organization of drumlins can thus be depicted as in Figure 13. We interpret that it is the higher density zones (drumlin patches) that misled some of the earlier statistical analyses of drumlin fields into concluding that clustering is the dominant signal. Well yes, clustering exists in drumlin fields but likely arises from antecedent variations in sediment thickness rather than from the drumlin-forming process. Once this is accounted for the drumlin-to-drumlin positioning is mostly regular.
The regularity of drumlin positioning demonstrates that they are a patterned phenomenon; a pattern being defined as a discernible regularity in which elements are repeated at specific length-scales. This finding is important as it should guide the quest for understanding the processes of drumlin formation and for evaluation of existing hypotheses. It says that a 'sitespecific' landform view that considers each drumlin as a separate individual formed at a specific site and essentially independent of others, is not appropriate for most drumlin fields (but see later for rock-cored drumlins, etc.). Arguably it might be better to consider the drumlin array or field as the primary phenomenon and to regard them as bedforms sensu stricto (cf. Menzies and Rose, 1987), characterized simply by amplitude, relief and orientation. In this view the definition of drumlin boundaries, as in most studies including this are somewhat arbitrary and unnecessary. Ripples on a beach for example are not mapped and analysed individually but treated as a patch of patterned and bumpy bedforms with a length-scale, Figure 12. Demonstration of regular patterning in drumlin arrangement. All 72 flow-sets analysed (with n > 30) displayed on a log-log plot illustrating that the main signal is of regular patterning. The majority of flow-sets (blue lines; 74% of total) are found to have a statistically significant demonstration of regularity, representing 93% of the drumlin sample (n = 42 488). The green lines are for those flow-sets that did not exhibit regularity (26% of the sample), approximating either complete spatial randomness or clustering. The black line is for flow-set 29, from Figure 11. In the middle panel, black dots record the position of the first prominent peak of each curve (main signal of regularity at this scale) and blue dots mark where the signals decline into complete spatial randomness (CSR). These provide the scale ranges for our conclusions in the panel beneath. [Colour figure can be viewed at wileyonlinelibrary.com] SPATIAL ORGANISATION OF DRUMLINS relief and orientation. Rather than speaking of drumlins perhaps it would be better to call the phenomenon drumlinized terrain.
When patterns are found in nature they are usually organized by a process (Ball, 1999). Some of the most striking natural patterns are geomorphological (Werner, 1999) or biological (Koch and Meinhardt, 1994). Exemplars include dunes (Kocurek and Ewing, 2005), freeze-thaw patterned ground (Kessler and Werner, 2003) river meander networks (Hallet, 1990), the fronds on a fern, and leaf positioning on a tree. It seems that all patterns in nature must come into being by interactions between elements as they evolve (grow or shrink). They are not 'printed' as a pattern but evolve into one. Elements of the pattern somehow 'know of' other elements to control their spacing, similar to geese in flight where the flock self-organizes into an en echelon pattern by each individual keeping a set distance and angle from just one neighbour. We suggest that such self-organization is intrinsic to the drumlinizing process with interactions occurring between bumps at a kilometre scale; they must interact with each other whilst evolving. The spatial organization of drumlins can therefore be regarded as an emergent property of the processes at work.
In Figure 14A we suggest that patterning of drumlins might arise from local self-organization, analogous to how regularly-spaced cusps form on a beach (Werner and Fink, 1993). It is thought that a slight depression encourages flow acceleration of water and increased erosion to create embayments with deposition at the periphery and with regularity arising from such competition between adjacent embayments and the smoothing imparted at the scale of the beach. The point being that there is a positive feedback between sediment mobility and relief. The drumlin model of differential erosion by regelation-infiltration with variations mediated by effective pressure (Iverson, 2000) seems to be an example of this type.
In other pattern-forming geomorphological systems, such as ripples and dunes, progress has been made with explanations rooted in non-linear dynamics and complex systems (Kocurek et al., 2010;Murray et al., 2014). Spontaneous development of a pattern from a system of mobile material with coupling between the fluid flow, sediment flux, and the bed shape, is usually taken to arise from instability in the system that amplifies tiny disturbances. Infinitesimally small perturbations at the bed grow at a preferred wavelength (spacing) and instability analysis predicts the wavelength at which bumps should grow, generating repeated culminations with regularity at this preferred scale ( Figure 14C). One possibility then is that the scale of regularity we have identified in our database is a Figure 13. Drumlins are not randomly distributed; rather they are a patterned phenomenon that we idealize in this cartoon. They mostly occur in patches of higher density with internal regularity and banding amongst more sparsely distributed drumlins and with numerous gaps evident in the overall drumlin field. The gaps might be due to thin sediment cover (or bedrock exposed in places) or inappropriate sediment properties. Figure 14. Suggestions of self-organization that might individually, or in some combination, explain patterning of drumlinized terrain. In (A) adjacent bumps interact and from this emerges a field-scale pattern with characteristic spacing, whereas as in (B), a perhaps random distribution of bumps self-organize into a pattern by collisions arising from migration and growth. In (C) instability occurs with the spacing controlled by the wavelength of the maximum growth rate. In this case simple stabilization occurs at this wavelength, but it could be that migration and coarsening also occurs (i.e. as per B) after the initial perturbation. [Colour figure can be viewed at wileyonlinelibrary.com] C. D. CLARK ET AL.
consequence of the wavelength of preferred growth in an unstable system. Examples of drumlin models of this type include the Hindmarsh-Fowler instability in the coupled flow of ice and sediment (Hindmarsh, 1998;Fowler, 2000), a thermomechanical instability in a partially frozen till substrate (Hooke and Medford, 2013) or as intrinsically required in the megaflood hypothesis (Shaw, 2002). The latest version of the Hindmarsh-Fowler theory (Fowler and Chapwanya, 2014) includes sub-glacial water flow within an upper layer of watersaturated sediments (i.e. more realistic than a thin water film across which stress cannot be transmitted) and this addition results in sediment erosion by water flow breaking up sub-glacial ribs to yield laterally-separated bumps resembling drumlins. Outputs from the numerical model estimate wavelengths consistent with measured drumlin dimensions (Clark et al., 2009) and the regularity that we report here. But there is some naivety in this view, in that it is not known if there is a mechanism that stabilizes the bumps at this wavelength, rather than allowing them to develop further. There are reasons for expecting further evolution to occur, and the prominent transverse banding found in drumlin fields might be an example of a coarsening process. From analysis of other geomorphological systems it is considered that finite-amplitude pattern formation could stabilize at the preferred wavelengths predicted by the instability analysis or, alternatively, as the bumps migrate they could collide and coalesce ( Figure 14B) with such interactions (coarsening) resulting in regularity at an entirely different scale (cf. Kocurek et al., 2010;Murray et al., 2014).
Field measurements of recently formed drumlins in Iceland (McCracken et al., 2016) appear to be inconsistent with the Hindmarsh-Fowler instability theory (Fowler and Chapwanya, 2014). If the identified instability mechanisms are not the underlying cause of drumlins it could be that individual drumlins initiate at randomly occurring locations, say triggered by stiff patches of sediment (although see Menzies et al., 2016;Hart, 1997), and then grow in dimensions over time such that they eventually meet and interact with their neighbours. Or they might migrate and interact by collisions ( Figure 14B). Regularity would then somehow emerge from these meso-scale interactions.
Our demonstration that patterning exists across extensive areas known to have variations in sediment characteristics imply that such variations are of second-order importance to the process of drumlin formation. The concept of equifinality generates a problem of terminology, in that whilst drumlin-arrays (local patches of self-organized drumlins) arise from interactions as they evolve, it is also true that shapes resembling drumlins (call them drumlins if you like) may form as 'site-specific landforms', say by sediment smeared around a bedrock protuberance. The latter do not falsify the finding of regularity, they just add some disorder (defects) to the pattern, and as proposed elsewhere (Clark, 2010) such elements (obstacle-drumlins and drumlin-clones) are likely to be different from the more ubiquitous emergent drumlins. We predict that emergent drumlins formed in thick drift (> 15 m) exhibit stronger regularity than those in thin drift (< 5 m) because in the latter case many of the forms will merely be obstacle-drumlins pinned to underlying bedrock protuberances, and which might be randomly distributed.
Ribbed moraine and mega-scale glacial lineations have already been found to exhibit strongly preferred wavelengths (Dunlop and Clark, 2006;Spagnolo et al., 2014) suggesting that they are also patterned phenomena requiring an explanation of their regular spacing. This patterning similarity between ribbed moraine, mega-scale glacial lineation and drumlins along with the finding that they vary in a continuum of shape and scale (Ely et al., 2016) adds credence to idea that they are genetically-related rather than separate manifestations formed by different processes.
In addition to the question of how drumlins form, the question of why they form has rarely been addressed (Smalley et al., 2000). In considering why sub-glacial basal shear stresses are typically found to occupy a narrow physical range, Lliboutry (1958) hinted at an answer to both questions. He suggested that ice flow seeks to modify the bed roughness, smoothing out overly rough beds to reduce friction, but when the bed is too smooth, drumlins are produced to provide resistance to flow. More work is required to explore interactions between roughness and flow. The ultimate goal is a combined understanding of how sub-glacial conditions yield a bumpy bed and the extent to which this roughness modulates ice flow velocity. That drumlins, along with mega-scale glacial lineation and ribbed moraine, are found to occupy a set scale and pattern highlights the potential for using these predictable emergent scales of roughness to help inform the next-generation of 'sliding laws' to improve modelling of ice sheet and ice stream flow (Ritz et al., 2015).

Conclusions
Recognizing that drumlin fields are patchy and contain gaps and choosing an appropriate statistical method for analysing such an inhomogeneous distribution we have shown that whilst drumlins are sometimes randomly placed, their main state is not random. They exhibit a strong and statistically significant signal of regularity across scales of 100 to 1200 m and possibly also up to 10 000 m in range. Because this finding holds across different ice sheets and a variety of substrate geologies we conclude that regularity is a ubiquitous tendency within drumlin fields. This patterning demonstrates spatial self-organization in the bedforming process with drumlins as a specific emergent manifestation of sub-glacial sediment mobility. Elements of the pattern must have some physical interplay with other elements to control their spacing and so kilometre-scale interactions between drumlins must occur as they evolve, or interactions may arise as a consequence of growth or migration. Hypotheses or models that can explain the spacing of drumlins are required. We highlight three suggestions for such self-organization that are not mutually exclusive: instability in the coupling of flow-sediment flux-bed shape; local feedback between sediment mobility and relief; coarsening by growth or migration. Numerical modelling tracking bump evolution beyond the initial instability phase (Chapwanya et al., 2011) to explore coarsening might prove to be fruitful, and further observations of bedforms evolving beneath ice sheets (e.g. King et al., 2009King et al., , 2016Johnson et al., 2010) should lead to improved understanding.
That drumlins are now known to be patterned landforms, along with mega-scale glacial lineation and ribbed moraine, and that these have been found to vary in a continuum of shape and scale (Ely et al., 2016) suggests that a unifying theory might be appropriate to explain these sub-glacial bedforms. In such a view these named types are merely emergent variants of the same phenomenon, whose differences arise from some key parameter, such as ice velocity, cumulative basal sliding distance, sediment thickness or from histories of interaction as they grow or migrate. for their pointers regarding spatial analysis and R software. The authors are grateful to the referees for helping improve the paper.