Geophysical Journal International

SummaryInduced polarization (IP) effects in airborne electromagnetic (AEM) surveys have commonly been investigated in helicopter-borne systems, leaving both a bibliographic and application gap for fixed-wing configurations. This gap partly reflects the large relative number of helicopter compared to fixed wing AEM systems, but also the geometric complexity of fixed wing platforms. In these platforms, nine geometric parameters come into play: the pitch, roll, and yaw of both transmitter and receiver, plus the three-axis offsets between the coils. Shifts in these factors can distort the measured data in ways that aren’t uniquely attributable, making it hard to pinpoint whether negative recordings truly arise from IP or from geometry-related effects. The non-fixed geometry also complicates removal of the primary field, often requiring iterative processing steps that may suppress or alter spectral content linked to IP. With advances in airborne IP understanding from helicopter-borne systems, revisiting fixed-wing platforms is both timely and necessary. Part A of this two-part study addresses this issue using the TEMPEST™ fixed-wing system connecting numerical modelling with field evidence. A suite of synthetic two-layer models with variable resistivity and chargeability parameters was developed to evaluate the system’s sensitivity to polarizable structures. The experiments demonstrate that IP effects, including negative secondary field responses, can be reliably detected in fixed-wing AEM data, both in X and Z magnetic field components. The capacity of these systems to detect IP phenomena is, however, strongly dependent on the electrical conductance of the environment. For instance, both fixed-wing and helicopter-borne systems, elevated near-surface conductance enhances the amplitude of purely electromagnetic induction currents, which in turn can dominate the recorded response and obscure the comparatively weaker polarization currents. More in general, IP detectability depends on the strength of the EM response generated by induction currents flowing elsewhere, which can dominate the small reverse current flow from a polarizable target. This highlights the critical role of near-surface conductivity in controlling the expression of IP responses and underscores the need to carefully account for these factors when interpreting survey data. The synthetic results are then connected with field-scale observations from a subset of the AusAEM dataset, over 470 000 line-km of TEMPESTTM data, where negative responses align with areas of low shallow conductance, confirming the simulation results. These finding open the way to the Part B of this study, where TEMPESTTM data are inverted taking into account IP and compared with helicopter-borne results and geological information.

SummaryAcoustic signals can couple to the ground, providing an opportunity to use seismic stations to investigate airborne sources. The study of Bishop et al. (2022) used wavefield simulations in a fluid-solid medium to quantify the role of topography on the seismic (ground) recordings of a monopole source in the air. We build upon this study by linking wavefield forward modeling with the source estimation code MTUQ, which can accommodate point forces or moment tensors in a solid medium, as well as sources in the air (or water) if they are enabled by the forward-modeling solver. We perform a series of synthetic numerical experiments to demonstrate that a dipole airborne source can be estimated using ground-based receivers, even within the presence of realistic topography. We investigate the influence of receiver coverage, topography, and assumed source location on the estimated results. The established capabilities raise the prospects for future efforts to estimate dipole sources in 3D models that include heterogeneity in the air and the earth in addition to topography.

AbstractSeismic simulation is fundamental for understanding earthquake physics and mitigating seismic hazards, but accurate seismic modeling requires fine computational grids, imposing severe memory and computational challenges. Traditional modeling solvers, relying on single-precision floating-point 32-bit (FP32) and scalar register-based computation, suffer from excessive memory consumption, low memory access efficiency, and limited computational efficiency. Compared with FP32, half-precision floating-point 16-bit (FP16) reduces memory consumption by 50% and improves memory access efficiency; relative to scalar registers, ARM’s Scalable Vector Extension (SVE) registers provide vectorized single-instruction multiple-data (SIMD) capabilities, significantly accelerating computation. However, leveraging the advantages of FP16 and SVE involves challenges such as FP16 overflow/underflow, SVE stencil adaptation, and SVE data misalignment from FP16 storage with FP32 arithmetic. Therefore, this study proposes three approaches on the ARM architecture: FP16-based, SVE-accelerated, and FP16-SVE hybrid; each is designed to tackle the respective challenges while exploiting FP16 memory efficiency and SVE computational acceleration. Correspondingly, the three solvers are implemented, validated, and benchmarked on both hypothetical models and real-world earthquake scenarios. The results of these solvers show near-perfect agreement with the reference solver, confirming their accuracy across diverse seismic scenarios. Moreover, the FP16-SVE hybrid solver halved memory usage and achieved up to 3× computational speedup, delivering more than 2.3× acceleration in the real-world earthquake simulation. The gains in high efficiency of memory and computation highlight the capability of the FP16-SVE hybrid solver to support large-scale, real-time seismic simulations and efficient earthquake hazard assessment on ARM platforms.

AbstractSurface waves such as Rayleigh, Love and Scholte waves can exhibit dispersion, i.e., variations in phase velocity with wavelength as a function of frequency. This property enables the inversion of 1D models of seismic velocity and density in the subsurface. Conventional deterministic and stochastic inversion schemes are widely applied to surface wave data but face two main challenges. The first is the identification of dispersion curves for fundamental and higher modes on wavefield-transformed images, which is often done manually. The second is the quantification of uncertainty, which can be computationally expensive in stochastic approaches or limited to data-propagated uncertainty in deterministic inversions. Our objectives are to (1) eliminate the need for manual or automatic dispersion curve picking, and (2) directly infer ensembles of 1D velocity models - and their associated uncertainties - from the full velocity spectrum, i.e., the complete dispersion image containing all modes. To this end, we employ Bayesian Evidential Learning, a predictive framework that reproduces experimental data from prior information while allowing prior falsification. In our application, ensembles of prior Earth models are sampled to predict 1D subsurface structures in terms of seismic velocity and, where applicable, attenuation from near-surface seismic wave data. This approach bypasses traditional inversion schemes and provides a computationally efficient tool for uncertainty quantification.

SummaryRepeating earthquakes are believed to result from recurring ruptures of a single asperity, driven by surrounding aseismic creep. However, their occurrence and behavior in intraplate regions remain poorly understood. This study investigates the repeating earthquakes in the Gyeryongsan region of the Korean Peninsula, a tectonically stable intraplate region, following the 2008 Mw 3.56 earthquake. We augmented the earthquake catalogue from 2007 to 2022 using template matching and identified one repeating earthquake family comprising ten events with irregular recurrence intervals. The repeating earthquakes, with a median magnitude of Mw 1.22, occurred within the rupture area of the Mw 3.56 mainshock, beginning in late 2010 and subsequently recurring intermittently between 2011 and 2019. Stress drops of nearby earthquakes increased gradually from 0.3-0.9 MPa to 8.6 MPa over a decade, indicating a fault strength recovery period substantially longer than that typically observed at plate boundaries. We interpret that the earthquakes occurred within a damaged fault zone, reflecting extremely low loading rates in the intraplate region. Our study provides insights into earthquake behaviour within intraplate damaged fault zones and documents a rare case of a repeating earthquake family that persisted over ∼12 years.

SummaryThe Palaeomagnetic Intensity (PINT) database documents variations in the full-vector of the ancient geomagnetic field that can be used to provide insights into the operation and evolution of the geodynamo. In this study, we report an update of PINT and the evolving behaviour of the palaeomagnetic field since 17 Ma. The update is the addition of 206 recently published site-mean data with ages between 0.06 and 2610 Ma that have been assessed using the palaeointensity quality criteria (QPI). Using this database, we analysed, for the first time, the distribution of values of the palaeosecular variation index (PSVi) in intervals drawn from the past 17 million years. Our results indicate that this index was enhanced prior to 5 Ma reflecting both lower average virtual dipole moments and higher angular deviations of the virtual geomagnetic pole (VGP) from the geographic pole. The present Brunhes chron is highlighted as being associated with especially high measurements of dipole moment which we hypothesise may be related to its already long duration relative to most other chrons of the last 17 Myr.

SummaryThe icy parts of the Earth, known as the cryosphere, are an integral part of the climate system. Comprehensively understanding the cryosphere requires dense observations, not only of its surface, but also of its internal structure and dynamics. Seismic methods play a central role in this endeavour. Fibre-optic sensing is emerging as a valuable complement and alternative to well-established inertial seismometers. Offering metre-scale channel spacing, interrogation distances of up to ∼100 km, and a bandwidth from mHz to kHz, it has enabled new seismological applications, for instance, under water, in cities and on volcanoes. Cryosphere research particularly benefits from fibre-optic sensing because long cables can be deployed with relative ease in icy environments where dense arrays of seismometers are difficult to install, including glaciers, ice sheets and deep boreholes. Intended to facilitate future fibre-optic seismology research in the cryosphere, this Expository Review combines a classical publication review with theoretical background, a practical field guide, a cryospheric signal gallery, and open-access data examples for hands-on training. Following a summary of recent findings about firn and ice structure, glacial seismicity, hydrology and avalanche dynamics, we derive the ideal instrument response of a distributed fibre-optic deformation sensor. To approach this ideal in field experiments, we propose numerous practical dos and don’ts concerning the choice and handling of fibre-optic cables, required equipment, splicing in the field at low temperatures, cable layout and trenching, and the deployment and coupling of cables in boreholes. A cryospheric signal gallery provides examples of data from a wide range of sources, such as explosions, land and air traffic, electricity generators, basal stick-slip icequakes, surface crevassing, englacial icequake cascades, floating ice shelf resonance, surface water flow and snow avalanches. Many of these data are enclosed as an open-access training resource, together with code for reading, visualisation and simple analyses. This review concludes with a discussion of grand open challenges in our understanding of cryosphere structure and dynamics, and how further advances in fibre-optic sensing may help to overcome them.

SummaryThis study proposes a dispersion-spectrum inversion method for improved estimation of shear-velocity (VS) profiles in marine sediments using underwater multichannel analysis of surface waves (UMASW). The method leverages an efficient forward modeling algorithm combined with a Markov chain Monte Carlo (McMC) global search to address the nonlinearity inherent in the inversion process. Comparisons with field and synthetic data demonstrate that the VS profiles inverted using the full dispersion spectrum (in terms of the frequency-phase velocity spectrum, FVS) exhibit greater stability and reliability than those obtained through traditional fundamental-mode (FM) inversion. Pseudo two-dimensional VS profiles are constructed by interpolating one-dimensional profiles obtained from the FVS-McMC inversion, showing more continuous subsurface interfaces that align with borehole data. Parametric studies further highlight the influence of Poisson’s ratio, water layer thickness, and VS contrast on the dispersion behavior, underscoring the robustness of the proposed approach for offshore site characterization.

SummaryAccurately mapping the distribution of natural methane hydrates is crucial for understanding their role in climate change and predicting the risks associated with hydrate dissociation. Attenuation shows great potential for remote hydrate detection, yet its behavior and underlying mechanisms are still not well understood. We conducted laboratory experiments to synthesize high-saturation methane hydrate in unconsolidated sands and measure attenuation based on ultrasonic waveforms. The resulting attenuation showed an unexpected decreasing trend during hydrate formation, contradicting previous studies in sands, where attenuation generally increases with hydrate saturation. Theoretical modeling suggests that attenuation is jointly controlled by hydrate and free gas. The gas reduction in pores due to hydrate formation substantially suppresses the attenuation induced by gas-bubble oscillation, and is therefore thought to be responsible for the observed attenuation reduction. By comparison, hydrate effects are relatively weak and strongly frequency-dependent. The discrepancy between our results and previous studies arises primarily from the distinct attenuation behavior across different ranges of gas content. Our samples fall within a relatively low gas content range, where attenuation is particularly sensitive to gas, highlighting its impact. These findings contribute new insights into the attenuation characteristics and mechanisms due to the coexistence of hydrate and gas in sediments.

SummaryDistributed Acoustic Sensing (DAS), a photonic technology that converts a fibre-optic cable into a long (tens of kilometres) high-linear-density (every few metres) array of seismo-acoustic sensors, can provide high-density, high-resolution strain measurements along the entire cable. The potential of such a distributed measurement has gained increasing attention in the seismology community for a wide range of applications. It has been shown that DAS has a sub-wavelength sensitivity to heterogeneities near the fibre-optic cable. This sensitivity is linked to the fact that the DAS measures deformation, as opposed to the displacements that seismometers measure. However, this sensitivity can create difficulties for many DAS applications, such as source location or distant imaging. Regardless, it can be advantageous in obtaining information about the subsurface near the cable. Here we present a method to locate small heterogeneities near the fibre-optic cable by inverting an indicator of the small-scale heterogeneities: the homogenised first-order corrector. We show that this first-order corrector can be used to locate heterogeneities near the fibre-optic cable at the gauge length precision, independent of the wavelength.

AbstractWe present an earthquake catalog in northeastern Tibetan plateau between September 2013 to April 2016 during the ChinArray-II deployment. Using continuous records from 676 transportable ChinArray-II stations and 172 permanent stations, the P/S phases are obtained using one deep learning phase picker. After associating these phases, the events are identified and located to establish the ChinArray-II Regional Earthquake Database (CARED-II). Benefiting from both improved station coverage and sensitive phase picker, CARED-II catalog has 156 057 events (around 3 million picks), about tenfold more than the manual routine catalog (15 967 events) using the permanent stations. The improved event catalog delineates the fault structures clearly. The deep structure of south-dipping north Qilian thrust faults is revealed, consisting with previous geology studies. The hidden faults and fault connectivity are revealed by improved seismicity, especially in the Alxa Block with sparse permanent stations and severe environments restricting geological field work. Moreover, small anthropogenic events are identified and related to highway tunnel construction across Qinling Mountain, forming a straight event cluster. The results demonstrate the high event detection ability of our procedure and reliability of the automatic catalog. Our array-based CARED-II catalog provides improved seismicity images in northeastern Tibet and could be used for further seismology and geotectonic studies.

SummaryThe global prevalence of organic pollutants presents a significant environmental challenge, necessitating sustainable remediation strategies. In situ biodegradation emerges as a cost-effective and eco-friendly solution. However, the real-time monitoring of in situ bacterial activities, particularly biodegradation processes, remains a challenge due to the limitations of traditional intrusive methods, including issues of representativeness, reproducibility, and high associated costs. Spectral induced polarization (SIP) has shown sensitivity to surface changes in subsurface environments, especially for biogeochemical reactivity monitoring including those associated with biodegradation. Despite this potential, advances have to be made to quantitatively link SIP parameters to in situ biodegradation processes. This study addresses this gap by conducting controlled biogeophysical experiments on a sand-packed column undergoing biodegradation facilitated by Rhodococcus wratislaviensis IFP 2006. SIP measurements were paired with bacterial growth kinetics to develop a quantitative model estimating bacterial growth. The results demonstrate that SIP, coupled with routine laboratory measurements, can effectively and quantitatively assess bacterial growth and the biodegradation of organic pollutants. These findings highlight the potential of SIP as a non-intrusive and reliable method for monitoring biodegradation in contaminated subsurface environments.

AbstractThe Surface Water and Ocean Topography mission (SWOT), equipped with the Ka-band Radar Interferometer (KaRIn), provides groundbreaking two-dimensional sea surface heights (SSHs), bringing new potential for optimizing the deflection of the vertical (DOVs). However, conventional DOV modeling—combining along- and cross-track geoid gradients with equal weights—fail to fully exploit the potential of SWOT/KaRIn observations and overlook the spatial variability in precision. We present a tailored method for optimizing DOVs estimation. The method combines geoid gradients in the along-track, cross-track, diagonal (forward and backward) directions with adaptive weighting. The refined weights are employed to exploit the potential of each geoid gradient based on the relationship between the standard deviation of SSHs and the significant wave height. To mitigate data gaps, prior and locally averaged geoid gradients are incorporated in the gaps and overlapping regions. SWOT/KaRIn-derived DOVs and gravity anomalies from the science-phase observations are validated against shipborne gravity in the Philippine Sea. Results indicate that the DOV model derived by the tailored method—particularly by combining triple-directional (along, cross, and diagonally forward) geoid gradients with refined weights—achieves a 7.3% improvement in accuracy over the conventional method. The supplement of additional geoid gradients is critical for mitigating leakage errors caused by missing or reduced observations in the gap regions. Furthermore, the gravity anomaly model recovered from DOVs by stacking 17-cycle observations achieved an accuracy of 2.97 mGal, representing a 7.2% improvement over single-cycle observations. The clear advantages of SWOT/KaRIn observations are gradually emerging in marine gravity recovery.

SummaryThe estimation of topographic gravity field models has attracted significant interest in recent years due to its growing relevance in Earth sciences. In this study, we present a robust methodology for the computation and comprehensive validation of global, complete spherical Bouguer and isostatic gravity anomalies that are essential for accurately interpreting subsurface mass distributions therefore geological structures. We synthesize these crucial gravitational functionals by leveraging spherical harmonic coefficients from high-resolution global gravity field models and various topographic/topographic-isostatic gravity field models. Our findings underscore the critical role of comprehensive terrain corrections in deriving physically meaningful, complete Bouguer gravity fields. The calculated global anomalies demonstrate strong coherence with established benchmark datasets, such as the World Gravity Map 2012. Residual differences are primarily attributed to variations in input Digital Terrain Models. Comparisons with regional Bouguer datasets reveal systematic biases that are largely explained by differing terrain correction methodologies. After removing this effect, there is a high level of consistency between the calculated global and published regional datasets, highlighting the utility of our global solutions, particularly in regions with sparse terrestrial data. Furthermore, the globally computed isostatic gravity anomalies exhibit significant agreement with both external global and diverse regional datasets, notably without the large systematic biases observed in Bouguer comparisons. This agreement reflects the effectiveness of the combined topographic and isostatic corrections in capturing Earth’s mass balance. This research provides valuable tools for new studies in the geoscience community by offering globally consistent and complete Bouguer and isostatic gravity field anomalies that have been rigorously validated for the ICGEM service.

SummaryDestabilization of volcanic edifices can generate debris avalanches with catastrophic impacts on their environment. We present the first high-resolution muography of Mount Unzen, Japan, conducted to characterize the structure of lava lobes formed on the volcano’s summit and flank during the 1990-1995 eruption. A multi-wire-proportional-chamber-based muon tracking system was operated for 203 days. The obtained high-resolution muographic image shows the internal density structure of Mount Unzen with a spatial resolution of 12 meters. Mean densities were respectively measured as 2,470 kg m−3 and 2,290 kg m−3 for the base rock and a fracture zone, and both were consistent with the results of prior drilling and sampling experiments. The mean density of lava lobes was measured significantly lower value of 1,570 kg m−3, indicating post-eruptive structural weakening. A comparison between the time-series of muographically measured density-lengths and daily precipitation records suggest that rainfall-induced gravitational destabilization did not occur during the observational period. This work demonstrates that long-term (multi-year) muon monitoring of the lava lobes can provide valuable complementary information for volcanic stability assessments.

SummaryThe asthenosphere is a weak layer in the upper mantle that supports the movement of the overriding tectonic plates and facilitates mantle convection. In this study, we compile a global dataset of SS precursors reflected at the base of the asthenosphere, also known as the 220-km discontinuity. The global dataset includes the oceanic SS precursors from Sun & Zhou (2025a) and new measurements with bounce points in continental regions. Similar to the oceanic dataset, the continental SS precursors are observed on about 45% of the SS waves, with bounce points distributed across all tectonic regions — from orogeny belts to stable cratons. We image the depth of the discontinuity at a global scale using finite-frequency tomography. In oceanic regions, the depth of the 220-km discontinuity agree well with the previous study, with discontinuity depth structure characterized by alternating linear bands of shallow and deep anomalies that roughly follow seafloor age contours. In continental regions, the variations are not spatially oscillatory but are instead much broader, with prominent perturbations associated with convergent plate boundaries. The base of the asthenosphere is shallow along the southern border of the Eurasian plate, from the Mediterranean region to Southeast Asia. Shallow discontinuity anomalies are also observed in the continental interiors – in Eurasia, from the northern Tian Shan through Mongolia to eastern Siberia, and in North America east of the Rocky Mountains. These anomalies form a linear structure roughly parallel to the Pacific subduction zones. The average depth of the discontinuity, as well as the velocity contrast across the interface, is globally consistent across both oceans and continents, with an average depth of approximately 251 km and a velocity increase of about 7%. Given that the continental lithosphere has been cooling for much longer than the oceanic lithosphere, the observed consistency in the average depth of the discontinuity implies that secular cooling does not significantly impact the thermal structure at the base of the asthenosphere.

SummaryA very frequent approach for studying lithospheric processes is to deploy temporary seismological networks in dedicated areas and to map the mantle structures with different approaches. One of them is the well-established relative travel time body wave tomography. Different circumstances often lead to a non-uniform deployment of stations both in space and time, and a wish to combine data which have been acquired asynchronously. This is the situation in Patagonia where two distinct seismic experiments provide complementary seismic data over the region covering the Patagonia slab window. Combining these data in one regional relative body wave tomography is however problematic as the two data sets are a priori with respect to two different reference models. In this contribution, we show that the number of finite-frequency relative travel time residuals varies very strongly from station to station for this data set, violating the assumption implicit in relative travel time tomography of a unique reference model due to an even data distribution for all events. We demonstrate the superiority of the inversion using relative sensitivity kernels compared with a traditional approach with absolute kernels and event terms. A resolution test proves how this is crucial for resolving the important issue of the eastern extent of the slab window. In addition, we discuss potential issues related to interference of the direct phases with core phases when measuring finite-frequency travel time residuals by cross-correlation of waveforms in necessarily relatively large time windows. We also briefly outline our preferred strategy for performing crustal correction, keeping in mind that finite-frequency residuals require frequency-dependent crustal corrections.

SummaryP-wave receiver functions (RFs), which utilize P-wave conversions to probe subsurface structures, face significant challenges in sedimentary environments. Specifically, strong reverberations generated by ultra-low-velocity sedimentary layers distort RF waveforms and obscure crustal signals, posing challenges for robust shallow crustal imaging. We develop a novel Bayesian joint inversion framework that simultaneously utilizes three complementary datasets—reverberant receiver functions, dereverberated receiver functions, and surface wave dispersion—to address this challenge. Our approach employs Unscented Kalman Inversion, a derivative-free method that efficiently handles nonlinear joint inversion problems. Synthetic tests demonstrate that our joint inversion recovers sediment thickness and Moho depth with uncertainties of ±0.50 km and ±1.0 km, respectively. Application to real data from the Songliao Basin verifies the approach, successfully reconstructing sediment thickness and Moho depth beneath sedimentary cover. This methodology demonstrates potential for advancing crustal investigations in complex sedimentary settings, such as continental rift basins and oceanic margins, where sedimentary sequences of variable thickness often obscure deeper structures.

SummaryThree-dimensional (3-D) forward modeling of magnetotelluric (MT) data remains a computationally challenging task, particularly when accurate broadband MT responses are simulated for real-world problems that often involve complex multi-scale bathymetry and/or topography. To overcome this challenge, we developed a new efficient numerical approach for 3-D MT forward modeling that combines high-order Nédélec-type finite elements and high-order meshes, allowing us to obtain superior accuracy and account for complex material boundaries and interfaces. Despite gains in accuracy, higher-order FE solvers are often considered impractical owing to higher memory consumption and a more ill-conditioned system. To overcome these limitations, we use an iterative solver accelerated by the Low-Order-Refined (LOR) preconditioner, which uses spectrally equivalent low-order operators, rendering the complexity independent of the polynomial degree. Another key novelty is a matrix-free implementation, where the action of the high-order operator is computed efficiently without explicit matrix assembly. The low-order system is solved using an Auxiliary Space Maxwell (AMS) solver based on a multigrid solver. We demonstrate the efficiency in a series of numerical experiments. Scalability analysis on a 3-D benchmark model demonstrates that the LOR preconditioner significantly outperforms the current state-of-the-art AMS preconditioner in terms of CPU time and memory usage, especially for higher polynomial degrees. Excellent scalability is confirmed by solving a problem with up to 1.5 × 109 degrees of freedom in less than 2 minutes using 16,384 CPU cores, which is, to our knowledge, the largest 3-D MT problem reported to date. We also illustrate that high-order hexahedral meshes allow for accurate discretization of complex geometries, such as topography, with substantially fewer elements than conventional linear meshes. Finally, the capability of the integrated approach is demonstrated on a real 3-D model crossing the ocean trench in the Aleutian subduction zone. The proposed methods pave the way for more efficient and accurate 3-D MT modeling that is crucial for the inversion of complex MT data sets.

SummaryAtmospheric models are based on various types of geophysical data, including lidar and radar. Infrasounds, acoustic waves that can propagate over large distances, have not yet been used in atmospheric models, although they provide valuable information. Besides their sensitivity to atmospheric phenomena such as gravity waves, infrasound also presents the advantage of being omnipresent. Previous studies explored the use of infrasound packet arrival properties for model estimation. However, properties such as arrival times present less information than full waveforms. We aim here to investigate, for the first time, the sensitivity of a full infrasound waveform to model parameters and to use these sensitivities in an inverse problem to recover atmospheric structure. For this purpose, infrasound propagation is modeled by Euler equations (i.e. Navier-Stokes equations in the absence of attenuation effects), and discretization is carried out here using the finite-differences method. Waveform sensitivity to atmospheric parameters is computed through the adjoint method via a novel and optimized double checkpointing-based procedure and validated by comparison with a small perturbation method. As an illustration, these sensitivity kernels are computed for the idealized case of an explosion in Finland, recorded by a CTBT station. These first results demonstrate the high sensitivity of infrasound waveforms to the atmospheric perturbations generated by gravity waves. Moreover, the sensitivity kernels of infrasound waveforms allow us to recover the variations of model parameters by solving an inverse problem. To demonstrate this capability, full waveform non-linear inversions are performed using the Limited Broyden-Fletcher-Goldfarb-Shanno method (L-BFGS): wind and sound speed profiles are inverted for a test case with idealized conditions and a synthetic dataset. These estimates of infrasound sensitivity kernels are closing a knowledge gap that allows the use of infrasound full waveforms to constrain atmospheric models.