\documentclass[11pt]{article}

\usepackage{graphicx}
\usepackage{amsmath,amssymb}
\usepackage[body={6in,8.5in},top=1.7in]{geometry}
\usepackage[T1]{fontenc}
\usepackage{aeguill}
\usepackage{multirow}
\usepackage{lastpage,xr}
\externaldocument[pr-]{../../pref06}
\bibliographystyle{siam}



\newcommand{\dd}{\mathrm{d}}
\newcommand{\totd}[2]{\frac{\dd #1}{\dd #2}}
\newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\beqn}{\begin{eqnarray}}
\newcommand{\eqn}{\end{eqnarray}}

\title{GFD 2006 Lecture 1: Introduction to Ice}
\author{Grae Worster; notes by Rachel Zammett and Devin Conroy}
\date{\today}


\begin{document}
\maketitle

\setcounter{page}{\pageref{pr-LastPage}}
\addtocounter{page}{1}

\section{Introduction}

Our aim in this course is to understand some of the processes
associated with ice in the natural environment. Figure \ref{cryo}
shows the location of some of Earth's ice during the
northern winter. These ice deposits may be categorized as sea ice,
ice sheets and shelves, and permafrost.

%Ice is relevant to Earth and its global circulations; it can flow,
%and there is convection associated with ice production.  Our aim
%in this course is to understand the interaction between a
%convective fluid flow and solidification.  Convection may be
%forced, by ocean circulation for example, or occur naturally, as
%is the case with buoyancy forces.
%Figure 1 shows examples of the three main types of
%ice that we will be looking at in this lecture course: sea ice,
%ice sheets and shelves, and permafrost.


\begin{figure}[h]
\begin{center}
 \includegraphics[height=60mm]{cryo.eps}
 \caption{Satellite image showing the ice cover in the northern hemisphere during northern winter, showing sea ice lying in the Arctic basin, the permanent ice sheet over Greenland and permafrost in the exposed land surface.}
\end{center}\label{cryo}
\end{figure}

\section{Ice sheets}

Firstly, figure \ref{cryo} shows the ice sheet that covers approximately
80\% of Greenland. This is about $10^5$ years
old and reaches depths of 2--3 kilometers.  On large scales, ice can be treated as a highly viscous, non-Newtonian fluid that can flow because it is a polycrystalline solid and contains a
percentage of unfrozen water (figure \ref{grains}).  Looking on a scale of about $100\mu
m$, we can see the ice grain junctions and the veins which lie
between them. The liquid water contained in
the veins between the ice crystals lubricates the flow, allowing
the ice to flow more easily.  This water can also transport
dissolved impurities, which will therefore move relative to the
ice crystals; this is important when analyzing ice cores, for
example.

\begin{figure}
\begin{center}
 \includegraphics[height=60mm]{grains.eps}
 \caption{Image of the intersection of four ice grains.   Between these grains lie the veins containing liquid water
 and dissolved impurities.  The scale bar on this picture is 100 $\mu$m.\label{grains}}
\end{center}
\end{figure}


Figure \ref{grains} also shows that there is a curvature to the
solid--liquid interface which is associated with the surface
energy of the phase boundary. We will see later that this surface
energy sets the scales for morphological instabilities of the
solid--liquid interface, such as those seen in snowflakes (figure
\ref{flake1}).   

\begin{figure}
\begin{center}
\includegraphics[scale=1,draft=false]{rachel2.eps}
 \caption{(left) Satellite image of the Larsen B ice shelf
on the coast of Antarctica.  Near the edge of the ice shelf, it is
possible to see the icebergs formed by calving. (right) An example
of an iceberg formed by calving at the edge of an ice sheet or ice
shelf. Here the vertical face is 30 m above the surface of the
ocean, meaning that approximately 300 m of ice lie below the
surface. \label{shelf}}
\end{center}
\end{figure}


The grounded ice cap
flows slowly towards the coast, sometimes flowing into floating
ice shelves, which ultimately break up to form icebergs. Projects such
as the Greenland Ice Core Project (GRIP) have obtained deep ice
cores from near Greenland's summit. Analyzing the properties of
the ice cores, such as oxygen isotope ratios, allow inferences
about the ancient climate to be drawn.

%Comparison of oxygen isotope records from the GISP2 and GRIP
%Greenland ice cores
%P. M. Grootes& ast;, M. Stuiver& ast;, J. W. C. White& dagger;, S.
%Johnsen& Dagger;& sect; & J. Jouzel& par



In figure \ref{shelf} we see the flow from the grounded ice sheet
to a floating ice shelf (Larsen B) in Antarctica.  At the edge of
the ice shelf we see the calving of icebergs;  this is responsible
for approximately 80\% of the mass lost from Antarctica.   The
icebergs are composed predominantly of freshwater ice, as the ice which comprises the ice sheets
first fell as snow.  Owing to the density difference between water
and ice, approximately 90\% of the volume of an iceberg is below
the surface of the ocean.

When these icebergs come into contact with the warm, salty ocean
they ablate, providing a freshwater flux to the ocean. This is
important as the production of deep ocean waters is sensitive to
changes in the freshwater budget.




\section{Sea ice}

Secondly, there is the sea ice which fills the Arctic basin and is
formed by direct freezing of the ocean.  It is typically 1--3 m
thick and less than 10 years old;  in its first year, sea ice
typically grows to a depth of 1 m.  This relative youth (in
comparison to ice sheets or glaciers, for example) is caused by
the movement of sea ice by polar winds and ocean currents to
warmer waters, where it melts.

%The main current responsible for this is the Beaufort Gyre, a
%clockwise ocean circulation which occurs north of Canada and
%Alaska, making one complete rotation every 4 years. The Transpolar
%Drift, which carries water and ice from east to west in Russian
%Siberia, then carries the sea ice to the warmer waters of the
%north Atlantic ocean.



%\begin{figure}[h!]
%\begin{minipage}[c]{.5\textwidth}
% \includegraphics[scale=0.75,draft=false]{icesheet.eps}
%  \label{shelf}
%\end{minipage}
%\begin{minipage}[c]{.5\textwidth}
% \includegraphics[scale=0.55,draft=false]{berg1.eps}
%  \label{berg1}
%\end{minipage}


\begin{figure}
\begin{minipage}[c]{.4\textwidth}
\includegraphics[height=50mm]{flake1.eps}
\end{minipage}
\begin{minipage}[c]{.4\textwidth}
\includegraphics[height=50mm]{surfacechannels.eps}
\end{minipage}
\caption{(left) Picture of a snowflake.  Here the smallest scale
at which instabilities occur is comparable to the radius of
the tip of one of the needles. \label{flake1}(right) Horizontal cross section  of sea
ice, showing both ice platelets and brine channels. The ice platelets are typically less than 1 mm wide and form
a porous matrix, which allows convection and the erosion of such
channels by the rejection of salt. These channels have a diameter
of a few millimetres.\label{platelets} }
\end{figure}

%\begin{figure}
%\begin{center}
% \includegraphics[height=60mm]{platelet.eps}
% \caption{Image of sea ice, showing both ice platelets and the brine between them. The ice platelets are typically less than 1 mm wide and form a porous matrix, which allows convection and the erosion of channels by the rejection of salt. These channels have a diameter of $O($mm$)$.\label{platelets}}
%\end{center}
%\end{figure}

Many of the structures and processes observed in sea ice develop
because the thermal diffusivity of heat is much larger than the
diffusivity of salt. In this course we shall see that sea ice can
be considered as an inhomogeneous porous medium.
 While
sheet ice only contains water in the veins between ice crystals,
sea ice has a much higher porosity (approximately 10 \% in old ice
and up to 40 \% in new ice).  The porous nature of sea ice means
that it can also be modified by internal convection.

We shall consider sea ice to be a mushy layer, which is a two-phase reactive porous medium.  We see in figure \ref{platelets}
that it is not macroscopically solid;  instead, it is composed of
ice platelets with salty brine between them.  The platelets which
form are composed of pure ice crystals, as the crystals reject the
salt contained in the ocean water. Some of this rejected salt
convects into the ocean below the sea ice, and the rest remains
between the crystals.

\begin{figure}
%\begin{minipage}[c]{.3\textwidth}
%\includegraphics[height=40mm]{threecmice.eps}
%\caption{Image of sea ice growing in the laboratory.  At this time, the layer of ice is 3 cm thick, and it is possible to see some convection occurring below it.  \label{threecmice}}
%\end{minipage}
%\begin{minipage}[c]{.4\textwidth}
%\includegraphics[height=50mm]{13cmice.eps}
%\caption{Image of sea ice growing in the laboratory where the layer of ice is now 13 cm thick, and it is possible to see the salty convective plumes
 %below it.  It is also possible to see the `pinching' instability at the base of these plumes. \label{13cmice}}
%\end{minipage}
\includegraphics[scale=1,draft=false]{rachel1.eps}
\caption{(left) Image  of sea ice growing in the laboratory. At
this time, the layer of ice (the dark upper region) is 3 cm thick,
and it is possible to see some convection occurring below it.
(right) Image of sea ice growing in the laboratory where the layer
of ice is now 13 cm thick, and it is possible to see the salty
convective plumes below it.  It is also possible to see the
`pinching' instability at the base of these plumes.
 Note that the scales of the images are different.}\label{13cmice}\label{3cmice}
\end{figure}

\begin{figure}
\begin{minipage}[c]{.4\textwidth}
\includegraphics[width=1\textwidth,trim= 0 0 0 0,clip]{frost1.eps}
\end{minipage}
\qquad
\begin{minipage}[c]{.5\textwidth}
\includegraphics[width=1\textwidth,trim= 0 0 0 0,clip]{heave.eps}
\end{minipage}
\caption{(left) An example of erosion caused to a rockery by winter frost. (right) Stone circles as an example of differential frost heave.}\label{frost1}\label{heave}
\end{figure}

\begin{figure}
\begin{minipage}[c]{.4\textwidth}
\includegraphics[width=.9\textwidth]{needles.eps}
\end{minipage}
\begin{minipage}[c]{.4\textwidth}
\includegraphics[width=.9\textwidth,trim= 160 160 160 160,clip]{colmn.eps}
\end{minipage}
\caption{(left) Ice needles protruding from soil.
\label{needles}(right) Photograph of a column of water-saturated
soil cooled at the top (Taber, 1929).  The black regions are ice
lenses, which contain no soil; between these are regions of
partially frozen soil.  There may also be ice between the soil
particles below the lowest lens.}
\end{figure}

This convection is also seen in the laboratory.  Figure
\ref{3cmice} shows shadowgraph pictures of sea ice growing in a
laboratory. In figure \ref{13cmice} (left) when the sea ice is only 3 cm
thick, it is possible to see some convection occurring in the salt
water below it, but it is small scale and has no obvious
structure. However, when the ice has grown to a thickness of about 13 cm
(figure \ref{13cmice}), it is possible to see strong convective plumes in
the water below. These have a high salt content and therefore
deliver a large flux of salt to the water below. We shall see that
there is a critical ice thickness at which such plumes occur; this
criterion for the onset of convection is determined by a form of
Rayleigh number.

\section{Permafrost}

The final type of ice we shall consider is permafrost, or
permanently frozen ground (defined as remaining below 0$^{\circ}$C
for more than two years). It occurs both on land and beneath
offshore Arctic continental shelves, and its thickness ranges from
less than one meter to greater than 1 kilometer. Permafrost
underlies about 15\% of the exposed land surface in the Northern
Hemisphere and causes deformation of the ground;  we shall be
looking at this and the associated flows.

Figure \ref{frost1} shows the effects of ice damage to rocks and
buildings if they are eroded by winter frost.  The ground may also
`heave', i.\,e.\,rise upwards due to water being pulled up from
the unfrozen ground below. Differential frost heave may form
patterned ground, such as hummocks and the stone circles seen in
figure \ref{heave}. Underlying this is the force of separation
between ice and other materials: in this context, we will consider
the other materials to be silicates. We will consider how ice
pushes on another material forming, for example, the ice needles
seen in figure \ref{needles}.

There are still some puzzles remaining.  Figure \ref{needles}
shows a laboratory experiment by Taber where a column of
water-saturated soil was frozen by cooling at its top.  It might
be expected that a freezing front which moves downwards is
observed, but instead a sequence of layers of alternating pure ice
and partially frozen soil forms.


\section{Student Problem}

If two identical ice cubes are placed in glasses of water and
whisky, where the liquids are at the same temperature, it is
observed that the ice cube in the whisky melts more quickly than
that in the water.  Why? (Hint: It is not because the melting
point of ice is lower in whisky than in water.)

\subsection*{Answer}

Initially when the ice cube is placed into a glass of whisky at
room temperature the ice melts, forming a layer of cold freshwater
adjacent to the phase boundary. Since water is denser than alcohol
and the the melted water is colder than the whisky, a plume forms
that convects the cool fresh water downwards and brings warmer
fluid with a higher alcohol concentration upwards.  This
convective mixing of the liquid below the ice cube supplies a heat
flux at the phase boundary; this flux is stronger than the
diffusive heat flux in the absence of convection.

%\begin{figure}
%\begin{center}
% \includegraphics[height=60mm]{berg2.jpg}
% \caption{Ice needles protruding from a log \lab{needles}}
%\end{center}
%\end{figure}



%The salt/ice ratio is important, as the removal of salt changes
%the density of the ice and convection is possible???  (figure
%\ref{rock}).

%\begin{figure}
%\begin{center}
% \includegraphics[height=60mm]{rock.jpg}
% \caption{\lab{rock}}
%\end{center}
%\end{figure}

\section{Stefan Condition}
The distinguishing feature of solidification or melting is the
evolution of a phase boundary which separates solid and liquid.
The speed of this interface can be determined by 
energy conservation, as illustrated in figure \ref{e:stefan},
which relates the rate of energy absorption or release to the
difference in heat fluxes across this boundary.  This 
\begin{figure}
\center{\includegraphics[width=.5\textwidth,trim= 0 0 0 0,
clip]{stefan.eps}} \caption{Illustration showing the control volume taken around the phase boundary and the energy fluxes into and out of it.}
\label{e:stefan}
\end{figure}
is formulated mathematically as follows by considering a
control volume around the phase boundary \beqn \mathbf{q_s} \cdot
\mathbf{n}-\mathbf{q_l} \cdot \mathbf{n} -\rho V_n H_s+\rho V_n
H_\ell=0 \Rightarrow \rho L V_n=\mathbf{n}\cdot
\left(\mathbf{q_l}-\mathbf{q_s}\right). \eqn Here
$\mathbf{q}=-k\nabla T$ is the heat flux from Fourier's law, $k$
is thermal conductivity, $\mathbf{n}$ is the unit normal vector
pointing from solid to liquid, $\rho$ is the density (assumed to
be the same in each phase), $V_n$ the interface velocity, $H$ is
the enthalpy, $L=H_\ell-H_s$ is the latent heat and subscripts $s$
and $\ell$ denote solid and liquid respectively. We assume that
the phase boundary is in equilibrium, implying that the
temperature is constant on either side of the interface. This
equation is known as the Stefan condition, attributed to Stefan
in 1891.


\section{Problem 1}

We consider a problem posed by Stefan in 1891, where solid ice is
growing into relatively warm ($T_m>T_B$) water from a cooled
boundary at $z=0$ (figure \ref{p1}). We assume that the liquid
portion is at the melting temperature $T_m$ initially and
therefore remains at this temperature. The governing equation is given by the
thermal diffusion equation \beqn \rho c_p\pd{T}{t}=\pd{}{z}\left(k
\pd{T}{z}\right)\Rightarrow \pd{T}{t}=\kappa\frac{\partial^2
T}{\partial t^2},\label{diff} \eqn where the conductivity
$k=\rho c_p \kappa$ is assumed to be constant and the thermal
diffusivity is represented by $\kappa$. The boundary conditions
for this equation are then \beqn T(t,z=0)=T_B, \qquad
T\left(t,z=a(t)\right)=T_m,\label{bc} \eqn where $z=a(t)$ denotes
the interface position and the unknown interface velocity is
determined by the Stefan condition \beqn \rho L
\totd{a}{t}=k\left.\pd{T}{z}\right|_{z=a(t)}. \eqn 

\subsection{Solution}
This problem can be solved using a similarity solution, as there is no intrinsic length scale in the problem. We can determine the form of this similarity variable using a scaling analysis and show that a fixed length scale cannot be formed. Let us define the appropriate scales, $\Delta T=T_m-T_B$, $D$ and $\tau$ for the temperature, length of the domain and time respectively. From the diffusion equation (\ref{diff}) and the Stefan condition (\ref{e:stefan}) we obtain the following

\beqn
\frac{\Delta T}{\tau} \sim \kappa \frac{\Delta T}{D^2}\Rightarrow D \sim \sqrt{\kappa\tau},\label{qq}\\
\rho L\frac{D}{\tau}\sim k\frac{\Delta T}{D}= \rho c_p \kappa
\frac{\Delta}{D}T\Rightarrow D\sim \sqrt{\kappa
\tau}S^{-1/2},\label{qw} \eqn where $S=L/(c_p \Delta T)$
is the Stefan number. Since the relationships between $D$ and $\tau$ in (\ref{qq}) and
(\ref{qw}) are the same, there is no intrinsic length scale, and
a similarity solution is possible. In addition there therefore is no
time scale so we choose $\tau=t$, the actual time, in which case $D \sim \sqrt {\kappa t}$. 

\begin{figure}
\center{\includegraphics[width=.5\textwidth,trim= 0 0 0 0,
clip]{p1.eps}} \caption{Growth of a planar solid into a liquid, maintained at the melting temperature $T_m$ from a cooled boundary maintained at temperature $T_B$. The position of the interface as a function of time is given by $a(t)$.}
\label{p1}
\end{figure}

We introduce the dimensionless variable $f$ such that
\beqn
T-T_B=\Delta T f\left(\frac{z}{D},\frac{t}{\tau}\right)=\Delta T\,f\left(\frac{z}{\sqrt{\kappa t}},1\right)=\Delta T\,f(\eta).
\eqn
We choose $\eta={z}/{2\sqrt{\kappa\,t}}$ for mathematical convenience, as the similarity variable may be multiplied by an arbitrary constant. In addition we know from scale analysis, using the diffusion equation, that length and therefore interface position can be assumed to have the form
\beqn
a=2\mu\sqrt{\kappa t},
\eqn
where the parameter $\mu$ must be determined as part of the solution. Rewriting the model in terms of the similarity variable, we arrive at the final set of non-dimensional equations
\beqn
f^{\prime \prime}=-2\eta\,f^{\prime}\label{f:sim},\\
f(\eta=0)=0,\label{1bc}\\
f(\eta=\mu)=1,\label{2bc}\\
2\,S\,\mu=f^{\prime}(\eta=\mu).\label{sbc} \eqn The solution to
equation (\ref{f:sim}) is determined using an integrating factor and determining 
\beqn
f^{\prime}=c_1\,e^{-\eta^2}\Rightarrow f=c_1 \int_0^\eta
e^{-y^2}\,dy+c_2\Rightarrow f=\hat{c_1}\mathrm{erf}(\eta) +c_2.
\eqn The error function, erf$(\eta)$, is defined by \beqn
\mathrm{erf}(x)=\frac{2}{\sqrt{\pi}}\int_0^x e^{-u^2}\,du, \eqn
with the following properties \beqn
\mathrm{erf}(0)=0,\\
\mathrm{erf}(\infty)=1. \eqn Boundary condition (\ref{1bc})
implies that $c_2=0$ and boundary condition (\ref{2bc}) gives
$\hat{c_1}=1/\mathrm{erf}(\mu)$. Finally the Stefan condition
(\ref{sbc}) gives us an parameter equation to be solved for the
parameter $\mu$. From equation (\ref{sbc}) we then have \beqn
\frac{1}{S}=\sqrt{\pi}\,\mu\,\mathrm{erf}(\mu)\,\mathrm{e}^{\mu^2}=F(\mu).\label{eig}
\eqn In figure \ref{stef} we plot the parameter $\mu$ as a function of  Stefan number. We
see that the growth speed increases as the Stefan number
decreases, which corresponds to increasing the driving temperature
difference or decreasing the amount of energy required to melt a
unit mass of solid. We should note that the interface position is
given by $a \propto \sqrt{t}$ and the interface velocity by
$\dot{a}\propto 1/\sqrt{t}$ so that the growth rate of a solid decreases with time.
\begin{figure}
\center{\includegraphics[width=.5\textwidth,trim= 0 0 0 0,clip]{fig1.10.eps}} \caption{(a) Solution to equation (\ref{eig}) for the eigenvalue $\mu$ as a function of the Stefan number, $S$. (b) Solution to equation (7) of Lecture 2.}
\label{stef}
\end{figure}


\subsection{Quasi Stationary Approximation}
For large Stefan numbers we have relatively small sensible heat compared to latent heat and the growth rate will be slow compared to the thermal diffusion rate. In this case, the temperature field will evolve more rapidly than the boundary position and we have a quasi-steady state regime for the the diffusion equation. This implies a linear profile of the solid temperature that is slowly decreasing in slope as the boundary moves. Integrating Laplace's equation twice and applying boundary conditions given in (\ref{bc}) we obtain the linear solution
\beqn
T=T_B+(T_m-T_B)\frac{z}{a}.
\eqn
Substituting this solution into the Stefan condition (\ref{sbc}) we arrive at the following expression for the interface position $a$
\beqn
a\totd{a}{t}=\frac{1}{S},\qquad a(0)=0 \qquad \Rightarrow a=\sqrt{\frac{2}{S}}\sqrt{\kappa t}.
\eqn


\section{Student Problem}

\subsection*{Question}
Solve the Stefan problem given in problem 1 for the case $\rho_s\neq\,\rho_{\ell}$.
\subsection*{Answer}
When the density between the solid and liquid differ by an appreciable amount there will be a normal velocity to the phase change surface due to the expansion or contraction of the liquid as it solidifies. To formulate this effect mathematically we draw a control volume around the moving surface and conserve mass and energy as follow
\beqn
\totd{}{t}\int_V \rho\,dV=\int_S  \rho_s\left(\dot{a}-V_s\right) dS_s-\int_S  \rho_{\ell}\left(\dot{a}-V_{\ell}\right) dS_{\ell},\\
\totd{}{t}\int_V \rho H\,dV=\int_S \rho_sH_s\left(\dot{a}-V_s\right) +\mathbf {n} \cdot \mathbf{q_s} \,\,dS_s-\int_S \rho_{\ell}H_{\ell}\left(\dot{a}-V_{\ell}\right)+\mathbf {n}\cdot \mathbf{q_{\ell}} \,\,dS_{\ell}.
\eqn
In the limit as $dx \rightarrow 0$ the amount of mass and energy stored within the control volume becomes negligible and we are left with the following relationships
\beqn
\rho_s(a-V_s)=\rho_{\ell}(a-V_{\ell}),\\
\rho_sH_s\left(\dot{a}-V_s\right) + \mathbf {n} \cdot\mathbf{q_s} =\rho_{\ell}H_{\ell}\left(\dot{a}-V_{\ell}\right)+ \mathbf{n}\cdot \mathbf{q_{\ell}}.
\eqn
Since there is no motion in the solid, $V_s=0$ and the mass conservation relation gives us a relationship for the fluid velocity.  Substituting this relationship into the energy conservation equation and noting that $L=H_s-H_{\ell}$ we obtain the modified Stefan condition. These conditions are

\beqn
V_{\ell}=\totd{a}{t}\left(\frac{\rho_{\ell}-\rho_s}{\rho_{\ell}}\right),\\
\rho_s L\totd{a}{t}=\left(\mathbf{q_{\ell}}-\mathbf{q_s}\right)\cdot\mathbf{n}.
\eqn

The addition of a fluid velocity on the liquid side adds an advective component to the governing temperature equation. The new model  can be solved by a similarity solution as a simple extension of the last section, but does not yield any new results. Since the temperature profile is homogeneous initially it must remain so for all time. On the other hand,  if the liquid were under-cooled (as we shall see in the third lecture) or if we were melting the solid, the advective component would give a small correction to the interface speed as long as $\rho_s\approx \rho_{\ell}$.

%\bibliography{lecture}


\begin{thebibliography}{5}

\bibitem{1} Taber, S. 1929. Frost heaving. {\it J. Geol.}, 37, 428-461.

\end{thebibliography}

\end{document}