Inversions

Introduction

Inversions are used to constrain poorly known model parameters such as basal friction. The method consists of finding a set of model inputs that minimizes the cost function \({\mathcal J}\) that measures the misfit between model and observations. For example, inverse methods are used to infer the basal friction \(k\):

\[\boldsymbol{\tau}_b = -k^2 N^r \|{\bf v}\|^{s-1} {\bf v}_b\]

and/or the depth-averaged ice hardness, \(B\), in Glen’s flow law:

\[\mu = \frac{B}{2\left( \dot{\varepsilon}_e^{1-\frac{1}{n}}\right) }\]

This section explains how to solve inverse problems and how optimization parameters must be tuned.

Cost functions

Absolute misfit

This is the classic way of calculating a misfit between a modeled and observed velocity field:

\[{\mathcal J\left({\bf v}\right)}=\int_{S} \dfrac{1}{2}\left(\left(v_x-v_x^{\text{obs}}\right)^{2}+\left(v_y-v_y^{\text{obs}}\right)^{2}\right) dS\]

where:

  • vx is the x component of the glacier modeled velocity
  • vy is the y component of the glacier modeled velocity
  • vxobs is the x component of the glacier observed velocity
  • vyobs is the y component of the glacier observed velocity

Relative misfit

The relative misfit is defined as follows:

\[{\mathcal J\left({\bf v}\right)}=\int_{S} \dfrac{1}{2}\left(\dfrac{\left(v_x-v_x^{\text{obs}}\right)^{2}}{\left(v_x^{\text{obs}}+\varepsilon\right)^{2}}+\dfrac{\left(v_y-v_y^{\text{obs}}\right)^{2}}{\left(v_y^{\text{obs}}+\varepsilon\right)^{2}}\right) dS\]

where:

  • \(\varepsilon\) is a minimum velocity used to avoid the observed velocity being equal to zero.

Logarithmic misfit

\[{\mathcal J\left({\bf v}\right)}=\int_{S} \left(\text{log}\left(\dfrac{\|{\bf v}\|+\varepsilon}{\|{\bf v}^{\text{obs}}\|+\varepsilon}\right) \right)^2 dS\]

where:

  • v is the glacier modeled velocity magnitude
  • vobs is the glacier observed velocity magnitude
  • \(\varepsilon\) is a minimum velocity used to avoid the observed velocity being equal to zero

Thickness misfit

\[{\mathcal J\left(H\right)}=\int_{\Omega} \dfrac{1}{2}\left(H-H^{\text{obs}}\right)^{2}d\Omega\]

where:

  • H is the ice thickness
  • Hobs is the measured ice thickness

Drag gradient

\[{\mathcal J\left(k\right)}=\int_{B} \gamma \dfrac{1}{2}\|\nabla k \|^{2}dB\]

where:

  • \(\gamma\) is a Tikhonov regularization parameter

Thickness gradient

\[{\mathcal J\left(k\right)}=\int_{\Omega} \gamma \dfrac{1}{2}\|\nabla H \|^{2}d\Omega\]

where:

  • \(\gamma\) is a Tikhonov regularization parameter

Model parameters

The parameters relevant to the stress balance solution can be displayed by typing:

>> md.inversion
  • md.inversion.iscontrol: 1 if inversion is activated, 0 for a forward run (default)
  • md.inversion.incomplete_adjoint: 1 linear viscosity, 0 non-linear viscosity
  • md.inversion.control_parameters: parameters that are inferred (ex: {'FrictionCoefficient'} or {'MaterialsRheologyBbar'}
  • md.inversion.cost_functions: list of individual cost functions that are summed to calculate the final cost function \(\mathcal J\) to be minimized (ex: [101, 501])
  • md.inversion.cost_functions_coefficients: weight of each individual cost function previously defined for each vertex (more/no weight can be put on certain regions)
  • md.inversion.min_parameters: minimum value for the inferred parameter
  • md.inversion.max_parameters: maximum value for the inferred parameter
  • md.inversion.vx_obs: x component of the surface velocity
  • md.inversion.vy_obs: y component of the surface velocity
  • md.inversion.vel_obs: surface velocity magnitude
  • md.inversion.thickness_obs: measured ice thickness

Minimization algorithms

Depending on the class of md.inversion, several optimization algorithm are available:

  • Brent search algorithm (md.inversion = inversion(), the default)
  • Toolkit for Advanced Optimization (TAO) (md.inversion = taoinversion())
  • M1QN3 algorithm (md.inversion = m1qn3inversion())

Each minimizer has its own optimization parameters described below.

ISSM has an interface to M1QN3 (Inria) [Gilbert1989]. This interface was largely based on [Nardi2009]. Here is a list of the relevant parameters:

  • md.inversion.maxsteps: maximum number of iterations (gradient computation)
  • md.inversion.maxiter: maximum number of Function evaluation (forward run)
  • md.inversion.dxmin: convergence criterion: two points less than dxmin from each other (sup-norm) are considered identical
  • md.inversion.gttol: gradient relative convergence criterion 2 (defined below)

Brent search minimizers

  • md.inversion.nsteps: number of optimization searches (gradient evaluations)
  • md.inversion.maxiter_per_step: maximum iterations during each optimization step
  • md.inversion.step_threshold: decrease threshold for next step (default is 30)
  • md.inversion.gradient_scaling: scaling factor on gradient direction during optimization, for each optimization step
\[\alpha\in\left[0,\mbox{gradient\_scaling} \right]\hspace{3em}p^{\text{new}}=p^{\text{old}}-\alpha\;\nabla_p {\mathcal J}/\|\nabla_p {\mathcal J}\|\]

Toolkit for Advanced Optimization (TAO)

ISSM has an interface to the Toolkit for Advanced Optimization (TAO) [Munson2012]. Here is a list of the relevant parameters:

  • md.inversion.maxsteps: maximum number of iterations (gradient computation)
  • md.inversion.maxiter: maximum number of Function evaluation (forward run)
  • md.inversion.algorithm: inimization algorithm. ex: 'tao_blmvm', 'tao_cg', 'tao_lmvm'
  • md.inversion.fatol: cost function absolute convergence criterion (defined below)
  • md.inversion.frtol: cost function relative convergence criterion (defined below)
  • md.inversion.gatol: gradient absolute convergence criterion (defined below)
  • md.inversion.grtol: gradient relative convergence criterion (defined below)
  • md.inversion.gttol: gradient relative convergence criterion 2 (defined below) with the following convergence criteria:
\[\begin{array}{lcl} f(X) - f(X^*) & < & \epsilon_{fatol} \\ \left|f(X) - f(X^*\right|/\left|f(X^*)\right| & < & \epsilon_{frtol} \\ \|g(X)\| & < & \epsilon_{gatol} \\ \|g(X)\|/\left|f(X)\right| & < & \epsilon_{grtol} \\ \|g(X)\|/\|g(X_0)\| & < & \epsilon_{gttol} \\ \end{array}\]

where:

  • \(f(X)\) is the cost function at \(X\)
  • \(g(X)\) is the cost function gradient with respect to \(X\)
  • \(X^*\) is the estimated “true” minimum
  • \(X_0\) is the initial guess

Running an inversion

To run an inversion, solve a stress balance solution with md.inversion.iscontrol = 1:

>> md = solve(md, 'Stressbalance');

References

  • Jean Charles Gilbert and Claude Lemarechal. Some numerical experiments with variable-storage quasi-Newton algorithms. Math. Program., 45(1-3):407-435, 1989.

  • Todd Munson, Jason Sarich, Stefan Wild, Steven Benson, and Lois Curfman McInnes. TAO 2.0 Users Manual. Technical Report ANL/MCS-TM-322, Mathematics and Computer Science Division, Argonne National Laboratory, 2012. http://www.mcs.anl.gov/tao.

  • Luigi Nardi, Charles Sorror, Fouad Badran, and Sylvie Thiria. YAO: A Software for Variational Data Assimilation Using Numerical Models. In Osvaldo Gervasi, David Taniar, Beniamino Murgante, Antonio Lagana, Youngsong Mun, and Marina L. Gavrilova, editors, LNCS 5593, Computational Science and Its Applications - ICCSA 2009, pages 621-636. Springer-Verlag, 2009.