Journal of Geodesy

Latest harmonic developments of the Earth tide-generating potential (TGP), e.g., HW95 (Hartmann and Wenzel in Geoph Res Lett 22:3553, 1995), RATGP95 (Roosbeek in Geophys J Int 126:197, 1996), KSM03 (Kudryavtsev in J Geodesy 77:829, 2004), include a number of terms caused by the joint effect of the Earth’s polar flattening (that can be numerically described by the \({J}_{2}\) geopotential coefficient) and the Moon/the Sun gravitational attraction. In the present study, we additionally consider the effect of the Earth’s equatorial flattening due to the Earth’s triaxiality. Explicit analytical expressions for the relevant part of the TGP are derived. New terms of the TGP development due to the Earth’s triaxial figure are found. Amplitudes of nineteen of them exceed the threshold level of 10–8 m2s−2 used by the modern tidal potential catalogs. Three of the new terms have the frequency sign opposite to that of the Earth rotation. It is not the case for any previously known term of the Earth TGP development. Every term has a new feature that an integer multiplier of the mean local lunar time used in its argument is not equal to the order of the spherical harmonic associated with the term. It necessitates a relevant modification of the standard HW95 format for representing the Earth TGP. The new terms are suggested for including in the current and future tidal potential catalogs.

Partial ambiguity resolution (PAR) has been widely adopted in real-time kinematic (RTK) and precise point positioning with augmentation from continuously operating reference station (PPP-RTK). However, current PAR methods, either in the position domain or the ambiguity domain, suffer from high false alarm and miss detection, particularly in challenging environments with poor satellite geometry and observations contaminated by non-line-of-sight (NLOS) effects, gross errors, biases, and high observation noise. To address these issues, a PAR method based on machine learning is proposed to significantly improve the correct fix rate and positioning accuracy of PAR in challenging environments. This method combines two support vector machine (SVM) classifiers to effectively identify and exclude ambiguities those are contaminated by bias sources from PAR without relying on satellite geometry. The proposed method is validated with three vehicle-based field tests covering open sky, suburban, and dense urban environments, and the results show significant improvements in terms of correct fix rate and positioning accuracy over the traditional PAR method that only utilizes ambiguity covariances. The fix rates achieved with the proposed method are 93.9%, 81.9%, and 93.1% with the three respective field tests, with no wrong fixes, compared to 72.8%, 20.9%, and 16.0% correct fix rates using the traditional method. The positioning error root mean square (RMS) is 0.020 m, 0.035 m, and 0.056 m in the east, north, and up directions for the first field test, 0.027 m, 0.080 m, and 0.126 m for the second field test, and 0.035 m, 0.042 m, and 0.071 m for the third field test. In contrast, only decimeter- to meter-level accuracy was obtainable with these datasets using the traditional method due to the high wrong fix rate. The proposed method provides a correct and fast time-to-first-fix (TTFF) of 3–5 s, even in challenging environments. Overall, the proposed method offers significant improvements in positioning accuracy and ambiguity fix rate with high reliability, making it a promising solution for PAR in challenging environments.

The development of optical clocks has experienced significant acceleration in recent years, positioning them as one of the most promising quantum optical sensors for next-generation gravimetric missions (NGGMs). This study investigates the feasibility of retrieving the temporal gravity field via improved optical clocks through a closed-loop simulation. It evaluates optical clock capabilities in temporal gravity field inversion by considering the clock noise characteristics, designing satellite formations, and simulating the performance of optical clocks. The results indicate that optical clocks exhibit higher sensitivity to low-degree gravity field signals. However, when the optical clock noise level drops below 1 × 10−19 \(/\sqrt{\uptau }\) (τ being the averaging time in seconds) in the satellite-to-ground (SG) mode or below 1 × 10−20 \(/\sqrt{\uptau }\) in the satellite-to-satellite (SS) mode, atmospheric and oceanic (AO) errors become the dominant source of error. At this noise level, optical clocks can detect time-variable gravity signals up to approximately degree 30. Compared to existing gravity measurement missions such as GRACE-FO, optical clocks exhibit consistent results in detecting signals below degree 20. If the orbital altitude is reduced to 250 km, the performance of optical clocks across all degrees aligns with the results of GRACE-FO. Furthermore, the study reveals that lowering the orbital altitude of satellite-based optical clocks from 485 to 250 km improves results by an average of 33%. Switching from the SS mode to the SG mode results in an average improvement of 51%, while each order-of-magnitude improvement in clock precision enhances results by an average of 60%. In summary, these findings highlight the tremendous potential of optical clock technology in determining Earth’s temporal gravity field and provide crucial technological support for NGGMs.

Lateral inhomogeneity in the Earth’s mantle affects the tidal response. The current study reformulates the expressions for estimating the lateral inhomogeneity effects on tidal gravity with respect to the unperturbed Earth and supplements some critical derivation process to enhance the methodology. The effects of lateral inhomogeneity are calculated using several real Earth models. By considering the collective contributions of seismic wave velocity disturbances and density disturbance, the global theoretical changes of semidiurnal gravimetric factor are obtained, which vary from − 0.22 to 0.22% compared to those in a layered Earth model, about 1/2 of the ellipticity’s effect. The gravity changes caused by lateral-inhomogeneous disturbances are also computed and turn out to be up to 0.16% compared to the changes caused by tide-generating potential. The current study compares the influences of lateral inhomogeneity with rotation and ocean tide loading. The results indicate that the rotation and ellipticity on tidal gravity are the most dominant factors, the ocean tide loading is the moderate one, and the lateral inhomogeneity in the mantle has the least significant influence. Moreover, an anti-correlation between the effective elastic thickness and gravimetric factor change caused by lateral inhomogeneity is found, implying that it is difficult to generate tidal response at regions with high rigidity. We argue that the gravimetric factor change can be used as an effective indicator for lithospheric strength.

The regional sea level budget and interannual sea level changes around Taiwan and Philippines are studied using altimetry, GRACE, and in-situ hydrographic data during 1993‒2021. Results show that the average sea level trend around Taiwan and Philippines during 1993–2021 derived from the altimetric data is 3.6 ± 0.2 mm/yr. Over 2002–2021, the study shows closure of sea level budget in the eastern ocean of Taiwan and Philippines within the observed data uncertainties, and the ocean mass accounts for 88%–100% of the observed sea level rise. In contrast, the sea level budget is not closed in the western ocean of Taiwan and Philippines, probably due to the lack of complete coverage by in-situ ocean observing systems. In addition, both regional sea level anomalies and their steric component around Taiwan and Philippines exhibit pronounced interannual and decadal variabilities. The trade wind stress associated with El Niño–Southern Oscillation and Pacific Decadal Oscillation offers a compelling explanation for the interannual and decadal signals of sea level anomalies in the southern ocean of Taiwan, with negative correlations of − 0.78 to − 0.64, indicating that trade wind stress makes a negative contribution to interannual-to-decadal sea level variability. In the northwestern ocean of Taiwan, the sea level variation is strongly influenced by the local monsoon system and shallow bathymetry with an annual amplitude of 90.3 ± 2.9 mm, larger than those in other regions around Taiwan and Philippines, where ocean mass is dominant with a high correlation with the sea level (+ 0.75 to + 0.78).

Numerical methods, like the finite element method (FEM) or finite volume method (FVM), are widely used to provide solutions in many boundary value problems. In previous studies, these numerical methods have also been applied in geodesy but demanded extensive computations because the upper boundary condition was usually set up at the satellite orbit level, hundreds of kilometers above the Earth. The relatively large distances between the lower boundary of the Earth's surface and the upper boundary exacerbate the computation loads because of the required discretization in between. Considering that many areas, such as the US, have uniformly distributed airborne gravity data just a few kilometers above the topography, we adapt the upper boundary from the satellite orbit level to the mean flight level of the airborne gravimetry. The significant decrease in the domain of solution dramatically reduces the large computation demand for FEM or FVM. This paper demonstrates the advantages of using FVM in the decreased domain in simulated and actual field cases in study areas of interest. In the simulated case, the FVM numerical results show that precision improvement of about an order of magnitude can be obtained when moving the upper boundary from 250 to 10 km, the upper altitude of the GRAV-D flights. A 2–3 cm level of accurate quasi-geoid model can be obtained for the actual datasets depending on different schemes used to model the topographic mass. In flat areas, the FVM solution can reach to about 1 cm precision, which is comparable with the counterparts from classical methods. The paper also demonstrates how to find the upper boundary if no airborne data are available. Finally, the numerical method provides a 3D discrete representation of the entire local gravity field instead of a surface solution, a (quasi) geoid model.

The Kalman filter stands as one of the most widely used methods for recursive parameter estimation. However, its standard formulation typically assumes that all state parameters avail initial values and dynamic models, an assumption that may not always hold true in certain applications, particularly in global navigation satellite system (GNSS) data processing. To address this issue, Teunissen et al. (2021) introduced a generalized Kalman filter that eliminates the need for initial values and allows linear functions of parameters to have dynamic models. This work proposes a least-squares approach to reformulate the generalized Kalman filter, enhancing its applicability to GNSS data processing when the parameter dimension varies with satellite visibility changes. The reformulated filter, named generalized least-squares filter, is equivalent to the generalized Kalman filter when all state parameters are recursively estimated. In this case, we demonstrate how both the generalized Kalman filter and the generalized least-squares filter adaptatively manage newly introduced or removed parameters. Specifically, when estimation is limited to parameters with dynamic models, the generalized least-squares filter extends the generalized Kalman filter by allowing the dimension of estimated parameters to vary over time. Moreover, we introduce a new element of least-squares smoothing, creating a comprehensive system for prediction, filtering, and smoothing. To verify, we conduct simulated kinematic and vehicle-borne kinematic GNSS positioning using the proposed generalized least-squares filter and compare the results with those from the standard Kalman filter. Our findings show that the generalized least-squares filter delivers better results, maintaining the positioning errors at the centimeter level, whereas the Kalman-filter-based positioning errors exceed several decimeters in some epochs due to improper initial values and dynamic models. Moreover, the normal equation reduction strategy employed in the generalized least-squares filter improves computational efficiency by 23% and 32% in simulated kinematic and vehicle-borne kinematic positioning, respectively. The generalized least-squares filter also allows for the flexible adjustment of smoothing window lengths, facilitating successful ambiguity resolution in several epochs. In conclusion, the proposed generalized least-squares filter offers flexibility for various GNSS data processing scenarios, ensuring both theoretical rigor and computational efficiency.

Vertical land motion and the redistribution of masses within and on the surface of the Earth affect the Earth’s gravity field. Hence, studying the ratio between temporal changes of the surface gravity \(\left( {\dot{g}} \right)\) and height ( \(\dot{h}\) ) is important in geoscience, e.g., for reduction of gravity observations, assessing satellite gravimetry missions, and tuning vertical land motion models. Sjöberg and Bagherbandi (2020) estimated a combined ratio of \(\dot{g}/\dot{h}\) in Fennoscandia based on relative gravity observations along the 63 degree gravity line running from Vågstranda in Norway to Joensuu in Finland, 688 absolute gravity observations observed at 59 stations over Fennoscandia, monthly gravity data derived from the GRACE satellite mission between January 2003 and August 2016, as well as a land uplift model. The weighted least-squares solution of all these data was \(\dot{g}/\dot{h}\)  =  − 0.166 ± 0.011 μGal/mm, which corresponds to an upper mantle density of about 3402 ± 95 kg/m3. The present note includes additional GRACE data to June 2017 and GRACE Follow-on data from June 2018 to November 2023. The resulting weighted least-squares solution for all data is \(\dot{g}/\dot{h}\)  =  − 0.160 ± 0.011 μGal/mm, yielding an upper mantle density of about 3546 ± 71 kg/m3. The outcomes show the importance of satellite gravimetry data in Glacial Isostatic Adjustment (GIA) modeling and other parameters such as land uplift rate. Utilizing a longer time span of GRACE and GRACE Follow-on data allows us to capture fine variations and trends in the gravity-to-height ratio with better precision. This will be useful for constraining and adjusting GIA models and refining gravity observations.

LARES-2 is a new geodetic satellite designed for high-accuracy satellite laser ranging. The orbit altitude of LARES-2 is similar to that of LAGEOS-1, whereas the inclination angle of 70° complements the LAGEOS-1 inclination of 110°; hence, both satellites form the butterfly configuration for the verification of the Lense–Thirring effect. Although the major objective of LARES-2 is testing general relativity, LARES-2 substantially contributes to geodesy in terms of the realization of terrestrial reference frames, recovery of the geocenter motion, pole coordinates, length-of-day, and low-degree gravity field coefficients. We analyze the first 1.5 years of LARES-2 data and test different empirical orbit models for LARES-2 with and without co-estimating low-degree gravity field coefficients to find the best combination strategy with LAGEOS satellites. We found that LARES-2 orbit determination is more accurate than that of LAGEOS-1/2 due to a different satellite construction consisting of a solid sphere with no inner structure. Neither the correction for D0 nor the empirical once-per-revolution along-track accelerations SC/SS have to be estimated for LARES-2 when co-estimating gravity field coefficients. The only empirical parameter needed for LARES-2 is the constant along-track acceleration S0 to compensate for the Yarkovsky–Schach effect. On the contrary, for LAGEOS-1/2, the non-gravitational perturbations affect C30 and Z geocenter estimates when once-per-revolution parameters are not estimated. LARES-2 does not face this issue. LARES-2 improves the formal errors of the Z geocenter component by up to 59% and C20 by up to 40% compared to the combined LAGEOS-1/2 solutions and provides C30 estimates unaffected by thermal orbit modeling issues.

For precise orbit determination (POD) and precise applications with POD products, one of the critical issues is the modeling of non-conservative forces acting on satellites. Since the official publication of Galileo satellite metadata in 2017, analytical models including the box-wing model and thermal thrust models have been established to absorb a substantial amount of solar radiation pressure (SRP) and thermal thrust. These models serve as the foundation for the best overall modeling approach, combining the analytical box-wing model and thermal thrust model with parameterization of the remaining non-conservative perturbing forces using various optimized Empirical CODE Orbit Models (ECOMs) of the Center for Orbit Determination in Europe (CODE). Firstly, we have demonstrated the significance of the second-order signals in the D direction and the first-order signals in the B direction through spectral analyses of the pure box-wing model, which are consistent with the currently recommended 7-parameter Empirical CODE Orbit Model 2 (ECOM2). In spite of this, we still found that degradation in orbit accuracy frequently occurs during deep eclipse seasons when using the ECOM2 model. We confirm a high-frequency signal existing in the fluctuating orbit overlap differences through the spectral analysis. Considering this, the ECOM2 force model should be extended to higher order and adapted to absorb the remaining effects of potential perturbing forces. After extending the ECOM2 force model to the sixth order in the Sun direction, we demonstrated the significance of fourth- and sixth-order sine terms for deep eclipses. Due to the higher-order periodic terms, the averaged RMS values of orbit overlap difference over deep eclipses can be reduced from 5.3, 10.8, and 23.8 cm to 3.2, 3.9, and 9.9 cm for in-orbit validation (IOV) satellites, from 5.0, 8.6, and 17.7 cm to 3.0, 3.0, and 7.1 cm for the first generation of full operational capability (FOC-1) satellites, and from 5.4, 8.6, and 19.0 cm to 3.6, 3.6, and 7.4 cm for the second generation of FOC (FOC-2) satellites, in the radial, cross-track, and along-track directions, respectively. Fluctuations with a peak amplitude of approximately 0.4 nm/s2 in the bias in the solar panel axis (Y) direction (Y-bias) are effectively mitigated by the higher-order terms. Due to the higher-order terms, the vertical positioning errors during kinematic precise point positioning (PPP) convergence can be improved from 42.3 to 37.1 cm at the 95.5% confidence level. Meanwhile, a low correlation level of up to 0.02 is found between the newly introduced higher-order parameters and earth rotation parameters (ERPs).

Carrier phase integer ambiguities must be determined for BDS-3/inertial navigation system (INS) tightly coupled (TC) integration to achieve centimetre-level positioning accuracy. However, cycle slip breaks the consistency of the integer ambiguities. Conventional multi-frequency cycle slip methods use the pseudorange; thus, requiring improvement when applied to kinematic situations. Furthermore, a concise and nonprior information-dependent model is crucial for real-time processing. In this study, an inertial-aided BDS-3 cycle slip detection and repair (I-CDR) method was developed. First, a BDS-3/INS TC model with I-CDR was created. The ionospheric delays were modelled as part of the TC states; therefore, they could be estimated and eliminated. Investigations were conducted on the effects of carrier phase noise, residual ionosphere delay, and INS-predicted position error on combined cycle slip detection (CCD) accuracy. The optimal CCDs under various frequency available configurations were determined. The effectiveness of I-CDR was demonstrated using land vehicle test data. The false alarm ratio was less than 1.0%, and the missed detection ratio was almost zero even in situations with challenging abundant 1-cycle slips in random epochs. Furthermore, the right determination ratio reached 100%. In addition, BDS-3 signal loss-recovery cases were simulated, and all cycle slips for all satellites could be repaired within 40s. I-CDR exhibits outstanding cycle slip detection and repair performance for dense 1-cycle slip and signal loss-recovery cases, demonstrating its suitability for BDS-3/INS TC integration.

This paper focuses on the retrieval of refractivity fields from GNSS measurements by means of least-squares collocation. Collocation adjustment estimates parameters that relate delays and refractivity without relying on a grid. It contains functional and stochastic models that define the characteristics of the retrieved refractivity fields. This work aims at emphasizing the capabilities and limitations of the collocation method in modeling refractivity and to present it as a valuable alternative to GNSS tomography. Initially, we analyze the stochastic models in collocation and compare the theoretical errors of collocation with those of tomography. We emphasize the low variability of collocation formal variances/covariances compared to tomography and its lower dependence on a-priori fields. Then, based on real and simulated data, we investigate the importance of station resolution and station heights for collocation. Increasing the network resolution, for example, from 10 to 2 km, results in improved a-posteriori statistics, including a 10% reduction in the error statistic for the retrieved refractivity up to 6 km. In addition, using additional stations at higher altitudes has an impact on the retrieved refractivity fields of about 1 ppm in terms of standard deviation up to 6 km, and a bias reduction of more than 3 ppm up to 3 km. Furthermore, we compare refractivity fields retrieved through tomography and collocation, where data of the COSMO weather model are utilized in a closed-loop validation mode to simulate tropospheric delays and validate the retrieved profiles. While tomography estimates are less biased, collocation captures relative changes in refractivity more effectively among the voxels within one height level. Finally, we apply tomography and collocation to test their capabilities to detect an approaching weather front. Both methods can sense the weather front, but their atmospheric structures appear more similar when the GNSS network has a well-distributed height coverage.

Deformation phenomena on Earth are inherently three dimensional. With SAR interferometry (InSAR), in many practical situations the maximum number of observations is two (ascending and descending), resulting in an infinite number of possible displacement estimates. Here we propose a practical solution to this underdeterminancy problem in the form of the strapdown approach. With the strapdown approach, it is possible to obtain “3D-global/2D-local” solutions, by using minimal and largely undisputed contextual information, on the expected driving mechanisms and/or spatial geometry. It is a generic method that defines a local reference system with transversal, longitudinal, and normal (TLN) axes, with displacement occurring in the transversal-normal plane only. Since the orientation of the local frame is based on the physics of the problem at hand, the strapdown approach gives physically more relevant estimates compared to conventional approaches. Moreover, using an a-priori uncertainty approximation on the orientation of the local frame it is possible to assess the precision of the final estimates. As a result, appropriate cartographic visualization using a vector map with confidence ellipses enables an improved interpretation of the results.

Accurate information on satellite antenna phase center offsets (PCOs) and phase variations (PVs) is indispensable for high-precision geodetic applications. In the absence of consistent pre-flight calibrations, satellite antenna PCOs and PVs of global navigation satellite systems are commonly estimated based on observations from a global network, constraining the scale to a given reference frame. As part of this estimation, flatness and zero-mean conditions need to be applied to unambiguously separate PCOs, PVs, and constant phase ambiguities. Within this study, we analytically investigate the impact of different boresight-angle-dependent weighting functions for PV minimization, and we compare antenna models generated with different observation-based weighting schemes with those based on uniform weighting. For the case of the GPS IIR/-M and III satellites, systematic differences of 10 mm in the PVs and 65 cm in the corresponding PCOs are identified. In addition, new antenna models for the different blocks of BeiDou-3 satellites in medium Earth orbit are derived using different processing schemes. As a drawback of traditional approaches estimating PCOs and PVs consecutively in distinct steps, it is shown that different, albeit self-consistent, PCO/PV pairs may result depending on whether PCOs or PVs are estimated first. This apparent discrepancy can be attributed to potentially inconsistent weighting functions in the individual processing steps. Use of a single-step process is therefore proposed, in which a dedicated constraint for PCO-PV separation is applied in the solution of the normal equations. Finally, the impact of neglecting phase patterns in precise point positioning applications is investigated. In addition to an overall increase of the position scatter, the occurrence of systematic height biases is illustrated. While observation-based weighting in the pattern estimation can help to avoid such biases, the possible benefit depends critically on the specific elevation-dependent weighting applied in the user’s positioning model. As such, the practical advantage of such antenna models would remain limited, and uniform weighting is recommended as a lean and transparent approach for the pattern estimation of satellite antenna models from observations.

An attractive feature of PPP-RTK is the possibility of reducing the amount of data that needs to be transferred to users. By leveraging the state-space Representation (SSR) of the corrections, the correction provider (i.e., a GNSS network) can consider distinct transfer rates for each of the individual corrections according to their temporal characteristics. Reducing the transfer rates comes at the cost of delivering time-delayed corrections, urging the user to time predict the corrections to bridge the gap between the corrections’ generation time and the positioning time. Consequently, the user Kalman filter needs to be equipped with a strategy to account for the errors caused by such predictions, minimizing the precision loss of the user parameter solutions. In this contribution, we apply a processing strategy for both the network and user filters to handle the latency of corrections. This enables the network to update corrections over longer time-intervals. To have the strategy applicable to regional networks, an ionosphere-weighted model is adopted for the corresponding observations, delivering minimum-variance spatially predicted ionospheric corrections to users. It is shown that certain components of the network filter’s dynamic model are duplicated and should be excluded from processing. To illustrate the performance of the strategy at work, three globally distributed regional networks are employed, and maximum correction latencies to meet different positioning criteria are evaluated. In terms of both the positioning precision and time-to-first-fix (TTFF), the strategy is numerically shown to outperform the user processing case in which the uncertainty of corrections is discarded.

We propose a benchmark for comparing gap-filling techniques used on global time-variable gravity field time-series. The Gravity Recovery and Climate Experiment (GRACE) and the GRACE Follow-On missions provide products to study the Earth’s time-variable gravity field. However, the presence of missing months in the measurements poses challenges for understanding specific Earth processes through the gravity field. We reproduce, adapt, and compare satellite-monitoring and interpolation techniques for filling these missing months in GRACE and GRACE Follow-On products on a global scale. Satellite-monitoring techniques utilize solutions from Swarm and satellite laser ranging, while interpolation techniques rely on GRACE and/or Swarm solutions. We assess a wide range of interpolation techniques, including least-squares fitting, principal component analysis, singular spectrum analysis, multichannel singular spectrum analysis, auto-regressive models, and the incorporation of prior data in these techniques. To inter-compare these techniques, we employ a remove-and-restore approach, removing existing GRACE products and predicting missing months using interpolation techniques. We provide detailed comparisons of the techniques and discuss their strengths and limitations. The auto-regressive interpolation technique delivers the best score according to our evaluation metric. The interpolation based on a least-squares fitting of constant, trend, annual, and semi-annual cycles offers a simple and effective prediction with a good score. Through this assessment, we establish a starting benchmark for gap-filling techniques in Earth’s time-variable gravity field analysis.

Inverting fault geometry and slip distribution simultaneously with geodetic observations based on Bayesian theory is becoming increasingly prevalent. A widely used approach, proposed by (Fukuda and Johnson, Geophys J Int 181:1441–1458, 2010) (F-J method), employs the least-squares method to solve the linear parameters of slip distribution after sampling the nonlinear parameters, including fault geometry, data weights and smoothing factor. Here, we present a modified version of the F-J method (MF-J method), which treats data weights and the smoothing factor as hyperparameters not directly linked to surface deformation. Additionally, we introduce the variance component estimation (VCE) method to resolve these hyperparameters. To validate the effectiveness of the MF-J method, we conducted inversion tests using both synthetic data and a real earthquake case. In our comparison of the MF-J and F-J methods using synthetic experiments, we found that the F-J method's inversion results for fault geometry were highly sensitive to the initial values and step sizes of hyperparameters, whereas the MF-J method exhibited greater robustness and stability. The MF-J method also exhibited a higher and more stable acceptance rate, enabling convergence to simulated values and ensuring greater accuracy of the parameter estimation. Furthermore, treating the fault length and width as unknown parameters and solving them simultaneously with other fault geometry parameters and hyperparameters using the MF-J method successfully resolved the issue of non-uniqueness in fault location solutions caused by the excessively large no-slip areas. In the 2017 Mw 7.3 Sarpol-e Zahab earthquake case study, the MF-J method produced a fault slip distribution with low uncertainty that accurately fit surface observation data, aligning with results from other research institutions. This validated the method's applicability and robustness in real-world scenarios. Additionally, we inferred that the second slip asperity was caused by early afterslip.

The 3D ionosphere structure is of interest in many fields such as radio frequency communication and global navigation satellite system (GNSS) applications. However, the limited temporal and spatial coverage of measurements poses a challenge for 3D electron density modeling. To overcome this challenge, we explore the use of kriging interpolation technique. The kriging interpolation is performed to obtain 3D representation of the ionosphere over electron density measurements retrieved by GNSS radio-occultation (RO) data. RO measurements are first reduced to “shape function,” the ratio of electron density to vertical total electron content (VTEC), aiming to create a background model. Then, the empirical residual semivariogram is analyzed for variation characteristics of the shape functions under different solar geomagnetic conditions. Finally, 3D kriging is adopted for shape function interpolation. Compared to the modeling results without kriging, the maximum root mean square error (RMSE) reduction reaches \(3.4\times {10}^{-4}~\text {km}^{-1}\) , which amounts to \(3.4\times {10}^{11}~\text {el/m}^{3}\) of electron density when VTEC is assumed as 100 TECU. This improvement accounts for 17.8% of root mean square (RMS) of shape function.