Inverse Model Results for Filchner-Ronne Catchment

This page contains the results of the inversions for basal drag and drag coefficient in the Filchner-Ronne catchment presented in Wolovick et al., (2023), along with the code used to perform the inversions and L-curves, analyze the results, and produce the figures presented in that paper. 

This all looks very complicated.  There's so many files here.  The description is so long.  I just want to know the basal drag!

If you don't want to get into the weeds of inverse modeling and L-curve analysis, or if you are uninterested in wading through our collection of model structures and scripts, then you should use the file BestCombinedDragEstimate.nc.  That file contains our best weighted mean estimate of the ice sheet basal drag in our domain, along with the weighted standard deviation of the scatter of the different models about the mean.  As discussed in the paper, this combined estimate is constructed from the weighted mean of 24 individual inversions, representing 8 separate L-curve experiments on our highest-resolution mesh, with three regularization values per L-curve (best estimate regularization, along with minimum and maximum acceptable regularization levels).  Each inversion is weighted according to the inverse of its total variance ratio, which is a quality metric incorporating both observational misfit and inverted structure.  For ease of use, these results have been interpolated from the unstructured model mesh onto a 250 m regular grid.  If you only want to know the basal drag in the Filchner-Ronne region, that is the only file you should use.

For users who want to go further, we will now explain the remaining files in this release.  First we give a brief summary of all of the scripts included here and their functions, and then we will give an explanation of the matfiles that contain the actual inversion and L-curve results.  Note that the scripts presented here are the matlab scripts used to organize and set up model runs for ISSM.  The Ice-Sheet and Sea-level System Model (ISSM) is a highly versatile parallelized finite-element ice sheet model run in C but controlled using Matlab or Python front-ends.   We do not include the underlying code for ISSM here; users who are interested in installing ISSM should go to the ISSM home page.  We merely include the Matlab scripts we used to organize our ISSM front-end, set up model structures, and then analyze and visualize results.

Main Matlab scripts:

These are the main functional scripts used to set up and run the model.

• ISSMInversion_v3.m.  This is the primary script we used to set up the inversions and perform L-curve analysis.  It requires a model mesh as input along with some gridded data.  It also produces an L-curve figure (figure 3) after performing the L-curve analysis.  This script can be run in two modes: "setupandsend", which prepares model structures and sends them to the cluster to be solved, and "loadandanalyze", which loads the solutions from the cluster, saves them to matfiles, and performs analysis and visualization.  In addition to L-curve analysis and the L-curve figure, this script can also produce a variety of additional figures of model output that we did not show in the paper.
• MakeISSMMesh_v4.m.  This is the script we used to make our model meshes.  It requires a domain boundary as input along with some gridded data.
• ModelBoundaryPicker_v1.m.  This script opens a crude graphical interface for picking the domain outline.
• OrganizeInversionsForRelease_v2.m.  This script assembles L-curve and inverse model results and organizes them into the data release you see here.  Note that it doesn't compute the combined drag estimate itself, (that is done by CombinedDragFigure_v1.m), but it does interpolate the combined drag estimate from the model mesh to the grid, and it produces the output netcdf file.

Note that the gridded data files needed by some of the above scripts are not included in our release here.  Users interested in using these scripts for their own projects will need to provide their own gridded inputs, for instance from BedMachine or Measures. 

Figure-making scripts:

These scripts produced almost all of the figures we presented in the paper, and also computed the statistics we presented in the tables in the paper.

• CombinedDragFigure_v1.m.  This script computes the combined drag estimate on the highest-resolution mesh, and makes a figure displaying it (Figure 12 in the paper).
• InversionComparisonFigure_HOSSA_v1.m.  This makes figure 11 in the paper and also computes the statistics shown in table 3.
• InversionComparisonFigure_m_v1.m  This makes figures 9 and 10, and also computes the statistics shown in table 2.
• InversionComparisonFigure_N_v1.m.  This makes figure 8, and also computes the statistics shown in table 1.
• InversionComparisonFigure_v1a.m.  This makes figure 4 in the paper.
• InversionResConvergenceFigure_v2.m.  This makes figure 6.
• InversionResMisfitFigure_v1.m.  This makes figure 7.
• InversionSettingFigure_v1.m.  This makes figure 1.
• InversionSpectrumFigure_v1.m.  This performs spectral analysis and makes figure 5.
• InversionThermalSettingFigure_v1.m.  This makes figure A1.
• MeshSizeFigure_v1.m.  This makes figure A2.
• NComparisonFigure_v1.m.  This makes figure 2.

Other utility Matlab functions:

These miscellaneous function do various tasks.  Many of them are called as subroutines of the scripts above.  Additionally, many of them are generally useful in contexts beyond the inverse modeling presented here.

• FlattenModelStructure.m.  ISSM has the unfortunate convention of saving every variable in 3D meshes on every single 3D mesh node, which is quire wasteful for variables that are actually 2D (ie, most of the model variables).  This function flattens all uneccessarily 3D information, but unlike the built-in ISSM function flatten.m, this script preserves the 3D geometry of the mesh, along with 3D variables that actually are 3D (such as englacial temperature, for example). This function can also be run in reverse to expand variables back to full 3D before calling solve().
• intuitive_lowpass.m.  This function low-pass filters a 1D dataset using a gaussian filter.  It has several options for handling boundary conditions at the end points.
• LaplacianInterpolation.m.  This function fills in missing data values for gridded data products by solving Poisson's equation (Laplacian=0).
• LaplacianInterpolation_mesh.m.  This function does the same thing but on an unstructured mesh.
• loadnetcdf.m.  This function loads variables from netcdf files into the Matlab workspace using a similar syntax as load() for matfiles.
• MultiWavelengthInterpolator.m.  This function interpolates gridded data onto an unstructured mesh using a multi-grid approach.  The grid is smoothed at multiple wavelengths and each mesh element interpolates from the wavelength that is appropriate for its size.  This functionality is useful for preventing aliasing in coarse-resolution areas when interpolating onto a mesh with variable mesh size.  It also produces results that are approximately (but not precisely) conservative.
• ThreeByThree.m.  This function iteratively performs a 3x3 smoothing on gridded data.
• unpack.m.  This function takes a structure and "unpacks" it by making every field into a variable in the workspace.

Matfiles with L-curve data and model structures.

The results of our L-curve analyses and our actual inversion results are stored in matfiles.  We performed 21 experiments shown in the paper; for each one we performed an independent L-curve analysis using 25 individual inversions, for a total of 525 inversions.  However, for this data release we simplify matters by only presenting 3 inversions per experiment, corresponding to the best regularization value (LambdaBest) and the maximum and minimum acceptable regularization values (LambdaMax and LambdaMin).  In addition, for each experiment we also provide an LCurveFile that summarizes the L-curve analysis but does not contain any actual model results.  In total, we present 84 matfiles in this data release.

Naming convention:

All matfiles presented here have the following naming convention:

Mesh#_eqn_m#_Ntype_LambdaType.mat

• Mesh#:  this represents the mesh on which the inversions were performed, ranging from Mesh1 (highest resolution) to Mesh10 (lowest resolution).
• eqn:  this represents the type of equations solved in the inversion.  Values are "SSA" or "HO".
• m#:  exponent in the sliding law.  Values are m1, m3, and m5.
• Ntype:  effective pressure source in the sliding law.  Values are "noN" (ie, Weertman sliding), "Nop", "Nopc", and "Ncuas".
• LambdaType:  values of this string are "LCurveFile" (for the file summarizing the whole L-curve experiment), "LambdaMin", "LambdaBest", and "LambdaMax".

Variables in the model files:

Every file ending with "LambdaMin", "LambdaBest", or "LambdaMax" is a model file containing the same set of variables.  Those variables are:

• md.  This is a a model structure variable usable by any ISSM installation.  Note that if you do not have ISSM installed on your machine, Matlab will not recognize class "model" and you will not be able to load this variable.  The results of the inversion are stored in md.results.StressbalanceSolution.  Other important things for the inversion, such as cost functions, cost function coefficients, and velocity observations, are stored in md.inversion.  Note that the process of normalizing cost function components described in the paper was performed in practice by manipulating values of md.inversion.cost_function_coefficients.  Other important fields of the model structure that might be relevant are md.mesh (describing the numerical mesh), md.geometry (ice sheet geometry), md.mask (ice and ocean masks), md.friction (containing the basal sliding law, including the slip exponent and effective pressure field), md.materials (containing material properties, including rheology), md.stressbalance (controlling the stress balance solution), and md.solvers and md.toolkits (numerical solvers).  In addition, users should note that md.cluster is class "ollie", which is a custom class made for AWI's old HPC setup.  This will produce a warning that the class of md.cluster is unrecognized if you load md on another computer, and you will need to replace md.cluster with something appropriate to your own setup if you want to use md in the "solve()" command.  Also note that md.miscellaneous.dummy contains a copy of the variable "GoodData_obs_vertices", a data quality metric produced by the script ISSMInversion_v3.m.   Finally, users should note that md.miscellaneous.name is still set to the original file name we used internally, not the cleaner names we have used for this release.
• MyInversionParameters.  This is a structure containing the parameters used by ISSMInversion_v3.m to set up and L-curve and inverse model.
• oldfilename.  This is the original file name that we used when we generated these results internally.  The cleaner file names you see here were created by the script OrganizeInversionsForRelease_v2.m
• thislambda.  The regularization lambda value used for this inversion.

Variables in the L-curve files:

The files ending in "LCurveFiles" contain data summarizing the results of each L-curve experiment.  Those variables are:

• Lambda.  This variable contains the 25 values of the regularization parameter lambda tested in the L-curve.  Size: [25,1].
• DataCost, RegularizationCost.  These two variables contain the values of J_obs and J_reg for each of the 25 inversions in the L-curve.  Size: [25,1].
• IsOutlier.  This variable identifies whether any of the 25 inversions have been identified as an outlier and excluded from the curve-fitting process.  Size: [25,1].
• LambdaWavelengths.  This variable contains all 50 of the smoothing wavelengths tested during the curve-fitting process.  Note that these are wavelengths in log(lambda) space.  Size: [1,50].
• LogLambda_curves.  This variable contains log(lambda), interpolated to fine spacing.  All variables on the smooth curves are defined on a sampling defined by this variable.  Size: [1000,1].
• LogDataCost_curves, LogRegCost_curves.  These variables contain all of the candidate smoothed curves of log(J_obs) and log(J_reg). Size: [1000,50].
• CurvatureVariance_curves, ScatterVariance_curves.  Variance measures (Eqs A1 and A2 in the paper) used to quantify how close the smooth curves are to fitting the 25 inversion data points (ScatterVariance) and how much variability the smooth curves have (CurvatureVariance).  Size: [1,50].
• TotalVariance_curves.  Total variance metric used to select optimal smoothing wavelength.  The optimal smoothing wavelength is defined by the minimum of this variable.  Size: [1,50].
• LogDataCost_bestcurve, LogRegCost_bestcurve.  These variables contain the final best smoothed curve of log(J_obs) and log(J_reg). Size: [1000,1].
• LogDataGradient_bestcurve, LogRegGradient_bestcurve.  These variables contain the first derivatives, d(log(J))/d(log(lambda)), for both J_obs and J_reg, computed for the best smoothed curve.  Note that their size has been reduced by 1 relative to the original curves (these values should be placed at the midpoints between adjacent samples in LogLambda_curves).  Size: [999,1].
• LogDataCurvature_bestcurve, LogRegCurvature_bestcurve.  These variables contain the second derivatives, d^2(log(J))/d(log(lambda))^2, for both J_obs and J_reg, computed for the best smoothed curve.  Note that their size has been reduced by 2 relative to the original curves (these values should be placed at the samples in LogLambda_curves, excluding the first and last samples).  Size: [998,1].
• TotalCurvature_bestcurve.  This variable corresponds to the sum LogDataCurvature_bestcurve+LogRegCurvature_bestcurve.  The best lambda value is the peak of this variable and the min/max acceptable values are defined where this variable drops to half of its maximum value.  Note that, while the fine lambda sampling of the smoothed curve allows us to define these lambda values with some precision, we still need to select models from the original 25 lambda values, which have much coarser resolution in log(lambda) space.  The "LambdaBest" model is selected to be the model that is closest to the fine-resolution lambdabest in log(lambda) space, while the LambdaMin and LambdaMax models are selected to be the first and last models inside those thresholds.  Size: [998,1].
• bestwavelength, bestwavelengthind.  These variables describe the best smoothing wavelength (wrt log(lambda)) used to generate the best smooth curve.  bestwavelengthind is the index of this wavelength within LambdaWavelengths.  Size: [1,1].
• bestlambda, bestlambda_minallowable, bestlambda_maxallowable.  These three variables represent the continuously determined best, min, and max lambda values.  Size: [1,1].
• bestlambdaind, bestlambdaind_minallowable, bestlambdaind_maxallowable. As above, but these represent the index of the lambda values within LogLambda_curves.  Size: [1,1].
• bestmodelind, bestmodelind_minallowable, bestmodelind_maxallowable.  The index of the models selected to represent the corresponding best lambda values.  Note that for the bestlambda, this is the model with log(lambda) closest to the best value, but for the min/max acceptable, this is the first or last model within the relevant bound.  Indices refer to the list of 25 models in each L-curve.  Size: [1,1].
• bestlambda_loguncertainty.  The error bar around the best lambda value in log(lambda) space.  This is equal to half of bestwavelength.  The corresponding uncertainty ratio on a linear scale is exp(bestlambda_loguncertainty).  Note that this same uncertainty value also applies to minimum and maximum acceptable lambda.  Size: [1,1].
• MyInversionParameters.  A structure containing the parameters used by ISSMInversion_v3.m to set up and the L-curve and inversions.
• LCurveFittingParameters.  A structure containing the parameters used by ISSMInversion_v3.m to perform the curve-fitting analysis.
• oldfilename.  The original file name of this L-curve in our internal system before we gave the files neater names for public release.