Changeset 11031
- Timestamp:
- 2019-05-21T21:41:02+02:00 (6 years ago)
- Location:
- NEMO/trunk/doc/latex/SI3/subfiles
- Files:
-
- 1 deleted
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/doc/latex/SI3/subfiles/abstract_foreword.tex
r11015 r11031 10 10 \chapter*{Abstract} 11 11 12 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.12 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. 13 13 14 14 %\vspace{-40pt} -
NEMO/trunk/doc/latex/SI3/subfiles/chap_domain.tex
r11015 r11031 75 75 \nlst{namitd} 76 76 77 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$).77 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$). 78 78 79 79 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). 101 101 102 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.102 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. 103 103 104 104 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
r11015 r11031 28 28 29 29 \subsection{Scales, thermodynamics and dynamics} 30 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.30 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. 31 31 32 32 \subsection{Subgrid scale variations} 33 33 34 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}.34 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}. 35 35 36 36 %-------------------------------------------------------------------------------------------------------------------- … … 49 49 %-------------------------------------------------------------------------------------------------------------------- 50 50 51 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.52 53 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}.51 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. 52 53 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}. 54 54 55 55 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 71 71 $c_i$ (cpic) & Pure ice specific heat & 2067 & J/kg/K & ? \\ 72 $c_w$ (rcp) & Seawater specific heat & 3991 & J/kg/K & \cite{TEOS -10} \\73 $L$ (lfus) & Latent heat of fusion (0$^\circ$C) & 334000 & J/kg/K & \cite{ BitzLipscomb99} \\74 $\rho_i$ (rhoic) & Sea ice density & 917 & kg/m$^3$ & \cite{ BitzLipscomb99} \\75 $\rho_s$ (rhosn) & Snow density & 330 & kg/m$^3$ & \cite{ MaykutUntersteiner71} \\76 $\mu$ (tmut) & Linear liquidus coefficient & 0.054 & $^\circ$C/(g/kg) & \cite{ Assur58} \\72 $c_w$ (rcp) & Seawater specific heat & 3991 & J/kg/K & \cite{TEOS_2010} \\ 73 $L$ (lfus) & Latent heat of fusion (0$^\circ$C) & 334000 & J/kg/K & \cite{bitz_1999} \\ 74 $\rho_i$ (rhoic) & Sea ice density & 917 & kg/m$^3$ & \cite{bitz_1999} \\ 75 $\rho_s$ (rhosn) & Snow density & 330 & kg/m$^3$ & \cite{maykut_1971} \\ 76 $\mu$ (tmut) & Linear liquidus coefficient & 0.054 & $^\circ$C/(g/kg) & \cite{assur_1958} \\ 77 77 \hline 78 78 \end{tabular} … … 85 85 % 86 86 87 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}.87 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}. 88 88 89 89 … … 154 154 \subsection{Dynamic formulation} 155 155 156 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}.156 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}. 157 157 158 158 %------------------------------------------------------------------------------------------------------------------------- … … 162 162 %------------------------------------------------------------------------------------------------------------------------- 163 163 164 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:164 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: 165 165 \begin{linenomath} 166 166 \begin{align} … … 168 168 \end{align} 169 169 \end{linenomath} 170 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}:170 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}: 171 171 \begin{linenomath} 172 172 \begin{align} … … 199 199 \end{align} 200 200 \end{linenomath} 201 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:201 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: 202 202 \begin{linenomath} 203 203 \begin{align} … … 206 206 \end{align} 207 207 \end{linenomath} 208 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.208 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. 209 209 210 210 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} 233 233 \end{linenomath} 234 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:234 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: 235 235 \begin{linenomath} 236 236 \begin{align} … … 250 250 251 251 Dynamical processes include the conservation of momentum, rheology, transport and mechanical redistribution. 252 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}:252 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}: 253 253 \begin{linenomath} 254 254 \begin{align} … … 257 257 \end{align} 258 258 \end{linenomath} 259 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:259 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: 260 260 \begin{linenomath} 261 261 \begin{align} … … 294 294 % Parameters 295 295 $P*$ (rn\_pstar) & ice strength thickness param. & 20000 & N/m2 & - \\ 296 $C$ (rn\_crhg) & ice strength concentration param. & 20 & - & \citep{ Hibler79} \\297 $H^*$ (rn\_hstar) & maximum ridged ice thickness param. & 25 & m & \citep{ Lipscombetal07} \\298 $p$ (rn\_por\_rdg) & porosity of new ridges & 0.3 & - & \citep{ Lepparantaetal95} \\296 $C$ (rn\_crhg) & ice strength concentration param. & 20 & - & \citep{hibler_1979} \\ 297 $H^*$ (rn\_hstar) & maximum ridged ice thickness param. & 25 & m & \citep{lipscomb_2007} \\ 298 $p$ (rn\_por\_rdg) & porosity of new ridges & 0.3 & - & \citep{lepp_ranta_1995} \\ 299 299 $amax$ (rn\_amax) & maximum ice concentration & 0.999 & - & -\\ 300 300 $h_0$ (rn\_hnewice) & thickness of newly formed ice & 0.1 & m & - \\ … … 311 311 312 312 313 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.314 315 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.313 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. 314 315 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. 316 316 317 317 \section{Ice thermodynamics} … … 319 319 320 320 \subsection{Transport in thickness space} 321 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.321 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. 322 322 323 323 \subsection{Thermodynamic source and sink terms} … … 391 391 %-------------------------------------------------------------------------------------------------------------------- 392 392 393 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}:393 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}: 394 394 \begin{linenomath} 395 395 \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. 400 400 401 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:401 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: 402 402 \begin{linenomath} 403 403 \begin{align} … … 418 418 \item the sensible heat flux from the ocean to the sea ice ($A.F_w$) 419 419 \end{itemize} 420 Other contributions are not assumed not to contribute. The ocean-to-ice sensible heat flux is formulated the bulk formula of \citep{ McPhee92}.420 Other contributions are not assumed not to contribute. The ocean-to-ice sensible heat flux is formulated the bulk formula of \citep{mcphee_1992}. 421 421 422 422 % Ice growth 423 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.423 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. 424 424 425 425 $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} 436 436 \end{linenomath} 437 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.437 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. 438 438 439 439 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. 448 448 449 \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.449 \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. 450 450 451 451 \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. 452 452 453 \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}.454 455 \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.453 \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}. 454 455 \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. 456 456 457 457 \end{document} -
NEMO/trunk/doc/latex/SI3/subfiles/chap_output_diagnostics.tex
r11015 r11031 16 16 \section{SIMIP diagnostics} 17 17 18 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.18 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. 19 19 20 20 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
r11015 r11031 36 36 %-------------------------------------------------------------------------------------------------------------------- 37 37 38 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.38 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. 39 39 40 40 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. 51 51 52 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}.53 54 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).52 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}. 53 54 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). 55 55 56 56 \nlst{namalb} … … 64 64 \vspace{0cm} 65 65 \includegraphics[height=10cm,angle=-00]{albedo_cloud_correction} 66 \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).}66 \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).} 67 67 % ocean uses 0.06 for overcast sky (Payne 74) and Briegleb and Ramanathan parameterization 68 68 \label{fig_albedo_cloud_correction} … … 72 72 %-------------------------------------------------------------------------------------------------------------------- 73 73 74 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}:74 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}: 75 75 \begin{equation} 76 76 \alpha(i,j,l) = [ 1 - c(i,j) ] \cdot \alpha_{cs} (i,j,l) + c (i,j) \cdot \alpha_{os}(i,j,l). 77 77 \end{equation} 78 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 through78 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 79 79 \begin{equation} 80 80 \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. 111 111 112 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:112 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: 113 113 \begin{eqnarray} 114 114 \alpha_{ice} = … … 129 129 values that are to be specified from the namelist. 130 130 131 \cite{ GrenfellPerovich04} suggest that the dependence of surface albedo on snow depth is exponential,131 \cite{grenfell_2004} suggest that the dependence of surface albedo on snow depth is exponential, 132 132 \begin{eqnarray} 133 133 \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. 144 144 145 Based on ideas developed from melt ponds on continental ice \citep{ ZuoOerlemans96}, the albedo of ponded ice was proposed to follow \citep{Lecomteetal11}:145 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}: 146 146 \begin{eqnarray} 147 147 \alpha_{pnd} = \alpha_{dpnd} - ( \alpha_{dpnd} -\alpha_{ice} ) \cdot exp( -h_{pnd} / 0.05 ) 148 148 \end{eqnarray} 149 $\alpha_{dpnd}$ is a namelist parameter. \cite{ EbertCurry93} also use such dependency for their multi-spectral albedo.149 $\alpha_{dpnd}$ is a namelist parameter. \cite{ebert_1993} also use such dependency for their multi-spectral albedo. 150 150 151 151 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} 154 154 155 The transmitted solar radiation below the surface is represented following \cite{ FichefetMaqueda97} and \cite{MaykutUntersteiner71}:155 The transmitted solar radiation below the surface is represented following \cite{fichefet_1997} and \cite{maykut_1971}: 156 156 \begin{eqnarray} 157 157 qtr\_ice\_top(i,j,l) = i_o(i,j) qsr\_ice(i,j,l), 158 158 \end{eqnarray} 159 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.159 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. 160 160 161 161 \subsection{Attenuation and transmission below the ice/ocean interface} -
NEMO/trunk/doc/latex/SI3/subfiles/chap_ridging_rafting.tex
r11015 r11031 29 29 \end{equation} 30 30 31 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.31 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. 32 32 33 33 \section{Dynamical inputs} 34 34 35 A general expression of $\Psi^g$, the mechanical redistribution function associated to the ice concentration, was proposed by \cite{ Thorndikeetal75}:35 A general expression of $\Psi^g$, the mechanical redistribution function associated to the ice concentration, was proposed by \cite{thorndike_1975}: 36 36 \begin{equation} 37 37 \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. 46 46 \end{itemize} 47 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.47 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. 48 48 49 The \textbf{net closing rate} is written as a sum of two terms representing the energy dissipation by shear and convergence \citep{ FlatoHibler95}:49 The \textbf{net closing rate} is written as a sum of two terms representing the energy dissipation by shear and convergence \citep{flato_1995}: 50 50 \begin{equation} 51 51 \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} 69 69 70 \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.70 \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. 71 71 72 \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}.72 \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}. 73 73 74 74 The deformation modes are formulated using \textbf{participation} and \textbf{transfer} functions with specific contributions from ridging and rafting: … … 90 90 \label{eq:bri} 91 91 \end{equation} 92 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:92 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: 93 93 \begin{equation} 94 94 b(h) = \frac{\text{exp} [ - G(h) / a^\star ]}{a^\star[1-\text{exp}(-1/a^\star)]} , 95 95 \end{equation} 96 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}.96 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}. 97 97 98 $\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:98 $\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: 99 99 \begin{equation} 100 100 \beta(h)= \frac{ \textrm{tanh} [ -C_{ra} (h-h_P) ] + 1 } {2}, 101 101 \end{equation} 102 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.102 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. 103 103 104 104 \section{Transfer functions} … … 115 115 \label{eq:nri} 116 116 \end{equation} 117 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.117 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. 118 118 119 119 For the numerical computation of the integrals, we have to compute several temporary values: … … 152 152 \section{Mechanical redistribution for other global ice variables} 153 153 154 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.154 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. 155 155 156 156 \end{document} -
NEMO/trunk/doc/latex/SI3/subfiles/chap_thermo.tex
r11015 r11031 55 55 \end{align} 56 56 \end{linenomath} 57 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}57 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} 58 58 \begin{linenomath} 59 59 \begin{align} … … 64 64 65 65 % Ice growth 66 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.66 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. 67 67 68 68 % Ice melt 69 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.69 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. 70 70 71 71 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
r11015 r11031 19 19 \section{Second order moments conserving (Prather 1986) scheme (\textit{ln\_adv\_Pra})} 20 20 21 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.21 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. 22 22 23 23 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
r9974 r11031 8 8 \textcolor{red}{[ \textit{July 2018} ]} \\ 9 9 10 %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.10 %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. 11 11 12 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.12 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. 13 13 14 14 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). 15 15 16 16 % Limitations & scope 17 %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.17 %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. 18 18 19 19 In order to handle all the subsequent required subjective choices, we applied the following guidelines or principles:
Note: See TracChangeset
for help on using the changeset viewer.