Earth Model Research Articles

Щодо підвищення точності аналітичних та чисельних методів геодезичних та картометричних операцій

Purpose of this work is to determine and justify the use of analytical and numerical geodetic and cartometric methods on the reference ellipsoid, as well as their calculation accuracy in the geographic information environment. Methodology. The research compiled a register of geodetic and cartometric methods used in geodetic practice and implemented in modern geographic information systems. Standard tools in GIS often use approximate numerical methods, which affects the accuracy of models of geospatial objects in the GIS environment. Therefore, we have analysed and established for each operation of geodetic and cartometric methods a mathematical model that determines a particular cartometric property with maximum accuracy either by analytical or numerical methods with the number of terms in the binomial series of 6 or more. Results. The author proposed 10 operations of the geodetic method and 3 operations of the cartometric method, for which mathematical models and their accuracy were established and substantiated with their corresponding implementation in the MATLAB v. R2018a. The defined list of geodetic and cartometric operations made it possible to move away from the classification of distance lengths that influenced the further use of certain surfaces (map projection plane, sphere, spheroid, or reference ellipsoid) and mathematical models of operations. The presented mathematical models allow performing the relevant geodetic and cartometric methods with maximum accuracy using modern computer technologies. The mathematical models of geodetic and cartometric methods are investigated, which have practically no limitations for achieving the required accuracy, especially for large and ultra-large distances. The scientific novelty of the research is to define and justify a clear list of mathematical models of numerical and analytical geodetic and cartometric methods instead of cartometric methods on the map and standard methods of instrumental GIS; using the surface of the reference ellipsoid, and not just cartographic projections, spheroid or sphere. The practical significance of the research lies in the use of numerical and analytical geodetic and cartometric methods that significantly increase the accuracy of operations in these works, as well as in the creation/updating of digital topographic maps, navigation and route planning, etc. The research results can be concluded that well-founded mathematical models will ensure an increase in the accuracy of computational operations taking into account the curvature of the Earth in all sectors and areas of the economy, which will affect the quality of accounting and monitoring of relevant objects, integration and geospatial analysis of heterogeneous geospatial data, improve the quality (topological consistency) of geospatial data, etc

Equation of state and thermodynamic properties of liquid Fe-O in the Earth’s outer core

Equation of state (EoS) plays a crucial role in the prediction of the composition of the outer core. Here, we calculated pressure (P)-volume (V)-temperature (T) data of liquid iron-oxygen alloys (Fe-X wt.% O, X = 0, 2.8, 6.1, and 9.9) under the outer core conditions (∼136–330 GPa, 4000–6000 K) by first-principles molecular dynamics simulations. We established an EoS for liquid Fe-O alloys with parameters including P, T, V, and O concentrations. Consequently, thermodynamic properties of liquid Fe-O alloys such as density (ρ), thermal expansion coefficient, isothermal and adiabatic bulk modulus, and sound velocity (VP) are calculated. To constrain the O content, we predicted the ρ-P and VP-P profiles along the geotherm and compared them with data from the Preliminary Reference Earth Model (PREM). We conclude that the adiabatic T profile as a function of depth affects the prediction of O content dramatically. With several anchored TICB, the composition of Fe-6.1 wt% O matches the PREM data with an acceptable range of error. But strictly speaking, the distribution in the outer core is probably uneven. In such case, we state that the O content in the outer core cannot be higher than approximately 6.1 wt%.

Introduction to the Distributional Finite Difference Method for 3D Seismic Wave Propagation and Comparison With the Spectral Element Method

AbstractWe have extended the distributional finite difference method (DFDM) to simulate the seismic‐wave propagation in 3D regional earth models. DFDM shares similarities to the discontinuous finite element method on a global scale and to the finite difference method locally. Instead of using linear staggered finite‐difference operators, we employ DFDM operators based on B‐splines and a definition of derivatives in the sense of distributions, to obtain accurate spatial derivatives. The weighted residuals method used in DFDM's locally weak formalism of spatial derivatives accurately and naturally accounts for the free surface, curvilinear meshing, and solid‐fluid coupling, for which it only requires setting the shear modulus and the continuity condition to zero. The computational efficiency of DFDM is comparable to the spectral element method (SEM) due to the more accurate mass matrix and a new band‐diagonal mass factorization. Numerical examples show that the superior accuracy of the band‐diagonal mass and stiffness matrices in DFDM enables fewer points per wavelength, approaching the spectral limit, and compensates for the increased computational burden due to four Lebedev staggered grids. Specifically, DFDM needs 2.5 points per wavelength, compared to the five points per wavelength required in SEM for 0.5% waveform error in a homogeneous model. Notably, while maintaining the same accuracy in the solid domain, DFDM demonstrates superior accuracy in the fluid domain compared to SEM. To validate its accuracy and flexibility, we present various 3D benchmarks involving homogeneous and heterogeneous regional elastic models and solid‐fluid coupling in both Cartesian and spherical settings.

Coupled Inversion of Amplitudes and Traveltimes of Primaries and Multiples for Monochannel Seismic Surveys

Engineers need to know properties of shallow marine sediments to build piers, pipelines and even offshore windfarms. We present a method for estimating the density, P velocity and thickness of these sediments. The traveltime inversion of primary and multiple reflections enables their semiquantitative estimation in marine surveys when using a minimal acquisition system such as a monochannel Boomer. Picking errors, ambient noise and interfering events lead to significant errors in the estimates. Similar, albeit milder, instabilities occur when inverting the signal amplitudes to determine the reflectivity of the layer interfaces. In this paper, we introduce a coupling between the separate inversion of amplitudes and traveltimes to obtain a better Earth model. The P velocity shows up in two stable terms provided by the separate inversions: the acoustic impedance of shallow sediments (through the amplitudes) and the transit time across the sediment layer (through the traveltimes). We couple the two inversion engines by imposing a smoothness condition on velocity and density and thickness of the layer while keeping the impedance and traveltime constant. We thus exploit the ambiguity of the solution to introduce geological criteria and reduce the noise contribution. We validated the proposed method with synthetic and real data.

The exact implementation of a spherical harmonic model for gravimetric quantities

Abstract This note presents formulas to evaluate a spherical harmonic model of Earth’s gravitational potential for essential gravimetric quantities without spherical and linear approximation. Typically, 10–13 significant digits of numerical accuracy for such computations are obtained over the globe using EGM2008 with FORTRAN 77 code that is also provided.

George McCready Price and Seventh-Day Adventism’s Influence on the Rise of Modern Creationism

Abstract Self-taught geologist and early convert to Seventh-day Adventism, George McCready Price played a significant, yet largely overlooked role in the rise of the modern creationist movement—and, to some degree, can be seen as the lynchpin leading to mainstream adoption of literal readings of the Christian scriptures. Influenced greatly by the writings of Ellen White, the prophetess of the Adventist movement, Price developed his ideas on flood geology which eventually became a cornerstone of creationist thinking. Price’s commitment to biblical literalism, influenced by White’s visions, shaped his worldview and led him to reconcile his interest in geology with the Adventist model of Earth’s history. His central argument for creationism was based on his belief in a worldwide flood, as described in the biblical narrative, which he saw as having shaped most of the Earth’s geological features. Although Price’s contributions were initially overlooked, his ideas gained prominence through the work of John Whitcomb and Henry Morris, who popularized his concepts in their influential book, The Genesis Flood. Price’s effect on modern creationism highlights the intersection of fundamentalist Christianity, particularly Adventism, with the evolving tension between theology and science in the 20th Century.

Tidal Calibration of the Gladwin Tensor Strain Monitor (GTSM) Array in Taiwan

AbstractTo ensure the accuracy and reliability of crustal strain measures, sensors require a thorough calibration. In Taiwan, the complicated dynamics of surface and subsurface hydrological processes under semi-tropical climate conditions conjugated with the rough surface topography could have impacted strainmeter deployment, pushing the installation conditions astray from the optimal ones. Here, we analyze the complex response of 11 Gladwin Strain Monitor (GTSM) strainmeter type deployed in north and central Taiwan and we propose a novel calibration methodology which relies on waveform modeling of Earth and ocean tidal strain-related deformations. The approach is completely data-driven, starting from a simple calibration framework and progressively adding complexity in the model depending on the quality of the data. However, we show that a simple quasi-isotropic model (three calibration factors) is generally suitable to resolve the orientation and calibration of 8 instruments out of 11. We also highlight the difficulty of clearly defining the behavior of instruments that are highly affected by hydrological forcing.

Biodesign: An Innovative Technology for Sustainable Development Goals in Health and Environment

Biodesign is a tool that integrates diverse fields of biology with design for the development of a good lifestyle and, sustainable environment along with futuristic support to the society and ecosystem. Using the approach, Identify-Invent-Implement, the biodesigners identify and unravel the unmet needs of society in the field of healthcare, environment, and architecture sectors by developing novel products and improving technologies with the incorporation of design into it. With evolving technology, biodesign has created an impact globally in the ecology by preventing the ecological imbalance, creating an impact on the psychological behaviour of humans by boosting their confidence and thus influencing one’s behaviour. Biodesign crossed various technological barriers by developing devices such as hearing, and vision aids, prosthetics, implants, and eco-friendly cements and also used microorganisms to enhance the quality of the bricks. Biodesign has crossed the limits of earth curvature and reached space by innovation like future space colonization and protecting indigenous species.

CSDMS Data Components: data–model integration tools for Earth surface processes modeling

Abstract. Progress in better understanding and modeling Earth surface systems requires an ongoing integration of data and numerical models. Advances are currently hampered by technical barriers that inhibit finding, accessing, and executing modeling software with related datasets. We propose a design framework for Data Components, which are software packages that provide access to particular research datasets or types of data. Because they use a standard interface based on the Basic Model Interface (BMI), Data Components can function as plug-and-play components within modeling frameworks to facilitate seamless data–model integration. To illustrate the design and potential applications of Data Components and their advantages, we present several case studies in Earth surface processes analysis and modeling. The results demonstrate that the Data Component design provides a consistent and efficient way to access heterogeneous datasets from multiple sources and to seamlessly integrate them with various models. This design supports the creation of open data–model integration workflows that can be discovered, accessed, and reproduced through online data sharing platforms, which promotes data reuse and improves research transparency and reproducibility.

Accurate Height Determination in Uneven Terrains with Integration of Global Navigation Satellite System Technology and Geometric Levelling: A Case Study in Lebanon

The technology for determining a point’s coordinates on the earth’s surface using the global navigation satellite system (GNSS) is becoming the norm along with ground-based methods. In this case, determining coordinates does not cause any particular difficulties. However, to identify normal heights using this technology with a given accuracy, special research is required. The fact is that satellite determinations of geodetic heights (h) over an ellipsoid surface differ from ground-based measurements of normal height (HN) over a quasi-geoid surface by a certain value called quasi-geoid height or height anomaly (ζ). In relation to determining heights of a certain territory, the concept of geoid height (N) is usually operated when dealing with a geoid model. In this work, geodetic and normal heights are determined for five control points in three different regions in Lebanon, where measurements are carried out using GNSS technology and geometric levelling. The obtained quasi-geoid heights are compared with geoid heights derived from the global Earth model EGM2008. The results obtained showed that, in the absence of gravimetric data, the combination of global Earth model data, geometric levelling for selected areas, and satellite determinations allows for the creation of a highly accurate altitude network for mountainous areas.

Analysis of the seismicity related to the tidal stress prior to the M w7.8 Kahramanmaras, Türkiye earthquake

Based on the preliminary reference Earth model, we calculated the temporal variation of tidal shear and normal stresses on the focal fault plane of the M w7.8 Kahramanmaras, Türkiye earthquake. For earthquakes (M ≥ 4.0) occurring along and around the East Anatolian fault zone (EAFZ) and the Dead Sea fault zone (DSFZ) between January 2002 and December 2022, we analysed tidal stress-related seismicity from two viewpoints: occurrence time and seismic strain. Seismicity was found to be significantly correlated with semi-diurnal tidal shear stress (TSS) prior to the Kahramanmaras M w7.8 main shock. An area with significant correlation was ∼140 km northeast of the epicentre. For three and more years prior to the Kahramanmaras event, the temporal rate of seismic strain increased during periods of increasing TSS and did not change much during periods of decreasing TSS. From the above results, the increase in TSS may have facilitated the occurrence of earthquakes in the years preceding the Kahramanmaras earthquake. The reason for these tide-triggered earthquakes may be that the seismic region is loaded to a critical state. This observation may shed light on the seismogenic processes that led to the Kahramanmaras earthquake and its precursors.

Timescales of glacial isostatic adjustment in Greenland: is transient rheology required?

SUMMARY The possibility of a transient rheological response to ice age loading, first discussed in the literature of the 1980s, has received renewed attention. Transient behaviour across centennial to millennial timescales has been invoked to reconcile apparently contradictory inferences of steady-state (Maxwell) viscosity based on two distinct data sets from Greenland: Holocene sea-level curves and Global Navigation Satellite System (GNSS) derived modern crustal uplift data. To revisit this issue, we first compute depth-dependent Fréchet kernels using 1-D Maxwell viscoelastic Earth models and demonstrate that the mantle resolving power of the two Greenland data sets is highly distinct, reflecting the differing spatial scale of the associated surface loading: the sea-level records are sensitive to viscosity structure across the entire upper mantle while uplift rates associated with post-1000 CE fluctuations of the Greenland Ice Sheet have a dominant sensitivity to shallow asthenosphere viscosity. Guided by these results, we present forward models which demonstrate that a moderate low viscosity zone beneath the lithosphere in Maxwell Earth models provides a simple route to simultaneously reconciling both data sets by significantly increasing predictions of present-day uplift rates in Greenland whilst having negligible impact on predictions of Holocene relative sea-level curves from the region. Our analysis does not rule out the possibility of transient deformation, but it suggests that it is not required to simultaneously explain these two data sets. A definitive demonstration of transient behaviour requires that one account for the resolving power of the data sets in modelling the glacial isostatic adjustment process.

Navigation by magnetic signatures in a realistic model of Earth's magnetic field.

Certain animal species use the Earth's magnetic field (i.e. magnetoreception) alongside their other sensory modalities to navigate long distances that include continents and oceans. It is hypothesized that several animals use geomagnetic parameters, such as field intensity and inclination, to recognize specific locations or regions, potentially enabling migration without a pre-surveyed map. However, it is unknown how animals use geomagnetic information to generate guidance commands, or where in the world this type of strategy would maximize an animal's fitness. While animal experiments have been invaluable in advancing this area, the phenomenon is difficult to studyin vivoorin situ, especially on the global scale where the spatial layout of the geomagnetic field is not constant. Alongside empirical animal experiments, mathematical modeling and simulation are complementary tools that can be used to investigate animal navigation on a global scale, providing insights that can be informative across a number of species. In this study, we present a model in which a simulated animal (i.e. agent) navigates via an algorithm which determines travel heading based on local and goal magnetic signatures (here, combinations of geomagnetic intensity and inclination) in a realistic model of Earth's magnetic field. By varying parameters of the navigation algorithm, different regions of the world can be made more or less reliable to navigate. We present a mathematical analysis of the system. Our results show that certain regions can be navigated effectively using this strategy when these parameters are properly tuned, while other regions may require more complex navigational strategies. In a real animal, parameters such as these could be tuned by evolution for successful navigation in the animal's natural range. These results could also help with developing engineered navigation systems that are less reliant on satellite-based methods.

Horizontal transverse isotropy anisotropic parameter inversion via azimuthal seismic velocity anisotropy and its application to anisotropic 3D in-situ stress estimation

Anisotropic parameter inversion in horizontal transverse isotropy (HTI) media plays an important role in predicting the fracture density and the anisotropic in-situ stress for unconventional reservoirs. The current industry practice is to use the azimuthal PP-wave reflection coefficient to estimate the HTI anisotropic parameters. Based on the linear slip theory, we develop an innovative approach that uses azimuthal P-wave phase velocity to calculate HTI anisotropic parameters, which presents superiority over the conventional azimuthal PP-wave reflection coefficient inversion. Specifically, we first verify that the azimuthal P-wave phase velocity is for the HTI elliptical fitting than the azimuthal PP-wave reflection coefficient using the analytical formulations. Second, we sort the prestack wide-azimuth data into offset vector tile sectors and perform amplitude-variation-with-offset inversion at each azimuth. Third, we apply an elliptical fitting to the obtained azimuthal P-wave phase velocities to estimate the HTI anisotropic parameters, fracture density, and fracture direction. Fourth, based on the HTI mechanical earth model, we formulate a cost-effective 3D in-situ stress estimation method using the obtained elastic parameters and fracture compliance. Finally, field examples from the Zhaotong area, China, demonstrate that the estimated fracture density and anisotropic in-situ stress are more accurate and have higher resolution compared with conventional methods. The dominant stress regime in the study area is a strike-slip faulting regime with a governing orientation of northeast–southwest and presents good alignments with well logs, which demonstrates the reliability and accuracy of our method for predicting fracture density and anisotropic in-situ stress.

A Digital Twin of the terrestrial water cycle: a glimpse into the future through high-resolution Earth observations

Climate change is profoundly affecting the global water cycle, increasing the likelihood and severity of extreme water-related events. Better decision-support systems are vital to accurately predict and monitor water-related environmental disasters and optimally manage water resources. These must integrate advances in remote sensing, in situ, and citizen observations with high-resolution Earth system modeling, artificial intelligence (AI), information and communication technologies, and high-performance computing. Digital Twin Earth (DTE) models are a ground-breaking solution offering digital replicas to monitor and simulate Earth processes with unprecedented spatiotemporal resolution. Advances in Earth observation (EO) satellite technology are pivotal, and here we provide a roadmap for the exploitation of these methods in a DTE for hydrology. The 4-dimensional DTE Hydrology datacube now fuses high-resolution EO data and advanced modeling of soil moisture, precipitation, evaporation, and river discharge, and here we report the latest validation data in the Mediterranean Basin. This system can now be explored to forecast flooding and landslides and to manage irrigation for precision agriculture. Large-scale implementation of such methods will require further advances to assess high-resolution products across different regions and climates; create and integrate compatible multidimensional datacubes, EO data retrieval algorithms, and models that are suitable across multiple scales; manage uncertainty both in EO data and models; enhance computational capacity via an interoperable, cloud-based processing environment embodying open data principles; and harness AI/machine learning. We outline how various planned satellite missions will further facilitate a DTE for hydrology toward global benefit if the scientific and technological challenges we identify are addressed.

Global marine gravity gradient tensor inverted from altimetry-derived deflections of the vertical: CUGB2023GRAD

Abstract. Geodetic applications of altimetry have largely been inversions of gravity anomaly. Previous studies of Earth's gravity gradient tensor mostly presented only the vertical gravity gradient (VGG). However, there are six unique signals that constitute the gravity gradient tensor. Gravity gradients are signals suitable for detecting short-wavelength topographic and tectonic features. They are derived from double differentiation of the disturbing potential and hence are susceptible to noise amplification which was exacerbated by low across-track resolution of altimetry data in the past. However, current generation of altimetry observations have improved spatial resolutions, with some better than 5 km. Therefore, this study takes advantage of current high-resolution altimetry datasets to present CUGB2023GRAD, a global (latitudinal limits of ±80°) 1 arcmin model of Earth's gravity gradient tensor over the oceans using deflections of the vertical as inputs in the wavenumber domain. The results are first assessed via Laplace's equation, whereby the resultant residual gradient is virtually zero everywhere. Further analysis at local regions in the Arctic and south Indian oceans showed that Txy, Txz and Tyz are the most dominant gravity gradients for bathymetric studies. This proves that bathymetric signatures in the non-diagonal tensor components are worth exploiting. Bathymetric coherence analysis of Tzz over the Tonga Trench showed strong correlation with multibeam shipboard depths. This study proves that current generation of altimetry geodetic missions can effectively resolve Earth's gravity gradient tensor. The CUGB2023GRAD model data can be freely accessed at https://doi.org/10.5281/zenodo.10511125 (Annan et al., 2024).

Full-waveform tomography reveals iron spin crossover in Earth’s lower mantle

Three-dimensional models of Earth’s seismic structure can be used to identify temperature-dependent phenomena, including mineralogical phase and spin transformations, that are obscured in 1-D spherical averages. Full-waveform tomography maps seismic wave-speeds inside the Earth in three dimensions, at a higher resolution than classical methods. By providing absolute wave speeds (rather than perturbations) and simultaneously constraining bulk and shear wave speeds over the same frequency range, it becomes feasible to distinguish variations in temperature from changes in composition or spin state. We present a quantitative joint interpretation of bulk and shear wave speeds in the lower mantle, using a recently published full-waveform tomography model. At all depths the diversity of wave speeds cannot be explained by an isochemical mantle. Between 1000 and 2500 km depth, hypothetical mantle models containing an electronic spin crossover in ferropericlase provide a significantly better fit to the wave-speed distributions, as well as more realistic temperatures and silica contents, than models without a spin crossover. Below 2500 km, wave speed distributions are explained by an enrichment in silica towards the core-mantle boundary. This silica enrichment may represent the fractionated remains of an ancient basal magma ocean.

3D Stresses and Velocities Caused by Continental Plateaus: Scaling Analysis and Numerical Calculations With Application to the Tibetan Plateau

AbstractUnderstanding stresses is crucial for geodynamics since they govern rock deformation and metamorphic reactions. However, the magnitudes and distribution of crustal stresses are still uncertain. Here, we use a 3D numerical model in spherical coordinates to investigate stresses and velocities that result from lateral crustal thickness variations around continental plateaus like those observed for the Tibetan plateau. We do not consider any far‐field deformation so that the plateau deforms by horizontal dilatation and vertical thinning. We assume viscous creep, a simplified plateau geometry, and simplified viscosity and density distributions to couple the numerical results with a scaling analysis. Specifically, we study the impact of the viscosity ratio between crust and lithospheric mantle, a rectangular plateau corner, a stress‐dependent power‐law flow law and Earth's double curvature on the crustal stress field and horizontal velocities. Two orders of magnitude variation in crustal and lithospheric mantle viscosities change the maximum crustal differential stress only by a factor of ≈2. We derive simple analytical estimates for the crustal deviatoric stress and horizontal velocity which agree to first order with 3D numerical calculations. We apply these estimates to calculate the average crustal viscosity in the eastern Tibetan plateau as ≈1022 Pa · s. Furthermore, our results show that a corner strongly affects the stress distribution, particularly the shear stresses, while Earth's curvature has a minor impact on the stresses. We further discuss potential implications of our results to strike‐slip faulting and fast exhumation around the Tibetan plateau's syntaxes.

Sensitivity of GNSS‐Derived Estimates of Terrestrial Water Storage to Assumed Earth Structure

AbstractGeodetic methods can monitor changes in terrestrial water storage (TWS) across large regions in near real‐time. Here, we investigate the effect of assumed Earth structure on TWS estimates derived from Global Navigation Satellite System (GNSS) displacement time series. Through a series of synthetic tests, we systematically explore how the spatial wavelength of water load affects the error of TWS estimates. Large loads (e.g., >1,000 km) are well recovered regardless of the assumed Earth model. For small loads (e.g., <10 km), however, errors can exceed 75% when an incorrect model for the Earth is chosen. As a case study, we consider the sensitivity of seasonal TWS estimates within mountainous watersheds of the western U.S., finding estimates that differ by over 13% for a collection of common global and regional structural models. Errors in the recovered water load generally scale with the total weight of the load; thus, long‐term changes in storage can produce significant uplift (subsidence), enhancing errors. We demonstrate that regions experiencing systematic and large‐scale variations in water storage, such as the Greenland ice sheet, exhibit significant differences in predicted displacement (over 20 mm) depending on the choice of Earth model. Since the discrepancies exceed GNSS observational precision, an appropriate Earth model must be adopted when inverting GNSS observations for mass changes in these regions. Furthermore, regions with large‐scale mass changes that can be quantified using independent data (e.g., altimetry, gravity) present opportunities to use geodetic observations to refine structural properties of seismologically derived models for the Earth's interior structure.

Development and validation of a general-purpose ReaxFF reactive force field for earth material modeling.

ReaxFF reactive force field bridges the gap between nonreactive molecular simulations and quantum mechanical calculations and has been widely applied during the past two decades. However, its application to earth materials, especially those under high T-P conditions relevant to Earth's interior, is still limited due to the lack of available parameters. Here, we present the development and validation of a ReaxFF force field containing several of the most common elements in Earth's crust, i.e., Si/Al/O/H/Na/K. The force field was trained against a large data set obtained from density functional theory (DFT) calculations, including charges, bond/angle distortion curves, equation of states, ion migration energy profiles, and condensation reaction energies. Different coordination environments were considered in the training set. The fitting results showed that the current force field can well reproduce the DFT data (the Pearson correlation coefficient, Rp, is 0.95). We validated the force field on mineral-water interfaces, hydrous melts/supercritical geofluids, and bulk crystals. It was found that the current force field performed excellently in predicting the structural, thermodynamic, and transport properties of various systems (Rp = 0.95). Moreover, possible applications and future development have been discussed. The results obtained in this study suggest that the current force field holds good promise to model a wide range of processes and thus open opportunities to advance the application of ReaxFF in earth material modeling.