2.1. Atmosphere and ocean boundary forcing

External forcing data that are relevant to Icepack

Variable	Description	External Interactions
\(z_o\)	Atmosphere level height	From atmosphere model to sea ice model
\(z_{o,s}\)	Atmosphere level height (scalar quantities) (optional)	From atmosphere model to sea ice model
\(\vec{U}_a\)	Wind velocity	From atmosphere model to sea ice model
\(Q_a\)	Specific humidity	From atmosphere model to sea ice model
\(\rho_a\)	Air density	From atmosphere model to sea ice model
\(\Theta_a\)	Air potential temperature	From atmosphere model to sea ice model
\(T_a\)	Air temperature	From atmosphere model to sea ice model
\(F_{sw\downarrow}\)	Incoming shortwave radiation (4 bands)	From atmosphere model to sea ice model
\(F_{L\downarrow}\)	Incoming longwave radiation	From atmosphere model to sea ice model
\(F_{rain}\)	Rainfall rate	From atmosphere model to sea ice model
\(F_{snow}\)	Snowfall rate	From atmosphere model to sea ice model
\(F_{frzmlt}\)	Freezing/melting potential	From ocean model to sea ice model
\(T_w\)	Sea surface temperature	From ocean model to sea ice model
\(S\)	Sea surface salinity	From ocean model to sea ice model
\(\nabla H_o\)	Sea surface slope	From ocean model via flux coupler to sea ice model
\(h_1\)	Thickness of first ocean level (optional)	From ocean model to sea ice model
\(\vec{U}_w\)	Surface ocean currents	From ocean model to sea ice model (available in Icepack driver, not used directly in column physics)
\(\vec{\tau}_a\)	Wind stress	From sea ice model to atmosphere model
\(F_s\)	Sensible heat flux	From sea ice model to atmosphere model
\(F_l\)	Latent heat flux	From sea ice model to atmosphere model
\(F_{L\uparrow}\)	Outgoing longwave radiation	From sea ice model to atmosphere model
\(F_{evap}\)	Evaporated water	From sea ice model to atmosphere model
\(\alpha\)	Surface albedo (4 bands)	From sea ice model to atmosphere model
\(T_{sfc}\)	Surface temperature	From sea ice model to atmosphere model
\(F_{sw\Downarrow}\)	Penetrating shortwave radiation	From sea ice model to ocean model
\(F_{water}\)	Fresh water flux	From sea ice model to ocean model
\(F_{hocn}\)	Net heat flux to ocean	From sea ice model to ocean model
\(F_{salt}\)	Salt flux	From sea ice model to ocean model
\(\vec{\tau}_w\)	Ice-ocean stress	From sea ice model to ocean model
\(F_{bio}\)	Biogeochemical fluxes	From sea ice model to ocean model
\(a_{i}\)	Ice fraction	From sea ice model to both ocean and atmosphere models
\(T^{ref}_{a}\)	2m reference temperature (diagnostic)	From sea ice model to both ocean and atmosphere models
\(Q^{ref}_{a}\)	2m reference humidity (diagnostic)	From sea ice model to both ocean and atmosphere models
\(F_{swabs}\)	Absorbed shortwave (diagnostic)	From sea ice model to both ocean and atmosphere models
\(E(f)\)	Wave spectrum as a function of frequency	From ocean surface wave model to sea ice model

The ice fraction \(a_i\) (aice) is the total fractional ice coverage of a grid cell. That is, in each cell,

\[\begin{split}\begin{array}{cl} a_{i}=0 & \mbox{if there is no ice} \\ a_{i}=1 & \mbox{if there is no open water} \\ 0<a_{i}<1 & \mbox{if there is both ice and open water,} \end{array}\end{split}\]

where \(a_{i}\) is the sum of fractional ice areas for each category of ice. The ice fraction is used by the flux coupler to merge fluxes from the sea ice model with fluxes from the other earth system components. For example, the penetrating shortwave radiation flux, weighted by \(a_i\), is combined with the net shortwave radiation flux through ice-free leads, weighted by (\(1-a_i\)), to obtain the net shortwave flux into the ocean over the entire grid cell. The CESM flux coupler requires the fluxes to be divided by the total ice area so that the ice and land models are treated identically (land also may occupy less than 100% of an atmospheric grid cell). These fluxes are “per unit ice area” rather than “per unit grid cell area.”

In some coupled climate models (for example, recent versions of the U.K. Hadley Centre model) the surface air temperature and fluxes are computed within the atmosphere model and are passed to CICE for use in the column physics. In this case the logical parameter calc_Tsfc in ice_therm_vertical is set to false. The fields fsurfn (the net surface heat flux from the atmosphere), flatn (the surface latent heat flux), and fcondtopn (the conductive flux at the top surface) for each ice thickness category are copied or derived from the input coupler fluxes and are passed to the thermodynamic driver subroutine, thermo_vertical. At the end of the time step, the surface temperature and effective conductivity (i.e., thermal conductivity divided by thickness) of the top ice/snow layer in each category are returned to the atmosphere model via the coupler. Since the ice surface temperature is treated explicitly, the effective conductivity may need to be limited to ensure stability. As a result, accuracy may be significantly reduced, especially for thin ice or snow layers. A more stable and accurate procedure would be to compute the temperature profiles for both the atmosphere and ice, together with the surface fluxes, in a single implicit calculation. This was judged impractical, however, given that the atmosphere and sea ice models generally exist on different grids and/or processor sets.

2.1.1. Atmosphere

The wind velocity, specific humidity, air density and potential temperature at the given level height \(z_\circ\) (optionally \(z_{\circ,s}\), see below) are used to compute transfer coefficients used in formulas for the surface wind stress and turbulent heat fluxes \(\vec\tau_a\), \(F_s\), and \(F_l\), as described below. The sensible and latent heat fluxes, \(F_s\) and \(F_l\), along with shortwave and longwave radiation, \(F_{sw\downarrow}\), \(F_{L\downarrow}\) and \(F_{L\uparrow}\), are included in the flux balance that determines the ice or snow surface temperature when calc_Tsfc is true. As described in the Thermodynamics section, these fluxes depend nonlinearly on the ice surface temperature \(T_{sfc}\). The balance equation is iterated until convergence, and the resulting fluxes and \(T_{sfc}\) are then passed to the flux coupler.

The snowfall precipitation rate (provided as liquid water equivalent and converted by the ice model to snow depth) also contributes to the heat and water mass budgets of the ice layer. Melt ponds generally form on the ice surface in the Arctic and refreeze later in the fall, reducing the total amount of fresh water that reaches the ocean and altering the heat budget of the ice; this version includes two new melt pond parameterizations. Rain and all melted snow end up in the ocean.

Wind stress and transfer coefficients for the turbulent heat fluxes are computed in subroutine atmo_boundary_layer following [32], with additions and changes as detailed in Appendix A of [55] for high frequency coupling (namelist variable highfreq). The resulting equations are provided here for the default boundary layer scheme, which is based on Monin-Obukhov theory (atmbndy = ‘stability’). Alternatively, atmbndy = ‘constant’ provides constant coefficients for wind stress, sensible heat and latent heat calculations (computed in subroutine atmo_boundary_const); atmbndy = ‘mixed’ uses the stability based calculation for wind stress and constant coefficients for sensible and latent heat fluxes.

The wind stress and turbulent heat flux calculation accounts for both stable and unstable atmosphere–ice boundary layers. We first define the “stability”

(1)\[\Upsilon = {\kappa g z_\circ\over u^{*2}} \left({\Theta^*\over\Theta_a\left(1+0.606Q_a\right)} + {Q^*\over 1/0.606 + Q_a}\right),\]

where \(\kappa\) is the von Karman constant, \(g\) is gravitational acceleration, and \(u^*\), \(\Theta^*\) and \(Q^*\) are turbulent scales for velocity difference, temperature, and humidity, respectively, computed as (given the ice velocity \(\vec{U}_i\)):

(2)\[\begin{split}\begin{aligned} u^*&=&c_u\;\textrm{max}\left(U_{\Delta\textrm{min}}, \left|\vec{U}_a - \vec{U}_i \right|\right), \\ \Theta^*&=& c_\theta\left(\Theta_a-T_{sfc}\right), \\ Q^*&=&c_q\left(Q_a-Q_{sfc}\right). \end{aligned}\end{split}\]

Note that atmo_boundary_layer also accepts an optional argument, zlvs, to support staggered atmospheric levels, i.e. receiving scalar quantities from the atmospheric model (humidity and temperature) at a different vertical level than the winds. In that case a separate stability \(\Upsilon_s\) is computed using the same formula as above but substituting \(z_o\) by \(z_{o,s}\).

Within the \(u^*\) expression, \(U_{\Delta\textrm{min}}\) is the minimum allowable value of \(|\vec{U}_{a} - \vec{U}_{i}|\) , which is set to of 0.5 m/s for high frequency coupling (highfreq =.true.). When high frequency coupling is turned off (highfreq =.false.), it is assumed in equation (2) that:

(3)\[\vec{U}_{a} - \vec{U}_{i} \approx \vec{U}_{a}\]

and a higher threshold is taken for \(U_{\Delta\textrm{min}}\) of 1m/s. Equation (3) is a poor assumption when resolving inertial oscillations in ice-ocean configurations where the ice velocity vector may make a complete rotation over a period of \(\ge\) 11.96 hours, as discussed in [55]. However, (3) is acceptable for low frequency ice-ocean coupling on the order of a day or more, when transient ice-ocean Ekman transport is effectively filtered from the model solution. For the \(\Theta^*\) and \(Q^*\) terms in (2), \(T_{sfc}\) and \(Q_{sfc}\) are the surface temperature and specific humidity, respectively. The latter is calculated by assuming a saturated surface, as described in the Thermodynamic surface forcing balance section.

Neglecting form drag, the exchange coefficients \(c_u\), \(c_\theta\) and \(c_q\) are initialized as

(4)\[\kappa\over \ln(z_{ref}/z_{ice})\]

and updated during a short iteration, as they depend upon the turbulent scales. The number of iterations is set by the namelist variable natmiter, nominally set to five but sometimes increased by users employing the highfreq option. A convergence tolerance atmiter_conv on ustar can be set to exit the natmiter loop early if desired. Here, \(z_{ref}\) is a reference height of 10m and \(z_{ice}\) is the roughness length scale for the given sea ice category. \(\Upsilon\) is constrained to have magnitude less than 10. Further, defining \(\chi = \left(1-16\Upsilon\right)^{0.25}\) and \(\chi \geq 1\), the “integrated flux profiles” for momentum and stability in the unstable (\(\Upsilon <0\)) case are given by

(5)\[\begin{split}\begin{aligned} \psi_m = &\mbox{}&2\ln\left[0.5(1+\chi)\right] + \ln\left[0.5(1+\chi^2)\right] -2\tan^{-1}\chi + {\pi\over 2}, \\ \psi_s = &\mbox{}&2\ln\left[0.5(1+\chi^2)\right].\end{aligned}\end{split}\]

In a departure from the parameterization used in [32], we use profiles for the stable case following [31],

(6)\[\psi_m = \psi_s = -\left[0.7\Upsilon + 0.75\left(\Upsilon-14.3\right) \exp\left(-0.35\Upsilon\right) + 10.7\right].\]

The coefficients are then updated as

(7)\[\begin{split}\begin{aligned} c_u^\prime&=&{c_u\over 1+c_u\left(\lambda-\psi_m\right)/\kappa} \\ c_\theta^\prime&=& {c_\theta\over 1+c_\theta\left(\lambda_s-\psi_s\right)/\kappa}\\ c_q^\prime&=&c_\theta^\prime\end{aligned}\end{split}\]

where \(\lambda = \ln\left(z_\circ/z_{ref}\right)\) and \(\lambda_s = \ln\left(z_{\circ,s}/z_{ref}\right)\) if staggered atmospheric levels are used, else \(\lambda_s=\lambda\). The first iteration ends with new turbulent scales from equations (2). After natmiter iterations the latent and sensible heat flux coefficients are computed, along with the wind stress:

(8)\[\begin{split}\begin{aligned} C_l&=&\rho_a \left(L_{vap}+L_{ice}\right) u^* c_q \\ C_s&=&\rho_a c_p u^* c_\theta^* + 1 \\ \vec{\tau}_a&=&{\rho_a (u^{*})^2 \left( \vec{U}_{a} - \vec{U}_{i} \right) \over \left| \vec{U}_{a} - \vec{U}_{i} \right|} \end{aligned}\end{split}\]

where \(L_{vap}\) and \(L_{ice}\) are latent heats of vaporization and fusion, \(\rho_a\) is the density of air and \(c_p\) is its specific heat. Again following [31], we have added a constant to the sensible heat flux coefficient in order to allow some heat to pass between the atmosphere and the ice surface in stable, calm conditions. For the atmospheric stress term in (8), we make the assumption in (3) when highfreq =.false..

The atmospheric reference temperature \(T_a^{ref}\) is computed from \(T_a\) and \(T_{sfc}\) using the coefficients \(c_u\), \(c_\theta\) and \(c_q\). Although the sea ice model does not use this quantity, it is convenient for the ice model to perform this calculation. The atmospheric reference temperature is returned to the flux coupler as a climate diagnostic. The same is true for the reference humidity, \(Q_a^{ref}\).

Additional details about the latent and sensible heat fluxes and other quantities referred to here can be found in the Thermodynamic surface forcing balance section.

2.1.2. Ocean

New sea ice forms when the ocean temperature drops below its freezing temperature. In the Bitz and Lipscomb thermodynamics, [6] \(T_f=-\mu S\), where \(S\) is the seawater salinity and \(\mu=0.054^\circ\)/ppt is the ratio of the freezing temperature of brine to its salinity (linear liquidus approximation). For the mushy thermodynamics, \(T_f\) is given by a piecewise linear liquidus relation. The ocean model calculates the new ice formation; if the freezing/melting potential \(F_{frzmlt}\) is positive, its value represents a certain amount of frazil ice that has formed in one or more layers of the ocean and floated to the surface. (The ocean model assumes that the amount of new ice implied by the freezing potential actually forms.)

If \(F_{frzmlt}\) is negative, it is used to heat already existing ice from below. In particular, the sea surface temperature and salinity are used to compute an oceanic heat flux \(F_w\) (\(\left|F_w\right| \leq \left|F_{frzmlt}\right|\)) which is applied at the bottom of the ice. The portion of the melting potential actually used to melt ice is returned to the coupler in \(F_{hocn}\). The ocean model adjusts its own heat budget with this quantity, assuming that the rest of the flux remained in the ocean.

In addition to runoff from rain and melted snow, the fresh water flux \(F_{water}\) includes ice melt water from the top surface and water frozen (a negative flux) or melted at the bottom surface of the ice. This flux is computed as the net change of fresh water in the ice and snow volume over the coupling time step, excluding frazil ice formation and newly accumulated snow.

Setting the namelist option update_ocn_f to true causes frazil ice to be included in the fresh water and salt fluxes. Some ocean models compute the frazil ice fluxes, which then might need to be corrected for consistency with mushy physics. This behavior is controlled using a combination of update_ocn_f, cpl_frazil and ktherm. In particular, cpl_frazil = 'external' assumes that the frazil ice fluxes are handled entirely outside of Icepack. When ktherm=2, cpl_frazil = 'fresh_ice_correction' sends coupling fluxes representing the difference between the mushy frazil fluxes and fluxes computed assuming the frazil is purely fresh ice. Otherwise the internally computed frazil fluxes are sent to the coupler.

There is a flux of salt into the ocean under melting conditions, and a (negative) flux when sea water is freezing. However, melting sea ice ultimately freshens the top ocean layer, since the ocean is much more saline than the ice. The ice model passes the net flux of salt \(F_{salt}\) to the flux coupler, based on the net change in salt for ice in all categories. In the present configuration, ice_ref_salinity is used for computing the salt flux, although the ice salinity used in the thermodynamic calculation has differing values in the ice layers.

A fraction of the incoming shortwave \(F_{sw\Downarrow}\) penetrates the snow and ice layers and passes into the ocean, as described in the Thermodynamic surface forcing balance section.

A thermodynamic slab ocean mixed-layer parameterization is available in icepack_ocean.F90 and can be run in the full CICE configuration. The turbulent fluxes are computed above the water surface using the same parameterizations as for sea ice, but with parameters appropriate for the ocean. The surface flux balance takes into account the turbulent fluxes, oceanic heat fluxes from below the mixed layer, and shortwave and longwave radiation, including that passing through the sea ice into the ocean. If the resulting sea surface temperature falls below the salinity-dependent freezing point, then new ice (frazil) forms. Otherwise, heat is made available for melting the ice.

When the namelist option calc_dragio is set to true, the ice-ocean drag coefficient, \(c_w\) (dragio), is computed from the thickness of the first ocean level, \(h_1\) (thickness_ocn_layer1), and an under-ice roughness length, \(z_{io}\) (iceruf_ocn). The computation follows [58] :

(9)\[c_w = c_w^* \lambda^2\]

where

(10)\[\begin{split}\begin{aligned} c_w^* &= \frac{\kappa^2} {\ln^2\left( h_1 / z_{io} \right)}, \\ \lambda &= \frac{h_1 - z_{io}} {h_1 \left[ \sqrt{c_w^*} \kappa^{-1} \left( \ln(2) - 1 + z_{io} / h_1 \right) + 1 \right] } \end{aligned}\end{split}\]

2.1.3. Variable exchange coefficients

In the default configuration, atmospheric and oceanic neutral drag coefficients (\(c_u\) and \(c_w\)) are assumed constant in time and space. These constants are chosen to reflect friction associated with an effective sea ice surface roughness at the ice–atmosphere and ice–ocean interfaces. Sea ice (in both Arctic and Antarctic) contains pressure ridges as well as floe and melt pond edges that act as discrete obstructions to the flow of air or water past the ice, and are a source of form drag. Following [70] and based on recent theoretical developments [41][40], the neutral drag coefficients can now be estimated from properties of the ice cover such as ice concentration, vertical extent and area of the ridges, freeboard and floe draft, and size of floes and melt ponds. The new parameterization allows the drag coefficients to be coupled to the sea ice state and therefore to evolve spatially and temporally. This parameterization is contained in the subroutine neutral_drag_coeffs and is accessed by setting formdrag = true in the namelist. (Note: see also Known bugs and other issues.)

Following [70], consider the general case of fluid flow obstructed by N randomly oriented obstacles of height \(H\) and transverse length \(L_y\), distributed on a domain surface area \(S_T\). Under the assumption of a logarithmic fluid velocity profile, the general formulation of the form drag coefficient can be expressed as

(11)\[C_d=\frac{N c S_c^2 \gamma L_y H}{2 S_T}\left[\frac{\ln(H/z_0)}{\ln(z_{ref}/z_0)}\right]^2,\]

where \(z_0\) is a roughness length parameter at the top or bottom surface of the ice, \(\gamma\) is a geometric factor, \(c\) is the resistance coefficient of a single obstacle, and \(S_c\) is a sheltering function that takes into account the shielding effect of the obstacle,

(12)\[S_{c}=\left(1-\exp(-s_l D/H)\right)^{1/2},\]

with \(D\) the distance between two obstacles and \(s_l\) an attenuation parameter.

As in the original drag formulation in CICE (Atmosphere and Ocean sections), \(c_u\) and \(c_w\) along with the transfer coefficients for sensible heat, \(c_{\theta}\), and latent heat, \(c_{q}\), are initialized to a situation corresponding to neutral atmosphere–ice and ocean–ice boundary layers. The corresponding neutral exchange coefficients are then replaced by coefficients that explicitly account for form drag, expressed in terms of various contributions as

(13)\[\tt{Cdn\_atm} = \tt{Cdn\_atm\_rdg} + \tt{Cdn\_atm\_floe} + \tt{Cdn\_atm\_skin} + \tt{Cdn\_atm\_pond} ,\]

(14)\[\tt{Cdn\_ocn} = \tt{Cdn\_ocn\_rdg} + \tt{Cdn\_ocn\_floe} + \tt{Cdn\_ocn\_skin}.\]

The contributions to form drag from ridges (and keels underneath the ice), floe edges and melt pond edges can be expressed using the general formulation of equation (11) (see [70] for details). Individual terms in equation (14) are fully described in [70]. Following [3] the skin drag coefficient is parametrized as

(15)\[{ \tt{Cdn\_(atm/ocn)\_skin}}=a_{i} \left(1-m_{(s/k)} \frac{H_{(s/k)}}{D_{(s/k)}}\right)c_{s(s/k)}, \mbox{ if $\displaystyle\frac{H_{(s/k)}}{D_{(s/k)}}\ge\frac{1}{m_{(s/k)}}$,}\]

where \(m_s\) (\(m_k\)) is a sheltering parameter that depends on the average sail (keel) height, \(H_s\) (\(H_k\)), but is often assumed constant, \(D_s\) (\(D_k\)) is the average distance between sails (keels), and \(c_{ss}\) (\(c_{sk}\)) is the unobstructed atmospheric (oceanic) skin drag that would be attained in the absence of sails (keels) and with complete ice coverage, \(a_{ice}=1\).

Calculation of equations (11) – (15) requires that small-scale geometrical properties of the ice cover be related to average grid cell quantities already computed in the sea ice model. These intermediate quantities are briefly presented here and described in more detail in [70]. The sail height is given by

(16)\[H_{s} = \displaystyle 2\frac{v_{rdg}}{a_{rdg}}\left(\frac{\alpha\tan \alpha_{k} R_d+\beta \tan \alpha_{s} R_h}{\phi_r\tan \alpha_{k} R_d+\phi_k \tan \alpha_{s} R_h^2}\right),\]

and the distance between sails

(17)\[D_{s} = \displaystyle 2 H_s\frac{a_{i}}{a_{rdg}} \left(\frac{\alpha}{\tan \alpha_s}+\frac{\beta}{\tan \alpha_k}\frac{R_h}{R_d}\right),\]

where \(0<\alpha<1\) and \(0<\beta<1\) are weight functions, \(\alpha_{s}\) and \(\alpha_{k}\) are the sail and keel slope, \(\phi_s\) and \(\phi_k\) are constant porosities for the sails and keels, and we assume constant ratios for the average keel depth and sail height (\(H_k/H_s=R_h\)) and for the average distances between keels and between sails (\(D_k/D_s=R_d\)). With the assumption of hydrostatic equilibrium, the effective ice plus snow freeboard is \(H_{f}=\bar{h_i}(1-\rho_i/\rho_w)+\bar{h_s}(1-\rho_s/\rho_w)\), where \(\rho_i\), \(\rho_w\) and \(\rho_s\) are respectively the densities of sea ice, water and snow, \(\bar{h_i}\) is the mean ice thickness and \(\bar{h_s}\) is the mean snow thickness (means taken over the ice covered regions). For the melt pond edge elevation we assume that the melt pond surface is at the same level as the ocean surface surrounding the floes [16][17][18] and use the simplification \(H_p = H_f\). Finally to estimate the typical floe size \(L_A\), distance between floes, \(D_F\), and melt pond size, \(L_P\) we use the parameterizations of [41] to relate these quantities to the ice and pond concentrations. All of these intermediate quantities are available for output, along with Cdn_atm, Cdn_ocn and the ratio Cdn_atm_ratio_n between the total atmospheric drag and the atmospheric neutral drag coefficient.

We assume that the total neutral drag coefficients are thickness category independent, but through their dependance on the diagnostic variables described above, they vary both spatially and temporally. The total drag coefficients and heat transfer coefficients will also depend on the type of stratification of the atmosphere and the ocean, and we use the parameterization described in the Atmosphere section that accounts for both stable and unstable atmosphere–ice boundary layers. In contrast to the neutral drag coefficients the stability effect of the atmospheric boundary layer is calculated separately for each ice thickness category.

The transfer coefficient for oceanic heat flux to the bottom of the ice may be varied based on form drag considerations by setting the namelist variable fbot_xfer_type to Cdn_ocn; this is recommended when using the form drag parameterization. The default value of the transfer coefficient is 0.006 (fbot_xfer_type = ’constant’).