\documentclass[../main/NEMO_manual]{subfiles}

\begin{document}

\chapter{Land Ice Ocean interactions (ISF and ICB)}
\label{chap:LIO}

\chaptertoc

\paragraph{Changes record} ~\\

{\footnotesize
  \begin{tabularx}{\textwidth}{l||X|X}
    Release & Author(s) & Modifications \\
    \hline
    {\em   5.0} & {\em Pierre Mathiot and Katherine Hutchinson} & {\em update of the ice shelf iceberg}\\[2mm] 
    {\em   4.2} & {\em Pierre Mathiot} & {\em update of the ice shelf section (2019 developments)}\\[2mm]  
    {\em   4.0} & {\em ...} & {\em ...} \\
    {\em   3.6} & {\em ...} & {\em ...}
  \end{tabularx}
}

\clearpage

Land ice includes ice sheets, icebergs, and ice-shelves. Land ice / ocean interactions encompasses ice-shelves, glacier termini and icebergs melting, as well as surface and sub-glacial runoff from the ice sheet.\\
Land ice builds up through the accumulation of snowfall over Greenland and Antarctica. It influences the ocean through the melting of ice-shelves or glacier termini at the edge of the continents. 
This can also be via the calving of icebergs that slowly drift at the ocean surface and also by surface and subsurface subglacial runoff induced by ice sheet surface melting (seasonal variations).\\

Although land ice and sea ice bear some physical similarities, the modelling components to handle them are usually drastically different due to the differing scales of the problem and the distinctly separate processes at play. 
In NEMO, basal ice-shelf melting is handled through the ISF module \citep{mathiot.jenkins.ea_GMD17}, icebergs are handled through the ICB module \citep{marsh.ivchenko.ea_GMD15},
and surface runoff is handled through the runoff module. Subglacial runoff is not implemented mostly because there is no consolidated data set to prescribed it. \\

\section[Ice Shelf (ISF)]{Interaction with ice shelves (ISF)}
\label{sec:LIO_isf}

\begin{listing}
  \nlst{namisf}
  \caption{\forcode{&namisf}}
  \label{lst:namisf}
\end{listing}

The activation or not of the ice shelf, and the accompanying ocean interaction, is controled by the namelist variable \np{ln_isf}{ln\_isf} in \nam{isf}{isf}. 
The following interactions modes are available:
\begin{description}
   \item $\bullet$ representation of the ice shelf/ocean melting/freezing for open cavities (cav, \np{ln_isfcav_mlt}{ln\_isfcav\_mlt}, fig. \ref{fig:ISF}b,c).
   \item $\bullet$ parametrisation of the ice shelf/ocean melting/freezing for closed cavities (par, \np{ln_isfpar_mlt}{ln\_isfpar\_mlt}, fig. \ref{fig:ISF}d,e).
   \item $\bullet$ coupling with an ice sheet model (\np{ln_isfcpl}{ln\_isfcpl}).
\end{description}

The outcomes of the ISF module are the fresh water flux and the associated heat flux. A description and result of sensitivity tests to \np{ln_isfcav_mlt}{ln\_isfcav\_mlt} and \np{ln_isfpar_mlt}{ln\_isfpar\_mlt} are presented in \citet{mathiot.jenkins.ea_GMD17}. 

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.66\textwidth]{SBC_isf_v4.2}
  \caption[Ice shelf location and fresh water flux definition]{
    Illustration of the location where the fwf is injected and
    whether or not the fwf is interactive or not.}
  \label{fig:ISF}
\end{figure}


  \subsection[Ocean/Ice shelf fluxes in opened cavities (\textit{isfcav.F90})]{Ocean/Ice shelf fluxes in opened cavities (\protect\mdl{isfcav})}
 \label{sec:LIO_isfcav}
 
     \np{ln_isfcav_mlt}{ln\_isfcav\_mlt}\forcode{ = .true.} activates the ocean/ice shelf thermodynamic interactions at the ice shelf/ocean interface. 
     If \np{ln_isfcav_mlt}{ln\_isfcav\_mlt}\forcode{ = .false.}, thermodynamic interactions are desactivated but the ocean dynamics inside the cavity are still active.
     If  \np{ln_isfcav_mlt}{ln\_isfcav\_mlt} is activated, the user needs to make sure that the domain\_cfg.nc input file includes ice shelf cavities (\autoref{subsec:ISF_dom}).\\

     As part of the \mdl{isfcavmlt} module, 3 options are available to represent the ice-shelf/ocean fluxes at the interface:
     \begin{description}
        \item[\np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = 'spe'}]:
        The fresh water flux is specified by forcing fields \np{sn_isfcav_fwf}{sn\_isfcav\_fwf}. Convention of the input file is: positive toward the ocean (i.e. positive for melting and negative for freezing).
        The latent heat flux is derived from the fresh water flux. 
        The heat content flux is derived from the fresh water flux assuming a temperature set to freezing point in the top boundary layer (\np{rn_htbl}{rn\_htbl}).

        \item[\np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = 'oasis'}]:
        The \forcode{'oasis'} is a prototype of what could be a method to spread precipitation landing on the Antarctic ice sheet as ice shelf melt inside the open cavities when a coupled model (Atmosphere/Ocean) is used. 
        It has been tested in the IPSL model only. Therefore, it is highly recommended that other coupled models use it with caution. Feedbacks are very welcome. 
        
        \item[\np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = '2eq'}]:
        The heat flux and the fresh water flux due to ice shelf melting/freezing are parameterized following \citet{grosfeld1997}. 
        This formulation is based on a balance between the vertical diffusive heat flux across the ocean top boundary layer (\autoref{eq:ISOMIP1}) 
        and the latent heat due to melting/freezing (\autoref{eq:ISOMIP2}):

        \begin{equation}
        \label{eq:ISOMIP1}
        \mathcal{Q}_{h} = \rho c_{p} \gamma (T_{w} - T_{f})
        \end{equation}
        \begin{equation}
        \label{eq:ISOMIP2}
        q = \frac{-\mathcal{Q}_{h}}{L_{f}}
        \end{equation}
        
        where $\mathcal{Q}_{h}$($W.m^{-2}$) is the heat flux,q($kg.s^{-1}m^{-2}$) the fresh-water flux, 
        $L_{f}$ the specific latent heat, $T_w$ the temperature averaged over a boundary layer below the ice shelf (explained below), 
        $T_{f}$ the freezing point using  the  pressure  at  the  ice  shelf  base  and  the  salinity  of the water in the boundary layer, 
        and $\gamma$ the thermal exchange coefficient.

        \item[\np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = '3eq'}]:
        For realistic studies, the heat and freshwater fluxes are parameterized following \citep{jenkins2001}. This formulation is based on three equations: 
        a balance between the vertical diffusive heat flux across the boundary layer with variations in latent heat due to melting/freezing of ice and the vertical diffusive heat flux into the ice shelf (\autoref{eq:3eq1}); 
        a balance between the vertical diffusive salt flux across the boundary layer with the salt source or sink represented by melting/freezing (\autoref{eq:3eq2}); 
        and a linear equation for the freezing temperature of sea water \cite[\autoref{eq:3eq3} detailed in the linearisation coefficient of ]{asaydavis2016}:

        \begin{equation}
        \label{eq:3eq1}
        c_{p} \rho \gamma_{T} (T_{w}-T_{b}) = -L_{f} q - \rho_{i} c_{p,i} \kappa \frac{T_{s} - T_{b}}{h_{isf}}
        \end{equation}
        \begin{equation}
        \label{eq:3eq2}
        \rho \gamma_S (S_{w} - S_{b}) = (S_{i} - S_{b})q
        \end{equation}
        \begin{equation}
        \label{eq:3eq3}
        T_{b} = \lambda_{1} S_{b} + \lambda_{2} +\lambda_{3} z_{isf}
        \end{equation}

        where $T_{b}$ is the temperature at the interface, $S_b$ the salinity at the interface, $\gamma_{T}$ and $\gamma_{S}$ the exchange velocities for temperature and salt respectively, 
        $S_{i}$ the salinity of the ice (assumed to be 0), $h_{isf}$ the ice shelf thickness, $z_{isf}$ the ice shelf draft, $\rho_{i}$ the density of the iceshelf, 
        $c_{p,i}$ the specific heat capacity of the ice, $\kappa$ the thermal diffusivity of the ice 
        and $T_{s}$ is the atmospheric surface temperature (at the ice/air interface, assumed to be -20 \deg{C}). 
        The Liquidus slope ($\lambda_1$), the liquidus intercept ($\lambda_2$) and the Liquidus pressure coefficient ($\lambda_{3}$) 
        for TEOS80 and TEOS10 are described in \citet{asaydavis2016} and in \citet{jourdain2017} respectively.
        The linear system formed by \autoref{eq:3eq1}, \autoref{eq:3eq2} and the linearised equation for the freezing temperature of sea water (\autoref{eq:3eq3}) can be solved for $S_{b}$ or $T_{b}$. 
        Afterward, the freshwater flux ($q$) and the heat flux ($\mathcal{Q}_{h}$) can be computed.

     \end{description}
%
     \begin{table}[h]
        \centering
        \caption{Description of the parameters hard coded into the ISF module}
        \label{tab:isf}
        \begin{tabular}{|l|l|l|l|}
        \hline
        Symbol    & Description               & Value              & Unit               \\
        \hline
        $C_p$     & Ocean specific heat       & 3992               & $J.kg^{-1}.K^{-1}$ \\
        $L_f$     & Ice latent heat of fusion & $3.34 \times 10^5$ & $J.kg^{-1}$        \\
        $C_{p,i}$ & Ice specific heat         & 2000               & $J.kg^{-1}.K^{-1}$ \\
        $\kappa$  & Heat diffusivity          & $1.54 \times 10^{-6}$& $m^2.s^{-1}$     \\
        $\rho_i$  & Ice density               & 920                & $kg.m^3$           \\
        \hline
        \end{tabular}
     \end{table}

     Temperature and salinity used to compute the fluxes in \autoref{eq:ISOMIP1}, \autoref{eq:3eq1} and \autoref{eq:3eq2} are the average values for the top boundary layer \citep{losch_JGR08}. 
     The boundary layer thickness is defined by \np{rn_htbl}{rn\_htbl}.
     The fluxes and friction velocity are computed using the mean temperature, salinity and velocity in the first \np{rn_htbl}{rn\_htbl} m (see \mdl{isftbl}).
     Then, the fluxes are spread over the same thickness (ie over one or several cells).
     If \np{rn_htbl}{rn\_htbl} is larger than the top $e3t$, then there is no more direct feedback between the freezing point at the interface and the temperature of the top cell.
     This can therefore lead to a super-cool temperature in the top cell under melting conditions.
     If \np{rn_htbl}{rn\_htbl} is smaller than the top $e3t$, then the top boundary layer thickness is set to the top cell thickness.\\

     For each melt formula \forcode{'3eq'} or \forcode{'2eq'}, the turbulent transfer velocities ($\gamma^{T,S}$) between the ocean and the ice can be computed in any of the following 3 ways:

     \begin{description}
        \item[\np{cn_gammablk}{cn\_gammablk}\forcode{='spe'}]:
        The salt and heat transfert exchange velocities are constant and defined by:
\[
\gamma^{T} = \Gamma^{T}
\]
\[
\gamma^{S} = \Gamma^{S}
\] 
        This is the recommended formulation for ISOMIP.

	\item[\np{cn_gammablk}{cn\_gammablk}\forcode{='vel'}]:
        The salt and heat exchange coefficients are velocity dependent and defined as
\[
\gamma^{T} = \Gamma^{T} \times u_{*} 
\]
\[
\gamma^{S} = \Gamma^{S} \times u_{*}
\]
        where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_htbl}{rn\_htbl} meters). 
         
        See \citet{jenkins.nicholls.ea_JPO10} for all the details on this formulation. This is the recommended formulation for realistic applications and ISOMIP+/MISOMIP configurations.

	\item[\np{cn_gammablk}{cn\_gammablk}\forcode{'vel_stab'}]:
        The salt and heat exchange coefficients are velocity and stability dependent, and are defined as:
\[
\gamma^{T,S} = \frac{u_{*}}{\Gamma_{Turb} + \Gamma^{T,S}_{Mole}} 
\]
        where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_tbl}{rn\_htbl} meters),
        $\Gamma_{Turb}$ the contribution of the ocean stability and
        $\Gamma^{T,S}_{Mole}$ the contribution of the molecular diffusion.
        See \citet{holland.jenkins_JPO99} for all the details on this formulation. 
        This formulation has not been extensively tested in NEMO (and is thus not recommended).
     \end{description}
     
     In the formulation presented above, the transfert coeficient $\Gamma^{T}$ and $\Gamma^{S}$ are respectively defined by \np{rn_gammat0}{rn\_gammat0} and \np{rn_gammas0}{rn\_gammas0}. The definition of the exchange velocities $\gamma^{T,S}$ is done in the \mdl{isfcavgam} module.\\
     
     The ice shelf fresh water fluxes are implemented as a volume flux, same as for the runoff.
The fresh water flux addition due to the ice shelf melting is distributed uniformly over the top boundary layer and added to
the horizontal divergence (\textit{hdivn}) at each relevant depth level in the subroutine \rou{isf\_hdiv}, called from \mdl{divhor}.
As for the volume flux, the associated heat fluxes (latent heat and heat content) are distributed uniformly vertically over the top boundary layer in \mdl{traisf}.
See the runoff section \autoref{sec:SBC_rnf} for all the details about the divergence correction.\\
     
\subsection[Ocean/Ice shelf fluxes in parametrised cavities  (\textit{isfpar.F90})]{Ocean/Ice shelf fluxes in parametrised cavities (\protect\mdl{isfpar})}
     For a low resolution model, many ice shelves are too small to be explicitly represented. 
     It is therfore necessary to parametrise or specify the melt from these unresolved cavities (\np[=.true.]{ln_isfpar_mlt}{ln\_isfpar\_mlt}) in addition of some explicit cavities like Ross ice shelf and Filchner Ronne ice shelf, for exemple.
     The  \mdl{isfpar} module offers two options defined in the \mdl{isfparmlt} module:

  \begin{description}

     \item[\np{cn_isfpar_mlt}{cn\_isfpar\_mlt}\forcode{ = 'bg03'}]:
     The fwf and heat flux are computed using the \citet{beckmann.goosse_OM03} parameterisation of isf melting.
     The fluxes are distributed along the ice shelf front between the depth of the average grounding line (GL)
     (\np{sn_isfpar_zmax}{sn\_isfpar\_zmax}) and the base of the ice shelf along the calving front
     (\np{sn_isfpar_zmin}{sn\_isfpar\_zmin}) as in (\np{cn_isfpar_mlt}{cn\_isfpar\_mlt}\forcode{ = 'spe'}). The exchange velocity is specified by \np{rn_isfpar_bg03_gt0}{rn\_isfpar\_bg03\_gt0} and
     the effective melting length (\np{sn_isfpar_Leff}{sn\_isfpar\_Leff}) is read from a file.
     This parametrisation, however, has not been tested for some time and based on \citet{Favier2019} and more recently \citet{burgard_2022}, 
     the advice is that this parametrisation should not be used.

     \item[\np{cn_isfpar_mlt}{cn\_isfpar\_mlt}\forcode{ = 'spe'}]:
     The fwf ($q$) is read from \np{sn_isfpar_fwf}{sn\_isfpar\_fwf} and distributed along the ice shelf front between
     the average grounding line (GL) depth (\np{sn_isfpar_zmax}{sn\_isfpar\_zmax}) and
     the  depth of the base of the ice shelf at the front edge (\np{sn_isfpar_zmin}{sn\_isfpar\_zmin}). Convention of the input file is positive toward the ocean (i.e. positive for melting and negative for freezing).
     The heat flux ($Q_h$) is computed as $Q_h = q \times L_f$.

     \item[\np{cn_isfpar_mlt}{cn\_isfpar\_mlt}\forcode{ = 'oasis'}]:
     The \forcode{'oasis'} method is a prototype. It spreads the precipitation received by the Antarctic ice sheet as ice shelf melt inside the cavity when a coupled model (Atmosphere/Ocean) is used. 
        It has only been tested in the IPSL model and it is therefore recommended that other coupled models use it with caution. Feedbacks are very welcome. 

  \end{description}

As for the open cavity case, the ice shelf fresh water fluxes are implemented as a volume flux. 
In this case the volume and heat fluxes are uniformly distributed between \np{sn_isfpar_zmin}{sn\_isfpar\_zmin} and \np{sn_isfpar_zmax}{sn\_isfpar\_zmax}.\\

To conclude these two sections on ice shelf melt, it is worth noting the following:
\begin{description}

     \item \np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = '2eq'}, \np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = '3eq'} and \np{cn_isfpar_mlt}{cn\_isfpar\_mlt}\forcode{ = 'bg03'} computes a melt rate based on
the water mass properties, ocean velocities and depth.
The resulting fluxes are thus highly dependent of the model resolution (horizontal and vertical) and 
realism of the water masses onto the shelf.\\

     \item \np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = 'spe'} and \np{cn_isfpar_mlt}{cn\_isfpar\_mlt}\forcode{ = 'spe'} read the melt rate from a file located in the \np{cn_isfdir}{cn\_isfdir} directory.
Via these files the user has total control of the fwf forcing.
This can be useful if the water masses on the shelf are not realistic,
the resolution (horizontal/vertical) is too coarse to have realistic melting or
for studies where one needs to control heat and fw input. 
Note, however, that if the forcing is not consistent with the dynamics below, an unrealistic low water temperature can manifest.\\

  \end{description}

\subsection[Available outputs (\textit{isfdiags.F90})]{Available outputs (\protect\mdl{isfdiags})}
The following outputs are availables via XIOS:
\begin{description}
   \item[for parametrised cavities]:
      \begin{xmllines}
 <field id="isftfrz_par"     long_name="freezing point temperature in the parametrization boundary layer" unit="degC"     />
 <field id="fwfisf_par"      long_name="Ice shelf melt rate"                           unit="kg/m2/s"  />
 <field id="qoceisf_par"     long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     />
 <field id="qlatisf_par"     long_name="Ice shelf latent heat flux"                    unit="W/m2"     />
 <field id="qhcisf_par"      long_name="Ice shelf heat content flux of injected water" unit="W/m2"     />
 <field id="fwfisf3d_par"    long_name="Ice shelf melt rate"                           unit="kg/m2/s"  grid_ref="grid_T_3D" />
 <field id="qoceisf3d_par"   long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" />
 <field id="qlatisf3d_par"   long_name="Ice shelf latent heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" />
 <field id="qhcisf3d_par"    long_name="Ice shelf heat content flux of injected water" unit="W/m2"     grid_ref="grid_T_3D" />
 <field id="ttbl_par"        long_name="temperature in the parametrisation boundary layer" unit="degC" />
 <field id="isfthermald_par" long_name="thermal driving of ice shelf melting"          unit="degC"     />
      \end{xmllines}
   \item[for open cavities]:
      \begin{xmllines}
 <field id="isftfrz_cav"     long_name="freezing point temperature at ocean/isf interface"                unit="degC"     />
 <field id="fwfisf_cav"      long_name="Ice shelf melt rate"                           unit="kg/m2/s"  />
 <field id="qoceisf_cav"     long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     />
 <field id="qlatisf_cav"     long_name="Ice shelf latent heat flux"                    unit="W/m2"     />
 <field id="qhcisf_cav"      long_name="Ice shelf heat content flux of injected water" unit="W/m2"     />
 <field id="fwfisf3d_cav"    long_name="Ice shelf melt rate"                           unit="kg/m2/s"  grid_ref="grid_T_3D" />
 <field id="qoceisf3d_cav"   long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" />
 <field id="qlatisf3d_cav"   long_name="Ice shelf latent heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" />
 <field id="qhcisf3d_cav"    long_name="Ice shelf heat content flux of injected water" unit="W/m2"     grid_ref="grid_T_3D" />
 <field id="ttbl_cav"        long_name="temperature in Losch tbl"                      unit="degC"     />
 <field id="isfthermald_cav" long_name="thermal driving of ice shelf melting"          unit="degC"     />
 <field id="isfgammat"       long_name="Ice shelf heat-transfert velocity"             unit="m/s"      />
 <field id="isfgammas"       long_name="Ice shelf salt-transfert velocity"             unit="m/s"      />
 <field id="stbl"            long_name="salinity in the Losh tbl"                      unit="1e-3"     />
 <field id="utbl"            long_name="zonal current in the Losh tbl at T point"      unit="m/s"      />
 <field id="vtbl"            long_name="merid current in the Losh tbl at T point"      unit="m/s"      />
 <field id="isfustar"        long_name="ustar at T point used in ice shelf melting"    unit="m/s"      />
 <field id="qconisf"         long_name="Conductive heat flux through the ice shelf"    unit="W/m2"     />
      \end{xmllines}
\end{description}

%% =================================================================================================
\subsection[Ice sheet coupling (\textit{isfcpl.F90})]{Ice sheet coupling (\protect\mdl{isfcpl})}
\label{subsec:ISF_iscpl}
A coupling interface with any ice sheet model is available. The coupling is done offline through a file exchange at the restart step.
At every coupling step, the user needs to build a new domain\_cfg.nc file using the latest ice shelf draft provided by the ice sheet model.
The coupling frequency is usually one year. It is not shown yet in the literature that higher coupling frequency between an ocean model and an ice sheet model is needed.

At each restart step, the procedure is the following:

\begin{description}
\item[Step 1]: the ice sheet model sends a new bathymetry and ice shelf draft netcdf file to NEMO.
\item[Step 2]: a new domcfg.nc file is built using the DOMAINcfg tools.
\item[Step 3]: NEMO runs for a specific period and outputs the average melt rate over that period.
\item[Step 4]: the ice sheet model runs using the melt rate outputed in step 3.
\item[Step 5]: The process is repeated from step 1.
\end{description}

When the coupling with an ice sheet model is activated (\np[= .true.]{ln_iscpl}{ln\_iscpl}), the isf draft is assumed to be different at each restart step with
potentially some new wet/dry cells due to the ice sheet dynamics/thermodynamics.
The wetting and drying scheme, applied at the restart phase, is very simple. The 6 different possible cases for the tracers and ssh are:

\begin{description}
   \item[Thin a cell]: T/S/ssh are unchanged.
   \item[Enlarge  a cell]: T/S/ssh are unchanged.
   \item[Dry a cell]: Mask, T/S, U/V and ssh are set to 0.
   \item[Wet a cell]: Mask is set to 1, T/S is extrapolated from neighbour points, ssh is unchanged. If no neighbours exist, then T/S is extrapolated from old top cell values. 
   If no neighbours exist along i, j and k directions (both previous tests failed), T/S/ssh and mask are set to 0.
   \item[Dry a column]: mask, T/S, U/Vand ssh are set to 0.
   \item[Wet a column]: Mask is set to 1, T/S/ssh are extrapolated from neighbours. If no neighbour exists, then  T/S/ssh and mask set to 0.
\end{description}

The horizontal extrapolation to fill new cells with realistic values is called \np{nn_drown}{nn\_drown} times.
It means that if the grounding line retreats by more than \np{nn_drown}{nn\_drown} cells between 2 coupling steps,
the code will be unable to fill all the new wet cells properly and the model is likely to blow up at initialisation.
The default number is set up for the MISOMIP idealised experiments.

The coupling method described above will strongly affect the barotropic transport under an ice shelf when the geometry changes.
In order to keep the model stable, an adjustment of the dynamics at initialisation, after the coupling step, is needed. 
The idea behind this is to keep $\pd[\eta]{t}$ as it should be without a change in geometry at initialisation. 
This will prevent any strong velocities due to large pressure gradients. 
To do so, the local horizontal divergence is corrected to match what it was before the change in geometry. This is done before $\pd[\eta]{t}$ is computed in the first time step.\\

This coupling procedure is able to take into account grounding line and calving front migration.
Note, however, that it is a non-conservative process. 
This could thus lead to a trend in heat/salt content and volume.\\

In order to remove this possible trend and keep the conservation level as close to 0 as possible,
a simple conservation scheme is available with \np{ln_isfcpl_cons}{ln\_isfcpl\_cons}\forcode{ = .true.}.
The heat/salt/vol. gain/loss are diagnosed, along with the location.
A correction increment is computed and applied each time step during the model run.
The corrective increment is applied into the cells themselves (if it is a wet cell), the neigbouring cells or the closest wet cell (if the cell is now dry).

 \subsection[Ice shelf load (\textit{isfload.F90})]{Ice shelf load (\protect\mdl{isfload})}
 \label{subsec:ISF_hpg}

Ice shelf impose a load on the ocean. The main hypothesis to compute the ice shelf load is that the ice shelf is in an hydrostatic equilibrium. Therefore, beneath an ice shelf, the total pressure is the sum of the pressure due to the ice shelf load and
the pressure due to the ocean load. The ice shelf pressure is computed by integrating a reference density profile $\rho_{isf}$ from the surface to the base of the ice shelf. The local ocean pressure can thenThis can be formulated like this:

\begin{equation}
  \label{eq:ISF_hpg_isf_1}
    p(z)=\int_{0}^{z_{isf}}{\rho_{isf}gdz} + \int_{z_{isf}}^{z}{\rho g dz}  \\
\end{equation}

where $p(z)$ is the pressure at depth z, $\rho_{isf}$ a reference density profile, $\rho$ is the water density at depth $z$ and $z_{isf}$ is the ice shelf draft. $\rho_{isf}$ is assume to be the density of a water parcel at \np{rn_isfload_s}{rn\_isfload\_s} g/kg and \np{rn_isfload_t}{rn\_isfload\_t} \deg{C} that
corresponds to the water replaced by the ice shelf. These values are, for the time being, uniform in space and time (\np[='uniform']{cn_isfload}{cn\_isfload}).
A detailed description of this method is described in \citet{losch_JGR08}.\\

In a z∗ coordinate framework, the pressure gradient is formulated as suggested by \citet{adcroft.campin_OM04}:

\begin{equation}
  \label{eq:ISF_hpg_isf_2}
    \nabla_{z} p(z) = \nabla_{z^*} p(z) + \rho g \nabla_{z*} z \\
\end{equation}

The hydrostatic pressure gradient at a given level, k, (first term in \autoref{eq:ISF_hpg_isf_2}) is computed by adding the pressure gradient due to the ice shelf load (defined as the first term of \autoref{eq:ISF_hpg_isf_1}) to the vertical integral of the in situ density gradient along the model level from the surface to that level.\\
The only pressure gradient scheme compatible with the under ice shelf seas is (\np[=.true.]{ln_hpg_isf}{ln\_hpg\_isf}).
Outside of the ice shelf cavities, this pressure gradient scheme is exactly the same as traditional sco scheme (\np{ln_hpg_sco}{ln\_hpg\_sco}, \autoref{subsec:DYN_hpg_sco}).\\
 
 \subsection[Under ice shelf cavity geometry (\textit{domisf.F90})]{Under ice shelf cavity geometry (\protect\mdl{domisf})} 
 \label{subsec:ISF_dom}
 
%\begin{listing}
%  \nlst{namzgr_isf}
%  \caption{\forcode{&namzgr_isf}}
%  \label{lst:namzgr_isf}
%\end{listing}

To produce a simulation with ice shelf cavities opened, the user needs to make sure the bathymetry in the domain\_cfg.nc input file includes ice shelf cavities (and ice shelf draft information). The logical flag \np{ln_isfcav}{ln\_isfcav} in the \texttt{DOMAINcfg} namelist controls whether or not the ice shelf cavities are closed. All the options available to define the shape of the under ice shelf cavities are listed in \nam{zgr_isf}{zgr\_isf} (\texttt{DOMAINcfg} only, \autoref{lst:namzgr_isf}).\\

First of all, the tool makes sure that the ice shelf draft ($h_{isf}$) is compatible with the bathymetry (ice shelf draft being defined by \np{cn_draft}{cn\_draft} and bathymetry by \np{cn_bathy}{cn\_bathy} in \autoref{apdx:DOMCFG}). This is done in 3 steps :

\begin{description}

\item{First:} the position of the grounding line (separation between the floating ice shelf and the grounded ice sheet) is enforced. 
Where the difference between the bathymetry and the ice shelf draft is smaller than \np{rn_glhw_min}{rn\_glhw\_min}, the cell are grounded (ie masked). 
  This step is needed in order to take into account possible small mismatches between the ice shelf draft value and the bathymetry value (sources are coming from different grid, different data processes, rounding errors, ...).

\item{Secondly:} unrealistically thin ice shelves are corrected and $h_{isf}$ is set to a minimum value \np{rn_isfdep_min}{rn\_isfdep\_min}. If \np{rn_isfdep_min}{rn\_isfdep\_min} is smaller than the surface level thickness, \np{rn_isfdep_min}{rn\_isfdep\_min} is set to $e3t\_1d(1)$. 

\item{Finally:} The water column thickness  \np{rn_isfhw_min}{rn\_isfhw\_min} is enforced. 
  Where the water column thickness is lower than \np{rn_isfhw_min}{rn\_isfhw\_min}, the ice shelf draft is adjusted to match this criterion. 
  If for any reason this adjustment defies the minimum ice shelf draft allowed (\np{rn_isfdep_min}{rn\_isfdep\_min}), the cell is masked.
  
\end{description}

Once all these adjustments are made, the ice shelf draft and bathymetry are compatible. The user can, in addition, also remove the wet one-cell-wide channels under an ice shelf (\np{ln_isfchannel}{ln\_isfchannel}).  \\
 
After the definition of the ice shelf draft, the tool defines the top level. 
The compulsory criterion is that the water column needs at least 2 wet cells in the water column at U- and V-points.
To do so, if there is a one-cell-thick water column, the tools adjust the ice shelf draft to fulfil the requirement by digging into the ice shelf draft.\\

The process is the following:
\begin{description}
\item{step 1:} The top level is defined in the same way as the bottom level is defined.
\item{step 2:} The isolated grid points in the bathymetry are filled (as it is done in a domain without an ice shelf)
\item{step 3:} If the water column at a U- or V- point is one wet cell thick, the ice shelf draft is adjusted. So the actual top cell becomes fully open and the new
  top cell thickness is set to the minimum cell thickness allowed (following the same logic as for the bottom partial cell). This step is iterated 4 times to ensure the condition is fulfilled along the 4 lateral sides of the water column.
In the case of a steep slope and a shallow water column, it is likely that 2 adjacent wet water columns under the same ice shelf are disconnected (bathymetry lies above its neigbourging ice shelf draft). 
The option \np{ln_isfconnect}{ln\_isfconnect} allow the tool to force the connection between these 2 cells.
Some limits in meters or levels on the digging allowed by the tool are available (respectively, \np{rn_zisfmax}{rn\_zisfmax} or \np{rn_kisfmax}{rn\_kisfmax}).
This will force the connection of two parts of an ice shelf at the expense of a long vertical pipe to connect the cells at very different levels.
\item{step 5:} If after the previous step, the ice shelf draft is shallower than \np{rn_isfdep_min}{rn\_isfdep\_min}, it means that there is an incompatibility and the water column is closed.
\end{description}

In addition to this process, two extra options are available to the user for easthetic or stability reason. 
Despite careful setting of your ice shelf draft and bathymetry input file as well as setting described in \autoref{subsec:zgrisf_isfd}, some situations are unavoidable.
For exemple, chimneys within the ice shelf draft (ie wet cell surrounded horizontally by land cell) can be present. The circulation in these cells is not resolved. The user can decide to remove those features by activating \np{ln_isfcheminey}{ln\_isfcheminey}. \\
An other exemple is, if you setup your ice shelf draft and bathymetry to do ocean/ice sheet coupling, 
you may decide to fill the whole antarctic with a bathymetry and an ice shelf draft value (ice/bedrock interface depth when grounded). 
If you do so, the subglacial lakes will show up (Vostock for example). Another possibility is with coarse vertical resolution, some ice shelves could be cut in 2 parts: 
one connected to the main ocean and the other one closed which can be considered as a subglacial sea be the model. Keeping these features can lead to these situations:
\begin{description}
\item $\bullet$ For subglacial lakes in the case of very weak circulation (often the case), the only heat flux is the conductive heat flux through the ice sheet. 
  This will lead to constant freezing until water reaches -20C. 
  This is one of the deficiencies of the 3 equation melt formulation (for details on this formulation, see: \autoref{sec:LIO_isfcav}).
\item $\bullet$ In the case of coupling with an ice sheet model, 
  the ssh in the subglacial lakes and the main ocean could be very different (ssh initial adjustment for example), 
  and so if for any reason both a connected at some point, the model is likely to fail.\\
\end{description}
The namelist option \np{ln_isfsubgl}{ln\_isfsubgl} allow you to remove these subglacial lakes.


% =================================================================================================
\section[Icebergs (ICB)]{Interaction with icebergs (ICB)}
\label{sec:LIO_icb}

\begin{listing}
  \nlst{namberg}
  \caption{\forcode{&namberg}}
  \label{lst:namberg}
\end{listing}

Icebergs are modelled as lagrangian particles in \NEMO\ \citep{marsh.ivchenko.ea_GMD15}. The activation of the model is controled by \np[=.true.]{ln_icebergs}{ln\_icebergs}.
The physical behaviour of icebergs is controlled by equations as described in \citet{martin.adcroft_OM10} 
(Note that the authors kindly provided a copy of their code to act as a basis for the implementation thereof in \NEMO).

\subsection[Iceberg initialisation (\textit{icbclv.F90})]{Iceberg initialisation (\protect\mdl{icbclv})}

\subsubsection{Iceberg properties}
Icebergs are initially spawned into one of ten classes which have specific mass and thickness as
described in the \nam{berg}{berg} namelist: \np{rn_initial_mass}{rn\_initial\_mass} and \np{rn_initial_thickness}{rn\_initial\_thickness}.
Each class has an associated scaling (\np{rn_mass_scaling}{rn\_mass\_scaling}),
which is an integer representing how many icebergs of this class are being described as one lagrangian point
(this reduces the numerical problem of tracking every single iceberg). Each iceberg is generated with a pre determined density (\np{rn_rho_bergs}{rn\_rho\_bergs}). The iceberg area is determined based on the iceberg thickness, mass and density. The length and width are determined with a pre defined length over width ratio \np{rn_LoW_ratio}{rn\_LoW\_ratio}.

\subsubsection{Iceberg generation }

Two initialisation schemes are possible to generate the iceberg particules.
\begin{description}
\item {\np{nn_test_icebergs}{nn\_test\_icebergs}~$>$~0 (\protect\mdl{icbini})} This scheme is mostly used for developing, testing and debugging. 
  In this scheme, the value of \np{nn_test_icebergs}{nn\_test\_icebergs} represents the class of iceberg to generate
  (so between 1 and 10). \np{rn_test_box}{rn\_test\_box} defines a box using four numbers representing the corners of
  the geographical box: lonmin, lonmax, latmin, latmax. One test icebergs is released in every the cells of this box at the beginning of the run.
  This happens each time the timestep equals \np{nn_nit000}{nn\_nit000}.
\item {\np[=-1]{nn_test_icebergs}{nn\_test\_icebergs} (\protect\mdl{icbclv})} In this scheme, the model reads a calving file supplied in the \np{sn_icb}{sn\_icb} parameter.
  This should be a file with a field on the configuration grid
  representing ice accumulation rate at each model point in $km^3$ of ice / y.
  These should be ocean points adjacent to land where icebergs are known to calve.
  Most points in this input grid will have a value of zero.
  When the model runs, ice is accumulated at each grid point which has a non-zero source term.
  At each time step, a test is performed to see if there is enough ice mass to
  calve an iceberg of each class, in order (1 to 10). How much of the calving rate is then used to fill each iceberg class before releasing an iceberg is controled by the array \np{rn_distribution}{rn\_distribution}.
  The sum of the \np{rn_distribution}{rn\_distribution} array needs to be 1.
  Note that this is the initial mass multiplied by the number each particle represents (\ie\ the scaling).
  If there is enough ice, a new iceberg is spawned and the total available ice to generate a new iceberg of its category is reduced accordingly.
\end{description}

\subsection[Iceberg dynamics (\textit{icbdyn.F90})]{Iceberg dynamics (\protect\mdl{icbdyn})}
\label{subsec:ICB_icbdyn}

Iceberg dynamics are affected by multiple factors. The Coriolis force and the ocean, atmosphere and sea ice drags ($\vec{\tau_{o/a/i}}$) are the main drivers of the iceberg dynamics. In addition, the icebergs are further driven by the wave radiation force ($\vec{F_r}$) and a pressure gradient force ($\vec{F_p}$) proportional to the horizontal ssh gradient. 

The momentum balance for an iceberg of mass M is given by:

\begin{equation}
   \label{eq:icbmomentum}
        M\frac{d\vec{v}}{dt}=-Mf\times\vec{v}+\vec{\tau_a}+\vec{\tau_o}+\vec{\tau_i}+\vec{F_p}+\vec{F_r} 
\end{equation}

See \citet{martin.adcroft_OM10} for the full details on each term presented here. 
By default, the iceberg dynamics are computed using ocean surface variable (ssu, ssv) and the icebergs are not sensible to the bottom topography (only to land), whatever the iceberg draft may be. 
\citet{Merino_OM2016} developed an option to use vertical profiles of ocean currents instead (\np{ln_M2016}{ln\_M2016}). Full details on the sensitivity to this parameter are described in \citet{Merino_OM2016}. \\

Icebergs are considered as virtual and levitating over the ocean. There are therefore no interactions between icebergs, and therefore there is no limitation regarding the maximum number of icebergs within a grid cell. Furthermore, generally the simulated icebergs are not affected by a shallow bank. If, however, \np{ln_M2016}{ln\_M2016} is activated, \np{ln_icb_grd}{ln\_icb\_grd} engages an option to prevent thick icebergs to move across shallow banks (ie shallower than the iceberg draft). 
In the case of an interaction between one iceberg and the bedrock or the coastline, the iceberg bounces back to its original position.
Therefore, when using \np{ln_M2016}{ln\_M2016} and  \np{ln_icb_grd}{ln\_icb\_grd}, special care is needed in the calving file to be sure that the bathymetry is deep enough at the calving location so as to prevent the accumulation of grounded icebergs. 
If this is not possible, then the user needs to build a file describing the \forcode{maxclass} to prevent NEMO filling the icebergs classes that are too thick for the local bathymetry.\\

\subsection[Iceberg thermodynamics (\textit{icbthm.F90})]{Iceberg thermodynamics (\protect\mdl{icbthm})}
\label{subsec:ICB_icbthm}
The total iceberg melt and the associated latent heat flux are the only fields that feedback to NEMO. The iceberg mass balance is driven by three components (bottom melt ($M_{b}$), lateral melt ($M_{v}$) and erosion by surface waves ($M_{e}$) ):
\begin{equation}
  \label{eq:icbM}
        \rho\frac{d(LWH)}{dt}=\rho(-LWM_b-H(L+W)(M_e+M_v)) 
\end{equation}
where L, W and H are respectively the iceberg length, width and thickness.
The various melt formulations are from \citet{bigg.wadley_CRST97} and \citet{gladstone.bigg_JGR01}. A complete description is available in Appendix A of \citet{martin.adcroft_OM10}. Only the main characteristics are described below.\\

The wave erosion depends on the sea state ($S_{s}$) defined by a fit to the Beaufort scale, the ice concentration ($A_{i}$) and the sea surface temperature ($T_{surf}$).
\begin{equation}
  \label{eq:icbMe}
        M_{e} = \frac{1}{12}S_{s} ( 1 + \cos(\pi A^{3}_{i}) ) ( T_{surf} + 2)    \\
\end{equation}
By creating a notch at the floatation line by melting, the wave erosion acts to disintegrate the iceberg by breaking the overhanging slab of ice, thus creating small bits of ice which are assumed to propagate with their larger parent and thus delay flux into the ocean. 
The proportion of the wave erosion that should be represented as bits is controlled by \np{rn_bits_erosion_fraction}{rn\_bits\_erosion\_fraction}.\\

The lateral melt ($M_{v}$) is mainly driven by the buoyant convection along the side of the icebergs. \citet{eltahan.venkatesh_JOMAE97} empirically estimated to be:
\begin{equation}
   \label{eq:icbMv}
        M_{v} = 7.62 \times 10^{-3} T_{w}  + 1.29 \times 10^{-3} T_{w}^{2}   \\
\end{equation}
with $T_{w}$ being the mean lateral temperature. By default NEMO uses the sea surface temperature for $T_{w}$. \\

The bottom melt ($M_{b}$) is mainly driven by the turbulence created by the relative motion of water with respect to the iceberg velocity. Its empirical estimation is :
\begin{equation}
    \label{eq:icbMb}
        M_{b} = 0.58 |\vec{v}_{icb}-\vec{v}_{w}|^{0.8} \frac{T_{w} - T_{icb}}{L^{0.2}}   \\
\end{equation}
Where $\vec{v}_{icb}$ is the iceberg velocity, $\vec{v}_{w}$ the basal ocean velocity, $T_{w}$ the bottom temperature and $T_{icb}$ the estimated internal temperature of the icebergs \cite[$-4^{\circ}C$]{loset_JGR93}. By default, NEMO uses the ocean surface properties for $T_{w}$. To avoid melting even under freezing condition, the basal melt is set to 0 is the surface temperature is below the surface freezing point.\\

As for the iceberg dynamics, ocean vertical temperature profiles and ocean basal properties can be used to compute the lateral and basal melt respectively, instead of the surface temperature, by using \np{ln_M2016}{ln\_M2016}.\\

If the iceberg width is too small compared to its draft, the iceberg becomes unstable and capsizes. The iceberg thickness and width are then swapped. Such events are detected using the empirical criterion of Weeks and Mellor (1978) :
\begin{equation}
  \label{eq:icbrollover}
        W< \sqrt{0.92D^2 + 58.32D}  \\
\end{equation}
where W and D are the iceberg width and draft.

Finally, as mentioned above, iceberg melt rate is the only NEMO input from the iceberg model. The user can, however, disable it 
by using \np{ln_passive_mode}{ln\_passive\_mode}. In this case, NEMO will affect the icebergs, the melt will be computed and outputted, but the icebergs themselves will not have any impact on the NEMO solution.

\subsection[Iceberg mpp (\textit{icblbc.F90})]{Iceberg mpp (\protect\mdl{icblbc})}
\label{subsec:ICB_icblbc}
The management of the communication between the subdomains of the iceberg model is very different to the other components of NEMO, as particules (icebergs) are exchanged between subdomains.
This is done every iceberg time step, after the iceberg dynamics. 
The computation of the iceberg thermodynamics and dynamics are done in the inner domain and over the first half of the overlap region (this is necessary as otherwise the interpolation is not possible due to missing iceberg data). For now, NEMO does not make the most of the halos (2 by default now). 
This leads to the limitation on the cfl to be 0.4 (\np{rn_speed_limit}{rn\_speed\_limit}) because of interpolation of the ssh at 0.5 in this case. This is planned to be improve upon for the next version.

\subsection[Iceberg outputs (\textit{icbdia.F90})]{Iceberg diagnostics (\protect\mdl{icbdia})}
\label{subsec:ICB_icbdia}

There are three types of diagnostics available: a text diagnostics for the budget, xios diagnostics for the gridded output and the trajectory for diagnostics per iceberg.
\np{ln_bergdia}{ln\_bergdia} activates one (or more) type(s) of diagnostics. If set to false, no diagnostics will be outputted.

\subsubsection{Iceberg text diagnostics}
Extensive text diagnostics can be produced mostly for debug purposes.
Separate output files are maintained for human-readable iceberg information.
A separate file is produced for each processor (independent of \np{ln_ctl}{ln\_ctl}).
The amount of information is controlled by two integer parameters:
\begin{description}
\item [{\np{nn_verbose_level}{nn\_verbose\_level}}] takes a value between one and four. It
  represents an increasing level of verbosity.
\item [{\np{nn_verbose_write}{nn\_verbose\_write}}] is the number of time steps between each write.
\end{description}

\subsubsection{XIOS diagnostics}
The following outputs are available via XIOS:
\begin{xmllines}
  <field_group id="icbvar" domain_ref="grid_T"  >
    <field id="berg_melt"          long_name="icb melt rate of icebergs"                       unit="kg/m2/s"                    />
    <field id="berg_melt_hcflx"    long_name="icb heat flux to ocean due to melting heat content"   unit="J/m2/s"                />
    <field id="berg_melt_qlat"     long_name="icb heat flux to ocean due to melting latent heat"    unit="J/m2/s"                />
    <field id="berg_buoy_melt"     long_name="icb buoyancy component of iceberg melt rate"     unit="kg/m2/s"                    />
    <field id="berg_eros_melt"     long_name="icb erosion component of iceberg melt rate"      unit="kg/m2/s"                    />
    <field id="berg_conv_melt"     long_name="icb convective component of iceberg melt rate"   unit="kg/m2/s"                    />
    <field id="berg_virtual_area"  long_name="icb virtual coverage by icebergs"                unit="m2"                         />
    <field id="bits_src"           long_name="icb mass source of bergy bits"                   unit="kg/m2/s"                    />
    <field id="bits_melt"          long_name="icb melt rate of bergy bits"                     unit="kg/m2/s"                    />
    <field id="bits_mass"          long_name="icb bergy bit density field"                     unit="kg/m2"                      />
    <field id="berg_mass"          long_name="icb iceberg density field"                       unit="kg/m2"                      />
    <field id="calving"            long_name="icb calving mass input"                          unit="kg/s"                       />
    <field id="berg_floating_melt" long_name="icb melt rate of icebergs + bits"                unit="kg/m2/s"                    />
    <field id="berg_real_calving"  long_name="icb calving into iceberg class"                  unit="kg/s"     axis_ref="icbcla" />
    <field id="berg_stored_ice"    long_name="icb accumulated ice mass by class"               unit="kg"       axis_ref="icbcla" />
  </field_group>
\end{xmllines}


\subsubsection[Iceberg trajectory (\textit{icbtrj.F90})]{Iceberg trajectory (\protect\mdl{icbtrj})}
Iceberg trajectories (position as well as iceberg shapes and ocean properties at iceberg location) can also be outputted, and this is enabled by setting \np{nn_sample_rate}{nn\_sample\_rate}~$>$~0.
A non-zero value represents the frequency of timesteps that information is written to the output file.
These output files are in NETCDF format.
When running with multiple processes, each output file contains only the icebergs in the corresponding sub-domain.
Trajectory points are written in the order of their parent iceberg in the model's "linked list" of icebergs.
So care is needed to recreate data for individual icebergs,
since their trajectory data may be spread across multiple files. To rebuild it, the user can use a specific python tools located in \textit{tools/MISCELLANEOUS/icb\_pp.py}.

\section[Land Ice Runoff]{Land Ice Runoff}
\label{sec:LIO_rnf}

Runoff from the ice sheet is produced by the melt of snow and ice at the ice sheet surface and by the basal melt under the grounded area by friction and geothermal heating. It can be classified in two categories: the surface runoff and the sub-glacial runoff.\\

The surface melt is originated from the melt at the ice sheet surface that runs directly off of the edge of the ice shelf or glacier termini. NEMO is able to represent it as any other surface river runoff. The usage of this module is described in \autoref{sec:SBC_rnf}.\\

The subglacial runoff is produced by the melting at the interface bedrock-ice by geothermal heat and friction of the ice, and also by the percolation of the ice sheet surface melt from the surface to the bedrock through moulins. This kind of runoff joins the ocean at the grounding line of the ice sheet (line where ocean, bedrock and the grounded ice meet). 
Because of the lack of observations of such runoff, subglacial runoff is not represented in NEMO.\\

\subinc{\input{../../global/epilogue}}

\end{document}
