%% Use the "review" option when submitting for review.
%% Removing the "review" option will switch the
%% manuscript to a two-column layout.
%%
%% Use option "jog" for submissions to the Journal of
%% Glaciology; use "aog" for submission to the Annals
%% of Glaciology.
% \documentclass[jog]{igs}
\documentclass[review,jog]{igs}

\usepackage[utf8]{inputenc}

\usepackage{mathabx}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{blindtext}

% for variance et expected value
\def\Var{{\textrm{Var}}\,}
\def\E{{\textrm{E}}\,}

% for having Fig. S1 instead of Fig. 1
\renewcommand{\thefigure}{S\arabic{figure}} % Redéfinition globale

% the default is for unnumbered section heads
% if you really must have numbered sections, remove
% the % from the beginning of the following command
% and insert the level of sections you wish to be
% numbered (up to 4):

\setcounter{secnumdepth}{2}

\jourvolume{12}
\jourissue{3}
\jourpubyear{2016}

\begin{document}

% \title[Validating an ISMIP model with observations]{Validating an ISMIP model with observations}

\title[Validating ensemble historical simulations]{Validating ensemble historical simulations of Upernavik
Isstrøm (1985-2019) with different friction parametrisation using observations of surface velocity and elevation}

\author[Jager and others]{Eliot JAGER$^1$,
Fabien GILLET-CHAULET$^1$, 
Jérémie MOUGINOT$^{1,2}$,
Romain MILLAN$^{1,3}$}

\affiliation{%
$^1$IGE, Univ. Grenoble Alpes, CNRS, IRD, 38000 Grenoble, France\\
$^2$Department of Earth System Science, University of California, Irvine, USA\\
$^3$Department of Geosciences and Natural Resource Management, University 
of Copenhagen, Copenhagen K, Denmark.\\
  Correspondence: Eliot JAGER 
  \email{eliot.jager@univ-grenoble-alpes.fr}}



\section{Introduction}

This document is a supplementary to \textit{"Validating ensemble historical simulations of Upernavik Isstrøm (1985-2019) using observations of surface velocity and elevation"} and provides detailed descriptions of the model equations and the data-set.


\section{Supplementary model description}

Here we give the detailed equations used in the model and described in sec 1.1 Model description.

\label{sec:IceFlowModel}

\subsection{Model equations}

For the force balance, we solve the Shelfy-Stream Approximation where the horizontal velocity $\mathbf{u}=(u,v)$ is solution of

\begin{equation}
\begin{cases}
2 \dfrac{\uppartial T_{xx}}{\uppartial x} + \dfrac{\uppartial T_{yy}}{\uppartial x} + \dfrac{\uppartial T_{xy}}{\uppartial y} - \tau_{bx} = \rho_i g H \dfrac{\uppartial z_s}{\uppartial x} \\
2 \dfrac{\uppartial T_{yy}}{\uppartial y} + \dfrac{\uppartial T_{xx}}{\uppartial y} + \dfrac{\uppartial T_{xy}}{\uppartial x} - \tau_{by} = \rho_i g H \dfrac{\uppartial z_s}{\uppartial y}
\end{cases}
\label{ssa}
\end{equation}
with
\begin{itemize}[label=\textbullet]
\item $\mathbf{T} $  the membrane stresses defined as  $T_{xx} = 2 H \bar{\eta} \uppartial u/ \uppartial x$, $T_{yy} = 2 H \bar{\eta} \uppartial v/ \uppartial y$, $T_{xy} = H \bar{\eta} (\uppartial u/ \uppartial y + \uppartial v/ \uppartial x)$
\item $\mathbf{\tau_b} = (\tau_{bx},\tau_{by})$ the basal shear stress,
\item $H=z_s-z_b$ is the ice thickness, with $z_s$ and $z_b$ the top and bottom surface elevations, respectively,
\item $\rho_i$ the ice density,
\item  $g$ the gravity.
\end{itemize}

The non-linear vertically-averaged effective viscosity $\bar{\eta}$ is given by 
\begin{equation}
\bar{\eta} =  \dfrac{1}{H} \int_{z_b}^{z_s} \dfrac{1}{2} (EA)^{-1/n} d_e^{-(1-n)/n} \text{d}z
\label{eta}
\end{equation}
with $A$ the rate factor , $E$ the enhancement factor, $n$ the Glen's exponent and $d_e$ the effective strain-rate given by
\begin{equation}
d_e^2 = \left( \dfrac{\uppartial u}{\uppartial x}\right)^2 + \left( \dfrac{\uppartial v}{\uppartial y}\right)^2 + \dfrac{\uppartial u}{\uppartial x}\dfrac{\uppartial v}{\uppartial y} + \dfrac{1}{4} \left( \dfrac{\uppartial u}{\uppartial y}+\dfrac{\uppartial v}{\uppartial x}\right)^2
 \label{d_e}
\end{equation}


The lateral sides of the model domain correspond to flow lines and we set a no-flow condition in the normal direction and a free-slip condition in the tangential direction. A Neumann condition that takes into account the sea water pressure is applied at the front boundary. 

The vertically integrated mass conservation equation for the evolution of the ice thickness is given by : 
\begin{equation}
 \dfrac{\uppartial H}{\uppartial t} + \mathbf{\nabla} \cdot (H\mathbf{u}) = \dot{a}_s - \dot{a}_b
\label{thickness}
\end{equation}
with $\dot{a}_s$ the surface mass balance and $\dot{a}_b$ the basal melt rate. The two free surfaces $z_s$ and $z_b$ are obtained from the floating condition, which for a constant sea level $z_{sl}=0$ m gives :
\begin{equation}
\label{Floating}
\begin{cases}
z_b=b & \text{for } H \geq -b \dfrac{\rho_w}{\rho_i}  \\
z_b=-H\dfrac{\rho_i}{\rho_w} & \text{otherwise} \\
\end{cases}
\end{equation}

\begin{equation}
\label{FloatingZs}
z_s = z_b + H
\end{equation}
with $b$ the bed elevation and $\rho_w$ the sea water density.


\subsection{Initial friction field inversion}

\label{sec:slcInv}

The initial friction coefficient field $\upbeta_W$ (section 2.1 Model parameters) is found using a classical inverse method by minimising the following cost function : 

\begin{equation}
\label{Jtot}
J(\upalpha) = J_{0}(\upalpha) + \lambda_{reg} J_{reg}(\upalpha) + \lambda_{div} J_{div}(\upalpha) 
\end{equation}{} 

where $\upalpha = \log_{10} (\upbeta_W)$, $\lambda_{reg}$ and $\lambda_{div}$ are different weights to the regularisation, $J_{reg}$, and divergence, $J_{div}$, cost functions. 


$J_0$ measures the misfit between an observed velocity field $\mathbf{u_{obs}}$ and the model velocity field $\mathbf{u^{mod}}$ at the $N_{obs}$ observation points: 

\begin{equation}
J_0 = \sum_{i=1}^{N_{obs}} \dfrac{1}{2}(||\mathbf{u_{mod}^i} - \mathbf{u_{obs}^i}||)^2 
\label{J0}
\end{equation}{} 

$J_{reg}$ is a Tikhonov regularisation term that imposes a smoothness constraint on $\upalpha$, to avoid over-fitting : 

\begin{equation}
J_{reg} = \underset{\Sigma}{\int} \dfrac{1}{2} ||\mathbf{\nabla} \upalpha||^2 d\Sigma
    \label{Jreg}
\end{equation}{}

with $\Sigma$ the surface contact between ice and bed.

According to Eq. \ref{thickness}, $J_{div}$ measures the misfit between the model and observed thickness rate of change $(\text{d} H/\text{d} t)_{obs}$ :

\begin{equation}
J_{div} = \underset{\Sigma}{\int} \dfrac{1}{2} \left(\mathbf{\nabla} \cdot (H\mathbf{u}) - \dot{a} - \left(\dfrac{\text{d} H}{\text{d} t}\right)_{obs} \right)^2 d\Sigma
    \label{Jdiv}
\end{equation}{}

For $(\text{d} H/\text{d} t)_{obs}$, we use an average surface elevation change with data collected between 1992 and 2018 \cite{Gourmelen2018}.


\subsection{parameterisation of $\upbeta_{RC}$}

\label{sec:b_f_d}

Here we give more details on the parameterisation Eq. 3 used for the friction coefficient $\upbeta_{RC}$ in the regularised-Coulomb friction law (Eq. 2). This parametrisation is meant to represent the dependence of $\upbeta_{RC}$ on the effective pressure $N$ in the vicinity of the ice front.

It is based on the observation that the friction coefficients determined during the inversion under the central flow lines of UI-N, UI-C and UI-S vary linearly with the distance to the front in the first few kilometers as shown for one of the initial state in Fig. \ref{fig:Beta_RC_f_d}. This dependence is found only under the central flow lines where the coefficients are the smallest. This result can be linked to the fact that the height above flotation $h_{af}$ also varies nearly linearly with the distance to the front along this flow lines. In \cite{tc-7-1679-2013}, they show that a linear relationship between $h_{af}$ and $\upbeta_{RC}$ ($\tau_c$ in their paper) can be obtained for the fastest parts of Jakobshavn. They also find that the time variations of $\upbeta_{RC}$ and $h_{af}$ are similar, although the variations of $\upbeta$ are more localised in the stream. However, introducing a direct dependence of $\upbeta_{RC}$ to the height above flotation would also affect the friction coefficients on the margins where they are higher, while our parameterisation introduces a dependence mostly under the fast flowing parts where the friction is the smallest. As indicated in the paper, making $\upbeta_{RC}$ depend on $h_{af}$ leads to results that are really sensitive to the surface elevation, which is poorly known in 1985. So we use distance to the front as a proxy of $h_{af}$.

\begin{figure}[h!]
    \centering
    \includegraphics[width=86mm]{Appendix/beta_fit_paper_end.png}
    \caption{Friction coefficient $\upbeta_{RC}$ for one of RCE member in function of the distance to the front during the inversion $d_{inv}$ for the points situated under the central flow-lines of UI-N (cyan), UI-C (blue) and UI-S (purple). The red dashed lines represent the two asymptotical behaviours, the vertical black dotted line corresponds to $d_{inv}=d_{lim}$ and in black/grey solid line $\upbeta_{fit}$ for two different front position in 2015/1985.}
    \label{fig:Beta_RC_f_d}
\end{figure}


$\upbeta_{RC}$ in Eq. (3) is made of two terms :

\begin{itemize}[label=\textbullet]
    \item $\upbeta_{ref}$ which represent the dependence of $\upbeta_{RC}$ on time invariant parameters
    \item $\upbeta_{fit}(d) = \upbeta_{lim} \dfrac{d}{d+d_{lim}}$ which introduces a dependence to the distance to the front d (Fig. \ref{fig:Beta_RC_f_d}).
\end{itemize}

$\upbeta_{fit}$ exhibits two asymptotic behaviours (Fig. \ref{fig:Beta_RC_f_d}) :

\begin{itemize}[label=\textbullet]
    \item when $d \ll d_{lim}$, $\upbeta_{fit}(d) \sim \frac{\upbeta_{lim}}{d_{lim}} \cdot d$ and changes linearly with $d$.
    \item when $d \gg d_{lim}$, $\upbeta_{fit}(d) \sim \upbeta_{lim}$ and is independent of d.
\end{itemize}

First, $\upbeta_{lim}$ is determined as the average of $\upbeta_{RC,inv}$, obtain from the inversion, over the validation area. Then, the value of $\frac{\upbeta_{lim}}{d_{lim}}$ is determined using a linear regression in the form $\beta_{RC} = \frac{\upbeta_{lim}}{d_{lim}} \cdot d$ (figure \ref{fig:Beta_RC_f_d}) for the points situated under the central flow-lines within the first 30 km from the front (see figure \ref{fig:finalSlc}).



Then, $\upbeta_{ref}$, which must be time invariant, is computed with the following relation:  

\begin{equation}
    \upbeta_{ref} = \upbeta_{RC} - \upbeta_{fit}(d_{inv}) 
\label{beta_rc_ref}
\end{equation}{}

Our parametrisation reproduce the behaviour observed in \cite{tc-7-1679-2013} (Fig. 13 in their paper, Fig. \ref{fig:finalSlc} here) with strong changes of the friction coefficient in the streams near the front, up to -60\% between the beginning and the end of the simulation, but no major changes outside (0\%).


\begin{figure*}[ht!]
    \centering
    \includegraphics[width=140mm]{Appendix/slc_change_fin.png}
    \caption{Changes of friction coefficient between the initial and final states for one member of RCE; In white, central flow-lines of the three main glaciers UI-N, UI-C and UI-S used to obtain $\frac{\upbeta_{lim}}{d_{lim}}$ in $\upbeta_{fit}$.}
    \label{fig:finalSlc}
\end{figure*}


\clearpage

\section{Supplementary data description}

\subsection{Averaging method}

As mentioned in the Validation Data section, we compile the collected surface velocity and surface elevation data into an annual time series. Here we provide plots of our data-set (Figures \ref{fig:vobs} and \ref{fig:hobs}) and further details on the spurious value filtering.

The individual velocity maps are filtered out for remaining spurious values when the relative difference compared to the median of all values is larger than 150\% for velocity larger than 1.0 kma$^{-1}$ and 100\% for velocity between 0.2 and 1.0 kma$^{-1}$, and when the absolute difference is larger than 50 ma$^{-1}$ for velocity lower than 200 ma$^{-1}$.

For the surface elevation and ASTER observations between 2000 and 2020, we used the MMASTER processing chain which is based on the MICMAC photogrammetric software \cite{rs9070704}. Photogrammetry based DEMs were coregistered to yearly reference elevation datasets. The reference datasets were assembled using all altimetry data on stable ground, and elevation measurements from the same year as the DEMs on ice \citep{tc-13-2633-2019}. We discarded the outliers pixels that have a difference with GIMP that is greater than 200-m. For elevations above 1000-m (i.e in the accumulation area), we also filtered pixels with a difference larger than 75 meters with the GIMP DEM \citep{Howat_2022}, as dh/dt are expected to be much lower in this region. The threshold values were determined arbitrarily after numerous iterative filtering tests specific to the Upernavik region. 

\clearpage
\begin{figure*}[h!]
    \centering
    \includegraphics[width=150mm]{Appendix/vobs.png}
    \caption{Time-serie of ice surface velocity collected for this work}
    \label{fig:vobs}
\end{figure*}

\begin{figure*}[ht!]
    \centering
    \includegraphics[width=150mm]{Appendix/hobs.png}
    \caption{Time-serie of ice surface elevation collected for this work}
    \label{fig:hobs}
\end{figure*}

\clearpage
\subsection{Uncertainty estimation}

The estimation of the uncertainties of surface velocities and surface elevations in figures 5 and 6 of the paper is done as follows.

For the surface flow velocity the error is determined using a three step procedure. First, we selected the year 2019, which has a high temporal resolution in 2019, and calculated the standard deviation of these seasonal velocities at each grid cell. Secondly, we evaluate an error function using a polynomial regression fitting the standard deviation relative to the magnitude of the surface flow velocity (Figure \ref{fig:stdv_f_v}). This error function can then be applied to any surface flow velocity field, and is weighted by the number of observations available within that pixel.

For the annual mean surface elevation, we assume an error of 0.1 m for satellite and airborne altimeters (ICESat, ICESat-2, ATM, LVIS). For the GIMP and GLISTIN-A DEMs we used the standard deviation grid provided along with the dataset \citep{GLISTIN_DataCitation,GIMP_DataCitation}. Error on the surface elevation from ASTER DEMs is set to 7m, and was calculated as the mean difference with altimetry data. We use a conservative error of 10 m on the 1985 historical DEM \citep{Korsgaard2016}.

\begin{figure}[h!]
    \centering
    \includegraphics[width=86mm]{Appendix/final_vobs_std.png}
    \caption{Evolution of the standard deviation of the velocity in 2019 as a function of its annual average. Each point represents a pixel and the orange line the polynomial regression.}
    \label{fig:stdv_f_v}
\end{figure}

% \clearpage

\section{Additional plots}


Additional figures are mentioned several times in the main text. These are listed below. In addition to the two ensembles WE and RCE presented in the main text, we run a third ensemble WE85. It has the same characteristics as WE with a constant friction field, but this field is corrected to reproduce the initial velocities.


To achieve this, we start from $\beta_W$ inversions using 2010s velocities, and compute $\upbeta_{W}(1985)$ by reusing Éq. 3 of the main text to obtain higher friction coefficients: 

\begin{equation}
    \upbeta_{W}(1985) = \upbeta_{W}(2010) + \upbeta_{lim} \dfrac{d(2010)}{d(2010)+d_{lim}} - \upbeta_{lim} \dfrac{d(1985)}{d(1985)+d_{lim}}
\label{paramBeta}
\end{equation}{}

with $d(1985)$ the distance to the ice front of 1985, $d(2010)$ the distance to the ice front of bedmachine v3, $\upbeta_{W}(2010)$ the friction field obtain by inverse method and $\upbeta_{W}(1985)$ the friction field used for WE85. To compute $d_{lim}$ and $\upbeta_{lim}$, we used the same method as explained in the main text. In the front retreat areas, i.e. where there is no ice in 2010 but where there was ice in 1985, we have used a $\upbeta_{W}(2010)$ value of 0. This extrapolation method is the same as that used for RCE.

In contrast to the WE friction field, this ensemble exhibits a linear dependence of friction on the distance from the 1985 front within the front retreat areas. Additionally, higher friction coefficients are observed upstream of the 2010 front for glaciers UI-N and UI-C, both of which have undergone front retreat.



This new ensemble correctly reproduces the initial velocity field (1985-2005), like RCE. We then show below its effect compared to RCE and WE on the velocities and transient elevations of points A and B, on the velocity field for the end of the period (2014-2019) and on the global variables of ice discharge and ice mass changes.


\begin{figure}[h!]
    \centering
    \includegraphics[width=86mm]{Appendix/Up_N_pointP_v_zsWE85.jpg}
    \caption{Surface velocity (top) and surface elevation (bottom) at point A (see figure \ref{fig:2DmapVelocity}) of WE (in blue), RCE (in red), WE85 (in purple) and observations (black dots with an estimate of the uncertainty in grey). For WE and RCE the mean is represented in solid line and the shading include 95\% of the ensemble members.}
    \label{fig:time_pointA}
\end{figure}

\begin{figure}[h!]
    \centering
    \includegraphics[width=86mm]{Appendix/Up_S_pointP_v_zsWE85.jpg}
    \caption{Same as figure \ref{fig:time_pointA} for point B (see figure \ref{fig:2DmapVelocity})}
    \label{fig:time_pointB}
\end{figure}


\begin{figure*}[h!]
    \centering
    \includegraphics[width=178mm]{Appendix/vvstatR15vsR18InitStateAB_WE85.jpg}
    \caption{Surface velocity bias (top), MAE (middle) of the ensemble mean and CRPS (bottom) for WE85 (left) and RCE (right) during the period 2014-2019. Point A ($+$) and Point B ($\times$) are used in Figures \ref{fig:time_pointA} and \ref{fig:time_pointB} as representative of UI-N and UI-S respectively. The grey and black lines in the first row are the 200 and 1000~m,a$^{-1}$ velocity contours computed from RCE 2014-2019 average.}
    \label{fig:2DmapVelocity}
\end{figure*}

\begin{figure*}[h!]
    \centering
    \includegraphics[width=140mm]{Appendix/IceDischargeVolumeWeVsWe85vsRCHist.jpg}
    \caption{Ice discharge (top graph) and cumulative ice mass change (bottom graph) for RCE (red), WE (blue) and WE85 (purple) between 1986 and 2019, with mean in solid line and the shading include 95\% of the ensemble members, against different observation : Mouginot (+), Mankoff ({\scriptsize Y}) and King ($\times$). On the right, histograms of ice Discharge and ice mass change in 2019}
    \label{fig:IceDischargeMean}
\end{figure*}


\clearpage

\bibliography{igsrefs}
\bibliographystyle{igs}


\end{document}

 