Author: Aaron Stubblefield (Columbia University, LDEO).
This program inverts for (1) basal velocity anomaly "w" and/or (2) the basal drag coefficient "beta" given the observed surface elevation change data "h_obs" (and possible horizontal surface velocity data) by solving a least-squares minimization problem
The main model assumptions are (1) Newtonian viscous ice flow, (2) a linear basal sliding law, and (3) that all fields are small perturbations of a simple background flow. These assumptions allow for efficient solution of the forward model: 2D (map-plane) Fourier transforms and convolution in time are the main operations. The model and numerical method are described in a forthcoming manuscript.
As of this commit, this code runs with the latest SciPy (https://www.scipy.org/) release. Plotting relies on Matplotlib (https://matplotlib.org/).
FFmpeg (https://www.ffmpeg.org/) can be used to create a video of the results.
The model is organized in 9 python files in the code directory as follows.
-
main.py is the main file that calls the inverse problem solver and then plots the solution.
-
inversion.py is the inverse problem solver: this defines the normal equations and calls the conjugate gradient solver.
-
conj_grad.py defines the conjugate gradient method that is used to solve the normal equations.
-
operators.py defines the forward and adjoint operators that appear in the normal equations.
-
kernel_fcns.py defines the relaxation and transfer functions that the forward and adjoint operators depend on.
-
regularizations.py defines the regularization options (L2 and H1).
-
params.py defines all of the model options and parameters.
-
synthetic_data.py defines synthetic elevation anomalies for the test problems.
-
plotting.py creates png images of the inversion.
To run the test problems, see the notebook files in the notebooks directory. These examples can be run in Binder:
The primary inversion options are passed to the main function in main.py.
The main options are whether to invert for the basal vertical velocity anomaly "w" or the drag coefficient "beta" (or both), and whether surface velocity data is included. The regularization parameters are also passed to the main function.
Regularization type (H1 or L2), numerical parameters, and physical parameters are set in params.py.
Synthetic data for the test problems can be set/modified in synthetic_data.py.
The synthetic data is created by solving the forward problem (given w or beta)
for the elevation anomaly, and then adding
some noise. The added noise is proportional to the norm of the elevation anomaly,
scaled by the noise_level
parameter.