The nearshore total freeboard (NTF) of the Great Lakes serves as a sentinel of regional climate change, yet their role in capturing short-term climate signals has been largely overlooked. Here, we employ a Random Forestbased model to reconstruct NTF heights from 1950 to 2025, leveraging observations derived from groundbased GNSS Interferometric Reflectometry (GNSS-IR) signal of opportunity technique. Our reconstruction reveals a significant long-term thinning trend of -2.5 mm/decade since 1978. More importantly, we identify a dominant 4-year periodic oscillation in NTF heights, accounting for over 60% of the data variance. Strikingly, modulated by the synergistic enhancement of the East Pacific-North Pacific, North Atlantic Oscillation, and Pacific Decadal Oscillation patterns, the amplitude of this oscillation increased by 400% in 2011-2024 (12.6 +/- 1.7 mm) relative to that in 1950-1980 (2.5 +/- 0.7 mm). Our findings provide robust evidence that the winter lake ice-snow system is becoming increasingly vulnerable to climate change.

We extend our previous formulation of the gravitational potential of a homogeneous cuboid (right rectangular prism) by deriving uniform analytical expressions for the gravitational acceleration and the full gravitational tensor. The formulation is valid in the interior and exterior of the cuboid and also on its boundary. Classical rectangular–prism formulas express the potential and its derivatives through logarithmic and arctangent kernels whose arguments vanish at faces, edges, and vertices. Although these singularities are removable analytically, their direct numerical evaluation introduces instabilities near the boundaries of the body. We introduce a unified closed-form expressions based on oriented face–distance coordinates combined with analytically regularized sp_log and sp_arctan functions and an alternating vertex-sign summation. This representation produces continuous expressions for the gravitational potential, acceleration, and tensor that remain well defined at interior, boundary, and exterior field points without case distinctions or limiting procedures. Validation using acceleration-field analysis and an analytic Laplacian test confirms exact agreement with Poisson’s equation inside the cuboid and Laplace’s equation outside. The resulting expressions provide stable analytic first and second derivatives suitable for high-accuracy modeling of gravity and gravity-gradient measurements.

We construct an estimator for gravitational potential on an irregular closed surface as a finite spherical polynomial, based on a finite number of sample observations. This estimator is a natural discrete adaptation of an approximation for gravitational potential suggested by Sanso and Sideris. We consider the following discrete analogue of their conjecture: that the error in our estimator is a contraction operator, implying that the estimation process can be iterated to approximate gravitational potential to arbitrarily levels of accuracy. We then offer evidence for the conjecture through a series of numerical simulations and conclude with a simulation using a model of the Earth’s topography and real-world data.

We present a uniform closed-form boundary-continuous formulation for the gravitational potential, acceleration, and gravitational gradient tensor of a homogeneous polyhedron. Classical polyhedral gravitational models, though exact in form, become numerically unstable at faces, edges, and vertices and are conventionally restricted to exterior points. Our method eliminates these limitations through analytic regularization of the logarithmic and arctangent terms, ensuring continuity and stability across interior, boundary, and exterior domains. A dyadic tensor formulation provides compact expressions for accelerations and gravitational gradients, guaranteeing compliance with Poisson’s and Laplace’s equations. The implementation is fully vectorized and supports multi-threaded parallelism, delivering high computational efficiency and precision. Validation on convex, concave, and multiply connected geometries demonstrates strict physical consistency and numerical robustness. The proposed framework establishes a universal, high-precision computational standard for polyhedral gravitational modelling, applicable to planetary dynamics, spacecraft navigation, solid Earth geophysics, and physical geodesy.

We derive closed-form expressions for the Newtonian potential, gravitational acceleration, and gravitational gradient tensor generated by a homogeneous, planar, zero-thickness polygon embedded in three-dimensional space. The analysis is formulated for a general simple polygon with constant areal mass density σ and is implemented by decomposing the polygon into oriented triangular elements. For each triangle, the formulation expresses the potential and acceleration using edge log-factors and the signed solid angle, and provides an explicit dyadic representation for the tensor in terms of gradients of these quantities. For the important special case of a rectangle, compact corner-sum formulas provide efficient expressions for the potential, acceleration, and tensor. These formulas are numerically stable except at analytically singular boundary configurations, such as edges and corners of the ideal zero-thickness surface-mass model. Four independent benchmarks validate the theory: (i) a disk benchmark that demonstrates rapid convergence of the triangulated solution to the analytical axis field; (ii) the thin-cuboid limit, confirming that the rectangle model is the correct zero-thickness limit of a finite cuboid (or plate) with matching masses; (iii) reconstruction of the cuboid potential by stacking rectangles, confirming correct volume integration behavior; and (iv) agreement of rectangle and two-triangle decompositions at both acceleration and tensor levels. These results establish the correctness and consistency of the proposed polygonal gravitational primitives for geodetic and geophysical forward modeling.

Sea level is significantly modulated by various climate modes. Yet, the impacts of individual modes, their evolution, and links to present-day sea-level rise remain elusive. Here, we develop a data-driven, deep-learning-aided framework combining satellite altimetry and tide-gauge data to reconstruct seven decades (1952–2022) of global sea level. An Empirical Orthogonal Function decomposition of our reconstruction separates global sea-level rise from variability driven by the modulated annual cycle, El Niño–Southern Oscillation, Interdecadal Pacific Oscillation, and North Pacific Gyre Oscillation (NPGO), while quantifying each component. We reveal intensification across all quantified variability components. The intensity (standard deviation) of global-mean variability from summing these components nearly doubles, from 0.73 mm in 1952–1982 to 1.45 mm in 1992–2022. The corresponding global-mean sea level (GMSL) trend also nearly doubles, from 1.59 ± 0.26 mm yr⁻¹ to 3.10 ± 0.28 mm yr⁻¹. In moving-window estimates, both metrics are more pronounced since the late 1960s, and GMSL trends show relatively high correlations with the NPGO-associated variability intensity series (r ≈ 0.45–0.85; scale-dependent). We conclude that climate modes’ impacts are not static; they also strengthen alongside intensified sea-level rise, potentially reshaping the spatiotemporal characteristics of sea level under an increasingly warmer Earth.

The region of convergence of the spherical harmonic expansion is determined by the (generally complex) singularities of the gravitational potential. This complex analysis perspective is at the heart of recent rigorous results concerning the divergence properties of the spherical harmonic expansion. In this paper we build physical intuition for these general mathematical results using illustrative examples (some familiar and some new) of idealized planets for which the analysis is particularly explicit. This approach provides new methods to determine the region of convergence, without computing and analyzing expansion coefficients, and gives a novel geometric understanding of the divergence phenomenon. It also explains the fundamental origin of the numerical instabilities inherent to polyhedral models of planets. For the sake of clarity, we illustrate this new approach for the special case of axisymmetric planets of constant density, with explicit comparisons, but the key ideas do not rely on these restrictions.

We introduce a new analytical approach to computing the gravitational potential of a homogeneous tetrahedron, a fundamental building block in polyhedral modeling. This formalism is entirely free of numerical singularities. Unlike many existing solutions for polyhedral bodies, which are valid only outside the body, our expressions are uniformly applicable across the interior, boundary, and exterior of the tetrahedron. The method eliminates singularities by explicitly resolving geometric and analytical irregularities near edges, vertices, and face planes. We provide implementations in Python, MATLAB, and Julia, together with algorithmic pseudocode to facilitate adoption. Numerical validation via continuity tests and a Laplacian test confirms the stability, precision, and consistency of the formulation across all spatial regimes. This work broadens the class of analytically tractable gravitational primitives, complementing our prior solution for the cuboid. Since any constant-density polyhedron can be decomposed into tetrahedra, the presented formulation establishes a general pathway for modeling the gravitational fields of bodies with arbitrary geometry, providing robust tools for geodetic and geophysical modeling.

Rapidly rising sea level is one of the major adverse consequences of anthropogenic climate change. Sea level rise poses an existential threat to coastal populations, particularly for urban settlements with accelerating growth rates. Contemporary empirical sea level reconstructions have been used to conflate short-term (∼ 3 decades) satellite altimetry geocentric sea level data and long-term (50 years or longer) tide gauge records to better estimate reliable sea level rise towards multi-decadal to centennial time scales. However, adequate separations and quantifications of low-frequency climate patterns and sea level trends globally at regional scales remain elusive. Here, we propose a new sea level reconstruction framework that incorporates Empirical Orthogonal Function (EOF) into the contemporary Cyclostationary EOF with Reduced Space Optimal Interpolation (CSEOF-OI) algorithm to better reconstruct sea level fields. Using 225 selected long-term gap-filled tide gauge records with vertical land motion adjusted and satellite altimetry, our global reconstructed monthly sea level time series, January 1950–January 2022, exhibits distinct delineations between modeled climate patterns and sea level trends at 1°×1° regional scales. The separated sea level patterns include trends, modulated annual cycles, the El Niño Southern Oscillation (ENSO), and the Pacific Decadal Oscillation (PDO). The third principal component of the reconstructed sea level exhibits a Pearson correlation coefficient of 0.87 with the Niño 3.4 ENSO index, and the fourth principal component correlates at 0.75 with the PDO index, indicating good agreement. The global mean sea level trend, accounting for the predominant climate periodicities, is 1.9 ± 0.2 mm yr−1 (95 % confidence), and the estimate during the satellite altimetry era (January 1993–December 2021) is 3.2 ± 0.3 mm yr−1 (95 % confidence). Compared with previous studies, we conclude that our 72-year sea-level reconstruction allows us to better separate the ENSO and PDO climate patterns, as well as the sea level they induced. Finally, we show that the short-term (5-year) rates of ENSO and PDO patterns significantly affect sea level both on a global and regional scale, altering global mean sea level trends by up to 1.1 ± 0.5 mm yr−1 (January 2011–January 2016). Over the past seven decades, the climate patterns exerted a minor impact on sea level trends, but substantially modulated apparent regional sea level accelerations, particularly in the western Pacific (e.g., 0.09 ± 0.05 mm yr−2 at the Kuroshio Current), and in the east and central equatorial Pacific Ocean (e.g., −0.04 ± 0.03 mm yr−2 near Costa Rica). The reconstructed sea level and analysis results datasets are available at https://doi.org/10.5281/zenodo.15288816 (Wang, 2025).

The Central and South-Central Andes form a “two-sided” mountain belt bounded by distinct zones of convergence in the western forearc and eastern foreland flanks. Previous geodetic studies of interseismic deformation in the Bolivian Subandes and the Argentine Precordillera found that the forearc to foreland velocity field decayed too slowly to be explained purely by elastic shortening driven by locking of the Nazca megathrust. The velocity field is more precisely explained if elastic deformation is augmented by eastward displacement of the entire Andes. Here, we extend the earlier interpretation of interseismic motion and argue that foreland décollements can participate in the co- and postseismic phases of the earthquake deformation cycle associated with the Nazca megathrust. These findings have direct implications in estimating recurrence interval, slip rate, and probabilistic seismic hazard analysis on both sides of the orogen.

Rapid melting of the Greenland Ice Sheet (GrIS) in response to global warming has been a major contributor to global sea level rise in the last 20 years. The ability of the Gravity Recovery and Climate Experiment (GRACE) to estimate GrIS mass changes is limited by its coarse (similar to 330 x 330 km(2)) spatial resolution. The Greenland Geodetic Network (GNET) senses the solid Earth's elastic responses to changing ice mass loads, as well as glacial isostatic adjustment (GIA). The GNET stations are sensitive to local ice mass changes at the scale of tens of kilometers, but have poor spatial coverage compared to GRACE, since all bedrock Global Navigation Satellite System (GNSS) stations are located near the margins of the GrIS. GRACE gravity observations and GNSS measurements of crustal displacement provide complementary constraints on GrIS mass changes. Here, we exploit this complementarity, by developing a joint inversion method that combines the norms of gradients and makes judicious use of the Lcurve to estimate GrIS mass changes at 0.25 degrees- grids. We modify the Laplacian operator, commonly used in previous studies, to make it suitable for the irregular inversion area and to avoid unrealistic inversion results at the land-sea boundary in Greenland. We have adopted a new weight allocation strategy to ensure that GRACE and GNSS data make similar contributions to the joint inversion results, avoiding the loss of information contained in GNSS due to the difference in spatial coverage between the two types of data. The joint inversion results are compared to satellite altimetry-derived GrIS mass changes and two GNET verification sites not involved in the inversion. This joint inversion method most strongly improves the spatial resolution of ice mass change estimates in low-altitude areas of Greenland. We recovered the melting signal of the GrIS leaking into the nonice-covered land, and joint inversion indicates during January 2008 to December 2020 an ice mass trend (-254.0 Gt/yr) which is slightly slower than that inferred using GRACE mascon methods (-261.3 Gt/yr), but the annual amplitude of ice mass change (152.6 Gt) is significantly higher than GRACE (135.0 Gt). We identified areas with significant changes in ice mass that were not resolved by GRACE, and found (not surprisingly) that mass fluctuations were greater at outlet glacier locations than in adjacent areas of the ice sheet. Ice mass changes inferred from vertical land motion are sensitive to GIA corrections, and the difference in the rate of ice mass loss evaluated using different GIA models can reach around 20 Gt/yr, similar to GRACE-inferred mass estimates. This study successfully applies inverting GNSS and GRACE data for ice mass change at the edge of Greenland.

Abstract Greenland's bedrock responds to ongoing ice loss with an elastic vertical land motion (VLM) that is measured by Greenland's Global Navigation Satellite System (GNSS) Network (GNET). The measured VLM also contains other contributions, including the long‐term viscoelastic response of the Earth to the deglaciation of the last glacial period. Greenland's ice sheet (GrIS) produces the most significant contribution to the total VLM. The contribution of peripheral glaciers (PGs) from both Greenland (GrPGs) and Arctic Canada (CanPGs) has not carefully been accounted for in previous GNSS analyses. This is a significant concern, since GNET stations are often closer to PGs than to the ice sheet. We find that, PGs produce significant elastic rebound, especially in North and East Greenland. Across these regions, the PGs produce up to 32% of the elastic rebound. For a few stations in the North, the VLM from PGs is larger than that due to the GrIS.

SUMMARY In this paper, we propose a fully coupled two-phase poroelastic deformation theory for a spherically layered and self-gravitating Earth. The earth model consists of a solid inner core, a fluid outer core and poroelastic mantle and crust. The boundary-value problems are posed in the Laplace-transformed domain using the spherical system of vector functions, and analytical solutions are obtained in each layer using the dual-variable and position matrix method. As an application, the surface loading problem is considered. The undrained and drained limits are discussed with the corresponding governing equations, which are different from the conventional ones where body force is ignored. Numerical examples show that the poroelastic effect can be significant since the difference between the results for undrained and drained limits are large with obvious time-dependent deformation. This newly derived theory for the coupled boundary-value problem will have broad applications, such as displacement and gravity change due to groundwater depletion, poroelastic rebound after an earthquake and risk evaluation of earthquakes induced by water injection.

This study investigates the impact of observation session duration and station latitude on the precision of Precise Point Positioning (PPP). Global Positioning System-only Receiver Independent Exchange files from 516 continuous stations spanning latitudes from 90°N to 90°S across the Americas (30°–130°W) were binned into sessions of 1, 2, 3, 4, 6, and 12 h and processed using the GPSPACE PPP software. These sub-daily solutions, along with the original 24 h ones, were compared against a reference coordinate derived from an extended linear trajectory model for each station. Analysis of the results, considering both accuracy and precision, was conducted across the entire latitude range of the study. We found that the precision of a PPP solution approximately follows a power-law relationship with observation duration. Parameters for these power-law relationships were determined for all latitude ranges that allow users to predict a result’s uncertainty as a function of session length. Findings indicate that longer observation sessions lead to reduced positioning errors, with vertical scatter decreasing with increasing latitude. Since all these stations are characterized by good to excellent sky view, our power-law rule-of-thumb provides a lower bound on the occupation time needed to achieve target positioning precision at locations with poorer sky visibility.

The Brillouin sphere is defined as the smallest sphere, centered at the origin of the geocentric coordinate system, that incorporates all the condensed matter composing the planet. The Brillouin sphere touches the Earth at a single point, and the radial line that begins at the origin and passes through that point is called the singular radial line. For about 60 years there has been a persistent anxiety about whether or not a spherical harmonic (SH) expansion of the external gravitational potential, V, will converge beneath the Brillouin sphere. Recently, it was proven that the probability of such convergence is zero. One of these proofs provided an asymptotic relation, called Costin's formula, for the upper bound, E-N, on the absolute value of the prediction error, e(N), of a SH series model, V-N(theta,lambda,r), truncated at some maximum degree, N = n(max). When the SH series is restricted to (or projected onto) a particular radial line, it reduces to a Taylor series (TS) in 1/r. Costin's formula is E-N similar or equal to BN-b(R/r)(N), where R is the radius of the Brillouin sphere. This formula depends on two positive parameters: b, which controls the decay of error amplitude as a function of N when r is fixed, and a scale factor B. We show here that Costin's formula derives from a similar asymptotic relation for the upper bound, A(n) on the absolute value of the TS coefficients, a(n), for the same radial line. This formula, A(n) similar or equal to Kn(-k), depends on degree, n, and two positive parameters, k and K, that are analogous to b and B. We use synthetic planets, for which we can compute the potential, V, and also the radial component of gravitational acceleration, g(r) = partial derivative V/partial derivative r, to hundreds of significant digits, to validate both of these asymptotic formulas. Let superscript V refer to asymptotic parameters associated with the coefficients and prediction errors for gravitational potential, and superscript g to the coefficients and predictions errors associated with g(r). For polyhedral planets of uniform density we show that b(V) = k(V) = 7/2 and b(g) = k(g) = 5/2 almost everywhere. We show that the frequency of oscillation (around zero) of the TS coefficients and the series prediction errors, for a given radial line, is controlled by the geocentric angle, alpha, between that radial line and the singular radial line. We also derive useful identities connecting K-V, B-V, K-g, and B-g. These identities are expressed in terms of quotients of the various scale factors. The only other quantities involved in these identities are alpha and R. The phenomenology of 'series divergence' and prediction error (when r < R) can be described as a function of the truncation degree, N, or the depth, d, beneath the Brillouin sphere. For a fixed r <= R, as N increases from very low values, the upper error bound E-N shrinks until it reaches its minimum (best) value when N reaches some particular or optimum value, N-opt. When N > N-opt, prediction error grows as N continues to increase. Eventually, when N >> N-opt, prediction errors increase exponentially with rising N. If we fix the value of N and allow R/r to vary, then we find that prediction error in free space beneath the Brillouin sphere increases exponentially with depth, d, beneath the Brillouin sphere. Because b(g) = b(V) - 1 everywhere, divergence driven prediction error intensifies more rapidly for g(r) than for V, both in terms of its dependence on N and d. If we fix both N and d, and focus on the 'lateral' variations in prediction error, we observe that divergence and prediction error tend to increase (as does B) as we approach high-amplitude topography.

Unmodeled displacements in GNSS times series, induced by instrumental artifacts or geophysical events, create significant biases in station trajectory parameters that can propagate into the reference frame itself. While non-tectonic ‘jumps’, such as equipment changes, affect only a specific GNSS station, seismically-induced displacements can affect large numbers of sites, severely threatening the frame’s stability. Manually reviewing individual GNSS time series for such effects is highly impractical because there can be thousands of GNSS stations in a frame, and the total number of earthquakes Mw ≥ 6.0 since GPS became fully operational is + 5100. To avoid this time-consuming task, automated methods rely on empirical power-law functions to determine which earthquake-station pairs require coseismic displacement parameters. Still, ‘conservative’ power-law functions tend to add coseismic offsets to stations that do not need them, which can also threaten the stability of the frame. In this work, we present an empirical formulation that was obtained using 809 global seismic events to fit power-law parameters that do not overestimate the region of influence of earthquakes. Our method is based on a two-level selection process: level 1 is isotropic and only considers the epicentral distance between the stations and the earthquake, and level 2 uses the geophysical parameters of the earthquake to predict a ‘tighter’ displacement pattern to select which stations require coseismic trajectory parameters. We applied our level 2 method to a database of 4700 event-station pairs and showed that it removed 55

The Greenland ice sheet (GrIS) is at present the largest single contributor to global-mass-induced sea-level rise, primarily because of Arctic amplification on an increasingly warmer Earth1-5. However, the processes of englacial water accumulation, storage and ultimate release remain poorly constrained. Here we show that a noticeable amount of the summertime meltwater mass is temporally buffered along the entire GrIS periphery, peaking in July and gradually reducing thereafter. Our results arise from quantifying the spatiotemporal behaviour of the total mass of water leaving the GrIS by analysing bedrock elastic deformation measured by Global Navigation Satellite System (GNSS) stations. The buffered meltwater causes a subsidence of the bedrock close to GNSS stations of at most approximately 5 mm during the melt season. Regionally, the duration of meltwater storage ranges from 4.5 weeks in the southeast to 9 weeks elsewhere. We also show that the meltwater runoff modelled from regional climate models may contain systematic errors, requiring further scaling of up to about 20% for the warmest years. These results reveal a high potential for GNSS data to constrain poorly known hydrological processes in Greenland, forming the basis for improved projections of future GrIS melt behaviour and the associated sea-level rise6.

Like many geophysical observations, relative gravity (RG) measurements are affected by random errors, systematic errors, and occasional blunders. When RG measurements are used to build large gravity networks in remote areas under adverse environmental or logistical conditions (such as extreme temperatures, heavy precipitation, rugged terrain, difficult or dangerous roads, and high altitudes), it is more likely for significant errors to occur and accumulate. Therefore, obtaining accurate gravity estimates at regional gravity networks largely depends on defensive data collection protocols and robust adjustment techniques. In this work, we present a measurement field protocol based on highly redundant observation patterns, and a two-step least squares adjustment scheme implemented as a MATLAB package. This software helps us identify blunders, mitigates the impact of random errors, and downweights or removes outlier observations. The methodology also guarantees that adjusted gravity values have well-constrained standard error estimates. We illustrate the capabilities of our approach through the case study of the Bolivian gravity network, where we determined the acceleration due to gravity at 2548 stations that spread over difficult and sometimes extreme environments, with a typical level of uncertainty of 0.10–0.15 mGal.

ANET-POLENET (Antarctic Network of the Polar Earth Observing Network) bedrock GNSS sites in the Ross Sea region of Antarctica surround an LGM load center in the Siple region of the Ross Embayment and record crustal motion due to GIA. Rather than a radial pattern of horizontal motion away from the former load, we instead observe three primary patterns of deformation; 1) motions are reversed towards the load in the southern region of the Transantarctic Mountains (TAM), 2) motions are radially away from the load in the Marie Byrd Land (MBL) region, and 3) an overall gradient in motion is present, with magnitudes progressively increasing from East to West Antarctica. We investigate the effects of alternative Earth model and ice loading scenarios, with the goal of understanding these distinct patterns of horizontal bedrock motion and their drivers. Using GIA models with a range of 1D Earth models, alternative ice loading scenarios for the Wilkes Subglacial Basin (LGM time scale) and the Siple Coast (centennial and millennial time scales) are explored. We find that no 1D model, regardless of the Earth model and ice loading scenario used, reproduces all three distinct patterns of observed motion at the same time. For select ice loading scenarios we also examine the influence of more complex rheology by invoking a boundary in Earth properties beneath the Transantarctic Mountains. This approach accounts for the strong lateral gradient in Earth properties across the continent by effectively separating East and West Antarctica into two different Earth model profiles. Some of our GIA models utilizing 3D Earth structure reproduce predicted motions that match all three observed patterns of deformation, and we find that a multiple order magnitude of change in upper mantle viscosity between East and West Antarctica is required to fit the observations.

We present a new, physically motivated triaxial reference ellipsoid for the Earth. It is an equipotential surface in the gravity field and closely approximates the geoid, akin to the conventional reference ellipsoid of revolution. According to Burša and Fialová (Studia Geophysica et Geodaetica 37(1):1–13, 1993), the triaxial reference ellipsoid is uniquely, but not exclusively, specified by the body’s total mass, the dynamic form factors of polar and equatorial flattening, the longitude of the equatorial major axis, the rotation rate, and the designated surface potential. We model the gravity field using triaxial ellipsoidal harmonics. While they are rarely considered practical for near-spherical planets, we leverage an intrinsic property that ellipsoidal harmonics yield an exact expression for the constant potential on a triaxial ellipsoid. A practical procedure is proposed to solve for the ellipsoidal parameters that converge iteratively to fulfill the exact condition of equipotentiality. We present the solution for the Earth Gravitational Model 2008.