Sea ice growth and analytical models

Sea ice covers about 7% of the Earth’s surface and its presence greatly affects the exchange of heat and momentum between the atmosphere and the ocean. The sea ice growth, and hence its associated thermodynamics, is a hot topic in Earth sciences given the ongoing changes in the Arctic and the consequent climate implications of that. For those who are interested, take a look at this introductory video made by NASA.

Everyone has produced ice cubes at least once in their life by just putting water in an ice tray and letting it cool down in the fridge. The idea behind that is quite intuitive: you decrease the temperature of a certain volume of water until it reaches the freezing point. This same concept can also be applied for sea ice, but the presence of the atmosphere with variable temperatures, winds and waves, snow, solar radiation, ocean and its internal heat content and salinity makes the problem way more difficult than just making ice cubes. How can we present this phenomenon in simple physical terms?

You will learn in this longer problemset the basic concepts of sea ice growth and investigate how its thickness can change according to different environmental conditions.

Mathematical Theory

Let’s start with the equation governing the heat transfer within the sea ice:

\[\frac{\partial}{\partial t}(\rho_i c_i T) = \nabla \cdot (\kappa_i \nabla T) + q \tag{1}\]

where \(t\) is time, \(\rho_{i}\) is the ice density, \(c_i\) is the specific heat of ice, \(\textit{T}\) is the ice temperature, \(\kappa_i\) is the heat conductivity of ice and \(q\) is an internal source term.

The left side of Eq. 1 represents how heat in sea ice changes in time, and the right side of the equation tells us about the sources and sinks of heat that promote such changes. As one can imagine, solving Equation 1 analytically is impossible, but we can simplify this.

The first simplification considered here is that vertical temperature gradients are usually much larger than the horizontal ones. In other words, Equation 1 may be rewritten as:

\[\frac{\partial}{\partial t}(\rho_i c_i T) = \frac{\partial}{\partial z}(\kappa_i \frac{\partial T}{\partial z}) + q \tag{2}\]

Since we are dealing with a partial differential equation, boundary conditions must be defined on the top of the sea ice, where vertical changes in temperature is determined by the heat flux (\(Q_T\)) at the surface

\[\kappa_i \frac{\partial T}{\partial z} = Q_T \tag{3}\]

and at the bottom, where the temperature is fixed to the freezing temperature of the sea water (\(T_f\)):

\[T = T_f \tag{4}\]

Variations in time of sea ice thickness (\(H\)), like freezing and melting, take place in the bottom and thus we have also to define how these changes happen:

\[\rho_i L \frac{dH}{dt} = \kappa_i \frac{\partial T}{\partial z} \vert \ _{bottom} - Q_w \tag{5}\]

where \(L\) is the latent heat of freezing, \(H\) is the ice thickness and \(Q_w\) is the heat flux from the water to the sea ice.

Stefan’s model

This set of equations (2 - 5) describes the evolution of sea ice thickness and they are the core of this task. However, we are still dealing with many variables that still make our search for an analytical model difficult.

Josef Stefan proposed in 1891 several simplifications to this set of equations to overcome this issue:

  1. No thermal inertia (\(c_i = 0\))

  2. No internal heat source (\(q = 0\))

  3. A known temperature at the top (\(T_0 = T_0(t)\))

  4. No heat flux from the water (\(Q_w = 0\))

  5. No snow

Assumptions (3) and (4) simplifies Equation 5 to:

\[\rho_i L \frac{dH}{dt} = \kappa_i \frac{T_f - T_0}{H} \tag{6}\]

where the vertical change in temperature (\(\frac{\partial T}{\partial z}\)) in Equation 5 is replaced by the initial conditions (\(T_0 \rightarrow\) top and \(T_f \rightarrow\) bottom)

We can then integrate Equation 6 with respect to time (dt) from \(t_0\) to \(t\):

\[\int_{t_0}^{t} H \frac{dH}{dt} dt = \kappa_i \frac{T_f - T_0}{\rho_i L} dt \tag{7}\]

This leads to:

\[\frac{H(t)^2 - H_0^2}{2} = \kappa_i \frac{T_f - T_0}{\rho_i L} * t \tag{8}\]

where \(H_0\) is the ice thickness at the initial step \(t_0\). Rearranging the terms, we find that:

\[H(t) = \sqrt{H_0^2 + 2 * \kappa_i \frac{T_f - T_0}{\rho_i L} * t} \tag{9}\]

Where the units of each variable are as follows:

  • \(H\), \(H_{0}\): meters (m)

  • \(T{_f}\), \(T{_0}\): \(^{\circ}C\)

  • \(\kappa_i\): \(\frac{W}{m ^{\circ}C}\)

  • \(\rho_i\): \(\frac{kg}{m^3}\)

  • \(L\): \(\frac{kJ}{kg}\)

  • \(t\): time (days)

Equation 9 is the one that will be investigated here.

Task 1

Implement Equation 9 in the function stefan_ice(), having as inputs \(H_{0}\), \(\kappa_i\), \(T_{f}\), \(T_{0}\), \(L\), \(\rho_i\) and time \(t\).

def stefan_ice(h0, ki, tf, t0, L, pi, t):
    h = []
    for i in range(0, len(t)):
        dh = (h0**2 + ((2*ki*(tf-t0))/(pi*L))*t[i])**0.5
        if i == 0:
            h.append(h0)
        else:
            h.append(dh)
    return h

Task 2

Using the following constants and stefan_ice(), estimate the sea ice growth over a period of 30 days.

\[ \begin{align}\begin{aligned}H_{0} = 0 m\\\kappa_i = 2 \frac{W}{m^{\circ}C}\\\rho_{i} = 917 \frac{kg}{m^3}\\L = 3.3 \times 10^{5} \frac{J}{kg}\\T_{f} = -2^{\circ}C\\T_{0} = -15^{\circ}C\end{aligned}\end{align} \]
import numpy as np
import matplotlib.pyplot as plt

h0 = 0 # initial thickness
ki = 2 # thermal conductivity
pi = 917 # sea ice density
L = 3.3e5 # latent heat of fusion
tf = -2 # freezing point of water
t0 = -15 # Temperature at the top

t = np.arange(0, 86400*30+24, 24)

h = stefan_ice(h0, ki, tf, t0, L, pi, t)

Task 3

Redo the sea ice growth estimations considering now an initial thickness of \(1 m\).

h1 = 1

h_1 = stefan_ice(h1, ki, tf, t0, L, pi, t)

Why, in physical terms, the sea ice grows rate is smaller here than in Task 2?

Task 4

It is not uncommon to spot snow on the top of sea ice in the high north. Consider now the presence of a snow pack with thickness \(h_{s}\) proportional to the sea ice thickness \(h_{s} \simeq \lambda H\). Equation 9 can be rewritten as:

\[H(t) = \sqrt{H_0^2 + 2*\kappa_i \frac{T_f - T_0}{\rho_i L*(1 + \frac{\lambda \kappa_i}{\kappa_s})} * \textit{t}} \tag{10}\]

Where \(\kappa_s\) is the heat conductivity of snow.

Implement Equation 10 in a function called stefan_ice_snow(). Use this function to reproduce an experiment with the same values as in Task 2 but using now the values \(\kappa_s\) and \(\lambda\) as additional input variables (for a total of 9 input variables to your function). Plot the three curves in the same plot.

def stefan_ice_snow(h0, ki, tf, t0, L, pi, lambd, ks, t):
    h = []
    for i in range(0, len(t)):
        dh = (h0**2 + ((2*ki*(tf-t0))/(pi*L*(1+((lambd*ki)/(ks)))))*t[i])**0.5
        if i == 0:
            h.append(h0)
        else:
            h.append(dh)
    return h

Now use stefan_ice_snow() to find snow pack thickness with for the values \(\kappa_s = 0.2 \frac{W}{m^{\circ}C}\) and \(\lambda = 0.3\). Plot the three curves in the same plot.

ks = 0.2
lambd = 0.3

hs = stefan_ice_snow(h1, ki, tf, t0, L, pi, lambd, ks, t)

Why does the presence of a snow cap makes the sea ice grow slower?

Task 5

When sea ice forms, salt is expelled to the underlying ocean through cavities called brine channels. An approximation of this process is given by the following equation:

\[\Delta S = (S - S_i) \frac{h}{D} \frac{\rho_i}{\rho_w} \tag{11}\]

Where \(\Delta S\) is the variation of salinity in the water column, \(S\) is the salinity in the underlying ocean, \(S_{i}\) is the salinity in the formed sea ice, \(h\) is the sea ice thickness (m),\(D\) is the ocean depth (m) and \(\rho_i\) and \(\rho_w\) are sea ice and water densities, respectively.

Use Equation 11 to estimate the increase of salinity in a 100 m deep ocean with a salinity \(S = 35\). Consider the salinity of sea ice \(S_i = 5\), \(\rho_i = 917 \frac{kg}{m^3}\) and \(\rho_w = 1025 \frac{kg}{m^3}\).

Use as \(h\) the values that you obtained at the last time step in the previous three tasks.

H = [h, h_1, hs] # ice thickness found previously
S = 35 # ocean salinity
Si = 5 # sea ice salinity
D = 100 # water depth
pw = 1025 # water density

for i in H:
    ice_t = i[-1] # take the last value

    d_S = (S-Si)*(ice_t/D)*(pi/pw)

    print(d_S)
0.17912007303045568
0.32267216115501274
0.2829387268696522

These values, ranging from 0.17 to 0.32, might look almost insignificant but salt injection in Arctic waters due to sea ice growth in the Barent Sea is the key component in water formation that insulates the sea ice from warmer waters coming from the North Atlantic (AW). Altering this balance might facilitate the contact of the AW with sea ice, promoting basal melting and thus increase Arctic sea ice reduction.