Study region: Lhasa Superconducting Gravimeter (SG) Observatory, Lhasa River alluvial plain Study focus: The Lhasa SG Observatory is the only continuously operating SG station on the Tibetan Plateau. Surrounded by thick Quaternary sediments, this site provides a critical window into the interactions between tectonic processes and hydrological mass redistribution. Precisely isolating gravity interference caused by local groundwater storage changes is essential for detecting subtle geodynamic signals, such as crustal thickening. We integrated high-precision SG observations, meteorological forcing, and in-situ groundwater levels (2010-2020) into a 1D physically-based Richards equation framework. We reconstructed the spatiotemporal evolution of soil moisture within the 3-meter unsaturated zone to accurately quantify the gravity effects induced by localized hydrological dynamics. New hydrological insights for the region: The physical model's reconstruction exhibits strong consistency with SG residuals at an hourly scale (cross-correlation coefficient: 0.62), significantly outperforming global hydrological products like ERA5 (0.18) and GLDAS (0.55). Groundwater-induced gravity fluctuations reach an amplitude of 10.62 mu Gal (1 mu Gal = 1 & sdot;10-8 m s-2), sufficient to mask contemporaneous tectonic signatures. Crucially, long-term regression identifies a persistent gravity decline of approximately-0.27 +/- 0.002 mu Gal & sdot;a-1 driven by continuous groundwater depletion. This trend accounts for nearly 14 %-40 % of the observed absolute gravity variation rate. Neglecting station-scale hydrological corrections can thus lead to substantial misjudgments of crustal thickening rates and Moho subsidence magnitudes on the Tibetan Plateau.

The detection and analysis of the Earth’s free oscillations in the form of normal modes contribute to our understanding of the planet’s structure, dynamics, and seismicity. However, false and dismissed recognitions of the focused normal mode signals prevent the detailed exploration of the Earth’s interior due to the complex patterns of the normal mode spectrum. Here, we develop a deep-learning-based neural network, named ModeNet, which is capable of precisely and efficiently selecting the frequency windows to cover the target normal modes on noisy spectra. ModeNet achieves a remarkable precision rate in the discrimination between normal modes and noises of 98.1

Identifying precursors of large earthquakes is critical for minimizing the losses of life and property. Recently, Bletery and Nocquet (2023) captured a similar to 2-h-long exponential acceleration of slip using the high-rate (5-min) Global Navigation Satellite System (GNSS) time series from the 48 hr before the 2011 M-W 9.0 Tohoku-oki earthquake, which was obtained by simply concatenating daily kinematic results together. Here, we apply their method to sum the horizontal displacements of 24 high-rate GNSS stations in the direction predicted by fault slip at the hypocenter of the 2023 M-W 7.8 Kahramanmaras earthquake to characterize its precursory phase. Results demonstrate a several-hour accelerating exponential slip before the mainshock. However, considering that single-day processing would lead to discontinuities at the day boundary, we process the multi-day GNSS data in continuous mode, repeat the experiment, and find that the observed acceleration-like signals vanish. Our work shows that inadequate data processing may lead to the detection of false precursory signals, highlighting the need to develop robust processing techniques to identify reliable precursory signals before large earthquakes.

The seafloor topography (ST) is essential for oceanographic research and geophysical applications, and its modeling relies on satellite altimetry-derived gravity data. The ST inversion commonly applies linear regression method that account primarily for shipborne bathymetry errors, while often neglecting the potential uncertainties in altimetry-derived gravity data. To address this question, we draw on the idea of total least squares, which has been widely applied in numerical analysis. In this study, we propose an improved linear regression framework for ST inversion: we use high-precision multibeam shipborne bathymetry data (MSB1) as a constraint to iteratively determine the optimal weight ratio between bathymetry and gravity data, and subsequently construct initial weight matrices; then the robust weighted total least squares (RWTLS) method is applied to estimate regression parameters, and the corresponding ST model is constructed. A case study is conducted in a region of the South China Sea (113 degrees-118 degrees E, 13 degrees-19 degrees N). Comparisons with the check data (MSB2) and existing models (topo_27.1, ETOPO2022, and SDUST2023BCO) indicate that the overall accuracy of the constructed ST model is approximately 150 m. And compared with the traditional method, the application of RWTLS improves ST inversion accuracy, with power spectral density analysis further revealing a significant enhancement in topographic energy within the 12-30 km wavelengths. Moreover, we find that the improvement in ST accuracy is mainly concentrated in offshore waters; in these areas, altimetry-derived gravity data tend to have greater uncertainty, nevertheless, the improved method effectively enhances the reliability of ST modeling by jointly accounting for the errors in both bathymetry and gravity data. This study proposes a robust inversion framework that considers multisource data uncertainties, offering a novel strategy for robust linear ST modeling.

Why is Earth, among the eight planets in our solar system, the only habitable one? Over the 4.6-billion-year evolution of the solar system, why did Mars and Venus evolve so differently? Where did life originate, and how will Earth evolve in the future? These questions are not only central to planetary science in the 21st century but are also deeply connected to humanity’s fundamental understanding of its own existence and planetary habitability. Addressing these grand scientific challenges demands a systematic, multidimensional research approach. In the temporal dimension, we need to trace the early formation and evolution of terrestrial planets; in the spatial dimension, we need to analyze the layered structure of planets and the coupling between these layers; from a comparative perspective, we also need to explore the atmospheric characteristics of exoplanets and the influence of their host stars on habitability. Supported by the Chinese Academy of Sciences’ Strategic Priority Research Program on “Formation, Evolution, and Habitability of Terrestrial Planets,” we have taken planetary habitability as the main research theme and conducted systematic, in-depth studies on terrestrial planets by integrating multiple approaches, including extraterrestrial sample analysis, deep-space exploration data processing, and numerical and experimental simulation. This paper comprehensively summarizes the significant advancements made by the project over the past five years. It covers topics ranging from the early processes and environmental evolution of terrestrial planets to open planetary systems linked to the external space environment, and from Earth’s Moon to exoplanets. Key achievements include the first confirmation of a solid inner core on Mars, revealing its core-mantle differentiation under high-pressure and high-temperature conditions, and the discovery that the youngest lunar basalts originated from a non-KREEP, volatile-poor mantle source region, challenging the long-held traditional hypothesis that “volatile-rich material drives late-stage volcanism”—a textbook-level achievement. While summarizing the latest research advances, this paper also looks toward future directions. Significantly improving the capability for multi-layered, multi-parameter detection of planets, especially global planetary survey capabilities, and vigorously developing related techniques and research methods, particularly the application of cutting-edge technologies such as quantum technology and artificial intelligence, will be the key to achieving further major breakthroughs in planetary science. By sharing these research findings and insights, we aim to inspire more young people to pursue careers in planetary science—a field full of opportunities and challenges—thereby promoting the sustainable development of planetary science in China and over the world.

Earth tidal perturbations affecting laser-ranged satellites are critical for refining satellite orbital dynamics modeling, and their accurate computation represents a prerequisite for high-precision fundamental physical effects and geodetic investigations based on satellite orbit analysis. This study focuses on the tidal perturbations induced by the asymmetric responses of LARES 2 and LAGEOS on their orbital nodes and inclinations. Perturbations induced by a total of 402 (392 2nd and 10 3rd-order) earth tide constituents on the two satellites were calculated, based on Kaula's orbital perturbation theory and Lagrange's planetary equations for satellites, considering the frequency dependence of Love numbers. The asymmetric characteristics of tidal perturbations between the two satellites were quantitatively analyzed. The minimum resolutions of orbital inclinations and nodes, used as screening thresholds for significant constituents, were derived from the RMS of overlapping orbit differences using orbital geometry and error propagation law. With these thresholds, 83 significant constituents were identified from the 402. The cumulative effect of the 319 minor constituents was further evaluated, and it was found that their total impact, from coherent superposition, noticeably exceeds the thresholds, thus becoming non-negligible. The results of this study provide accurate tidal perturbation parameters for LARES 2 and LAGEOS, and offer methodological references for the screening of Earth tide constituents in high-precision satellite orbital dynamics research, laying a foundation for subsequent studies on inverting geophysical parameters from satellite orbits and verifying fundamental physical effects, particularly the relativistic Lense-Thirring effect.

It remains unclear how the spatial variability in interseismic coupling contributes to assessing future seismic rupture patterns, such as the slip magnitude and rupture speed. The 2025 Mw 7.8 Myanmar earthquake provides a rare opportunity to address this question. Herein, we systematically integrate the viscoelastic deformation model, regional slowness-enhanced backprojection, finite fault inversion, and seismic waveform analysis to characterize the interseismic fault coupling and coseismic rupture history associated with the Myanmar event in detail. Our results reveal that the Myanmar event ruptured multiple segments of the Sagaing fault, involving two distinct phases: an initial bilateral rupture followed by a persistent supershear rupture. Notably, this stable supershear rupture spatially coincides with the seismic gap. The interseismic coupling model indicates that the gap is characterized by high stress buildup and low dissipated-to-potential energy ratio Gc/G0 (i.e., fracture energy over static energy release rate), indicating a high readiness for supershear rupture. Moreover, the coseismic high-slip regions correlate closely with the strongly coupled zones. The accumulated slip deficit on each fault section since the last major earthquake scales with its respective average coseismic slip, revealing loosely slip-predictable behavior. Our study thus reveals the possibility of assessing the coseismic slip and supershear rupture potential of future earthquakes based on the interseismic coupling and rupture histories of fault systems.

Whether and when earthquakes of different sizes can be distinguished early in their rupture process is critical to improving earthquake early warning (EEW) systems. To address this question, we develop a transformer-based framework called source time function Magnitude Network (STF-MgNet), which leverages earthquake STFs with M w 5.5 to investigate whether the initial stages of the rupture process can predict the earthquake's final magnitudes. The proposed STF-MgNet utilizes transformer blocks in conjunction with a U-Net architecture to effectively manage long-range dependencies and capture essential information in data sets, boosting its performance. Training on 2,126 global M w ≥ 5.5 earthquakes reveals three distinct nucleation-phase diagnostic regimes: (a) For 5.5 M w ≤ 7 events, subsecond analysis of the STF achieves 80.0% magnitude estimation accuracy; (b) For 7 < M w ≤ 8 earthquakes, 3–5 s of observations capture key features of STFs, yielding 86.9% accuracy; (c) For M w > 8 events, approximately 5 s of monitoring resolves the interplay between determinism and stochasticity, maintaining 83.7% accuracy despite greater complexity. Thus, no matter how an earthquake begins, earthquakes of different sizes can be distinguished at some point. Moreover, the earthquake magnitude prediction accuracy of STF-MgNet improves with increasing input sequences of STFs. This study supports the hypothesis that rupture onset at nucleation differs for earthquakes of different final magnitudes.

Accurate determination of the parameters of the Earth’s free core nutation (FCN) provides insights into the core-mantle coupling mechanism and helps refine modeling of celestial pole offsets. Recent studies suggest that the FCN period exhibits time-varying characteristics, potentially related to core-mantle interactions, but precise detection remains challenging. Traditionally, FCN periods estimated from superconducting gravimeter (SG) observations, based on resonance phenomena in diurnal tidal waves, rely mostly on the Ψ_1 wave. However, the low signal-to-noise ratio (SNR) of the Ψ_1 wave leads to significant uncertainties in the precise detection of temporal variations in the FCN period. In this study, we propose a method that relies solely on K_1 wave with a higher SNR, based on the sensitivity of the K_1 wave to the time-varying FCN parameter. The simulation results show that the new method can overcome the limitations of using the Ψ_1 wave and can effectively capture the variations in the FCN period within a few days under current SG observational precision. The method is applied to observations from nine SG stations in the International Geodynamic and Earth Tide Service (IGETS) network. The results indicate that the proposed method can significantly improve the detection of temporal variation of the FCN period and can obtain results close to those of the very long baseline interferometry (VLBI) technique. This method can enhance the effective utilization of SG observations in detecting weak signals related to the dynamics of the Earth’s core.

Geophysical inversion plays a pivotal role in understanding the Earth's internal structure. Recently generative neural networks (GNNs), such as normalizing flows models (NFMs), have gained popularity for solving Bayesian inversion problems. However, the posterior probability density functions (PDFs) obtained by amortized GNN‐based methods often deviates from the target distribution. This discrepancy arises because traditional amortized methods use joint PDFs as the objective in loss functions, rather than the conditional PDFs of the observed data. To address this, we propose the Iterative Normalizing Flows Model (INFM), a novel approach that mitigates loss function bias by progressively narrowing the prior distribution's support set in each iteration, while ensuring that the posterior distribution accurately converges to the target distribution. Our experiment, validated on high‐dimensional Bayesian inversion tasks, shows that INFM significantly enhances inversion accuracy without increasing network complexity or computational cost. When applied to the Earth's 1‐D structure model inversion, our method revealed key insights, such as a lower core density compared to the Preliminary Reference Earth Model (PREM) model and the presence of anisotropy in both the mantle and core, consistent with previous studies. These findings suggest that the INFM method offer high computational efficiency and accuracy, making it well‐suited for large‐scale geophysical inversion problems.

Interplay between seismic and aseismic slip could shed light on the frictional properties and seismic potential of faults. The well-recorded 2023 Kahramanmaraş earthquake doublet provides an excellent opportunity to understand their partitioning on strike-slip faults. Here, we utilize InSAR and strong motion data to derive the coseismic rupture during the doublet, ~4-month postseismic afterslip, and slip distributions of two Mw>6.0 aftershocks. Our results show that afterslip appears to be complementary to coseismic slip and aftershocks, accounting for ~11.3% of the coseismic moment. Aftershocks mainly fall within the regions of positive Coulomb stresses caused by afterslip and follow a temporal decay similar to that of afterslip, indicating that aftershock production is the failure of small asperities loaded by the afterslip. The early postseismic afterslip is released ~93.7% aseismically and ~6.3% seismically by aftershocks. Our modeling results thus depict a complex fault system with highly variable slip patterns and stresses. The study examines interplay between seismic and aseismic slip associated with the 2023 Kahramanmaraş earthquake doublet, constrained by InSAR and strong motion data, shedding light on the frictional properties and seismic potential of faults.

Accurate theoretical simulations of Earth deformation play a crucial role in understanding the Earth's internal structure, interpreting observational data, and understanding its deformation processes. This paper develops a novel error control algorithm based on adaptive Runge-Kutta theory that can handle multiple initial solution conditions. The application of the algorithm focuses on solving the deformation problem caused by surface mass loading and the normal modes of the Earth. By implementing matrix transformations to achieve global truncation error control, the proposed method achieves precise numerical solutions, and the results show significant improvements in accuracy compared to traditional methods. The results show that when calculating the 1000 degree load Love numbers for a homogeneous model, the adaptive Runge-Kutta methods RK8(7) and RK5(4) achieve two orders of magnitude higher accuracy than the fourth-order Runge-Kutta method (RK4), with computational speeds increased by 20 times and 5 times, respectively. For higher degree load Love numbers of the PREM model, the convergence of the RK4 method is poor, resulting in absolute errors of up to 1x10-4. Conversely, the RK8(7) method shows improved convergence with increasing degree and maintains relative stability. In addition, by setting the new initial solution and integration position, the RK method can calculate up to 10 million degrees for the load Love numbers, with an accuracy of 1x10(-7). When the tolerance of 1x10(-6) is set, the absolute error of the eigenperiod calculated by RK8(7) is generally about 1x10-6 second compared with the analytical solution. In addition, the root-finding speed for a degree of 250 is 17.42 times faster than that of RK4. Finally, cross-validation with the Mineos software confirms that the eigenperiods of the PREM model can achieve an accuracy of at least 1x10(-3) second. Specifically, the absolute error of the Slichter mode eigenperiod compared to other numerical methods and observed periods is 0.0055 hour and 0.0156 hour, respectively, further demonstrating the accuracy and reliability of this method.

The Surface Water and Ocean Topography (SWOT) wide-swath altimetry satellite was launched in December 2022. The performance of novel wide-swath altimetry in seafloor topography modeling needs to be evaluated. This study utilized 15 cycles of SWOT Level-3 product to construct seafloor topography model of the South China Sea by linear regression analysis. The root mean square error of the difference between the model and shipborne bathymetry at checkpoints is about 120 m, which is 20 m better than topo_27.1 and DTU18BAT, and 40 m better than ETOPO1. First, the effects of the shipborne bathymetry at control points and priori bathymetry model in different topography-gravity scaling factor estimation strategies (A: using robust least squares (RBLSQ) to estimate regional scaling factor; B: using ratio method to calculate scaling factors at control points; C: using the moving window method and RBLSQ to obtain scaling factor grids.) on SWOT seafloor topography modeling are explored. We find that the control point number barely affects strategy A/C but significantly affects strategy B, while the priori bathymetry model mainly affects strategy C. Then, the three strategies are applied to the traditional radar altimetry gravity anomaly, and the results are compared with the SWOT-derived seafloor topography. The results show that incorporating SWOT data can improve the accuracy of seafloor topography estimation by about 7 m, and improve the power spectral density in the wavelength range about 10∼20 km, which can help to reveal more detailed topography information.

Since 2000, eastern Taiwan has experienced 38 Mw >= 5:5 earthquakes, leaving three seismic gaps along the Longitudinal Valley fault (LVF). In April 2024, the Mw 7.3 Hualien earthquake occurred near the LVF. Herein, we first apply comprehensive geodetic data including Interferometric Synthetic Aperture Radar and Global Navigation Satellite System to estimate two potential fault geometries and invert for the coseismic slip. Our results suggest that a transpressive WNW-dipping low-angle fault related to the Central range fault is responsible for the Mw 7.3 Hualien earthquake. We then perform the Coulomb stress analysis to probe earthquake interaction in eastern Taiwan. The increased stress of similar to 2.6 bars due to the preceding major earthquakes at the hypo- center of the 2024 event significantly pushes this fault toward failure. Moreover, the conjugate LVF and the Milun fault are activated, and some aftershocks are promoted here. Finally, we note that the Coulomb stress changes from historical earthquakes and the 2024 Hualien earthquake exert positive stress on the seismic gaps in the northern LVF, potentially influencing future ruptures.

It is commonly believed that the atmosphere is decoupled from the solid Earth. Thus, it is difficult for the seismic wave energy inside the Earth to propagate into the atmosphere, and atmospheric pressure wave signals excited by earthquakes are unlikely to exist in atmospheric observations. An increasing number of studies have shown that earthquakes, volcanoes, and tsunamis can perturb the Earth's atmosphere due to various coupling effects. However, the observations mainly focus on acoustic waves with periods of less than 10 min and inertial gravity waves with periods of greater than 1 h. There are almost no clear observations of gravity waves that coincide with observations of low-frequency signals of the Earth's free oscillation frequency band within 1 h. This paper investigates atmospheric gravity wave signals within 1 h of surface-atmosphere observations using the periodogram method based on seismometer and microbarometer observations from the global seismic network before and after the July 29, 2021 MW8.2 Alaska earthquake in the United States. The numerical results show that the atmospheric gravity wave signals with frequencies similar to those of the Earth's free oscillations 0S2 and 0T2 can be detected in the microbarometer observations. The results confirm the existence of atmospheric gravity waves, indicating that the atmosphere and the solid Earth are not decoupled within this frequency band and that seismic wave energy excited by earthquakes can propagate from the interior of the Earth to the atmosphere and enhance the atmospheric gravity wave signals within 1 h.

On December 18, 2023, the Ms 6.2 Jishishan earthquake occurred in the northeastern region of the Tibetan Plateau, causing heavy casualties and property damage in Gansu and Qinghai provinces. In this study, we integrate space imaging geodesy, finite fault inversion, and back-projection method to decipher its rupture property, including fault geometry, coseismic slip distribution, rupture direction, and propagation speed. The results reveal that the seismogenic fault dips to the southwest at an angle of 29°. The major slip asperity is dominated by reverse slip and is concentrated within a depth range of 7–16 km, which explains the significant uplift near the epicenter observed by both the Sentinel-1 ascending and descending InSAR data. Moreover, the teleseismic array waveforms indicate a northwest propagating rupture with an overall slow rupture velocity of ∼1.91 km/s (AK array) or 1.01 km/s (AU array).

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.

Diurnal tidal oscillations in the coupled atmosphere–ocean system generate important contributions to the Earth’s free core nutation (FCN) and annual and sub-annual components of forced nutation in the celestial pole offsets. The determination of FCN parameters cannot avoid the influence of geophysical fluid excitation neither with the direct analysis of FCN signal (direct approaches) nor with the resonance analysis of forced nutation (resonance approaches). There is a significant difference in the FCN parameters obtained with resonance and direct approaches from celestial pole offsets observed through very long baseline interferometry (VLBI). The source of the difference between the two lacks quantitative analysis, which causes difficulties in interpreting the validity of the derived FCN parameters. Using both approaches, we conducted a simulation of celestial pole offsets to quantitatively demonstrate how geophysical fluid excitation affects the determination of FCN parameters from VLBI observations. Using the same excitation source, the FCN period obtained by the direct approach deviated from the set value (430.21 d) by more than 10 d, while the FCN period obtained by the resonance approach showed no deviation from the set value by more than 1 d. The results indicate that the resonance approach more accurately reflects the intrinsic period of the FCN. The impact of atmospheric and oceanic contributions on the determination of the FCN period with the resonance approach was within 2 d. Numerical simulation shows that discrepancies in FCN parameters caused by geophysical excitation were nonnegligible in constructing accurate FCN models.

The Earth's Free Core Nutation (FCN) causes Earth tides and forced nutation with frequencies close to the FCN that exhibit resonance effects. High-precision superconducting gravimeter (SG) and very long baseline interferometry (VLBI) provide good observation techniques for detecting the FCN parameters. However, some choices in data processing and solution procedures increase the uncertainty of the FCN parameters. In this study, we analyzed the differences and the effectiveness of weight function and ocean tide corrections in the FCN parameter detection using synthetic data, SG data from thirty-one stations, and the 10 celestial pole offset (CPO) series. The results show that significant discrepancies are caused by different computing options for a single SG station. The stacking method, which results in a variation of 0.24–5 sidereal days (SDs) in the FCN period (T) and 103-104 in the quality factor (Q) due to the selection of the weighting function and the ocean tide model (OTM), can effectively suppress this influence. The statistical analysis results of synthetic data shows that although different weight choices, while adjusting the proportion of diurnal tidal waves involved, do not significantly improve the accuracy of fitted FCN parameters from gravity observations. The study evaluated a series of OTMs using the loading correction efficiency. The fitting of FCN parameters can be improved by selecting the mean of appropriate OTMs based on the evaluation results. Through the estimation of the FCN parameters based on the forced nutation, it was found that the weight function P1 is more suitable than others, and different CPO series (after 2009) resulted in a difference of 0.4 SDs in the T and of 103 in the Q. We estimated the FCN parameters for SG: (T = 430.4 ± 1.5 SDs and Q = 1.52 × 104 ± 2.5 × 103) and for VLBI: (T = 429.8 ± 0.7 SDs, Q = 1.88 × 104 ± 2.1 × 103).

On December 18, 2023, the M-S 6.2 Jishishan earthquake occurred in the northeastern region of the Qinghai-Xizang Plateau, causing heavy casualties and property damage in Gansu and Qinghai Provinces. In this study, we integrate space imaging geodesy, finite fault inversion, and back-projection methods to decipher its rupture property, including fault geometry, coseismic slip distribution, rupture direction, and propagation speed. The results reveal that the seismogenic fault dips to the southwest at an angle of 29 degrees. The major slip asperity is dominated by reverse slip and is concentrated within a depth range of 7-16 km, which explains the significant uplift near the epicenter observed by both the Sentinel-1 ascending and descending InSAR data. Moreover, the teleseismic array waveforms indicate a northwest propagating rupture with an overall slow rupture velocity of similar to 1.91 km/s (AK array) or 1.01 km/s (AU array).