Crack patterns over uneven substrates

Cracks in thin layers are influenced by what lies beneath them. From buried craters to crocodile skin, crack patterns are found over an enormous range of length scales. Regardless of absolute size, their substrates can dramatically influence how cracks form, guiding them in some cases, or shielding regions from them in others. Here we investigate how a substrate's shape affects the appearance of cracks above it, by preparing mud cracks over sinusoidally varying surfaces. We find that as the thickness of the cracking layer increases, the observed crack patterns change from wavy to ladder-like to isotropic. Two order parameters are introduced to measure the relative alignment of these crack networks, and, along with Fourier methods, are used to characterise the transitions between crack pattern types. Finally, we explain these results with a model, based on the Griffith criteria of fracture, that identifies the conditions for which straight or wavy cracks will be seen, and predicts how well-ordered the cracks will be. Our metrics and results can be applied to any situation where connected networks of cracks are expected, or found.


I. INTRODUCTION
From networks of cracks across the surfaces of planets [1][2][3][4] to artistic craquelure in pottery glazes and paintings [5][6][7], or bespoke fracture textures in thin films [8,9], crack patterns naturally occur on a wide variety of scales. They are even implicated in the formation of scales on crocodiles [10]. All these structures represent the outcome of failure processes that occur in a brittle medium supported by an uneven substrate. Here we show how the shape of the substrate plays an important role in determining the geometry of the resulting cracks. We will demonstrate how to use the substrate to template fractures, explore the conditions necessary to generate a variety of patterns, and introduce general metrics to quantify the orientational ordering of any set of connected cracks.
To introduce the concerns of our work we will first briefly present three examples of cracks on the scales of kilometres, centimetres, and micrometers. These motivational examples, and those shown additionally in Fig. 1, also aim to highlight the need for quantitative approaches to the analysis of fracture patterns.
On large length scales, radial graben are associated with buried craters on Mercury [3,4] and Mars [1,2]. These features are due to the repeated filling and cooling of lava in impact craters. The uneven stresses caused by the edge of the crater are believed to lead to wrinkle ridges above the crater rim [4]. Contraction cracks also form within the interior of the crater, leading to a mixture of polygonal and circular graben patterns [3,4]. The radial symmetry of the substrate imposes a characteristic structure to the such cracks, as in Fig. 1(a).
On a more intermediate length scale, craquelure ( Fig. 1(b,c)) is the result of differential strains between a coating and its substrate [5][6][7]14]. This may be due to tension in the framing of a painted canvas [5], for * lucas.goehring@ds.mpg.de  example, or changes in the relative humidity at which a painting has been stored [14]. Alternatively, it may result from an artist's choice of a pottery glaze with a higher coefficient of thermal contraction than the clay it covers. Such craquelure patterns are extremely hard to forge, and can thus aid in the attribution of paintings [6]. They have been shown to be influenced by local painting style, choice of substrate, including, for paint on wood, even the species of wood used [6,15]. The analysis of similar patterns, in dried blood ( Fig. 1(d)), is also being developed as a forensic tool, for example to help determine the relative humidity when the blood was spilt [12,16]. On the micrometer scale, recent attention has focussed on controlling cracks. For example, Nam et al. [8] showed how to selectively initiate, refract, and stop both straight and wavy cracks in silicon nitride thin films, vapour-deposited on silicon wafers, by adding notches, cutting stairs, or introducing interfacial layers between the wafer and film. Kim et al. [9] instead looked at soft substrates, with crack formation guided by strategically placed notches along grooves; their work aims to control failure in situations, such as flexible electronics, where it must be expected. They find the most regular results when they design for crack spacings that are near the natural crack spacing of the brittle film.
A common feature of all these situations is that the cracking layer is thin, compared to its lateral extent. Nonetheless, the layer thickness imposes the natural length-scale of such problems, and the spacing of the cracks in each case reflects this thickness [9,[17][18][19]. Here we study the model problem of contraction cracks in 'thin' elastic layers adhered to rigid substrates and show how the shape of the substrate also plays an important role in determining the geometry of the crack pattern.
For this end we choose the system of desiccation cracks in drying clays as a simple, yet representative, experiment. As clays dry they tend to shrink, uniformly. Socalled channel cracks open in the film, to relieve its growing strain energy (see e.g. [20][21][22]). The appropriate physical model for stress in a wet clay, poroelasticity, has an exact mathematical analogy with the thermoelastic model of how most solids shrink as they cool [23][24][25]. In biology, comparable conditions can also result from the differential growth of layered tissues [26][27][28]. This is known to lead to the formation of scales on the heads of Nile crocodiles [10] (Fig. 1(e)), for example, while a similar mechanism has been proposed to account for leaf venation patterns [26,27].

II. MATERIALS AND METHODS
Briefly, bentonite clay was dried over rigid patterned substrates, where it cracked; a schematic summary of the experiment is given in Fig. 2(a). Compared to the wet clay, the substrate is undeformable and unbreakable.
Each substrate consisted of a 20×20 cm 2 plate. Five sinusoidal plates were cut from acrylic blocks by computernumerical-control (CNC) milling, with a resolution of 200-400 µm. These plates were corrugated; their surface z varied sinusoidally along one direction, x, with a constant profile along the other, y, direction such that z = A sin(2πx/λ). The wavelengths λ and amplitudes A of the plates are given in Table I. Two pairs of plates had the same dimensionless ratios a = A/λ, in order to check the scale-independence of the fracture patterns. Two additional plates were prepared by 3D printing from    acrylic photopolymer, with a vertical resolution of 30 µm, and a contour resolution of 5 µm (4D Concepts). These plates had radial wave patterns, with surface profiles (see Table I) that varied periodically in the radial direction, r, away from the centre of the plate, such that z = A cos(2πr/λ). Finally, a flat acrylic plate was used for control experiments. All plates were smooth to the touch. For each plate 10 cm high acrylic walls were glued to each edge, to make a waterproof container. Slurries of potassium bentonite (Acros Organics K-10 montmorillonite) were prepared by mixing a set weight of bentonite powder with twice as much water, by mass. Slurries were then stirred until smooth and immediately poured evenly into a container, over a substrate. At this concentration the slurries behaved as fluids, rather than pastes. These conditions avoid the memory effect [29,30], whereby the flow or vibration of a yield-stress paste can affect the way it cracks (see Fig. 1(f)). We checked that this effect was absent by drying slurries over flat substrates, and finding isotropic crack patterns.
After pouring, the bentonite settled slightly, separating into a dense clay covered by a thin water layer. We measured the thicknesses of these submerged clay layers in flat-bottomed containers after 6 hours of settling, without significant evaporation, and found that the consolidated wet clay consistently contained 0.49 grams of bentonite per cubic centimetre. For slurries poured over the sinusoidal substrates we therefore used the total mass of the bentonite as a proxy for the average layer thickness H of the wet clay layer. All experiments were performed in the range of H from 3 to 15 mm.
Once a slurry had been poured, it was left to dry in a quiet, warm environment for 6-48 hours, depending on the amount of clay used. When drying is complete, the bentonite becomes lighter in colour; no further cracks appear after this time. For the linear wave plates, drying was accelerated by halogen heat lamps, which raised the surface temperature of the clay to ∼ 50 • C. As elevated temperatures would melt the printed plates, drying over the radial plates was enhanced by raising the room temperature to about 30 • C. During drying, the appearance and evolution of the cracks were imaged every 5 minutes using a digital SLR camera (Nikon D5000 series), mounted above the drying container. By comparing the initial wet height H to the dry height of various layers, we estimate a total vertical strain of 0.27, during drying.

A. Image processing
Images were pre-processed using ImageJ and Matlab, as shown in Fig. 3. Each image was cropped to isolate the cracked clay, then converted to grayscale by averaging across all colour channels, and contrast-enhanced. Next, a band-pass filter was applied to remove high-frequency noise (3-pixel cutoff to reduce single-pixel noise) and lowfrequency variations in lighting intensity (40-pixel cutoff, comparable to the width of the largest cracks). Finally, each image was thresholded to give a black-and-white representation of the crack network, as in Fig. 3 To characterise the crack patterns we will focus on the shape of the crack network, and on the relative positions and orientations of features such as the vertices between cracks. For some measurements, we thus chose to eliminate unconnected cracks and spurs from the binary images. Isolated cracks, and other noise, were removed by deleting objects below a threshold size. Every image was then thickened, by 3 pixels, to ensure that small breaks in the crack network, or en-passant cracks [31], were filled in. The result was then eroded to leave a skeleton consisting of one-pixel-wide cracks that connected various vertices or nodes, and which delimited isolated patches, or peds. Cracks that branched off the main network, but stopped without connecting to anything, were then removed. Both the binary images and the skeletons were used as the basis for further analyses.
Finally, we note here that the error bars given throughout this work are estimated by subdividing the processed    images into quarters, analysing each subsection separately, and then calculating the standard deviation of the four different results.

III. RESULTS AND DISCUSSION
We found that as the thickness of the brittle layer increases, its contraction cracks slowly lose information about their substrate. However, in order to understand the effect of an uneven substrate, we must first briefly summarise the base case of cracks in flat bentonite layers. For the flat control plate, isotropic crack patterns were seen for all layer thicknesses (H from 4-15 mm). Cracks appeared one by one, sequentially breaking the clay into smaller and smaller pieces, as in [32,33]. This process stops at some finite spacing, when the pattern saturates [18,25] -further drying increases the widths of the cracks in these 'mature' networks, as in [34,35], rather than create new ones. We confirmed that the natural crack spacing over the flat plate is between 2-3 times the wet layer thickness H, using the line-dropping method [36]. The results are compatible with the crack spacings reported in [36], but consistently smaller, as the dry heights of bentonite layers were reported in that study.
For clays dried over the sinusoidally varying substrates the typical sequence of crack patterns, after drying is complete, is shown in Fig. 4. A wavy pattern developed in thin layers, with the first cracks meandering back and forth as they travelled along the troughs of the underlying substrate. Various shorter cracks then filled in the pattern, as in Fig. 4(a). For intermediate thicknesses of clay regular straight cracks appeared over each ridge, followed by cracks between these primary cracks, and at right angles to them, like the rungs of a ladder. As in Fig. 4(b), additional short cracks also often formed along the troughs of the pattern, especially in the later stages of drying. For thick layers, as in Fig. 4(c), there was no obvious effect of the substrate on the crack pattern, which was essentially isotropic.
Linear elastic fracture mechanics has no inherent length-scales, other than those set by the geometry of the system under study. If this assumption holds, then the dependence of the fracture pattern on the shape of our substrates can be entirely captured by considering two dimensionless parameters, such as the dimensionless layer thickness, h = H/λ, and a dimensionless amplitude, a = A/λ, of the substrate's relief. Since in many situations, including ours, the spacing of thin-film cracks scales with the thickness of the cracking layer [9,[17][18][19], the parameter h should describe the relative match between the natural spacing of any cracks and the wavelength of the perturbations affecting their pattern. The dimensionless amplitude a describes the sharpness of the peaks on the substrate.
For the linear wave plates we categorised each fullydried crack pattern by visual inspection as either wavy, straight or isotropic. As summarised in Fig. 5, the different patterns divide up the expected dimensionless phase space into well-defined regions. At the boundaries between the different patterns were small areas of coexistence, for example between wavy and straight cracks. This coexistence represents two scenarios -either both patterns were observed within different regions of a single experiment, or were seen for the same h, but on different plates. For dimensionless layer heights in excess of about h = 1, only isotropic patterns were ever seen.
Although there are clear visual differences between the patterns shown in Fig. 4, there are no well-established methods by which they can be characterised. In order to analyse our patterns quantitatively we therefore need to develop new order parameters similar to, for example, the orientational order parameter of liquid crystals [37]. We present three such measurements next, based on Fourier methods, the relative orientation of the regions outlined by the crack network and the positions of neighbouring crack vertices.

A. Fourier analysis of crack patterns
The presence of evenly spaced cracks that are wellaligned with the substrate's features suggests the use of Fourier methods to quantify our crack patterns. Although a two-dimensional transform could be used to measure crack spacings, we focus here on a onedimensional approach, which allows us to measure the relative alignment and periodicity of the cracks with respect to the linear substrates. To this end, we first collapse each binary image to find the average density of cracks at each position along the x-direction. This is done by summing over the y-direction of the images, where cracks are represented by ones and uncracked material by zeros. We then find the power spectrum of the crack density as the absolute value of the square of its Fourier transform, normalised to have an integrated power of one.
In cases where cracks are aligned with the substrate there are sharp peaks in the power spectrum of the crack density when the wavenumber ξ x is an integer multiple of the wavelength λ of the substrate, as in Fig. 6. In other words, the substrate prefers there to be a integer number of cracks, per wavelength. The responses of all five linear wave plates were very similar: typically, most power was found in the first or the second peak, i.e. when ξ x λ = 1 or 2. Furthermore, there was a change in the dominant mode of the power spectra around h = 0.5−0.75. Below this, most power was in the ξ x λ = 2 peak, while above this the ξ x λ = 1 peak was also strong. Visually, the point where these intensities cross over coincides with the transition between patterns with cracks along every ridge ( Fig. 6(a), ξ x λ = 1), and patterns with cracks along both the ridges and troughs ( Fig.  6(b), ξ x λ = 2). The appearance of wavy cracks, at very small h, are associated with a decline in the power of the ξ x λ = 2 mode. Since they meander, their crack density is not as strongly localised as the ladder-like cracks, and the Fourier analysis reflects this. The isotropic cracks, around h = 1 or thicker, show no significant coupling to their substrate's periodicity.
An example of all of the above behaviour is given in Fig. 7 for plate 2, where a = 0.25. The results for plate 3 are consistent with this, and extend the range of suppressed ξ x λ = 2 down to h = 0.2. For plates 4 and 5, where the substrate is relatively sharper (a = 0.5), the sequence of responses is identical, but shifted to slightly higher values of h: the transition from one to two cracks per wavelength is at h = 0.75, for example. For plate 1, with shallow relief, this transition is shifted down to about 0.4-0.5, but the Fourier analysis does not produce very clear results, otherwise. These variations in the conditions leading to the different types of crack patterns are consistent with the visual classifications given in Fig. 5.
Thus, Fourier methods capture the regularity of the observed crack patterns, and quantify how the transition from one to two cracks per wavelength occurs when the periodicity of the substrate is about 1.5 to 2 times the thickness of the cracking layer -i.e. approaching the natural crack spacing of the clay over a flat substrate. However, it does not provide a good indication of the alignment of the various crack features with the substrate. The absolute intensity of the power spectral peaks, which could be such a measure, is too easily affected by noise to serve in this way. The next two sections aim to provide more robust measures of alignment.
B. An orientational order parameter As the clay dries in our experiments, its cracks subdivide and isolate individual regions of material, whose shapes are then determined by their sets of final bounding cracks. For example, in Fig. 4(b) the cracks around the dried peds are mostly parallel or perpendicular to the substrate relief. The cracked regions therefore have rectangular shapes, and tile the plane like rows of bricks. In contrast, the cracks in Fig. 4(c) are only weakly influenced by the substrate and define cracked regions with more arbitrary shapes and orientations. This section aims to measure the degree of alignment of these regions, with respect to the ripples along their substrate.
From each image we extracted a binary skeleton of the crack network (see Sect. II A). Each connected region in this image is identified and an ellipse is fit to cover its area. The major axis of this ellipse goes through the region's centroid and its orientation, θ, is chosen to minimise the second moment of area of the region, around the major axis. Figure 8 shows examples of elliptical fits of cracked regions with three different patterns.
To characterise the crack-bounded regions, through their best-fit ellipses, we use an orientational order parameter that measures the relative alignment of rod-like particles, or any objects that are invariant after rotation by 180 • . This parameter is well-known in the physics of liquid crystals (e.g [37]), where it is used to describe the transition from an isotropic to a nematic (orientational order without positional order) phase, for example. In our case, for dimension d = 2, the general definition [38] of this nematic order parameter simplifies to  where the angled brackets represent averaging over all ellipses. The orientation, θ, of the major axis of each ellipse is measured with respect to the y-axis for the linear wave plates, and the radial direction for the radially-symmetric substrates.
The orientational ordering for all our experiments is shown in Fig. 9. For each dimensionless amplitude a there is good agreement between the different experiments, although the radial plates give slightly higher values of S 1 than the linear plates. For small enough h, in particular around h = 0.5 and below, the peds are very strongly aligned with the relief of the substrate, and this alignment slowly decreases as the layer thickness increases, until S 1 becomes statistically equivalent to its value over the flat plate above a relative thickness of at most h = 1.25. Furthermore, the intensity of the ordering decreases with the relative sharpness of the plates, a, for experiments with the same h. However, unlike the Fourier analysis, or the order parameter we introduce in the next section to describe the relative alignment of the individual cracks, S 1 remains high for very low h. It is a sensitive measure of the straight-to-isotropic transition, but because it averages over the entire shape of a cracked region, it does not detect the meandering of the cracks.

C. Crack alignment order parameter
A second order parameter is introduced to describe the relative alignment of the cracks with their substrate. Working from the skeletons of the cracks, we first find the coordinates of all vertices, or intersections between cracks. These skeletons form closed networks, where all isolated features and spurs have been removed. The vertices are therefore points along the skeleton that touch at least three different cracked regions, and this property is used to identify them. Crack segments are then defined by the paths between any two neighbouring, connected vertices, and their orientation φ by the vectors connecting these two endpoints. Since one end of a crack cannot be effectively distinguished from the other, the distribution of crack angles is collapsed onto the range from 0-180 • . These angles are referenced to the y or radial directions, for the linear or radial wave plates, respectively.
It should be emphasised that this method creates a very local definition of the orientation of a crack, and splits long cracks into many individual segments, between adjacent intersections with other cracks. This is different from some other network analyses of cracks [33,39,40], which focus instead on the entire length of a crack, regardless of any intersections it has with other features. An example of a crack pattern, and the relative abundance of crack orientations in it, is given in Fig. 10.
In our more well-ordered experiments, as in Fig. 10, there are typically prominent cracks along the ridges of the substrate, and shorter perpendicular cracks that join these together. Their measure must reflect this symmetry. In analogy with S 1 , it should take a value of +1 for cracks that are either parallel or perpendicular to some reference direction, 0 for random patterns, and −1 for cracks that are all completely mis-aligned with the reference (here, cracks at ±45 • to the substrate features). The simplest such definition is for a crack alignment order parameter S 2 = cos(4φ) . (2) We calculated this parameter for each drying experiment, and summarise the results in Fig. 11. We found that the crack alignment order parameter, S 2 , matches the orientational order parameter of the cracked regions, S 1 , in most situations. There is also good data collapse between the linear and radial wave plates, and between all plates of the same relative amplitude a, confirming that our results do not depend on the absolute scale of the drying experiments. There is a loss of ordering around h = 1, again at slightly higher values of h for the sharper substrate ridges, and slightly lower values for the shallower relief. However, unlike S 1 , there is a clear downturn in S 2 at low layer thicknesses -the cracks become less aligned for very thin layers. As it is a local measure, the crack orientation is more sensitive than our other metrics to any curvature of the cracks, and some of the decrease in S 2 for very low h is due to the wavy cracks. Nevertheless, in all cases S 2 starts to   fall well before we can see wavy cracks by eye. In Section IV we will show how this effect, and the other transitions in the crack pattern, can be explained by the relative forcing of the uneven substrate, which will be shown to be strongest at the intermediate layer thicknesses where both order parameters are largest.

D. Extension to other crack networks
The two order parameters S 1 and S 2 introduced above describe, respectively, the relative alignment of crackbounded areas, or individual cracks, with respect to some known angle. In our experiments there was an obvious direction to choose for this angle, based on the substrates used -either the local radial vector, or one of the cartesian axes. However these metrics can also be used to determine whether a crack network is aligned or not without knowing, in advance, any such preferred direction. By extending the liquid-crystal analogy of the order parameter S 1 , one may use the local director, or average local orientation of the relevant objects, as the reference direction. In other words, one may first calculate the average θ or φ, either over all cracks or areas in a region of interest, or locally over neighbourhoods that are large enough to provide good statistical averaging. Then, in calculating the order parameters through Eqs. 1 and 2, all angles can be referenced to this local director orientation. This would allow one to compare the relative degree of alignment of natural crack patterns like polygonal terrain [41,42] and the hydraulic channels of heavily fractured rock [39,43], for example. It would also allow one to better measure manufactured patterns such as craquelure [6,15], or custom-designed micro-cracks [8,9].

IV. MODELLING THE CRACK PATTERN TRANSITIONS
When a crack opens, it releases stored strain energy from its environment. In Griffith's model of fracture this energy is used to create the new surface area of the crack, as it grows [44,45]. In this section we consider a finite element model of a channel crack opening over a sinusoidal surface, and use the Griffith criteria for fracture to predict where the first cracks will appear. We will show how the observed transitions from wavy to straight to isotropic cracks can be explained in terms of the relative strain energy released by a crack opening at different locations in a thin elastic layer drying -or cooling, or with some other misfit strain -over such a rippled substrate.
In particular, we consider a linear poroelastic [46] model for cracks in drying clay. The total stressσ in the clay is a sum of the elastic stress σ supported by the network of clay particles and the capillary pressure p of the interstitial liquid, such that The negative sign here accounts for the convention that a positive stress is tensile, while a positive pressure is compressive. The pressure is generated by the tiny curved menisci of the water still trapped in the clay, and for these conditions p < 0. The network stress can be given in terms of the strain by the linear elastic constitutive . 12: Numerical solution of a channel crack opening over (a) a trough or (b) a peak of the substrate. The deformations caused by the cracks are exaggerated, so as to be visible. The greyscale shows how much of the pre-crack stress σ * xx is released by the crack -with black shown where all the stress is released, and white where the crack does not affect the stress at all.
where E is the Young's modulus and ν the Poisson's ratio of the clay. The strain tensor is, in turn, derived from the elastic displacements, u, of the clay by The pressure term in Eq. 3 provides the driving force for fracture. As the clay dries p becomes more and more negative (its magnitude increases), requiring the clay to compensate with a tensile strain that can lead to crack growth. In order to instead model cooling, or tissue growth (e.g. for crocodile skin [10] or leaves [27]), this pressure can be replaced by an equivalent thermal contraction term [23,24], or the effects of differential growth in a brittle layer and its substrate. We solved Eqs. 3-5 using Matlab's finite-element routines, on domains like those shown in Fig. 12, namely semi-rectangular regions with straight edges at the top (z = 1, dimensionless units, scaled so that the average film thickness H = 1), left (x = 0) and right (x = 10) sides, and a corrugated base where z = 1 − = A cos(2πx/λ + ψ). Here ψ is an offset phase, allowing us to examine the strain energy released by fracture at various positions along the substrate, and is the local film thickness. Plane strain conditions, with no displacements allowed in the y-direction, were used to model a channel crack advancing with a constant profile along the sinusoidal plate. The upper surface was taken to be traction free (σ xz = σ zz = 0 on z = 1), while the lower surface was treated as no-slip, such that u = 0 there, and only vertical displacements were allowed (u x = 0) along the right boundary. The left-face boundary was either treated like the right one, for an intact film, or as a traction-free boundary, for a cracked film. Finally, the pressure p was taken to be constant everywhere, as appropriate for quasi-static fracture.
The work done by the growth of the channeling crack can be measured by integrating the product of the precrack stress, σ * , and the opening displacements of the crack, δ, as in [21,45]. This is, essentially, the sum of the forces on the incipient crack surfaces, acting over the distance they travel as the crack opens. The strain energy released, per unit area of new crack, is thus For a simulation with a specific A and λ we calculated the pre-crack stress field once, using ψ = 0, and used this solution to find σ * at different positions along the substrate. We then calculated the opening displacements of a crack at different ψ. Accounting for the opening displacement to either side of a crack can be done by letting δ x = u x (ψ) + u x (−ψ), and δ y = u y (ψ) − u y (−ψ). This gives us the values of G for cracks growing along the substrate at different ψ, ranging from the crest of the sine wave, to its trough. For a flat substrate, this problem has been studied in depth by Beuth [21], Xia and Hutchinson [22] and Yin et al. [47,48]. They showed that G scales as where g is a dimensionless parameter that only depends on the elastic mismatch between the cracking film and its substrate, and whereĒ = E/(1 − ν 2 ) and σ 0 = −p(1 − 2ν)/(1 − ν). For example [21], for ν = 1/3, g = 0.71. We confirmed their solutions with respect to the exponential decay of the height-averaged stress away from the crack [22] and the dependence of g on ν [21] and the depth of fracture [21,47], and checked that g was independent of the average film thickness, the film's elastic modulus E and the capillary pressure p. The mesh used, and the horizontal domain length, were found to be sufficient to agree with these known solutions to within 1%. For an uneven substrate G depends on location, measured here by the offset ψ. As the clay dries the relative shape of G(ψ) will not change, but rather its magnitude will simply increase along with the capillary pressure. The Griffith criteria predict that fractures can only grow once G reaches some critical value, G c , which is a material parameter (see e.g. [45]). The first cracks in the brittle drying layer should therefore form at the positions where the strain energy released by these cracks, i.e G, would be a maximum, as these are the first places where G = G c . For a Poisson ratio of ν = 1/3, we calculated how G varied with position in a wide range of experimental conditions, from A = 0.01-0.90 and λ = 0.5-8.0. To focus only on geometric effects we used Eq. 7 to rescale G by the relative stiffness 2Ē/πσ 2 0 , so thatḠ = g. The results are summarised in Fig. 13(a), which shows the values of ψ that maximise G for various A and λ. Depending on the shape of the substrate and the thickness of the clay, the preferred positions of the earliest cracks  can be either (i) over the peaks of the substrate, (ii) over the troughs or (iii) at pairs of points displaced to either side of each peak. Figure 13(a) also replots the observations from Fig. 5, corresponding to the conditions where wavy and ladderlike cracks are found. For all experiments where straight cracks were seen the model correctly predicts that the first cracks should appear over the ridges of the substrate. Where wavy cracks are seen, the model instead generally predicts two stable positions for cracks, along the sides of the ridges. Although our model has no dynamic component, this empirical agreement of the division of phase-space suggests that the wavy cracks result from cracks moving back and forth between their two bistable positions. The boundary between the straight and wavy cases corresponds to a forward pitchfork bifurcation, as sketched in Fig. 13(b).
To consider the more gradual transition from ladderlike to isotropic cracks, we looked at the relative variations inḠ, defining ∆Ḡ as the difference between the maximum and minimum values ofḠ(ψ), for each value of A and λ studied. If ∆Ḡ is small, then there is only weak guidance for the cracks, and they should appear more random than in cases of larger variation of energy release rate. We tested this hypothesis by comparing ∆Ḡ with  the crack alignment order parameter of our experiments. First, we calculated ∆Ḡ for each point in our simulation, mapping the results in Fig. 14 alongside the data showing S 2 for our various experiments with linear wave plates. This figure shows that there is a neighbourhood of driving parameters, centred around λ/H = 2.2 and A/H = 0.6, where ∆Ḡ is largest. Our experiments pass through this region, and show their most well-ordered patterns as they do so. To quantify this agreement, we then calculated ∆Ḡ for the conditions corresponding to each experiment where ladder-like or isotropic cracks were seen. As shown in Fig. 14(b), there is indeed a good correlation between the range of energy release rates allowed for any particular geometry, and the relative alignment of the cracks that formed over it.
Finally, we note that some insight into the limits of our results can be gained through a qualitative comparison to channel cracks in flat films and coatings. Here the stresses relaxed by the crack will decay away over a length comparable to the thickness of the film [20,22]. This effect can also be seen in Fig. 12. For a wavy substrate this limit is approached when the amplitude of the relief goes to zero, and under such conditions isotropic cracks with a spacing ∼ H are clearly expected. However, for finite A/H the limit can still offer insight. For very long wavelength ripples, large λ/H, each crack will see a local environment, of size ∼ H, that is approximately flat. Although these conditions were not probed by our experiments, we would also expect an isotropic pattern in this case, but one with a local modulation of the crack spacing reflecting the changes in local thickness. For the opposite limit, of very sharp features where λ/H 1, the region around any crack will sample many ridges and troughs, and these local variations will tend to cancel each other out (see Saint-Venant's principle e.g. [49]); again, isotropic cracks would be expected. In other words, for finite-amplitude relief a brittle layer must filter out the high and low-frequency information of its substrate, and any contraction cracks that form should reflect this.

V. CONCLUSIONS
Our results explore the effects of an uneven substrate on the development of contraction cracks in thin films, and are inspired by natural patterns of cracks around buried craters, paintings, and vapour-deposited thin films. We showed that the wavy, periodic undulation of a substrate can cause the ordering of a crack pattern that forms above it, and has the strongest influence on this pattern when the thickness of the cracking layer is comparable to, or slightly smaller than, the periodicity of the substrate.
Generally we found that for increasing layer thickness the crack pattern changed from one with wavy cracks running along the troughs of the substrate, to straight cracks appearing over each peak and trough of the substrate, to straight cracks on only the peaks of the substrate, and finally to a disordered or isotropic pattern. The wavy-to-straight transition was the most abrupt change, as evidenced by examples of single experiments that were wavy on one side of the plate, and straight on the other. All these patterns were quantified by Fourier methods and orientational order parameters, similar to the nematic order parameter of liquid crystals, which are easily adaptable to other types of crack networks. Finally, we showed how the transitions between pattern type could be explained by a simple finite-element model of channeling cracks advancing over uneven relief, regardless of whether those cracks were driven by drying, cooling or the differential growth of crocodilian skin.