Development of ice sheet model using COMSOL Multiphysics®: Comparison with Ice Sheet System Model (ISSM)

Glacier mass loss related to ice flow is a quantitative factor that controls sea level rise and global warming. Ice flow is attributed to the complicated interaction between grounding line migration, calving basal friction, and topography. Previous studies using massive geophysical observation, including IPR (Ice Penetrating Radar), InSAR (Interferometric Synthetic Aperture Rada), and heat flow exploration, highlight the need to develop a sophisticated numerical model. Because it is difficult to numerically solve the full-stokes equation obtaining ice velocities, the simplified governing equation (i.e., HO, Higher Order model; SSA, Shallow Shelf Approximation; SIA, Shallow Ice Approximation) is widely used in the ice sheet dynamics community. The SSA approach, which assumes a depth-independent velocity model, has computational cost reduction based on simplified descriptions of the full-stokes equation. Here we developed the two-dimensional SSA numerical model to better understand ice dynamics using COMSOL Multiphysics^® (hereafter COMSOL), which is a user-friendly finite element software providing convenient GUI, mesh generation, post-processing tools, various types of elements, and the order of shape function. To verify the application of COMSOL to ice sheet dynamics, we compared it with an open-source finite element package, Ice Sheet System Model (ISSM). The concise toy model successfully simulated the distribution of viscosity and velocity and the evolution of surface topography in both COMSOL and ISSM. Furthermore, we applied realistic bathymetry data (i.e., Bedmap2) of Pine Island Glacier, where the large ice mass loss occurs, to clarify more similarity of each model. The surface velocity patterns calculated by COMSOL and ISSM are significantly similar for various physical properties (e.g., ice viscosity, friction, and rate of accumulation and melt). We propose that COMSOL, which efficiently handles mesh generation and visualization, and various weak forms, can sufficiently be applied in ice dynamics.