Changeset 11031

Timestamp:

2019-05-21T21:41:02+02:00 (5 years ago)

Author:

nicolasmartin

Message:

Update the chapters following the new keys for BibTeX entries

Location:

NEMO/trunk/doc/latex/SI3/subfiles

Files:

• abstract_foreword.tex (modified) (1 diff)
• ametsoc.bst (deleted)
• chap_domain.tex (modified) (2 diffs)
• chap_model_basics.tex (modified) (20 diffs)
• chap_output_diagnostics.tex (modified) (1 diff)
• chap_radiative_transfer.tex (modified) (8 diffs)
• chap_ridging_rafting.tex (modified) (6 diffs)
• chap_thermo.tex (modified) (2 diffs)
• chap_transport.tex (modified) (1 diff)
• introduction.tex (modified) (1 diff)

NEMO/trunk/doc/latex/SI3/subfiles/abstract_foreword.tex

@@ -10 +10 @@
 \chapter*{Abstract}
 
-SI$^3$ (Sea Ice modelling Integrated Initiative) is the sea ice engine of NEMO (Nucleus for European Modelling of the Ocean). It is adapted to regional and global sea ice and climate problems. It is intended to be a flexible tool for studying sea ice and its interactions with the other components of the Earth System over a wide range of space and time scales. SI$^3$ is based on the Arctic Ice Dynamics Joint EXperiment (AIDJEX) framework \citep{Coonetal74}, combining the ice thickness distribution framework, the conservation of horizontal momentum, an elastic-viscous plastic rheology, and energy-conserving halo-thermodynamics. Prognostic variables are the two-dimensional horizontal velocity field, ice volume, area, enthalpy, salt content, snow volume and enthalpy. In the horizontal direction, the model uses a curvilinear orthogonal grid. In the vertical direction, the model uses equally-spaced layers. In thickness space, the model uses thickness categories with prescribed boundaries. Various physical and numerical choices are available to describe sea ice physics. SI$^3$ is interfaced with the NEMO ocean engine, and, via the OASIS coupler, with several atmospheric general circulation models. It also supports two-way grid embedding via the AGRIF software.
+SI$^3$ (Sea Ice modelling Integrated Initiative) is the sea ice engine of NEMO (Nucleus for European Modelling of the Ocean). It is adapted to regional and global sea ice and climate problems. It is intended to be a flexible tool for studying sea ice and its interactions with the other components of the Earth System over a wide range of space and time scales. SI$^3$ is based on the Arctic Ice Dynamics Joint EXperiment (AIDJEX) framework \citep{coon_1974}, combining the ice thickness distribution framework, the conservation of horizontal momentum, an elastic-viscous plastic rheology, and energy-conserving halo-thermodynamics. Prognostic variables are the two-dimensional horizontal velocity field, ice volume, area, enthalpy, salt content, snow volume and enthalpy. In the horizontal direction, the model uses a curvilinear orthogonal grid. In the vertical direction, the model uses equally-spaced layers. In thickness space, the model uses thickness categories with prescribed boundaries. Various physical and numerical choices are available to describe sea ice physics. SI$^3$ is interfaced with the NEMO ocean engine, and, via the OASIS coupler, with several atmospheric general circulation models. It also supports two-way grid embedding via the AGRIF software.
 
 %\vspace{-40pt}

NEMO/trunk/doc/latex/SI3/subfiles/chap_domain.tex

@@ -75 +75 @@
 \nlst{namitd}
 
-Thickness space is discretized using $jl=1, ..., jpl$ thickness categories, with prescribed boundaries $hi\_max(jl-1),hi\_max(jl)$. Following \cite{Lipscomb01}, ice thickness can freely evolve between these boundaries. The number of ice categories $jpl$ can be adjusted from the namelist ($nampar$).
+Thickness space is discretized using $jl=1, ..., jpl$ thickness categories, with prescribed boundaries $hi\_max(jl-1),hi\_max(jl)$. Following \cite{lipscomb_2001}, ice thickness can freely evolve between these boundaries. The number of ice categories $jpl$ can be adjusted from the namelist ($nampar$).
 
 There are two means to specify the position of the thickness boundaries of ice categories. The first option (ln\_cat\_hfn) is to use a fitting function that places the category boundaries between 0 and 3$\overline h$, with $\overline h$ the expected mean ice thickness over the domain (namelist parameter rn\_himean), and with a greater resolution for thin ice (Fig. \ref{fig_dom_icecats}). More specifically, the upper limits for ice in category $jl=1, ..., jpl-1$ are:
@@ -100 +100 @@
 The other option (ln\_cat\_usr) is to specify category boundaries by hand using rn\_catbnd. The first category must always be thickner than rn\_himin (0.1 m by default).
 
-The choice of ice categories is important, because it constraints the ability of the model to resolve the ice thickness distribution. The latest study \citep{Massonnetetal18b} recommends to use at least 5 categories, which should include one thick ice with lower bounds at $\sim$4 m and $\sim$2 m for the Arctic and Antarctic, respectively, for allowing the storage of deformed ice.
+The choice of ice categories is important, because it constraints the ability of the model to resolve the ice thickness distribution. The latest study \citep{massonnet_2018} recommends to use at least 5 categories, which should include one thick ice with lower bounds at $\sim$4 m and $\sim$2 m for the Arctic and Antarctic, respectively, for allowing the storage of deformed ice.
 
 With a fixed number of cores, the cost of the model linearly increases with the number of ice categories. Using $jpl=1$ single ice category is also much cheaper than with 5 categories, but seriously deteriorates the ability of the model to grow and melt ice. Indeed, thin ice thicknes faster than thick ice, and shrinks more rapidly as well. When nn\_virtual\_itd=1 ($jpl$ = 1 only), two parameterizations are activated to compensate for these shortcomings. Heat conduction and areal decay of melting ice are adjusted to closely approach the 5 categories case.

NEMO/trunk/doc/latex/SI3/subfiles/chap_model_basics.tex

@@ -28 +28 @@
 
 \subsection{Scales, thermodynamics and dynamics}
-Because sea ice is much wider -- $\mathcal{O}$(100-1000 km) -- than thick -- $\mathcal{O}$(1 m) -- ice drift can be considered as purely horizontal: vertical motions around the hydrostatic equilibrium position are negligible. The same scaling argument justifies the assumption that heat exchanges are purely vertical\footnote{The latter assumption is probably less valid, because the horizontal scales of temperature variations are $\mathcal{O}$(10-100 m)}. It is on this basis that thermodynamics and dynamics are separated and rely upon different frameworks and sets of hypotheses: thermodynamics use the ice thickness distribution \citep{Thorndikeetal75} and the mushy-layer \citep{Worster92} frameworks, whereas dynamics assume continuum mechanics \citep[e.g.,][]{Lepparanta05}. Thermodynamics and dynamics interact by two means: first, advection impacts state variables; second, the horizontal momentum equation depends, among other things, on the ice state.
+Because sea ice is much wider -- $\mathcal{O}$(100-1000 km) -- than thick -- $\mathcal{O}$(1 m) -- ice drift can be considered as purely horizontal: vertical motions around the hydrostatic equilibrium position are negligible. The same scaling argument justifies the assumption that heat exchanges are purely vertical\footnote{The latter assumption is probably less valid, because the horizontal scales of temperature variations are $\mathcal{O}$(10-100 m)}. It is on this basis that thermodynamics and dynamics are separated and rely upon different frameworks and sets of hypotheses: thermodynamics use the ice thickness distribution \citep{thorndike_1975} and the mushy-layer \citep{worster_1992} frameworks, whereas dynamics assume continuum mechanics \citep[e.g.,][]{lepp_ranta_2011}. Thermodynamics and dynamics interact by two means: first, advection impacts state variables; second, the horizontal momentum equation depends, among other things, on the ice state.
 
 \subsection{Subgrid scale variations}
 
-Sea ice properties -- in particular ice thickness -- feature important changes at horizontal scales $\mathcal{O}$(1m) \citep{Thorndikeetal75}. An explicit representation of these variations is not and will not be -- at least in the next twenty years or so -- accessible to large-scale sea ice models. Yet important features, such as energy exchanges through the ice, quite non-linearly depend on ice thickness \citep{Maykut86}; whereas ice motion depends on the presence of open water, thin and thick ice at the very least, suggesting that subgrid-scale variations in ice properties must be accounted for, at least in a statistical fashion \citep{MaykutThorndike73}. 
+Sea ice properties -- in particular ice thickness -- feature important changes at horizontal scales $\mathcal{O}$(1m) \citep{thorndike_1975}. An explicit representation of these variations is not and will not be -- at least in the next twenty years or so -- accessible to large-scale sea ice models. Yet important features, such as energy exchanges through the ice, quite non-linearly depend on ice thickness \citep{maykut_1986}; whereas ice motion depends on the presence of open water, thin and thick ice at the very least, suggesting that subgrid-scale variations in ice properties must be accounted for, at least in a statistical fashion \citep{maykut_1973}. 
 
 %--------------------------------------------------------------------------------------------------------------------
@@ -49 +49 @@
 %--------------------------------------------------------------------------------------------------------------------
 
-The \textit{multi-category} framework \citep{MaykutThorndike73} addresses this issue by treating the ice thickness as an independent variable next to spatial coordinates and time, and introducing a thickness distribution\footnote{$g(h)$, termed the \textit{ice thickness distribution} is the density of probability of ice thickness \citep{Thorndikeetal75}.} $g(h)$ as the main prognostic model field. In the discrete world, the thickness distribution is converted into $L$ thickness categories. Ice thickness categories occupy a fraction of each grid cell, termed ice concentration ($a_l, l=1, 2, ..., L$), with specific thickness and properties. 
-
-The \textit{single-category} framework \citep{Hibler79} tackles the subgrid-scale issue by drastically simplifying the ice thickness distribution. The grid cell is divided into open water and sea ice characterized by a single ice concentration $A$ and mean thickness $H$. Single-category models (in particular LIM2) typically add parameterizations to represent the effects of unresolved ice thickness distribution on ice growth and melt \citep[see, e.g.][]{MellorKantha89,FichefetMaqueda97}.
+The \textit{multi-category} framework \citep{maykut_1973} addresses this issue by treating the ice thickness as an independent variable next to spatial coordinates and time, and introducing a thickness distribution\footnote{$g(h)$, termed the \textit{ice thickness distribution} is the density of probability of ice thickness \citep{thorndike_1975}.} $g(h)$ as the main prognostic model field. In the discrete world, the thickness distribution is converted into $L$ thickness categories. Ice thickness categories occupy a fraction of each grid cell, termed ice concentration ($a_l, l=1, 2, ..., L$), with specific thickness and properties. 
+
+The \textit{single-category} framework \citep{hibler_1979} tackles the subgrid-scale issue by drastically simplifying the ice thickness distribution. The grid cell is divided into open water and sea ice characterized by a single ice concentration $A$ and mean thickness $H$. Single-category models (in particular LIM2) typically add parameterizations to represent the effects of unresolved ice thickness distribution on ice growth and melt \citep[see, e.g.][]{mellor_1989,fichefet_1997}.
 
 SI$^3$ provides the choice between either a multi- or a single-category framework. The default mode is multi-category. The single-category mode can be activated by setting the number of categories ($jpl=1$) and by activating the virtual thickness distribution parameterizations (nn\_monocat=1).
@@ -70 +70 @@
  & Description & Value & Units & Ref \\ \hline
 $c_i$ (cpic) & Pure ice specific heat & 2067 & J/kg/K & ? \\
-$c_w$ (rcp) & Seawater specific heat & 3991 & J/kg/K & \cite{TEOS-10} \\
-$L$ (lfus) & Latent heat of fusion (0$^\circ$C) & 334000 & J/kg/K & \cite{BitzLipscomb99} \\
-$\rho_i$ (rhoic) & Sea ice density & 917 & kg/m$^3$ & \cite{BitzLipscomb99} \\
-$\rho_s$ (rhosn) & Snow density & 330 & kg/m$^3$ & \cite{MaykutUntersteiner71} \\
-$\mu$ (tmut) & Linear liquidus coefficient & 0.054 & $^\circ$C/(g/kg) & \cite{Assur58} \\
+$c_w$ (rcp) & Seawater specific heat & 3991 & J/kg/K & \cite{TEOS_2010} \\
+$L$ (lfus) & Latent heat of fusion (0$^\circ$C) & 334000 & J/kg/K & \cite{bitz_1999} \\
+$\rho_i$ (rhoic) & Sea ice density & 917 & kg/m$^3$ & \cite{bitz_1999} \\
+$\rho_s$ (rhosn) & Snow density & 330 & kg/m$^3$ & \cite{maykut_1971} \\
+$\mu$ (tmut) & Linear liquidus coefficient & 0.054 & $^\circ$C/(g/kg) & \cite{assur_1958} \\
 \hline
 \end{tabular}
@@ -85 +85 @@
 %
 
-Ice thermodynamics are formulated assuming that sea ice is covered by snow. Within each thickness category, both snow and sea ice are horizontally uniform, hence each thickness category has a specific ice thickness ($h_i^l$) and snow depth ($h_s^l$). Snow is assumed to be fresh, with constant density and thermal conductivity. Sea ice is assumed to be a \textit{mushy layer}\footnote{Mushy layers are two-phase, two-component porous media.} \citep{Worster92} of constant density, made of pure ice and brine in thermal equilibrium, related by a linear liquidus relationship \citep{BitzLipscomb99}. A vertically-averaged bulk salinity $S_l$ uniquely characterizes brine fraction for each thickness category, and changes through  time from a simple parametrization of brine drainage. The linear vertical salinity profile ($S^*_{kl}$) is reconstructed from the vertical mean \citep{Vancoppenolleetal09b}. The diffusion of heat affects the vertical temperature profile, discretized on a unique layer of snow and multiple ice layers (typically 2-5) for each category, whereas thermal properties depend on local brine fraction. Growth and melt rates are computed, also for each ice category. The choice of the main thermodynamic constants  is described in Tab. \ref{PhysicalConstants_table}.
+Ice thermodynamics are formulated assuming that sea ice is covered by snow. Within each thickness category, both snow and sea ice are horizontally uniform, hence each thickness category has a specific ice thickness ($h_i^l$) and snow depth ($h_s^l$). Snow is assumed to be fresh, with constant density and thermal conductivity. Sea ice is assumed to be a \textit{mushy layer}\footnote{Mushy layers are two-phase, two-component porous media.} \citep{worster_1992} of constant density, made of pure ice and brine in thermal equilibrium, related by a linear liquidus relationship \citep{bitz_1999}. A vertically-averaged bulk salinity $S_l$ uniquely characterizes brine fraction for each thickness category, and changes through  time from a simple parametrization of brine drainage. The linear vertical salinity profile ($S^*_{kl}$) is reconstructed from the vertical mean \citep{vancoppenolle_2009}. The diffusion of heat affects the vertical temperature profile, discretized on a unique layer of snow and multiple ice layers (typically 2-5) for each category, whereas thermal properties depend on local brine fraction. Growth and melt rates are computed, also for each ice category. The choice of the main thermodynamic constants  is described in Tab. \ref{PhysicalConstants_table}.
 
 
@@ -154 +154 @@
 \subsection{Dynamic formulation}
 
-The formulation of ice dynamics is based on the continuum approach. The latter holds provided the drift ice particles are much larger than single ice floes, and much smaller than typical gradient scales. This compromise is rarely achieved in practice \citep{Lepparanta05}. Yet the continuum approach generates a convenient momentum equation for the horizontal ice velocity vector $\mathbf{u}=(u,v)$, which can be solved with classical numerical methods (here, finite differences on the NEMO C-grid). The most important term in the momentum equation is internal stress. We follow the viscous-plastic (VP) rheological framework \citep{Hibler79}, assuming that sea ice has no tensile strength but responds to compressive and shear deformations in a plastic way. In practice, the elastic-viscous-plastic (EVP) technique of  \citep{Bouillonetal13} is used, more convient numerically than VP.  It is well accepted that the VP rheology and its relatives are the minimum complexity to get reasonable ice drift patterns \citep{Kreyscheretal00}, but fail at generating the observed deformation patterns \citep{Girardetal09}. This is a long-lasting problem: what is the ideal rheological model for sea ice and how it should be applied are still being debated \citep[see, e.g.][]{Weiss13}. 
+The formulation of ice dynamics is based on the continuum approach. The latter holds provided the drift ice particles are much larger than single ice floes, and much smaller than typical gradient scales. This compromise is rarely achieved in practice \citep{lepp_ranta_2011}. Yet the continuum approach generates a convenient momentum equation for the horizontal ice velocity vector $\mathbf{u}=(u,v)$, which can be solved with classical numerical methods (here, finite differences on the NEMO C-grid). The most important term in the momentum equation is internal stress. We follow the viscous-plastic (VP) rheological framework \citep{hibler_1979}, assuming that sea ice has no tensile strength but responds to compressive and shear deformations in a plastic way. In practice, the elastic-viscous-plastic (EVP) technique of  \citep{bouillon_2013} is used, more convient numerically than VP.  It is well accepted that the VP rheology and its relatives are the minimum complexity to get reasonable ice drift patterns \citep{kreyscher_2000}, but fail at generating the observed deformation patterns \citep{girard_2009}. This is a long-lasting problem: what is the ideal rheological model for sea ice and how it should be applied are still being debated \citep[see, e.g.][]{weiss_2013}. 
 
 %-------------------------------------------------------------------------------------------------------------------------
@@ -162 +162 @@
 %-------------------------------------------------------------------------------------------------------------------------
 
-We first present the essentials of the thickness distribution framework \citep{Thorndikeetal75}. Consider a given region of area $R$ centered at spatial coordinates $(\mathbf{x})$ at a given time $t$. $R$ could be e.g. a model grid cell. The ice thickness distribution $g(\mathbf{x},t, h)$ is introduced as follows:
+We first present the essentials of the thickness distribution framework \citep{thorndike_1975}. Consider a given region of area $R$ centered at spatial coordinates $(\mathbf{x})$ at a given time $t$. $R$ could be e.g. a model grid cell. The ice thickness distribution $g(\mathbf{x},t, h)$ is introduced as follows:
 \begin{linenomath}
 \begin{align}
@@ -168 +168 @@
 \end{align}
 \end{linenomath}
-where $dA(h,h+dh)$ is the surface fraction of all parts of $R$ with ice thickness between $h$ and $h+dh$. Using this definition, the spatial structure of ice thickness is lost (see Fig. \ref{fig_g_h}), and $h$ becomes an extra independent variable, next to spatial coordinates and time, that can be thought as random. $g$ is by definition normalized to 1. The conservation of area, expressed in terms of $g(h)$, is given by \citep{Thorndikeetal75}: 
+where $dA(h,h+dh)$ is the surface fraction of all parts of $R$ with ice thickness between $h$ and $h+dh$. Using this definition, the spatial structure of ice thickness is lost (see Fig. \ref{fig_g_h}), and $h$ becomes an extra independent variable, next to spatial coordinates and time, that can be thought as random. $g$ is by definition normalized to 1. The conservation of area, expressed in terms of $g(h)$, is given by \citep{thorndike_1975}: 
 \begin{linenomath}
 \begin{align}
@@ -199 +199 @@
 \end{align}
 \end{linenomath}
-Ice volume per area is the extensive counterpart for ice thickness, connected with volume through $h_l^i = v_l^i / a_l$. Evolution equations for extensive variables can be readily derived from equation \ref{eq:gt} by integration between thickness boundaries of the $l^{th}$ category \citep{Bitzetal01}. This applies to all model extensive variables (see Table \ref{GVariables_table}). For ice area, this reads:
+Ice volume per area is the extensive counterpart for ice thickness, connected with volume through $h_l^i = v_l^i / a_l$. Evolution equations for extensive variables can be readily derived from equation \ref{eq:gt} by integration between thickness boundaries of the $l^{th}$ category \citep{bitz_2001}. This applies to all model extensive variables (see Table \ref{GVariables_table}). For ice area, this reads:
 \begin{linenomath}
 \begin{align}
@@ -206 +206 @@
 \end{align}
 \end{linenomath}
-wher $\Theta^a_l$ refers to the effect of thermodynamics. Enthalpy is a particular case because it also has a vertical depth dependence $z$, which corresponds to $K$ vertical layers of equal thickness. The solution adopted here, following from \cite{ZhangRothrock01}, is that enthalpy from the individual layers are conserved separately. This is a practical solution, for lack of better.
+wher $\Theta^a_l$ refers to the effect of thermodynamics. Enthalpy is a particular case because it also has a vertical depth dependence $z$, which corresponds to $K$ vertical layers of equal thickness. The solution adopted here, following from \cite{zhang_2001}, is that enthalpy from the individual layers are conserved separately. This is a practical solution, for lack of better.
 
 One of the major actions of LIM is to resolve conservation equations for all extensive variables that characterize the ice state.  Let us now connect this detailed information with classical sea ice fields. The ice concentration $A$ and the ice volume per area\footnote{Ice volume per area is equivalent to the grid-cell averaged ice thickness.} $V_i$ (m) directly derive from $g$: 
@@ -232 +232 @@
 \end{align}
 \end{linenomath}
-where $m=\rho_i V_i + \rho_s V_s $ is the ice and snow mass per unit area, $\mathbf{u}$ is the ice velocity, $\mathbf{\sigma}$ is the internal stress tensor, $\mathbf{\tau}_a$ and $\mathbf{\tau}_w$ are the air and ocean stresses, respectively, $f$ is the Coriolis parameter, $\mathbf{k}$ is a unit vector pointing upwards, $g$ is the gravity acceleration and $\eta$ is the ocean surface elevation. The EVP approach used in LIM \citep{Bouillonetal13} gives the stress tensor as a function of the strain rate tensor $\dot{\mathbf{\epsilon}}$ and some of the sea ice state variables:
+where $m=\rho_i V_i + \rho_s V_s $ is the ice and snow mass per unit area, $\mathbf{u}$ is the ice velocity, $\mathbf{\sigma}$ is the internal stress tensor, $\mathbf{\tau}_a$ and $\mathbf{\tau}_w$ are the air and ocean stresses, respectively, $f$ is the Coriolis parameter, $\mathbf{k}$ is a unit vector pointing upwards, $g$ is the gravity acceleration and $\eta$ is the ocean surface elevation. The EVP approach used in LIM \citep{bouillon_2013} gives the stress tensor as a function of the strain rate tensor $\dot{\mathbf{\epsilon}}$ and some of the sea ice state variables:
 \begin{linenomath}
 \begin{align}
@@ -250 +250 @@
 
 Dynamical processes include the conservation of momentum, rheology, transport and mechanical redistribution.
-To resolve the momentum equation, atmospheric stress is taken either as forcing or from an atmospheric model, oceanic stress and sea surface elevation from the ocean model, the Coriolis term is trivial. The last term, the divergence of the internal stress tensor $\sigma$, is the most critical term in the momentum equation and requires a rheological formulation. The EVP approach used in LIM gives the stress tensor components as \citep{Bouillonetal13}:
+To resolve the momentum equation, atmospheric stress is taken either as forcing or from an atmospheric model, oceanic stress and sea surface elevation from the ocean model, the Coriolis term is trivial. The last term, the divergence of the internal stress tensor $\sigma$, is the most critical term in the momentum equation and requires a rheological formulation. The EVP approach used in LIM gives the stress tensor components as \citep{bouillon_2013}:
 \begin{linenomath}
 \begin{align}
@@ -257 +257 @@
 \end{align}
 \end{linenomath}
-where $\Delta$ is a particular measure of the deformation rate, $\Delta_{min}$ a parameter determining a smooth transition from pure viscous vlow ($\Delta<<\Delta_{min}$) to pure plastic flow ($\Delta >> \Delta_{min}$), and $e$ is a parameter giving the ratio between the maximum compressive stress and twice the maximum shear stress. In the pure plastic regime, the stress principal components should lie on the edge of an elliptical yield curve (Fig. \ref{fig_yield}). In the viscous regime, they are within the ellipse. The ice strength $P$ determines the plastic failure criterion and connects the momentum equation with the state of the sea ice. $P$ is not well constrained and must be parameterized. The heuristic option of \cite{Hibler79} was here adopted as a reference formulation:
+where $\Delta$ is a particular measure of the deformation rate, $\Delta_{min}$ a parameter determining a smooth transition from pure viscous vlow ($\Delta<<\Delta_{min}$) to pure plastic flow ($\Delta >> \Delta_{min}$), and $e$ is a parameter giving the ratio between the maximum compressive stress and twice the maximum shear stress. In the pure plastic regime, the stress principal components should lie on the edge of an elliptical yield curve (Fig. \ref{fig_yield}). In the viscous regime, they are within the ellipse. The ice strength $P$ determines the plastic failure criterion and connects the momentum equation with the state of the sea ice. $P$ is not well constrained and must be parameterized. The heuristic option of \cite{hibler_1979} was here adopted as a reference formulation:
 \begin{linenomath}
 \begin{align}
@@ -294 +294 @@
 % Parameters
 $P*$ (rn\_pstar) & ice strength thickness param. & 20000 & N/m2 & - \\
-$C$ (rn\_crhg) & ice strength concentration param. & 20 & - & \citep{Hibler79} \\
-$H^*$ (rn\_hstar) & maximum ridged ice thickness param. & 25 & m & \citep{Lipscombetal07} \\
-$p$ (rn\_por\_rdg) & porosity of new ridges & 0.3 & - & \citep{Lepparantaetal95} \\
+$C$ (rn\_crhg) & ice strength concentration param. & 20 & - & \citep{hibler_1979} \\
+$H^*$ (rn\_hstar) & maximum ridged ice thickness param. & 25 & m & \citep{lipscomb_2007} \\
+$p$ (rn\_por\_rdg) & porosity of new ridges & 0.3 & - & \citep{lepp_ranta_1995} \\
 $amax$ (rn\_amax) & maximum ice concentration & 0.999 & - & -\\
 $h_0$ (rn\_hnewice) & thickness of newly formed ice & 0.1 & m & - \\
@@ -311 +311 @@
 
 
-Transport connects the horizontal velocity fields and the rest of the ice properties. LIM assumes that the ice properties in the different thickness categories are transported at the same velocity. The scheme of \cite{Prather86}, based on the conservation of 0, 1$^{st}$ and 2$^{nd}$ order moments in $x-$ and $y-$directions,  is used, with some numerical diffusion if desired. Whereas this scheme is accurate, nearly conservative, it is also quite expensive since, for each advected field, five moments need to be advected, which proves CPU consuming, in particular when multiple categories are used. Other solutions are currently explored.
-
-The dissipation of energy associated with plastic failure under convergence and shear is accomplished by rafting (overriding of two ice plates) and ridging (breaking of an ice plate and subsequent piling of the broken ice blocks into pressure ridges). Thin ice preferentially rafts whereas thick ice preferentially ridges \citep{TuhkuriLensu02}. Because observations of these processes are limited, their representation in LIM is rather heuristic. The amount of ice that rafts/ridges depends on the strain rate tensor invariants (shear and divergence) as in \citep{FlatoHibler95}, while the ice categories involved are determined by a participation function favouring thin ice \citep{Lipscombetal07}. The thickness of ice being deformed ($h'$) determines whether ice rafts ($h'<$ 0.75 m) or ridges ($h'>$ 0.75 m), following \cite{Haapala00}. The deformed ice thickness is $2h'$ after rafting, and is distributed between $2h'$ and $2 \sqrt{H^*h'}$ after ridging, where $H^* = 25$ m \citep{Lipscombetal07}. Newly ridged ice is highly porous, effectively trapping seawater. To represent this, a prescribed volume fraction (30\%) of newly ridged ice \citep{Lepparantaetal95} incorporates mass, salt and heat are extracted from the ocean. Hence, in contrast with other models, the net thermodynamic ice production during convergence is not zero in LIM, since mass is added to sea ice during ridging. Consequently, simulated new ridges have high temperature and salinity as observed \citep{Hoyland02}. A fraction of snow (50 \%) falls into the ocean during deformation.
+Transport connects the horizontal velocity fields and the rest of the ice properties. LIM assumes that the ice properties in the different thickness categories are transported at the same velocity. The scheme of \cite{prather_1986}, based on the conservation of 0, 1$^{st}$ and 2$^{nd}$ order moments in $x-$ and $y-$directions,  is used, with some numerical diffusion if desired. Whereas this scheme is accurate, nearly conservative, it is also quite expensive since, for each advected field, five moments need to be advected, which proves CPU consuming, in particular when multiple categories are used. Other solutions are currently explored.
+
+The dissipation of energy associated with plastic failure under convergence and shear is accomplished by rafting (overriding of two ice plates) and ridging (breaking of an ice plate and subsequent piling of the broken ice blocks into pressure ridges). Thin ice preferentially rafts whereas thick ice preferentially ridges \citep{tuhkuri_2002}. Because observations of these processes are limited, their representation in LIM is rather heuristic. The amount of ice that rafts/ridges depends on the strain rate tensor invariants (shear and divergence) as in \citep{flato_1995}, while the ice categories involved are determined by a participation function favouring thin ice \citep{lipscomb_2007}. The thickness of ice being deformed ($h'$) determines whether ice rafts ($h'<$ 0.75 m) or ridges ($h'>$ 0.75 m), following \cite{haapala_2000}. The deformed ice thickness is $2h'$ after rafting, and is distributed between $2h'$ and $2 \sqrt{H^*h'}$ after ridging, where $H^* = 25$ m \citep{lipscomb_2007}. Newly ridged ice is highly porous, effectively trapping seawater. To represent this, a prescribed volume fraction (30\%) of newly ridged ice \citep{lepp_ranta_1995} incorporates mass, salt and heat are extracted from the ocean. Hence, in contrast with other models, the net thermodynamic ice production during convergence is not zero in LIM, since mass is added to sea ice during ridging. Consequently, simulated new ridges have high temperature and salinity as observed \citep{h_yland_2002}. A fraction of snow (50 \%) falls into the ocean during deformation.
 
 \section{Ice thermodynamics}
@@ -319 +319 @@
 
 \subsection{Transport in thickness space}
-Transport in thickness space describes how vertical growth and melt moves ice state variables among the different thicknesses at a velocity $f(h)$, the net ice growth/melt rate, which needs to be first computed.  In discretized form, this term moves ice properties between neighbouring categories. The linear remapping scheme of \cite{Lipscomb01} is used. This scheme is semi-lagrangian, second-order, is less diffusive and converges faster than other options.
+Transport in thickness space describes how vertical growth and melt moves ice state variables among the different thicknesses at a velocity $f(h)$, the net ice growth/melt rate, which needs to be first computed.  In discretized form, this term moves ice properties between neighbouring categories. The linear remapping scheme of \cite{lipscomb_2001} is used. This scheme is semi-lagrangian, second-order, is less diffusive and converges faster than other options.
 
 \subsection{Thermodynamic source and sink terms}
@@ -391 +391 @@
 %--------------------------------------------------------------------------------------------------------------------
 
-Based on these, brine fraction reduces to $\phi = -\mu S/T$ (see Fig. \ref{fig_thermal_properties}), where $\mu$ relates the freezing point of brine to salinity, and one can derive the specific enthalpy $q_m(S,T)$, defined as the energy required to warm and melt a unit control volume of sea ice at temperature $T$ (in Celsius) and salinity $S$ until 0$^\circ$ C, taken as a reference zero-energy level \citep{BitzLipscomb99,Schmidtetal04}:
+Based on these, brine fraction reduces to $\phi = -\mu S/T$ (see Fig. \ref{fig_thermal_properties}), where $\mu$ relates the freezing point of brine to salinity, and one can derive the specific enthalpy $q_m(S,T)$, defined as the energy required to warm and melt a unit control volume of sea ice at temperature $T$ (in Celsius) and salinity $S$ until 0$^\circ$ C, taken as a reference zero-energy level \citep{bitz_1999,schmidt_2004}:
 \begin{linenomath}
 \begin{align}
@@ -399 +399 @@
 where $c_i$ is pure ice specific heat, $L$ is latent heat of fusion at 0$^\circ$C, and $c_w$ is water specific heat. The first term expresses the warming of solid ice. The second term expresses internal change in brine fraction, which is often the largest because the Stefan number  ($c_i T/L$) is generally small. The last term gives the warming of the remaining water from $T_{fr} = -\mu S$ until 0$^\circ$C. Similar, but simpler and linear expressions for snow and water can be derived.
 
-The second overarching aspect is that all growth and melt processes must be calculated consistently with the enthalpy formulation. Energetics of phase transitions are handled using the formalism of \cite{Schmidtetal04}. For each phase transition, initial and final states (temperature and salinity) are defined, and the ice-to-ocean mass flux to the ice $F_m$ (kg/s) relates to the energy gain or loss $\Delta Q$ through:
+The second overarching aspect is that all growth and melt processes must be calculated consistently with the enthalpy formulation. Energetics of phase transitions are handled using the formalism of \cite{schmidt_2004}. For each phase transition, initial and final states (temperature and salinity) are defined, and the ice-to-ocean mass flux to the ice $F_m$ (kg/s) relates to the energy gain or loss $\Delta Q$ through:
 \begin{linenomath}
 \begin{align}
@@ -418 +418 @@
 \item the sensible heat flux from the ocean to the sea ice ($A.F_w$)
 \end{itemize}
-Other contributions are not assumed not to contribute. The ocean-to-ice sensible heat flux is formulated the bulk formula of \citep{McPhee92}.
+Other contributions are not assumed not to contribute. The ocean-to-ice sensible heat flux is formulated the bulk formula of \citep{mcphee_1992}.
 
 % Ice growth
-If $B^{opw}$ is such that the SST would decrease below the freezing point, the remainder of the heat is used to form new ice. The heat loss is converted into a volume of new ice $v_0$. The thickness $h_0$ of the new ice grown during a sea ice time step depends on unresolved small-scale currents and waves and is prescribed. The fraction $a_0=v_0/h_0$ is computed accordingly. The salinity of this new ice $S_0$ is given by the salinity-thickness empirical relationship of \cite{Kovacs96}. The temperature assumed for this new ice is the local freezing point. If by contrast $B^{opw}$ is positive and there still is ice in the grid cell, then $B^{opw}$ is directly redirected to bottom melting. This argument follows from \cite{MaykutMcPhee95}, who found that most of solar heat absorbed in the surface waters is converted into melting. In practise, this prevents the SST to be above freezing as long ice is present.
+If $B^{opw}$ is such that the SST would decrease below the freezing point, the remainder of the heat is used to form new ice. The heat loss is converted into a volume of new ice $v_0$. The thickness $h_0$ of the new ice grown during a sea ice time step depends on unresolved small-scale currents and waves and is prescribed. The fraction $a_0=v_0/h_0$ is computed accordingly. The salinity of this new ice $S_0$ is given by the salinity-thickness empirical relationship of \cite{kovacs_1996}. The temperature assumed for this new ice is the local freezing point. If by contrast $B^{opw}$ is positive and there still is ice in the grid cell, then $B^{opw}$ is directly redirected to bottom melting. This argument follows from \cite{maykut_1995}, who found that most of solar heat absorbed in the surface waters is converted into melting. In practise, this prevents the SST to be above freezing as long ice is present.
 
 $B^{opw}$ can be seen as a predictor of the heat budget of the first ocean level. As such, it only helps to compute new ice formation and the extra bottom melt in summer, but is not part of the conservation of heat in the model. To ensure heat conservation, the heat effectively contributing to changing sea ice is removed from the non-solar flux sent to the ocean. This includes: (i) the heat loss used for ice formation, (ii) the heat gain used to melt ice, and (iii) the sensible heat given by the ocean to the ice. Finally, because ice dynamics are not able to maintain the small amount of open water that is observed, a maximum ice fraction ($amax, <1$) is prescribed.
@@ -435 +435 @@
 \end{align}
 \end{linenomath}
-which state that the local change in enthalpy is given by the divergence of the vertical conduction ($F_c=-k(S,T) \partial T / \partial z$) and radiation ($F_r$) fluxes. $\rho$ is the density of ice or snow. Re-expressed as a function of temperature, this becomes the heat diffusion equation.  This equation is non-linear in $T$, because of $q$ and $k$, and its main specificity is that internal melting requires large amounts of energy near the freezing point. The thermal conductivity is formulated following \cite{Pringleetal07}, empirically accounting for the reduction of thermal conductivity at large brine fractions. 
+which state that the local change in enthalpy is given by the divergence of the vertical conduction ($F_c=-k(S,T) \partial T / \partial z$) and radiation ($F_r$) fluxes. $\rho$ is the density of ice or snow. Re-expressed as a function of temperature, this becomes the heat diffusion equation.  This equation is non-linear in $T$, because of $q$ and $k$, and its main specificity is that internal melting requires large amounts of energy near the freezing point. The thermal conductivity is formulated following \cite{pringle_2007}, empirically accounting for the reduction of thermal conductivity at large brine fractions. 
 
 At the ice base, we assume that the temperature is at the local freezing point.  Ice grows or melt if the heat balance between the oceanic sensible heat flux ($F_w$) and internal conduction is negative or positive.
@@ -447 +447 @@
 where $Q^{sr}$ and $Q^{ns}$ are the net downwelling atmospheric solar and non-solar flux components. If the solution of this equation without melting gives a surface temperature ($T_{su}$) below 0$^\circ$ C, then there is no melting and the heat available for surface melting $Q_{sum}=0$.  Otherwise $T_{su}$ is capped at 0$^\circ$ C and $Q_{sum}$ is calculated as a residual.
 
-\textbf{Radiation}. Radiation contributes to the surface and internal heat budget. The radiative transfer scheme is currently basic, composed of surface albedo, transmission through the ice interior and attenuation with vertical depth. The albedo is computed empirically as a function of ice thickness, snow depth and surface temperature, using a reformulation of the parameterization of \cite{ShineHenderson85}. When snow is present, all the absorbed radiation is transformed into sensible heat available for conduction or melting. Over snow-free ice, a fraction of solar radiation is transmitted below the surface and attenuates exponentially with depth, until it reaches the base of the ice.
+\textbf{Radiation}. Radiation contributes to the surface and internal heat budget. The radiative transfer scheme is currently basic, composed of surface albedo, transmission through the ice interior and attenuation with vertical depth. The albedo is computed empirically as a function of ice thickness, snow depth and surface temperature, using a reformulation of the parameterization of \cite{shine_1985}. When snow is present, all the absorbed radiation is transformed into sensible heat available for conduction or melting. Over snow-free ice, a fraction of solar radiation is transmitted below the surface and attenuates exponentially with depth, until it reaches the base of the ice.
 
 \textbf{Growth and melt processes}.  Snow grows from precipitation and loses mass from melting and snow-ice conversion once the snow base is below sea level. Sea ice grows and melts by various means. Ice forms by congelation or melt at the base, can melt at the surface and form from snow-to-ice conversion at the snow-ice interface if the latter is below sea level. Some new ice is also added to the system when seawater is trapped into newly formed pressure ridges.
 
-\textbf{Salt dynamics}. Bulk salinity is empirically parameterized, as a function of salt uptake during growth, gravity drainage and flushing. The shape of the vertical profile depends on the bulk salinity \citep{Vancoppenolleetal09b}.
-
-\textbf{Single-category parameterizations}. If the single-category representation is adopted, then two parameterizations can be activated, following \citep{FichefetMaqueda97}. First, the thermal conductivity of both ice and snow is multiplied by a factor $>1$ accounting for the unresolved thin ice, effectively increasing the ice growth rate. Second, to account for the loss of thin ice in summer, the ice concentration is reduced in proportion to the loss of ice thickness. Both parameterizations have been tuned to match the results in multi-category mode.
+\textbf{Salt dynamics}. Bulk salinity is empirically parameterized, as a function of salt uptake during growth, gravity drainage and flushing. The shape of the vertical profile depends on the bulk salinity \citep{vancoppenolle_2009}.
+
+\textbf{Single-category parameterizations}. If the single-category representation is adopted, then two parameterizations can be activated, following \citep{fichefet_1997}. First, the thermal conductivity of both ice and snow is multiplied by a factor $>1$ accounting for the unresolved thin ice, effectively increasing the ice growth rate. Second, to account for the loss of thin ice in summer, the ice concentration is reduced in proportion to the loss of ice thickness. Both parameterizations have been tuned to match the results in multi-category mode.
 
 \end{document}

NEMO/trunk/doc/latex/SI3/subfiles/chap_output_diagnostics.tex

@@ -16 +16 @@
 \section{SIMIP diagnostics}
 
-The SIMIP protocol \citep{Notzetal16} was designed for CMIP6, to standardize sea ice model outputs in climate simulations. We tried to follow the data request as closely as possible. Outputs are in most cases directly managed with XIOS2 in \textbf{limwri.F90}, but not always. In the code, output fields keep their native LIM reference name. 
+The SIMIP protocol \citep{notz_2016} was designed for CMIP6, to standardize sea ice model outputs in climate simulations. We tried to follow the data request as closely as possible. Outputs are in most cases directly managed with XIOS2 in \textbf{limwri.F90}, but not always. In the code, output fields keep their native LIM reference name. 
 
 A corresponding entry exists in \textbf{field\_def\_nemo-lim.xml}, where fields are given their SIMIP specifications (standard name, long name, units). At the end of the file the fields are gathered in the field groups \textbf{SIday\_fields}, \textbf{SImon\_fields} and \textbf{SImon\_scalar} for separation of the daily (SIday) and monthly (SImon) requests. 

NEMO/trunk/doc/latex/SI3/subfiles/chap_radiative_transfer.tex

@@ -36 +36 @@
 %--------------------------------------------------------------------------------------------------------------------
 
-Solar radiation in the snow-ice system is represented following the principles of \cite{MaykutUntersteiner71}, see Fig.\ref{fig_radiative_transfer}, using a unique band of solar radiation. Incident solar radiation (W/m$^2$, counted per unit ice area - not per grid cell area) is specified in the SBC routines and is a priori category dependent, because multiple atmosphere-surface reflexions are frequent in polar regions imply that incident radiation depends on the surface albedo and therefore surface state. 
+Solar radiation in the snow-ice system is represented following the principles of \cite{maykut_1971}, see Fig.\ref{fig_radiative_transfer}, using a unique band of solar radiation. Incident solar radiation (W/m$^2$, counted per unit ice area - not per grid cell area) is specified in the SBC routines and is a priori category dependent, because multiple atmosphere-surface reflexions are frequent in polar regions imply that incident radiation depends on the surface albedo and therefore surface state. 
 
 Net solar radiation qsr\_ice(i,j,l) is obtained by substracting the reflected part of the incident radiation using the surface albedo $\alpha(i,j,l)$, parameterized as a function of environmental conditions. 
@@ -50 +50 @@
 The surface albedo determines the amount of solar radiation that is reflected by the ice surface, hence also net solar radiation. The philosophy of the parameterization of surface albedo is the following: each ice category has its own albedo value $\alpha(i,j,l)$, determined as a function of cloud fraction, ice thickness, snow depth, melt pond fraction and depth, using observation-based empirical fits. 
 
-The original \cite{ShineHenderson85} parameterization had a few inconsistencies and flaws that the revisited parameterization described hereafter fixes. In particular, the dependencies of albedos on ice thickness, snow depth and cloud fraction have been revised in the light of recent observational constraints \citep{Brandtetal05,GrenfellPerovich04}. In addition, the asymptotic properties of albedo are better specified and now fully consistent with oceanic values. Finally, the effect of melt ponds has been included \citep{Lecomteetal15}.
-
-The user has control on 5 reference namelist values, which describe the asymptotic values of albedo of snow and ice for dry and wet conditions, as well as the deep ponded-ice albedo. Observational surveys, in particular during SHEBA in the Arctic \citep{Perovichetal02alb} and further additional experiments \citep{GrenfellPerovich04}, as well as by \cite{Brandtetal05} in the Antarctic, have provided relatively strong constraints on the surface albedo. In this context, the albedo can hardly be used as the main model tuning parameter, at least outside of these observation-based bounds (see namalb for reference values).
+The original \cite{shine_1985} parameterization had a few inconsistencies and flaws that the revisited parameterization described hereafter fixes. In particular, the dependencies of albedos on ice thickness, snow depth and cloud fraction have been revised in the light of recent observational constraints \citep{brandt_2005,grenfell_2004}. In addition, the asymptotic properties of albedo are better specified and now fully consistent with oceanic values. Finally, the effect of melt ponds has been included \citep{lecomte_2015}.
+
+The user has control on 5 reference namelist values, which describe the asymptotic values of albedo of snow and ice for dry and wet conditions, as well as the deep ponded-ice albedo. Observational surveys, in particular during SHEBA in the Arctic \citep{perovich_2002} and further additional experiments \citep{grenfell_2004}, as well as by \cite{brandt_2005} in the Antarctic, have provided relatively strong constraints on the surface albedo. In this context, the albedo can hardly be used as the main model tuning parameter, at least outside of these observation-based bounds (see namalb for reference values).
 
 \nlst{namalb}
@@ -64 +64 @@
 \vspace{0cm}
 \includegraphics[height=10cm,angle=-00]{albedo_cloud_correction}
-\caption{Albedo correction $\Delta \alpha$ as a function of overcast sky (diffuse light) albedo $\alpha_os$, from field observations \cite[][their Table 3]{GrenfellPerovich04} (squares) and 2nd-order fit (Eq. \ref{eq_albedo_cloud_correction}). Red squares represent the irrelevant data points excluded from the fit. For indication, the amplitude of the correction used in the ocean component is also depicted (blue circle).}
+\caption{Albedo correction $\Delta \alpha$ as a function of overcast sky (diffuse light) albedo $\alpha_os$, from field observations \cite[][their Table 3]{grenfell_2004} (squares) and 2nd-order fit (Eq. \ref{eq_albedo_cloud_correction}). Red squares represent the irrelevant data points excluded from the fit. For indication, the amplitude of the correction used in the ocean component is also depicted (blue circle).}
 % ocean uses 0.06 for overcast sky (Payne 74) and Briegleb and Ramanathan parameterization
 \label{fig_albedo_cloud_correction}
@@ -72 +72 @@
 %--------------------------------------------------------------------------------------------------------------------
 
-Because the albedo is not an intrinsic optical property, it depends on the type of light (diffuse of direct), which is practically handled by weighting the clear (cs) and overcast (os) skies values by cloud fraction $c(i,j)$ \citep{FichefetMaqueda97}:
+Because the albedo is not an intrinsic optical property, it depends on the type of light (diffuse of direct), which is practically handled by weighting the clear (cs) and overcast (os) skies values by cloud fraction $c(i,j)$ \citep{fichefet_1997}:
 \begin{equation}
 \alpha(i,j,l) = [ 1 - c(i,j) ] \cdot \alpha_{cs} (i,j,l) + c (i,j) \cdot \alpha_{os}(i,j,l).
 \end{equation}
-For concision, we drop the spatial and category indices hereafter. \cite{GrenfellPerovich04} observations at Point Barrow, on the Alaskan Coast, suggest that clear and overcast sky albedos are directly related through
+For concision, we drop the spatial and category indices hereafter. \cite{grenfell_2004} observations at Point Barrow, on the Alaskan Coast, suggest that clear and overcast sky albedos are directly related through
 \begin{equation}
 \alpha_{cs} = \alpha_{os} - \Delta \alpha(\alpha_{os}).
@@ -110 +110 @@
 The surface fractions $f_{ice}$, $f_{snw}$ and $f_{pnd}$ are currently crudely parameterized: if snow is present ($h_s>0$), then $f_{snw}=1$ and $f_{ice}=f_{pnd}=0$. In the absence of snow, $f_{pnd}$ is either specified or calculated (depending on melt pond options in nampnd), and $f_{ice}=1.-f_{pnd}$. Admittedly, more refined parameterizations of $f_{snw}$ could improve the realism of the model. Note finally that the dependence of surface albedo on the presence of melt ponds can be included or not (namelist parameter ln\_pnd\_alb). If the latter is set to false, $f_{pnd}$ is always assumed zero in the albedo computations.
 
-Works by \cite{Brandtetal05} and references therein, indicate that the dependence of the albedo of bare ice on ice thickness depends is linear/logarithmic/constant from thin to thick ice. Hence, the following expressions capture  the essence of their works:
+Works by \cite{brandt_2005} and references therein, indicate that the dependence of the albedo of bare ice on ice thickness depends is linear/logarithmic/constant from thin to thick ice. Hence, the following expressions capture  the essence of their works:
 \begin{eqnarray}
 \alpha_{ice} = 
@@ -129 +129 @@
 values that are to be specified from the namelist.
 
-\cite{GrenfellPerovich04} suggest that the dependence of surface albedo on snow depth is exponential,
+\cite{grenfell_2004} suggest that the dependence of surface albedo on snow depth is exponential,
 \begin{eqnarray}
 \alpha_{snw} = \alpha_{snw}^{\infty} - ( \alpha_{snw}^{\infty} - \alpha_{ice} ) * exp( -h_s / h_s^{ref} ),
@@ -143 +143 @@
 values that are to be specified from the namelist.
 
-Based on ideas developed from melt ponds on continental ice \citep{ZuoOerlemans96}, the albedo of ponded ice was proposed to follow \citep{Lecomteetal11}:
+Based on ideas developed from melt ponds on continental ice \citep{zuo_1996}, the albedo of ponded ice was proposed to follow \citep{lecomte_2011}:
 \begin{eqnarray}
 \alpha_{pnd} = \alpha_{dpnd} - ( \alpha_{dpnd} -\alpha_{ice} ) \cdot exp( -h_{pnd} / 0.05 )
 \end{eqnarray}
-$\alpha_{dpnd}$ is a namelist parameter. \cite{EbertCurry93} also use such dependency for their multi-spectral albedo.
+$\alpha_{dpnd}$ is a namelist parameter. \cite{ebert_1993} also use such dependency for their multi-spectral albedo.
 
 The dependencies of surface albedo on ice thickness, snow depth and pond depth are illustrated in Fig. \ref{fig_albedo_dependencies}.
@@ -153 +153 @@
 \subsection{Transmission below the snow/ice surface}
 
-The transmitted solar radiation below the surface is represented following \cite{FichefetMaqueda97} and \cite{MaykutUntersteiner71}:
+The transmitted solar radiation below the surface is represented following \cite{fichefet_1997} and \cite{maykut_1971}:
 \begin{eqnarray}
 qtr\_ice\_top(i,j,l) = i_o(i,j) qsr\_ice(i,j,l),
 \end{eqnarray}
-where $i_o=0$ in presence of snow, and depends on cloud fraction otherwise, based on works of \cite{GrenfellMaykut77}. This parameterization needs to be re-evaluated and likely updated.
+where $i_o=0$ in presence of snow, and depends on cloud fraction otherwise, based on works of \cite{grenfell_1977}. This parameterization needs to be re-evaluated and likely updated.
 
 \subsection{Attenuation and transmission below the ice/ocean interface}

NEMO/trunk/doc/latex/SI3/subfiles/chap_ridging_rafting.tex

@@ -29 +29 @@
 \end{equation}
 
-Second, we must say something about volume conservation, and this will be done more specifically later. Following \cite{Thorndikeetal75}, we separate the $\Psi^X$'s into \textbf{(i)} \textit{dynamical inputs}, \textbf{(ii)} \textit{participation functions}, i.e., how much area of ice with a given thickness participates to mechanical deformation \textbf{(iii)} \textit{transfer functions}, i.e., where in thickness space the ice is transferred after deformation.
+Second, we must say something about volume conservation, and this will be done more specifically later. Following \cite{thorndike_1975}, we separate the $\Psi^X$'s into \textbf{(i)} \textit{dynamical inputs}, \textbf{(ii)} \textit{participation functions}, i.e., how much area of ice with a given thickness participates to mechanical deformation \textbf{(iii)} \textit{transfer functions}, i.e., where in thickness space the ice is transferred after deformation.
 
 \section{Dynamical inputs}
 
-A general expression of $\Psi^g$, the mechanical redistribution function associated to the ice concentration, was proposed by \cite{Thorndikeetal75}:
+A general expression of $\Psi^g$, the mechanical redistribution function associated to the ice concentration, was proposed by \cite{thorndike_1975}:
 \begin{equation}
 \Psi^g = \vert \dot{\epsilon} \vert [ \alpha_o(\theta) \delta(h) + \alpha_d(\theta) w_d(h,g) ],
@@ -45 +45 @@
 \item $\vert \dot{\epsilon} \vert \alpha_d$, the net closing rate.
 \end{itemize}
-Following \cite{Thorndikeetal75}, we choose $\int_0^\infty w_d(h,g) = - 1$. In order to satisfy area conservation, the relation $\vert \dot{\epsilon} \vert \alpha_o - \vert \dot{\epsilon} \vert \alpha_d = \nabla \cdot \mathbf{u}$ must be verified. In the model, there are two ways to compute the divergence of the velocity field. A first way is to use the velocity components ($\dot{\epsilon}_I = \nabla \cdot \mathbf{u}\vert^{rhg}$) as computed after the rheology (superscript $rhg$). Another way is to derive it from the horizontal transport of ice concentration and open water fraction. In principle, the equality $A^o + \sum_{l=1}^L g^i_L = 1$ should always be verified. However, after ice transport (superscript $trp$), this is not the case, and one can diagnose a velocity divergence using the departure from this equality: $\nabla \cdot \mathbf{u}\vert^{trp} = ( 1 - A^o - \sum_{l=1}^L g^i_L )/ \Delta t$. In general, we will use $\dot{\epsilon}_I$ unless otherwise stated.
+Following \cite{thorndike_1975}, we choose $\int_0^\infty w_d(h,g) = - 1$. In order to satisfy area conservation, the relation $\vert \dot{\epsilon} \vert \alpha_o - \vert \dot{\epsilon} \vert \alpha_d = \nabla \cdot \mathbf{u}$ must be verified. In the model, there are two ways to compute the divergence of the velocity field. A first way is to use the velocity components ($\dot{\epsilon}_I = \nabla \cdot \mathbf{u}\vert^{rhg}$) as computed after the rheology (superscript $rhg$). Another way is to derive it from the horizontal transport of ice concentration and open water fraction. In principle, the equality $A^o + \sum_{l=1}^L g^i_L = 1$ should always be verified. However, after ice transport (superscript $trp$), this is not the case, and one can diagnose a velocity divergence using the departure from this equality: $\nabla \cdot \mathbf{u}\vert^{trp} = ( 1 - A^o - \sum_{l=1}^L g^i_L )/ \Delta t$. In general, we will use $\dot{\epsilon}_I$ unless otherwise stated.
 
-The \textbf{net closing rate} is written as a sum of two terms representing the energy dissipation by shear and convergence \citep{FlatoHibler95}:
+The \textbf{net closing rate} is written as a sum of two terms representing the energy dissipation by shear and convergence \citep{flato_1995}:
 \begin{equation}
 \vert \dot{\epsilon} \vert \alpha_d (\theta) = C_s \frac{1}{2} ( \Delta - \vert \dot{\epsilon}_I \vert) - \textrm{min} (\dot{\epsilon}_I,0),
@@ -68 +68 @@
 \end{equation}
 
-\textbf{Rafting} is the piling of two ice sheets on top of each other. Rafting doubles the participating ice thickness and is a volume-conserving process. \cite{Babkoetal02} concluded that rafting plays a significant role during initial ice growth in fall, therefore we included it into the model. 
+\textbf{Rafting} is the piling of two ice sheets on top of each other. Rafting doubles the participating ice thickness and is a volume-conserving process. \cite{babko_2002} concluded that rafting plays a significant role during initial ice growth in fall, therefore we included it into the model. 
 
-\textbf{Ridging} is the piling of a series of broken ice blocks into pressure ridges. Ridging redistributes participating ice on a various range of thicknesses. Ridging does not conserve ice volume, as pressure ridges are porous. Therefore, the volume of ridged ice is larger than the volume of new ice being ridged. In the model, newly ridged is has a prescribed porosity $p=30\%$ (\textit{ridge\_por} in \textit{namelist\_ice}), following observations \citep{Lepparantaetal95,Hoyland02}. The importance of ridging is now since the early works of \citep{Thorndikeetal75}.
+\textbf{Ridging} is the piling of a series of broken ice blocks into pressure ridges. Ridging redistributes participating ice on a various range of thicknesses. Ridging does not conserve ice volume, as pressure ridges are porous. Therefore, the volume of ridged ice is larger than the volume of new ice being ridged. In the model, newly ridged is has a prescribed porosity $p=30\%$ (\textit{ridge\_por} in \textit{namelist\_ice}), following observations \citep{lepp_ranta_1995,h_yland_2002}. The importance of ridging is now since the early works of \citep{thorndike_1975}.
 
 The deformation modes are formulated using \textbf{participation} and \textbf{transfer} functions with specific contributions from ridging and rafting:
@@ -90 +90 @@
 \label{eq:bri}
 \end{equation}
-where $b(h)$ is an exponential weighting function with an e-folding scale $a*$ \citep{Lipscombetal07} (\textit{astar} in \textit{namelist\_ice}) which preferentially apportions the thinnest available ice to ice deformation:
+where $b(h)$ is an exponential weighting function with an e-folding scale $a*$ \citep{lipscomb_2007} (\textit{astar} in \textit{namelist\_ice}) which preferentially apportions the thinnest available ice to ice deformation:
 \begin{equation}
 b(h) = \frac{\text{exp} [ - G(h) / a^\star ]}{a^\star[1-\text{exp}(-1/a^\star)]} ,
 \end{equation}
-It is numerically more stable than the original version of \cite{Thorndikeetal75}. This scheme is still present in the code and can be activated using \textit{partfun\_swi} from \textit{namelist\_ice}, with the associated parameter \textit{Gstar}.
+It is numerically more stable than the original version of \cite{thorndike_1975}. This scheme is still present in the code and can be activated using \textit{partfun\_swi} from \textit{namelist\_ice}, with the associated parameter \textit{Gstar}.
 
-$\beta(h)$ partitions deformation ice between rafted and ridged ice. $\beta(h)$ is formulated following \cite{Haapala00}, using the \cite{Parmerter75} law, which states that, under a critical participating ice thickness $h_P$, ice rafts, otherwise it ridges:
+$\beta(h)$ partitions deformation ice between rafted and ridged ice. $\beta(h)$ is formulated following \cite{haapala_2000}, using the \cite{parmerter_1975} law, which states that, under a critical participating ice thickness $h_P$, ice rafts, otherwise it ridges:
 \begin{equation}
 \beta(h)= \frac{ \textrm{tanh} [ -C_{ra} (h-h_P) ] + 1 } {2},
 \end{equation}
-where $C_{ra}=5$ m$^{-1}$ (\textit{Craft} in \textit{namelist\_ice})  and $h_P=0.75$ m (\textit{hparmeter} in \textit{namelist\_ice}) \citep{Haapala00,Babkoetal02}. The $tanh$ function is used to smooth the transition between ridging and rafting.  If namelist parameter \textit{raftswi} is set to 0, ice only ridges and does not raft. 
+where $C_{ra}=5$ m$^{-1}$ (\textit{Craft} in \textit{namelist\_ice})  and $h_P=0.75$ m (\textit{hparmeter} in \textit{namelist\_ice}) \citep{haapala_2000,babko_2002}. The $tanh$ function is used to smooth the transition between ridging and rafting.  If namelist parameter \textit{raftswi} is set to 0, ice only ridges and does not raft. 
 
 \section{Transfer functions}
@@ -115 +115 @@
 \label{eq:nri}
 \end{equation}
-The redistributor $\gamma(h',h)$ specifies how area of thickness $h'$ is redistributed on area of thickness $h$. We follow \citep{Hibler80} who constructed a rule, based on observations, that forces all ice participating in ridging with thickness $h'$ to be linearly distributed between ice that is between $2h'$ and $2\sqrt{H^*h'}$ thick, where $H^\star=100$ m (\textit{Hstar} in \textit{namelist\_ice}). This in turn determines how to construct the ice volume redistribution function $\Psi^v$. Volumes equal to participating area times thickness are removed from thin ice. They are redistributed following Hibler's rule. The factor $(1+p)$ accounts for initial ridge porosity $p$ (\textit{ridge\_por} in \textit{namelist\_ice}, defined as the fractional volume of seawater initially included into ridges. In many previous models, the initial ridge porosity has been assumed to be 0, which is not the case in reality since newly formed ridges are porous, as indicated by in-situ observations \citep{Lepparantaetal95,Hoyland02}. In other words, LIM3 creates a higher volume of ridged ice with the same participating ice.
+The redistributor $\gamma(h',h)$ specifies how area of thickness $h'$ is redistributed on area of thickness $h$. We follow \citep{hibler_1980} who constructed a rule, based on observations, that forces all ice participating in ridging with thickness $h'$ to be linearly distributed between ice that is between $2h'$ and $2\sqrt{H^*h'}$ thick, where $H^\star=100$ m (\textit{Hstar} in \textit{namelist\_ice}). This in turn determines how to construct the ice volume redistribution function $\Psi^v$. Volumes equal to participating area times thickness are removed from thin ice. They are redistributed following Hibler's rule. The factor $(1+p)$ accounts for initial ridge porosity $p$ (\textit{ridge\_por} in \textit{namelist\_ice}, defined as the fractional volume of seawater initially included into ridges. In many previous models, the initial ridge porosity has been assumed to be 0, which is not the case in reality since newly formed ridges are porous, as indicated by in-situ observations \citep{lepp_ranta_1995,h_yland_2002}. In other words, LIM3 creates a higher volume of ridged ice with the same participating ice.
 
 For the numerical computation of the integrals, we have to compute several temporary values:
@@ -152 +152 @@
 \section{Mechanical redistribution for other global ice variables}
 
-The other global ice state variables redistribution functions $\Psi^X$ are computed based on $\Psi^g$ for the ice age content and on $\Psi^{v^i}$ for the remainder (ice enthalpy and salt content, snow volume and enthalpy). The general principles behind this derivation are described in Appendix A of \cite{Bitzetal01}. A fraction $f_s=0.5$ (\textit{fsnowrdg} and \textit{fsnowrft} in \textit{namelist\_ice}) of the snow volume and enthalpy is assumed to be lost during ridging and rafting and transferred to the ocean. The contribution of the seawater trapped into the porous ridges is included in the computation of the redistribution of ice enthalpy and salt content (i.e., $\Psi^{e^i}$ and $\Psi^{M^s}$). During this computation, seawater is supposed to be in thermal equilibrium with the surrounding ice blocks. Ridged ice desalination induces an implicit decrease in internal brine volume, and heat supply to the ocean, which accounts for ridge consolidation as described by \cite{Hoyland02}. The inclusion of seawater in ridges does not imply any net change in ocean salinity. The energy used to cool down the seawater trapped in porous ridges until the seawater freezing point is rejected into the ocean.
+The other global ice state variables redistribution functions $\Psi^X$ are computed based on $\Psi^g$ for the ice age content and on $\Psi^{v^i}$ for the remainder (ice enthalpy and salt content, snow volume and enthalpy). The general principles behind this derivation are described in Appendix A of \cite{bitz_2001}. A fraction $f_s=0.5$ (\textit{fsnowrdg} and \textit{fsnowrft} in \textit{namelist\_ice}) of the snow volume and enthalpy is assumed to be lost during ridging and rafting and transferred to the ocean. The contribution of the seawater trapped into the porous ridges is included in the computation of the redistribution of ice enthalpy and salt content (i.e., $\Psi^{e^i}$ and $\Psi^{M^s}$). During this computation, seawater is supposed to be in thermal equilibrium with the surrounding ice blocks. Ridged ice desalination induces an implicit decrease in internal brine volume, and heat supply to the ocean, which accounts for ridge consolidation as described by \cite{h_yland_2002}. The inclusion of seawater in ridges does not imply any net change in ocean salinity. The energy used to cool down the seawater trapped in porous ridges until the seawater freezing point is rejected into the ocean.
 
 \end{document}

NEMO/trunk/doc/latex/SI3/subfiles/chap_thermo.tex

@@ -55 +55 @@
 \end{align}
 \end{linenomath}
-Hence there is no consideration of the entrainment of heat at the base of the first ocean level, or of solar radiation transmitted below the ice. The ocean-to-ice turbulent sensible heat flux is formulated following \citep{McPhee92}
+Hence there is no consideration of the entrainment of heat at the base of the first ocean level, or of solar radiation transmitted below the ice. The ocean-to-ice turbulent sensible heat flux is formulated following \citep{mcphee_1992}
 \begin{linenomath}
 \begin{align}
@@ -64 +64 @@
 
 % Ice growth
-If the $B^{opw}$ is such that the SST would decrease below the freezing point, the remainder of the heat is used to form new ice. The heat loss is converted into mass through [\ref{eq_phasechange}], giving by multiplication by density a volume of new ice $v_0$. The thickness $h_0$ of the new ice grown during a sea ice time step depends on unresolved small currents and waves and is prescribed. The fraction $a_0=v_0/h_0$ is computed accordingly. The salinity of this new ice $S_0$ is given by the S-h empirical relationship of \cite{Kovacs96}. The temperature assumed for this new ice is the local freezing point.
+If the $B^{opw}$ is such that the SST would decrease below the freezing point, the remainder of the heat is used to form new ice. The heat loss is converted into mass through [\ref{eq_phasechange}], giving by multiplication by density a volume of new ice $v_0$. The thickness $h_0$ of the new ice grown during a sea ice time step depends on unresolved small currents and waves and is prescribed. The fraction $a_0=v_0/h_0$ is computed accordingly. The salinity of this new ice $S_0$ is given by the S-h empirical relationship of \cite{kovacs_1996}. The temperature assumed for this new ice is the local freezing point.
 
 % Ice melt
-If there is ice in the grid cell and that the $B^{opw}$ is positive, it is directly given attributed to the heat available for bottom melting. This argument follows from \cite{MaykutMcPhee95}, who found that most of solar heat absorbed in the surface waters is converted into melting. In practise, this means that the SST cannot go above freezing as long ice is present.
+If there is ice in the grid cell and that the $B^{opw}$ is positive, it is directly given attributed to the heat available for bottom melting. This argument follows from \cite{maykut_1995}, who found that most of solar heat absorbed in the surface waters is converted into melting. In practise, this means that the SST cannot go above freezing as long ice is present.
 
 The heat loss used for ice formation, heat gain used to melt ice and the sensible heat given by the ocean to the ice, are all removed from the non-solar heat flux transmitted to the ocean.

NEMO/trunk/doc/latex/SI3/subfiles/chap_transport.tex

@@ -19 +19 @@
 \section{Second order moments conserving (Prather 1986) scheme (\textit{ln\_adv\_Pra})}
 
-The scheme of  \cite{Prather86} explicitly computes the conservation of second-order moments of the spatial distribution of global sea ice state variables. This scheme preserves positivity of the transported variables and is practically non-diffusive. It is also computationally expensive, however it allows to localize the ice edge quite accurately. As the scheme is conditionally stable, the time step is split into two parts if the ice drift is too fast, based on the CFL criterion. 
+The scheme of  \cite{prather_1986} explicitly computes the conservation of second-order moments of the spatial distribution of global sea ice state variables. This scheme preserves positivity of the transported variables and is practically non-diffusive. It is also computationally expensive, however it allows to localize the ice edge quite accurately. As the scheme is conditionally stable, the time step is split into two parts if the ice drift is too fast, based on the CFL criterion. 
 
 State variables per unit grid cell area are first multiplied by grid cell area. Then, for each state variable, the 0$^{th}$ (mean), 1$^{st}$ (x, y) and 2$^{nd}$ (xx, xy, yy) order moments of the spatial distribution are transported. At 1st time step, all moments are zero (if prescribed initial state); or read from a restart file, and then evolve through the course of the run. Therefore, for each global variable, 5 additional tracers have to be kept into memory and written in the restart file, which significantly increases the required memory. Advection following x and y are computed independently. The succession order of x- and y- advection is reversed every day.

NEMO/trunk/doc/latex/SI3/subfiles/introduction.tex

@@ -8 +8 @@
 \textcolor{red}{[ \textit{July 2018} ]} \\
 
-%Near the poles of the Earth, the seas and oceans freeze when seawater at the freezing point loses heat. The resulting forms of saline ice are collectively called \textit{sea ice} \citep{WMO70}, reaching up to a few meters in thickness, where as sea ice coverage is about 5\% of the global ocean, about 30 million square kilometers. All sea ice characteristics vary over a wide range of spatio-temporal scales, reflecting changes in heat, mass and momentum exchanges with the atmosphere and the ocean; the clearest temporal signal being an ample seasonal cycle. Sea ice formation and melting affects water mass formation in the ocean \citep{GoosseFichefet99}. It not only impacts, but also reflects the state of the climate system \citep{Budyko69,NotzStroeve16}.  Sea ice also affects marine life, water chemistry and human activities in polar regions. Local populations use sea ice for travelling and hunting, whereas navigation and resource exploitation are dependent on sea ice conditions. For such reasons, ocean modelling systems, including NEMO, must include a sea ice component.
+%Near the poles of the Earth, the seas and oceans freeze when seawater at the freezing point loses heat. The resulting forms of saline ice are collectively called \textit{sea ice} \citep{WMO70}, reaching up to a few meters in thickness, where as sea ice coverage is about 5\% of the global ocean, about 30 million square kilometers. All sea ice characteristics vary over a wide range of spatio-temporal scales, reflecting changes in heat, mass and momentum exchanges with the atmosphere and the ocean; the clearest temporal signal being an ample seasonal cycle. Sea ice formation and melting affects water mass formation in the ocean \citep{goosse_1999}. It not only impacts, but also reflects the state of the climate system \citep{budyko_1969,notz_2016}.  Sea ice also affects marine life, water chemistry and human activities in polar regions. Local populations use sea ice for travelling and hunting, whereas navigation and resource exploitation are dependent on sea ice conditions. For such reasons, ocean modelling systems, including NEMO, must include a sea ice component.
 
-The sea Ice Modelling Integrated Initiative (SI$^3$) is the sea ice engine of the Nucleus for European Modelling of the Ocean (NEMO). It is intended to be a flexible tool for studying sea ice and its interactions with the other components of the Earth System over a wide range of space and time scales. SI$^3$ is a curvilinear grid, finite-difference implementation of the classical AIDJEX\footnote{AIDJEX=\textbf{A}rctic \textbf{I}ce \textbf{D}ynamics \textbf{J}oint \textbf{EX}periment} model \citep{Coonetal74}, combining the conservation of momentum for viscous-plastic continuum, energy and salt-conserving halo-thermodynamics, an explicit representation of subgrid-scale ice thickness variations, snow and melt ponds. An option to switch back to the \textit{single-category} (or \textit{2-level}) framework of \cite{Hibler79} provides a cheap sea ice modelling solution.
+The sea Ice Modelling Integrated Initiative (SI$^3$) is the sea ice engine of the Nucleus for European Modelling of the Ocean (NEMO). It is intended to be a flexible tool for studying sea ice and its interactions with the other components of the Earth System over a wide range of space and time scales. SI$^3$ is a curvilinear grid, finite-difference implementation of the classical AIDJEX\footnote{AIDJEX=\textbf{A}rctic \textbf{I}ce \textbf{D}ynamics \textbf{J}oint \textbf{EX}periment} model \citep{coon_1974}, combining the conservation of momentum for viscous-plastic continuum, energy and salt-conserving halo-thermodynamics, an explicit representation of subgrid-scale ice thickness variations, snow and melt ponds. An option to switch back to the \textit{single-category} (or \textit{2-level}) framework of \cite{hibler_1979} provides a cheap sea ice modelling solution.
 
 SI$^3$ is the result of the recommendation of the Sea Ice Working Group (SIWG) to reduce duplication and better use development resources. SI$^3$ merges the capabilities of the 3 formerly used NEMO sea ice models (CICE, GELATO and LIM). The \textbf{3} in SI$^3$ refers to the three formerly used sea ice models. It also refers to linkages between 3 different media (ocean, ice, snow). The model can be spelt 'SI3' in situations where the superscript could be problematic (i.e., within code and svn repository etc.) The model name would be pronounced as 'si-cube' for short (or 'sea ice cubed' for slightly longer). 
 
 % Limitations & scope
-%There are limitations to the applicability of models such as SI$^3$. The continuum approach is not invalid for grid cell size above at least 1 km, below which sea ice particles may include just a few floes, which is not sufficient \citep{Lepparanta05}. Second, one must remember that our current knowledge of sea ice is not as complete as for the ocean: there are no fundamental equations such as Navier Stokes equations for sea ice. Besides, important features and processes span widely different scales, such as brine inclusions (1 $\mu$m-1 mm) \citep{PerovichGow96}, horizontal thickness variations (1 m-100 km) \citep{Percivaletal08}, deformation and fracturing (10 m-1000 km) \citep{Marsanetal04}. These impose complicated and often subjective subgrid-scale treatments. All in all, there is more empirism in sea ice models than in ocean models. 
+%There are limitations to the applicability of models such as SI$^3$. The continuum approach is not invalid for grid cell size above at least 1 km, below which sea ice particles may include just a few floes, which is not sufficient \citep{lepp_ranta_2011}. Second, one must remember that our current knowledge of sea ice is not as complete as for the ocean: there are no fundamental equations such as Navier Stokes equations for sea ice. Besides, important features and processes span widely different scales, such as brine inclusions (1 $\mu$m-1 mm) \citep{perovich_1996}, horizontal thickness variations (1 m-100 km) \citep{percival_2008}, deformation and fracturing (10 m-1000 km) \citep{marsan_2004}. These impose complicated and often subjective subgrid-scale treatments. All in all, there is more empirism in sea ice models than in ocean models. 
 
 In order to handle all the subsequent required subjective choices, we applied the following guidelines or principles: