UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Chaotic pattern dynamics on sun-melted snow Mitchell, Kevin A. 2008

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2008_fall_mitchell_kevin.pdf [ 1.76MB ]
JSON: 24-1.0066733.json
JSON-LD: 24-1.0066733-ld.json
RDF/XML (Pretty): 24-1.0066733-rdf.xml
RDF/JSON: 24-1.0066733-rdf.json
Turtle: 24-1.0066733-turtle.txt
N-Triples: 24-1.0066733-rdf-ntriples.txt
Original Record: 24-1.0066733-source.json
Full Text

Full Text

Chaotic Pattern Dynamics on Sun-melted Snow by Kevin A. Mitchell  B.Sc., The University of British Columbia, 2005  A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in The Faculty of Graduate Studies (Physics)  The University Of British Columbia (Vancouver) October, 2008 c Kevin A. Mitchell 2008  Abstract This thesis describes the comparison of time-lapse field observations of suncups on alpine snow with numerical simulations. The simulations consist of solutions to a nonlinear partial differential equation which exhibits spontaneous pattern formation from a low amplitude, random initial surface. Both the field observations and the numerical solutions are found to saturate at a characteristic height and fluctuate chaotically with time. The timescale of these fluctuations is found to be instrumental in determining the full set of parameters for the numerical model such that it mimics the nonlinear dynamics of suncups. These parameters in turn are related to the change in albedo of the snow surface caused by the presence of suncups. This suggests the more general importance of dynamical behaviour in gaining an understanding of pattern formation phenomena.  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  v  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  vi  Abstract  List of Tables  Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . viii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Non-dimensionalisation . . . . . . . . . . . . . . . . . . . . 2.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . 2.3 Linear Stability . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Roughening and Pattern Formation . . . . . . . . . . . . . 2.5 The Early Nonlinear Regime . . . . . . . . . . . . . . . . . 2.6 The Long-time Nonlinear Regime . . . . . . . . . . . . . . 2.7 Power Spectral Density in the Long-Time Nonlinear Regime 2.8 Saturation Height . . . . . . . . . . . . . . . . . . . . . . . 2.9 Effect of the |∇η|2 Term on the Mean . . . . . . . . . . . . 2.10 Diffusion of the Simulated Suncups . . . . . . . . . . . . .  . . .  4 4 5 6 6 12 13 15 27 29 32  3 Field Observations . . . . . . . . . . . . . . . . 3.1 Measuring the Local Peak-to-Peak Height . . 3.2 Roughening: Time and Height scale . . . . . 3.3 Characteristic Length . . . . . . . . . . . . . 3.4 Diffusion of Real Suncups . . . . . . . . . . . 3.5 Contribution of Suncups to the Ablation Rate  . . . . . .  40 41 44 47 54 60  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  64  4 Conclusion  . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . . .  1  iii  Table of Contents Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  66  Appendices A Numerical Methods . . . . . . . A.1 Spectral Differentiation . . . . A.2 Exponential Time Differencing A.3 Adaptive Time Stepping . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  68 68 69 73  B Perspective Transformation . . . . . . . . . . . . . . . . . . B.1 Homogeneous Coordinates . . . . . . . . . . . . . . . . . . B.2 Matrix Formulation of the Perspective Transformation . . B.3 Inversion of the Perspective Projection . . . . . . . . . . . B.4 Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . B.5 Changes to the Parameters Over the Course of Observation B.6 Setting the Scale . . . . . . . . . . . . . . . . . . . . . . . . B.7 Cross-sections . . . . . . . . . . . . . . . . . . . . . . . . .  . . . . . . . .  77 78 78 80 81 82 84 86  iv  List of Tables 3.1  Measured scales and the coefficients for Eq. 1.1 derived from them . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  63  v  List of Figures 1.1  Top view of suncups . . . . . . . . . . . . . . . . . . . . . . .  2  2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10 2.11  8 9 10 11 14 15 16 17 18 19  2.25  Simulated RMS Roughening . . . . . . . . . . . . . . . . . . . PSDs of surface in the linear regime . . . . . . . . . . . . . . Pattern formation top view . . . . . . . . . . . . . . . . . . . Pattern formation cross section . . . . . . . . . . . . . . . . . Early nonlinear saturation cross section . . . . . . . . . . . . Early nonlinear saturation PSD . . . . . . . . . . . . . . . . . Early nonlinear saturation top view . . . . . . . . . . . . . . . Nonlinear height saturation . . . . . . . . . . . . . . . . . . . Long-time nonlinear solutions top view . . . . . . . . . . . . . Long-time nonlinear cross-sections . . . . . . . . . . . . . . . Time and ensemble averaged radial power spectral density in the nonlinear regime . . . . . . . . . . . . . . . . . . . . . . . Determination of the characteristic length . . . . . . . . . . . Time and ensemble averaged cross-sectional power spectral density in the nonlinear regime . . . . . . . . . . . . . . . . . Modified cross-sectional PDS’s . . . . . . . . . . . . . . . . . Modified radial PSDs . . . . . . . . . . . . . . . . . . . . . . . Characteristic length as a function of θ . . . . . . . . . . . . . Mean RMS saturation height as a function of θ . . . . . . . . Local peak-to-peak height as a function of θ . . . . . . . . . . Time evolution the local peak-to-peak height . . . . . . . . . |∇η|2 as a function of θ . . . . . . . . . . . . . . . . . . . . Trajectories of the minima of the numerically simulated suncups Density of suncups per characteristic length . . . . . . . . . . Creation/Destruction rate of suncups as a function of θ . . . RMS displacement of numerical suncup minima as a function of time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Diffusion coefficient for suncup minima as a function of θ . .  3.1  Map of the snow observation sites . . . . . . . . . . . . . . . .  41  2.12 2.13 2.14 2.15 2.16 2.17 2.18 2.19 2.20 2.21 2.22 2.23 2.24  21 22 23 25 26 27 28 30 31 33 34 35 36 38 39  vi  List of Figures 3.2 3.3 3.4 3.5 3.6  Selkirk observation site . . . . . . . . . . . . . . . . . . . . . . Whistler webcam observation site . . . . . . . . . . . . . . . . Experimental cross-sectional surface profiles . . . . . . . . . . Numerical height roughening scaled to experimental data . . Height and time scales as a function of θ from fits to experimental data . . . . . . . . . . . . . . . . . . . . . . . . . . . . PSDs of snow surface cross-sections . . . . . . . . . . . . . . . Modified PSDs of the snow surface cross-sections . . . . . . . Cross-sectional numerical PSDs at 5 days . . . . . . . . . . . Length scale of numerical solution with respect to experimental cross-sections . . . . . . . . . . . . . . . . . . . . . . . . . Experimental suncup trajectories in isotropic “snow coordinates” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . RMS displacement of the experimental and numerical suncup minima . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Variation in the experimental determination of θ over the observed range of ∆t . . . . . . . . . . . . . . . . . . . . . . . Net surface displacement as a function of time . . . . . . . .  58 61  B.1 Perspective transformation . . . . . . . . . . . . . . . . . . . B.2 Snow surface image PSD . . . . . . . . . . . . . . . . . . . . . B.3 Snow surface cross-section photo . . . . . . . . . . . . . . . .  83 85 87  3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14  42 43 45 46 48 50 51 53 54 56 57  vii  Acknowledgements I would like to thank my colleagues in the MBE lab for providing support and encouragement in the completion of my programme, especially Anders Ballestad who showed me the way of the theorist in an experimental research group and Arvin Yazdi, Dan Beaton, and Mike Whitwick who in particular assisted in the development of my hands-on experimental skills. I would also like to thank Matt Choptuik and Ralf Wittenberg who taught and guided me in successfully implementing my numerical solutions. My second reader, Joerg Rottler also deserves his share of credit for finding the time in his busy schedule on short notice to read this thesis. Most of all however, I would like to thank Tom Tiedje for believing in my ability to do scientific research, fighting for me at the administrative level and supporting me financially even before my admission to the Masters programme.  viii  Chapter 1  Introduction Suncups are a type of spontaneous pattern formation found on snow surfaces exposed to extended periods of direct sunlight. They are characterised as quasi-periodic parabolic depressions in the snow surface with rounded bowllike valleys that can range 20−80cm in diameter [1] as shown in Fig. 1.1. The bowls are joined at sharp peaks ranging between 2 − 50cm above the bottom of the valleys [17]. As their formation requires an undisturbed snow pack that is deep enough to persist under prolonged melting, they are typically found in temperate alpine regions during the summer. Suncups are important because they are an effect of the snow melting process. They play an active role in this process since they reduce the albedo of snow thereby increasing its melt rate. Snow melt is the source of many of the worlds rivers, supplying fresh water for drinking, industry and agriculture. At the same time, snow melt plays a significant role in the shrinkage of temperate glaciers. In a climate of uncertainty about fresh water resources and rising global temperatures, a better understanding of suncups and their role in melting is of particular contemporary interest. Suncups are also interesting in their own right as an example of spontaneous pattern formation in nature. Many such examples exist. Penitentes are a qualitatively different snow ablation feature with much taller and sharper peaks that occur when humidity is low and sublimation is important [3]. Aeolean sand dunes are caused by the transport of individual sand grains by wind [12]. Sastrugi are their analogue in snow, typically found in the polar regions where the snow pack is persistent and dry [22]. Frost boils are caused by frost heaving [14]. Regmaglypts are suncup-like surface patterns that form on meteorites under atmospheric ablation [11]. Similar pattern formation can also be observed in materials science on a vastly different length scale. Submicron sputter ripples are produced on sputtering targets in ion beam erosion [13]. Understanding how suncups form and evolve can hopefully shed some light on these other physical processes as well. By considering the interaction of light with a granular medium whose local decrease in height is proportional to absorption of sunlight, the following 1  Chapter 1. Introduction  Figure 1.1: Top view of suncups as seen from the chairlift to the summit of Whistler Mountain (Photo: Tom Tiedje, July 2007) PDE model for the formation and evolution of suncups was derived [21] ∂t h(x, t) = −c0 − c1 ∇2 h − c2 ∇4 h − c3 |∇h|2 + c4 ∇2 |∇h|2  (1.1)  Here, h is the local height of the surface at time t and position x in the plane of the surface. The constants c0,1,2,3,4 > 0 are related to the spectrally averaged diffusion length of light immediately beneath the snow surface which can be determined from the optical properties of snow as documented in [23]. The c0 term is the trivial rate of recession of the mean surface height for a flat snow surface. The remaining differential terms describe how this rate is locally increased or decreased by the surface shape. As will be shown in Section 2.3, the first linear term ∇2 h is unstable and dominant at long length scales while the second linear ∇4 h is stable and dominant at short length scales. The crossover point between the two terms defines a characteristic length with the fastest exponential growth rate in the linear regime. Setting c4 = 0 in Eq. 1.1 results in the Kuramoto-Sivashinsky equation originally derived to describe laminar flame fronts [9] and reaction diffusion systems [20]. The nonlinear term |∇h|2 results in a saturation of the RMS height h 2 to an approximately constant characteristic height scale which is consistent with the behaviour of suncups. After this saturation takes place, 2  Chapter 1. Introduction solutions fluctuate chaotically with time. The nonlinear term |∇h|2 with the negative sign as in Eq. 1.1 also gives solutions a suncup-like shape with rounded valleys and sharp peaks. The addition of the ∇2 |∇h|2 term to the Kuramoto-Sivashinsky equation maintains the height saturation and suncup shape of the solutions, though it does lead to an increase in the characteristic length [16]. More importantly however, it is found that this additional nonlinear term causes the chaotic fluctuations to slow down, or their timescale to be increased. Since the chaotic fluctuations only take place in the height saturated nonlinear regime, the strength of this nonlinear term provides a way of adjusting the nonlinear timescale of the solutions independent of the growth rate predicted by the linear terms. An additional distinction between the effects of ∇2 |∇h|2 and |∇h|2 is also noteworthy. All the differential terms in Eq. 1.1 with the exception of |∇h|2 have a zero mean. Since |∇h|2 is the square of a real value however, it is everywhere greater than or equal to zero and so its mean is nonzero and positive. Given a negative sign as in Eq. 1.1, this leads to an increase in the melt rate relative to the c0 observed for a flat surface. Some of this increased ablation rate can be accounted for by a corresponding increase in absorbed energy (or coalbedo). Determination of c3 as well as |∇h|2 then constitutes a measurement of the degree to which suncups increase the energy absorption of the snow. This thesis investigates the growth and fluctuation of suncups through numerical solution of Eq. 1.1 described in Chapter 2 and time-lapse observations of snow surfaces described in Chapter 3. The methods used to solve Eq. 1.1 are described in Appendix A while the techniques used to process the raw experimental data are described in Appendix B. By comparing simulations and observations, the coefficients c0,1,2,3,4 are uniquely determined. This determination demonstrates the necessity of a nonzero c4 and c3 to adjust the timescale of the fluctuations observed in real suncups. Thus Eq. 1.1 is shown to be the minimal PDE model necessary to quantitatively reproduce the nonlinear behaviour of suncups. The value of c3 is also used to estimate the increase in the coalbedo due the presence of suncups.  3  Chapter 2  Numerical Results In the interests of better understanding the properties of Eq. 1.1, numerical simulations are carried out. The methods used for these simulations are described in Appendix A. This chapter describes the most salient findings of these simulations particularly with respect to the observations of snow which will be discussed in Chapter 3.  2.1  Non-dimensionalisation  In the interests of reducing the parameter space that must be explored in the simulations, Eq. 2 may be cast in dimensionless form by measuring time, length and height respectively in units of t0 , x0 and h0 where t0 = x0 = h0 =  c2 c21  (2.1) c2 c1 c3 c1  (2.2) 2  +  c4 c2  2  −1/2  (2.3) (2.4)  This is be achieved by substituting t = τ t0 x = χx0 and h = ηh0 so that  ∂τ η = −  c0 t 0 − ∇2χ η − ∇4χ η + cos(θ)|∇χ η|2 + sin(θ)∇2χ |∇χ η|2 h0  where θ = cot−1  c2 c3 c1 c4  (2.5)  (2.6)  is the only free parameter which is taken to lie in the range θ ∈ [π/2, π] for consistency with the shape suncups. This is the only parameter that must be 4  Chapter 2. Numerical Results varied to obtain qualitatively different solutions to Eq. 1.1. If solutions for all values of the parameter θ are somehow known, then solutions to Eq. 1.1 may be obtained for all values of c1,2,3,4 simply by scaling the time, length and height appropriately. Additionally, the coordinate system may be shifted to one that moves with the recession of the surface at the rate c0 t0 /h0 . This is achieved via the transformation η → η − c0 t0 /h0 t which eliminates the first term in Eq. 2.5. ∂τ η = −∇2χ η − ∇4χ η + cos θ|∇χ η|2 + sin θ∇2χ |∇χ η|2  (2.7)  The dependence of the three scales in Eqs. 2.1-2.3 on the original four parameters c1,2,3,4 is not a unique recipe for reducing the number of free parameters in Eq. 2.7 to 1. The present scaling relations however are selected so that the relative strength of the nonlinear terms may be varied to both extremes while the sum of squares of the coefficients is constrained to unity. Throughout the rest of this chapter the χ subscripts will be dropped from ∇ for simplicity of notation.  2.2  Boundary Conditions  The physical boundaries of the snow surface are beyond the scope of the optical model for suncups. Mathematically, however there is no sidestepping the issue. In many ways, it seems natural to take χ ∈ (−∞, ∞) ⊗ (−∞, ∞), (i.e., χ extends infinitely in all directions along the snow surface). For large snow fields where the physical boundaries are far from a typical point on the snow surface, this would seem reasonable. Such an infinite domain however is not possible to simulate numerically. The closest that might be achieved is a large periodic domain χ ∈ [0, L) ⊗ [0, L). The periodicity is trivially achieved with the spectral method described in Appendix A. This might in some ways even be more appropriate a model for the actual snow surface as it prevents the fluctuations of too large a length scale from occurring. As will be seen below, solutions to Eq. 2.7 indeed have significant fluctuations that extend to length scales as large as the domain size will allow. It is unlikely that such long wavelength fluctuations represent any physical property of the snow surface and at the very least they are not of interest for the present work. Thus, the size of the domain L sets an upper limit on the maximum resolvable length scale of the numerical simulation just as the mesh spacing ∆χ sets a lower limit on the minimum resolvable length scale. This analogy suggests that in the same way as numerical results should not significantly 5  Chapter 2. Numerical Results dependent on ∆χ, neither should they depend on L. In the present work, care is taken to ensure that both of these conditions are met.  2.3  Linear Stability  Though general solution of the nonlinear Eq. 2.7 is not feasible, it can easily be seen that the equation does admit the trivial solution η(χ, τ ) = 0. It is of interest then to examine the stability of this steady state by perturbing it by a small amount δη with low slope |∇δη| ≪ 1. In this case, the nonlinear terms become vanishingly small compared to the linear terms ∂τ δη = −∇2 δη − ∇4 δη + O(|∇δη|2 )  (2.8)  This equation may be solved by taking the Fourier transform ∂τ δ ηˆ(κ, t)) ≈ κ2 δ ηˆ − κ4 δ ηˆ  (2.9)  where κ = {2πL 1/n 1/m : m, n ∈ Z} is the wave number and κ = |κ|. The Fourier space solution to this linearised equation is given by δ ηˆ(κ, t) ≈ δ ηˆ(κ, 0)e(κ  2 −κ4  )τ  (2.10)  which remains a good approximation until the condition |∇δη| ≪ 1 is no longer met. From this, it is clear that in the small slope regime where the approximation is valid, there is exponential growth at length scales for which κ < 1 (ℓ > 2π). Length scales smaller than this on the other hand, are exponentially damped. This loosely suggests that even if only the weaker condition |δη| ≪ 1 is satisfied initially, the damping of small scale fluctuations will quickly ensure that |∇δη| ≪ 1 is also satisfied. Perhaps most significant√is that it can also be seen that there is a length scale at κ = 1/2 (ℓ = 2π 2) at which the exponential growth is a maximum with an exponential growth time constant of τc = 4.  2.4  Roughening and Pattern Formation  To simulate the spontaneous pattern formation observed in snow, Eq. 2.7 is solved numerically with low amplitude random initial conditions. It is clear from the previous section that such perturbation from a totally flat surface is necessary since η(χ, τ ) = 0 is a steady state solution. This is  6  Chapter 2. Numerical Results reasonable given that even a relatively flat snow surface will have some small fluctuations. The initial condition grid is initialised using the GNU libc implementation of the POSIX standard rand() function to randomly compute each grid point in real space. The resulting solution is normalised to a maximum peak-to-peak amplitude η(χ, 0) ∞ = a0 ≪ 1 and set to zero mean η(χ, 0) = 0. The rand() function deterministically returns pseudo-random numbers based on a given random seed. There exist pseudo random number generators with far stronger randomness (i.e., less correlation between successive numbers generated), but these require linking additional libraries, reducing code portability. This is deemed unnecessary for the present case since there is no reason to believe that adjacent points on the snow surface are truly randomly uncorrelated. Indeed, this is almost certainly not the case. This of course, is not to say that any correlation between these weakly pseudo random numbers corresponds to that of a snow surface. The goal here is not to demonstrate that pattern formation arises from a surface matching that of snow in great detail, but that it is essentially independent of the particular initial condition used. This might be predicted given Eq. 2.10, which says that any low amplitude initial condition will quickly develop a strongly dominant characteristic length regardless of its particular form. Thus, for the present work it is deemed sufficient that the initial condition contain a statistically uniform distribution of fluctuations at all length scales (i.e., white noise). The GNU rand() function is found sufficient to provide this as seen in the τ = 0 trace of Fig. 2.2. The roughening of the RMS surface amplitude η 2 from an initial condition generated as described above with a0 = 1 × 10−6 is shown in Fig. 2.1. Solutions for three values of θ ∈ [π/2, π] are shown demonstrating that the RMS height during the initial exponential roughening stage is essentially independent of the nonlinear coefficients. For this reason, it is called the linear regime. The radial power spectral densities (PSD) for θ = 0.75π at key points in the linear regime are shown in Fig. 2.2. The radial power spectral density at κ is defined as the average of the 2D PSD on a circle of radius κ centred at κ = 0. P SDr =  L2 ˆ k|)|2 |h(| (2π)2 N 4  φ  (2.11)  Colour coded top views of the real space surface are shown at these same points in Fig. 2.3 with cross-sections in Fig. 2.4 for reference.  7  Chapter 2. Numerical Results  100  η  2  10−2 10−4 10−6 10−8  θ = 0.50π θ = 0.75π θ = 1.00π 0  20  40  60  80  100  τ (a) Roughening of the RMS surface amplitude as a function of time  η  2  10−6  10−7 0  0.5  1  1.5 τ  2  2.5  3  (b) Fast initial smoothing. Axis limits are indicated by the boxed area in (a)  Figure 2.1: Roughening in the RMS surface amplitude for three values of θ at two timescales: (a) medium timescale at which the RMS amplitude shows exponential growth indicated by the fitted line which has time constant τc ≈ 4.1 (b) short timescale at which rapid smoothing takes place for short length scales. 8  Chapter 2. Numerical Results  100 τ = 0.0 τ = 0.7 τ = 40.0  PSDr  10−10  10−20  10−30  0  0.5  1  1.5  2 κ  2.5  3  3.5  4  Figure 2.2: Power Spectral Density of the surface for θ = 0.75π at key points in the linear regime. Solutions for different values of θ are qualitatively similar, particularly at long length scales κ < 1. The τ = 0 trace shows the PSD of the white noise initial condition. The PSD of the surface after the fast initial smoothing is shown at τ = 0.7. Finally, the τ = 40 trace shows the PSD as the surface roughens predominantly at the linear characteristic length κ = 1/2 indicated by the vertical line. Small nonlinear effects can already be seen in this last trace in the form of small peaks at harmonics of the fundamental characteristic length.  9  Chapter 2. Numerical Results  τ = 0.0 5×10−7  χy  0  100  0  200  -5  τ = 0.7  χy  0  3×10−7 2 1 0 -1 -2 -3  100  200  τ = 40.0  0  2×10−3  χy  1 100  0 -1  200  -2 0  100 χx  200  Figure 2.3: Top view of the simulated surface at key times during the pattern formation process. The white noise initial condition is shown at τ = 0. The surface at τ = 0.7 has undergone rapid smoothing. At τ = 40 a clear pattern has formed with a characteristic length.  10  Chapter 2. Numerical Results  τ = 0.0 6 4 2 0 -2 -4 -6 −7 6 ×10  τ = 0.7  4  η  2 0 -2 -4 -6  τ = 40.0  3 2 1 0 -1 -2 -3 0  100 χy  200  Figure 2.4: Cross-sections along the line χx = 0 of the solutions shown in Fig. 2.3. It is noteworthy that in all three cases, the surface retains updown symmetry, that is the same solutions with η → −η would be equally plausible. 11  Chapter 2. Numerical Results It can be seen from Fig. 2.1(b) that prior to the exponential growth, the RMS height undergoes an initial rapid decrease. Since the initial condition contains fluctuations at all length scales as seen by the τ = 0 trace in Fig. 2.2, a large fraction of the surface fluctuations are at length scales κ > 1 which are exponentially damped according to Eq. 2.10. The initial decrease thus corresponds to this initial decay which is far more rapid than the growth for longer length scales κ < 1. This is confirmed by the fact that spatial frequencies κ > 1 have been drastically damped in the τ = 0.7 trace of Fig. 2.2 while frequencies κ < 1 are essentially unchanged from their τ = 0 state. Likewise, the image of the surface at τ = 0.7 in Fig. 2.3 and its crosssection in Fig. 2.4 are significantly smoother than the surface at τ = 0, but still appear relatively random at low frequencies. Once these short wavelength fluctuations have subsided, the exponential growth at the longer wavelengths begins to dominate and the RMS height increases exponentially. The τ = 40 trace in Fig. 2.2 demonstrates that this growth is predominantly at the linear characteristic length κ = 1/2 where κ2 − κ4 achieves its maximum. This characteristic length can clearly be identified in the surface image for τ = 40 in Fig. 2.3 though it may not be as apparent from the corresponding cross-section in Fig. 2.4. According to Eq. 2.10, the time constant of growth at this dominant length scale should be τc = 4. It is not surprising then that the fit to the exponential growth regime of the RMS height in Fig. 2.1 has a time constant of τc ≈ 4.1. This is naturally slightly larger than for the characteristic length alone. The slower growth of other length scales κ < 1 also contribute to the growth in the RMS height. Nevertheless this growth is strongly dominated by the characteristic length.  2.5  The Early Nonlinear Regime  The previous section described the behaviour of Eq. 2.7 in the linear regime. In this regime, the time evolution of the solution is easily understood and predicted through analysis of the linear terms alone. It is possible to neglect the nonlinear terms due to the condition |∇η| ≪ 1 being met. Clearly however, this inequality will quickly be violated given the exponential growth in the surface amplitude seen above. A clear sign of the transition from the linear to the nonlinear regime can be seen in Fig. 2.1 where the RMS height begins to level off to η 2 ∼ 1 at τ ∼ 70. This saturation of the height cannot be understood by considering only the linear terms of Eq. 2.7. Moreover, it is at this point that the RMS 12  Chapter 2. Numerical Results heights for solutions with different values of θ begin to diverge indicating that the nonlinear terms, and more specifically, their ratio, have become important. Solutions themselves also look qualitatively different once this nonlinear saturation has taken place. Fig. 2.5 shows cross sections of solutions at varying θ at τ = 70. Unlike the cross sections in Fig. 2.4, these no longer have up-down symmetry. The tops of peaks have become notably sharper, while valleys are rounded. Again, it appears that unlike in the linear regime, the value of θ has clear effects on the solution. Fig. 2.5 shows that the solution for θ = 0.5π has deeper and more defined parabolic valleys, whereas the solution for θ = π is more disordered with more low frequency fluctuations. These observations are even more evident from the PSD shown in Fig. 2.6 where the low frequences are significantly damped for θ = π/2, and the peak at the characteristic length is much more distinct. The top view of the solutions Fig. 2.7 additionally show an increase hexagonal ordering as θ decreases from 1.00π.  2.6  The Long-time Nonlinear Regime  The persistence of the height saturation hinted at in Fig. 2.1 is demonstrated in Fig. 2.8. Here, each curve represents the ensemble average of the RMS height over 4 separate solutions with different initial conditions but the same value of θ. It is likely that a larger ensemble average would result in a smoother time dependence. Nevertheless, it is clear that the RMS height indeed levels off and remains at a relatively constant value for the length of time simulated for θ = π, θ = 0.75π and θ = 0.6π. The RMS height for θ = 0.55π however doesn’t appear to have reached its full saturation value within this time. This is complicated by the fact that the adaptive timestep algorithm described in Section A.3 finds that a much smaller timestep is required as θ → π/2. The reasons for this are discussed toward the end of Section A.2 as are some ideas for mitigating this problem. The present work however sidesteps this issue by narrowing the range investigated to θ ∈ [0.6π, π]. This will prove to be more than sufficient for comparison to experimental data for snow discussed in Chapter 3. Within this narrowed range, solutions for all values of θ appear to saturate to an approximately constant RMS height. As seen in Fig. 2.8 the solutions for θ near to π find this saturation height very quickly, while θ = 0.6π doesn’t fully reach its saturation height until τ ∼ 2000. Since there is sig13  Chapter 2. Numerical Results  θ = 0.50π 4 2 0 -2 -4 θ = 0.75π 4  η  2 0 -2 -4 θ = 1.00π 4 2 0 -2 -4 0  100 χy  200  Figure 2.5: Cross-section of the surface at τ = 70 along the line χx = 0 just as the RMS height begins to saturate. The solutions no longer have up-down symmetry and clear distinction can be made between the solutions of different θ. 14  Chapter 2. Numerical Results  PSDr  100  θ = 0.50π θ = 0.75π θ = 1.00π  10−5  10−10 0.5  1  1.5  2 κ  2.5  3  3.5  Figure 2.6: Early nonlinear saturation PSD. Low frequency fluctuations are reduced as θ decreases while the exponential fall off of higher frequencies becomes slower. Harmonics are present in the θ = π/2 solution. nificant interest in the intrinsic properties of the equation independent of particular initial conditions, care is taken to ensure any time averages or results concerning the nonlinear dynamics of Eq. 2.7 are taken only from the fully saturated regime. Top views for solutions in this regime are shown in Fig. 2.9 for τ = 4000 while cross-sections at this time are shown in Fig. 2.10. Here the differences across different values of θ are even clearer. Additionally, it is quite clear that the characteristic length gets longer as θ increases. Additionally, significant low frequency fluctuations now appear for all values of θ shown.  2.7  Power Spectral Density in the Long-Time Nonlinear Regime  The difference in length scale and the presence of significant low frequency fluctuations motivates a closer study of the power spectral densities for the long-time nonlinear regime. Since it appears that solutions do not vary qualitatively in this regime, the time average of the radial power spectral density is taken in addition to the small ensemble average over 4 solutions 15  Chapter 2. Numerical Results  θ = 0.50π 0  χy  5  100  0  -5  200  θ = 0.75π  0  χy  5  100  0  -5  200  θ = 1.00π  0  χy  5  100  200  0  -5 0  100 χx  200  Figure 2.7: Top view of the surface at τ = 70 just as the RMS height begins to saturate. Differences in the degree of hexagonal ordering and in the amount of low frequency noise for different θ can be distinguished.  16  Chapter 2. Numerical Results  5 θ θ θ θ  4  = 0.55π = 0.60π = 0.80π = 1.00π  h  2  3 2 1 0 0  500  1000  1500  2000 t  2500  3000  3500  4000  Figure 2.8: Long timescale behaviour of the RMS height for solutions on a domain of length L = 64π demonstrating nonlinear saturation near η 2 = O(1). The values plotted are the ensemble average over four simulations with different random seeds for their initial conditions.  17  Chapter 2. Numerical Results  θ = 0.60π 0  y  5 100  0  -5 200 θ = 0.80π 0  y  5 100  0  -5 200 θ = 1.00π 0  y  5 100  0  -5 200  0  100 x  200  Figure 2.9: Top view of solutions in the long-time nonlinear regime (τ = 4000) well after the solutions shown have fully saturated. All solutions have significant low frequency fluctuations, but the characteristic length increases as θ decreases in the range shown. 18  Chapter 2. Numerical Results  θ = 0.60π 6 4 2 0 -2 θ = 0.80π 6  h  4 2 0 -2 θ = 1.00π 6 4 2 0 -2 0  100 x  200  Figure 2.10: Cross-sections of the solutions in the long-time nonlinear regime (τ = 4000) along the line χx = 0.  19  Chapter 2. Numerical Results with different initial conditions. This has the desirable effect of reducing stochastic fluctuations. The time averages are taken after the solutions have sufficiently annealed. For the range θ ∈ [0.6π, π], it is found adequate to take the time average from τ ∈ [2000, 4000]. By this time, the amplitude of all the individual Fourier modes of the solutions in a domain of size L = 64π as well as L = 32π appear to have saturated such that they fluctuate about a constant mean value. Together with the ensemble average, this gives an average over a total of total of 8000 dimensionless time units for each pair of θ and L. The ensemble and time averaged radial power spectral density is shown for three values of θ in Fig. 2.11 on loglog, linear and semilog axes. These show clear peaks near, but at slightly longer wavelength than the fastest growing mode in the linear regime κ = 1/2. Furthermore, the position of the peak shifts toward lower frequency as the θ decreases. Below this peak, the PSD of all θ have power law behaviour very close to κ−2 as can be seen by comparison to the dashed line in the loglog plot Fig. 2.11. The slight divergence of the θ = 0.6π trace from this trend is likely due random fluctuations. A longer time average, or more practically, a large ensemble average is needed to comment more conclusively on this feature. Finally, at high frequencies, the power spectral density is found to fall off exponentially for all θ simulated, which is consistent with the solutions being analytic (i.e., all derivatives exist). Interestingly however, the rate of exponential falloff does decrease as θ decreases. This is likely a result of the fact that the pointed ridges of the solutions are sharper at lower θ. The observation that the position of the peak in the power spectral density shifts with θ leads naturally to the question of just how the two are related. To answer this, a method is needed for determining the position of the peak. The most straight forward method involves finding the local maximum and minimum sample points of the radial PSD. These are points that are higher or lower respectively than their two adjacent points. The local minimum with the lowest κ is then defined as κmin . The local maximum κmax with the lowest κ such that κmin < κmax is taken to correspond to the characteristic length. Given that the largest domain size simulated is L = 64π, the resolution at which the PSD is sampled is ∆κ = 1/32. Merely taking κc = κmax would then yield an uncertainty of 1/64. This can be improved upon by interpolating the PSD about κmax in the peak region. For this purpose, a Newton divided difference interpolation polynomial [5] is used to interpolate 20  Chapter 2. Numerical Results  102 100 κ−2  10−2 10−4  100  10−1  3 θ = 0.60π θ = 0.80π θ = 1.00π  2.5 PSDr  2 1/2  1.5 1 0.5 0 0.5  0  1  1.5  2  102 100 10−2 10−4 10−6 10−8 0  0.5  1  1.5  2 κ  2.5  3  3.5  Figure 2.11: Time and ensemble averaged power spectral density of the solution for three values of θ in the fully saturated nonlinear regime. The same data is plotted in all three plots above on different axis scales to highlight different features. The loglog plot highlights the low frequency κ−2 behaviour, the linear plot highlights the peak at the characteristic length compared to κ = 1/2 (dashed line). Finally, the semi-log plot shows the exponential falloff of the PSD at high frequencies. 21  Chapter 2. Numerical Results  1 κc  PSDr  0.8  κmax  0.6 0.4  κmin  0.2 0 0  0.2  0.4  0.6  0.8  1  κ Figure 2.12: Demonstration of the determination of the characteristic by finding the peak in the polynomial interpolated PSD. The location of the peak is indicated at κc = 0.69. The example shown is for θ = 0.9π and L = 64π. the PSD in the region κ ∈ [κmax − 1/8, κmax + 1/8]. For the L = 64π solutions, this range spans 9 sample points while for L = 32π it spans 5. The maximum of the interpolation polynomial is found by solving for the zeros of its derivative. It is expected that κc will be near the centre of the domain of the interpolation polynomial where it is the most accurate. Having determined the position of the power spectral peak κc , the characteristic length is then ℓc = 2π/κc . An other method of finding the characteristic length is to examine the cross-sectional PSDs. These are derived by averaging the one dimensional PSDs of all vertical and horizontal cross-sections of the solution. ˆ y )|2 k ˆ x )|2 k + |h(k L |h(k x y (2.12) 2 2πN 2 Like the radial PSDs, the time and ensemble average is then taken as well. The resulting PSDs are shown in Fig. 2.13. It may be noted that unlike the κ−2 dependence of the radial PSDs at low frequencies, the cross-sectional PSDs exhibit a κ−1 dependence. The reason for this may be easily understood by first realising the 1D PSD in P SDcs =  22  Chapter 2. Numerical Results  104  102  κ−1  100  100  10−1  PSDcs  1500 θ = 0.60π θ = 0.80π θ = 1.00π  1000  1/2  500  0 0  0.2  0.4  0.6  0.8  1  104 102 100 10−2 10−4 10−6 10−8  0  0.5  1  1.5  2 κ  2.5  3  3.5  4  Figure 2.13: Time and ensemble averaged cross-sectional power spectral density of the solution for three values of θ in the fully saturated nonlinear regime. This plot may be compared with the radial PSDs in Fig. 2.11  23  Chapter 2. Numerical Results the y direction averaged along the x direction is equivalent to the average of the 2D PSD along the x direction to within a constant. If the 2D PSD expressed verbosely as F[F[f (χx , χy )](κy )](κx ) (to within a constant) goes −1 for low frequencies as shown in Fig. 2.11, then as |κ|−2 = κ2x + κ2y (F[f (χx , χy )](κy ))2 dχx = C = C′  (F[F[f (χx , χy )](κy )](κx ))2 dκx (κ2x + κ2y )−1 dκx = C ′  π κy  (2.13)  Other than the different power law, the cross-sectional PSDs differ from the radial PSDs in the fact that the peaks corresponding to the characteristic length are much more difficult to distinguish. They appear as more of a shoulder than a proper peak. They can however be made into peaks by multiplying the cross-sectional PSDs by a factor of κ which removes the κ−1 dependence as shown in Fig. 2.14. A second measurement of the characteristic length can then be determined from the cross-sectional PSDs using the same polynomial interpolation method as above. The procedure of removing the background power law from the crosssectional PSDs may be applied also to the radial PSDs. In this case, we multiply them by κ2 as shown in Fig. 2.15. Though the peaks were distinguishable before, this makes them even clearer. In fact, they are now the absolute maxima of the PSDs which makes locating them algorithmically much easier. A third value for the characteristic length is computed from these modified radial PSDs using the now familiar polynomial interpolation method. Having computed the characteristic length using the three methods outlined above, the relation between these length scales and the ratio of the nonlinear terms tan θ is shown in Fig. 2.16. It can be seen found that the characteristic length has an approximately linear relationship of slope ∼ −0.13 with the ratio between the nonlinear coefficients tan θ. The characteristic length measured from the raw radial PSD is slightly longer than that measured from the modified cross-sectional PSD which in turn, is slightly longer than the value measured from the modified radial PSD. The direct radial method gives a characteristic length that is always greater than the linear value and has the noteworthy property that unlike the other two measurements, it achieves its minimum not at the endpoint tan θ = tan π = 0, but at tan θ = tan 0.95π ≈ −0.16. The other two measurements give a characteristic length that is slightly smaller than the linear characteristic length when θ ≥ 0.95π. The three measurements are in 24  Chapter 2. Numerical Results  103 102 101 100  10−1  100  400 θ = 0.60π θ = 0.80π θ = 1.00π  κ PSDcs  300  1/2  200 100 0 0  0.2  0.4  0.6  0.8  1  102 100 10−2 10−4 10−6  0  0.5  1  1.5  2 κ  2.5  3  3.5  4  Figure 2.14: Cross sectional PSDs multiplied by κ. This flattens out the κ−1 dependence at low frequencies and results in a much clearer peak at the characteristic length.  25  Chapter 2. Numerical Results  100 10−1 10−2 10−3  10−1 1/2  0.6 κ2 PSDr  100 θ = 0.60π θ = 0.80π θ = 1.00π  0.4 0.2 0 0  0.5  1.5  1  2  100 10−2 10−4 10−6 10−8 0  0.5  1  1.5  2 κ  2.5  3  3.5  Figure 2.15: Radial PSDs as shown in Fig. 2.11 multiplied by κ2 .Though the peaks were mostly distinguishable before, they become even more prominent under this modification.  26  Chapter 2. Numerical Results 1.5  √ lc /(2π 2)  1.4 1.3 1.2 1.1 PSDr κ PSDcs κ2 PSDr  1 0.9 -3.5  -3  -2.5  -1.5  -2  -1  -0.5  0  tan θ  Figure 2.16: Characteristic length as a function of θ measured from the peaks in three different PDSs. The characteric length is expressed in mul√ tiples of the linear characteristic length 2π 2. Error bars are calculated by taking the difference in the values obtained from domains of length L = 32π and L = 64π. best agreement (∼ 4%) in the range tan(θ) ∈ [−1.5, −0.5] (θ ∈ [0.7π, 0.85π]). This is likely due to the fact that the PSD peak is narrowest and most clearly defined in this range as seen by the θ = 0.8π trace in Figs. 2.11, 2.14, 2.15.  2.8  Saturation Height  It was noted in passing in Section 2.5 that the RMS height of the solution varied for different values of θ. As it was with the characteric length ℓc , it is natural to ask about the relationship between the saturation height and θ. With this motivation, the fully saturated RMS height η 2 is time averaged over τ ∈ [2000, 4000] and ensemble averaged over the 4 solutions available for each L and θ. The result for two domains sizes L is plotted as a function of θ in Fig. 2.17. Unfortunately, this figure demonstrates that the saturation value of this height measurement is not independent of the domain length L. The reason for this is easily understood given the PSDs plotted in Fig. 2.11. The loglog scale in particular shows that a non27  Chapter 2. Numerical Results  2.6 L = 32π L = 64π  2.4  2  η  2  2.2  1.8 1.6 1.4 0.6  0.65  0.7  0.75  0.8 θ/π  0.85  0.9  0.95  1  Figure 2.17: Time average of the RMS height after full saturation as a function of θ for L = 64π and L = 32π. The dependence on domain size demonstrates that the RMS height is not a suitable parameter to compare with experimental observation. negligible portion of the surface fluctuations are at length scales on the order of the domain size. If the power law behaviour of the PSD persists to larger domains, then the RMS height will steadily continue to grow with L no matter how large the domain is made. What is needed is a more local measure of surface height and furthermore, one that is more readily compared to experimental data where the full surface topography is not necessarily known. To this end, the crosssectional local peak-to-peak amplitude η lptp1 is defined with a gross abuse of notation. This is calculated by taking a cross-section of the surface and locating the local minima and maxima. The cross-sectional local peak-topeak amplitude is then defined as the mean difference in height between the minima and their two adjacent maxima. Naturally, given the full two dimensional grid of a numerical solution, this can be computed for every possible cross section in both the x and y direction and the average value taken to represent the surface as a whole. A similar, more intuitive measure of two dimensional suncup depth may also be defined. Rather than taking the mean difference between minima  28  Chapter 2. Numerical Results and the adjacent maxima in cross-sections, the 2-D local peak-to-peak height η lptp2 is defined as the mean difference between 2-D minima and their nearest local maximum. This value will in general be larger than η lptp1 because a minimum in a given cross-section will not generally be the minimum of the full 2-D suncup. Though this value is easier to estimate by visual inspection, it is more difficult to experimentally measure in a precise manner. The time and ensemble average of both of the above described local peak-to-peak height measurements averaged over the same time domain as the RMS height above is shown in Fig. 2.18. Both values are found to be independent of domain size as shown by the error bars which indicate the discrepancy between the values measured from L = 64π and L = 32π solutions. Interestingly, the 2-D value is consistently found to be very close √ to a factor of 2 larger than the cross-sectional value. The ensemble averaged time dependence of both local peak-to-peak height values is shown in Fig. 2.19 demonstrating that both are significantly more robust against fluctuations than the ensemble averaged RMS height in Fig. 2.8. This offers an additional advantage when measuring this value given limited experimental data.  2.9  Effect of the |∇η|2 Term on the Mean  Being the square of a real-valued number, the term |∇η|2 in Eq. 2.7 is always positive. This implies that for a surface with nonzero slope anywhere, its mean is positive definite. It is the only term in Eq. 2.7 to have this property. All others have average to zero over the whole surface. Recalling that in the range θ ∈ [π/2, π], cos θ is always negative, increasing the strength of this term by increasing θ towards π will increase the rate at which the mean height recedes according to ∂τ η = cos θ |∇η|2  (2.14)  It should be noted that this recession is in addition to the −c0 t0 /h0 recession that has been subtracted out of Eq. 2.5 by shifting the coordinates of the surface into the moving frame. For this reason, the recession associated with Eq. 2.14 it is sometimes referred to as “excess velocity”[2]. Given a value for θ, determination of the excess velocity of the surface also requires knowledge of the value |∇η|2 . This value depends on the particular shape of the surface and, not surprisingly varies with θ much like the characteristic length and the characteristic height discussed above. The time and ensemble averaged θ dependence of |∇η|2 is shown in Fig 2.20. 29  Chapter 2. Numerical Results  5 lptp2 lptp1  4.5  η  4 3.5 3 2.5 2 0.6  0.65  0.7  0.75  0.8 θ/π  0.85  0.9  0.95  1  Figure 2.18: Local peak-to-peak height time averaged over the domain τ ∈ [2000, 4000] as a function of θ. Both the cross-sectional and 2-D values are shown. The error bars represent the difference between the measurements from the L = 32π and L = 64π domains  30  Chapter 2. Numerical Results  3.5 3  lptp1  2  η  2.5  1.5 1 0.5 0 0  500  1000  1500  2000 τ  2500  3000  3500  4000  (a) Time dependence of the cross-sectional local peak-to-peak height  5  η  lptp2  4 3 2 θ = 0.60π θ = 0.80π θ = 1.00π  1 0 0  500  1000  1500  2000 τ  2500  3000  3500  4000  (b) Time dependence of the 2-D local peak-to-peak height  Figure 2.19: Time evolution of both the cross-sectional η η lptp2 local peak-to-peak height for three values of θ.  lptp1  and 2-D  31  Chapter 2. Numerical Results The relation has a roughly parabolic shape with a minimum beteween θ = 0.80 and θ = 0.75. Discrepancy between |∇η|2 computed for L = 32π and L = 64π is essentially nonexistent as indicated by the error bars. The values for |∇η|2 for the range of θ shown have a relative standard deviation of 34%. It is worth noting that if the height and length are measured by their characteristic values as determined above for each θ, then the rescaled value of |∇η|2 has a relative standard deviation of only 7% as shown in Fig. 2.20(b).  2.10  Diffusion of the Simulated Suncups  Having measured the height and length scale of the solution in the nonlinear regime, the next natural measurement is the timescale. This is perhaps both the most interesting and difficult measurement to make, particularly in a way that might be easily measured experimentally. One possibility would be to measure the correlation time either of real space grid points or Fourier modes. Such a measurement may be worthwhile for a more mathematical investigation of the properties of Eq. 2.7, however it is rather impractical to measure experimentally. The strategy adopted is to measure the dynamic feature which is most visually striking when viewing the solution in motion. This feature is arguably the movement of individual suncups around the surface. This may be quantified by tracking the position of local minima as the solution evolves. For these purposes, local minima are defined as grid points in the solution whose neighbours are all greater in height. This is verified to accurately and comprehensively identify the locations of individual suncups as shown in Fig. 2.21. As a further, more quantitative sanity check, the mean density of the minima per characteristic length squared in the range τ ∈ [2000, 4000] is plotted in Fig. 2.22. As might be expected, the density is within 20% of one minima per characteristic length squared. Lower θ do however appear to have a slightly lower relative density of minima. The individual minima are linked between snapshots of the surface at time intervals of ∆τ = 1 by associating the nearest points with each other. Minima are regularly observed to be created and destroyed. This occurs when there are no available maxima or minima nearby to link a given point at either the previous or next time snapshot respectively. Visually, this typically appears as the joining or splitting of existing minima. The mean creation/destruction rate relative to the total number of minima is shown as a function of θ in Fig. 2.23. Both the creation and destruction rate 32  replacemen Chapter 2. Numerical Results  1.4 1.2  |∇η|2  1 0.8 0.6 0.4 0.2 0 0.6  0.65  0.7  0.75  0.8 θ/π  0.85  0.9  0.95  1  0.85  0.9  0.95  1  (a) unscaled  14  |∇η|2 (lc / η  2 lptp1 )  12 10 8 6 4 2 0 0.6  0.65  0.7  0.75  0.8 θ/π  (b) scaled by the characteristic height and length  Figure 2.20: |∇η|2 as a function of θ. Both (a) the raw value as well as (b) the value measured in units of characteristic lengths ℓc and local peak-to-peak height η lptp1 are shown. Time averages are taken on the domain τ ∈ [2000, 4000] and the ensemble average is taken over 4 solutions with differing random seed for the initial condition. Error bars show the miniscule discrepancy between L = 32π and L = 64π domains. 33  Chapter 2. Numerical Results  0  100  200 0  100  200  Figure 2.21: Trajectories of the minima of the numerically simulated suncups. The trajectories shown represent the motion of the minima of the presently shown solution at τ = 4000 for which θ = 0.80π. The trajectories shown are arbitrarily limited to the past 250 dimensionless time units for optimal esthetic effect. Naturally, this is only an upper bound on the age of displayed trajectories due to the creation of new minima during this time.  34  Chapter 2. Numerical Results  1.1 1.05  #/ℓ2c  1 0.95 0.9 0.85 0.8 0.75 0.6  0.65  0.7  0.75  0.8 θ/π  0.85  0.9  0.95  1  Figure 2.22: Density of suncups per characteristic length as a function of θ. The characteristic length measured from the modified radial power spectral density κ2 PSDr is used. As expected, this is close to unity for all values of θ. Error bars indicate the discrepancy between L = 64π and L = 32π.  35  relative creation/destruction rate  Chapter 2. Numerical Results  10−1  10−2  10−3 0.6  0.65  0.7  0.75  0.8 θ/π  0.85  0.9  0.95  1  Figure 2.23: Time averaged Creation/Destruction rate of suncups in the time interval τ ∈ [2000, 4000] as a function of θ. The relation is found to be surprisingly close to exponential in the range θ ∈ [0.6π, π]. Error bars indicate the discrepancy between the values from L = 32π and L = 64π. are essentially identical meaning that any coarsening that took place in the early nonlinear regime has ceased before the sampled time. This is a further sanity check that the solution is in a statistically steady state since it indicates that the number of minima is not changing on average. The relation between θ and the relative creation/destruction rate is surprisingly close to exponential in the range of θ simulated with θ = π having the greatest relative creation/destruction rate. The time continuity of the minima that results from this linking permits the tracking of trajectories as shown in Fig. 2.21. These trajectories may then be used to calculate the displacement of individual minima enumerated by i as a function of time |χi (τj + ∆τ ) − χi (τj )| where τj is some reference time. The root mean square of this value over all minima is computed to give the time dependent diffusion length ℓd (∆τ ) (whose nomenclature will soon become obvious). ℓd (∆τ ) =  |χi (τj + ∆τ ) − χi (τj )|2  i  (2.15)  The average can however be improved by averaging over all possible 36  Chapter 2. Numerical Results reference times τj , ℓd (∆τ ) =  |χi (τj + ∆τ ) − χi (τj )|2  j i  (2.16)  Note that the average over each reference time j is done first for each minimum i. This is because the allowable set of τj for the ith minimum is restricted to the times for which that minimum exists (i.e., between its creation and destruction). Ensembles of 4 solutions for each value of θ ∈ [0.6π, π] with L = 64π and L = 32π are used to assemble a set of trajectories χi (τ ) restricted to τ ∈ [2000, 4000] from which an ℓd (∆τ ) is computed using Eq. 2.16. The ℓd (∆τ ) resulting from the L = 64π solutions at several values of θ are plotted in Fig. 2.24. A striking feature is that for all values of θ, the diffusion length follows a power law ℓd (∆t) = (D∆t)b  (2.17)  with approximately the same exponent b ∼ 1/2 which is the standard result for diffusion of a random walker. The value of D on the other hand, which is equivalent to the diffusion coefficient, increases with θ as seen by the increasing vertical offset of curves on the log scale y-axis. To quantify this relationship, Eq. 2.17 is fitted to each ℓd (∆τ ) to determine D and b. Since there appears to be a slight distortion of the displacement curves at the scale of the solution grid spacing, the fit only takes into account data points with ℓd larger than this size. Each data point is weighted in the fit according to the square root of the number of trajectories used in the average at that ∆τ . This number naturally decreases at longer time since there are fewer trajectories that live long enough to contribute to the ℓd at larger ∆τ . The mean value of the exponent is found to be b = 0.58 with a standard deviation of δb = 0.04. The diffusion coefficient D as a function of θ is plotted in Fig. 2.25. The monotonic increase of D with θ provides a clear correspondence between the value of θ from an experimentally measurable quantity. This will be used to determine θ from experimental data of snow surface dynamics in Chapter 3.  37  Chapter 2. Numerical Results  ∆τ 1/2  ℓd  101  100  100  101  102  103  ∆τ (a) Unscaled (dotted line indicates mesh size)  ∆τ 1/2  ℓd (ch.len.)  100  θ θ θ θ θ  10−1  100  101  102  = 1.0π = 0.9π = 0.8π = 0.7π = 0.6π 103  ∆τ (b) Scaled to characteristic length  Figure 2.24: RMS displacement of the numerical suncup minima as a function of time for different values of θ. The two plots above show this relation in both (a) natural units and (b) with each θ trace scaled to its characteristic length as determined from κ2 PSDr . The grid spacing size is marked by the dashed line in (a) which suggests the reason for the slight bowing of trace near this point. 38  Chapter 2. Numerical Results  D  100  10−1  0.6  0.65  0.7  0.75  0.8 θ/π  0.85  0.9  0.95  1  0.85  0.9  0.95  1  D (ch.len.) b  1  (a) Unscaled  10−2  10−3 0.6  0.65  0.7  0.75  0.8 θ/π  (b) Scaled to characteristic length  Figure 2.25: Diffusion coefficient for suncup minima as a function of θ. The two plots above show that both the (a) unscaled and (b) scaled relations are quite close to exponential spanning slightly more than a decade in the range θ ∈ [0.6π, π]. Error bars are computed from the difference between L = 64π and L = 32π. 39  Chapter 3  Field Observations So that the results of the previous section may be compared with physical snow surfaces, two field observation sites are established in Southern British Columbia as shown in Fig. 3.1. The first site is in a valley in the Selkirk Mountains at an elevation of 2240m and coordinates 51◦ 3′ 52′′ − 117◦ 32′ 0′′ . Two areas of snow 2−3.5 m in diameter where suncups had formed naturally are flattened and the reformation of suncups is observed during the period July 24-30, 2007. Cross-sectional photographs such as the one shown in Fig. 3.2 are taken of trenches cut in these flattened areas at different stages of suncup development. An additional trench is also cut out of an undisturbed area of snow. The edges of these trenches are digitised and transformed into true cross-sectional traces using the methods described in Section B.7. This allows for the measurement of the roughening rate as well as the saturation height and characteristic length of suncups. The second site established is on top of Whistler Mountain located at ◦ 50 3′ 36” − 122◦ 57′ 38” and 2160m above sea level. An overview of the site is shown in Fig. 3.3. Here, a webcam is set up which makes use of a preexisting wireless infrastructure servicing the Whistler Blackcomb ski resort. The webcam is aimed at a snowbank several metres deep where suncups had developed. Images are downloaded every minute during the period July 21September 10, 2007 and night time images are automatically removed from the dataset. One image from each day is semi-manually selected for optimal contrast on the snow surface. The images selected in this way range in time from 1:30pm-9:30pm due to variable lighting conditions. Selected images are then transformed to “snow coordinates” using the methods described in Appendix B so that the snow surface appears as though viewed from directly above. These data are used to measure the dynamics of suncups in the nonlinear regime. In some cases, the Whistler data also serves to corroborate some of the measurements made at the Selkirk observation site.  40  Chapter 3. Field Observations  Figure 3.1: Map of the snow observation sites. Whistler Mountain is the leftmost Marker. The Selkirk site is marked on the right.  3.1  Measuring the Local Peak-to-Peak Height  The edges of the trenches in the photographs taken at the Selkirk site are digitised by hand and transformed so that they trace out the intersection of a vertical plane with the horizontal snow surface. Because of this, the y-axis of the resulting traces shown in Fig. 3.4 may be interpreted as the vertical deviation of the snow surface, while the x-axis may be interpreted as the lateral position along a line on the snow surface. Mathematically, the trace can be interpreted as the position dependent height h(x, y0 ) where y0 is a constant. Some of the noise associated with the digitisation process is removed from these cross-sections. To achieve this without fundamentally altering the topography, a Savitzky-Golay[19] smoothing algorithm of order 8 with a span of approximately 300 pixels or ∼ 30cm is applied to the data. This smoothing method has the advantage of preserving the height of minima and maxima which are particularly important in measuring the height scale for comparison to the numerical simulations. Local minima and maxima are located within the smoothed cross-sections as shown in Fig. 3.4. Since some of the lower frequency noise in the digitised cross-sections cannot be removed by the relatively conservative SavitzkyGolay smoothing, the local minima and maxima are located in the con41  Chapter 3. Field Observations  Figure 3.2: One of the flattened patches at the Selkirk observation site. This area was flattened two days before the photo was taken. The edges of the trenches are digitised and transformed into vertical height profiles as described in Appendix B.  42  Chapter 3. Field Observations  Figure 3.3: Overview of the Whistler site where dynamic observations of the snow surface in the foreground were made. The camera was mounted on the pole to the right. Whistler village can be seen in the background in the upper right.  43  Chapter 3. Field Observations volution of the cross-sections with a Gaussian of width ∼ 2 − 5cm. This reduces false positive identification of the minima and maxima. Once the approximate locations of the extrema in the convolved cross-sections are determined, their location in the cross-sections with only Savitzky-Golay smoothing are refined by searching for the largest or smallest nearby points to within the width of the Gaussian. The resulting maxima and minima are indicated by the dots in Fig. 3.4. By taking the mean difference in height between minima and adjacent maxima, an experimental measurement of the cross-sectional local peak-to-peak height h lptp1 as defined for the numerical simulations in Section 2.8 is made. Since, this is the only measure of topography height significantly discussed in this chapter, the lptp1 subscript will be omitted in further discussion so that unless otherwise specified · indicates the cross-sectional local peak-to-peak height.  3.2  Roughening: Time and Height scale  The three cross-sections shown in Fig. 3.4 each represent one of the stages of suncup development for which cross-sections were taken, namely 2 days, 5 days and the long-time limit (i.e., the unmodified surface). This last crosssection is taken to be equivalent to a surface 40 days after flattening. This number is arbitrary except that it is significantly longer than the 5 days required for the suncups to achieve their saturation height, but not so long as to exceed the amount of time that the suncups are likely to have been present. Outside of these broad requirements, its placement is ensured not to significantly affect the subsequent results. Since multiple cross-sections are taken at each of the above mentioned times, they may be averaged together at each time. An error estimate on this average may be computed in the form of the standard error δ h(t) =  h(t)  2  − N  h(t)  2  (3.1)  where N is the number of data points included in the average. For clarity, it should be stressed that · applies only to the spatial, and not the time component of h, so that h(t) is a time dependent function. Similarly · applies only to the set of h available at a given time so that again, h(t) is also a time dependent function. In addition to these three data points for which the cross-sectional local peak-to-peak height is directly measured, more approximate “by eye” estimates of the surface height are also made after 0, 1, 3, and 4 days. 44  Chapter 3. Field Observations t = 2 days 4 2 0 -2 -4 -6 t = 5 days 4 h (cm)  2 0 -2 -4 -6 t ≫ 5 days 4 2 0 -2 -4 -6 0  50  100  150  200 x (cm)  250  300  350  Figure 3.4: Experimental cross-sectional surface profiles from the Selkirk observation site. Dots indicate maxima and minima. 45  Chapter 3. Field Observations  h (cm)  101  100  0  5  10  15  20 t (days)  25  30  35  40  Figure 3.5: Cross-sectional local peak-to-peak height from numerical simulations (solid line) with θ = 0.8π scaled to fit experimental data (errorbars). The error on these is assumed to be 20%. These height estimates however more closely resemble the 2-D local peak-to-peak height measurements. In order to meaningfully include these measurements with those from the crosssections, they are scaled such that the height on the day 4 is equal to the local peak-to-peak height measured directly from the cross-section on day 5. The net result is a time series of 7 data points in time h(ti ) ≡ hi with error estimates δ h(ti ) ≡ δ hi corresponding to the roughening and saturation of the snow surface height as shown in Fig. 3.5. This shows that suncups grow in amplitude by a little less than an order of magnitude in 5 days. This is consistent with less precise observations at the Whistler site where a flattened area was observed to increase in roughness by approximately half an order of magnitude in ∼ 3 days. The Selkirk data shown in Fig. 3.5 are compared with the cross-sectional local peak-to-peak height from the numerical solutions of Eq. 2.7. Solutions for 9 values of θ ∈ [0.6π, π] spaced at equal intervals of 0.05π are considered. The time series of the local peak-to-peak height ηθ (τ ) for these numerical solutions are scaled in height and time by factors of h0 and t0 respectively. These factors are the same as those defined in Eq. 2.1 and Eq. 2.3 used to  46  Chapter 3. Field Observations arrive at the dimensionless formulation Eq. 2.7. Together with a time offset τ0 , these scales are adjusted so as to minimise the sum N i=1  hi − h0 ηθ (ti /t0 + τ0 ) δ hi  2  (3.2)  Cubic interpolation is used to evaluate ηθ at times between regularly sampled intervals of 1 dimensionless time unit and Nelder-Mead simplex optimisation[10] is used to minimise Eq. 3.2. This procedure is applied for each θ to obtain h0 (θ) and t0 (θ). Values of τ0 (θ) are also obtained in this process, but are not of interest for further calculations. They merely indicate the time in the numerical simulation defined by h0 η(τ0 ) ≈ h(0days). In order to obtain error estimates on the values of h0 (θ) and t0 (θ), the fit is repeatedly performed for each θ with the experimental points hi perturbed by a set of random values ρji δ hi where ρji is a set of normally distributed random variables with unit standard deviation and 1 ≤ j ≤ 10, 000. A set of 10, 000 hj0 (θ), tj0 (θ) and τ0j (θ) are then implicitly defined my minimising the 10, 000 sums N i=1    hi + ρji δ hi    − hj0 ηθ (ti /tj0 + τ0j ) δ hi  2   (3.3)  The resulting sets {hj0 (θ)} {tj0 (θ)} for each θ can be used to estimate the errors δh0 (θ) and δt0 (θ) by taking their standard deviation. The scales h0 and t0 along with their uncertainties are shown as a function of θ in Fig. 3.6.  3.3  Characteristic Length  The length scale x0 is also measured from the cross-sections. The digitisation and transformation process results in gridpoints that are in general, unevenly spaced. Thus, the cross-sections are first cubically interpolated to a uniform grid xi = i∆x with i = 0, 1 . . . N − 1 with spacing ∆x equal to the mean spacing of the raw cross-sections. To reduce spectral leakage when taking taking the Fourier transform, the height functions h(xi ) are then multiplied element-wise by the Hanning window function [4] wi =  1 − cos( 2π(i+1) N +1 ) 2  (3.4)  47  Chapter 3. Field Observations  0.6 0.5  t0 (days)  0.4 0.3 0.2 0.1 0 0.6  0.65  0.7  0.75  0.8 θ/π  0.85  0.9  0.95  1  0.85  0.9  0.95  1  (a) times scale  2  h0 (cm)  1.5  1  0.5  0 0.6  0.65  0.7  0.75  0.8 θ/π  (b) height scale  Figure 3.6: Height and time scales determined by fitting numerical crosssectional peak-to-peak height to experimental measurements. Error bars are computed using the Monte Carlo method described in the text.  48  Chapter 3. Field Observations The cross-sections are grouped according to their time since flattening, and the photo from which they were digitised. This results in two groups at t = 2 days and t = 5 days corresponding to photos taken from the front and back of the same pair of trenches. Only one group is available for the lone t ≫ 5 days trench because only one photo was taken. The length of the longest cross-section Lmax in physical units is determined for each group and all cross-sections in the group are zero-paded ˆ i ) with to that length. In this way the discrete Fourier transforms h(k k = {2πi/Lmax : i ∈ Z} for all the cross-sections so that they have equal grid spacing ∆k = 2π/Lmax . The power spectral density of each cross-section is computed and the PSDs in each group are averaged together. The resulting PSDs for one of the groups at each time are shown in Fig. 3.7. Unlike the numerical PSDs, a low frequency power law is not readily evident. This is consistent with a more qualitative observation of Fig.1.1 or Fig. 3.3 where significant low frequency fluctuations do not appear to be present. Nevertheless, for consistency with the analysis of the numerical cross-sectional PDSs, the experimental cross-sectional PSDs are also multiplied by k. As shown in Fig. 3.8, this in fact, aids in the unambiguous determination of the power spectral peak particularly for the t = 5 days PSD. The highest peak in these modified PSDs is taken to correspond to the characteristic length. The same polynomial interpolation method as described in Section 2.7 is used to determine the location of the peak to subgrid precision. The peak for the 2 day cross-section is found to be too ill-defined to obtain a meaningful measurement of the characteristic length and so is omitted from further calculations. The peaks from the undisturbed snow on the other hand, are the most clearly defined with a characteristic length of ℓ = 66 ± 5cm. The peaks after 5 days are discernible at the slightly shorter length scale of ℓ = 53 ± 3cm, though there is also a strong sub-harmonic peak. These limited data seems to suggest that suncups coarsen over time similar to the numerical solutions. More data is however needed to better backup such a claim. Again, less precise observations from the Whistler site support these figures for the characteristic length. The square used to calibrate the perspective transform discussed in Section B.4 is known to be approximately 10′ × 10′ . When the snow coordinate grid is calibrated to the characteristic length of the snow surface as discussed in Section B.6, this square is found to be approximately 6.5 characteristic lengths in width and height. This implies a characteristic length of ∼ 47 cm for the surface, which is consistent 49  Chapter 3. Field Observations  t = 2 days 25 20 15 10 5 0 t = 5 days  PSDcs (cm3 )  200 150 100 50 0 t ≫ 5 days  150  100  50  0 0.05  0.1  0.2 0.15 −1 k (cm )  0.25  0.3  Figure 3.7: Average PSDs of groups of snow surface cross-sections.  50  Chapter 3. Field Observations  t = 2 days 1.5  1  0.5  0 t = 5 days  k PSDcs (cm2 )  8 6 4 2 0 t ≫ 5 days  10 8 6 4 2 0 0.05  0.1  0.15 0.2 −1 k (cm )  0.25  0.3  Figure 3.8: Modified PSDs of the snow surface cross-sections. These are multiplied by a factor of k for consistency with the analysis of the numerical data. This also allows a more unambiguous determination of the power spectral peak particularly for the 5 day PSD. 51  Chapter 3. Field Observations with the above values. The error for the above characteristic lengths is estimated by taking the difference between values determined for photos taken from opposite sides of the trenches. Since there is only a single photo taken of the cross-sections in the undisturbed snow, their characteristic length is assumed to be measured with the same relative error as the 5 day cross-sections. There are thus two data points for the characteristic length. One is at t = 5 days, while the other is at significantly larger unknown time. Unlike with the height scale, this poses a problem as the transition from the linear characteristic length to the slightly longer nonlinear characteristic length shown in Fig. 2.16 is a gradual process. It is not known at what stage the t ≫ 5 day trace is and so an arbitrary, but plausible time cannot be selected. Indeed it is not even known for certain that the snow surface coarsens in the same way as the numerical solutions. Thus, only the characteristic length measured 5 days after flattening is used in the computation of the length scale x0 as defined in Eq. 2.2. To determine x0 the experimental length scale ℓexp = 53 ± 3cm must be compared to the relevant numerical length scale. To this end, the crosssectional PSDs are computed from numerical solutions of Eq. 2.7 for θ ∈ [0.6π, π] at the equivalent of 5 days after flattening as determined via the time scale t0 (θ). Unlike the cross-sectional PSDs from the long time nonlinear regime shown in Fig. 2.13, the time average is not taken in the present case. This is not done since 5 days is very shortly after nonlinear saturation has taken place and the characteristic length still has to grow to its full nonlinear value. The ensemble average however is taken over the PSDs of 4 solutions as in Section 2.7. PSDs are again multiplied by k to remove the k −1 behaviour at low frequency and bring out the peak at the characteristic length. Additionally, it is found that since the time average is not taken, it is necessary to remove some noise in the PSDs by convolving them with a narrow Gaussian of width 1/16 (2∆k for the L = 64π solutions and ∆k for the L = 32π solutions). The PSDs before and after multiplication by k and smoothing are shown in Fig. 3.9 for three values of θ. The position of the peaks in the PSDs for each θ and each L are determined again using polynomial interpolation. This gives the characteristic length at 5 days as a function of θ: ℓc (θ, 5 days/t0 (θ)). The length scale x0 that converts between the numerical and experimental data is computed by taking the ratio of the numerical characteristic length with the experimental  52  Chapter 3. Field Observations  1200 θ = 0.6π θ = 0.8π θ = 1.0π  1000  PSDcs  800 600 400 200 0  0  0.5  1.5  1 k  5  (a) unsmoothed raw PSD  200 θ = 0.6π θ = 0.8π θ = 1.0π  k PSDcs  150  100  50  0  0  1  0.5  1.5  k (b) smoothed PSD with low frequencies attenuated  53  Figure 3.9: Cross-sectional numerical PSDs at three values of θ at τ t0 (θ) = 5 days.  Chapter 3. Field Observations  6 5  x0 (cm)  4 3 2 1 0 0.6  0.65  0.7  0.75  0.8 θ/π  0.85  0.9  0.95  1  Figure 3.10: Length scale x0 of the numerical solution with respect to the experimental cross-sections. The value is only weakly dependent on θ since the cross-sections are taken from the early nonlinear √ regime where the characteristic length is uniformly very near to the 2π 2 predicted from the linear stability analysis of Section 2.3. value: ℓexp = 54 ± 6 cm. x0 (θ) =  ℓexp ℓc (θ,  5days t0 (θ) )  =  53 ± 3  ℓc (θ, 5days t0 (θ) )  cm  (3.5)  This scale factor x0 is plotted as a function of θ in Fig. 3.10.  3.4  Diffusion of Real Suncups  The webcam at the Whistler site provides long term time-lapse observations of the suncup dynamics. This allows the motion of the suncups to be measured experimentally in much the same way that the minima of individual suncups were tracked for the numerical solutions in Section 2.10. Of course, locating the position of individual suncups is more difficult because the data is limited to an image of the surface rather than full topographical data. The spatial variation in the brightness of the image is  54  Chapter 3. Field Observations related in a non trivial way to the actual topography via the lighting conditions and optical properties of the snow at the time of observation. To add to this difficulty, said lighting conditions were unsurprisingly not constant throughout the observation period. In the absence of a sophisticated machine vision algorithm which is beyond the scope of this thesis, the suncup minima are identified manually. A GUI is designed to step sequentially through images of the snow surface on each of the 49 days of observation. This is used to select and track the trajectory of 455 individual suncups in close proximity to the camera. The suncups selected are all tracked for the extent of their existence within the observation time and area. Once the trajectories are digitised, it is possible to transform them along with the surface images into isotropic coordinates so that both appear to be viewed from above as shown in Fig. 3.11. The details of this procedure including calibration of the coordinates to units of the characteristic length scale of the surface are found in Appendix B. The importance of this is that it allows meaningful measurements of distance to be made from the images in units of the characteristic length. In particular, this makes it possible to compute the RMS displacement of the suncups as a function of time in exactly the same way as described in Section. 2.10 for the numerical simulations. Since the experimental displacement is measured in units of characteristic lengths, it is already in the same units as the data plotted in Fig. 2.24(b). All that is required for a meaningful comparison between the numerical and experimental data is a scaling of the time coordinate of the numerical displacement curves for each θ into time units of days via t0 (θ). The scaled numerical data are plotted together with the experimental displacement in Fig. 3.12. A weighted fit of Eq. 2.17 to the experimental data results in an exponent of b = 0.52. The fits to the numerical displacement curves done in Section 2.10 are also redone, this time restricting the range of the fit to the length of the experimental data for more realistic comparison. This turns out to slightly increase the power law exponent for the numerical data from b = 0.58 as found in Section 2.10 to b = 0.64 with a standard deviation of 0.05 over different values of θ. If the exponents of the experimental and numerical data were identical, a comparison of the fitted values for D alone would suffice for an experimental estimate of θ. Since b differs however, the fitted lines on a loglog plot are not parallel. This means that the relative placement of the experimental fit with respect to the numerical curves will be weakly time dependent. In general, the numerical displacements may be expressed as a function 55  Chapter 3. Field Observations  Figure 3.11: Experimental suncup trajectories (red lines) in isotropic “snow coordinates”. The termination point of the trajectories shown (black dots) as well as the image onto which they are superimposed are on the 37th day of observation (August 29). Thus, the oldest trajectories are a maximum of 36 days old. Tick marks correspond to units of the surface characteristic length.  56  Chapter 3. Field Observations  θ = 1.0π θ = 0.9π θ = 0.8π θ = 0.7π θ = 0.6π expt.  ℓd (ch.len.)  100  ∆t1/2  10−1  100  101 ∆t(days)  Figure 3.12: RMS displacement of experimental and numerical suncup minima scaled to comparable units. The y-axis is measured in units of characteristic lengths determined in both cases from the peak of the radial power spectral density. The x-axis is measured in units of days where the numerical data have been scaled according by t0 (θ), which is shown in Fig. 3.6(a).  57  Chapter 3. Field Observations  1  θexp (∆t)/π  0.9 0.8 0.7 0.6 0.5 100  101 ∆t (days)  Figure 3.13: Variation in the experimental determination of θ over the observed range of ∆t. If numerical and experimental exponents b were exactly equal, this would be a constant. of time and θ: ℓd (∆t, θ) = (D(θ)∆t)b(θ)  (3.6)  Though the form of D(θ) and b(θ) is not known analytically, they are sampled at several values of θ ∈ [0.6π, π]. This allows Eq. 3.6 to be “solved” for θ(ℓd , ∆t) via numerical interpolation provided that the diffusion length lies in the range ℓd (θ, ∆t) ∈ [ℓd (0.60π, ∆t), ℓd (π, ∆t)], which is seen to be the case for the experimental trace in Fig. 3.12. In this way, θexpt (∆t) is determined as shown in Fig. 3.13. To understand how the error bars in Fig. 3.13 are calculated, recall that the error on the time scale t0 (θ) was computed using a Monte Carlo method which resulted in the set {t0 (θ)} for each θ. Since each of the numerical displacement curves ℓd (θ, t) is scaled according to t0 (θ), this generalises to a set of such curves {ℓd (θ, t)} each scaled by a slightly different t0 within the set {t0 }. Additionally, since the least squares fits used to determine b(θ) and D(θ) also have a small uncertainty associated with them, these are also incorporated into {ℓd (θ, t)}. This is achieved by perturbing b and D by normally distributed sets {δb(θ)} and {δD(θ)} whose standard deviation 58  Chapter 3. Field Observations is equal to the error on the fit. These sets, by construction, have the same number of elements (10, 000) as {t0 }. Hence they may be combined on an element-wise basis {ℓd (θ, t)} =  (D(θ) + {δD(θ)})  t{t0 } t0  b(θ)+{δb(θ)}  (3.7)  In other words, {ℓd (θ, t)} also has 10000 elements just like each of the sets used to compute it. Since the sets are all generated from different random seeds, they may be considered uncorrelated so that this combination is a reasonable thing to do. Additionally, the experimental data have a relative uncertainty of 6% associated with determining the characteristic length of the snow surface as discussed in Section B.6. This may be introduced as an element-wise multiplication of {ℓd (θexpt , t)} by 1 + {δX} where {δX} is a normally distributed set with standard deviation 0.06. Similarly, there are uncertainties in the fit parameters D(θexpt ) and b(θexpt ). The analogous set of experimental displacement curves is {ℓd (θexpt , t)} = (1 + {δX}) ((D(θexpt ) + {δD(θexpt )}) t)b(θexpt )+{δb(θexpt )} (3.8) The interpolation of θexpt (t) is repeated for each pair of elements from Eqs. 3.7-3.8 to resulting in a set of curves {θexpt (t)} whose standard deviation is used as the size of the error bars in Fig. 3.13. The t dependence of θexpt is to be interpreted as a result of experimental uncertainty rather than a true variation of the θexpt parameter throughout the course of observation. As such, the final value for the experimental measurement of θexpt is taken as the average of θexpt (t) over the observation time. The average is weighted by a factor of 1/t so that it has uniform weighting on a log scale in time. An error on θexpt is found by taking the standard deviation of the entire set {θexpt (t)} over the entire course of experimental observation. This absorbs the variation of θexpt over t into its error estimate. The resulting value for θexpt may be used to compute the scales in Eqs. 2.1-2.3 by interpolating the curves in Figs. 3.6,3.10. Error estimates on the interpolated values of h0 (θexpt ) and t0 (θexpt ) are computed by interpolating the set of curves {h0 (θ)} and {t0 (θ)} from the original Monte Carlo fitting described in Section 3.2 at the set of points {θexpt }. On the other hand {x0 (θ)} must be generated by perturbing the points on x0 (θ) by normal 59  Chapter 3. Field Observations distributions with standard deviation equal to the error. The standard deviations of the resulting sets {h0 ({θexpt })}, {t0 ({θexpt })} and {x0 ({θexpt })} are taken as the errors for h0 (θexpt ), t0 (θexpt ) and x0 (θexpt ) respectively. The values of h0 (θexpt ), t0 (θexpt ) and x0 (θexpt ) together with estimates on their error are used to compute the coefficients c1,2,3,4 of the fully dimensional Eq. 1.1 by inverting Eqs. 2.1-2.3:  c1 = c2 = c3 = c4 =  x20 t0 x40 t0 x20 cos(θ) t0 h0 x40 sin(θ) t0 h0  (3.9) (3.10) (3.11) (3.12)  The resulting measured values are given in Table 3.1  3.5  Contribution of Suncups to the Ablation Rate  Recalling the observation in Section 2.9 that |∇h|2 is the only differential term that contributes to the height recession, this implies that ∂t h = −c0 − c3 |∇h|2  (3.13)  This is significant because the increase in ablation that results from |∇h|2 means that a rough snow surface has a larger coalbedo than a flat snow surface. The coalbedo is defined as the fraction of incident sunlight that is absorbed by the surface. Given the numerical measurement of |∇χ η|2 (θ) in dimensionless units in Fig. 2.14, |∇χ η|2 (θexpt ) may be interpolated together with an error estimate via the standard deviation of |∇χ η|2 ({θexpt }). This may be converted into experimental units with a factor of (h0 /x0 )2 . Together with Eq. 3.11 this allows Eq. 3.13 to be evaluated as ∂t h = −c0 − cos(θexpt )  h0 t0  |∇χ η|2 (θexpt )  (3.14)  60  Chapter 3. Field Observations  250  |∆h| cm  200 150  100 50 0 0  10  20  30  40  50  t (days) Figure 3.14: The net displacement of the snow surface at the Whistler site determined from the dynamic perspective transform parameters. The y-axis is scaled to units of centimetres using the value of x0 in Table 3.1and the numerically measured radial characteristic length. The fitted line has a slope of 4.0 ± 0.3cm/day which measures the rate of surface recession assuming the surface does not move laterally. Using the values computed above, the second term contributes 2.0 ± 0.5cm/day to the ablation rate. In order to put this measurement in context, the rate at which the snow surface recedes is measured at the Selkirk site to be ∂t h = 4.5cm/day. Additionally, the method by which the location of the snow surface is visually tracked at the Whistler site as described in Section B.5 results in a net surface displacement as a function of time R0 (t) = X0 (t)2 + Y0 (t)2 + Z0 (t)2 ) which is plotted in Fig. 3.14. The slope of this line is found to be 4.0 ± 0.3cm/day using the value for x0 in Table 3.1 corroborating the Selkirk measurement. However, since it is known that the camera at the Whistler site moved during the course of observations and since the Selkirk measurement is made by more direct means the Selkirk value is considered the final measurement of the recession rate. This would suggest that suncups are responsible for 45% of the surface recession rate. This is a surprisingly large number for surface topography 61  Chapter 3. Field Observations with such a small aspect ratio of h0 /x0 = 0.1. This would suggest that the measurements of one or both of θexpt and h0 are too large and/or t0 is too small. Additionally, it is possible that additional effects need to be incorporated into the model. For example, the diffusion rate of the suncups is augmented by the motion of the sun, or by stochastic processes such as rain altering the snow surface.  62  Chapter 3. Field Observations  x0 h0 t0 θ c0 c1 c2 c3 c4  5.2 ± 0.4cm 2.9 ± 0.2cm 0.39 ± 0.08days (0.79 ± 0.05)π 2.5 ± 0.5 84.2 ± 29cm2 /day (2.6 ± 1.3) × 103 cm4 /day 28 ± 12cm/day (1.0 ± 0.5) × 103 cm3 /day  Table 3.1: them  Measured scales and the coefficients for Eq. 1.1 derived from  63  Chapter 4  Conclusion This thesis has compared numerical simulations of the dimensionless Eq. 2.7 to field observations of suncups in order to determine length, time and height scales x0 , t0 and h0 in addition to the dimensionless nonlinear parameter θ. The scales x0 , t0 and h0 were determined by straightforward comparison of characteristic feature sizes and growth time from a relatively flat surface. Since these scales are of the same order of magnitude for all θ simulated, they are sufficient to give a rough estimate of the coefficients of the linear terms c1,2 as well as the combined strength of the nonlinear terms c23 + c24 as shown by Eqs 3.9-3.12. Knowledge of θ however is necessary to calculate the individual nonlinear coefficients c3,4 and to get improved estimates of the linear coefficients c1,2 . Experimental measurement of θ requires long term observation of the nonlinear dynamics of the snow surface. From such observations, the rate of “diffusion” of suncups on the surface was determined and compared to that of numerical solutions to Eq. 2.7 for a range of θ. It was thus found that an intermediate value of θ was required to reproduce the nonlinear dynamics of the snow surface. This intermediate value indicates that both c3 and c4 must be nonzero (i.e., that both nonlinear terms are required to minimally describe the formation and dynamics of suncups). The value of c3 is of particular interest as it indicates the degree to which the |∇h|2 term adds to the rate of recession of a flat snow surface c0 . Since this additional recession is most likely caused by an increase in energy absorption, it implies a corresponding increase in the coalbedo. It is a striking result that information about the energy balance of the snow surface can be determined from the motion of individual suncups formed on it. In particular, a faster rate of suncup motion is associated with a higher absorption of solar energy. It appears however that the full picture is not yet understood. The measured value of c3 |∇h|2 is surprisingly large compared to a directly measured ∂t h . It is possible that additional non-optical effects such as mechanical collapse of the surface or water transport significantly contribute to the measured size of c3 |∇h|2 . 64  Chapter 4. Conclusion An other discrepancy is the relative absence of low frequency fluctuations in the observed snow surface when compared to the inverse square power law in the PSD of the equation. It is likely that low frequencies play a significant role in the dynamic behaviour of the equation. Their relative absence in the physical system may indicate a significant shortcoming of the model that needs to be addressed to make further progress in understanding how it relates to snow. In addition to the specific snow problem, much has been learned about the PDE used to model it. It has been shown that the characteristic length measured in several different ways increases approximately linearly with the ratio tan θ of the dimensionless nonlinear terms at long times. It is also shown that the over the range of θ investigated, the solutions, given time, consistently adopt an approximately inverse square power-law behaviour in the spatial frequency k which dominates at low frequencies. Dynamic features of the equation, particularly with respect to the parameter θ have also been discovered. Minima in the solution undergo a random walk on the surface and the RMS displacement as a function of time goes as approximately the square root of time for all θ. An increased time-scale associated with decreasing θ (or increasing the coefficient of ∇2 |∇h|2 ) is quantified in the form of a “diffusion coefficient” which decreases exponentially with θ. Similarly, the steady state rate at which minima are created or destroyed also decreases exponentially with θ. These last two results are surprising since the dependence is on the relatively arbitrary parameter θ. It would be interesting to investigate if these relations persisted into the θ ∈ [π/2, 0.6π] range not investigated in this thesis. The analysis described here indicates the importance of dynamic time dependence in pattern formation phenomena. While the Kuramoto-Sivashinsky equation (with c4 = 0) is sufficient to describe the characteristic length, height and time scale, the addition of a higher order nonlinear term is required to mimic the full nonlinear dynamics of suncups. Furthermore, the equation is found to have interesting properties in its own right, particularly as changes are made to the ratio of the nonlinear terms. Since equations describing pattern formation frequently arise out of some type of low order expansion, this suggests the importance of understanding the long term nonlinear dynamics of both physical system and the equations used to model it in order to gain some physical insight  65  Bibliography [1] Glacier Ice. University of Washington Press, Seattle, revised edition, 2006. [2] Albert-Laszlo Barabsi and Harry Eugene Stanley. Fractal Concepts in Surface Growth. Cambridge University Press, 1995. [3] V. Bergeron, C. Berger, and M. D. Betterton. Controlled irradiative formation of penitentes. Phys. Rev. Lett., 96:098502, 2006. [4] Ronald N. Bracewell. The Fourier Transform and Its Applications. McGraw-Hill, 3 edition, 1999. [5] Richard L. Burden and J. Douglas Faires. Numerical Analysis. Thomson, 8 edition, 2005. [6] Ingrid Carlbom and Joseph Paciorek. Planar geometric projections and viewing transformations. Comp. Sur., 10:465, 1978. [7] S. M. Cox and P. C. Matthews. Exponential time differencing for stiff systems. J. Comp. Phys., 176:430, 2002. [8] A. K. Kassam and L. N. Trefethen. Fourth order time-stepping for stiff PDEs. SIAM J. Sci. Comp., 26:12141233, 2005. [9] Yoshiki Kuramoto and Toshio Tsuzuki. Persistent propagation of concentration waves in media far from thermal equilibrium. Prog. of Th. Phys., 55:356, 1976. [10] J.C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright. Convergence properties of the nelder-mead simplex method in low dimensions. SIAM J. Opt., 9(1):112, 1998. [11] T.C. Lin and Pu Qun. On the formation of regmaglypts on meteorites. Fluid Dynamics Research, 1:191, 1987. [12] C. Misbah and A. Valance. Sand ripple dynamics in the case of out-ofequilibrium aeolian regimes. Euro. Phys. J. E, 12:523, 2004. 66  Bibliography [13] J. Muoz-Garcia, M. Castro, and R. Cuerno. Nonlinear ripple dynamics on amorphous surfaces patterned by ion beam sputtering. Phys. Rev. Lett., 96:086101, 2006. [14] R. A. Peterson and W. B. Krantz. Differential frost heave model for patterned ground formation: Corroboration with observations along a north american arctic transect. J. of Geophys. Res., 113:G03S04, 2008. [15] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes in C: The Art of Scientific Computing. Cambridge University Press, 2nd edition, 1988. [16] Martin Raible, Stefan J. Linz, and Peter Hanggi. Amorphous thin film growth: Effects of density inhomogeneities. Phys. Rev. E, 64:031506, 2001. [17] J. J. Rhodes, R. T. Armstrong, and S. G. Warren. Mode of formation of ”ablation hollows” controlled by dirt content of snow. J. Glaciology, 33:135, 1987. [18] E. B. Saff and A. D. Snider. Fundamentals of Complex Analysis for Mathematics, Science, and Engineering. Prentic Hall, Upper Saddle River, New Jersey, 2nd edition, 1993. [19] A. Savitzky and Marcel J.E. Golay. Smoothing and differentiation of data by simplified least squares procedures. An. Chem., 36:1627, 1964. [20] G. I. Sivashinsky. Non-linear analysis of hydrodynamic instability in laminar flames–I. derivation of basic equations. Act. Astro., 4:1177, 1977. [21] Tom Tiedje, K. A. Mitchell, B. Lau, A. Ballestad, and E. Nodwell. Radiation transport model for ablation hollows on snowfields. J. of Geophys. Res., 111:F02015, 2006. [22] S. G. Warren, R. E. Brandt, and P. O. Hinton. Effect of surface roughness on bidirectional reflectance of antarctic snow. J. of Geophys. Res., 103:25789, 1998. [23] W.J. Wiscombe and S.G. Warren. A model for the spectral albedo of snow. J. of Atm. Sci., 1921.  67  Appendix A  Numerical Methods As with most non-linear PDEs, general closed form solutions of Eq. 2.7 are not feasible. This makes the numerical approximation of the solutions an important tool in understanding their properties and behaviour. Even such numerical solutions however must be implemented carefully to ensure both accurate and efficiently computed solutions. To this end, a spectral method is used for computing spatial derivatives and a fourth order Runge-Kutta exponential time differencing method with adaptive time step size is used to advance the solution in time. This chapter describes the implementation and advantages of such an algorithm.  A.1  Spectral Differentiation  Spectral differentiation takes advantage of the diagonalisation of differential operators in Fourier space and the efficiency of the fast Fourier transform in transferring to and from real space to carry out non-linear multiplication. This can be seen by taking the discrete Fourier transform F of Eq. 2.7 on an N × N 2D computational grid of size L × L ∂τ ηˆi,j  =  |κi,j |2 − |κi,j |4 ηˆi,j + −cos (θ) + sin (θ) |κi,j |2 F  F −1 {Iκˆ η}  2 i,j  (A.1)  √ where I = −1, ηˆi,j = ηˆ(κi,j , τ ) = F{η(χ, τ )}(κi,j ) and κi,j = 2π/L(iˆ x+ j yˆ) given perpendicular unit vectors x ˆ and yˆ and indicies i, j = −N/2 . . . N/2− 1 for even N . The linear terms are now greatly simplified to a straightforward wave number multiplication. This multiplication can be carried out in O(N 2 ) time. Additionally, though the use of successive forward and backward Fourier transforms to compute the nonlinear term may at first seem cumbersome,the efficiency of the fast Fourier transform allows it to be completed in O(N 2 log(N )) time. This is in contrast to computing the equivalent  68  Appendix A. Numerical Methods convolution directly in Fourier space, which would have a time complexity of O N 3 . The time complexity of computing the non-linear terms using the fast Fourier transform in this way still however dominates the numerical evaluation of the right hand side of Eq. A.1 such that the overall time complexity is O(N 2 log(N )). This might be compared to an equivalent finite difference method for which the time complexity is O(N 2 ). This would suggest that a finite difference method is more favourable in terms of speed of execution, however the merits of using a spectral method become clear if its accuracy is considered. While a typical finite difference scheme would have an error which goes to zero with increasing gridpoints N as O((L/N )3 ), the error of the spectral method above scales as O((L/N )N ). That is finite difference methods have slow power law convergence, while spectral methods have exponential convergence.  A.2  Exponential Time Differencing  Exponential time differencing takes advantage of the fact that the solution of just the linear terms of a nonlinear PDE can be arrived at analytically. To better understand this, we first rewrite Eq. A.1 as ∂τ ηˆi,j (τ ) = λi,j ηˆi,j (τ ) + γi,j N [ˆ η (τ )]i,j  (A.2)  where λi,j = |κi,j |2 −|κi,j |4 , γi,j = −cos(θ)+sin(θ)|κi,j |2 , and N [ˆ η (κ, t)]i,j = −1 2 F{|F {iκˆ η (κ, τ )}| }i,j . Taking only the linear terms in Eq. A.2 ∂τ ηˆi,j (τ ) = λi,j ηˆi,j  (A.3)  The analytical solution is simply ηˆi,j (τ ) = eλi,j τ ηˆi,j (0)  (A.4)  Firstly, writing the solution to the linear terms in this way draws attention to a difficulty that arises for large κ. The smallest timescale in Eq. A.4, τ ∼ κ−4 places a limit on the largest integration timestep that can be used to faithfully integrate the solution. This is in spite of the fact that the Fourier modes in which we are primarily interested are of order κ ∼ 1 and hence of timescale τ ∼ 1. Thus, a numerical method that mitigates this problem is desirable.  69  Appendix A. Numerical Methods Secondly, it also might be hoped that there exists a numerical method for time stepping Eq. A.2 which preserves the exact solution given in Eq. A.4 in the limit that the N [ˆ η ] → 0. Indeed, these goals are complimentary and can be realised by multiplying Eq. A.2 by the integration factor e−λi,j τ e−λi,j τ ∂τ ηˆi,j = λi,j e−λi,j τ ηˆi,j + γi,j e−λi,j τ N [ˆ η ]i,j  (A.5)  which can be rearranged η ]i,j ∂τ e−λi,j τ ηˆi,j = γi,j e−λi,j τ N [ˆ  (A.6)  Both sides of Eq. A.6 can then be integrated from time τ to τ + ∆τ . τ +∆τ  τ +∆τ  ∂s e−λi,j s ηˆi,j (s) ds = γi,j  e−λi,j s N [ˆ η (s)]i,j ds  (A.7)  τ  τ  Using the fundamental theorem of calculus to evaluate the integral on the left hand side and shifting the integration variable in the remaining integral on the right hand side, we may isolate the advanced time solution ηˆi,j (τ + ∆τ ) to obtain ∆τ  ηˆi,j (τ + ∆τ ) = eλi,j ∆τ ηˆi,j (τ ) + γi,j eλi,j ∆τ  e−λi,j s N [ˆ η (τ + s)]i,j ds (A.8)  0  We may explicitly discretise the time variable defining τn = n∆τ and ηˆi,j,n = ηˆi,j (τn ) ∆τ  ηˆi,j,n+1 = eλi,j ∆τ ηˆi,j,n + γi,j eλi,j ∆τ  e−λi,j s N [ˆ η (τn + s)]i,j ds  (A.9)  0  The first term on the right hand side corresponds to the linear analytical solution mentioned previously while the second term corresponds to the remaining non-linear portion of the equation that must be integrated numerically. It can easily be seen that this separation of the linear and nonlinear time stepping in Eq. A.9 achieves the desired end of preserving the exact solution in the limit N [ˆ η ] → 0. Additionally, because the integration of the small timescales resulting from large κ in the linear terms is handled analytically, this will not be a limiting factor on the size of the numerical time step. 70  Appendix A. Numerical Methods It should also be noted that Eq. A.9 is completely integration-method agnostic. Any numerical algorithm for evaluating the integral in the second term may be used. It is found however, that for the particular case of Eq. A.1, a fourth order Runge-Kutta method is an ideal choice due to its relative stability. Such a method has been derived for the case of Eq. A.9 which takes account of the particular form of the integral, but is non-trivial to derive without the aid of a computer algebra system [7]. It is given by, ai,j,n = eλi,j ∆τ /2 ηˆi,j,n + αi,j N [ˆ ηn ]i,j λi,j ∆τ /2  bi,j,n = e  ηˆi,j,n + αi,j N [an ]i,j  ci,j,n = eλi,j ∆τ /2 ai,j,n + αi,j (2N [bn ]i,j − N [ˆ ηn ]i,j ) λi,j ∆τ  ηˆi,j,n+1 = e  (A.10) (A.11) (A.12)  ηˆi,j,n + βi,j N [ˆ ηn ]i,j  +σi,j (N [an ]i,j + N [bn ]i,j ) + ρi,j N [cn ]i,j  (A.13)  where eλi,j ∆τ /2 − 1 λi,j  αi,j  = γi,j  βi,j  = γi,j  σi,j  = 2γi,j  ρi,j  = γi,j  (A.14)  (4 − 3λi,j ∆τ + λ2i,j ∆τ 2 )eλi,j ∆τ − 4 − λi,j ∆τ ∆τ 2 λ3i,j  (−2 + λi,j ∆τ )eλi,j ∆τ + 2 + λi,j ∆τ ∆τ 2 λ3i,j  (4 − λi,j ∆τ )eλi,j ∆τ − 4 − 3λi,j ∆τ − λ2i,j ∆τ 2 ∆τ 2 λ3i,j  (A.15) (A.16) (A.17)  This however is not the final word, as it might be noted that α, β, σ, and ρ become difficult to evaluate numerically near the value λi,j = 0, which is a perfectly reasonable value for it to take. Indeed, this occurs in the case of Eq. A.1 at |κ| = 1 which is the same order as the characteristic length. A robust solution has been proposed [8] which makes use of the Cauchy’s integral formula [18]. f (z) 1 dz (A.18) 2πI Γ z − z0 √ Here, I is again used to denote −1 and Γ is a contour in the complex plane enclosing z0 . Eqs A.14-A.17 are analytically continued into the complex plane and substituted for f into Eq. A.18 with λi,j corresponding to f (z0 ) =  71  Appendix A. Numerical Methods the variable z0 . The integral is evaluated over the unit circle centred at λi,j which is parameterised as Γ = {z(φ) = λi,j + eiφ : φ ∈ [0, 2π]}. For example, αi,j is evaluated as αi,j = γi,j  2π  1 2π  0  Iφ  e(λi,j +e )∆τ /2 − 1 dφ λi,j + eIφ  (A.19)  This however can be simplified owing to the fact that the imaginary part of the integrand cancels in the upper and lower portions of the circle. αi,j = γi,j  1 π  Iφ  π 0  ℜ  e(λi,j +e )∆τ /2 − 1 λi,j + eIφ  dφ  (A.20)  This integral is evaluated numerically using the midpoint rule αi,j = γi,j  1 M  M −1 k=0  Iπ(k+0.5)/M  ℜ  )∆τ /2 − 1 e(λi,j +e λi,j + eIπ(k+0.5)/M  (A.21)  It is found that M = 16 is sufficient to achieve convergence to within machine precision. Eqs. A.15-A.17 are all evaluated in an analogous way. This method is found to be excellent at integrating Eq. A.1 for θ ∈ [0.6.1]π, however as θ → π/2, smaller and smaller timesteps are found to be necessary to achieve the same accuracy. The reason for this likely has to do with the prefactor γi,j to the integral in Eq. A.9. As θ → π/2, the sin(θ)κ2i,j component of γi,j gets stronger thereby amplifying the high κi,j components of the numerically evaluated integral, and with them, their error. This is essentially the same small time scale problem that was solved for the linear terms by utilising the exponential time differencing method. It however is of lesser severity, as we are now dealing with time scales τ ∼ 1/κ2 rather than τ ∼ 1/κ4 . In future work, it may be worth investigating ways to mitigate this τ ∼ 1/κ2 problem as well. A first step would be to utilise a black box implicit integrator such as LSODA to evaluate the integral and prefactor in Eq. A.9. This would still retain the advantages of the exponential time differencing method, however unlike the fourth order Runge-Kutta scheme outlined above, it would not make use of the special form of the integral. For this, it might be advantageous to investigate adapting an implicit timestepping scheme directly as has been done for the explicit RK4.  72  Appendix A. Numerical Methods  A.3  Adaptive Time Stepping  It is found empirically that the time step size necessary to integrate Eq. A.1 to a given accuracy can vary widely over the course of integration and as noted above, over different values of θ. This is especially the case when white noise random initial conditions are used which, by definition, contain equal amplitudes of all resolvable spatial frequencies. This is in contrast to typical solutions of Eq. A.1 for which such transients have been allowed to die out and high frequency modes decrease in amplitude exponentially with |κ|. As long as such high frequency modes have significant amplitude, the time step will be restricted to a smaller value particularly for θ ← π/2. When the high frequency modes have been damped out however, the timestep can typically be significantly larger. In order to achieve a relatively constant accuracy over the entire integration time, as well as over the range of θ investigated, an adaptive time step algorithm is implemented. This involves periodically generating an error estimate. If the error is above a certain tolerance, the timestep is decreased and the integration is repeated from the last verified time using the smaller stepsize. On the other hand, when the error is below the threshold the timestep is considered successful and integration proceeds. The error together with a knowledge of the order of convergence for the method can be used to estimate the amount by which the timestep can safely be increased or decreased to achieve the desired accuracy in a minimum computation time. Due to the generality with which it can be implemented, a step doubling scheme is used to generate the error estimate necessary for the adaptive timestep scheme. Every two timesteps are accompanied by single timesteps of double the size and the difference in the results is taken as the error. This naturally, gives an error for each Fourier mode κi,j . ˆ i,j,n (∆τ ) ∆ˆ ηi,j,n = ηˆi,j,n (2∆τ ) − h  (A.22)  Neglecting the time saved from increasing the timestep when it is permissible, it at first appears that such a scheme will take 1.5 times as long as a scheme which simply directly applies Eqs. A.10-A.13 at a constant time step. This however can be improved upon by noting that we may eliminate an evaluation of N [ˆ ηn ]i,j in Eq. A.10 since the same result can be used for both the first ∆τ timestep as well as the 2∆τ step. This results in the improved execution speed of 1.38 times a constant timestep scheme [15]. As mentioned above, the typical amplitude of the Fourier modes falls off 73  Appendix A. Numerical Methods exponentially, thus the amplitudes, and hence their errors span many orders of magnitude. This makes it desirable to enforce a relative tolerance on ∆ηi,j,n proportional to ηi,j,n by a factor Tr . Below some absolute error however, this becomes impractical. Due to machine roundoff error, we simply cannot expect the same relative tolerance to be enforced across more than a few orders of magnitude. It is also simply not necessary ensure relative accuracy of modes smaller than some cutoff as they simply don’t have a significant effect on the solution. Thus, we would like also to implement the enforcement of an absolute tolerance Ta bellow which the relative tolerance is relaxed. This is achieved by enforcing the inequality ǫn ≡ max i,j  |∆ˆ ηi,j,n | 2 Ta N + Tr |ˆ ηi,j,n |  ≤1  (A.23)  The factor of the number of gridpoints over the whole domain N 2 is necessary due to the fact that N −1 i,j=0  |ηi,j |2 =  1 N2  N −1 i,j=0  |ˆ ηi,j |2  (A.24)  where ηi,j enumerate the real space gridpoints corresponding to the Fourier space solution ηˆi,j . This ensures, for example, that ǫn is not directly dependent on the number of grid points used to represent the solution. As stated previously, not only can the error estimate ǫn be used to evaluate the quality of the current solution, it can also be used to estimate the step size that should be used for the next timestep. Given that the Runga-Kutta time stepping method used is fourth order (i.e., has an error of O(∆τ 5 )), if the current time step resulted in an error of ǫn , we can make the reasonable assumption that the step size we should use for the next timestep to achieve ǫn+1 ∼ 1 ∆τn+1 = s∆τn  1 ǫn  1 5  (A.25)  where s is a “safety factor” taken (somewhat arbitrarily) to be 0.9 in the present case. The fourth order Runge-Kutta scheme above however does not lend well to blindly implementing Eq. A.25 since this means that Eqs. A.14-A.17 must all be recalculated for a new ∆τ at each time step. Instead, the time step  74  Appendix A. Numerical Methods is halved whenever the following condition more stringent condition than Eq. A.23 is not met ǫn < s5  (A.26)  conversely, the stepsize is doubled when the condition ǫn <  s 2  5  (A.27)  Thus ǫn should typically obey s 5 < ǫn < s5 (A.28) 2 The inequality on the right is strict, in the sense that any time step not satisfying it is redone at a smaller stepsize. The inequality on the left on the other hand, will by definition not hold for timesteps immediately preceding a step doubling. Of course this is not a problem, since it simply makes the given step slightly more accurate than desired. This halfing/doubling scheme has the advantage that ∆τ is changed and hence Eqs. A.14-A.17 are recalculated only when strictly necessary. Furthermore, since it is necessary at any given time to have coefficients Eqs. A.14-A.17 for both ∆τ and 2∆τ in order to compute Eq. A.22, it is only necessary to recalculate one such set of coefficients, while the other may be reused. This works rather well in practice for integrating Eq. A.1, as there are extended periods (particularly after transients have died out) during which the timestep does not change. For situations in which it changes more frequently, this method may not be as practical. On the other hand, given that all possible timesteps are within a power of two, it is not unlikely that the same timestep may be revisited. Thus, computed values of the coefficients Eqs. A.14-A.17 for various ∆τ might be retained for later use. This wouldn’t be terribly impractical given the typically large amounts of memory available in most modern computers. Computations utilising this scheme without retaining previously computed coefficients had ample memory left over to do so. Indeed, execution time, not available memory was the primary factor limiting the size of computation that could be practically undertaken. ˆ i,j,n at each A final note should be made regarding the computation of h timestep. The value computed using the timestep of 2∆τ need not simply be discarded. It may be used to improve the order of convergence of the  75  Appendix A. Numerical Methods value obtained using ∆τ by making use of “local extrapolation” [15]. (16ˆ ηi,j,n (∆t) − ηˆi,j,n (2∆τ ) + O(η 6 ) (A.29) 15 Note that this changes the truncation error of ηˆi,j,n from fourth to fifth order in ∆τ . This may be derived from the Taylor series of both ηˆi,j,n (∆τ ) and ηˆi,j,n (2∆τ ). In future work, it may be worth while investigating the existence and potential superiority of an embedded Runge-Kutta method to integrate and estimate the error on Eqs. A.10-A.13. Runge-Kutta methods are not unique in that there can exist multiple possible choices of the coefficients Eqs. A.14A.17 to achieve accuracy of a given order. It has been shown [15] that there exist less specialised Runga-Kutta schemes that make use of two different sets of coefficients to produce estimates of ηˆi,j,n+1 to fourth and fifth order using the same function evaluations (analogous to N in Eqs. A.10-A.13). This could, in principle mean an even more efficient adaptive time step method than outlined above. ηˆi,j,n =  76  Appendix B  Perspective Transformation Experimental data discussed in previous chapters come predominantly from photographs. This poses a challenge to making quantitative measurements, as the scale is not constant throughout the photographically depicted threedimensional space. In order to derive physical measurements, calibration is necessary together with some assumptions about geometry of the object being imaged (in this case, the snow surface). Lacking complete information, some calibration adjustments are must be made “by eye”. While this is admittedly not ideal from a scientific method perspective, such augmentation of more objective methods is found to be both necessary and surprisingly powerful. As usual, hindsight suggests methods by which the data acquisition might be improved to yield better results. The fundamental theoretical tool used in analysing the photographs is the perspective projection [6]. This is a member of the class of projections known as planar geometric projections. Planar geometric projections consist of a non-unique mapping between points in three-dimensional space and a two-dimensional plane. Lines passing through points on the plane and the centre of perspective map all points they pass through back to their originating point on the plane. All planar projections have the property that they map straight lines in three dimensional space to straight lines on the plane. A planar projection for which the centre of perspective is at a finite distance is, by definition, a perspective transformation. The perspective transformation is used as a simplified model for the creation of a photograph from a three-dimensional scene. This neglects any distortion that arises from passing the image through a camera lens which is found to be negligible for present purposes. This chapter details how such a model may be used to gain usable spatial information from the photographic image of a two-dimensional snow surface embedded in three-dimensional space.  77  Appendix B. Perspective Transformation  B.1  Homogeneous Coordinates  As stated above, the scale at a point in a photographic image varies depending on the position of the corresponding point within the original threedimensional scene. This can be seen as the direct result of applying the perspective transformation in the process of capturing the image. It is convenient to encode this spatially dependent scale directly into the coordinate vectors used to describe positions in the image. Where x, y and z are the conventional Cartesian coordinates of the position vector x, a fourth s may be added for the scale such that   x y   (B.1) xh =  z  s  This set of four coordinates for representing points in three-dimensional space are referred to in the literature as homogeneous coordinates [6]. Interpretation of s as the scale suggests s = 1 for a coordinate system that has not yet been subject to the perspective transformation. The first three elements x, y and z may be interpreted similarly to their Cartesian counterparts with the exception that the scale has not yet been applied. To completely convert a vector of homogeneous coordinates back to standard three-dimensional coordinates its first three elements must be divided by s   x 1  y x= (B.2) s z  B.2  Matrix Formulation of the Perspective Transformation  The effect of the perspective projection on a point in three-dimensional space can be described as a composition of linear transformations acting on its homogeneous coordinate vector xh . The first transformation applied is a translation, or a uniform shift of all points within the space by specified amounts X0 , Y0 and Z0 in each of the three spatial directions. Use of homogeneous coordinates conveniently  78  Appendix B. Perspective Transformation allows the expression of this transformation in matrix form as   1 0 0 X0 0 1 0 Y0   T (X0 ) =  0 0 1 Z0  0 0 0 1  (B.3)  A rotation of angle φ1 is then performed around the z-axis followed by a rotation of θ around the new x-axis. The z-axis after applying this second rotation points in the direction of the camera. Equivalently, this axis is parallel to the normal of the plane of projection. The final rotation of φ2 about this axis leaves its orientation unchanged. In matrix form, all three rotations are denoted by   cos(φ2 ) −sin(φ2 ) 0 0 0 0 1 0 sin(φ2 ) cos(φ2 ) 0 0 cos(θ) −sin(θ) 0 0   R(φ1 , θ, φ2 ) =   0 0 1 0 sin(θ) cos(θ) 0 0 0 0 0 1 0 0 0 1   cos(φ1 ) −sin(φ1 ) 0 0 sin(φ1 ) cos(φ1 ) 0 0   (B.4)  0 0 1 0 0 0 0 1   It might be noted that application of the above linear operators has no effect on the fourth coordinate s which remains equal to unity up to this point. Now however, having oriented the camera within the threedimensional space, the perspective transformation is applied, distorting the scale of the scene by making s depend on the spatial coordinate. If the centre of perspective is located at the Cartesian coordinates U0 V0 W0 in the translated and rotated coordinate system, the perspective transformation is denoted by the matrix   1 0 0 U0 0 1 0 V0   (B.5) P (U0 ) =  0 0 1 0 0 0 1/W0 1  The next step is to project the image onto the plane of projection. This amounts to reducing the number of coordinates by one which is achieved by  79  Appendix B. Perspective Transformation the matrix   1 0 0 0 Z = 0 1 0 0 0 0 0 1  (B.6)  In summary, the homogeneous coordinate uh = u v s of a point on the plane of projection may be expressed in terms of the homogeneous coordinate vector xh of a point in three-dimensional space via matrix multiplication uh = ZP (U0 )R(φ1 , θ, φ2 )T (X0 )xh  (B.7)  It is then possible to convert uh into Cartesian coordinates via Eq. B.2. An additional constant scaling factor of s0 can also included at this stage to convert from the natural units of measurement in the three-dimensional coordinate system to units of pixels within the projected image. u=  s0 uh u = v s vh  (B.8)  For convenience, the origin of the image coordinate system u = 0 is taken to be the centre pixel of the image.  B.3  Inversion of the Perspective Projection  The previous section described the procedure by which a point in three dimensional space may be mapped onto the plane of projection using a composition of linear transformations. What is needed however is the reverse, namely to map points which have been projected onto the image back to their origin in three-dimensional space. Unfortunately, this is not generally possible since the perspective transformation does not map a given x uniquely to a given u. Information is lost in the application of Eq. B.6 and Eq. B.8. Luckily however, the present problem is less general. The points of interest in three-dimensional space are not located arbitrarily, but are constrained to the plane of the snow surface. Such a point may be described by just a two-coordinate vector whose coordinate axes lie in the plane of the snow. x=  x y  (B.9)  80  Appendix B. Perspective Transformation This is trivially converted to a full three-dimensional homogeneous coordinate vector x by setting s = 1 and z = 0 so that the snow surface lies in the x-y plane.   x y   (B.10) xh =  0 1 In this case, coordinates taken from the image of the snow surface u are the left hand side of a system of two equations in the two variables x and y which identify the location of the corresponding point on the snow surface. s0 u(x, y) uh (x, y) = v(x, y) s(x, y) vh (x, y)  (B.11)  where     x uh (x, y) y   vh (x, y)  = ZP (U0 )R(φ1 , θ, φ2 )T (X0 )   0 s(x, y) 1  (B.12)  The expressions for u(x, y) and v(x, y) are too complex to reproduce here to any benefit, however they are perfectly calculable. This system of equations may be solved for x(u, v) and y(u, v) whose expressions are algebraically derivable though again, unwieldy. Having established that they may be practically computed, it is convenient to define the forward u = Qp (x) (B.13) and backward x = Qp−1 (u)  (B.14)  p = X0 φ1 θ φ2 U0 s0  (B.15)  perspective transforms where  B.4  Calibration  Given Q and Q−1 , there remains the task of finding their parameters p. This is somewhat simplified by placing the centre of perspective directly behind the origin of the plane of projection. Since the centre of perspective may 81  Appendix B. Perspective Transformation be interpreted as the position of the observer [6], this most closely models the creation of the image. A convenient result of choosing the centre of the image as the origin is that this effectively sets   0 U0 =  0  (B.16) W0  Starting with the minimal parameters X0 = 0, φ1 = θ = φ2 = 0, W0 = ∞ and (arbitrarily) s0 = 50, a portion of the image is transformed into the snow coordinates. Likewise, a unit grid spanning the extent of the transformed image xgrid i,j =  i j  (B.17)  is transformed back into image coordinates and superimposed on the original image. Adjustments are then made to θ, φ2 , and W0 so that the transformed snow surface appears at roughly uniform scale as though it is viewed from directly above and the grid superimposed on the original image appears to lie on the snow surface. Note that changing X0 and φ1 are not necessary to achieve this as they merely change the origin of the snow coordinate system or rotate it about its z-axis. Further refinement of p is achieved by using a 10′ × 10′ square marked off in rope on the snow surface as reference. The parameters θ, φ2 and W0 are tweaked so that opposing sides of the square are the same length and parallel in the snow coordinates (i.e., that the square is indeed a square) as shown in Fig. B.4. This is made easier to judge by orienting φ1 such that the sides of the square are parallel to the snow coordinate axes. In retrospect, the determination of θ,φ2 , and W0 could have been significantly improved if there were more a priori information about the snow surface. In particular, more reference points than just the square would improve the estimate of the parameters.  B.5  Changes to the Parameters Over the Course of Observation  Since the snow surface recedes as it moves, and there is some movement of the camera throughout the course of observation, the parameters of the perspective transformation are not exactly the same for images from different 82  Appendix B. Perspective Transformation  (a) Original image with grid indicating snow coordinates  (b) Transformed image in snow coordinates  83  Figure B.1: Demonstration of the forward (a) and backward (b) perspective transformation of the snow surface. The scale s0 has been set so that grid spacing is equal to the characteristic length using the methods described in Section B.6.  Appendix B. Perspective Transformation days. Repeating the above procedure to determine pd for each day d is not however possible since the rope that marked off the square did not remain in a square shape throughout the observation and eventually disappeared. Even if this were possible, the process would be painstaking and the error in the resulting parameters would be uncorrelated between days. This would artificially enhance the diffusion of the suncup minima. Instead, use is made of the suncup minima data points ui,d which link points on the snow surface (enumerated by i) between days d and d + 1. The mean square magnitude of the displacement vector over all points existing on both these days Qp−1 (ui,d ) − Qp−1 (ui,d+1 ) d d+1  (B.18)  should be zero since it is equally likely to point in all directions. Given a set of minima ui,d and ui,d+1 that exist on both days d and d + 1, together with known values pd for the parameters of the perspective transform on day d, the values for pd+1 that minimise the sum χ(pd+1 )2 = i  Qp−1 (ui,d ) − Qp−1 (ui,d+1 ) d d+1  2  (B.19)  are taken as the parameters for day d + 1. The minimum of χ(pd+1 ) is found using a Nelder-Mead simplex optimisation algorithm[10]. This procedure is applied recursively to determine the pd for each day of the observation period. It should be noted that only the X0 , φ1 , θ, φ2 , and W0 components of pd are allowed to vary in this minimisation. The parameters U0 , V0 and s0 are functions of the camera optics which are assumed constant.  B.6  Setting the Scale  Having established values of pd for the entire observation period, it is possible to transform each image into snow coordinates such that relative distances within the image may be reasonably estimated. Absolute distance measurements however, are not yet possible since s0 was set arbitrarily in Section B.4. Finding a more meaningful value for s0 requires measuring the characteristic length of the snow surface. To this end, rectangular sub-images containing exclusively the snow surface are extracted from the transformed  84  Appendix B. Perspective Transformation  ×10−4 5  k PSDr  4 3 2 1 0 0.5  1  1.5 2 −1 k/2π (ch.len. )  2.5  Figure B.2: Power spectral density of the transformed image of the snow surface. The characteristic length kc is indicated by the red dot which is scaled to 2π so that length is measured in units of the characteristic length. images for each day. Each of these sub-images is pixel-wise multiplied by a two-dimensional Hanning window defined for an N xM pixel image Ii,j as 2π(j+1) (1 − cos( 2π(i+1) N +1 ))(1 − cos( M +1 )) (B.20) 4 The windowed images are then zeropaded to the size of the largest one so that they all have the same dimensions. The radial power spectral density is computed for each image and averaged over all images. The result is found to have k −1 dependence at low frequencies which is removed by multiplying by a factor of k as shown in Fig. B.6. The resulting power spectrum has a peak at kc which corresponds to the characteristic length of the pattern on the snow surface. If this characteristic length is used for the unit of measure in the snow coordinates, this peak will occur at kc = 2π. This may be achieved simply by rescaling k → 2πk/kc . The snow coordinates will correspondingly be rescaled as x → xkc /2π. This scaling may be built in to the transformations Q and Q−1 using the scaling relations  wi,j =  85  Appendix B. Perspective Transformation  QX0 ,φ1 ,θ,φ2 ,U0 ,s0 (nx) = QnX0 ,φ1 ,θ,φ2 ,nU0 ,s0 /n (x) nQ−1  X0 ,φ1 ,θ,φ2 ,U0 ,s0  (u) = Q−1  nX0 ,φ1 ,θ,φ2 ,nU0 ,s0 /n  (u)  (B.21) (B.22)  Thus, given the unscaled set of parameters pd = X0 φ1 θ φ2 U0 s0  (B.23)  the rescaled values will be pd =  kc 2π X0  φ1 θ φ2  kc 2π U0  2π kc s0  (B.24)  This rescaling is uniformly applied to the parameters pd for all days. This ensures that all lengths in the snow coordinates are measured in units of the characteristic length.  B.7  Cross-sections  The cross-sections of the snow surface discussed in Section 3.1 are obtained from digitised photos of trenches dug in the snow as shown in Fig. B.7. Unlike the webcam photos whose analysis is described above, these crosssectional photos contain an object of precisely known length. An ice axe which can be seen in Fig. B.6 is present in all cross-sectional photos analysed and has a handle measured to be 60cm in length. This allows accurate measurements of the snow topography to be made in physical units. The first task is to determine the orientation of the snow surface in a way similar to that described in Section B.4. This involves setting values for the parameters p. In particular, these parameters will include an s0 calibrated by the ice axe which converts between pixels and centimetres. Determining the other parameters however is more difficult in the present case. Even though the dimensions of the ice axe handle are precisely known, this gives only two reference points as opposed to the four used for the webcam. There is thus a greater reliance on “by eye” adjustments. This could have been avoided had a square object of known length been used for caibration. The result is shown by the grid in Fig. B.6. The grid can be thought of as lying at the mean height of the snow surface. The next step is to determine the line lying within the plane of the grid that the cross-section intersects. This is achieved by performing a linear least squares fit to the digitised cross-sectional trace. All deviation of the 86  Appendix B. Perspective Transformation  Figure B.3: An example of a trench from which cross-sections of the snow surface were obtained. This trench is in an otherwise undisturbed patch of snow. The edge of the trench is digitised as shown. The grid shown indicates the plane of the snow surface and is scaled to 30cm/grid spacing. The 60cm handle of the ice axe on the far side of the trench is used as a scale for measurements. trace from this line of best fit is assumed to be in the vertical direction (i.e., parallel to the z-axis in the snow coordinates). The line of best fit is transformed into snow coordinates to determine its angle of rotation and offset from the x-axis. It is sufficient to transform the end points of the line ua and ub . xa,b = Qp−1 (ua,b )  (B.25)  The angle between this line and the x-axis is given by ∆φ1 = tan−1  yb − ya xb − xa  (B.26)  To make the line parallel to the x-axis, this can be built into the transform parameters p determined above by adding it to φ1 φ′1 = φ1 + ∆φ1  (B.27)  Note that, φ1 only rotates the snow coordinates about the z-axis so that the x-y plane of the rotated coordinate system still lies on the snow surface. The new set of parameters p′ with the rotation built in is then applied to the line of best fit so that its y-offset in the rotated snow coordinates may  87  Appendix B. Perspective Transformation be determined. x′a,b = Qp−1 ′ (ua,b )  (B.28)  Both x′a and x′b have the same y component ∆y since they are parallel to the x-axis thanks to the rotation. ∆y = xa,b · yˆ  (B.29)  It is possible to shift the origin of the snow coordinates such that ∆y = 0 by modifying p′ to form a new set of parameters p′′ . This places the line of best fit on the x-axis. The required p′′ is computed from p′ simply by adding ∆y to Y0 . Y0′ = Y0 + ∆y  (B.30)  The above procedure is carried out to transform each cross-sectional trace onto the x-axis in the snow coordinates. Since any deviation from the line of best fit is in the z direction in snow coordinates, and having constrained this line to the x-axis by shifting p appropriately, it follows that the y-coordinate of the cross-section is uniformly 0 for all points on the trace. Thus, we may modify Eqs. B.11-B.12 to read s0 u(x, z) uh (x, z) = v(x, z) s(x, z) vh (x, z)  (B.31)  where     x uh (x, z) 0  vh (x, z)  = ZP (U0 )R(φ1 , θ, φ2 )T (X0 )   z  s(x, z) 1  (B.32)  Which may be similarly solved to derive transformations P and P −1 that convert between image and snow-coordinates constrained to the x-z plane instead of the x-y plane as for Q and Q−1 . u = Pp (x)  (B.33)  Pp−1 (u)  (B.34)  x=  The x and y components computed by applying P −1 (x) to the digitised images points represent the horizontal and vertical positions respectively. 88  Appendix B. Perspective Transformation Furthermore, the s0 calibrated using the ice axe conveniently gives these values as measured in centimetres.  89  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items