Publications
Global mapping of lake-terminating glaciers
Jakob Steiner, William Armstrong, Will Kochtitzky, Robert McNabb, Rodrigo Aguayo, Tobias Bolch, Fabien Maussion, Vibhor Agarwal, Iestyn Barr, Nathaniel R. Baurley, Mike Cloutier, Katelyn DeWater, Frank Donachie, Yoann Drocourt, Siddhi Garg, Gunjan Joshi, Byron Guzman, Stanislav Kutuzov, Thomas Loriaux, Caleb Mathias, Brian Menounos, Evan Miles, Aleksandra Osika, Kaleigh Potter, Adina Racoviteanu, Brianna Rick, Miles Sterner, Guy D. Tallentire, Levan Tielidze, Rebecca White, Kunpeng Wu, and Whyjay Zheng
Earth System Science Data, 18, 1665–1681 | 2026
Journal Paper

Proglacial lakes at glacier termini have received widespread attention in the literature for their role in accelerating melt, velocity and contributing to cryospheric hazards. Although global and regional inventories for both glaciers and lakes exist, lake-terminating glaciers have not been consistently identified at the global scale. Based on the most recent global glacier inventory (RGI 7.0), which so far identifies marine-terminating glaciers in most regions, but not lake-terminating glaciers, we present a global inventory of lake-terminating glaciers, differentiating between three categories based on the degree of contact between the glacier terminus and any proglacial lakes. Contributors manually assigned categories to glaciers using satellite imagery from as close to the target date of 2000 as possible, aided by regional lake inventories where available. The resulting dataset corresponds to the year 2000 (±1.5), matching to the timestamp of RGI 7.0 outlines (2001 ± 6.2). We find that of 274 531 glaciers worldwide, 1.4 % terminate in lakes, with regional percentages varying between 0.5 % and 6.7 % across the 19 RGI regions. These glaciers account for 11.4 % of global glacier area (0.2 % to 41.8 % across regions). With multiple submissions available for 1260 individual glaciers, we find mapping conflicts between contributors to be low (6.7 %). The lake termini data set is available at https://doi.org/10.5281/zenodo.15524733 (Steiner et al., 2025) as well as at https://github.com/GLIMS-RGI/lake_terminating (last access: 2 February 2026). This dataset is integrated into the forthcoming update to the RGI, v7.1.

Investigating seasonal velocity variations of selected glaciers in high mountain asia
Francesca Baldacchino, Whyjay Zheng, Kunpeng Wu, Vassiliy Kapitsa, Alexandr Yegorov, and Tobias Bolch
Science of Remote Sensing, 12, 100266 | 2025
Journal Paper

Glacier velocity is a sensitive indicator of mass balance and is key to understanding how glaciers respond to climate change. Monitoring glacier velocity at high temporal resolutions enables a better understanding of the drivers of glacier dynamics. Previous studies have found that the glaciers in High Mountain Asia (HMA) tend to slow down concomitant to losing mass at an accelerating rate on decadal timescales. However, few studies have explored seasonal variations in glacier velocities and have typically focused on large, fast-flowing glaciers. We select one debris-covered glacier, and four clean-ice glaciers in HMA. Sentinel-1 and -2 images are used to calculate the glacial velocities using the feature tracking module provided by the Cryosphere And Remote Sensing Toolkit (CARST). We develop a novel, regularised linear inverse model to extract the seasonally resolved glacial velocity time series (6-day intervals) with rigorous uncertainty estimates. Our results show that three of the five glaciers have strong seasonal signals in velocities, with faster velocities in spring and/or summer compared to winter. We also find an up-glacier propagation of the late spring and/or summer accelera-tions and a down-glacier propagation of the autumn accelerations. We suggest that changes in the subglacial hydrology efficiency drive the observed seasonal variations in velocities. We also highlight that icefalls may alter the glacier flow response by blocking the development of subglacial drainage channels and thus the seasonal propagation of velocities. Our methodol-ogy enables us to successfully extract seasonal signals in glaciers that flow slowly and provide a further understanding of glacier dynamics.

Community estimate of global glacier mass changes from 2000 to 2023
The GlaMBIE Team (incl. Whyjay Zheng)
Nature, 639(8054), 382–288 | 2025
Journal Paper

Abstract

Glaciers are indicators of ongoing anthropogenic climate change. Their melting leads to increased local geohazards, and impacts marine and terrestrial ecosystems, regional freshwater resources, and both global water and energy cycles. Together with the Greenland and Antarctic ice sheets, glaciers are essential drivers of present and future sea-level rise. Previous assessments of global glacier mass changes have been hampered by spatial and temporal limitations and the heterogeneity of existing data series. Here we show in an intercomparison exercise that glaciers worldwide lost 273 ± 16 gigatonnes in mass annually from 2000 to 2023, with an increase of 36 ± 10% from the first (2000–2011) to the second (2012–2023) half of the period. Since 2000, glaciers have lost between 2% and 39% of their ice regionally and about 5% globally. Glacier mass loss is about 18% larger than the loss from the Greenland Ice Sheet and more than twice that from the Antarctic Ice Sheet. Our results arise from a scientific community effort to collect, homogenize, combine and analyse glacier mass changes from in situ and remote-sensing observations. Although our estimates are in agreement with findings from previous assessments at a global scale, we found some large regional deviations owing to systematic differences among observation methods. Our results provide a refined baseline for better understanding observational differences and for calibrating model ensembles, which will help to narrow projection uncertainty for the twenty-first century.

GLAcier Feature Tracking testkit (GLAFT): a statistically and physically based framework for evaluating glacier velocity products derived from optical satellite image feature tracking
Whyjay Zheng, Shashank Bhushan, Maximillian Van Wyk De Vries, William Kochtitzky, David Shean, Luke Copland, Christine Dow, Renette Jones-Ivey, and Fernando Pérez
The Cryosphere, 17, 4063–4078 | 2023
Journal Paper
Software

Abstract

Glacier velocity measurements are essential to understand ice flow mechanics, monitor natural hazards, and make accurate projections of future sea-level rise. Despite these important applications, the method most commonly used to derive glacier velocity maps, feature tracking, relies on empirical parameter choices that rarely account for glacier physics or uncertainty. Here we test two statistics- and physics-based metrics to evaluate velocity maps derived from optical satellite images of Kaskawulsh Glacier, Yukon, Canada, using a range of existing feature-tracking workflows. Based on inter-comparisons with ground truth data, velocity maps with metrics falling within our recommended ranges contain fewer erroneous measurements and more spatially correlated noise than velocity maps with metrics that deviate from those ranges. Thus, these metric ranges are suitable for refining feature-tracking workflows and evaluating the resulting velocity products. We have released an open-source software package for computing and visualizing these metrics, the GLAcier Feature Tracking testkit (GLAFT).

Short summary

We design and propose a method that can evaluate the quality of glacier velocity maps. The method includes two numbers that we can calculate for each velocity map. Based on statistics and ice flow physics, velocity maps with numbers close to the recommended values are considered to have good quality. We test the method using the data from Kaskawulsh Glacier, Canada, and release an open-sourced software tool called GLAcier Feature Tracking testkit (GLAFT) to help users assess their velocity maps.

Glacier geometry and flow speed determine how Arctic marine-terminating glaciers respond to lubricated beds
Whyjay Zheng
The Cryosphere, 16, 1431–1445 | 2022
Journal Paper

Abstract

Basal conditions directly control the glacier sliding rate and the dynamic discharge of ice. Recent glacier destabilization events indicate that some marine-terminating glaciers quickly respond to lubricated beds with increased flow speed, but the underlying physics, especially how this vulnerability relates to glacier geometry and flow characteristics, remains unclear. This paper presents a 1D physical framework for glacier dynamic vulnerability assuming sudden basal lubrication as an initial perturbation. In this new model, two quantities determine the scale and the areal extent of the subsequent thinning and acceleration after the bed is lubricated: Péclet number (Pe) and the product of glacier speed and thickness gradient (dubbed J0 in this study). To validate the model, this paper calculates Pe and J0 using multi-sourced data from 1996 to 1998 for outlet glaciers in the Greenland ice sheet and Austfonna ice cap, Svalbard, and compares the results with the glacier speed change during 1996/1998–2018. Glaciers with lower Pe and J0 are more likely to accelerate during this 20-year span than those with higher Pe and J0, which matches the model prediction. A combined factor of ice thickness, surface slope, and initial flow speed physically determines how much and how fast glaciers respond to lubricated beds in terms of speed, elevation, and terminus change.

Short Summary

A glacier can speed up when surface water reaches the glacier's bottom via crevasses and reduces sliding friction. This paper builds up a physical model and finds that thick and fast-flowing glaciers are sensitive to this friction disruption. The data from Greenland and Austfonna (Svalbard) glaciers over 20 years support the model prediction. To estimate the projected sea-level rise better, these sensitive glaciers should be frequently monitored for potential future instabilities.

Mapping ice flow velocity using an easy and interactive feature tracking workflow
Whyjay Zheng, Shane Grigsby, Facundo Sapienza, Jonathan Taylor, Tasha Snow, Fernando Pérez, and Matthew Siegfried
2021 EarthCube Annual Meeting | Jun 15-17, 2021
Conference Presentation
Other

An open-source, interactive, and highly customizable workflow for glacier velocity mapping using feature tracking technique and satellite images.

Aseismic Deformation During the 2014 Mw 5.2 Karonga Earthquake, Malawi, From Satellite Interferometry and Earthquake Source Mechanisms
Whyjay Zheng, Sarah Jaye Oliva, Cynthia Ebinger, and Matthew Pritchard
Geophysical Research Letters, 47, e2020GL090930 | 2020
Journal Paper

Abstract

Aseismic deformation has been suggested as a mechanism to release the accumulated strain in rifts. However, the fraction and the spatial distribution of the aseismic strain are poorly constrained during amagmatic episodes. Using Sentinel-1 interferograms, we identify the surface deformation associated with the 2014 Mw5.2 Karonga earthquake, Malawi, and perform inversions for fault geometry. We also analyze aftershocks and find a variety of source mechanisms within short timescales. A significant discrepancy in the earthquake depth determined by geodesy (3–6 km) and seismology (11–13 km) exists, although both methods indicate Mw5.2. We propose that the surface deformation is caused by aseismic slip from a shallow depth. This vertical partitioning from seismic to aseismic strain is accommodated by intersecting dilatational faults in the shallow upper crust and sedimentary basin, highlighting the importance of considering aseismic deformation in active tectonics and time-averaged strain patterns, even in rifts with little volcanism.

Plain Language Summary

During the early stages of continental rifting, some accumulated tectonic energy is released with little or no earthquake activity (i.e., aseismic strain). This is most often observed in places where magma intrusion occurs. However, in rift basins lacking evidence of magma intrusion, how much and where the aseismic strain is released remains largely unknown. We investigate a magnitude 5.2 earthquake that occurred in 2014 near Karonga, Malawi, using both local seismic data and satellite ground deformation data. Our analyses show that although the main earthquake and other aftershocks ruptured at ∼12 km deep, the surface deformation is caused by an aseismic source from a shallower region (3–6 km deep). The aftershocks occur seconds to minutes apart and have different rupture sources, suggesting a highly damaged zone. The depth discrepancy between earthquake locations and the aseismic energy source might be caused by a weak, damaged layer where many faults intersect. Since intersecting faults are common in rift systems, our study suggests that this depth discrepancy might be common as well, indicating a large portion of energy released as aseismic strain in a continental rifting system.

Investigating Mass Loss and Changing Ice Dynamics of Arctic Ice Caps Using Remote Sensing
Whyjay Zheng
Doctoral Dissertation, Cornell University | 2020
Other

Abstract

Glacier thinning and retreat have accelerated globally in the last century and are the largest contributor to rising sea levels. For the Arctic region, observations and modeling results have shown that extensive warming is taking place. However, the recent glacier dynamics (mass balance and ice discharge) in many Arctic regions have not been well studied due to the remote nature of these glaciers. This thesis uses multiple types of satellite data to quantify the mass balance and ice discharge for three Arctic regions showing dramatic glacier change in recent decades possibly due to Arctic warming. The objective is to resolve the mass budget and velocity pattern on a per glacier basis and understand the mechanisms driving recent changes. To facilitate the entire workflow, our research team has developed the Cryosphere and Remote Sensing Toolkit (CARST) software, and I am the lead author. CARST provides useful python and bash scripts that use satellite imagery, particularly SAR and optical images, to monitor changes of glaciers and ice caps through time. The first study area is Franz Josef Land (FJL), Russia, which is currently subjected to a rapidly-warming climate in the Arctic. I combine surface elevation data derived from different sources and times, including the WorldView satellite series and the ArcticDEM data set (2011–2015), SPOT-5 (2007), CryoSat-2 (2011–2015), and a digitized cartographic map (1953). I calculate elevation change rate (dh/dt) in two different periods, and the results show a two-fold rate of ice loss over the past 60 years, from -2.18 ± 0.72 Gt/yr (1953–2011/2015 average) to -4.43 ± 0.78 Gt/yr (2011–2015). Despite being spatially variable, a trend of increased thinning from NE towards SW is discovered, suggesting a link to the local gradient in temperature and precipitation. Ice loss is mostly focused on marine-terminating glaciers probably due to the interaction between glaciers and warming ocean water. These retreating glaciers generated a new island in 2016 and more islands are likely to emerge in the foreseeable future as FJL’s ice loss has reached an unprecedented rate. The research focus in the following chapter shifts to the neighboring archipelago called Severnaya Zemlya, Russia. A surge-like collapse initiated in 2013 in Vavilov Ice Cap, one of many ice caps in this region. By spring 2019, this ice cap had lost 9.5 Gt of ice. Using time series of surface elevation and glacier velocity derived from multiple satellite data sets such as WorldView (elevation), ArcticDEM (elevation), ASTER (elevation), Landsat 8 (velocity), Sentinel-1, (velocity), Sentinel-2 (velocity), Radarsat-2 (velocity), and ALOS-2 (velocity), I identify a shift of flow pattern starting in 2017 when shear margins formed within the grounded marine piedmont fan. Multiple summer speedups occurred after the new flow pattern formed, possibly with the aid of basal lubrication due to surface melt. With the analysis using multiple physical models, it is suggested that the collapsed ice cap has entered a new ice stream-like regime in which diffusion of surface thinning controls the glacier dynamics. This is the first documented case of an ice stream-like feature ever being formed, and this glacier now flows at a higher speed and drains the ice cap more efficiently. To publicize the findings and their scientific implications, I made two videos showing the temporal changes of the terminus position and speed pattern, which are available on Youtube. In the last chapter, I further develop a physical framework for the glacier perturbation model to understand how different glaciers respond to basal lubrication. The modified 1-D flowline model suggests two physical quantities, Péclet number (Pe) and a value dubbed J0, governing glacier vulnerability to basal lubrication. To test the model, I use the Ice Thickness Models Intercomparison eXperiment (ITMIX) data set and the NASA MEaSUREs ITS_LIVE data set. ITMIX contains velocity, elevation, and ice thickness data from Austfonna Ice Cap, Svalbard, where multiple glacier collapse events occurred within the past 10 years. I calculate Pe and |J0| using the data from ITMIX and compare them with the speed change revealed by the ITS_LIVE data set. The results show that a low Pe and a high |J0| correspond to the high magnitude of glacier speedup during 1995–2018, as suggested by the model prediction. My analysis implies that basal lubrication can lead to a prolonged or even permanent change of glacier dynamics for some glaciers. These "weak" glaciers might be able to waste ice more rapidly than we thought, posing a warning of an underestimated sea level rise projection.

The possible transition from glacial surge to ice stream on Vavilov Ice Cap
Whyjay Zheng, Matthew Pritchard, Michael Willis, and Leigh Stearns
Geophysical Research Letters, 46, 13892–13902 | 2019
Journal Paper

Abstract

Surge-type glaciers typically undergo cyclical flow instability due to mass accumulation; however, some recent glacier surges have caused irreversible ice loss in a short period. At Vavilov Ice Cap, Russia, surge-like behavior initiated in 2013 and by spring 2019 the ice cap had lost 9.5 Gt of ice (11% mass of the entire basin). Using time series of surface elevation and glacier velocity derived from satellite optical and synthetic-aperture radar imagery, we identify a shift of flow pattern starting in 2017 when shear margins formed within the grounded marine piedmont fan. Multiple summer speedups correlate with warmer summers during 2015–2019 and suggest that surface melt may access the subglacial environment. Force balance analysis and examination of the Péclet number show that glacier thinning propagated upstream in 2016–2017, and diffusion became a significant dynamic response to thinning perturbations. Our results suggest that the glacier has entered a new ice stream-like regime.

Plain Language Summary

A glacier surge is a sudden speedup of glacier flow coinciding with a large advance of the ice front. Some glaciers surge periodically every 10–100 years, and so surge mechanism is thought to be independent of climate change. However, some recent surges have evacuated so much ice that another surge is unlikely to occur in the same place again. A glacier surge at the Vavilov Ice Cap, Russia, is one of these cases. Since 2013, the glacier has drained more than 11% of the ice basin (9.5 Gt) into the ocean. After the initial surge in 2013, the glacier still retains fast flow at around 1.8 km/year, an unusually high and long-lasting speed for a glacier surge. To understand the future of the surge, we use satellite images to calculate surface elevations and ice speeds, and analyze their change over time for the glacier. The results reveal that the glacier now physically behaves more like an "ice stream," a stream of fast-flowing ice within an ice sheet, which can probably flow at high speed for a long time and drain ice efficiently. This is the first documented case of an ice stream-like feature ever being formed.

Faulting processes during early-stage rifting: seismic and geodetic analysis of the 2009–2010 Northern Malawi earthquake sequence
James Gaherty, Whyjay Zheng, Donna Shillington, Matthew Pritchard, Scott Henderson, Patrick Chindandali, Hassan Mdala, A. Shuler, N. Lindsey, Sarah Jaye Oliva, Scott Nooner, C. Scholz, D. Schaff, G. Ekström, and Meredith Nettles
Geophysical Journal International, 217(3), 1767–1782 | 2019
Journal Paper

Summary

In December, 2009, a rare sequence of earthquakes initiated within the weakly extended Western Rift of the East African Rift system in the Karonga province of northern Malawi, providing a unique opportunity to characterize active deformation associated with intrabasinal faults in an early-stage rift. We combine teleseismic and regional seismic recordings of the largest events, InSAR imagery of the primary sequence, and recordings of aftershocks from a temporary (4-month) local network of six seismometers to delineate the extent and geometry of faulting. The locations of ∼1900 aftershocks recorded between January and May 2010 are largely consistent with a west-dipping normal fault directly beneath Karonga as constrained by InSAR and CMT fault solutions. However, a substantial number of epicentres cluster in an east-dipping geometry in the central part of the study area, and additional west-dipping clusters can be discerned near the shore of Lake Malawi, particularly in the southern part of the study area. Given the extensive network of hanging wall faults mapped in the Karonga region on the surface and in seismic reflection images, the distribution of events is strongly suggestive of multiple faults interacting to produce the observed deformation, and the InSAR data permit this but do not require it. We propose that fault interaction contributed to the seismic moment release as a series of Mw 5-to-6 events instead of a normal main shock–aftershock sequence. We find the depth of fault slip during the main shocks constrained by InSAR peaks at less than 6 km, while the majority of recorded aftershocks are deeper than 6 km. This depth discrepancy appears to be robust and may be explained by fault interaction. Structural complexities associated with fault interaction may have limited the extent of coseismic slip during the main shocks, which increased stress deeper than the coseismic slip zone on the primary fault and synthetic faults to the east, causing the energetic aftershock series. There is no evidence of deformation at the Rungwe volcanic province ∼50 km north of the earthquake sequence between 2007 and 2010, consistent with previous interpretations of no significant magmatic contribution during the sequence.

Cryosphere And Remote Sensing Toolkit (CARST) v1.0.1
Whyjay Zheng, William Durkin, Andrew Melkonian, and Matthew Pritchard
Zenodo | 2019
Software
Other

Cryosphere And Remote Sensing Toolkit (CARST) is a toolkit that provides useful python and bash scripts that use satellite imagery, particularly SAR and optical images, to monitor changes of glaciers and ice caps through time. The toolkit has two main approaches: ice elevation changes (also known as dh/dt) and pixel tracking. The version 1.0.1 has been released and a few small bugs have been fixed from v1.0. Please see README file for more information.

Massive destabilization of an Arctic ice cap
Michael Willis, Whyjay Zheng, William Durkin, Matthew Pritchard, Joan Ramage, Julian Dowdeswell, Toby Benham, Robin Bassford, Leigh Stearns, Andrey Glazovsky, Yuri Macheret, and Claire Porter
Earth and Planetary Science Letters, 502, 146–155 | 2018
Journal Paper

Abstract

Ice caps that are mostly frozen at the bedrock-ice interface are thought to be stable and respond slowly to changes in climate. We use remote sensing to measure velocity and thickness changes that occur when the margin of the largely cold-based Vavilov Ice Cap in the Russian High Arctic advances over weak marine sediments. We show that cold-based to polythermal glacier systems with no previous history of surging may evolve with unexpected and unprecedented speed when their basal boundary conditions change, resulting in very large dynamic ice mass losses (an increase in annual mass loss by a factor of ∼100) over a few years. We question the future long-term stability of cold and polythermal polar ice caps, many of which terminate in marine waters as the climate becomes warmer and wetter in the polar regions.

Accelerating glacier mass loss on Franz Josef Land, Russian Arctic
Whyjay Zheng, Matthew Pritchard, Michael Willis, Paul Tepes, Noel Gourmelen, Toby Benham, and Julian Dowdeswell
Remote Sensing of Environment, 211, 357–375 | 2018
Journal Paper

Abstract

The glaciers of the Franz Josef Land (FJL) archipelago in the Russian Arctic are subjected to rapidly-warming temperatures but are small contributors to sea level. We analyze ice surface elevation data derived from satellite stereo imagery (WorldView and SPOT), radar altimetry (CryoSat-2), and a digitized 1953 cartographic map to calculate elevation change rates (dh/dt). Mass loss from FJL doubled between 2011 and 2015 compared to 1953–2011/2015, increasing from a rate of −2.18 ± 0.72 Gt yr−1 to −4.43 ± 0.78 Gt yr−1. This 2011−2015 rate indicates an acceleration in ice loss from that observed in 2003–2009 by multiple studies using ICESat and GRACE. Glacier thinning rates are spatially highly variable. We observe glacier thinning rates of up to 10 m per year, and in general we see a trend of increased thinning from the NE towards the SW. Glacier retreat is widespread and has led to the creation of at least one new island. Historically, ice wastage from FJL is thought to have been relatively small, but accelerating ice loss may be the new normal for this archipelago in a warming Arctic.

Elastic Flexure Model of Iapetus' Equatorial Ridge
Whyjay Zheng
Master Thesis, National Taiwan University | 2013
Other

Abstract

Iapetus may be the most peculiar satellite in the Solar System. This Saturnian moon has a mean radius of 735 km, but an averagely 10-kilometer-high mountain ridge lies precisely on its 75% equatorial circumference. The ridge is so high that Iapetus appears walnut shaped, and it is named “equatorial ridge” after this amazing truth. The ridge was discovered by the Cassini spacecraft in 2005, but the formation theory is still under debate because of the lack of observational data. Several hypotheses, which are roughly divided into endogenic (tectonic buckling) and exogenic (ring remnant) processes, are attributed to explain its origin. Previous studies also noted that the shape of Iapetus is an oblate spheroid related to a hydrostatic spin period of 16 h, but Iapetus now is tidally synchronized with a 79-day period. Because the surface of Iapetus is old and heavily cratered, the formation of the ridge and the oblate spheroid had finished in the early stage of Iapetus (> 4000 Ma). Thus, it’s plausible to assume that Iapetus had a high thermal flux when the equatorial ridge formed. The assumption leads to a result that the surface would bend when the applying load like the ridge exerted. Therefore, upon calculating the deflection of the surface, we could obtain some constraints for the thermal history of Iapetus, and the proper origin model of the equatorial ridge. According to this idea, we attempted to construct analytical and numerical flexural models of the equatorial ridge by utilizing elastic lithosphere theory. The equatorial ridge is treated as a perfectly linear load on Iapetus’ hard shell (i.e. elastic layer of Iapetus). The Digital Terrain Model (DTM) data are inputted and transformed to a vertical load function, and also reveals that large deflection exists in some foothills area. This few-kilometre deflection implies a very thin elastic layer enough to regard it as a flat plate. Moreover, there are no tectonic signals on Iapetus, so the flat-Earth and one-plate condition could adapt to the flexure model. To obtain an analytical solution, the equatorial ridge is simplified to a central loading point. This can be rearranged into an explicit deflecting function in the 1-D coordinate system, so the deflection can be computed if the elastic thickness is given. In the numerical model, the point vertical force is replaced by a loading map. The finite difference method is used to solve the ODE flexural function. Consider the elastic thickness may vary with different areas; we also set a variable-thickness program for the numerical modelling. The modelling results illustrate that an over 100-km elastic layer would not cause any significant deflection; it coincides with the previous suggested. However, a deflecting curve with a range of 5-10 km elastic thickness well fits the terrain data, especially for the distance between a bulge and the ridge. Numerical solution also shows that there are 2 factors contributing the geomorphological changes: cratering and the flexure. Cratering created a deep hole and a thinner elastic layer. These new results seem controversial to the previous studies, but the modelled surface profile is highly consistent with numerical ridge DTM profile except the plateau regions which are suspected to be caused by cratering end load pressure. Such a thin shell implies that the ridge formed when the heat flux stayed high (~18 mW m-2). Therefore, the formation of the ridge probably happened before the despin (oblate shaping) event. The thin-layer flexure model also solves the problem of the angle of response because the ridge sank in the deflected surface, lowered the slope from the angle of response to the observed slope of the ridge. Since there is no evidence relating to endogenic processes, the exogenic origin is in favour. In conclusion, the flexural model of Iapetus' equatorial ridge reveals the possibility of thinner (5-10 km) hard shell, fits the surface profile and thermal history, and supplies more clues to the origin of Iapetus, the interesting satellite in the Solar System.