An Empirical Estimate of GIA-Induced Vertical Motion in the Great Lakes Basin Derived from an Ensemble of GIA Models

Vertical land motion in the Great Lakes Basin (GLB) arises from the combined effects of ongoing Glacial Isostatic Adjustment (GIA) and shorter-term environmental and hydrological loadings. Because the present-day GIA hinge line crosses the region, even small errors in separating long-term uplift from elastic responses can strongly bias geophysical interpretations. Over the past two decades, the GLB has experienced pronounced lake-level fluctuations. An analysis of GRACE/FO data indicates minimal Total Water Storage (TWS) changes across the GLB during 2002-2012, a period during which lake levels were relatively stable and vertical motions should therefore reflect GIA alone. In contrast, from 2012 to 2019 lake levels rose to record highs, and since 2020 they have been falling at a comparable rate. The area is densely instrumented with continuous GNSS stations, providing an exceptional opportunity to investigate how long-term GIA and short-term hydrological forcing interact. Our goal is to develop a robust, precise and accurate estimate of the GIA signal so that we can accurately remove GIA from observations and constrain surface/groundwater storage changes.

We compared the predictions of many GIA models with the pre-2012 observations, which should reflect the GIA signal alone, but none of the existing models adequately reproduce the observed data. Despite differences in viscosity structure or ice history, every model produces the same systematic bias: the hinge line (zero uplift) is positioned too far south. However, the shape of the modelled profiles matches the GNSS curvature extremely well.  Therefore, we developed a spatial optimization framework to minimize geometric misalignments between GIA model predictions and GNSS vertical velocities across the Great Lakes (2002–2012.5). Seventy-two GIA realizations based on diverse ice histories (ICE-6G_C, ICE-6G_D, ICE-7G_NA, ANU-ICE, NAICE, etc.) and Earth rheologies were subjected to systematic horizontal translations (via a grid search with limits ranging from ±2° to ±8°), with and without allowing for small planar rotations, yielding 576 model-configuration combinations evaluated using RMS misfit, concordance correlation, and 5-fold cross-validation. The best-fitting models achieve the lowest misfits (approximately 0.40 mm/yr), and highest concordance (ccc > 0.90). The models that fit well give very consistent hinge line predictions across the core of our region but are more variable toward the edges of the model domain.

We introduced a hierarchical set of model ensembles constructed by ranking all 576 optimized configurations by post-alignment RMS and grouping them into four tiers: ELITE (RMS ≤ 0.49 mm/yr), GOOD (≤ 0.59 mm/yr), MEDIUM (≤0.79 mm/yr), and ALL (>0.80 mm/yr). These hybrid fields reveal a systematic progression, with the ELITE and GOOD ensembles capturing the GNSS-derived deformation shape with narrow uncertainty bands, while the MEDIUM and ALL ensembles exhibit progressively larger uncertainties that grow with ensemble size. The mean models of the ELITE and GOOD ensembles are nearly identical and provide the most stable uplift geometry and the smallest GPS-calibrated uncertainties, with representative values below 0.18 mm/yr (ELITE) and 0.26 mm/yr (GOOD) for Michigan, demonstrating that tightly constrained multi-model ensembles can outperform any individual GIA realization.