The Formation of Lake Stars
The Formation of Lake Stars
Victor C. Tsai1∗ and J. S. Wettlaufer2†
Department of Earth & Planetary Sciences, Harvard University, Cambridge, Massachusetts, 02138 and
Departments of Geology & Geophysics and Physics, Yale University
New Haven, Connecticut 06520-8109
(Dated: 11 April, 2007)
Star patterns, reminiscent of a wide range of diffusively controlled growth forms from snowflakes
to Saffman-Taylor fingers, are ubiquitous features of ice covered lakes. Despite the commonality
and beauty of these “lake stars” the underlying physical processes that produce them have not been
explained in a coherent theoretical framework. Here we describe a simple mathematical model that
captures the principal features of lake-star formation; radial fingers of (relatively warm) water-rich
regions grow from a central source and evolve through a competition between thermal and porous
media flow effects in a saturated snow layer covering the lake. The number of star arms emerges
from a stability analysis of this competition and the qualitative features of this meter-scale natural
phenomena are captured in laboratory experiments.
I. INTRODUCTION
The scientific study of the problems of growth and form occupies an anomalously broad set of disciplines. Whether
the emergent patterns are physical or biological in origin, their quantitative description presents many challenging
and compelling issues in, for example, applied mathematics [1], biophysics [2], condensed matter [3] and geophysics
[4] wherein the motion of free boundaries is of central interest. In all such settings a principal goal is to predict the
evolution of a boundary that is often under the influence of an instability. Here we study a novel variant of such a
situation that occurs naturally on the frozen surfaces of lakes.
Lakes commonly freeze during a snowfall. When a hole forms in the ice cover, relatively warm lake water will flow
through it and hence through the snow layer. In the process of flowing through and melting the snow this warm water
creates dark regions. The pattern so produced looks star-like (see Figure 1) and we refer to it as a “lake star”. These
compelling features have been described qualitatively a number of times (e.g. [5, 6, 7]) but work on the formation
process itself has been solely heuristic. Knight [5] outlines a number of the physical ideas relevant to the process, but
does not translate them into a predictive framework to model field observations. Knight’s main idea is that locations
FIG. 1: Typical lake star patterns. The branched arms are approximately 1 m in length. Quonnipaug Lake, Guilford,
Connecticut, 8 March, 2006.
∗Electronic address: vtsai@fas.harvard.edu
†Electronic address: john.wettlaufer@yale.edu
http://arxiv.org/abs/0704.1504v1
mailto:vtsai@fas.harvard.edu
mailto:john.wettlaufer@yale.edu
FIG. 2: Schematic of the geometry of the model. The perspective is looking down on a nascent star. The equations (refer to
text for numbering) are shown in the domains of the system where they are applicable.
with faster flow rates melt preferentially, leading to even faster flow rates and therefore to an instability that results
in fingers. This idea has features that resemble those of many other instabilities such as, for example, those observed
during the growth of binary alloys [8], in flow of water through a rigid hot porous media [9], or in more complex
geomorphological settings [10], and we structure our model accordingly.
Katsaros [6] and Woodcock [7] attribute the holes from which the stars emanate and the patterns themselves to
thermal convection patterns within the lake, but do not measure or calculate their nature. However, often the holes
do not exhibit a characteristic distance between them but rather form from protrusions (e.g. sticks that poke through
the ice surface) [5] and stars follow thereby ruling out a convective mechanism as being necessary to explain the
phenomena. The paucity of literature on this topic provides little more than speculation regarding the puncturing
mechanism but lake stars are observed in all of these circumstances. Therefore, while hole formation is necessary for
lake star formation, its origin does not control the mechanism of pattern formation, which is the focus of the present
work.
II. THEORY
The water level in the hole is higher than that in the wet snow–slush–layer [5] and hence we treat this warm water
[11] region as having a constant height above the ice or equivalently a constant pressure head, which drives flow
of water through the slush layer, which we treat as a Darcy flow of water at 0◦C. We model the temperature field
within the liquid region with an advection-diffusion equation and impose an appropriate (Stefan) condition for energy
conservation at the water-slush interface. The water is everywhere incompressible. Finally, the model is closed with
an outer boundary condition at which the pressure head is assumed known.
Although we lack in-situ pressure measurements, circular water-saturated regions (a few meters in radius) are
observed around the lake stars. Hence, we assume that the differential pressure head falls to zero somewhere in the
vicinity of this circular boundary. The actual boundary at which the differential pressure head is zero is not likely
to be completely uniform (as in Figure 4 of Knight [5]) but treating it as uniform is a good approximation in the
linear regime of our analysis. Finally, we treat the flow as two-dimensional. Thus, although the water in direct
contact with ice must be at 0◦C, we consider the depth-averaged temperature, which is above freezing. Additionally,
the decreasing pressure head in the radial direction must be accompanied by a corresponding drop in water level.
Therefore, although the driving force is more accurately described as deriving from an axisymmetric gravity current,
the front whose stability we assess is controlled by the same essential physical processes that we model herein. Our
analysis could be extended to account for these three-dimensional effects.
The system is characterized by the temperature T , a Darcy fluid velocity u, pressure p, and an evolving liquid-
slush interface a. The liquid properties are κ (thermal diffusivity), CP (specific heat at constant pressure) and µ
(dynamic viscosity) and the slush properties are Π (permeability), ξ (solid fraction) and L (latent heat). We non-
dimensionalized the equations of motion by scaling the length, temperature, pressure and velocity with r0, T0, p0,
and Πp0
, respectively. Thus, our model consists of the following system of dimensionless equations:
+ u · ∇θ = ǫ∇2θ r0 < r < a(φ, t), (1)
θ = 0 a(φ, t) < r < 1, (2)
p = 1 ri < r < a(φ, t), (3)
∇2p = 0 a(φ, t) < r < 1, (4)
∇ · u = 0 ri < r < a(φ, t), (5)
= u |a+ r = a(φ, t), (6)
u = −∇p a(φ, t) < r < 1, (7)
with boundary conditions
ȧ = −
∇θ r = a(φ, t), (8)
1 r = ri
0 r = a(φ, t)
0 r = 1
, (9)
1 r = ri
1 r = a(φ, t)
0 r = 1
, (10)
where (1) describes the temperature evolution in the liquid, (4) and (5) describe mass conservation with a Darcy flow
(7) in the slush, (8) is the Stefan condition, and (9) and (10) are the temperature and pressure boundary conditions,
respectively (see Figure 2). Note that (3) and (5) can both be satisfied since the liquid region has an effectively infinite
permeability.
The dimensionless parameters ǫ and S of the system are given by
, and S ≡
, (11)
which describe an inverse Peclet number and a Stefan number respectively. Because the liquid must be less than
or equal to 4◦C, we make the conservative estimates that T0 < 4
◦C, ξ > 0.3, and use the fact that L/CP ≈ 80
from which we see S > 6 ≫ 1. Using κ ≈ 10−7m2s−1, and the field observations of Knight [5] to constrain u0
(1cm/hr < u0 < 10cm/hr) and r0 (0.3m < r0 < 3m), we find that ǫ < 0.1 ≪ 1. We therefore employ the quasi-
stationary (S ≫ 1) and large Peclet number (ǫ ≪ 1) approximations, and hence equations (1) - (10) are easily
solved for a purely radial flow with cylindrical symmetry (no φ dependence) and circular liquid-slush interface. This
(boundary layer) solution is
u = ur̂ = −
ln(a0)
r̂ ri < r < 1, (12)
ln(r)
ln(a0)
r > a0, (13)
θ0 = 1−
(−1/ ln(a0)+2ǫ)
r < a0, (14)
0 0.5 1 1.5 2 2.5 3
σ       
/ln(a
/ln2(a
 = 0.06
 = 0.5
approx
FIG. 3: Stability curve: Non-dimensional growth rate σ versus non-dimensional wavenumber k′. Scales for the axes are given
at the upper left (σ axis) and the lower right corners (k′ axis). σ is plotted for the range of plausible a0 (dot-dashed blue and
dashed red curves) and for the approximation (18) (solid green curve).
Sa0ȧ0
−1/ ln(a0) + 2ǫ
= 1, (15)
where equation (15) has an approximate implicit solution for a0 given by
a20 ln(a0) =
. (16)
We perform a linear stability analysis around this quasi-steady cylindrically symmetrical flow. Proceeding in the
usual way, we allow for scaled perturbations in θ and a with scaled wavenumber k′ = ǫk, non-dimensional growth rate
σ, and amplitudes f(r) and g respectively. Keeping only terms linear in ǫ, 1/S and g, we solve (4) subject to (10),
substitute into (6) and satisfy (5) and (1). This gives the non-dimensional growth rate (σ) as a function of scaled
wave number (k′):
2a0 ln
2(a0)S
1 + 4k′2 ln2(a0)− 1
−k′ ln(a0)
. (17)
Equation (17) can be approximated in 0 ≤ x . 1 as
ln2(a0)S
x(1− x), (18)
where x ≡ −k′ ln(a0)/a0.
The stability curve (17) and the approximation (18) are plotted in Figure 3. The essential features of (17) are a
maximum in the range 0 < k′ < a0/ ln(a0), zero growth rate at k
′ = a0/ ln(a0) and a linear increase in stability with
k′ for large k′. The long-wavelength cut-off is typical of systems with a Peclet number, here with the added effect of
latent heat embodied in the Stefan number. This demonstrates the competition between the advection and diffusion
of heat and momentum (in a harmonic pressure field); the former driving the instability and the latter limiting its
extent. The maximum growth rate occurs at approximately
k′max ≈
−2 ln(a0)
, (19)
with (non-dimensional) growth rate
σmax ≈
4S ln
. (20)
FIG. 4: Schematic showing r0, r
, rLS and rℓ.
Translating (19) and (20) back into dimensional quantities, we find that the most unstable mode has angular size
given by
φdegrees =
720◦κ
, (21)
and has growth rate given by
σdim =
4Sr0 ln
2(r0/a0)
. (22)
III. EXTRACTING INFORMATION FROM FIELD OBSERVATIONS
Field observations of lake stars cannot be controlled. A reasonable estmate for r0 is the radius of the wetted
(snow) region around the lake stars, and observations [5, 6, 7] bound the value as 1.5m . r0 . 4m. This is simply
because if there were significant excess pressure at this point then the wetting front would have advanced further.
However, it is also possible that the effective value of r0, say r
0 , is less than this either because the wetted radius
is smaller earlier in the star formation process or because the ambient pressure level is reached at smaller radii.
Here, we take a0 to be the radius of the roughly circular liquid-filled region at the center of the lake star (rℓ) as
the best approximation during the initial stages of star formation (see Figure 4). Field observations show that
0.1m . rℓ . 0.5m, [5, 6, 7] and hence 0.07 . rℓ/r0 . 0.15. We note that equations (21) and (22) are more sensitive
to a0/r0 than a0 or r0 independently[12]. With this interpretation of r0 we find a reasonable estimate of u0 as
1.4 · 10−5m/s . u0 . 2.8 · 10
−5m/s. Using these parameter values, the most unstable mode should have wavelength
between 8◦ and 130◦. Letting the number of branches be N = 360◦/φdeg, then 3 < N < 45 and we clearly encompass
the observed values for lake stars (4 < N < 15), but note that values (N > 15) are never seen in the field.
Despite the dearth of field observations, many qualitative features embolden our interpretation. For example, the
stars with larger values of a0/r0 have a larger number of branches. Moreover, for any value of a0/r0, our analysis
predicts an increase in N with r0 and u0. Indeed, u0 increases with p0 (higher water height within the slush layer)
and Π (less well-packed snow). Therefore, we ascribe some of the variability among field observations to variations
in these quantities (which have not been measured in the field) and the remainder to nonlinear effects. Because the
dendritic arms are observed long after onset and are far from small perturbations to a radially symmetric pattern,
as one might see in the initial stages of the Saffman-Taylor instability, the process involves non-linear cooperative
phenomena. Hence, our model should only approximately agree with observations. Although a rigorous non-linear
analysis of the long term star evolution process (e.g. [3]) may more closely mirror field observations, the present state
of the latter does not warrant that level of detail. Instead, we examine the model physics through simple proof of
concept experimentation described presently.
FIG. 5: Typical experimental run where small- scale fingers are present. For scale, the nozzle head has diameter of 5 mm.
FIG. 6: Typical run where channels form. This picture is taken from the underside. Note: part of the slush broke off when it
was flipped to image it. The ruler scale is in cm.
IV. DEMONSTRATING LAKE STARS IN THE LABORATORY
A 30 cm diameter circular plate is maintained below freezing (≈ −0.5◦C), and on top of this we place a 0.5 to 1 cm
deep layer of slush through which we flow 1◦C water. Given the technical difficulties associated with its production,
the grain size, and hence the permeability, of the slush layer, is not a controlled variable. This fact influences our
results quantitatively. In fourteen runs we varied the initial size of the water-filled central hole (a0), that of the
circular slush layer (r0), and the flow rate (Q), which determines u0. The flow rate is adjusted manually so that the
water level (h0) in the central hole remains constant [13]. Fingering is observed in every experimental run and hence
we conclude that fingers are a robust feature of the system. Two distinct types of fingering are observed: small-scale
fingering (see Figure 5) that forms early in an experimental run, and larger channel-like fingers (see Figure 6) that are
ubiquitous at later times and often extend from the central hole to the outer edge of the slush. Since the channel-like
fingers provide a direct path for water to flow, effectively shorting Darcy flow within the slush, their subsequent
dynamics are not directly analogous to those in natural lake stars. However, in all runs, the initial small-scale fingers
have the characteristics of lake stars and hence we focus upon them. We note that because the larger channel-like
fingers emerge out of small-scale fingers, they likely represent the non-linear growth of the linear modes of instability,
a topic left for future study. Finally, we measure the distance between fingers (df ), so that for each experiment we
can calculate u0 = Q/(2πr0h0), φcalc ≡ φdegrees, from equation (21), and φobs = 180
◦df/(πa0), and we can thereby
compare experiment, theory and field observations.
In Figure 7 we plot φobs versus φcalc for the various field observations for which we have estimates of parameters,
the laboratory experiments described above, and the model [equation (21)]. There is a large amount of scatter in both
0 100 200 300 400 500
 (deg)
FIG. 7: Comparison of theory, experiment and field observations. Circles are field observations (cyan = best constrained field
observation, black = range of plausible field observations), triangles are experimental results (blue upward-pointing triangles
were unambiguous; red left-pointing triangles have channels but show no clear small-scale fingers, so channel spacing is taken
for df ; green right-pointing triangles were compromised by the quality of the images). Errors are approximately 0.3 cm, 0.5
cm, 2 mm, 5 ml/min and 0.2 cm (respectively) for the five measured quantities. All experimental results thus have error bars
of at least a factor of two in the x-coordinate and 30% in the y-coordinate. Typical error bars are shown on one measurement.
The solid red line is the theoretical prediction; the dotted green line is the best fit line to the blue triangles.
the experimental and observational data and the data does not lie on the one-to-one curve predicted by the model.
However, the experiments are meant to demonstrate the features of the model predictions, and the results have the
correct qualitative trend (having a best-fit slope of 0.34). We also attempt to find trends in the experimental data
not represented by the model by comparing y ≡ φobs/φcalc vs. various combinations of control parameters (≡ x)
including r0, a0, r0/a0, r0u0, r0/a0 ln(r0/a0) and ln(r0/a0)/(a0u0). For all plots of y vs. x, our model predicts a zero
slope (and y-intercept of 1). A non-random dependence of y on x would point to failure of some part of our model.
Thus, to test the validity of our model, we perform significance tests on all non-flagged data with the null hypothesis
being a non-zero slope. In all cases, the null hypothesis is accepted (not rejected) at the 95% confidence level. Thus,
although the agreement is far from perfect, the simple model captures all of the significant trends in the experimental
data.
V. CONCLUSIONS
By generalizing and quantifying the heuristic ideas of Knight [5], we have constructed a theory that is able to
explain the radiating finger-like patterns on lake ice that we call lake stars. The model yields a prediction for the
wavelength of the most unstable mode as a function of various physical parameters that agrees with field observations.
Proof of concept experiments revealed the robustness of the fingering pattern, and to leading order the results also
agree with the model. There is substantial scatter in the data, and the overall comparison between field observations,
model and experiment demonstrates the need for a comprehensive measurement program and a fully nonlinear theory
which will yield better quantitative comparisons. However, the general predictions of our theory capture the leading
order features of the system.
VI. ACKNOWLEDGEMENTS
We thank K. Bradley and J. A. Whitehead for laboratory and facilities support. This research, which began at the
Geophysical Fluid Dynamics summer program at the Woods Hole Oceanographic Institution, was partially funded by
National Science Foundation (NSF) grant OCE0325296, NSF Graduate Fellowship (VCT), NSF grant OPP0440841
(JSW), and Department of Energy grant DE-FG02-05ER15741 (JSW).
[1] T.Y. Hou, J.S. Lowengrub, and M.J. Shelley, J. Comp. Phys. 169, 302 (2001).
[2] M.P. Brenner, L.S. Levitov, and E.O. Budrene, Biophys. J. 74, 1677 (1998); H. Levine and E. Ben-Jacob, Phys. Biol. 1,
14 (2004).
[3] M.C. Cross and P.C. Hohenberg, Rev. Mod. Phys. 65, 851 (1993); I.S. Aranson and L.S. Tsimring, Rev. Mod. Phys. 78,
64 (2006)
[4] N. Goldenfeld, P.Y. Chan and J. Veysey, Phys. Rev. Lett. 96, 254501 (2006); M.B. Short, J.C. Baygents, and R.E.
Goldstein, Phys. Fluids 18, 083101 (2006).
[5] C.A. Knight, in Structure and dynamics of partially solidified systems, ed. D.E. Loper, (Martinus Nijhoff, Dordrecht, 1987),
pp. 453-465.
[6] K.B. Katsaros, Bull. Amer. Meteor. Soc. 64, 277 (1983).
[7] A.H. Woodcock, Limnol. Oceanogr. 10 R290 (1965).
[8] M.G. Worster, Annu. Rev. Fl. Mech. 29, 91 (1997).
[9] S. Fitzgerald and A.W. Woods, Nature 367, 450 (1994).
[10] N. Schorghofer, B. Jensen, A. Kudrolli and D.H. Rothman, J. Fluid Mech. 503, 357 (2004).
[11] A finite body of fresh water cooled from above will have a maximum below ice temperature below of 4 ◦C.
[12] For the later stages of growth, clearly in the nonlinear regime not treated presently, a0 may also be interpreted as the
radius of the lake star (rLS). Field observations show 1m . rLS . 2m [5, 6, 7] and hence 0.3 . rLS/r0 . 0.6.
[13] In many of the runs, we begin the experiment without the central hole. In practice, however, the first few drops of warm
water create a circular hole with radius one to three times the radius of the water nozzle (0.5cm < a0 < 1.0cm). It is
significantly more difficult to prepare a uniform permeability sample with a circular hole initially present; these runs are
therefore more difficult to interpret.
	Introduction
	Theory
	Extracting information from field observations
	Demonstrating Lake Stars in the Laboratory
	Conclusions
	Acknowledgements
	References
