2.7. Thermodynamics¶
The current Icepack version includes two thermodynamics
options, the Bitz and Lipscomb model [6]
(ktherm = 1) that assumes a fixed salinity profile, and a “mushy”
formulation (ktherm = 2) in which salinity evolves
[71]. For each thickness category, Icepack computes
changes in the ice and snow thickness and vertical temperature profile
resulting from radiative, turbulent, and conductive heat fluxes. The ice
has a temperature-dependent specific heat to simulate the effect of
brine pocket melting and freezing, for ktherm = 1 and 2.
Each thickness category \(n\) in each grid cell is treated as a
horizontally uniform column with ice thickness
\(h_{in} = v_{in}/a_{in}\) and snow thickness
\(h_{sn} = v_{sn}/a_{in}\). (Henceforth we omit the category
index \(n\).) Each column is divided into \(N_i\) ice layers
of thickness \(\Delta h_i = h_i/N_i\) and \(N_s\) snow layers of
thickness \(\Delta h_s = h_s/N_s\). Minimum ice and snow thickness
is specified by namelist parameters hi_min and hs_min.
The surface temperature (i.e., the temperature of ice or snow at the interface with the atmosphere) is \(T_{sf}\), which cannot exceed \(0^{\circ}C\). The temperature at the midpoint of the snow layer is \(T_s\), and the midpoint ice layer temperatures are \(T_{ik}\), where \(k\) ranges from 1 to \(N_i\). The temperature at the bottom of the ice is held at \(T_f\), the freezing temperature of the ocean mixed layer. All temperatures are in degrees Celsius unless stated otherwise.
The tfrz_option namelist specifies the freezing temperature formulation.
minus1p8 fixes the freezing temperature at -1.8C. constant fixes
the freeing point at whatever value is specified by the parameter Tocnfrz.
linear_salt sets the freezing temperature based on salinity,
\(Tf = -depressT * sss\). And mushy uses the mushy formulation for setting
the freezing temperature.
Each ice layer has an enthalpy \(q_{ik}\), defined as the negative of the energy required to melt a unit volume of ice and raise its temperature to \(0^{\circ}C\). Because of internal melting and freezing in brine pockets, the ice enthalpy depends on the brine pocket volume and is a function of temperature and salinity. We can also define a snow enthalpy \(q_s\), which depends on temperature alone.
Given surface forcing at the atmosphere–ice and ice–ocean interfaces
along with the ice and snow thicknesses and temperatures/enthalpies at
time \(m\), the thermodynamic model advances these quantities to
time \(m+1\) (ktherm = 2 also advances salinity). The calculation
proceeds in two steps. First we solve a set of equations for the new
temperatures, as discussed in the New temperatures section. Then we
compute the melting, if any, of ice or snow at the top surface, and the
growth or melting of ice at the bottom surface, as described in
the Growth and melting section. We begin by describing the melt ponds and surface
forcing parameterizations, which are closely related to the ice and snow
surface temperatures.
Sometimes instabilities can arise when the temperature is close to the melt point in
the snow and sea ice and there is abundant internal shorwave absorbed. One can choose
to “move” the excess internal shortwave in this case up to the top surface to be reabsorbed.
The namelist parameters for this option are sw_redist, sw_frac, and sw_dtemp.
By default, sw_redist is set to .false.
2.7.1. Snow fraction¶
In several places in the code, the snow fraction over ice (either sea ice or pond lids) varies
as a function of snow depth. That is, thin layers of snow are assumed to be patchy, which
allows the shortwave flux to increase gradually as the layer thins, preventing sudden changes
in the shortwave reaching the sea ice (which can cause the thermodynamics solver to not converge).
For example, the parameter snowpatch is used for the CCSM3 radiation scheme, with a default
value of 0.02:
The parameters hs0 and hs1 are used similarly for delta-Eddington radiation calculations with
meltponds, with hs0 over sea ice and hs1 over pond ice.
In the tests shown in [27], \(h_{s0}=0\) for all cases except with the cesm pond scheme; that pond scheme has now been deprecated. \(h_{s0}\) can be used with the topo pond scheme, although its impacts have not been documented. We enforce \(hs0=0\) for level-ice ponds because the infiltration of snow by pond water accomplishes the gradual radiative forcing transition for which the patchy-snow parameters were originally intended. When level-ice ponds are not used, then a typical value for hs0 is 0.03.
With level-ice ponds, the pond water is allowed to infiltrate snow over the level ice area, invisible to the radiation scheme, until the water becomes deep enough to show through the snow layer. The pond fraction is computed during this process and then used to set the snow fraction such that \(f_{snow}+f_{pond}=1\). The ponds are only on the level ice area, and so there is still snow on the ridges even if the entire level ice area becomes filled with ponds.
See [27] for a discussion of the impacts of varying hs1, whose default value is 0.03.
2.7.2. Melt ponds¶
Three explicit melt pond parameterizations are available in Icepack, and all must use the delta-Eddington radiation scheme, described below. The ccsm3 shortwave parameterization incorporates melt ponds implicitly by adjusting the albedo based on surface conditions.
For each of the three explicit parameterizations, a volume \(\Delta V_{melt}\) of melt water produced on a given category may be added to the melt pond liquid volume:
where
is the fraction of the total melt water available that is added to the
ponds, \(\rho_i\) and \(\rho_s\) are ice and snow densities,
\(\Delta h_i\) and \(\Delta h_s\) are the thicknesses of ice and
snow that melted, and \(F_{rain}\) is the rainfall rate. Namelist
parameters are set for the level-ice (tr_pond_lvl) parameterization;
in the cesm and topo pond schemes the standard values of \(r_{max}\)
and \(r_{min}\) are 0.7 and 0.15, respectively.
Radiatively, the surface of an ice category is divided into fractions of snow, pond and bare ice. In these melt pond schemes, the actual pond area and depth are maintained throughout the simulation according to the physical processes acting on it. However, snow on the sea ice and pond ice may shield the pond and ice below from solar radiation. These processes do not alter the actual pond volume; instead they are used to define an “effective pond fraction” (and likewise, effective pond depth, snow fraction and snow depth) used only for the shortwave radiation calculation.
In addition to the physical processes discussed below, tracer equations and definitions for melt ponds are also described in the Tracers section.
2.7.2.1. Topographic formulation (tr_pond_topo = true)¶
The principle concept of this scheme is that melt water runs downhill under the influence of gravity and collects on sea ice with increasing surface height starting at the lowest height [16][17][18]. Thus, the topography of the ice cover plays a crucial role in determining the melt pond cover. However, Icepack does not explicitly represent the topography of sea ice. Therefore, we split the existing ice thickness distribution function into a surface height and basal depth distribution assuming that each sea ice thickness category is in hydrostatic equilibrium at the beginning of the melt season. We then calculate the position of sea level assuming that the ice in the whole grid cell is rigid and in hydrostatic equilibrium.
Melt Ponds¶
Figure Melt Ponds illustrates (a) Schematic illustration of the relationship between the height of the pond surface \(h_{pnd,tot}\), the volume of water \(V_{Pk}\) required to completely fill up to category \(k\), the volume of water \(V_{P} - V_{Pk}\), and the depth to which this fills category \(k + 1\). Ice and snow areas \(a_i\) and \(a_s\) are also depicted. The volume calculation takes account of the presence of snow, which may be partially or completely saturated. (b) Schematic illustration indicating pond surface height \(h_{pnd,tot}\) and sea level \(h_{sl}\) measured with respect to the thinnest surface height category \(h_{i1}\), the submerged portion of the floe \(h_{sub}\), and hydraulic head \(\Delta H\) . A positive hydraulic head (pond surface above sea level) will flush melt water through the sea ice into the ocean; a negative hydraulic head can drive percolation of sea water onto the ice surface. Here, \(\alpha=0.6\) and \(\beta=0.4\) are the surface height and basal depth distribution fractions. The height of the steps is the height of the ice above the reference level, and the width of the steps is the area of ice of that height. The illustration does not imply a particular assumed topography, rather it is assumed that all thickness categories are present at the sub-grid scale so that water will always flow to the lowest surface height class.
Once a volume of water is produced from ice and snow melting, we calculate the number of ice categories covered by water. At each time step, we construct a list of volumes of water \(\{V_{P1}, V_{P2}, . . . V_{P,k-1}, V_{Pk},\) \(V_{P,k+1}, . . . \}\), where \(V_{Pk}\) is the volume of water required to completely cover the ice and snow in the surface height categories from \(i = 1\) to \(i = k\). The volume \(V_{Pk}\) is defined so that if the volume of water \(V_{P}\) is such that \(V_{Pk} < V_{P} < V_{P,k+1}\) then the snow and ice in categories \(n = 1\) to \(n = k + 1\) are covered in melt water. Figure Melt Ponds (a) depicts the areas covered in melt water and saturated snow on the surface height (rather than thickness) categories \(h_{top,k}\). Note in the code, we assume that \(h_{top,n}/h_{in} = 0.6\) (an arbitrary choice). The fractional area of the \(n\)th category covered in snow is \(a_{sn}\). The volume \(V_{P1}\), which is the region with vertical hatching, is the volume of water required to completely fill the first thickness category, so that any extra melt water must occupy the second thickness category, and it is given by the expression
where \(V_{sw}\) is the fraction of the snow volume that can be occupied by water, and \(h_{s1}\) is the snow depth on ice height class 1. In a similar way, the volume required to fill the first and second surface categories, \(V_{P2}\), is given by
The general expression for volume \(V_{Pk}\) is given by
(Note that we have implicitly assumed that \(h_{si} < h_{top,k+1} - h_{top,k}\) for all \(k\).) No melt water can be stored on the thickest ice thickness category. If the melt water volume exceeds the volume calculated above, the remaining melt water is released to the ocean.
At each time step, the pond height above the level of the thinnest surface height class, that is, the maximum pond depth, is diagnosed from the list of volumes \(V_{Pk}\). In particular, if the total volume of melt water \(V_{P}\) is such that \(V_{Pk} < V_{P} < V_{P,k+1}\) then the pond height \(h_{pnd,tot}\) is
where \(h_{par}\) is the height of the pond above the level of the ice in class \(k\) and partially fills the volume between \(V_{P,k}\) and \(V_{P,k+1}\). From Figure Melt Ponds (a) we see that \(h_{top,k} - h_{top,1}\) is the height of the melt water, which has volume \(V_{Pk}\), which completely fills the surface categories up to category \(k\). The remaining volume, \(V_{P} - V_{Pk}\), partially fills category \(k + 1\) to the height \(h_{par}\) and there are two cases to consider: either the snow cover on category \(k + 1\), with height \(h_{s,k+1}\), is completely covered in melt water (i.e., \(h_{par} > h_{s,k+1}\)), or it is not (i.e., \(h_{par} \le h_{s,k+1}\)). From conservation of volume, we see from Figure Melt Ponds (a) that for an incompletely to completely saturated snow cover on surface ice class \(k + 1\),
and for a saturated snow cover with water on top of the snow on surface ice class \(k + 1\),
As the melting season progresses, not only does melt water accumulate upon the upper surface of the sea ice, but the sea ice beneath the melt water becomes more porous owing to a reduction in solid fraction [13]. The hydraulic head of melt water on sea ice (i.e., its height above sea level) drives flushing of melt water through the porous sea ice and into the underlying ocean. The mushy thermodynamics scheme (ktherm = 2) handles flushing. For ktherm \(\ne 2\) we model the vertical flushing rate using Darcy’s law for flow through a porous medium
where \(w\) is the vertical mass flux per unit perpendicular cross-sectional area (i.e., the vertical component of the Darcy velocity), \(\Pi_v\) is the vertical component of the permeability tensor (assumed to be isotropic in the horizontal), \(\mu\) is the viscosity of water, \(\rho_o\) is the ocean density, \(g\) is gravitational acceleration, \(\Delta H\) is the the hydraulic head, and \(h_i\) is the thickness of the ice through which the pond flushes. As proposed by [20] the vertical permeability of sea ice can be calculated from the liquid fraction \(\phi\):
Since the solid fraction varies throughout the depth of the sea ice, so does the permeability. The rate of vertical drainage is determined by the lowest (least permeable) layer, corresponding to the highest solid fraction. From the equations describing sea ice as a mushy layer [14], the solid fraction is determined by:
where \(S\) is the bulk salinity of the ice, \(S_{br}(T)\) is the concentration of salt in the brine at temperature \(T\) and \(c_i\) is the concentration of salt in the ice crystals (set to zero).
The hydraulic head is given by the difference in height between the upper surface of the melt pond \(h_{pnd,tot}\) and the sea level \(h_{sl}\). The value of the sea level \(h_{sl}\) is calculated from
where \(0.4 \sum^{N}_{n=1} a_{in} h_{i,n}\) is the mean thickness of the basal depth classes, and \(h_{sub}\) is the depth of the submerged portion of the floe. Figure Melt Ponds (b) depicts the relationship between the hydraulic head and the depths and heights that appear in Equation (12). The depth of the submerged portion of the floe is determined from hydrostatic equilibrium to be
where \(\rho_m\) is the density of melt water, \(V_P\) is the total pond volume, \(V_s\) is the total snow volume, and \(V_i\) is the total ice volume.
When the surface energy balance is negative, a layer of ice is formed at the upper surface of the ponds. The rate of growth of the ice lid is given by the Stefan energy budget at the lid-pond interface
where \(L_0\) is the latent heat of fusion of pure ice per unit volume, \(T_i\) and \(T_p\) are the ice surface and pond temperatures, and \(k_i\) and \(k_p\) are the thermal conductivity of the ice lid and pond respectively. The second term on the right hand-side is close to zero since the pond is almost uniformly at the freezing temperature [67]. Approximating the temperature gradient in the ice lid as linear, the Stefan condition yields the classic Stefan solution for ice lid depth
where \(\Delta T\) is the temperature difference between the top and the bottom of the lid. Depending on the surface flux conditions the ice lid can grow, partially melt, or melt completely. Provided that the ice lid is thinner than a critical lid depth (1 cm is suggested) then the pond is regarded as effective, i.e. the pond affects the optical properties of the ice cover. Effective pond area and pond depth for each thickness category are passed to the radiation scheme for calculating albedo. Note that once the ice lid has exceeded the critical thickness, snow may accumulate on the lid causing a substantial increase in albedo. In the current CICE model, melt ponds only affect the thermodynamics of the ice through the albedo. To conserve energy, the ice lid is dismissed once the pond is completely refrozen.
As the sea ice area shrinks due to melting and ridging, the pond volume over the lost area is released to the ocean immediately. In [17], the pond volume was carried as an ice area tracer, but in [18] and here, pond area and thickness are carried as separate tracers, as in the Tracers section.
Unlike the cesm and level-ice melt pond schemes, the liquid pond water
in the topo parameterization is not necessarily virtual; it can be
withheld from being passed to the ocean model until the ponds drain by
setting the namelist variable l_mpond_fresh = .true. The refrozen pond
lids are still virtual. Extra code needed to track and enforce
conservation of water has been added to icepack_itd.F90 (subroutine
zap_small_areas), icepack_mechred.F90 (subroutine ridge_shift),
icepack_therm_itd.F90 (subroutines linear_itd and lateral_melt),
and icepack_therm_vertical.F90 (subroutine thermo_vertical), along
with global diagnostics in icedrv_diagnostics.F90.
2.7.2.2. Level-ice formulation (tr_pond_lvl = true)¶
This meltpond parameterization represents a combination of ideas from the empirical CESM melt pond scheme and the topo approach, and is documented in [27]. The ponds evolve according to physically based process descriptions, assuming a thickness-area ratio for changes in pond volume. A novel aspect of the new scheme is that the ponds are carried as tracers on the level (undeformed) ice area of each thickness category, thus limiting their spatial extent based on the simulated sea ice topography. This limiting is meant to approximate the horizontal drainage of melt water into depressions in ice floes. (The primary difference between the level-ice and topo meltpond parameterizations lies in how sea ice topography is taken into account when determining the areal coverage of ponds.) Infiltration of the snow by melt water postpones the appearance of ponds and the subsequent acceleration of melting through albedo feedback, while snow on top of refrozen pond ice also reduces the ponds’ effect on the radiation budget.
Melt pond processes, described in more detail below, include addition of liquid water from rain, melting snow and melting surface ice, drainage of pond water when its weight pushes the ice surface below sea level or when the ice interior becomes permeable, and refreezing of the pond water. If snow falls after a layer of ice has formed on the ponds, the snow may block sunlight from reaching the ponds below. When melt water forms with snow still on the ice, the water is assumed to infiltrate the snow. If there is enough water to fill the air spaces within the snowpack, then the pond becomes visible above the snow, thus decreasing the albedo and ultimately causing the snow to melt faster. The albedo also decreases as snow depth decreases, and thus a thin layer of snow remaining above a pond-saturated layer of snow will have a lower albedo than if the melt water were not present.
The level-ice formulation assumes a thickness-area ratio for changes
in pond volume, while the CESM scheme assumes this ratio for the total
pond volume. Pond volume changes are distributed as changes to the area
and to the depth of the ponds using an assumed aspect ratio, or shape,
given by the parameter \(\delta_p\) (pndaspect),
\(\delta_p = {\Delta h_p / \Delta a_{p}}\) and
\(\Delta V = \Delta h_p \Delta a_{p} = \delta_p\Delta a_p^2 = \Delta h_{p}^2/\delta_p\).
Here, \(a_{p} = a_{pnd} a_{lvl}\), the mean pond area over the ice.
Given the ice velocity \(\bf u\), conservation equations for level ice fraction \(a_{lvl}a_i\), pond area fraction \(a_{pnd}a_{lvl}a_i\), pond volume \(h_{pnd}a_{pnd}a_{lvl}a_i\) and pond ice volume \(h_{ipnd}a_{pnd}a_{lvl}a_i\) are
(We have dropped the category subscript here, for clarity.) Equations (18) and (19) express conservation of melt pond volume and pond ice volume, but in this form highlight that the quantities tracked in the code are the tracers \(h_{pnd}\) and \(h_{ipnd}\), pond depth and pond ice thickness. Likewise, the level ice fraction \(a_{lvl}\) is a tracer on ice area fraction (Equation (16)), and pond fraction \(a_{pnd}\) is a tracer on level ice (Equation (17)).
Pond ice. The ponds are assumed to be well mixed fresh water, and therefore their temperature is 0\(^\circ\)C. If the air temperature is cold enough, a layer of clear ice may form on top of the ponds. There are currently three options in the code for refreezing the pond ice. Only option A tracks the thickness of the lid ice using the tracer \(h_{ipnd}\) and includes the radiative effect of snow on top of the lid.
A. The frzpnd = ‘hlid’ option uses a Stefan approximation for growth of
fresh ice and is invoked only when \(\Delta V_{melt}=0\).
The basic thermodynamic equation governing ice growth is
assuming a linear temperature profile through the ice thickness \(h_i\). In discrete form, the solution is
where
When \(\Delta V_{melt}>0\), any existing pond ice may also melt. In this case,
where \(F_\circ\) is the net downward surface flux.
In either case, the change in pond volume associated with growth or melt of pond ice is
where \(\rho_0\) is the density of fresh water.
B. The frzpnd = ‘cesm’ option uses the same empirical function as in the
CESM melt pond parameterization.
Radiative effects. Freshwater ice that has formed on top of a melt pond is assumed to be perfectly clear. Snow may accumulate on top of the pond ice, however, shading the pond and ice below. The depth of the snow on the pond ice is initialized as \(h_{ps}^0 = F_{snow}\Delta t\) at the first snowfall after the pond ice forms. From that time until either the pond ice or the pond snow disappears, the pond snow depth tracks the depth of snow on sea ice (\(h_s\)) using a constant difference \(\Delta\). As \(h_s\) melts, \(h_{ps}=h_s-\Delta\) will be reduced to zero eventually, at which time the pond ice is fully uncovered and shortwave radiation passes through.
To prevent a sudden change in the shortwave reaching the sea ice (which can prevent the thermodynamics from converging), thin layers of snow on pond ice are assumed to be patchy, thus allowing the shortwave flux to increase gradually as the layer thins. This is done using the same parameterization for patchy snow as is used elsewhere in Icepack, but with its own parameter \(h_{s1}\):
If any of the pond ice melts, the radiative flux allowed to pass through the ice is reduced by the (roughly) equivalent flux required to melt that ice. This is accomplished (approximately) with \(a_{pnd}^{eff} = (1-f_{frac})a_{pnd}a_{lvl}\), where (see Equation (23))
Snow infiltration by pond water. If there is snow on top of the sea ice, melt water may infiltrate the snow. It is a “virtual process” that affects the model’s thermodynamics through the input parameters of the radiation scheme; it does not melt the snow or affect the snow heat content.
A snow pack is considered saturated when its percentage of liquid water content is greater or equal to 15% (Sturm and others, 2009). We assume that if the volume fraction of retained melt water to total liquid content
then effectively there are no meltponds present, that is, \(a_{pnd}^{eff}=h_{pnd}^{eff}=0\). Otherwise, we assume that the snowpack is saturated with liquid water.
We assume that all of the liquid water accumulates at the base of the snow pack and would eventually melt the surrounding snow. Two configurations are therefore possible, (1) the top of the liquid lies below the snow surface and (2) the liquid water volume overtops the snow, and all of the snow is assumed to have melted into the pond. The volume of void space within the snow that can be filled with liquid melt water is
and we compare \(V_p\) with \(V_{mx}\).
Case 1: For \(V_p < V_{mx}\), we define \(V_p^{eff}\) to be the volume of void space filled by the volume \(V_p\) of melt water: \(\rho_0 V_p = (\rho_0-\rho_s) V_p^{eff},\) or in terms of depths,
The liquid water under the snow layer is not visible and therefore the ponds themselves have no direct impact on the radiation (\(a_{pnd}^{eff}=h_{pnd}^{eff}=0\)), but the effective snow thickness used for the radiation scheme is reduced to
Here, the factor \(a_p=a_{pnd}a_{lvl}\) averages the reduced snow depth over the ponds with the full snow depth over the remainder of the ice; that is, \(h_s^{eff} = h_s(1-a_p) + (h_s -h_p^{eff})a_p.\)
Case 2: Similarly, for \(V_p \ge V_{mx}\), the total mass in the liquid is \(\rho_0 V_p + \rho_s V_s = \rho_0 V_p^{eff},\) or
Thus the effective depth of the pond is the depth of the whole slush layer \(h_p^{eff}\). In this case, \(a_{pnd}^{eff}=a_{pnd}a_{lvl}\).
Drainage. A portion \(1-r\) of the available melt water drains immediately into the ocean. Once the volume changes described above have been applied and the resulting pond area and depth calculated, the pond depth may be further reduced if the top surface of the ice would be below sea level or if the sea ice becomes permeable.
We require that the sea ice surface remain at or above sea level. If the weight of the pond water would push the mean ice–snow interface of a thickness category below sea level, some or all of the pond water is removed to bring the interface back to sea level via Archimedes’ Principle written in terms of the draft \(d\),
There is a separate freeboard calculation in the thermodynamics which considers only the ice and snow and converts flooded snow to sea ice. Because the current melt ponds are “virtual” in the sense that they only have a radiative influence, we do not allow the pond mass to change the sea ice and snow masses at this time, although this issue may need to be reconsidered in the future, especially for the Antarctic.
The mushy thermodynamics scheme (ktherm = 2) handles flushing. For ktherm \(\ne 2\), the permeability of the sea ice is calculated using the internal ice temperatures \(T_i\) (computed from the enthalpies as in the sea ice thermodynamics). The brine salinity and liquid fraction are given by [46] [eq 3.6] \(S_{br} = {1/ (10^{-3} - 0.054/T_i)}\) and \(\phi = S/S_{br}\), where \(S\) is the bulk salinity of the combined ice and brine. The ice is considered permeable if \(\phi \ge 0.05\) with a permeability of \(p=3\times 10^{-8}\min(\phi^3)\) (the minimum being taken over all of the ice layers). A hydraulic pressure head is computed as \(P=g\rho_w\Delta h\) where \(\Delta h\) is the height of the pond and sea ice above sea level. Then the volume of water drained is given by
where \(d_p\) is a scaling factor (dpscale), and \(\mu=1.79\times 10^{-3}\) kg m \(^{-1}\) s \(^{-1}\) is the dynamic viscosity.
Conservation elsewhere. When ice ridges and when new ice forms in open water, the level ice area changes and ponds must be handled appropriately. For example, when sea ice deforms, some of the level ice is transformed into ridged ice. We assume that pond water (and ice) on the portion of level ice that ridges is lost to the ocean. All of the tracer volumes are altered at this point in the code, even though \(h_{pnd}\) and \(h_{ipnd}\) should not change; compensating factors in the tracer volumes cancel out (subroutine ridge_shift in icepack_mechred.F90).
When new ice forms in open water, level ice is added to the existing sea ice, but the new level ice does not yet have ponds on top of it. Therefore the fractional coverage of ponds on level ice decreases (thicknesses are unchanged). This is accomplished in icepack_therm_itd.F90 (subroutine add_new_ice) by maintaining the same mean pond area in a grid cell after the addition of new ice,
and solving for the new pond area tracer \(a_{pnd}^\prime\) given the newly formed ice area \(\Delta a_i = \Delta a_{lvl}\).
2.7.3. Thermodynamic surface forcing balance¶
The net surface energy flux from the atmosphere to the ice (with all fluxes defined as positive downward) is
where \(F_s\) is the sensible heat flux, \(F_l\) is the latent heat flux, \(F_{L\downarrow}\) is the incoming longwave flux, \(F_{L\uparrow}\) is the outgoing longwave flux, \(F_{sw}\) is the incoming shortwave flux, \(\alpha\) is the shortwave albedo, and \(i_0\) is the fraction of absorbed shortwave flux that penetrates into the ice. The albedo may be altered by the presence of melt ponds. Each of the explicit melt pond parameterizations (CESM, topo and level-ice ponds) should be used in conjunction with the Delta-Eddington shortwave scheme, described below.
2.7.3.1. Shortwave radiation: Delta-Eddington¶
Two methods for computing albedo and shortwave fluxes are available, the
“ccsm3” method, described in the next section, and a multiple scattering
radiative transfer scheme that uses a Delta-Eddington approach
(shortwave = dEdd).
“Inherent” optical properties (IOPs) for snow and sea ice, such as extinction coefficient and single scattering albedo, are prescribed based on physical measurements; reflected, absorbed and transmitted shortwave radiation (“apparent” optical properties) are then computed for each snow and ice layer in a self-consistent manner. Absorptive effects of inclusions in the ice/snow matrix such as dust and algae can also be included, along with radiative treatment of melt ponds and other changes in physical properties, for example granularization associated with snow aging.
The Delta-Eddington formulation is described in detail in [8]. Since publication of this technical paper, a number of improvements have been made to the Delta-Eddington scheme, including a surface scattering layer and internal shortwave absorption for snow, generalization for multiple snow layers and more than four layers of ice, and updated IOP values.
In addition, a 5-band option for snow has
been added based on [10] using parameters derived from the
SNICAR snow model (shortwave = dEdd_snicar_ad). The 3-band
Delta-Eddington data is still used for non-snow-covered surfaces. The
5-band option calculates snow radiative transfer properties for 1 visible and
4 near-infrared bands, and the reflection, absorption and transmission of
direct and diffuse shorwave incidents are computed separately, thus removing
the snow grain adjustment used in the 3-band Delta-Eddington scheme. Also,
albedo and absorption of snow-covered sea ice are adjusted for solar zenith
angles greater than 75 degrees. Because the 5-band lookup tables are very
large, they can be slow to compile. The setting ICE_SNICARHC is false
for simulations not using the dEdd_snicar_ad option, and must be set
to true in order to use the hard-coded (HC) lookup tables generated from the
SNICAR model.
The namelist parameters R_ice and R_pnd adjust the albedo of bare or
ponded ice by the product of the namelist value and one standard
deviation. For example, if R_ice = 0.1, the albedo increases by
\(0.1\sigma\). Similarly, setting R_snw = 0.1 decreases the snow
grain radius by \(0.1\sigma\) (thus increasing the albedo). Two
additional tuning parameters are available for this scheme, dT_mlt and
rsnw_mlt. dT_mlt is the temperature change needed for a change in snow
grain radius from non-melting to melting, and rsnw_mlt is the maximum
snow grain radius when melting. An absorption coefficient for algae
(kalg) may also be set. See [8] for details; the
CESM melt pond and Delta-Eddington parameterizations are further
explained and validated in [23].
2.7.3.2. Shortwave radiation: CCSM3¶
In the parameterization used in the previous version of the Community Climate System Model (CCSM3), the albedo depends on the temperature and thickness of ice and snow and on the spectral distribution of the incoming solar radiation. Albedo parameters have been chosen to fit observations from the SHEBA field experiment. For \(T_{sf} < -1^{\circ}C\) and \(h_i >\) ahmax, the bare ice albedo is 0.78 for visible wavelengths (\(<700\) nm) and 0.36 for near IR wavelengths (\(>700\) nm). As \(h_i\) decreases from ahmax to zero, the ice albedo decreases smoothly (using an arctangent function) to the ocean albedo, 0.06. The ice albedo in both spectral bands decreases by 0.075 as \(T_{sf}\) rises from \(-1^{\circ}C\) to . The albedo of cold snow (\(T_{sf} < -1^{\circ}C\)) is 0.98 for visible wavelengths and 0.70 for near IR wavelengths. The visible snow albedo decreases by 0.10 and the near IR albedo by 0.15 as \(T_{sf}\) increases from \(-1^{\circ}C\) to \(0^{\circ}C\). The total albedo is an area-weighted average of the ice and snow albedos, where the fractional snow-covered area is
and \(h_{snowpatch} = 0.02 \ {\mathrm m}\). The envelope of albedo values is shown in Figure Albedo. This albedo formulation incorporates the effects of melt ponds implicitly; the explicit melt pond parameterization is not used in this case.
Albedo¶
Figure Albedo illustrates Albedo as a function of ice thickness and temperature, for the two extrema in snow depth, for the ccsm3 shortwave option. Maximum snow depth is computed based on Archimedes’ Principle for the given ice thickness. These curves represent the envelope of possible albedo values.
The net absorbed shortwave flux is \(F_{swabs} = \sum (1-\alpha_j) F_{sw\downarrow}\), where the summation is over four radiative categories (direct and diffuse visible, direct and diffuse near infrared). The flux penetrating into the ice is \(I_0 = i_0 \, F_{swabs}\), where \(i_0 = 0.70 \, (1-f_{snow})\) for visible radiation and \(i_0 = 0\) for near IR. Radiation penetrating into the ice is attenuated according to Beer’s Law:
where \(I(z)\) is the shortwave flux that reaches depth \(z\) beneath the surface without being absorbed, and \(\kappa_i\) is the bulk extinction coefficient for solar radiation in ice, set to \(1.4 \ {\mathrm m^{-1}}\) for visible wavelengths [12]. A fraction \(\exp(-\kappa_i h_i)\) of the penetrating solar radiation passes through the ice to the ocean (\(F_{sw\Downarrow}\)).
2.7.3.3. Longwave radiation, turbulent fluxes¶
While incoming shortwave and longwave radiation are obtained from the atmosphere, outgoing longwave radiation and the turbulent heat fluxes are derived quantities. Outgoing longwave takes the standard blackbody form, \(F_{L\uparrow}=\epsilon\sigma \left(T_{sf}^{K}\right)^4\), where \(\epsilon=0.985\) is the emissivity of snow or ice, \(\sigma\) is the Stefan-Boltzmann constant and \(T_{sf}^{K}\) is the surface temperature in Kelvin. (The longwave fluxes are partitioned such that \(\epsilon F_{L\downarrow}\) is absorbed at the surface, the remaining \(\left(1-\epsilon\right)F_{L\downarrow}\) being returned to the atmosphere via \(F_{L\uparrow}\).) The sensible heat flux is proportional to the difference between air potential temperature and the surface temperature of the snow or snow-free ice,
\(C_s\) and \(C_l\) (below) are nonlinear turbulent heat transfer coefficients described in the Atmosphere section. Similarly, the latent heat flux is proportional to the difference between \(Q_a\) and the surface saturation specific humidity \(Q_{sf}\):
where \(q_1 = 1.16378 \times 10^7 \, \mathrm{kg/m^3}\), \(q_2 = 5897.8 \, \mathrm{K}\), \(T_{sf}^K\) is the surface temperature in Kelvin, and \(\rho_a\) is the surface air density.
The net downward heat flux from the ice to the ocean is given by [43]:
where \(\rho_w\) is the density of seawater, \(c_w\) is the specific heat of seawater, \(c_h = 0.006\) is a heat transfer coefficient, \(u_*=\sqrt{\left|\vec{\tau}_w\right|/\rho_w}\) is the friction velocity, and \(T_w\) is the sea surface temperature. A minimum value of \(u_*\) is available; we recommend \(u_{*\min} = 5\times 10^{-4}\) m/s, but the optimal value may depend on the ocean forcing used and can be as low as 0.
\(F_{bot}\) is limited by the total amount of heat available from the ocean, \(F_{frzmlt}\). Additional heat, \(F_{side}\), is used to melt the ice laterally following [44] and [64]. \(F_{bot}\) and the fraction of ice melting laterally are scaled so that \(F_{bot} + F_{side} \ge F_{frzmlt}\) in the case that \(F_{frzmlt}<0\) (melting; see Growth and melting).
2.7.4. New temperatures¶
2.7.4.1. Bitz and Lipscomb thermodynamics (ktherm = 1)¶
The “Bitz99” thermodynamic sea ice model is based on [45] and [6], and is described more fully in [37]. The vertical salinity profile is prescribed and is unchanging in time. The snow is assumed to be fresh, and the midpoint salinity \(S_{ik}\) in each ice layer is given by
where \(z \equiv (k-1/2)/N_i\), \(S_{\max} = 3.2\) ppt, and \(a=0.407\) and \(b=0.573\) are determined from a least-squares fit to the salinity profile observed in multiyear sea ice by [60]. This profile varies from \(S=0\) at the top surface (\(z = 0\)) to \(S=S_{\max}\) at the bottom surface (\(z=1\)) and is similar to that used by [45]. Equation (40) is fairly accurate for ice that has drained at the top surface due to summer melting. It is not a good approximation for cold first-year ice, which has a more vertically uniform salinity because it has not yet drained. However, the effects of salinity on heat capacity are small for temperatures well below freezing, so the salinity error does not lead to significant temperature errors.
Temperature updates
Given the temperatures \(T_{sf}^m\), \(T_s^m\), and \(T_{ik}^m\) at time \(m\), we solve a set of finite-difference equations to obtain the new temperatures at time \(m+1\). Each temperature is coupled to the temperatures of the layers immediately above and below by heat conduction terms that are treated implicitly. For example, the rate of change of \(T_{ik}\) depends on the new temperatures in layers \(k-1\), \(k\), and \(k+1\). Thus we have a set of equations of the form
where \({\bf A}\) is a tridiagonal matrix, \({\bf x}\) is a column vector whose components are the unknown new temperatures, and \({\bf b}\) is another column vector. Given \({\bf A}\) and \({\bf b}\), we can compute \({\bf x}\) with a standard tridiagonal solver.
There are four general cases: (1) \(T_{sf} < 0^{\circ}C\), snow
present; (2) \(T_{sf} = 0^{\circ}C\), snow present;
(3) \(T_{sf} < 0^{\circ}C\), snow absent; and
(4) \(T_{sf} = 0^{\circ}C\), snow absent. For case 1 we have
one equation (the top row of the matrix) for the new surface
temperature, \(N_s\) equations for the new snow temperatures, and
\(N_i\) equations for the new ice temperatures. For cases 2 and 4 we
omit the equation for the surface temperature, which is held at \(0^{\circ}C\), and
for cases 3 and 4 we omit the snow temperature equations. Snow is
considered absent if the snow depth is less than a user-specified
minimum value, hs_min. (Very thin snow layers are still transported
conservatively by the transport modules; they are simply ignored by the
thermodynamics.)
The rate of temperature change in the ice interior is given by [45]:
where \(\rho_i = 917 \ \mathrm {kg/m^{3}}\) is the sea ice density
(assumed to be uniform), \(c_i(T,S)\) is the specific heat of sea
ice, \(K_i(T,S)\) is the thermal conductivity of sea ice,
\(I_{pen}\) is the flux of penetrating solar radiation at
depth \(z\), and \(z\) is the vertical coordinate, defined to be
positive downward with \(z = 0\) at the top surface. If shortwave =
‘ccsm3’, the penetrating radiation is given by Beer’s Law:
where \(I_0\) is the penetrating solar flux at the top ice surface
and \(\kappa_i\) is an extinction coefficient. If shortwave =
‘dEdd’, then solar absorption is computed by the Delta-Eddington scheme.
The specific heat of sea ice is given to an excellent approximation by [48]
where \(c_0 = 2106\) J/kg/deg is the specific heat of fresh ice at , \(L_0 = 3.34 \times 10^5\) J/kg is the latent heat of fusion of fresh ice at , and \(\mu = 0.054\) deg/ppt is the (liquidus) ratio between the freezing temperature and salinity of brine.
Following [72] and
[45], the standard thermal conductivity
(conduct = ‘Maykut71’) is given by
where \(K_0 = 2.03\) W/m/deg is the conductivity of fresh ice and \(\beta = 0.13\) W/m/ppt is an empirical constant. Experimental results [69] suggest that Equation (44) may not be a good description of the thermal conductivity of sea ice. In particular, the measured conductivity does not markedly decrease as \(T\) approaches \(0^{\circ}C\), but does decrease near the top surface (regardless of temperature).
An alternative parameterization based on the “bubbly brine” model of
[50] for conductivity is available
(conduct = ‘bubbly’):
where \(\rho_i\) and \(\rho_0=917\) kg/m\(^3\) are densities of sea ice and pure ice. Whereas the parameterization in Equation (44) asymptotes to a constant conductivity of 2.03 W m\(^{-1}\) K \(^{-1}\) with decreasing \(T\), \(K_i\) in Equation (45) continues to increase with colder temperatures.
The equation for temperature changes in snow is analogous to Equation (42), with \(\rho_s = 330\) kg/m\(^3\), \(c_s = c_0\), and \(K_s = 0.30\) W/m/deg replacing the corresponding ice values. If shortwave = ‘ccsm3’, then the penetrating solar radiation is equal to zero for snow-covered ice, since most of the incoming sunlight is absorbed near the top surface. If shortwave = ‘dEdd’, however, then \(I_{pen}\) is nonzero in snow layers.
It is possible that more shortwave penetrates into an ice layer than is
needed to completely melt the layer, or else it causes the computed
temperature to be greater than the melting temperature, which until now
has caused the vertical thermodynamics code to abort. A parameter
frac = 0.9 sets the fraction of the ice layer than can be melted through.
A minimum temperature difference for absorption of radiation is also
set, currently dTemp = 0.02 (K). The limiting occurs in
icepack_therm_vertical.F90, for both the ccsm3 and delta Eddington
radiation schemes. If the available energy would melt through a layer,
then penetrating shortwave is first reduced, possibly to zero, and if
that is insufficient then the local conductivity is also reduced to
bring the layer temperature just to the melting point.
We now convert Equation (42) to finite-difference form. The resulting equations are second-order accurate in space, except possibly at material boundaries, and first-order accurate in time. Before writing the equations in full we give finite-difference expressions for some of the terms.
First consider the terms on the left-hand side of Equation (42). We write the time derivatives as
where \(T\) is the temperature of either ice or snow and \(m\) is a time index. The specific heat of ice layer \(k\) is approximated as
which ensures that energy is conserved during a change in temperature. This can be shown by using Equation (43) to integrate \(c_i \, dT\) from \(T_{ik}^m\) to \(T_{ik}^{m+1}\); the result is \(c_{ik}(T_{ik}^{m+1} - T_{ik}^m)\), where \(c_{ik}\) is given by Equation (46). The specific heat is a nonlinear function of \(T_{ik}^{m+1}\), the unknown new temperature. We can retain a set of linear equations, however, by initially guessing \(T_{ik}^{m+1} = T_{ik}^m\) and then iterating the solution, updating \(T_{ik}^{m+1}\) in Equation (46) with each iteration until the solution converges.
Next consider the first term on the right-hand side of Equation (42). The first term describes heat diffusion and is discretized for a given ice or snow layer \(k\) as
where \(\Delta h\) is the layer thickness and \(K_{k}\) is the effective conductivity at the upper boundary of layer \(k\). This discretization is centered and second-order accurate in space, except at the boundaries. The flux terms on the right-hand side (RHS) are treated implicitly; i.e., they depend on the temperatures at the new time \(m+1\). The resulting scheme is first-order accurate in time and unconditionally stable. The effective conductivity \(K^*\) at the interface of layers \(k-1\) and \(k\) is defined as
which reduces to the appropriate values in the limits \(K_k \gg K_{k-1}\) (or vice versa) and \(h_k \gg h_{k-1}\) (or vice versa). The effective conductivity at the top (bottom) interface of the ice-snow column is given by \(K^*=2K/\Delta h\), where \(K\) and \(\Delta h\) are the thermal conductivity and thickness of the top (bottom) layer. The second term on the RHS of Equation (42) is discretized as
where \(\tau_k\) is the fraction of the penetrating solar radiation \(I_0\) that is transmitted through layer \(k\) without being absorbed.
We now construct a system of equations for the new temperatures. For \(T_{sf} < 0^{\circ}C\) we require
where \(F_{ct}\) is the conductive flux from the top surface to the ice interior, and both fluxes are evaluated at time \(m+1\). Although \(F_0\) is a nonlinear function of \(T_{sf}\), we can make the linear approximation
where \(T_{sf}^*\) is the surface temperature from the most recent iteration, and \(F_0^*\) and \((dF_0/dT_{sf})^*\) are functions of \(T_{sf}^*\). We initialize \(T_{sf}^* = T_{sf}^m\) and update it with each iteration. Thus we can rewrite Equation (48) as
Rearranging terms, we obtain
the first equation in the set of equations (41). The temperature change in ice/snow layer \(k\) is
where \(T_0 = T_{sf}\) in the equation for layer 1. In tridiagonal matrix form, Equation (50) becomes
where \(\eta_k = \Delta t/(\rho_k c_k \Delta h_k)\). In the equation for the bottom ice layer, the temperature at the ice–ocean interface is held fixed at \(T_f\), the freezing temperature of the mixed layer; thus the last term on the LHS is known and is moved to the RHS. If \(T_{sf} = 0^{\circ}C\) , then there is no surface flux equation. In this case the first equation in Equation (41) is similar to Equation (51), but with the first term on the LHS moved to the RHS.
These equations are modified if \(T_{sf}\) and
\(F_{ct}\) are computed within the atmospheric model and
passed to the host sea ice model (calc_Tsfc = false; see Atmosphere). In this case there
is no surface flux equation. The top layer temperature is computed by an
equation similar to Equation (51) but with the first term on the
LHS replaced by \(\eta_1 F_{ct}\) and moved to the RHS. The
main drawback of treating the surface temperature and fluxes explicitly
is that the solution scheme is no longer unconditionally stable.
Instead, the effective conductivity in the top layer must satisfy a
diffusive CFL condition:
For thin layers and typical coupling intervals (\(\sim 1\) hr),
\(K^*\) may need to be limited before being passed to the atmosphere
via the coupler. Otherwise, the fluxes that are returned to the host sea ice model may
result in oscillating, highly inaccurate temperatures. The effect of
limiting is to treat the ice as a poor heat conductor. As a result,
winter growth rates are reduced, and the ice is likely to be too thin
(other things being equal). The values of hs_min and \(\Delta t\)
must therefore be chosen with care. If hs_min is too small, frequent
limiting is required, but if hs_min is too large, snow will be ignored
when its thermodynamic effects are significant. Likewise, infrequent
coupling requires more limiting, whereas frequent coupling is
computationally expensive.
This completes the specification of the matrix equations for the four cases. We compute the new temperatures using a tridiagonal solver. After each iteration we check to see whether the following conditions hold:
\(T_{sf} \leq 0^{\circ}C\).
The change in \(T_{sf}\) since the previous iteration is less than a prescribed limit, \(\Delta T_{\max}\).
\(F_0 \geq F_{ct}\). (If \(F_0 < F_{ct}\), ice would be growing at the top surface, which is not allowed.)
The rate at which energy is added to the ice by the external fluxes equals the rate at which the internal ice energy is changing, to within a prescribed limit \(\Delta F_{\max}\).
We also check the convergence rate of \(T_{sf}\). If \(T_{sf}\) is oscillating and failing to converge, we average temperatures from successive iterations to improve convergence. When all these conditions are satisfied—usually within two to four iterations for \(\Delta T_{\max} \approx 0.01^{\circ}C\) and \(\Delta F_{max} \approx 0.01 \ \mathrm{W/m^2}\)—the calculation is complete.
To compute growth and melt rates (Growth and melting), we derive expressions for the enthalpy \(q\). The enthalpy of snow (or fresh ice) is given by
Sea ice enthalpy is more complex, because of brine pockets whose salinity varies inversely with temperature. Since the salinity is prescribed, there is a one-to-one relationship between temperature and enthalpy. The specific heat of sea ice, given by Equation (43), includes not only the energy needed to warm or cool ice, but also the energy used to freeze or melt ice adjacent to brine pockets. Equation (43) can be integrated to give the energy \(\delta_e\) required to raise the temperature of a unit mass of sea ice of salinity \(S\) from \(T\) to \(T^\prime\):
If we let \(T^\prime = T_{m} \equiv -\mu S\), the temperature at which the ice is completely melted, we have
Multiplying by \(\rho_i\) to change the units from \(\mathrm {J/kg}\) to \(\mathrm {J/m^{3}}\) and adding a term for the energy needed to raise the meltwater temperature to , we obtain the sea ice enthalpy:
Note that Equation (52) is a quadratic equation in \(T\). Given the layer enthalpies we can compute the temperatures using the quadratic formula:
where
The other root is unphysical.
2.7.4.2. Mushy thermodynamics (ktherm = 2)¶
The “mushy” thermodynamics option treats the sea ice as a mushy layer [14] in which the ice is assumed to be composed of microscopic brine inclusions surrounded by a matrix of pure water ice. Both enthalpy and salinity are prognostic variables. The size of the brine inclusions is assumed to be much smaller than the size of the ice layers, allowing a continuum approximation: a bulk sea-ice quantity is taken to be the liquid-fraction-weighted average of that quantity in the ice and in the brine.
Enthalpy and mushy physics
The mush enthalpy, \(q\), is related to the temperature, \(T\), and the brine volume, \(\phi\), by
where \(q_{br}\) is the brine enthalpy, \(q_i\) is the pure ice enthalpy, \(\rho_i\) and \(c_i\) are density and heat capacity of the ice, \(\rho_{w}\) and \(c_{w}\) are density and heat capacity of the brine and \(L_0\) is the latent heat of melting of pure ice. We assume that the specific heats of the ice and brine are fixed at the values of cp_ice and cp_ocn, respectively. The enthalpy is the energy required to raise the temperature of the sea ice to \(0^{\circ}C\), including both sensible and latent heat changes. Since the sea ice contains salt, it usually will be fully melted at a temperature below \(0^{\circ}C\). Equations (52) and (53) are equivalent except for the density used in the term representing the energy required to bring the melt water temperature to \(0^{\circ}C\) (\(\rho_i\) and \(\rho_w\) in equations (52) and (53), respectively).
The liquid fraction, \(\phi\), of sea ice is given by
where the brine salinity, \(S_{br}\), is given by the liquidus relation using the ice temperature.
Within the parameterizations of brine drainage the brine density is a function of brine salinity [46]:
Outside the parameterizations of brine drainage the densities of brine and ice are fixed at the values of \(\rho_w\) and \(\rho_i\), respectively.
The permeability of ice is computed from the liquid fraction as in [20]:
where \(\phi_\Pi\) is 0.05.
The liquidus relation used in the mushy layer module is based on observations of [4]. A piecewise linear relation can be fitted to observations of Z (the ratio of mass of salt (in g) to mass of pure water (in kg) in brine) to the melting temperature: \(Z = aT + b\). Salinity is the mass of salt (in g) per mass of brine (in kg) so is related to Z by
The data is well fitted with two linear regions,
where
\(T_0\) is the temperature at which the two linear regions meet. Fitting to the data, \(T_0=-7.636^\circ\)C, \(a_1=-18.48 \;\mathrm{g} \;\mathrm{kg}^{-1} \;\mathrm{K}^{-1}\), \(a_2=-10.3085\;\mathrm{g} \;\mathrm{kg}^{-1} \;\mathrm{K}^{-1}\), \(b_1=0\) and \(b_2=62.4 \;\mathrm{g}\;\mathrm{kg}^{-1}\).
Two-stage outer iteration
As for the Bitz99 thermodynamics [6] there are two qualitatively different situations that must be considered when solving for the vertical thermodynamics: the surface can be melting and at the melting temperature, or the surface can be colder than the melting temperature and not melting. In the Bitz99 thermodynamics these two situations were treated within the same iterative loop, but here they are dealt with separately. If at the beginning of the time step the ice surface is cold and not melting, we solve the ice temperatures assuming that this is also true at the end of the time step. Once we have solved for the new temperatures we test to see if the answer is consistent with this assumption. If the surface temperature is below the melting temperature then we have found the appropriate consistent solution. If the surface is above the melting temperature at the end of the initial solution attempt, we recalculate the new temperatures assuming the surface temperature is fixed at the melting temperature. Alternatively if the surface is at the melting temperature at the start of a time step, we assume initially that this is also the case at the end of the time step, solve for the new temperatures and then check that the surface conductive heat flux is less than the surface atmospheric heat flux as is required for a melting surface. If this is not the case, the temperatures are recalculated assuming the surface is colder than melting. We have found that solutions of the temperature equations that only treat one of the two qualitatively different solutions at a time are more numerically robust than if both are solved together. The surface state rarely changes qualitatively during the solution so the method is also numerically efficient.
Temperature updates
During the calculation of the new temperatures and salinities, the liquid fraction is held fixed at the value from the previous time step. Updating the liquid fraction during the Picard iteration described below was found to be numerically unstable. Keeping the liquid fraction fixed drastically improves the numerical stability of the method without significantly changing the solution.
Temperatures are calculated in a similar way to Bitz99 with an outer Picard iteration of an inner tridiagonal matrix solve. The conservation equation for the internal ice temperatures is
where \(q\) is the sea ice enthalpy, \(K\) is the bulk thermal conductivity of the ice, \(w\) is the vertical Darcy velocity of the brine, \(q_{br}\) is the brine enthalpy and \(F\) is the internally absorbed shortwave radiation. The first term on the right represents heat conduction and the second term represents the vertical advection of heat by gravity drainage and flushing.
The conductivity of the mush is given by
where \(K_i = 2.3 \mathrm{Wm}^{-1}\mathrm{K}^{-1}\) is the conductivity of pure ice and \(K_{br}=0.5375 \mathrm{Wm}^{-1}\mathrm{K}^{-1}\) is the conductivity of the brine. The thermal conductivity of brine is a function of temperature and salinity, but here we take it as a constant value for the middle of the temperature range experienced by sea ice, \(-10^\circ\)C [63], assuming the brine liquidus salinity at \(-10^\circ\)C.
We discretize the terms that include temperature in the heat conservation equation as
where the superscript signifies whether the quantity is evaluated at the start (\(t_0\)) or the end (\(t\)) of the time step and the subscript indicates the vertical layer. Writing out the temperature dependence of the enthalpy term we have
The mush thermal conductivities are fixed at the start of the timestep. For the lowest ice layer \(T_{k+1}\) is replaced with \(T_{bot}\), the temperature of the ice base. \(\Delta h\) is the layer thickness and \(z^\prime_k\) is the distance between the \(k-1\) and \(k\) layer centers.
Similarly, for the snow layer temperatures we have the following discretized equation:
For the upper-most layer (either ice layer or snow layer if it present) \(T_{k-1}\) is replaced with \(T_{sf}\), the temperature of the surface.
If the surface is colder than the melting temperature then we also have to solve for the surface temperature, \(T_{sf}\). Here we follow the methodology of Bitz99 described above.
These discretized temperature equations form a tridiagional matrix for
the new temperatures and are solved with a standard tridiagonal solver.
A Picard iteration is used to incorporate nonlinearity in the equations.
The surface heat flux is a function of surface temperature and with each
iteration, the surface heat flux is calculated with the new surface
temperature until convergence is achieved. Convergence normally occurs
after a few iterations once the temperature changes during an iteration
fall below \(5\times10^{-4}\;^\circ\mathrm{C}\) and the energy
conservation error falls below 0.9 ferrmax.
Salinity updates
Several physical processes alter the sea ice bulk salinity. New ice forms with the salinity of the sea water from which it formed. Gravity drainage reduces the bulk salinity of newly formed sea ice, while flushing of melt water through the ice also alters the salinity profile.
The salinity equation takes the form
where \(w\) is a vertical Darcy velocity and \(G\) is a source term. The right-hand side depends indirectly on the bulk salinity through the liquid fraction (\(S = \phi S_{br}\)). Since \(\phi\) is fixed for the time step, we solve the salinity equation explicitly after the temperature equation is solved.
A. Gravity drainage. Sea ice initially retains all the salt present in the sea water from which it formed. Cold temperatures near the top surface of forming sea ice result in higher brine salinities there, because the brine is always at its melting temperature. This colder, saltier brine is denser than the underlying sea water and the brine undergoes convective overturning with the ocean. As the dense, cold brine drains out of the ice, it is replaced by fresher seawater, lowering the bulk salinity of the ice. Following [71], gravity drainage is assumed to occur as two simultaneously operating modes: a rapid mode operating principally near the ice base and a slow mode occurring everywhere.
Rapid drainage takes the form of a vertically varying upward Darcy flow. The contribution to the bulk salinity equation for the rapid mode is
where \(S\) is the bulk salinity and \(B_{br}\) is the brine salinity, specified by the liquidus relation with ice temperature. This equation is discretized using an upwind advection scheme,
The upward advective flow also carries heat, contributing a term to the heat conservation Equation (54),
where \(q_{br}\) is the brine enthalpy. This term is discretized as
where the maximum is taken over all the ice layers between layer \(k\) and the ice base. \(\tilde{w}_j\) is given by
where \(Ra_c\) is a critical Rayleigh number and \(Ra(z)\) is the local Rayleigh number at a particular level,
where \(\Delta \rho\) is the difference in density between the brine at \(z\) and the ocean, \(\Pi\) is the minimum permeability between \(z\) and the ocean, \(h\) is the ice thickness, \(\kappa\) is the brine thermal diffusivity and \(\eta\) is the brine dynamic viscosity. Equation (55) reduces the flow rate for Rayleigh numbers below the critical Rayleigh number.
The unmodified flow rate, \(w\), is determined from a hydraulic pressure balance argument for upward flow through the mush and returning downward flow through ice free channels:
where
There are three tunable parameters in the above parameterization,
\(a\), the diameter of the channel, \(\Delta x\), the horizontal
size of the mush draining through each channel, and \(Ra_c\), the
critical Rayleigh number. \(\rho_p\) is the density of brine in the
channel which we take to be the density of brine in the mush at the
level that the brine is draining from. \(l\) is the thickness of
mush from the ice base to the top of the layer in question. We assume
that \(\Delta x\) is proportional to \(l\) so that
\(\Delta x = 2 \beta l\). \(a\) (a_rapid_mode), \(\beta\)
(aspect_rapid_mode) and \(Ra_c\) (Ra_c_rapid_mode) are all
namelist parameters with default values of \(0.5\;\mathrm{mm}\), 1
and 10, respectively. The value \(\beta=1\) gives a square aspect
ratio for the convective flow in the mush.
The slow drainage mode takes the form of a simple relaxation of bulk salinity:
The decay constant, \(\lambda\), is modeled as
where \(S^\ast\) is a tuning parameter for the drainage strength, \(T_{bot}\) is the basal ice temperature, \(T_{sf}\) is the upper surface temperature and \(h\) is the ice thickness. The bulk salinity relaxes to a value, \(S_c(z)\), given by
where \(S_{br}(z)\) is the brine salinity at depth \(z\) and
\(\phi_c\) is a critical liquid fraction. Both \(S^\ast\) and
\(\phi_c\) are namelist parameters,
dSdt_slow_mode \(=1.5\times10^{-7}\;\mathrm{m}\;\mathrm{s}^{-1}\;\mathrm{K}^{-1}\)
and phi_c_slow_mode \(=0.05\).
B. Downwards flushing. Melt pond water drains through sea ice and flushes out brine, reducing the bulk salinity of the sea ice. This is modeled with the mushy physics option as a vertical Darcy flow through the ice that affects both the enthalpy and bulk salinity of the sea ice:
These equations are discretized with an upwind advection scheme. The flushing Darcy flow, \(w_f\), is given by
where \(\overline{\Pi}\) is the harmonic mean of the ice layer permeabilities and \(\Delta h\) is the hydraulic head driving melt water through the sea ice. It is the difference in height between the top of the melt pond and sea level.
Basal boundary condition
In traditional Stefan problems the ice growth rate is calculated by determining the difference in heat flux on either side of the ice/ocean interface and equating this energy difference to the latent heat of new ice formed. Thus,
where \((1-\phi_i)\) is the solid fraction of new ice formed and
the right hand is the difference in heat flux at the ice–ocean interface
between the ice side and the ocean side of the interface. However, with
mushy layers there is usually no discontinuity in solid fraction across
the interface, so \(\phi_i=1\) and Equation (56)
cannot be used explicitly. To circumvent this problem we set the
interface solid fraction to be 0.15, a value that reproduces
observations. \(\phi_i\) is a namelist parameter (phi_i_mushy =
0.85). The basal ice temperature is set to the liquidus temperature
\(T_f\) of the ocean surface salinity.
Tracer consistency
In order to ensure conservation of energy and salt
content, the advection routines will occasionally limit changes to
either enthalpy or bulk salinity. The mushy thermodynamics routine
determines temperature from both enthalpy and bulk salinity. Since the
limiting changes performed in the advection routine are not applied
consistently (from a mushy physics point of view) to both enthalpy and
bulk salinity, the resulting temperature may be changed to be greater
than the limit allowed in the thermodynamics routines. If this situation
is detected, the code corrects the enthalpy so the temperature is below
the limiting value. The limiting value, Tliquidus_max can be specified
in namelist. Conservation of energy is ensured by placing the
excess energy in the ocean, and the code writes a warning (see Error Messages and Aborts)
that this has
occurred to the diagnostics file. This situation only occurs with the
mushy thermodynamics, and it should only occur very infrequently and
have a minimal effect on results. The addition of the heat to the ocean
may reduce ice formation by a small amount afterwards.
2.7.5. Growth and melting¶
Melting at the top surface is given by
where \(q\) is the enthalpy of the surface ice or snow layer [1] (recall that \(q < 0\)) and \(\delta h\) is the change in thickness. If the layer melts completely, the remaining flux is used to melt the layers beneath. Any energy left over when the ice and snow are gone is added to the ocean mixed layer. Ice cannot grow at the top surface due to conductive fluxes; however, snow–ice can form. New snowfall is added at the end of the thermodynamic time step.
Growth and melting at the bottom ice surface are governed by
where \(F_{bot}\) is given by Equation (39) and \(F_{cb}\) is the conductive heat flux at the bottom surface:
If ice is melting at the bottom surface, \(q\)
in Equation (58) is the enthalpy of the bottom ice layer. If
ice is growing, \(q\) is the enthalpy of new ice with temperature
\(T_f\) and salinity \(S_{max}\) (ktherm = 1) or ocean surface
salinity (ktherm = 2). This ice is added to the bottom layer.
In general, frazil ice formed in the ocean is added to the thinnest ice category. The new ice is grown in the open water area of the grid cell to a specified minimum thickness; if the open water area is nearly zero or if there is more new ice than will fit into the thinnest ice category, then the new ice is spread over the entire cell.
If tr_fsd=true, a floe size must be assigned to the new frazil ice.
If spectral ocean surface wave forcing is provided (and set using the
namelist option wave_spec_type), this will be used
to calculate a tensile stress on new floes that determines their maximum
possible size [62][53]. If no ocean surface wave forcing
is provided, all floes are assumed to grow as pancakes, at the smallest
possible floe size.
If tr_fsd=true, lateral growth at the edges of exisiting floes may
also occur, calculated using the prognostic floe size distribution as
described in [24] and [51]. The lateral growth
that occurs is a portion of the total new ice growth, depending on the
area of open water close to floe edges. Lateral growth
modifies the ITD and the FSD.
If tr_fsd=true, floes may weld together thermodynamically during
freezing conditions according to the probability that they overlap,
assuming they are replaced randomly on the domain. Evolution of the
FSD is described using a coagulation equation. The total number of floes
that weld with another, per square meter, per unit time, in the
case of a fully covered ice surface was estimated from observations in
[52]. In its original model implementation, with 12 floe size categories,
the tendency term for floe welding was divided by a
constant equal to the area of the largest floe, (approx 2 km^2), with this choice made
as the product of sensitivity studies to balance the climatological tendencies of
wave fracture and welding. So that results do not vary as the number or range of
floe size categories varies, we fix this scaling coefficient, c_weld.
If the latent heat flux is negative (i.e., latent heat is transferred from the ice to the atmosphere), snow or snow-free ice sublimates at the top surface. If the latent heat flux is positive, vapor from the atmosphere is deposited at the surface as snow or ice. The thickness change of the surface layer is given by
where \(\rho\) is the density of the surface material (snow or ice), and \(L_v = 2.501 \times 10^6 \ \mathrm{J/kg}\) is the latent heat of vaporization of liquid water at \(0^{\circ}C\). Note that \(\rho L_v\) is nearly an order of magnitude larger than typical values of \(q\). For positive latent heat fluxes, the deposited snow or ice is assumed to have the same enthalpy as the existing surface layer.
After growth and melting, the various ice layers no longer have equal thicknesses. We therefore adjust the layer interfaces, conserving energy, so as to restore layers of equal thickness \(\Delta h_i = h_i / N_i\). This is done by computing the overlap \(\eta_{km}\) of each new layer \(k\) with each old layer \(m\):
where \(z_m\) and \(z_k\) are the vertical coordinates of the old and new layers, respectively. The enthalpies of the new layers are
If tr_fsd=false, lateral melting is accomplished by multiplying the state variables by
\(1-r_{side}\), where \(r_{side}\) is the fraction of ice melted
laterally [44][64], and adjusting the ice
energy and fluxes as appropriate. We assume a floe diameter of 300 m.
If tr_fsd=true, lateral melting is accomplished using the [44]
lateral heat flux, but applied to the ice using the prognostic floe size distribution
as described in [24] and [51]. Lateral melt modifies
the ITD and the FSD.
2.7.6. Snow-ice formation¶
At the end of the time step we check whether the snow is deep enough to lie partially below the surface of the ocean (freeboard). From Archimedes’ principle, the base of the snow is at sea level when
Thus the snow base lies below sea level when
In this case, for ktherm = 1 (Bitz99) we raise the snow base to sea level
by converting some snow to ice:
In rare cases this process can increase the ice thickness substantially. For this reason snow–ice conversions are postponed until after the remapping in thickness space (Transport in thickness space), which assumes that ice growth during a single time step is fairly small.
For ktherm = 2 (mushy), we model the snow–ice formation process as
follows: If the ice surface is below sea level then we replace some snow
with the same thickness of sea ice. The thickness change chosen is that
which brings the ice surface to sea level. The new ice has a porosity of
the snow, which is calculated as
where \(\rho_s\) is the density of snow and \(\rho_i\) is the density of fresh ice. The salinity of the brine occupying the above porosity within the new ice is taken as the sea surface salinity. Once the new ice is formed, the vertical ice and snow layers are regridded into equal thicknesses while conserving energy and salt.