ALBERT

Description: 〈span〉〈div〉Summary〈/div〉This study presents for the first time an integrated image of the crust and lithospheric mantle of Alaska and its adjacent western shelves of the Chukchi and Bering seas based on joint modeling of potential field data constrained by thermal analysis and seismic data. We also perform 3D forward modelling and inversion of Bouguer anomalies to analyze density heterogeneities at the crustal level. The obtained crustal model shows northwest-directed long wavelength thickening (32 to 36 km), with additional localized trends of thicker crust in the Brooks Range (40 km) and in the Alaska and St. Elias ranges (50 km). Offshore, 28–30 km thick crust is predicted near the Bearing slope break and 36–38 km in the northern Chukchi Shelf. In interior Alaska, the crustal thickness changes abruptly across the Denali fault, from 34–36 to the north to above 30 km to the south. This sharp crustal thickness gradient agrees with the presence of a crustal tectonic buttress guiding block motion west and south towards the subduction zone. The average crustal density is 2810 kg∙m〈sup〉−3〈/sup〉. The denser crust, up to 2910 kg∙m〈sup〉−3〈/sup〉, is found south of the Denali Fault likely related to the oceanic nature of the Wrangellia Composite Terrane rocks. Offshore, less dense crust (〈 2800 kg∙m〈sup〉−3〈/sup〉) is found along the sedimentary basins of the Chukchi and Beaufort shelves. At LAB levels, there is a regional SE–NW trend that coincides with the current Pacific plate motion, with a lithospheric root underneath the Brooks Range, Northern Slope, and Chuckchi Sea, that may correspond to a relic of the Chukotka-Artic Alaska microplate. The obtained lithospheric root (above 180 km) agrees with the presence of a boundary of cold, strong lithosphere that deflects the strain towards the South. South of the Denali Fault the LAB topography is quite complex. East of 150 ° W, below Wrangellia and the eastern side of Chugach terranes, the LAB is much shallower than it is west of this meridian. The NW trending limit separating thinner lithosphere in the east and thicker in the west agrees with the two–tiered slab shape of the subducting Pacific Plate.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Earthquake focal mechanisms put primary control on the distribution of ground motion, and also bear on the stress state of the crust. Most routine focal mechanism catalogs still use 1D velocity models in inversions, which may introduce large uncertainties in regions with strong lateral velocity heterogeneities. In this study, we develop an automated waveform-based inversion approach to determine the moment tensors of small-to-medium-sized earthquakes using 3D velocity models. We apply our approach in the Los Angeles region to produce a new moment tensor catalog with a completeness of M〈sub〉L 〈/sub〉≥ 3.5. The inversions using the Southern California Earthquake Center Community Velocity Model (3D CVM-S4.26) significantly reduces the moment tensor uncertainties, mainly owing to the accuracy of the 3D velocity model in predicting both the phases and the amplitudes of the observed seismograms. By comparing the full moment tensor solutions obtained using 1D and 3D velocity models, we show that the percentages of non-double-couple components decrease dramatically with the usage of 3D velocity model, suggesting that large fractions of non-double-couple components from 1D inversions are artifacts caused by unmodeled 3D velocity structures. The new catalog also features more accurate focal depths and moment magnitudes. Our highly accurate, efficient, and automatic inversion approach can be expanded in other regions, and can be easily implemented in near real-time system.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The far-field assumption is widely used and suitable for the moment-tensor inversion, in which the source-receiver distance is quite long. However, the description of far field is uncertain and an explicit far-field range is missing. In this study, the explicit far-field range is determined and the errors of moment-tensor solutions produced by the far-field approximation are analyzed. The distance, for which the far-field assumption is satisfied and the effect of the near-field term can be ignored, is directionally dependent. For the shear dislocation, in the directions near the nodal lines of the far-field P waves, the far-field distance is heavily dependent on the displacement component used to invert moment tensors. The radial component of displacement, which is parallel to the wave propagation direction, is recommended for the inversion and the corresponding far-field distance is quite short. In the directions far from the nodal lines, the selection of displacement components has little influence on the far-field distance. The maximum far-field distance appears in the directions of the pressure and tensile axes of the source and the value is about 30 wavelengths of radiated waves. Using more receivers (〉6) in the moment-tensor inversion can shorten the far-field distance. The effect of the near-field term on the moment-tensor inversion for tensile dislocations and isotropic sources (explosion or implosion) can be ignored. The conclusions obtained in this study are helpful for determining the positions of receivers and evaluating the accuracy of moment-tensor solutions, with far-field assumption being applied in the inversion.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We present 2-D attenuation images of the Mt. Etna volcanic region on the basis of separation of intrinsic and scattering effects. The analysis presented here exploits a large active seismic database that fully covers the area under study. We observe that scattering effects dominate over intrinsic attenuation, suggesting that the region is very heterogeneous. Comparison with analyses conducted at other volcanoes reveals that the Mt. Etna region is characterised by high intrinsic attenuation, resulting from the presence of large volcanoclastic deposits at shallow depth. The 2-D distributions of intrinsic and scattering anomalies show the presence of regions characterised by high and low attenuation effects, corresponding to several tectonic and volcanic features. In particular, we identify a high attenuation region in the SW sector of the Mt. Etna volcanic complex, which is correlated with high seismicity rates and volcanism. This work supports the hypothesis of a link between the dynamics of the SW flank and the recharge of the volcano in the last decades, occurring under the summit crater and, secondarily, the upper South rift zone.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Coda-Q is used to estimate the attenuation and scattering properties of the Earth (Aki & Chouet 1975). So far focus has been on earthquake data at frequencies above 1 Hz, as the high noise level in the first and second microseismic peak, and possibly lower scattering coefficient, hinder stable measurements at lower frequencies. In this work, we measure and map coda-Q in the period bands 2.5 s–5 s, 5 s–10 s and 10 s–20 s in the greater Alpine region using noise cross-correlations between station pairs, based on data from permanent seismic stations and from the temporary AlpArray experiment. The observed coda-Q for short interstation distances is independent of azimuth so there is no indication of influence of the directivity of the incoming noise field on our measurements. In the 2.5 s–5 s and 5 s–10 s period bands, our measurements are self-consistent, and we observe stable geographic patterns of low and high coda-Q in the period bands 2.5 s–5 s and 5 s–10 s. In the period band 10 s–20 s, the dispersion of our measurements increases and geographic patterns become speculative. The coda-Q maps show that major features are observed with high resolution, with a very good geographical resolution of for example low coda-Q in the Po Plain. There is a sharp contrast between the Po Plain and the Alps and Apennines where coda-Q is high, with the exception a small area in the Swiss Alps which may be contaminated by the low coda-Q of the Po Plain. The coda of the correlations is too short to make independent measurements at different times within the coda, so we cannot distinguish between intrinsic and scattering Q. Measurements on more severely selected datasets and longer timeseries result in identical geographical patterns but lower numerical values. Therefore, high coda-Q values may be overestimated, but the geographic distribution between high and low coda-Q areas is respected. Our results demonstrate that noise correlations are a promising tool for extending coda-Q measurements to frequencies lower than those analysed with earthquake data.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Salt structures are high potential targets for oil and gas exploration. However, large-scale salt domes with irregular surfaces pose significant challenges for velocity model building. For full waveform inversion, in the absence of a high-fidelity initial model, the success of the inversion depends on low-frequency seismic data, which are scarce in the exploration data sets. This paper presents a new idea to solve the problem of salt structure velocity modeling. Firstly, we propose an envelope-based full-band seismic data reconstruction algorithm. The smoothness of envelope is used to segment the events in seismic data, and the phase independence of envelope is used for the identification of the seismic event's arrival-time to obtain the apparent reflection sequences of the subsurface. Full-band seismic data are obtained by convolving the apparent reflection sequence with full-band source. Window averaging function and threshold strategy are used to ensure the accuracy of seismic event segmentation and the stability of the algorithm when dealing with noisy data. Then the multiscale reflection waveform inversion based on reconstructed data is proposed for salt structure velocity building. The numerical experiment results of the Sigbee2A model demonstrate the performance of the inversion algorithm in the case where the seismic data lack low-frequency components and contain noise. The limitations of the algorithm have also been analyzed and studied.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Teleseismic waveforms contain information on fault slip evolution during an earthquake, as well as on the fault geometry. A linear finite-fault inversion method is a tool for solving the slip-rate function distribution under an assumption of fault geometry as a single or multiple-fault-plane model. An inappropriate assumption of fault geometry would tend to distort the solution due to Green’s function modelling errors. We developed a new inversion method to extract information on fault geometry along with the slip-rate function from observed teleseismic waveforms. In this method, as in most previous studies, we assumed a flat fault plane, but we allowed arbitrary directions of slip not necessarily parallel to the assumed fault plane. More precisely, the method represents fault slip on the assumed fault by the superposition of five basis components of potency-density tensor, which can express arbitrary fault slip that occurs underground. We tested the developed method by applying it to real teleseismic 〈span〉P〈/span〉 waveforms of the 〈span〉M〈/span〉〈sub〉W〈/sub〉 7.7 2013 Balochistan, Pakistan, earthquake, which is thought to have occurred along a curved fault system. The obtained spatiotemporal distribution of potency-density tensors showed that the focal mechanism at each source knot was dominated by a strike-slip component with successive strike angle rotation from 205° to 240° as the rupture propagated unilaterally towards the south-west from the epicentre. This result is consistent with Earth’s surface deformation observed in optical satellite images. The success of the developed method is attributable to the fact that teleseismic body waves are not very sensitive to the spatial location of fault slip, whereas they are very sensitive to the direction of fault slip. The method may be a powerful tool to extract information on fault geometry along with the slip-rate function without requiring detailed assumptions about fault geometry.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉It is essential to pick 〈span〉P〈/span〉-wave and 〈span〉S〈/span〉-wave arrival times rapidly and accurately for the microseismic monitoring systems. Meanwhile, it is not easy to identify the arrivals at a true phase automatically using traditional picking method. This is one of the reasons that many researchers are trying to introduce deep neural networks to solve these problems. Convolutional neural networks (CNNs) are very attractive for designing automatic phase pickers especially after introducing the fundamental network structure from semantic segmentation field, which can give the probability outputs for every labelled phase at every sample in the recordings. The typical segmentation architecture consists of two main parts: (1) an encoder part trained to extracting coarse semantic features; (2) a decoder part responsible not only for recovering the input resolution at the output but also for obtaining sparse representation of the objects. The fundamental segmentation structure performs well; however, the influence of the parameters in the structure on the pickers has not been investigated. It means that the structure design just depends on experience and tests. In this paper, we solve two main questions to give some guidance on network design. First, we show what sparse features will learn from the three-component microseismic recordings using CNNs. Second, the influence of two key parameters in the network on pickers, namely, the depth of decoder and activation functions, is analysed. Increasing the number of levels for a certain layer in the decoder will increase the burden of demand on trainable parameters, but it is beneficial to the accuracy of the model. Reasonable depth of the decoder can balance prediction accuracy and the demand of labelled data, which is important for microseismic systems because manual labelling process will decrease the real-time performance in monitoring tasks. Standard rectified linear unit (ReLU) and leaky rectified linear unit (Leaky ReLU) with different negative slopes are compared for the analysis. Leaky ReLU with a small negative slope can improve the performance of a given model than ReLU activation function by keeping some information about the negative parts.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Accurate synthetic seismic wavefields can now be computed in 3-D earth models using the spectral element method (SEM), which helps improve resolution in full waveform global tomography. However, computational costs are still a challenge. These costs can be reduced by implementing a source stacking method, in which multiple earthquake sources are simultaneously triggered in only one teleseismic SEM simulation. One drawback of this approach is the perceived loss of resolution at depth, in particular because high-amplitude fundamental mode surface waves dominate the summed waveforms, without the possibility of windowing and weighting as in conventional waveform tomography.This can be addressed by redefining the cost-function and computing the cross-correlation wavefield between pairs of stations before each inversion iteration. While the Green’s function between the two stations is not reconstructed as well as in the case of ambient noise tomography, where sources are distributed more uniformly around the globe, this is not a drawback, since the same processing is applied to the 3-D synthetics and to the data, and the source parameters are known to a good approximation. By doing so, we can separate time windows with large energy arrivals corresponding to fundamental mode surface waves. This opens the possibility of designing a weighting scheme to bring out the contribution of overtones and body waves. It also makes it possible to balance the contributions of frequently sampled paths versus rarely sampled ones, as in more conventional tomography.Here we present the results of proof of concept testing of such an approach for a synthetic 3-component long period waveform data set (periods longer than 60 s), computed for 273 globally distributed events in a simple toy 3-D radially anisotropic upper mantle model which contains shear wave anomalies at different scales. We compare the results of inversion of 10 000 s long stacked time-series, starting from a 1-D model, using source stacked waveforms and station-pair cross-correlations of these stacked waveforms in the definition of the cost function. We compute the gradient and the Hessian using normal mode perturbation theory, which avoids the problem of cross-talk encountered when forming the gradient using an adjoint approach. We perform inversions with and without realistic noise added and show that the model can be recovered equally well using one or the other cost function.The proposed approach is computationally very efficient. While application to more realistic synthetic data sets is beyond the scope of this paper, as well as to real data, since that requires additional steps to account for such issues as missing data, we illustrate how this methodology can help inform first order questions such as model resolution in the presence of noise, and trade-offs between different physical parameters (anisotropy, attenuation, crustal structure, etc.) that would be computationally very costly to address adequately, when using conventional full waveform tomography based on single-event wavefield computations.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We present a method for measuring the dispersion of higher mode surface wave phase velocities from a single seismogram using a hierarchical transdimensional Bayesian approach. The 1-D shear velocity profiles down to 800km depth between sources and seismic stations are regarded as the controlling parameters to tune the phase velocities of fundamental and higher modes. The misfits between synthetics and real waveforms indicate whether the phase velocities are recovered well from the data. We use Monte Carlo Markov chains (MCMC) to approximate the posterior distribution of each model parameters, and assess the uncertainties from these probability density functions. These techniques can test models of varying dimensions while being parsimonious, thereby letting the data themselves control the complexity of the solution. Another advantage is that the algorithm can decide how much data noise is needed in order to explain the data without overfitting them. The data noise can be treated as an unknown and different noise levels can be applied to the different time windows considered. The posterior noise distributions can then be used as an indicator of the quality of the waveform fit within each frequency-time window. We considered phase velocities between 50s and 200s for each mode, and performed a reliability analysis to determine which modes and periods are reliably constrained. In this paper, we first present the method and demonstrate its feasibility with synthetic tests, which show that the technique is robust. We then illustrate it with applications to real data. We applied the method to two paths sampling Australia using earthquakes at regional distances, and obtained results that agree well with previous studies. The new method can be used in regional and global tomographic studies to obtain phase velocity maps and 3-D models of seismic velocities and anisotropy at depths that are not well resolved by fundamental mode surface waves or body waves.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉It is difficult to separate additive random noise from spatially coherent signal using a rank-reduction method that is based on the truncated singular value decomposition (TSVD) operation. This problem is due to the mixture of the signal and the noise subspaces after the TSVD operation. This drawback can be partially conquered using a damped rank reduction method, where the singular values corresponding to effective signals are adjusted via a carefully designed damping operator. The damping operator works most powerfully in the case of a small rank and a small damping factor. However, for complicated seismic data, e.g., multi-channel reflection seismic data containing highly curved events, the rank should be large enough to preserve the details in the data, which makes the damped rank reduction method less effective. In this paper, we develop an optimal damping strategy for adjusting the singular values when a large rank parameter is selected so that the estimated signal can best approximate the exact signal. We first weight the singular values using optimally calculated weights. The weights are theoretically derived by solving an optimization problem that minimizes the Frobenius-norm difference between the approximated signal components and the exact signal components. The damping operator is then derived based on the initial weighting operator to further reduce the residual noise after the optimal weighting. The resulted optimally damped rank reduction method is nearly an adaptive method, i.e., insensitive to the rank parameter. We demonstrate the performance of the proposed method on a group of synthetic and real five-dimensional seismic data.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We present a rheological model of the continental lithosphere during the rift–drift transition, founded on theoretical considerations, numerical modelling, and supported by geophysical and geological evidence from the Red Sea region. The model is based on the fundamental assumption that during the rifting phase the necking lithospheric mantle of the conjugate continental margins retains and accumulates elastic strain. This hypothesis is tested numerically, thereby it is shown that during the phase of extension the uppermost mantle has the capability to store and maintain recoverable elastic strain over geological times, differently from the upper crust, which exhibits a continuous short–term irreversible deformation by seismic release. After the onset of sea–floor spreading, strain recovery and release of the strain energy accumulated in the lithospheric mantle occurs through a phase of non–linear anelastic relaxation. During this phase, the upper crust of the conjugate continental margins experiences post–rift deformation with tectonic inversion of former extensional structures, while the extra–space created along the axial zone as a consequence of the rapid contraction of the margins triggers a rapid upwelling of asthenosphere that induces an initial pulse of fast spreading followed by a steady phase of oceanic crust accretion. We present geophysical evidence supporting this model, including: 1) The observed pattern of oceanic magnetic anomalies in the Red Sea, 2) The distribution of finite crustal strains across the continental margins of Nubia and Arabia, and 3) The distribution of earthquake epicentres along the western margin of the Arabian plate. We also present new structural data acquired during three geological campaigns performed in 2015 and 2016 along the western Arabian margin, which are consistent with a post–rift phase of compression and inversion of the rift structures. We will also show that a selection of realistic rheological parameters supports non–linear viscoelastic behaviour of the continental lithosphere during the rift–drift transition. Finally, we will show that in the case of the Red Sea ∼40 per cent of the total extensional strain accumulated during the rifting stage has been recovered in the southernmost part of the Arabian margin conjugate to the Nubian plate (∼19° N), while this percentage decreases to ∼14 per cent around 23.8° N, where the continental margin faces the youngest spreading segment, and it is zero north of this area, where the Red Sea is still in the rifting stage. Taking into account of the age of oceanization along the central and northern Red Sea this implies an average recovery of ∼10 per cent Myr〈sup〉–1〈/sup〉.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉In December, 2009, a rare sequence of earthquakes initiated within the weakly extended Western Rift of the East African Rift system in the Karonga province of northern Malawi, providing a unique opportunity to characterize active deformation associated with intrabasinal faults in an early-stage rift. We combine teleseismic and regional seismic recordings of the largest events, InSAR imagery of the primary sequence, and recordings of aftershocks from a temporary (4-month) local network of six seismometers to delineate the extent and geometry of faulting. The locations of ∼1900 aftershocks recorded between Jan-May 2010 are largely consistent with a west-dipping normal fault directly beneath Karonga as constrained by InSAR and CMT fault solutions. However, a substantial number of epicenters cluster in an east-dipping geometry in the central part of the study area, and additional west-dipping clusters can be discerned near the shore of Lake Malawi, particularly in the southern part of the study area. Given the extensive network of hanging wall faults mapped in the Karonga region on the surface and in seismic reflection images, the distribution of events is strongly suggestive of multiple faults interacting to produce the observed deformation, and the InSAR data permit this but do not require it. We propose that fault interaction contributed to the seismic moment release as a series of M〈sub〉w〈/sub〉 5 and 6 events instead of a normal mainshock-aftershock sequence. We find the depth of fault slip during the mainshocks constrained by InSAR peaks at less than 6 km, while the majority of recorded aftershocks are deeper than 6 km. This depth discrepancy appears to be robust and may be explained by fault interaction. Structural complexities associated with fault interaction may have limited the extent of co-seismic slip during the mainshocks, which increased stress deeper than the co-seismic slip zone on the primary fault and synthetic faults to the east, causing the energetic aftershock series. There is no evidence of deformation at the Rungwe volcanic province ∼50 km north of the earthquake sequence between 2007-2010, consistent with previous interpretations of no significant magmatic contribution during the sequence.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉In co-seismic slip distribution inversion, the Laplacian smoothness constraint is used to avoid the rank deficiency of the coefficient Green matrix and ensure the smoothness among the patches of the fault plane. However, by introducing the classic Laplacian smoothness constraint, the inversion of the maximum slip is commonly underestimated. To better solve this problem, here we propose a new method of weighting to the Laplacian second-order smoothness matrices by the slip. The Laplacian smoothness constraint with unequal weights and the classic Laplacian constraint are used to carry out the simulation experiments. The simulation experiments show that the accuracy of the estimated maximum slip from our method is outperforming the classical Laplacian method by a 12 per cent–19 per cent. Moreover, the estimations of the average slip and moment magnitude have been enhanced, which improved ranging from 4 per cent -12.5 and 0.4 per cent -9 per cent over the inversion results of classic Laplacian smoothness constraint, respectively. In order to further validate the general applicability of the proposed method, we conduct the experiments in terms of the inversion of L'Aquila and Taiwan Meinong earthquakes. The maximum slip inversion results show that the inversion of Laplacian smoothness constraint with unequal weights are larger than that of the classic Laplacian smoothness constraint, which is consistent with the conclusion of simulation experiments. In addition, the parameters inversion results of the actual earthquake from the proposed method are in agreement with other studies which also indicate the feasibility and effectiveness of our method.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Estimating subsurface density is important for imaging various geologic structures such as volcanic edifices, reservoirs, and aquifers. Muon tomography has recently been used to complement traditional gravity measurements as a powerful method for probing shallow subsurface density structure beneath volcanoes. Gravity and muon data have markedly different spatial sensitivities and, as a result, the combination is useful for imaging structures on spatial scales that are larger than the area encompassed by crossing muon trajectories. Here we explore and test a joint inversion of gravity and muon data in a study area where there is an independently-characterised target anomaly: a regionally extensive, high-density layer beneath Los Alamos, New Mexico, USA. We resolve the nearly flat-lying structure using a unique experimental set-up wherein surface and subsurface gravity and muon measurements are obtained above and below the target volume. Our results show that with minimal geologic (prior) constraints, the joint inversion correctly recovers salient features of the expected density structure. The results of our study illustrate the potential of combining surface and subsurface (e.g., borehole) gravity and muon measurements to invert for shallow geologic structures.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The age of the inner core nucleation is a first-order problem in the thermal evolution of the Earth that can be addressed with paleomagnetism. We conducted a paleointensity study on the 1.3 Ga Gardar basalts from southern Greenland to investigate previously reported high ancient geomagnetic field intensities. Biggin 〈span〉et al.〈/span〉 (2015) used the earlier result to identify nucleation of Earth's solid inner core at 1.3 Ga. We collected 106 samples from 39 flows from the lavas of the Eriksfjord Formation, sampling 17 of the lower flows, 8 of the middle flows and 14 of the upper flows. Rock magnetic analyses, including magnetic hysteresis, first order reversal curves (FORCs), and magnetic susceptibility versus temperature measurements, suggest that the predominate magnetic mineral in the lower basalts is low Ti titanomagnetite, whereas the middle and upper flows have varying amounts of hematite. The magnetic hysteresis data suggest magnetic grains range from multi-domain to single domain in character, with an apparent dominance of pseudo-single behavior. Thellier-Thellier double heating experiments using the IZZI methodology yielded vector endpoint diagrams and Arai plots showing two components of magnetization, one up to approximately 450˚C and the higher temperature component typically from 450˚C up to 580˚C, but sometimes to as high as 680˚C. We attribute the lower temperature component, to partial overprinting by the nearby Ilimaussaq intrusion, and acquisition of viscous remanent magnetization. We use the Thellier auto-interpreter (Shaar & Tauxe, 2013) assigning standard selection criteria vetted by cumulative distribution plots. This approach yields a paleointensity of 6.5 μT ± 5.9 μT (1 SD) based on 27 samples from 13 flows and a nominal virtual dipole moment (VDM) of 1.72 × 10〈sup〉22〈/sup〉 Am〈sup〉2〈/sup〉. However, we cannot exclude the possibility of bias in this value related to chemical remanent magnetization (CRM) and multi-domain effects. We isolate a conservative upper bound on paleointensity as the highest paleointensity result that is free of CRM effects. This yields a paleointensity of ∼18 μT, and a VDM of ∼4.5 × 10〈sup〉22〈/sup〉 Am〈sup〉2〈/sup〉, which is a field strength similar to many other Proterozoic values. Thus, our analysis of the Gardar basalts supports the conclusion of Smirnov 〈span〉et al.〈/span〉 (2016) that there is no paleointensity signature of inner core growth 1.3 billion years ago.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We present new shear wave anisotropy measurements in the continental crust along the Mexican subduction zone obtained from tectonic tremor. The new measurements were made in the states of Jalisco, Colima, Michoacán, and Oaxaca. To make a complete analysis of the anisotropic crustal structure, we also include previous measurements reported in Guerrero using tremor signals. Since tectonic tremor is abundant along the Mexican subduction zone, it offers an opportunity to determine anisotropy parameters in this region. Polarization and splitting analyses were performed using broadband, three-component seismograms. Results show that splitting times range between 0.07 and 0.34 s. These values are similar to the splitting magnitudes typically observed in the continental crust. Fast polarization azimuths are variable in Jalisco, Colima, and Michoacán, but some of them tend to align with the regional stress field (margin-normal, maximum horizontal compressive stresses). On the other hand, the fast axes at the remaining stations are margin-parallel, suggesting that in this case anisotropy could be controlled by active crustal faulting or geological structures. In Oaxaca, fast polarization directions tend to align with Tertiary inactive faults and are oblique with respect to the local stress field, which suggest that anisotropic geological structures are the source of anisotropy.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉An indirect boundary integral equation method (IBIEM) is proposed to investigate the scattering and diffraction of the qP-, qSV- and SH-waves by a three-dimensional (3-D) canyon embedded in a multi-layered transversely isotropic (TI) half-space. First, by employing the Fourier expansion in the circumferential direction and Hankel transform in the radial direction in a cylindrical coordinate system, the 3-D dynamic Green's functions for concentrated loads are derived semi-analytically. The scattered fields, which arise from the presence of the canyon irregularity, are constructed by applying a set of fictitious concentrated loads on the auxiliary surface and calculating the Green's functions. The forces are applied on the auxiliary surface instead of on the real canyon boundary to avoid singularity. Then, the direct stiffness method is adopted to solve the dynamic responses of the displacements and stresses due to the body waves incident from the underlying TI half-space (free fields). Finally, the boundary condition is introduced to determine the densities of the concentrated loads in a least-squares manner. And the total dynamic responses are acquired by adding the free fields to the scattered fields. The numerical stability of the method is examined, and the validity of the proposed formulations is achieved by comparing its results with the existing numerical solutions. By taking a semi-spherical canyon cut in a homogeneous TI half-space, in a single TI layered half-space and in a three layered TI half-space as examples, numerical calculations are performed in the frequency domain. And the effects of material anisotropy, type of incident wave, incident angle, frequency of excitation and sedimentary sequence of the TI layers on the displacement amplitudes are studied. Numerical results illustrate that the dynamic responses for the TI medium can be significantly different from those for the isotropic case, and the amplification spectral is strongly dependent on the type of incident wave, incident angle, frequency of excitation and location of the observation points. The variation of the TI parameters alters the resonant frequencies of the layer and further alters the dynamic interaction between the layer and the canyon. The sedimentary sequence of the layers is crucial for accurate assessment of the surface ground motion around the canyon.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We present a new global model for the Earth’s lithosphere and upper mantle (〈span〉LithoRef18〈/span〉) obtained through a formal joint inversion of 3D gravity anomalies, geoid height, satellite-derived gravity gradients and absolute elevation complemented with seismic, thermal and petrological prior information . The model includes crustal thickness, average crustal density, lithospheric thickness, depth-dependent density of the lithospheric mantle, lithospheric geotherms, and average density of the sublithospheric mantle down to 410 km depth with a surface discretization of 2°×2°. Our results for lithospheric thickness and sublithospheric density structure are in excellent agreement with estimates from recent seismic tomography models. A comparison with higher resolution regional studies in a number of regions around the world indicates that our values of crustal thickness and density are an improvement over a number of previous global crustal models. Given the strong similarity with recent tomography models down to 410 km depth, 〈span〉LithoRef18〈/span〉 can be readily merged with these seismic models to include seismic velocities as part of the reference model. We include several analyses of robustness and reliability of input data, method and results. We also provide easy-to-use codes to interrogate the model and use its predictions for the development of higher-resolution models.Considering the model‘s features and data fitting statistics, 〈span〉LithoRef18〈/span〉 will be useful in a wide range of geophysical and geochemical applications by serving as a reference or initial lithospheric model for i) higher-resolution gravity, seismological and/or integrated geophysical studies of the lithosphere and upper mantle, ii) including far-field effects in gravity-based regional studies, iii) global circulation/convection models that link the lithosphere with the deep Earth, iv) estimating residual, static and dynamic topography, v) thermal modelling of sedimentary basins, v) studying the links between the lithosphere and the deep Earth, among others. Several avenues for improving the reliability of 〈span〉LithoRef18〈/span〉’s predictions are also discussed. Finally, the inversion methodology presented in this work can be applied in other planets for which potential field data sets are either the only or major constraints to their internal structures (e.g. Moon, Venus, etc).〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The accuracy and efficiency of numerical simulations of seismic wave propagation and earthquake ground motion in realistic models strongly depend on discrete grid representation of the material heterogeneity and attenuation. We present a generalization of the orthorhombic representation of the elastic medium (Kristek 〈span〉et al.〈/span〉 2017) to the viscoelastic medium to make it possible to account for a realistic attenuation in a heterogeneous viscoelastic medium with material interfaces. An interface is represented by an averaged orthorhombic medium with rheology of the Generalized Maxwell body (GMB-EK, equivalent to the Generalized Zener body). The representation is important for the possibility of applying one explicit finite-difference scheme to all interior grid points (points not lying on a grid border) no matter what their positions are with respect to the material interface. This is one of the key factors of the computational efficiency of the finite-difference modelling. Smooth or discontinuous heterogeneity of the medium is accounted for only by values of the effective (i.e. representing reasonably averaged medium) grid moduli and densities. Accuracy of modelling thus very much depends on how the medium heterogeneity is represented/averaged. We numerically demonstrate accuracy of the developed orthorhombic representation. The orthorhombic representation neither changes the structure of calculating stress-tensor components nor increases the number of arithmetic operations compared to a smooth weakly heterogeneous viscoelastic medium. It is applicable to the velocity-stress, displacement-stress and displacement FD schemes on staggered, partly-staggered, Lebedev and collocated grids. We also present an optimal procedure for a joint determination of the relaxation frequencies and anelastic coefficients.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉This paper develops an analytical solution for the transient response of a semi-infinite solid with a fluid surface subjected to an arbitrary line source in the solid. The line source is first decomposed into P- and S-pulses, and the wave fields for transient responses of the fluid and semi-infinite solid are presented. Applying the Fourier and Laplace transform methods, the analytical solution in the transform domain, which is an integral solution, is obtained. Using de Hoop's method, a suitable distortion of the contour is provided to change the path of integration and the positions of the singularities and branch points, and a new form of integration is obtained. As the new integration has a form similar to that of the Laplace transform pair, the inverse Laplace transform is directly obtained, and the analytical solution in the time domain is developed. In practice, the compressional wave velocity in the solid is usually larger than that in the fluid. For the cylindrical S-pulse line source, the head wave can always be observed in the fluid-solid system. For the cylindrical P-pulse line source, the head wave can be observed in the fluid-solid system only when the wave is refracted from the solid to the fluid. Numerical examples are provided to discuss the behaviour of the fluid-solid system.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Continuous long-term sound sources are recorded at hydroacoustic station H03S, a three-element hydrophone array south of Robinson Crusoe Island, between 23 April 2014 and 20 August 2017. The origin of the signal between 3-20Hz is investigated by using cross-correlation, array processing using plane wave beamforming, and spectral analysis. One-bit normalization is successfully applied as a cross-correlation preprocessing step in order to suppress undesired earthquake events in the data. Travel times retrieved from averaged cross-correlations do not yield a coherent array direction of arrival. Averaged envelopes of the cross-correlations, however, indicate a coherent signal approaching H03S from a south-southwest direction. Beamforming indicates two dominant back azimuth directions: 172°-224° (Antarctica) and 242° (Monowai Volcanic Seamount). This continuous source field creates possibilities to investigate the applicability of acoustic thermometry at hydrophones H03 S1-S2. Cross-correlation and array processing indicate significant directional variation in local modal propagation, most likely related to the steep slope in the bathymetry near H03S. In addition, it is demonstrated that the ambient noise field is not sufficiently equipartitioned. It is shown that this causes a large error in the estimated temperature, primarily due to the short receiver spacing. These large errors have not been addressed in previous studies on deep-ocean acoustic thermometry. Hence, it is shown that acoustic thermometry does not perform well on small arrays such as H03S. The power spectral density yields a strong broadband signal in January - March, most likely related to iceberg noise. A narrow banded signal around 17Hz between April and September corresponds to whale calls. The best-beam sound pressure levels towards Antarctica are compared to ERA5 climatologies for sea ice cover and normalized stress into the ocean, supporting the hypothesis of iceberg noise.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The Ground-source Airborne Time-domain Electromagnetic (GATEM) method is an efficient electromagnetic prospecting technique for engineering, geological surveying and mineral exploration. The electromagnetic response and inversion are obtained via the electrical conductivity of the geological medium, which is conventionally characterized by a homogeneous half-space or a piecewise smooth spatial distribution. However, these approximations hardly reflect the conditions of real geological media because of their rough, self-similar characteristics. This paper primarily considers the GATEM response and inversion with a long grounded wire source in a rough geological medium based on the fractional diffusion. To manage the rough geological medium, the fractional Maxwell equation is solved in the Laplace domain, and then the solution is transformed to the time domain using Guptasarma¡®s digital filter method. According to the mapping relationship between the resistivity and the induced voltage theoretically calculated in a rough geological medium, the neural network method is used to invert the resistivity. The characteristics of the response in a rough medium are analysed. Then, the inversion results are verified for a layered rough geological medium model; these results are closer to those of the real model than the results of the method derived using a smooth spatial distribution. Finally, the method is applied to GATEM field data in Jincheng, Shanxi Province, China. The inversions performed using a rough geological medium result in a reliable conductivity distribution.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We introduce a novel block rational Krylov method to accelerate three-dimensional time-domain marine controlled-source electromagnetic modelling with multiple sources. This method approximates the time-varying electric solutions explicitly in terms of matrix exponential functions. A main attraction is that no time stepping is required, while most of the computational costs are concentrated in constructing a rational Krylov basis. We optimize the shift parameters defining the rational Krylov space with a positive exponential weight function, thereby producing smaller approximation errors at later times and reducing iteration numbers. Furthermore, for multi-source modelling problems, we adopt block Krylov techniques to incorporate all source vectors in a single approximation space. The method is tested on two examples: a layered seafloor model and a 3D hydrocarbon reservoir model with seafloor bathymetry. The modelling results are found to agree very well with those from 1D semi-analytic solutions and finite-element time-domain solutions using a backward Euler scheme, respectively. Numerical benchmarks show that the block method benefits from better memory and cache efficiency, resulting in about 1.26 to 1.48-fold speedup compared to non-block methods. Further efficiency gains are achieved through optimized rational Krylov techniques, allowing our approach to outperform classical time stepping schemes.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Seismic inversion is essentially ill-posed because of inaccurate and insufficient data. Existing methods for improving stability and reducing the size of solution space are usually model-driven. These methods often rely on sophisticated mathematical models associated with prior knowledge or expectation of the structure, distribution, and correlation of subsurface parameters. In general, model-driven methods are established only for particular usage situations, which are known to have poor adaptivity. Also, because of their abstraction and simplification, these methods often have difficulty achieving satisfactory accuracy and resolution for complex geology. In this paper, we present a data-driven inversion framework using cooperative sparse representation (CSR) for multi-parameter simultaneous inversion. We first let the system learn a joint dictionary from well-log data using multi-component dictionary learning algorithm. Such a learned dictionary captures structural features and correlations of the elastic parameters. Then, in the inversion process, we replace the constraints of mathematical models with that of cooperative sparse representation over the learned dictionary on model parameters. This multi-parameter inversion framework incorporates the correlation among different parameters and is data-driven. It only relies upon the widely accepted assumption that subsurface parameters have a certain degree of consistency and lateral continuity within the same geological formation. Tests on a complex synthetic example and field data demonstrate that our algorithm is effective in improving the resolution and accuracy of solutions, especially for the density parameter that is insensitive to amplitude but extremely important for reservoir identification.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The structure and processes of collisional orogens at the edges of cratons are complex and the extent to which the crust is affected by subduction of cratonic material, slab break-off and post-collisional magmatism is poorly constrained. The Carpathian arc is an Alpine orogenic system, involving the collision and subduction of the East European Craton and the closure of the Tethys ocean. Unusual Neogene back-arc magmatism developed after subduction ended ∼9 Ma ago and upper-mantle Mw〉7 earthquakes are associated with a relic slab actively breaking-off in the Vrancea Zone. To investigate the crustal and uppermost mantle structure, we analyse teleseismic earthquakes recorded at 56 broadband seismic stations in Romania and Moldova. Using phase-weighted 〈span〉H〈/span〉 − κ stacking of receiver functions, we estimate the bulk crustal thickness and 〈span〉Vp〈/span〉/〈span〉Vs〈/span〉 ratio, a good discriminant of crustal composition. This technique was used on all stations and assumes a single-layer crust with a discontinuity at its base. Additionally, by jointly inverting receiver functions and ambient noise we obtain shear wavespeed (〈span〉Vs〈/span〉) models to ∼60 km depth beneath 34 stations and provide another independent measure of the Moho depth and its apparent seismic sharpness. Crustal thickness estimates from the two methods are inconsistent in regions where the Moho is gradational or an intra-crustal discontinuity is more dominant. We support this interpretation with a range of synthetic tests. Above the actively detaching Vrancea slab, multiple intracrustal discontinuities and a gradational Moho are identified on 〈span〉Vs〈/span〉 profiles, and high 〈span〉Vp〈/span〉/〈span〉Vs〈/span〉 (〉1.83) are consistent with magmatic underplating. On the back-arc side of the slab, magmatism-pervaded crust has low 〈span〉Vp〈/span〉/〈span〉Vs〈/span〉 signature (〈1.73), consistent with felsic composition or the presence of water without extensive partial melt. In the back-arc Transylvanian Basin and Apuseni Mountains, the upper mantle has anomalously low 〈span〉Vs〈/span〉 (〈4 km/s to 60 km depth), the crust is dominantly felsic and thin (〈35 km), with a Moho probably inherited from the obducted Neo-Tethys ophiolites. On the fore-arc side, sedimentary basins have low 〈span〉Vs〈/span〉 (〈2.5 km/s) and high 〈span〉Vp〈/span〉/〈span〉Vs〈/span〉 signatures (〉1.8), indicating a high fluid-pressure regime. Finally, Precambrian units in the foreland are seismically heterogenous, with Moho depths of 30 km-50 km, and extend ∼100 km beneath the South and Southeast Carpathians, implying that the Vrancea slab underlies Precambrian-aged crust.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The formation of fold-thrust belts at convergent margins is a dynamic process. Accretion of weak sediments to the front of the overriding plate results in crustal thickening and continued flexural subsidence of the underthrusting plate. Fold-thrust belts are often treated as a Coulomb wedge having self-similar geometries with a critical taper, and either a rigid or isostatically compensated base. In this paper we build upon this work by developing a new dynamic model to investigate both the role of the thickness and material properties of the incoming sediment, and the flexure in the underthrusting plate in controlling the behaviour and evolution of fold-thrust belts. Our analysis shows that the evolution of fold-thrust belts can be dominated by either gravitational spreading or vertical thickening, depending on the relative importance of sediment flux, material properties and flexure. We apply our model to the Makran accretionary prism and the Indo-Burman Ranges, and show that for the Makran flexure must be considered in order to explain the dip of the sediment-basement interface from seismic reflection profiles. In the Indo-Burman Ranges, we show that incoming sediment thickness has a first-order control on the variations in the characteristics of the topography from north to south of the Shillong Plateau.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The converted wave data (P-to-s/S-to-p), commonly termed as receiver functions, contain noises of various origins. Such noises may influence the modeling and may sometimes lead to over interpretations of the data. In order to suppress noise, we use a robust sparsity enhancing tool, i.e. the Seislet Transform (ST), to process receiver function data by applying regularization in the seislet domain. The transform utilizes the multiscale orthogonal basis and the basis-functions are oriented along the dominant seismic phases following local linearity. The inversion results of both the synthetic and field examples from the Hi-CLIMB network and station HYB from the Indian shield show an excellent performance over the original data sets.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉In the absence of quantitative relationships, a cross-gradient function can be used to correlate unknown physical properties in a joint inversion of geophysical data sets. It introduces a structural correlation between properties. A commonly used, inexact approach adds a weighted cross-gradient term as a penalty to the cost function being minimized during inversion. This weighting factor needs to be tuned to balance the regularization and cross gradient terms. In this paper we propose non-linear weighting for the cross-gradient function which addresses the very different magnitudes of the cross gradient and regularization terms. This approach also couples the weighting factors for the regularization and correlation terms reducing the number of tuning parameters. The approach is investigated for a synthetic case. Results are also shown for the 3D joint inversion of high-resolution magnetic and gravity anomaly data from Southern Queensland in Australia with over 30 Million cells.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The goal of next generation gravity missions (NGGM) is to improve the monitoring of mass transport in the Earth system by an increased space–time sampling capability as well as higher accuracies of a new generation of instrumentation. They should be able to fulfil the scientific and societal needs of providing high-resolution short-time gravity field solutions for geophysical applications like for for example service applications such as flood and drought monitoring and forecast or applications in water management. To facilitate this need a near-real time (NRT) processing scheme based on a coparametrization of low-resolution daily and longer-term gravity field solution, combined with a sliding window averaging, was set up. In contrast to other strategies that are usually based on Kalman filtering, the proposed NRT concept is independent of any prior information about the temporal gravity field, and does not require any regularization. The enhanced spatial-temporal resolution opens the possibility to self-dealias high-frequency atmospheric and oceanic signals, and additionally provides gravity field solutions on short timescales. In order to quantify the capabilities of the proposed NRT approach, a numerical closed-loop simulation of a low-low satellite-to-satellite tracking (ll-sst) mission for a two-pair Bender-type constellation with realistic noise assumptions was performed. While for the daily parametrization a spherical harmonics degree and order of 15 turns out to be a favourable choice, by applying the sliding window NRT approach stable daily gravity field estimates up to degree/order 50 with latencies of down to 1 d could be achieved.〈/span〉

Description: 〈span〉The title of my recent article (Flament 〈a href="http://academic.oup.com/gji#bib1"〉2019〈/a〉) should be ‘Present-day dynamic topography and lower-mantle structure from palaeogeographically constrained mantle flow models’, not ‘Present-dayd dynamic topography and lower-mantle structure from palaeogeographically constrained mantle flow models’. This spelling mistake was introduced at proof stage and was present in the published version of the article.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Separating compressional and shear wavefields is an important step in elastic reverse-time migration, which can remove wave-mode crosstalk artefacts and improve imaging quality. In vertical (VTI) and titled (TTI) transversely isotropic media, the state-of-the-art techniques for wavefield separation are based on either non-stationary filter or low-rank approximation. They both require intensive Fourier transforms for models with strong heterogeneity. Based on the eigenform analysis, we develop an efficient pseudo-Helmholtz decomposition method for the VTI and TTI media, which produces vector 〈span〉P〈/span〉 and 〈span〉S〈/span〉 wavefields with the same amplitudes, phases and physical units as the input elastic wavefields. Starting from the elastic VTI wave equations, we first derive the analytical eigenvalues and eigenvectors, then use the Taylor expansion to approximate the square-root term in the eigenvalues, and finally obtain a zero-order and a first-order pseudo-Helmholtz decomposition operator. Because the zero-order operator is the true solution for the case of ϵ = δ, it produces accurate wavefield separation results for elliptical anisotropic media. The first-order separation operator is more accurate for non-elliptical anisotropy. Since the proposed pseudo-Helmholtz decomposition requires solving an anisotropic Poisson’s equation, we propose two fast numerical solvers. One is based on the sparse lower-upper (LU) factorization, which can be repeatedly applied to the input elastic wavefields once computing the lower and upper triangular matrices. The second solver assumes the model parameters are laterally homogeneous within a given migration aperture. This assumption allows us to efficiently solve the anisotropic Poisson’s equation in the 〈span〉z〈/span〉 − 〈span〉kx〈/span〉 domain, where 〈span〉kx〈/span〉 and 〈span〉z〈/span〉 denote the horizontal wavenumber and depth, respectively. Using the coordinate transform, we extend the pseudo-Helmholtz decomposition to the TTI media. The separated vector wavefields are used to produce 〈span〉PP〈/span〉 and 〈span〉PS〈/span〉 images by applying a dot-product imaging condition. Several numerical examples demonstrate the feasibility and applicability of the proposed methods.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉When performing a cooperative inversion using a structural constraint, extracting an accurate average direction from the high-resolution model is important because two models with different resolutions should be used in the procedure. However, when the average direction of the high-resolution model, which indicates the major change direction of the model at each inversion block location, is calculated, the conventional gradient methods such as cross-gradient have a limitation on components having opposite directions. Therefore, in this study, an effective average-direction extraction algorithm was developed by introducing the concept of the structure tensor to accurately calculate average-direction information. And finally, based on the extracted average-direction information, structure-tensor-constrained cooperative inversion algorithm was proposed. To verify the effectiveness of the proposed inversion method, the method was tested using data obtained from the synthetic model containing an anomalous body with a complicated shape and the result compared with results of individual EM inversion and the conventional cross-gradient-constrained cooperative inversion. Lastly, to evaluate the performance with a realistic mode, the proposed cooperative inversion was applied to the data acquired using the complex SEG Advanced Modelling Program model. In all experiments, the cooperative inversion with the structure-tensor constraint provided better location estimation results as well as better estimations of the shape of the anomaly. In addition, the resistivity distribution of the anomaly was estimated to be closer to the truth in the inversion result.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉A complex Meso-Cenozoic history characterizes the offshore Indus Basin Pakistan, including rifting, drifting, collision, transpression and recent inverted extension that is poorly understood. We use reflection seismic data to document the geometry of the structures that accommodated the deformation and offshore well data to correlate the sedimentary sequences which were affected by the different tectonic episodes. Synthetic seismograms have been prepared for calibrating seismic data with well information for identification of key reflections representative of major lithological horizons. Seismic stratigraphy helped in dividing the sedimentary successions into rifted margin and transpressive margin sequences. Additionally, the basement represented by chaotic and semitransparent reflections, and Moho discontinuity by moderate amplitude bright reflections are also identified on seismic records. Syn-depositional normal faulting in extensional regime (Cretaceous–Oligocene) are identified on seismic profiles that converted into strike slip faults in Miocene and listric faults in Recent sediments forming inverted extensional regime. We infer that the basin inversion from extension to inverted extension and transpression occurred in recent times. Number of depositional features including shelves, deltaic deposits, canyons and channel levee complexes have been identified that developed in different tectonic episodes. Five tectonic subsidence curves using backstripping procedure have been calculated in order to constrain the geodynamic evolution of the Meso-Cenozoic sedimentary succession of the offshore Indus Pakistan. The analysis of tectonic subsidence curves, covering a time span of 140 Ma, allowed us to quantify the total subsidence and to document the activities of syn-rift and post-rift faulting systems. Different stages, in terms of duration and magnitude of subsidence-uplift trends, have been identified in the evolution of the basin. Two uplift events are prominent in Eocene (∼40 Ma) and Miocene time (∼12 Ma). A rapid resumption of subsidence with non-conformity in Early Palaeocene is seen along the margin.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We investigate the physical processes that generate seismicity during hydraulic fracturing. Fluid processes (increases in pore pressure and poroelastic stress) are often considered to be the primary drivers. However, some recent studies have suggested that elastic stress interactions may significantly contribute to further seismicity. In this work we use a microseismic data set acquired during hydraulic fracturing to calculate elastic stress transfer during a period of fault activation and induced seismicity. We find that elastic stress changes may have weakly promoted initial failure, but at later times stress changes generally acted to inhibit further slip. Sources from within tight clusters are found to be the most significant contributor to the cumulative elastic stress changes. Given the estimated 〈span〉in situ〈/span〉 stress field, relatively large increases in pore pressure are required to reach the failure envelope for these faults—on the order of 10 MPa. This threshold is far greater than the reliable cumulative elastic stress changes found in this study, with the vast majority of events receiving no more than 0.1 MPa of positive Δ〈span〉CFS〈/span〉, further indicating that elastic stress changes were not a significant driver, and that interaction with the pressurized fluid was required to initiate failure. Thus, cumulative stress transfer from small events near the injection well does not appear to play a significant role in the reactivation of nearby faults.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The influence of cooling rate on the intensity of thermoremanent magnetization (TRM) and the necessity to correct archaeo/palaeointensities for this effect have long been recognized. However the reliability of the correction is still questioned. We studied 35 bricks baked in two modern kilns (SK and BK) in known experimental conditions and with measurements of the direction and intensity of the geomagnetic field at the site. The smallest kiln (SK, 0.2 m〈sup〉3〈/sup〉) cooled in around 12 hr and the biggest kiln (BK, 8 m〈sup〉3〈/sup〉) in around 40 hr. Thermomagnetic, hysteresis and backfield curves indicated that the main magnetic carriers were Ti-poor titanomagnetites and Ti-poor titanohematites. The fraction of the TRM carried by Ti-poor titanohematites is the main difference between the two sets of bricks. This fraction is around 5–10 per cent in bricks from BK kiln and up to 40 per cent in those from SK kiln. Intensities of the Earth's magnetic field were determined using the original Thellier-Thellier protocol with correction of TRM anisotropy. The average intensities overestimate the expected field intensity by 5 per cent (SK) and 6 per cent (BK). This result emphasizes the necessity of the cooling rate correction. In order to have a detailed evaluation of the cooling rate effect, we used several slow cooling rates: 0.8, 0.4, 0.2 and 0.1° C min〈sup〉–1〈/sup〉. The correction factors obtained with the 0.8° C min〈sup〉–1〈/sup〉 cooling ranged between –2 and 21 per cent and were proportional to the TRM fraction carried by Ti-poor titanohematite. The higher proportion of these grains in bricks from SK kiln led to an overestimate of the correction factor and an underestimate of the intensity by 7 per cent. However, the expected intensity is recovered when temperature steps higher than 580 °C (i.e. in the range of Ti-poor titanohematites unblocking temperatures) were excluded from the calculation of archaeointensity and cooling rate correction. In the case of the BK kiln bricks, for which Ti-poor titanohematites does not contribute significantly to the TRM, all tested cooling rates give average intensities close to the expected value. Incorrectly estimating the duration of the archaeological cooling has therefore a low impact on the accuracy of the archaeointensity data on these kinds of material.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Time-series of 〈span〉VP〈/span〉/〈span〉VS〈/span〉 ratio have been used to track local changes in elastic properties of rock volumes. Identifying such variations can provide information on the geophysical processes taking place inside a rock volume during the seismic cycle. A value of 〈span〉VP〈/span〉/〈span〉VS〈/span〉 ratio can be computed from traveltime of 〈span〉P〈/span〉 and 〈span〉S〈/span〉 waves generated from a single local event and it is representative of the value of the 〈span〉VP〈/span〉/〈span〉VS〈/span〉 ratio for the rocks traversed by the seismic ray, between the source and the receiver. It is straightforward, during a seismic sequence, to generate time-series of 〈span〉VP〈/span〉/〈span〉VS〈/span〉 ratio for events located close together and a single station. Such time-series should be able to monitor temporal variations of elastic parameters in the rock volume. Due to the very small nature of the expected changes in 〈span〉P〈/span〉- and 〈span〉S〈/span〉-wave velocity, the evaluation of 〈span〉VP〈/span〉/〈span〉VS〈/span〉 ratio time-series has been problematic in the past, and subjective choices about, for example the time-averaging scheme applied or event selection for constructing the time-series, have been proven to strongly affect the outcomes of the analysis. In this contribution, we present the application of a new methodology for a statistical evaluation of changes in 〈span〉VP〈/span〉/〈span〉VS〈/span〉 ratio time-series. The new methodology belongs to the wide class of ‘change-point analysis’ algorithms and is developed in the framework of Bayesian inference. The posterior probability distribution (PPD) of the change-point locations is obtained using a trans-dimensional Markov chain Monte Carlo (trans-D McMC) algorithm, where the existence and number of change-points is directly dictated by the data themselves. We apply the new algorithm to the seismic catalogue produced by the Alto Tiberina Near Fault Observatory seismic network (Northern Apennines, Italy). Here the high rate of background seismic release and the dense seismic network allow for a robust statistical analysis. The occurrence of change-points in 〈span〉VP〈/span〉/〈span〉VS〈/span〉 time-series identified with the proposed procedure is represented in space and time. The space–time distributions of change-points in the study area shows a clear peak of change-points following the occurrence of local main events, clustered along the main fault system activated. The robustness of the proposed approach makes it appropriate as an automatic, real-time tool for monitoring rock property changes related to seismic activity.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Seven years of 〈span〉GRACE〈/span〉 intersatellite range-rate measurements are used to test the new ocean tide model FES2014 and to compare against similar results obtained with earlier models. These qualitative assessments show that FES2014 represents a marked improvement in accuracy over its earlier incarnation, FES2012, with especially notable improvements in the Arctic Ocean for constituents K〈sub〉1〈/sub〉 and S〈sub〉2〈/sub〉. Degradation appears to have occurred in two anomalous regions: the Ross Sea for the O〈sub〉1〈/sub〉 constituent and the Weddell Sea for M〈sub〉2〈/sub〉.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Cratons are the oldest parts of the lithosphere, some of them surviving since Archean. Their long-term survival has sometimes been attributed to high viscosity and low density. In our study, we use a numerical model to examine how shear tractions exerted by mantle convection work to deform cratons by convective shearing. We find that although tractions at the base of the lithosphere increase with increasing lithosphere thickness, the associated strain-rates decrease. This inverse relationship between stress and strain-rate results from lateral viscosity variations along with the model’s free-slip condition imposed at the Earth’s surface, which enables strain to accumulate along weak zones at plate boundaries. Additionally, we show that resistance to lithosphere deformation by means of convective shearing, which we express as an apparent viscosity, scales with the square of lithosphere thickness. This suggests that the enhanced thickness of the cratons protects them from convective shear and allows them to survive as the least deformed areas of the lithosphere. Indeed, we show that the combination of a smaller asthenospheric viscosity drop and a larger cratonic viscosity, together with the excess thickness of cratons compared to the surrounding lithosphere, can explain their survival since Archean time.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉A promising way to perform seismological studies in the Arctic region is deploying seismic stations on ice floes. The pioneering works by the Alfred Wegener Institute (AWI) Bremerhaven (Schlindwein 〈span〉et al.〈/span〉 2007; Laderach and Schlindwein, 2011) have demonstrated the efficiency of such floating networks to explore local and regional seismicity and to build 3D seismic models. However, problems remain, related to the identification of different types of seismic waves, particularly S-waves. Here, we perform 2D and 3D numerical simulations of seismic waves emitted by an earthquake to explore the possibility of recording different phases on the sea surface. We use different types of simple shear source models, namely, strike-slip, vertical displacement and normal faults. In the calculated wave field, we obtain three major types of seismic waves recorded on the sea surface: 〈span〉Pw, Sw〈/span〉 and 〈span〉SPw〈/span〉 (〈span〉w〈/span〉 denotes an acoustic wave in the water layer) and numerous multiple waves. The clarity of the recorded phases strongly depends on the type of wave, source mechanism, epicentral distance, thickness of the water layer and depth of the source. For example, the 〈span〉Pw〈/span〉 phase is clearest for the strike-slip mechanism, less clear for the normal fault and almost invisible for the vertical displacement. The 〈span〉Sw〈/span〉 phase is observable in all of these cases; however, it can be confused with the 〈span〉SPw〈/span〉 phase that arrives earlier. In addition, at some distances, the 〈span〉Sw〈/span〉 wave interferes with the multiple 〈span〉Pw2〈/span〉 wave and therefore is hardly detectable. In summary, the numerical simulations in a model with a water layer have demonstrated several non-obvious features of wave propagation that should be taken into account when analysing experimental data recorded on ice floes.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉From a suite of 56 chemically-driven dynamo simulations with aspect ratio χ (inner to outer core radii) ranging from 0.10 to 0.44, we conduct the first systematic investigation of the impact of inner-core size on the reversing behaviour of dynamos. We show that the growth of the inner core leads to a transition between a “small inner-core” regime (χ ≤ 0.18), when the field produced is intermediately strong and dipolar, and a “large inner-core” regime (χ 〉 0.26), when the field is stronger and more dipolar. During that transition the field is weaker and slightly less dipolar. For aspect ratios 0.20 ≤ χ ≤ 0.22, reversal frequencies may be more sensitive to changes in the vigour of the convection, allowing high frequencies to be reached much more easily. Although other factors than the size of the inner core likely contribute to controlling the reversal frequency of the Earth’s dynamo, we hypothesise that the occurrence of such a transition for the Earth’s core between the end of the Precambrian and the end of the Devonian could possibly account for the manifestation of an unusual long-lasting episode of predominantly reversal hyperactivity and complex low intensity fields during that still poorly documented period of time.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Velocity macro-model building is an essential step of the seismic imaging workflow. Indeed, obtaining acceptable results through migration or full waveform inversion is highly dependent on the kinematic accuracy of the background/initial velocity model. Two decades ago, stereotomography was proposed as an alternative to reflection traveltime tomography, the first relying on semi-automatic picking of locally coherent events associated with small reflection or diffraction segments tied to scatterers in depth by a pair of rays, while the latter on interpretive picking of laterally-continuous reflections. The flexibility of stereotomography paved the way for many developments that have shown the efficiency of the method whilst emphasizing on the complementary information carried out by traveltimes and slopes of locally-coherent events. A recent formulation recast stereotomography under a matrix-free formulation based on eikonal solvers and the adjoint-state method. In the latter, like in the previous works, the scatterer positions and the velocity field are updated jointly to tackle the ill-famed velocity-position coupling in reflection tomography. Following on from this adjoint-state formulation, we propose a new parsimonious formulation of slope tomography that offers the chance to restrain the problem to minimizing the residuals of a single data class being a slope, in search of a sole parameter class being the subsurface velocity field. This parsimonious formulation results from a variable projection, which is implemented by enforcing a consistency between the scatterer coordinates and the velocity macro-model through migration of kinematic attributes. We explain why the resulting reduced-parametrization inversion is more suitable for tomographic problems than the most common joint inversion strategy. We benchmark our method against the complex Marmousi model along with a validation through time domain full waveform inversion and then present the results of a field data case study.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Stoneley waves propagating in a borehole will produce reflected waves when they encounter a porous formation. In previous studies, the simplified Biot-Rosenbaum theory is used to calculate the Stoneley wave reflection coefficient for porous formations. Such simplified theory ignores the effect of formation frame elasticity, thus can't obtain accurately the Stoneley wave reflection for the porous formation with small stiffness. In this study, to take the effect of formation frame elasticity into account, we use the Biot's theory in the low-frequency limit to simulate the Stoneley wave reflection by employing the velocity-stress finite-difference time-domain (FDTD) method. In addition, as permeability of the formation varies in the axial direction, and the viscosity of the pore fluid also changes in a radial direction due to mud invasion, it is difficult to use the simplified theory to calculate the reflection coefficient of the Stoneley wave in such heterogeneous cases. Therefore, this study first investigates the effect of the formation permeability heterogeneity on Stoneley wave reflection coefficient by this FDTD method.The FDTD method is verified by a comparison with the real axis integration (RAI) method with respect to the Stoneley wave propagation in a borehole surrounded by a homogeneous porous formation. The reflection coefficient obtained by the FDTD method is smaller than that using the simplified theory, which shows that elasticity of the formation frame affects Stoneley wave reflection: the effect of elasticity on the reflection coefficient is greater when the formation frame is less rigid or when the porosity and permeability of the formation are lower. According to the simulation results of the FDTD method, a modified simplified theory which can improve the calculation accuracy of Stoneley wave reflection coefficient is proposed. Furthermore, the effects of the permeability heterogeneity on the Stoneley wave reflection are investigated: the reflection coefficient peaks change when permeability alters in an axial direction, and the peak interval increases. For the mud invasion model, the reflection coefficient is almost identical to that of the homogeneous model, which has the same permeability as the borehole wall of the mud invasion model.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We propose a new Bayesian method to reveal the 〈span〉Vs〈/span〉 structure of the near surface of the earth using spatial autocorrelation (SPAC) functions and apply this new method to synthetic, broadband, and geophone datasets. The principle of SPAC is introduced, and an implementation of the Bayesian Monte Carlo inversion (BMCI) for modeling SPAC coherency functions is described. To demonstrate its effectiveness, BMCI is applied to synthetic tests, data from 14 SPAC array sites in the Salt Lake Valley (SLV), Utah, and two arrays (one broadband and one geophone) located in south central Utah. The 〈span〉Vs〈/span〉 models derived from previous SPAC analysis of the 14 SLV sites differ by 10 per cent at most from those determined by BMCI and lie within uncertainties determined for the BMCI models. These agreements demonstrate the effectiveness of the BMCI method. The synthetic tests and applications to the SLV SPAC data show BMCI has great potential to resolve 〈span〉Vs〈/span〉 structure down to at least 400 m. To achieve resolution for deeper 〈span〉Vs〈/span〉 structure, longer duration deployments, wider array apertures, and additional seismometers or geophones can be employed. Additionally, when the target frequencies are greater than 0.1 Hz, there is no apparent disadvantage in using geophone data for BMCI compared to broadband data. Most significantly, BMCI places a quantifiable constraint on the uncertainties of the 〈span〉Vs〈/span〉 models as well as 〈span〉Vs30〈/span〉.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Although volume 3-D modelling solutions has become widespread in recent time, thin sheet approximation of Earth's conductivity distribution can still serve as a useful tool when quasi-3-D conductivity structures in the heterogeneous subsurface are investigated and the available database of observations is limited to long-period electromagnetic induction data from large-scale deep sounding arrays. We present results of stochastic Monte Carlo Markov Chains (MCMC) inversion of long-period induction arrows based on the Bayesian statistics strategy.We concentrated on the different methodological aspects of MCMC for Gibbs sampling and for adaptive Metropolis algorithm together with convergence of these methods. The results are presented on a case study from the transition zone between the Bohemian Massif and the West Carpathians where a phantom effect caused by superposition of the prominent SW–NE trending Carpathian Conductivity Anomaly and NW–SE trending anomalous structure related to the fault system at the eastern boundary of the Bohemian Massif appears.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉In Magnetotellurics (MT) natural electromagnetic field variations are recorded to study the electrical conductivity structure of the subsurface. Thereby long time-series of electromagnetic data are subdivided into smaller segments, which are Fourier transformed and typically averaged in a statistically robust manner to obtain MT transfer functions. Unfortunately, nowadays the presence of man-made electromagnetic noise sources often deteriorates a significant fraction of the recorded time-series by overprinting the desired natural field variations. Available approaches to obtain undisturbed and high quality MT results include, for example robust statistics, remote reference or multi-station analyses which aim at the removal of outliers or uncorrelated noise. However, we have observed that intermittent noise often affects a certain time span resulting in a second cluster of transfer functions in addition to the expected true MT distribution. In this paper, we present a novel criterion for the detection and pre-selection of EM noise in form of outliers or additional clusters based on a distance measure of each data segment with regard to the centre of the data distribution. For this purpose, we utilize the Mahalanobis distance (MD) which computes the distance between two multivariate points considering the covariance matrix of the data that quantifies the shape and the size of multivariate data distributions. As the MD considers the covariance matrix, it corrects not only for different variances but also for any correlation between the data. The computation of both, the mean value and covariance matrix, is susceptible to ouliers (e.g. noise) and requires a statistically robust estimation. We tested several robust estimators, for example median absolute deviation or minimum covariance determinant algorithm and finally implemented an automatic criterion using a deterministic minimum covariance determinant algorithm. We will present results using MT data from various field experiments all over the world, which illustrate successfull data improvement. This approach is able to remove scattered data points as well as to reject complete data cluster originating from noise sources. However, like all purely statistical algorithms the criterion is limited to cases where the majority of the recorded data is well-behaved, that is noise content is below 50 per cent. If the majority of data points originates from noise sources, the new criterion will fail if used in an automatic way. In these cases, additional input by the user either manually or in an automated fashion can be utilized. We therefore suggest to use an add-on criterion to back the MD selection and subsequent robust stacking in form of a physically motivated constraint based on the magnetic incidence direction. This property indicates whether the magnetic field originates from various sources in the far field or from a strong and well defined source in the near field.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Numerical simulations of earthquake ground motions are used both to anticipate the effects of hypothetical earthquakes by forward simulation and to infer the behavior of the real earthquake source ruptures by inversion of recorded ground motions. In either application it is necessary to assume some Earth structure that is necessarily inaccurate and to use a computational method that is also inaccurate for simulating the wave field Green's functions. We refer to these two sources of error as ‘propagation inaccuracies,’ which might be considered to be epistemic. We show that the variance of the Fourier spectrum of the synthetic earthquake seismograms caused by propagation inaccuracies is related to the spatial covariance on the rupture surface of errors in the computed Green's functions, which we estimate for the case of the 2009 L'Aquila, Italy, earthquake by comparing erroneous computed Green's functions with observed L'Aquila aftershock seismograms (empirical Green's functions). We further show that the variance of the synthetic seismograms caused by rupture variability (aleatory uncertainty) is related to the spatial covariance on the rupture surface of aleatory variations in the rupture model, and we investigate the effect of correlated variations in Green's function errors and variations in rupture models. Thus, we completely characterize the variability of synthetic earthquake seismograms induced by errors in propagation and variability in rupture behavior. We calculate the spectra of the variance of the ground motions of the L'Aquila main shock caused by propagation inaccuracies for two specific broadband stations, the AQU and the FIAM stations. These variances are distressingly large, being comparable or in some cases exceeding the data amplitudes, suggesting that the best-fitting L'Aquila rupture model significantly over-fits the data and might be seriously in error. If these computed variances are typical, the accuracy of many other rupture models for past earthquakes may need to be reconsidered. The results of this work might be useful in seismic hazard estimation because the variability of the computed ground motion, caused both by propagation inaccuracies and variations in the rupture model, can be computed directly, not requiring laborious consideration of multiple Earth structures.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We present a parallel and high-order Nédélec finite element solution for the marine controlled-source electromagnetic (CSEM) forward problem in 3D media with isotropic conductivity. Our parallel Python code is implemented on unstructured tetrahedral meshes, which supports multiple-scale structures and bathymetry for general marine 3D CSEM modeling applications. Based on a primary/secondary field approach, we solve the diffusive form of Maxwell’s equations in the low-frequency domain. We investigate the accuracy and performance advantages of our new high-order algorithm against a low-order implementation proposed in our previous work. The numerical precision of our high-order method has been successfully verified by comparisons against previously published results that are relevant in terms of scale and geological properties. A convergence study confirms that high-order polynomials offer a better trade-off between accuracy and computation time. However, the optimum choice of the polynomial order depends on both the input model and the required accuracy as revealed by our tests. Also, we extend our adaptive-meshing strategy to high-order tetrahedral elements. By using adapted-meshes to both physical parameters and high-order schemes, we are able to achieve a significant reduction in computational cost without sacrificing accuracy in the modeling. Furthermore, we demonstrate the excellent performance and quasi-linear scaling of our implementation in a state-of-the-art HPC architecture.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Modelling marine controlled-source electromagnetic (CSEM) responses in the fictitious time domain is a novel approach, which facilitates the full exploration of EM diffusive properties in the fictitious wave domain (FWD). Concepts, such as reflections, refractions, diffractions and transmissions, which are used for the analysis of elastic wave propagation can thus be adopted in FWD for interpreting CSEM data. In this paper, we use a high-order finite difference time domain (FDTD) algorithm for modelling marine CSEM responses in both the fictitious time domain and the diffusive frequency domain. A complex frequency shifted perfectly matched layer (CFS-PML) boundary condition is adopted to the FDTD modelling. We demonstrate the performance of the CFS-PML boundary condition and validate the high-order FDTD code in the fictitious wave domain with the half-space seawater model and in the frequency domain with the 1D canonical reservoir model. We investigate and analyze the propagation characteristics of electromagnetic fields in the fictitious wave domain, where we apply wave propagation concepts to interpret marine CSEM data. Similarities between wave and field propagations relevant for marine CSEM data are demonstrated through several 1D to 3D numerical examples.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The Amatrice-Norcia-Visso sequence is characterized by complex behavior that is somewhat atypical of mainshock-aftershock sequences (Marzocchi et al. 2017), as there were multiple large mainshocks that continued for months.In this study we focus on the Amatrice sequence (mainshock August 24〈sup〉th〈/sup〉, 2016 Mw = 5.97) to evaluate the apparent stress values and magnitude dependent scaling in order to improve our knowledge of processes that control small and large earthquakes within this active region of Italy. Apparent stress is proportional to the ratio of radiated seismic energy and seismic moment, and as such, these stress parameters play an important role in hazard prediction as they have a strong effect on the observed and predicted ground shaking.We analyze 83 events of the sequence from August 24〈sup〉th〈/sup〉 to October 16〈sup〉th〈/sup〉 2016, within a radius of 20 km from the mainshock and with an Mw ranging between 5.97 and 2.72.Taking advantage of the averaging nature of coda waves, we analyze coda-envelope-based spectral ratios (Mayeda et al. 2007) between neighboring event pairs. We use equations proposed by Walter et al. 2017 to consider stable, low-frequency and high frequency spectral ratio levels (LFL and HFL, respectively) which provide measures of the corner frequency and apparent stress ratios of the events within the sequence.The results demonstrate non-self-similar behavior within the sequence, suggesting a change in dynamics between the largest events and the smaller aftershocks.The apparent stress and corner frequency estimates are compared to those obtained by Malagnini and Munafò (2018) who utilized hundreds of direct S-wave spectral ratio measurements to obtain their results. Although our analysis is based only on 83 events, our results are in very good agreement, demonstrating once more that the use of coda waves is very stable (Sato and Fehler, 1998; Mayeda et al. 2003, 2007) and provides lower variance measures than those using direct waves.A comparison with recent Central Apennines source scaling models derived from various seismic sequences (1997-1998 Colfiorito, 2002 San Giuliano di Puglia, 2009 L'Aquila), shows that the Amatrice sequence source scaling in this study is well represented by the models proposed by Pacor et al. 2016 and Malagnini and Mayeda 2008.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉With abundant seismic data for small repeating earthquakes, it is important to construct a dynamic model that can explain various aspects of related observations. In this work, we study small repeating earthquakes on a fault governed by rate- and state-dependent friction laws. The earthquakes occur on a velocity-weakening patch surrounded by a much larger velocity-strengthening region. The whole fault is subject to long-term tectonic loading. The model with a circular patch and the aging form of rate- and state-dependent friction laws has been shown to reproduce the scaling of recurrence time versus seismic moment for small repeating earthquakes in a previous study. Here we investigate the behaviour of small repeating earthquakes in related models under different scenarios, including several forms of the state evolution equations in rate- and state-dependent friction laws, rectangular velocity-weakening patch geometries, quasi-dynamic versus fully dynamic representation of inertial effects and 2-D versus 3-D simulations. We find that the simulated scalings between the recurrence time and seismic moment for these different scenarios is similar while differences do exist. We propose a theoretical consideration for the scaling between the recurrence time and seismic moment of small repeating earthquakes. For patch radii smaller than or comparable to the full nucleation size, the scaling is explained by the increase of seismic to aseismic slip ratio with magnitude. For patch radii larger than the full nucleation size, the scaling is explained by the model in which the recurrence time is determined by the earthquake nucleation time, which is in turn determined by the time for aseismic slip to penetrate the distance of the full nucleation size into the patch. The obtained theoretical insight is used to find the combinations of fault properties that allow the model to fit the observed scaling and range of the seismic moment and recurrence time.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We use a version of source encoding to facilitate the calculation of the gradient of a misfit function independent of the number of sources or receivers. The crosstalk-free method requires only two numerical simulations per iteration, namely, one source-encoded forward simulation and one source-encoded adjoint simulation. Importantly for practical applications, each source does not need to be recorded by all receivers. The method is implemented in two complementary ways. In the first approach the encoded forward and adjoint wavefields are run until they reach steady state, at which point they are ‘decoded’ to obtain their stationary parts. These stationary parts are combined to calculate the misfit gradient by summing their respective contributions. In the second approach the steady-state encoded forward and adjoint wavefields are convolved over a time period proportional to the inverse of the encoded frequency spacing. Using this strategy, the encoded forward and adjoint wavefields do not need to be decoded, nor is there a need to calculate or store intermediary stationary contributions to the gradient. We consider a wide variety of source-encoded misfit functions, including waveform differences, phase and amplitude measurements,‘double-difference’ phase and amplitude measurements, cross-correlation traveltime measurements, and a generic adjoint tomography misfit function. When measurements in specific time windows are involved in the construction of the source-encoded misfit function, as in adjoint tomography, the computational cost scales linearly with the number of seismic sources, because the necessary synthetic seismograms must be computed individually. In contrast, when using ‘super measurements’ based on source-encoded Fourier coefficients of entire observed and simulated seismograms, as in pure full waveform inversion, one iteration requires just two numerical simulations, independent of the number of sources and receivers. We illustrate the method based on examples from both earthquake and exploration seismology, highlighting inversion options and strategies involving frequency- and time-domain encoding, decoding (with more than 16 000 frequencies), encoded frequency randomization, encoding multiple frequencies per source, effects of noise, a variable number of receivers per event, various measurements and related misfit functions and attenuation.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Laboratory and field studies have demonstrated the applicability of nanoparticles (NP) for accelerated contaminant degradation. Beside other limitations (e.g. costs, delivery, longevity, non-target specific reactions), concerns of regulators arose regarding toxicity of injected NP and particles delivered off-target (i.e. renegade particles). Renegade particles also significantly reduce the efficiency of the remediation. The delivery of particles off-target is caused, mainly, by unintended fracking, where the fractures act then as preferential flow paths changing the trajectory of the particles. Hence, the real-time monitoring of particle injection is of major importance to verify correct particle delivery and thus help to optimize the remediation strategy. However, to date NP monitoring techniques rely on the analysis of soil and water samples, which cannot provide information about clogging or the formation of fractures away of the sampling points. To overcome these limitations, in this study we investigate the applicability of complex-conductivity imaging (CCI), a geophysical electrical method, to characterize possible pore clogging and fracking during NP injections. We hypothesize that both processes are related to different electrical footprints, considering the loss of porosity during clogging and the accumulation of NP in areas away of the target after fracking. Here, we present CCI results for data collected before and during the injection of Nano-Goethite particles (NGP) applied to enhance biodegradation of a BTEX (benzene, toluene, ethylbenzene and xylene) contaminant plume. Imaging results for background data revealed consistency with the known lithology, while overall high electrical conductivity values and a negligible induced-polarization magnitude correspond with the expected response of a mature hydrocarbon plume. Monitoring images revealed a general increase (∼15 per cent) in the electrical conductivity due to the injected NGP suspension in agreement with geochemical data. Furthermore, abrupt changes in this trend, shortly before daylighting events, show the sensitivity of the method to pore clogging. Such interpretation is in line with the larger variations in CCI resolved in the unsaturated zone, clearly indicating the accumulation of renegade NGP close to the surface due to fracking. Our results demonstrate the applicability of the CCI method for the assessment of pore clogging accompanying particles injection.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The work in this paper is motivated by the increasing use of electrical and electromagnetic methods in geoscience problems where steel-cased wells are present. Applications of interest include monitoring carbon capture and storage and hydraulic fracturing operations. Also of interest is detecting flaws or breaks in degrading steel-casings – such wells pose serious environmental hazards. The general principles of electrical methods with steel-cased wells are understood and several authors have demonstrated that the presence of steel-cased wells can be beneficial for detecting signal due to targets at depth. However, the success of a DC resistivity survey lies in the details. Secondary signals might only be a few percent of the primary signal. In designing a survey, the geometry of the source and receivers, and whether the source is at the top of the casing, inside of it, or beneath the casing will impact measured responses. Also the physical properties and geometry of the background geology, target, and casing will have a large impact on the measured data. Because of the small values of the diagnostic signals, it is important to understand the detailed physics of the problem and also to be able to carry out accurate simulations. This latter task is computationally challenging because of the extreme geometry of the wells, which extend kilometers in depth but have millimeter variations in the radial direction, and the extreme variation in the electrical conductivity which is typically 5-7 orders of magnitude larger than that of the background geology. In this paper, we adopt a cylindrical discretization for numerical simulations to investigate three important aspects of DC resistivity in settings with steel-cased wells. (1) We examine the feasibility of using a surface-based DC resistivity survey for diagnosing impairments along a well in a casing integrity experiment. This parameter study demonstrates the impact of the background conductivity, the conductivity of the casing, the depth of the flaw, and the proportion of the casing circumference that is compromised on amplitude of the secondary electric fields measured at the surface. (2) Next, we consider elements of survey design for exciting a conductive or resistive target at depth. We show that conductive targets generate stronger secondary responses than resistive targets, and that having an electrical connection between the target and well can significantly increase the measured secondary responses. (3) Finally, we examine common strategies for approximating the fine-scale structure of a steel cased well with a coarse-scale representation to reduce computational load. We show that for DC resistivity experiments, the product of the conductivity and the cross-sectional area of the casing is the important quantity for controlling the distribution of currents and charges along its length. To promote insight into the physics, we present results by plotting the currents, charges, and electric fields in each of the scenarios examined. All of the examples shown in this paper are built on open-source software and are available as Jupyter notebooks.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Dispersion curve inversion is one of the core components of Rayleigh wave surveys, which mainly include linear and nonlinear inversion theoretical systems. Damped least squares (DLS) is the most mature and commonly used method of linear optimization, but it relies heavily on more accurate initial models, otherwise it can easily fall into a local minimum or can even result in an incorrect inversion. As a representative method of nonlinear optimization, genetic algorithm (GA) may be more feasible to obtain a global optimal solution for the geophysical inversion in theory. However, the GA algorithm is less stable, as well as less efficient in the later period of the inversion. In the past, the above two systems have been used independently to perform inversion processing. Faced with complex seismic geological conditions, they often display poor adaptability and lack balance between speed and accuracy. For this reason, we made a reasonable and effective improvement to the generation of the initial population and the coding of the classic GA to overcome the time-consuming and memory-intensive computational issues. We redefined the selection and crossover function to prevent the ‘premature convergence’ phenomenon in genetic iterations. Simultaneously, the DLS and the steepest descent method (SD) are embedded in the GA inversion process to linearly optimize the dominant individuals in the population (i.e. local extrema) in the local space to guide the population to move quickly and stably advance towards the global optimal direction. Next, a robust DLS inversion is used to obtain the final S-wave velocity model using the adaptive GA inversion results as the input velocity model. The model and actual dataset processing results show that our proposed nested joint inversion can combine the advantages of linear and nonlinear inversion that can effectively suppress the multi-solution problem of the dispersion curve inversion and significantly improve the inversion efficiency and accuracy.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Despite their importance as a fundamental constraint on Earth properties, regional-scale measurements of body-wave seismic attenuation are scarce. This is partially a result of the difficulty in producing robust estimates of attenuation. In this paper, we focus on measuring differential attenuation on records of teleseismic P waves. We examine a unique dataset of 5 records of the North Korean nuclear test of 2017 measured at five broadband seismic stations deployed within a few meters of each other but using different installation procedures. Given their extreme proximity, we expect zero differential intrinsic attenuation between the different records. However, we find that different attenuation measurement methods and implementation parameters in fact produce significant apparent differential attenuation (Δt*). Frequency-domain methods yield a wide range of Δt* estimates between stations, depending on measurement bandwidth and nuances of signal processing. This measurement instability increases for longer time windows. Time domain methods are largely insensitive to the frequency band being considered but are sensitive to the time window that is chosen. We determine that signal-generated noise can affect measurements in both the frequency and time domain. In some cases, the range of results amounts to a significant fraction of the range of differential attenuation across the conterminous United States as determined by a recent study. We suggest some approaches to manage the inherent instability in these measurements and recommend best practices to confidently estimate body wave attenuation.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉One of the largest induced seismic events attributed to hydraulic-fracturing operations, with a moment magnitude (〈span〉Mw〈/span〉) of 4.1, occurred on 12 January 2016 in Alberta, Canada. The hydraulic-fracturing treatment was monitored using a dense array of 93 shallow wellbores, each containing three geophones, producing a high-resolution dataset that is presented here as a case study. Over 9700 smaller events were observed both before and after the event, with a magnitude of completeness of approximately 〈span〉Mw〈/span〉 -1.9. The use of perforation shots with known locations enabled refinement of the velocity model, such that event hypocentres could be located with high accuracy and precision. Most of the hypocentres of induced events were located within the sedimentary section, 100s of metres above the treatment zone. Event locations of the foreshocks and aftershocks reveal fault activation along both small and large structures. The mechanisms and locations suggest a complex north - south oriented strike-slip fault system that we interpret to be part of a regional flower structure. The 〈span〉b〈/span〉-value for the dataset is close to unity, as expected for fault activation events; however, the three largest magnitude events deviate significantly from the fit, suggesting that it would be difficult to forecast these larger events based on the Gutenberg-Richter relation. This has significant implications for the mitigation of hydraulic-fracturing induced seismicity. The three largest events appear to be located on fault segments that are less optimally oriented for slip than those imaged by other event clusters, suggesting that larger events require a larger perturbation to the stress state in order to nucleate. Monitoring setups with the capability of the system presented here, if operated in real time, could be of great value to discern fault activation due to seismicity aligning on planar structures prior to a mainshock. This approach may have considerable potential to improve short-term forecasts of induced seismicity.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The maximum depth of magnetization in the Earth’s crust is generally thought to coincide with the Curie temperature of magnetite (580°C), and is commonly called Curie depth. A popular approach to estimating Curie depth is based on comparing the power-spectral density (PSD) of total-field magnetic anomaly data with theoretical expressions that describe the statistical properties of crustal magnetization and its depth distribution. However, the self-affine nature of magnetization (characterized by a power-law exponent in the PSD) renders this problem difficult in practice and several approximate methods have been devised to obtain estimates of the depth distribution and infer Curie depth variations using curve fitting of limited wavenumber ranges of biased PSD estimates. In most studies, however, errors are not carefully propagated through the estimation process and uncertainties in estimated parameters are not reported. In addition, robust spectral estimation techniques are required to map the spectral properties in space at the highest resolution possible and obtain Curie depth maps. Here we investigate the combined use of the robust multitaper spectral estimation technique with a Bayesian formalism to evaluate the recovery of the depth and statistical properties of the buried magnetized layer using a systematic analysis of synthetically-generated data. Based on these tests, we come to the conclusion that most Curie depth estimates are not reliable, although it may be possible to obtain more reliable estimates through incorporation of prior information on one or all model parameters. Finally, we suggest that the spectral analysis of magnetic anomaly data may be better suitable to hotter settings, where shallow Curie depth estimates are more robust and in the oceans, where magnetization is likely to be uniform. These studies may also be more useful in retrieving the power-law exponent of magnetization as well as the depth to the top of the magnetized layer.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The engines of surface deformation in the Anatolia–Aegean region are a matter of debate, including the origin of the high elevations of the Anatolian plateau. Recent publications based on geological and thermomechanical modelling emphasize the role of dynamic topography in the plateau uplift. However, quantitative estimates of the contribution of dynamic topography are affected by large uncertainties due to insufficient knowledge of the crustal structure, in particular crustal thickness and density. To reduce these uncertainties, we provide a new accurate crustal thickness map of the Anatolia–Aegean domain computed from a large volume of broadband seismic data. In addition, we display high-resolution seismic sections of the internal structure of the crust in Western and Central Anatolia. Density contrasts are derived from the same seismic data set and Bouguer gravity anomaly computed from the EGM2008 model. Our crustal thickness model is highly correlated with the topography suggesting that the Anatolian plateau is close to isostatic equilibrium. The average density difference between crust and upper mantle computed from our crustal model and Bouguer gravity anomaly is low compared to the global average, ∼0.315 ×10〈sup〉3〈/sup〉 kg m〈sup〉−3〈/sup〉. The ratio of surface elevation to crustal thickness is lower than average, 1/9.4, which also indicates a low-density crust. Differences between isostatic topography and observed topography are overall small (〈500 m). The east-to-west gradients of crustal thickness and topography changes are nearly constant in between the Taurides and Pontides at the northern and southern borders of Anatolia. The observed constant crustal thickness gradient may indicate a low viscosity lower crust supported by the thin mantle lithosphere evidenced by seismic tomography beneath the Anatolian plateau. We propose that viscous flow in the lower crust has smoothed out lateral changes in the crustal structure expected for such a heterogeneous collage of continental fragments. This flow may originate from gravitational potential energy differences between Eastern Anatolia (thick crust, high elevations) and the Aegean Sea (thin crust, low elevations), suggesting that gravity plays an integral part in the westward escape of Anatolia.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉A promising way to perform seismological studies in the Arctic region is deploying seismic stations on ice floes. The pioneering works by the Alfred Wegener Institute Bremerhaven have demonstrated the efficiency of such floating networks to explore local and regional seismicity and to build 3-D seismic models. However, problems remain, related to the identification of different types of seismic waves, particularly 〈span〉S〈/span〉 waves. Here, we perform 2-D and 3-D numerical simulations of seismic waves emitted by an earthquake to explore the possibility of recording different phases on the sea surface. We use different types of simple shear source models, namely strike-slip, vertical displacement and normal faults. In the calculated wave field, we obtain three major types of seismic waves recorded on the sea surface—〈span〉Pw, Sw〈/span〉 and 〈span〉SPw〈/span〉 (〈span〉w〈/span〉 denotes an acoustic wave in the water layer)—and numerous multiple waves. The clarity of the recorded phases strongly depends on the type of wave, source mechanism, epicentral distance, thickness of the water layer and depth of the source. For example, the 〈span〉Pw〈/span〉 phase is clearest for the strike-slip mechanism, less clear for the normal fault and almost invisible for the vertical displacement. The 〈span〉Sw〈/span〉 phase is observable in all of these cases; however, it can be confused with the 〈span〉SPw〈/span〉 phase that arrives earlier. In addition, at some distances, the 〈span〉Sw〈/span〉 wave interferes with the multiple 〈span〉Pw2〈/span〉 wave and therefore is hardly detectable. In summary, the numerical simulations in a model with a water layer have demonstrated several non-obvious features of wave propagation that should be taken into account when analysing experimental data recorded on ice floes.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Numerical modelling of seismic waves is a method for simulating the propagation of waves in the Earth. The objective is to make a prediction of the seismogram when given an assumed structure of the subsurface. The diffusive-viscous theory can be used to describe the attenuation of seismic waves propagating in fluid-saturated rocks, and to study the relationship between the frequency dependence of reflections and fluid-saturation in a porous medium, since the generally used theories such as acoustic, elastic theory, among others, are unable to effectively characterize the subsurface of the earth. We derive the second-order Runge–Kutta finite-volume scheme for the diffusive-viscous wave equation and based on this scheme, we simulate the propagation of seismic waves in a fluid-saturated medium. The numerical results indicate that the propagating waves in the fluid-saturated media attenuate greatly when compared with those of the acoustic scenario.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We present a new approach to estimate 3-D seismic velocities along a target interface. This approach uses an artificial neural network trained with user-supplied geological and geophysical input features derived from both a 3-D seismic reflection volume and a 2-D wide-angle seismic profile that were acquired from the Galicia margin, offshore Spain. The S-reflector detachment fault was selected as the interface of interest. The neural network in the form of a multilayer perceptron was employed with an autoencoder and a regression layer. The autoencoder was trained using a set of input features from the 3-D reflection volume. This set of features included the reflection amplitude and instantaneous frequency at the interface of interest, time-thicknesses of overlying major layers and ratios of major layer time-thicknesses to the total time-depth of the interface. The regression model was trained to estimate the seismic velocities of the crystalline basement and mantle from these features. The ‘true’ velocities were obtained from an independent full-waveform inversion along a 2-D wide-angle seismic profile, contained within the 3-D data set. The autoencoder compressed the vector of inputs into a lower dimensional space, then the regression layer was trained in the lower dimensional space to estimate velocities above and below the targeted interface. This model was trained on 50 networks with different initializations. A total of 37 networks reached minimum achievable error of 2 per cent. The low standard deviation (〈300  m s〈sup〉−1〈/sup〉) between different networks and low errors on velocity estimations demonstrate that the input features were sufficient to capture variations in the velocity above and below the targeted S-reflector. This regression model was then applied to the 3-D reflection volume where velocities were predicted over an area of ∼400 km〈sup〉2〈/sup〉. This approach provides an alternative way to obtain velocities across a 3-D seismic survey from a deep non-reflective lithology (e.g. upper mantle) , where conventional reflection velocity estimations can be unreliable.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The strainmeter record observed at Isabella (ISA), California, for the 1960 Chilean earthquake (〈span〉M〈/span〉〈sub〉w〈/sub〉 = 9.5) is one of the most important historical records in seismology because it was one of the three records that provided the opportunity for the first definitive observations of free oscillations of the Earth. Because of the orientation of the strainmeter rod with respect to the back azimuth to Chile, the ISA strainmeter is relatively insensitive to G (Love) waves and higher order (order ≥ 6) toroidal modes, yet long-period G waves and toroidal modes were recorded with large amplitude on this record. This observation cannot be explained with the conventional low-angle thrust mechanism typical of great subduction-zone earthquakes and requires an oblique mechanism with half strike-slip and half thrust. The strain record at Ogdenburg, New Jersey, the Press–Ewing seismograms at Berkeley, California, and the ultra-long period displacement record at Pasadena, California, also support the oblique mechanism. We tested the performance of the ISA strainmeter using other events including the 1964 Alaskan earthquake and found no instrumental problems. Thus, the ISA observation of large G/R and toroidal/spheroidal ratios most likely reflects the real characteristics of the 1960 Chilean earthquake, rather than an observational artefact. The interpretation of the large strike-slip component is not unique, but it may represent release of the strike-slip strain that has accumulated along the plate boundary as a result of oblique convergence at the Nazca–South American plate boundary. The slip direction of the 2010 Chilean (Maule) earthquake ( 〈span〉M〈/span〉〈sub〉w〈/sub〉 = 8.8) is rotated by about 10° clockwise from the plate convergence direction suggesting that right-lateral strain comparable to that of an 〈span〉M〈/span〉〈sub〉w〈/sub〉 = 8.3 earthquake remained unreleased and accumulates near the plate boundary. One possible scenario is that the strike-slip strain accumulated over several great earthquakes like the 2010 Maule earthquake was released during the 1960 Chilean earthquake. If this is the case, we cannot always expect a similar behaviour for all the great earthquakes occurring in the same subduction zone and such variability needs to be considered in long-term hazard assessment of subduction-zone earthquakes.〈/span〉

Description: 〈span〉〈div〉Abstract〈/div〉Multi-phase reactive transport processes are ubiquitous in igneous systems. A challenging aspect of modelling igneous phenomena is that they range from solid-dominated porous to liquid-dominated suspension flows and therefore entail a wide spectrum of rheological conditions, flow speeds, and length scales. Most previous models have been restricted to the two-phase limits of porous melt transport in deforming, partially molten rock and crystal settling in convecting magma bodies. The goal of this paper is to develop a framework that can capture igneous system from source to surface at all phase proportions including not only rock and melt but also an exsolved volatile phase. Here, we derive an 〈span〉n〈/span〉-phase reactive transport model building on the concepts of Mixture Theory, along with principles of Rational Thermodynamics and procedures of Non-equilibrium Thermodynamics. Our model operates at the macroscopic system scale and requires constitutive relations for fluxes within and transfers between phases, which are the processes that together give rise to reactive transport phenomena. We introduce a phase- and process-wise symmetrical formulation for fluxes and transfers of entropy, mass, momentum, and volume, and propose phenomenological coefficient closures that determine how fluxes and transfers respond to mechanical and thermodynamic forces. Finally, we demonstrate that the known limits of two-phase porous and suspension flow emerge as special cases of our general model and discuss some ramifications for modelling pertinent two- and three-phase flow problems in igneous systems.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉An integral component of the surface nuclear magnetic resonance forward model involves predicting the magnitude of the transverse magnetization following excitation. To predict the transverse magnetization, the Bloch equation must be solved. Traditional surface NMR forward models solve a simplified version of the Bloch equation where the relaxation terms are neglected. A shortcoming of this approach is that it can struggle to accurately describe the impact of relaxation during pulse effects. To address this concern, an alternative forward model based on solution of the full-Bloch equation is proposed. The advantage of the proposed scheme is that it implicitly accounts for relaxation during pulse effects, increases the flexibility to implement alternative parametrizations of the inverse model, and can readily describe an arbitrary excitation protocol given that it no longer requires closed form expressions of the transverse magnetizations. To demonstrate the potential of the updated forward modelling scheme, a novel approach for the inversion of complex-valued free-induction decay (FID) data is presented. The inverse model is reparametrized in order to produce depth profiles of the water content, T〈sub〉2〈/sub〉* and T〈sub〉2〈/sub〉. This approach has great potential to enhance the ability of FID measurements to provide insights into pore size and permeability as it can provide direct sensitivity to T〈sub〉2〈/sub〉. In contrast, traditional approaches that employ a forward model based on the simplified Bloch equation and estimate only T〈sub〉2〈/sub〉* are plagued by uncertainty surrounding the link between T〈sub〉2〈/sub〉* and pore size/permeability. Synthetic and field results are presented to demonstrate the feasibility of the proposed forward model and FID inversion framework.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉This study tests three models for generating stochastic earthquake-tsunami scenarios on subduction zones by comparison with deep ocean observations from 18 tsunamis during the period 2006–2016. It focusses on the capacity of uncalibrated models to generate a realistic distribution of hypothetical tsunamis, assuming the earthquake location, magnitude and subduction interface geometry are approximately known, while details of the rupture area and slip distribution are unknown. Modelling problems like this arise in tsunami hazard assessment, and when using historical and palaeo-tsunami observations to study pre-instrumental earthquakes. Tsunamis show significant variability depending on their parent earthquake’s properties, and it is important that this is realistically represented in stochastic tsunami scenarios. To clarify which aspects of earthquake variability should be represented, three scenario generation approaches with increasing complexity are tested: a simple fixed-area-uniform-slip (FAUS) model with earthquake area and slip deterministically related to moment magnitude; a variable-area-uniform-slip (VAUS) model which accounts for earthquake area variability; and a heterogeneous-slip (HS) model which accounts for both earthquake area variability and slip heterogeneity. The models are tested using deep-ocean tsunami time-series from 18 events (2006–2016) with moment magnitude 〈span〉M〈/span〉〈sub〉w〈/sub〉 〉 7.7. For each model and each observed event a ‘corresponding family of model scenarios’ is generated which includes stochastic scenarios with earthquake location and magnitude similar to the observation, with no additional calibration. For an ideal model (which perfectly characterizes the variability of tsunamis) the 18 observed events should appear like a random sample of size 18, created by taking one draw from each of the 18 ‘corresponding family of model scenarios’. This idea forms the basis of statistical techniques to test the models. First a goodness-of-fit criterion is developed to identify stochastic scenarios ‘most similar’ to the observed tsunamis, and assess the capacity of different models to produce good-fitting scenarios. Both the HS and VAUS models show similar capacity to generate tsunamis matching observations, while the FAUS model performs much more poorly in some cases. Secondly, the observed tsunami stage ranges are tested for consistency with the null hypothesis that they were randomly generated by the model. The null hypothesis cannot be rejected for the HS model, whereas both uniform-slip models exhibit a statistically significant tendency to produce small tsunamis too often. Finally, the statistical properties of slip for stochastic earthquake scenarios are compared against earthquake scenarios that best fit the observations. For the VAUS models the best-fitting model scenarios have higher slip on average than the stochastic scenarios, highlighting biases in this model. The techniques developed in this study can be applied to test stochastic tsunami scenario generation techniques, identify and partially correct their biases, and provide better justification for their use in applications.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Joint inversion of multiple electromagnetic data sets, such as controlled source electromagnetic and magnetotelluric data, has the potential to significantly reduce uncertainty in the inverted electrical resistivity when the two data sets contain complementary information about the subsurface. However, evaluating quantitatively the model uncertainty reduction is made difficult by the fact that conventional inversion methods—using gradients and model regularization—typically produce just one model, with no associated estimate of model parameter uncertainty. Bayesian inverse methods can provide quantitative estimates of inverted model parameter uncertainty by generating an ensemble of models, sampled proportional to data fit. The resulting posterior distribution represents a combination of a priori assumptions about the model parameters and information contained in field data. Bayesian inversion is therefore able to quantify the impact of jointly inverting multiple data sets by using the statistical information contained in the posterior distribution. We illustrate, for synthetic data generated from a simple 1-D model, the shape of parameter space compatible with controlled source electromagnetic and magnetotelluric data, separately and jointly. We also demonstrate that when data sets contain complementary information about the model, the region of parameter space compatible with the joint data set is less than or equal to the intersection of the regions compatible with the individual data sets. We adapt a trans-dimensional Markov chain Monte Carlo algorithm for jointly inverting multiple electromagnetic data sets for 1-D earth models and apply it to surface-towed controlled source electromagnetic and magnetotelluric data collected offshore New Jersey, USA, to evaluate the extent of a low salinity aquifer within the continental shelf. Our inversion results identify a region of high resistivity of varying depth and thickness in the upper 500 m of the continental shelf, corroborating results from a previous study that used regularized, gradient-based inversion methods. We evaluate the joint model parameter uncertainty in comparison to the uncertainty obtained from the individual data sets and demonstrate quantitatively that joint inversion offers reduced uncertainty. In addition, we show how the Bayesian model ensemble can subsequently be used to derive uncertainty estimates of pore water salinity within the low salinity aquifer.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Earth observation satellites yield a wealth of data for scientific, operational and commercial exploitation. However, the redistribution of mass in the system Earth is not yet part of the standard inventory of Earth Observation (EO) data products to date. It is derived from the Gravity Recovery and Climate Experiment (GRACE) mission and its Follow-On mission (GRACE-FO). Among many other applications, mass redistribution provides fundamental insights into the global water cycle. Changes in continental water storage impact the regional water budget and can, in extreme cases, result in floods and droughts that often claim a high toll on infrastructure, economy and human lives. The initiative for a European Gravity Service for Improved Emergency Management (EGSIEM) established three different prototype services to promote the unique value of mass redistribution products for Earth Observation in general and for early-warning systems in particular. The first prototype service is a scientific combination service to derive improved mass redistribution products from the combined knowledge of the European GRACE analysis centres. Second, the timeliness and reliability of such products is a primary concern for any early-warning system and therefore EGSIEM established a prototype for a near real-time service that provides dedicated gravity field information with a maximum latency of 5 d. Third, EGSIEM established a prototype of a hydrological/early warning service that derives wetness indices as indicators of hydrological extremes and assessed their potential for timely scheduling of high-resolution optical/radar satellites for follow-up observations in case of evolving hydrological extreme events.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Coastal populations are impacted by relative sea level variations, which consist both of absolute sea level variations and of vertical land motions. This paper focuses on the Southwest and Central Pacific region, a recognized vulnerable region to sea level rise and where a large range of vertical land motion dynamics is observed. We analyse vertical displacement rates obtained from Global Navigation Satellite Systems (GNSS) by different analysis centres. We study the role played by modelled parameters, such as step discontinuities (due to equipment changes, earthquakes, etc.), in the position time-series analysis. We propose a new modelling approach based on a joint inversion of GNSS position time-series from different analysis centres. The final uncertainty on the vertical land motion rates is estimated as a combination of the uncertainty due to the GNSS data processing itself and the uncertainty due to the stability of the reference frame in which the GNSS data are expressed. We find that the dominant trend in the Southwest and Central Pacific is a moderate subsidence, with an average rate of −1.1 mm yr〈sup〉–1〈/sup〉, but significant variations are observed, with displacement rates varying from an uplift of 1.6 ± 0.3 mm yr〈sup〉–1〈/sup〉 to a subsidence of −5.4 ± 0.3 mm yr〈sup〉–1〈/sup〉. Taking into account the geodynamic context, we assess, for each station, the relevance of current estimates of linear vertical displacement rate and uncertainty for forecasting future coastal sea levels.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The Groningen gas field is one of the largest gas fields in Europe. The continuous gas extraction led to an induced seismic activity in the area. In order to monitor the seismic activity and study the gas field many permanent and temporary seismic arrays were deployed. In particular, the extraction of the shear wave velocity model is crucial in seismic hazard assessment. Local 〈span〉S〈/span〉-wave velocity-depth profiles allow us the estimation of a potential amplification due to soft sediments.Ambient seismic noise tomography is an interesting alternative to traditional methods that were used in modelling the 〈span〉S〈/span〉-wave velocity. The ambient noise field consists mostly of surface waves, which are sensitive to the 〈span〉S〈/span〉wave and if inverted, they reveal the corresponding 〈span〉S〈/span〉-wave structures.In this study, we present results of a depth inversion of surface waves obtained from the cross-correlation of 1 month of ambient noise data from four flexible networks located in the Groningen area. Each block consisted of about 400 3-C stations. We compute group velocity maps of Rayleigh and Love waves using a straight-ray surface wave tomography. We also extract clear higher modes of Love and Rayleigh waves.The 〈span〉S〈/span〉-wave velocity model is obtained with a joint inversion of Love and Rayleigh waves using the Neighbourhood Algorithm. In order to improve the depth inversion, we use the mean phase velocity curves and the higher modes of Rayleigh and Love waves. Moreover, we use the depth of the base of the North Sea formation as a hard constraint. This information provides an additional constraint for depth inversion, which reduces the 〈span〉S〈/span〉-wave velocity uncertainties.The final 〈span〉S〈/span〉-wave velocity models reflect the geological structures up to 1 km depth and in perspective can be used in seismic risk modelling.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We examine seismic anisotropy beneath eastern China by analysing shear wave splitting of teleseismic SKS, SKKS and PKS phases recorded at 402 permanent and temporary stations. In general, our fast directions are approximately oriented in the E–W direction. The delay times vary greatly between 0.5 and 1.9 s. The delay times at stations in the Sichuan basin are less than 0.8 s, which could explain the fossilized anisotropy in the lithospheric mantle and little deformation after the formation of stable block. The fast directions in the Qinling–Dabie orogenic belts orient NWW–SEE, E–W and NEE–SWW, and are roughly parallel to the strikes of orogens that may indicate that the anisotropy contribution comes mainly from the lithospheric mantle. The fast directions mostly oriented in the nearly E–W direction, are generally subparallel to the absolute plate motion directions beneath the North China Craton and South China Block, suggesting it is primarily related to drag induced by the asthenospheric flow.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Archaeointensity determinations on burnt archaeological material are complex and reliable data scarce, although this kind of material can be of great interest in archaeological investigations. With the goal of analysing the reliability of archaeointensity determinations, an interlaboratory comparison study has been performed combining different experimental protocols on present-day reproductions of Pre-Columbian Mesoamerican archaeological artefacts and two brick samples. Samples were baked in an original kiln from an artisan workshop in western Mexico. The ambient magnetic field at the site during the experiment was measured and continuous temperature data were recorded at four different positions in the kiln during the heating–cooling procedure.Archaeointensity determinations were carried out with four different methods at four different palaeomagnetic laboratories: Thellier–Coe (Burgos, Spain), microwave (Liverpool, UK), multispecimen (Morelia, Mexico) and multispecimen with the extended protocols for fraction and domain-state correction (Montpellier, France). 26 conventional resistive heating determinations with the Thellier–Coe protocol yielded a 100 per cent success rate, while 7 out of 8 microwave-heating determinations with the Thellier–Coe protocol also provided successful results. Also, two multispecimen determinations performed with both multispecimen methods provided statistically reliable results. In all cases, a good agreement between the determined archaeointensities and the ambient field at the production site could be observed.Highly reversible magnetization-versus-temperature curves yielded slightly Al, Mg or Ti-substituted magnetite as the main ferromagnetic (〈span〉s.l〈/span〉.) phase. In addition, in several samples, a thermally stable low Curie-temperature phase displaying a high coercivity behaviour could be observed in thermomagnetic curves and by thermal demagnetization of saturation isothermal remanent magnetization. This phase is interpreted as ε-Fe〈sub〉2〈/sub〉O〈sub〉3〈/sub〉. To our knowledge, its occurrence has never been reported through the experimental recreation of burnt archaeological materials. No correlation could be observed between the proxies of domain-state behaviour and deviation of palaeointensity determinations from the expected result.Results obtained on clay samples heated in this type of ancient kiln can be considered a good source for determining the geomagnetic field strength variation in the past. Matching palaeointensity results obtained with different methods based on different principles can be taken as a quality criterion for result reliability and consistency.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Recently an ambitious experiment combining deep seismic surveys from near-vertical and wide-angle acquisition methods was carried out in Brazil. The seismic lines are essentially coincident and crossed the Parnaíba Basin from west to east near latitude 5°S. Here, the wide-angle reflection and refraction (WARR) and deep seismic reflection (DSR) results, which were previously interpreted independently, are compared by directly correlating WARR interfaces converted to TWTT with the major reflective horizons identified in the zero-offset image and by considering coincident reflectivity patterns displayed in both data sets. This integrated WARR and DSR analysis allowed a spatial association of the apparently acoustically featureless crust imaged in the DSR profile to the high reflectivity observed in the WARR data. Numerical tests and elastic modelling show that variations of the elastic properties of the crust, particularly as they are characterized by low 〈span〉Vp〈/span〉 and 〈span〉Vs〈/span〉 contrasts with a possible increase of the 〈span〉Vp〈/span〉/〈span〉Vs〈/span〉 ratio, can only weakly explain the observed reflectivity patterns but that fine-scale lithological heterogeneity within the crust is capable of replicating the observed contrasting seismic responses. The segment of the Parnaíba Basin crust that is characterized by fine-scale lithological heterogeneity lies directly above a mafic crustal underplate defined by the WARR model and was named as the Grajaú domain on the basis of WARR-derived velocity model. The applied methodologies allow added value to be taken from the independent seismic data sets and provide new information about crustal structure that may have important implications for overlying intracontinental basin evolution.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The analysis of surface wave dispersion curves (DCs) is widely used for near-surface S-wave velocity (VS) reconstruction. However, a comprehensive characterization of the near-surface requires also the estimation of P-wave velocity (VP). We focus on the estimation of both VS and VP models from surface waves using a direct data transform approach. We estimate a relationship between the wavelength of the fundamental mode of surface waves and the investigation depth and we use it to directly transform DCs into VS and VP models in laterally varying sites. We apply the workflow to a real dataset acquired on a known test site. The accuracy of such reconstruction is validated by a waveform comparison between field data and synthetic data obtained by performing elastic numerical simulations on the estimated VP and VS models. The uncertainties on the estimated velocity models are also estimated.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We analyze 67 new conductive heat flow measurements on the southern flank of the Costa Rica Rift (CRR). Heat flow measurements cover five sites ranging in oceanic crustal age between approximately 1.6 and 5.7 Ma, and are co-located with a high-resolution multi-channel seismic line that extends from slightly north of the first heat flow site (1.6 Ma) to beyond ODP Hole 504B in 6.9 Ma crust. For the five heat flow sites, the mean observed conductive heat flow is ≈ 85 mWm〈sup〉−2〈/sup〉. This value is approximately 30 per cent of the mean lithospheric heat flux expected from a half-space conductive cooling model, indicating that hydrothermal processes account for about 70 per cent of the heat loss. The advective heat loss fraction varies from site to site and is explained by a combination of outcrop to outcrop circulation through exposed basement outcrops and discharge through faults. Super-critical convection in Layer 2A extrusives occurs between 1.6 and 3.5 Ma, and flow through a thinly-sedimented basement high occurs at 4.6 Ma. Advective heat loss diminishes rapidly between ≈ 4.5 and ≈ 5.7 Ma, which contrasts with plate cooling reference models that predict a significant deficit in conductive heat flow up to ages ≈ 65 ± 10 Ma. At ≈ 5.7 Ma the CRR topography is buried under sediment with an average thickness ≈ 150 m, and hydrothermal circulation in the basement becomes sub-critical or perhaps marginally critical. The absence of significant advective heat loss at ≈ 5.7 Ma at the CRR is thus a function of both burial of basement exposure under the sediment load and a reduction in basement permeability that possibly occurs as result of mineral precipitation and original permeability at the time of formation. Permeability is a non-monotonic function of age along the southern flank of the CRR, in general agreement with seismic velocity tomography interpretations that reflect variations in the degree of ridge-axis magma supply and tectonic extension. Hydrothermal circulation in the young oceanic crust at southern flank of CRR is affected by the interplay and complex interconnectedness of variations in permeability, sediment thickness, topographical structure, and tectonic and magmatic activities with age.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The application of the singular value decomposition method (SVD) for filtering of seismic data has become common in recent decades, as it promotes significant improvements of the signal-to-noise ratio, highlighting reflections in seismograms. One particular way to apply SVD in a single (or multivariate) time-series is the singular spectrum analysis (SSA) method, normally applied on constant-frequency slices in one or many spatial dimensions. We demonstrate that SSA method applied in the time domain corresponds to filtering the time-series with a symmetric zero-phase filters, which are the autocorrelations of the eigenvectors of the data covariance matrix, preserving the phase of the original data. In this paper, we explore the SSA method in the time domain, and we propose a new recursive-iterative SSA (RI-SSA) algorithm, which uses only the first eigenvector of the data covariance matrix to decompose a discrete time-series into signal components. From the analytic signal of each component we compute a time–frequency representation. By interpretation of the time signals and their time–frequency representations, groups with similar features are summed to produce a smaller number of signal components. The resulting RI-SSA signal decomposition is exact and phase-preserving, but non-unique. Applications to land seismic data for ground-roll removal and to two synthetic signals for time–frequency analysis give good results.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We derive two sets of new wave equations involving the fractional time derivatives or fractional Laplacians for simulating seismic wave propagation in viscoelastic anisotropic (VA) media based on the Kjartansson's constant-〈span〉Q〈/span〉 model. The approximate fractional Laplacian wave equation is developed under the assumptions of the small velocity anisotropy parameter $\delta$ and attenuation anisotropy parameter ${\delta _Q}$ (or weak velocity and attenuation anisotropy). Both the formulas have the advantages of simple form and can accurately describe the constant-〈span〉Q〈/span〉 (i.e. frequency-independent quality factor) attenuation and arbitrary attenuation anisotropy behaviours compared to the widely used VA theory based on the generalized standard linear solids (GSLS) model and memory variables. Under the assumption of homogeneous plane wave and the constant-〈span〉Q〈/span〉 attenuation mechanism, we further derive exact analytical expressions for the phase velocities, attenuation coefficients, quality factors of the 〈span〉P〈/span〉- and 〈span〉SV〈/span〉-waves, and analyse their direction dependence in different attenuation anisotropy cases. For numerical modelling, we implement the fractional finite-difference (FD) method with the Grünwald–Letnikov (GL) approximation to solve the fractional time wave equation, and the generalized Fourier pseudospectral (PS) method to solve the fractional Laplacian wave equation, respectively. The PS method is highly efficient compared with the fractional FD method since it avoids additional memory to store the past wavefields. Numerical results of the homogeneous VTI (transversely isotropic with a vertical symmetry axis) model validate the accuracy of the two proposed schemes and illustrate the influence of attenuation and attenuation anisotropy on seismic wavefields, which are consistent with theoretical analysis. Finally, the modelling of the Hess VTI model shows the applicability of our formulations and algorithms in heterogeneous media.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Different mechanisms have been proposed as explanations for seismic anisotropy at the base of the mantle, including crystallographic preferred orientation of various minerals (bridgmanite, post-perovskite and ferropericlase) and shape preferred orientation of elastically distinct materials such as partial melt. Investigations of the mechanism for D" anisotropy usually yield ambiguous results, as seismic observations rarely (if ever) uniquely constrain a mechanism or orientation and usually rely on significant assumptions to infer flow patterns in the deep mantle. Observations of shear wave splitting and polarities of SdS and PdP reflections off the D" discontinuity are among our best tools for probing D" anisotropy; however, currently available data sets cannot constrain one unique scenario among those suggested by the mineral physics literature. In this work, we determine via a forward modelling approach what combinations of body wave phases (e.g. SKS, SKKS and ScS) are required to uniquely constrain a mechanism for D" anisotropy. We test nine models based on single-crystal and polycrystalline elastic tensors provided by mineral physics studies. Our modelling predicts fast shear wave splitting directions for SKS, SKKS and ScS phases, as well as polarities of 〈span〉P-〈/span〉 and 〈span〉S-〈/span〉wave reflections off the D" interface, for a range of propagation directions, via solution of the Christoffel equation. We run tests using randomly selected synthetic data sets based on a given starting model, controlling the total number of measurements, the azimuthal distribution, and the type of seismic phases. For each synthetic data set, we search over all possible elastic tensors and orientations to determine which are consistent with the synthetic data. Overall, we find it difficult to uniquely constrain the mechanism for anisotropy with a typical number of seismic anisotropy measurements (based on currently available studies) with only one measurement technique (SKS, SKKS, ScS or reflection polarities). However, data sets that include SKS, SKKS and ScS measurements or a combination of shear wave splitting and reflection polarity measurements increase the probability of uniquely constraining the starting model and its orientation. Based on these findings, we identify specific regions (i.e. North America, northwestern Pacific and Australia) of the lowermost mantle with sufficient ray path coverage for a combination of measurement techniques.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Free air gravity anomalies are now available from Tibetan Plateau and surrounding regions. They show that the Plateau itself is isostatically compensated, and that its elevation is supported by its margins acting as barriers. Its interior has an elastic thickness of less than 4 km. The main features of the Plateau can be understood by using a simple lithostatic flexural model, which can account for the paired gravity anomalies and the thrust faulting on its margins. The resulting estimates of the magnitudes of flexural forces agree with those from other observations. The thick crust that underlies the Plateau generates radiogenic heat that greatly reduces its viscosity and accounts for its weakness. As Tibet flows south over cold Indian lithosphere, heat is conducted downwards, heating up the upper mantle beneath the Moho. Combined with isostatic compensation, the resulting thermal expansion then produces a gradient in crustal thickness, from about 75 km in the south to 65 km in the north. Though this lithostatic model provides a framework for understanding many of the features of Tibet, it does not account for the difference in the dynamical behaviour between the north and south of the Plateau.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The slip and static stress drop of the 2015 〈span〉M〈/span〉〈sub〉w〈/sub〉 7.8 Gorkha, Nepal earthquake has been studied using its excellent near-source static observations and a newly developed finite fault inversion algorithm with an average stress drop constraint. A series of nonlinear inversions with different target stress drops are conducted to search for the solution that not only most accurately fits the geodetic data, but also has an energy-based stress drop, $\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$, matching the target value. Our result reveals that the misfit to the geodetic data gradually decreases when the $\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$ of the inverted slip model increases from 2 to 7 MPa; but becomes nearly constant when $\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$ further increases. Hence, only the lower bound of $\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$, that is, ${\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}} ^{\mathrm{ min}}}$(∼7 MPa), of the Gorkha earthquake can be constrained with near-fault geodetic measurements, consistent with our previous study using just far-field seismic data. The upper bound of ${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$ can be constrained with an extra constraint on the maximum local stress drop. An artificial upper bound can also be reached if using a large subfault size to represent the source. The lower bound of $\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$ leads to the lower bound of the apparent available energy and the upper bound of the radiation efficiency (${\eta _\mathrm{ R}}$, 0.09–0.15), though the latter is also sensitive to the determination of the radiated seismic energy. We find that the lower ${\eta _\mathrm{ R}}$ and the reported high rupture velocity (80–93 per cent shear wave speed) can be reconciled by considering the aspect ratio of the dominant slip patch. We recommend using ${\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}} ^{\mathrm{ min}}}$ to replace various slip smoothing constraints to stabilize the finite fault inversion because of its clear physical meaning.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Study of induced, triggered, stimulated, and nuance earthquakes can provide a unique opportunity to probe into stress/frictional conditions on the subsurface fault. Although, mechanism of faulting is tectonically driven phenomenon, anthropogenic crustal (un-)loading process can also drive long-term fault slip or modulate seismicity in the nearby seismogenic crust. Using Gravity Recovery and Climate Experiment (GRACE) data, surface subsidence rate captured by cGPS, and three-dimensional hydrological stress models, we report that the November 12, 2017 M7.3 Iran-Iraq border earthquake in the Arabia-Eurasia collision boundary region of Zagros fold-and-thrust belt was significantly influenced by the anthropogenic groundwater unloading in the Euphrates and Tigris River plains and western Iran, a region which is probably the most intensely irrigated region in the Middle East. We find that the anthropogenic crustal unloading contributed about one tenth of the secular interseismic stress change due to tectonic process. Thus the resulting cumulative stress accumulation was significantly influenced by the extensive groundwater pumping.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Global mass redistribution between the Earth subsystems oceans, atmosphere, and continental hydrosphere cause a predominantly seasonal signal in Earth rotation excitation that superimposes the effects of each individual Earth subsystem. Especially for annual length-of-day variations a consistent consideration of the global mass balance among atmosphere, ocean, and continental water is necessary to compare the simulated effective angular momentum functions for Earth rotation from geophysical models with geodetic observations. In addition to atmospheric, oceanic, and hydrological contributions, we estimate the contributions due to the global mass balance effect using the new ESMGFZ SLAM product as well as estimates of the barystatic ocean bottom pressure anomalies from the GRACE Level 3 GravIS products. For the annual cycle the global mass balance effect overcompensates the contributions to length-of-day variations from terrestrial hydrology. Moreover, most of the atmospheric surface pressure contribution is also compensated. The global mass balance effect has to be calculated for each combination of geophysical Earth system models individually. Considering the global mass balance, model based mass induced excitation on seasonal length-of-day variations coincide well with estimates from satellite gravimetry. Moreover, the mass terms can be determined accurate enough to attribute the remaining gap in the length-of-day excitation budget between models and observation clearly to an underestimation of atmospheric wind speeds in the global European weather forecast model. Magnifying its wind speeds by +7% the sum of all ESMGFZ angular momentum functions can almost perfectly explain the total length-of-day excitation.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Mascon solutions, one of the recent advances in Gravity Recovery and Climate Experiment (GRACE) products, are being widely used due to their significant improvement in spatial accuracy and their ability to be directly employed without postprocessing. However, few studies have verified the performance of mascon solutions at different regional scales, especially those smaller than GRACE's native resolution. We investigated the consistency between mascon and spherical harmonic solutions (hereafter named mascon and SH solutions) as release-05 products from the Center for Space Research (CSR) and determined whether mascon solutions had a higher resolution than SH solutions (approximately 300 km). In addition, we used ICESat, LEGOS, and DAHITI observations of water levels as independent supplementary data to validate the comparison results. We considered three cases of mass signals with varying spatial scales to examine the performance of mascon solutions in specific regions: 1) the abrupt water-level change in Longyang Reservoir (an area smaller than 1° × 1°) driven by heavy rainfall; 2) mass increase in the middle and lower Yangtze River basin; and 3) mass loss of Tianshan Mountain glaciers and their subregions, including Bosten Lake at the foot of the glaciers. We concluded that mascon solutions have the potential to reveal gravity signals in areas larger than 3° × 3° from mass sources of Terrestrial Water Storage anomalies (TWSA), whereas they perform poorer than gridded SH products (GSH) or SH solutions with a decorrelation filter (DDK4) in some specific regions smaller than 3° × 3°. The results of the three different cases provide references to non-geodetic users to select suitable GRACE solutions for different applications.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Seismic waves radiated from an earthquake induce mass redistribution in the Earth, which in turn generates prompt gravity perturbations at all distances before the arrival of seismic waves. Here we obtain an analytic solution for the transient elastic response to prompt gravity perturbations induced by a portion of earthquake mass redistribution in an infinite homogeneous isotropic non-self-gravitation elastic medium. This solution reveals the fundamental nature of an unusual net-inner-force-free deformation regime that emerges between the onset of an earthquake and P-wave arrival, which we find inherent in action at a distance with the inverse-square law. These results provide physical insights into why a gravimeter installed in the Earth is characterized by significantly decreased sensitivity to earthquake-induced prompt gravity perturbations.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Microseismic data is usually of low signal-to-noise ratio (SNR), which makes it difficult to utilize the microseismic waveforms for imaging and inversion. We develop a useful denoising algorithm based on a non-stationary least-squares decomposition model to enhance the quality of microseismic signals. The microseismic signals are assumed to be represented by a superposition of several smoothly variable components. We construct a least-squares inverse problem to solve for the the smooth components. We constrain the least-squares inversion via both time and space constraints. The temporal smoothness constraint is applied to ensure the stability when calculating the non-stationary autoregression coefficients. The space-smoothness constraint is applied to extract the spatial correlation among multi-channel microseismic traces. The new algorithm is validated via several synthetic and real microseismic data and is proved to be effective. Comparison with the state-of-the-art algorithms demonstrates that the proposed method is more powerful in suppressing random noise of a wide range of levels than its competing methods.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉One of the largest induced seismic events attributed to hydraulic-fracturing operations, with a moment magnitude (〈span〉M〈/span〉〈sub〉w〈/sub〉) of 4.1, occurred on 12 January 2016 in Alberta, Canada. The hydraulic-fracturing treatment was monitored using a dense array of 93 shallow wellbores, each containing three geophones, producing a high-resolution data set that is presented here as a case study. Over 9700 smaller events were observed both before and after the event, with a magnitude of completeness of approximately 〈span〉M〈/span〉〈sub〉w〈/sub〉 −1.9. The use of perforation shots with known locations enabled refinement of the velocity model, such that event hypocentres could be located with high accuracy and precision. Most of the hypocentres of induced events were located within the sedimentary section, 100s of metres above the treatment zone. Event locations of the foreshocks and aftershocks reveal fault activation along both small and large structures. The mechanisms and locations suggest a complex north–south oriented strike-slip fault system that we interpret to be part of a regional flower structure. The 〈span〉b〈/span〉-value for the data set is close to unity, as expected for fault activation events; however, the three largest magnitude events deviate significantly from the fit, suggesting that it would be difficult to forecast these larger events based on the Gutenberg–Richter relation. This has significant implications for the mitigation of hydraulic-fracturing induced seismicity. The three largest events appear to be located on fault segments that are less optimally oriented for slip than those imaged by other event clusters, suggesting that larger events require a larger perturbation to the stress state in order to nucleate. Monitoring setups with the capability of the system presented here, if operated in real time, could be of great value to discern fault activation due to seismicity aligning on planar structures prior to a main shock. This approach may have considerable potential to improve short-term forecasts of induced seismicity.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Quantifying spatial and temporal changes in thermospheric neutral density is important for various applications such as precise orbit determination, estimating mission lifetime and re-entry prediction of Earth orbiting objects. It is also crucial for analysis of possible collisions between active satellite missions and space debris. Empirical models are frequently applied to estimate neutral densities at the position of satellites. But their accuracy is severely constrained by model simplifications and the sampling limitation of solar and geomagnetic indices used as inputs. In this study, we first estimate thermospheric neutral density by processing the high-accuracy accelerometer measurements on-board of the twin-satellite mission Gravity Recovery And Climate Experiment (GRACE). Daily density corrections (in terms of scales) are then computed for the commonly used NRLMSISE-00 empirical model. The importance of these daily scales is examined within an orbit determination practice. Finally, three data-driven prediction techniques based on Artificial Neural Network (ANN) are applied to forecast the daily density corrections for few days to months. Our numerical results indicate that GRACE derived scales are correlated with solar and geomagnetic indices and can improve the timing (from few hours to days) and magnitude of model simulations (up to 10–100 times) during high solar or geomagnetic activity when they usually perform poorly. We found that the Non-linear Autoregressive with Exogenous (External) Input (NARX) ANN technique performs well in predicting the corrections with an average fit of 0.8 or more in terms of squared correlation coefficients for time-scales of 7–90 days.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉It is difficult to separate additive random noise from spatially coherent signal using a rank-reduction (RR) method that is based on the truncated singular value decomposition (TSVD) operation. This problem is due to the mixture of the signal and the noise subspaces after the TSVD operation. This drawback can be partially conquered using a damped RR (DRR) method, where the singular values corresponding to effective signals are adjusted via a carefully designed damping operator. The damping operator works most powerfully in the case of a small rank and a small damping factor. However, for complicated seismic data, e.g. multichannel reflection seismic data containing highly curved events, the rank should be large enough to preserve the details in the data, which makes the DRR method less effective. In this paper, we develop an optimal damping strategy for adjusting the singular values when a large rank parameter is selected so that the estimated signal can best approximate the exact signal. We first weight the singular values using optimally calculated weights. The weights are theoretically derived by solving an optimization problem that minimizes the Frobenius-norm difference between the approximated and the exact signal components. The damping operator is then derived based on the initial weighting operator to further reduce the residual noise after the optimal weighting. The resulted optimally damped rank-reduction method is nearly an adaptive method, i.e. insensitive to the rank parameter. We demonstrate the performance of the proposed method on a group of synthetic and real 5-D seismic data.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Array-based techniques for multimode phase speed measurements of surface waves and modal waveform decomposition are developed and applied to USArray. Measuring phase speeds of multimode surface waves is not a straightforward issue since wavetrains of several mode-branches propagate with similar group speeds and overlap each other in an observed seismogram. In this study, we develop a series of array-based methods for (1) multimode phase speed measurements based on a classical 〈span〉f–k〈/span〉 analysis using a linear array, and (2) modal waveform decomposition for the centroid location of the array using the linear Radon transform. These methods are applied to synthetic and observed waveforms to check their validity and utility. A series of synthetic experiments revealed that the precise measurement of multimode phase speeds requires a long linear array over 2000–3000 km. Extracted dispersion curves reflect an average structure at the centroid of the array. Decomposed modal waveforms at the centroid location of an array match well with theoretically predicted waveforms. We then apply our method to seismograms observed at USArray stations. Decomposed mode-branch waveforms are used to measure interstation phase speeds through a nonlinear waveform fitting method, suggesting that we can achieve stable measurements of local phase speeds even for the fundamental mode Love wave that tends to be overlapped and interfered by higher modes. Our hybrid array-based approach can be of help in further enhancing the horizontal and vertical resolutions of regional tomography models by making the most of highly dense broad-band seismic networks.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉In coseismic slip distribution inversion, the Laplacian smoothness constraint is used to avoid the rank deficiency of the coefficient Green matrix and ensure the smoothness among the patches of the fault plane. However, by introducing the classic Laplacian smoothness constraint, the inversion of the maximum slip is commonly underestimated. To better solve this problem, here we propose a new method of weighting to the Laplacian second-order smoothness matrices by the slip. The Laplacian smoothness constraint with unequal weights and the classic Laplacian constraint are used to carry out the simulation experiments. The simulation experiments show that the accuracy of the estimated maximum slip from our method is outperforming the classical Laplacian method by 12–19 per cent. Moreover, the estimations of the average slip and moment magnitude have been enhanced, which improved ranging from 4 to 12.5 per cent and 0.4 to 9 per cent over the inversion results of classic Laplacian smoothness constraint, respectively. In order to further validate the general applicability of the proposed method, we conduct the experiments in terms of the inversion of L'Aquila and Taiwan Meinong earthquakes. The maximum slip inversion results show that the inversion of Laplacian smoothness constraint with unequal weights are larger than that of the classic Laplacian smoothness constraint, which is consistent with the conclusion of simulation experiments. In addition, the parameters inversion results of the actual earthquake from the proposed method are in agreement with other studies which also indicate the feasibility and effectiveness of our method.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Debate is ongoing as to which tectonic model is most consistent with the known geology of southeast Australia, formerly part of the eastern margin of Gondwana. In particular, numerous tectonic models have been proposed to explain the enigmatic geological relationship between Tasmania and the mainland, which is separated by Bass Strait. One of the primary reasons for the lack of certainty is the limited exposure of basement rocks, which are masked by the sea and thick Mesozoic–Cenozoic sedimentary and volcanic cover sequences. We use ambient noise tomography recorded across Bass Strait to generate a new shear wave velocity model in order to investigate crustal structure. Fundamental mode Rayleigh wave phase velocity dispersion data extracted from long-term cross-correlation of ambient noise data are inverted using a transdimensional, hierarchical, Bayesian inversion scheme to produce phase velocity maps in the period range 2–30 s. Subsequent inversion for depth-dependent shear wave velocity structure across a dense grid of points allows a composite 3-D shear wave velocity model to be produced. Benefits of the transdimensional scheme include a data-driven parametrization that allows the number and distribution of velocity unknowns to vary, and the data noise to also be treated as an unknown in the inversion. The new shear wave velocity model clearly reveals the primary sedimentary basins in Bass Strait as slow shear velocity zones which extend down to 14 km in depth. These failed rift basins, which formed during the early stages of Australia–Antarctica break-up, appear to be overlying thinned crust, where high velocities of 3.8–4.0 km s〈sup〉−1〈/sup〉 occur at depths greater than 20 km. Along the northern margin of Bass Strait, our new model is consistent with major tectonic boundaries mapped at the surface. In particular, we identify an east dipping velocity transition zone in the vicinity of the Moyston Fault, a major tectonic boundary between the Lachlan and Delamerian orogens, which are part of the Phanerozoic accretionary terrane that makes up eastern Australia. A pronounced lineament of high shear wave velocities (∼3.7–3.8 km s〈sup〉−1〈/sup〉) in the lower crust of our new model may represent the signature of relict intrusive magmatism from failed rifting in the early stages of Australia–Antarctica break-up along a crustal scale discontinuity in the Selwyn Block microcontinent which joins Tasmania and Victoria.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Finite-difference (FD) methods are widely used for numerical solution of acoustic and elastic wave equations. Temporal high-order FD methods exhibit better accuracy and stability than the methods with second-order differencing in time. Also, the implicit calculation of spatial derivatives can bring significant improvement in accuracy. The present implicit FD methods with high-order accuracy in time are based on centred grids. In this paper, we propose an implicit staggered-grid FD (SFD) scheme with a combined stencil, which is the combination of rhombus/pyramid and cross stencils in 2-D/3-D case, for modelling scalar wave propagation. Our scheme computes the temporal and spatial derivatives using high-order temporal and implicit spatial FD operators based on the combined stencil and the conventional stencil, respectively. We derive the dispersion relations of the FD scheme for 2-D and 3-D scalar wave equations and estimate temporal and implicit spatial FD coefficients by Taylor series expansion (TE) and least squares (LS). According to the kinds of FD coefficients, we formulate four implicit SFD operators: TE–TE, TE–LS, LS–TE and LS–LS operators. We carry out the comparison between our scheme and several existing SFD schemes: the conventional explicit and implicit, optimal explicit and implicit and explicit temporal high-order schemes. 2-D and 3-D dispersion analysis, stability analysis and modelling examples reveal that our implicit scheme has greater accuracy than other schemes and requires slightly stricter stability condition than the explicit temporal high-order SFD schemes. Owing to higher accuracy, our implicit SFD scheme with LS-LS operators allows for shorter FD operators and larger grid spacing, which can increase the computational efficiency.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Knowledge of deformation at plate boundaries has been improved greatly by the development of observational techniques in space geodesy. However, most theoretical and numerical models of coseismic deformation have remained very simple and do not include realistic Earth structure. 3-D material heterogeneity and topography are often neglected because simple models are assumed to be sufficient and available tools cannot easily accommodate complex heterogeneity. In this study, we demonstrate the importance of 3-D heterogeneity using a spectral element method that incorporates topography and 3-D material properties. Using a parabolic hill model and a topographic model of the 2010 Maule earthquake, we show that topographic features can alter the shape of observed surface deformation patterns. We also estimate the coseismic surface deformation due to a model of the slip distribution of the 2015 Gorkha earthquake using realistic topography and 3-D elastic structure, and find that the presence of topography causes changes in the shape of observed surface displacement patterns, while material heterogeneities primarily affect the magnitude of observed displacements. Our results show that the inclusion of topography in particular can affect predictions of coseismic deformation modelling.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Dynamo action in the Earth’s outer core is expected to be controlled by a balance between pressure, Coriolis, buoyancy and Lorentz forces, with marginal contributions from inertia and viscous forces. Current numerical simulations of the geodynamo, however, operate at much larger inertia and viscosity because of computational limitations. This casts some doubt on the physical relevance of these models. Our work aims at finding dynamo models in a moderate computational regime which reproduce the leading-order force balance of the Earth. By performing a systematic parameter space survey with Ekman numbers in the range 10〈sup〉−6〈/sup〉 ≤ 〈span〉E〈/span〉 ≤ 10〈sup〉−4〈/sup〉, we study the variations of the force balance when changing the forcing (Rayleigh number, 〈span〉Ra〈/span〉) and the ratio between viscous and magnetic diffusivities (magnetic Prandtl number, 〈span〉Pm〈/span〉). For dipole-dominated dynamos, we observe that the force balance is structurally robust throughout the investigated parameter space, exhibiting a quasi-geostrophic (QG) balance (balance between Coriolis and pressure forces) at zeroth order, followed by a first-order MAC balance between the ageostrophic Coriolis, buoyancy and Lorentz forces. At second order this balance is disturbed by contributions from inertia and viscous forces. Dynamos with a different sequence of the forces, where inertia and/or viscosity replace the Lorentz force in the first-order force balance, can only be found close to the onset of dynamo action and in the multipolar regime. To assess the agreement of the model force balance with that expected in the Earth’s core, we introduce a parameter quantifying the distance between the first- and second-order forces. Analysis of this parameter shows that the strongest-field dynamos can be obtained close to the onset of convection (〈span〉Ra〈/span〉 close to critical) and in situations of reduced magnetic diffusivity (high 〈span〉Pm〈/span〉). Decreasing the Ekman number gradually expands this regime towards higher supercriticalities and lower values of 〈span〉Pm〈/span〉. Our study illustrates that most classical numerical dynamos are controlled by a QG-MAC balance, while cases where viscosity and inertia play a dominant role are the exception rather than the norm.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The Amatrice–Norcia–Visso sequence is characterized by complex behaviour that is somewhat atypical of main-shock–aftershock sequences, as there were multiple large main shocks that continued for months.In this study we focus on the Amatrice sequence (main shock 2016 August 24, 〈span〉M〈/span〉〈sub〉w〈/sub〉 = 5.97) to evaluate the apparent stress values and magnitude-dependent scaling in order to improve our knowledge of processes that control small and large earthquakes within this active region of Italy. Apparent stress is proportional to the ratio of radiated seismic energy and seismic moment, and as such, these stress parameters play an important role in hazard prediction as they have a strong effect on the observed and predicted ground shaking.We analyse 83 events of the sequence from 2016 August 24 to October 16, within a radius of 20 km from the main shock and with an 〈span〉M〈/span〉〈sub〉w〈/sub〉 ranging between 5.97 and 2.72.Taking advantage of the averaging nature of coda waves, we analyse coda-envelope-based spectral ratios between neighbouring event pairs. We use equations proposed by Walter 〈span〉et al.〈/span〉 to consider stable, low-frequency and high-frequency spectral ratio levels which provide measures of the corner frequency and apparent stress ratios of the events within the sequence.The results demonstrate non-self-similar behaviour within the sequence, suggesting a change in dynamics between the largest events and the smaller aftershocks.The apparent stress and corner frequency estimates are compared to those obtained by Malagnini and Munafò who utilized hundreds of direct 〈span〉S〈/span〉-wave spectral ratio measurements to obtain their results. Although our analysis is based only on 83 events, our results are in very good agreement, demonstrating once more that the use of coda waves is very stable and provides lower variance measures than those using direct waves.A comparison with recent Central Apennines source scaling models derived from various seismic sequences (1997–1998 Colfiorito, 2002 San Giuliano di Puglia, 2009 L'Aquila) shows that the Amatrice sequence source scaling in this study is well represented by the models proposed by Pacor 〈span〉et al.〈/span〉 and Malagnini and Mayeda.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We applied a nonlinear teleseismic tomography algorithm to explore the 3-D structure of the upper mantle beneath the Arabia–Eurasia continental collision zone in western Iran. An unprecedented data set consisting of 32 738 teleseismic 〈span〉P-〈/span〉wave relative arrival time residuals from 129 permanent and temporary seismic stations was inverted for seismic imaging. The seismic images suggest a thick high-velocity lithospheric mantle beneath the Zagros. Other high-velocity domain observed below ∼300 km depth, which is not connected to the Zagros lithosphere, is interpreted as a slab segment in the upper mantle beneath the collision zone. The low-velocity anomalies beneath Central Iran and Alborz, consistent with a weak upper mantle, may be a result of upwelling asthenosphere and partial melting of the subducted lithosphere.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We found that 〈span〉SH〈/span〉〈sub〉diff〈/sub〉 phases generated by earthquakes in the Fiji–Tonga, recorded in India, are accompanied by secondary pulses. We interpreted them as a consequence of multipathing of 〈span〉S〈/span〉 waves caused by the Pacific large low-shear-velocity province (LLSVP). We analysed the differential traveltimes between 〈span〉SH〈/span〉〈sub〉diff〈/sub〉 and the secondary pulse, together with the absolute 〈span〉SH〈/span〉〈sub〉diff〈/sub〉 arrival times, to constrain the thickness and velocity perturbations in the western end of the Pacific LLSVP. Our preferred model shows a lateral variation in the thickness of the LLSVP; the southern part reveals a thicker (300 km) low-velocity region compared to the northern part (200 km). However, the velocity perturbations of the LLSVP appear to be comparable ($-1.5\hbox{ per cent}$). The results are consistent with a scenario that the LLSVP is a chemically distinct pile with significant surface topography.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Non-uniqueness in the geophysical inverse problem is well recognized and so too is the ability to obtain solutions with different character by altering the form of the regularization function. Of particular note is the use of ℓ〈sub〉〈span〉p〈/span〉〈/sub〉 norms with 〈span〉p〈/span〉 ∈ [0, 2] which gives rise to sparse or smooth models. Most algorithms are designed to implement a single ℓ〈sub〉〈span〉p〈/span〉〈/sub〉 norm for the entire model domain. This is not adequate when the fundamental character of the model changes throughout the volume of interest. In such cases we require a generalized regularization function where each sub-volume of the model domain has penalties on smallness and roughness and its own suite of ℓ〈sub〉〈span〉p〈/span〉〈/sub〉 parameters.Solving the inverse problem using mixed ℓ〈sub〉〈span〉p〈/span〉〈/sub〉 norms in the regularization (especially for 〈span〉p〈/span〉 〈 1) is computationally challenging. We use the Lawson formulation for the ℓ〈sub〉〈span〉p〈/span〉〈/sub〉 norm and solve the optimization problem with Iterative Reweighted Least Squares. The algorithm has two stages; we first solve the 〈span〉l〈/span〉〈sub〉2〈/sub〉-norm problem and then we switch to the desired suite of ℓ〈sub〉〈span〉p〈/span〉〈/sub〉 norms; there is one value of 〈span〉p〈/span〉 for each term in the objective function. To handle the large changes in numerical values of the regularization function when 〈span〉p〈/span〉 values are changed, and to ensure that each component of the regularization is contributing to the final solution, we successively rescale the gradients in our Gauss–Newton solution. An indicator function allows us to evaluate our success in finding a solution in which components of the objective function have been equally influential.We use our algorithm to generate an ensemble of solutions with mixed ℓ〈sub〉〈span〉p〈/span〉〈/sub〉 norms. This illuminates some of the non-uniqueness in the inverse problem and helps prevent overinterpretation that can occur by having only one solution. In addition, we use this ensemble to estimate the suite of 〈span〉p〈/span〉 values that can be used in a final inversion. First, the most common features of our ensemble are extracted using principal component analysis and edge detection procedures; this provides a reference model. A correlation of each member of the ensemble with the reference model, carried out in a windowed domain, then yields a set of 〈span〉p〈/span〉 values for each model cell. The efficacy of our technique is illustrated on a synthetic 2-D cross-well example. We then apply our technique to the field example that motivated this research, the 3-D inversion of magnetic data at a kimberlite site in Canada. Since the final regularization terms have different sets of 〈span〉p〈/span〉 values in different regions of model space we are able to recover compact regions associated with the kimberlite intrusions, continuous linear features with sharp edges that are associated with dykes and a background that is relatively smooth. The result has a geologic character that would not have been achievable without the use of spatially variable mixed norms.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The ice cap covering Antarctica has long limited our understanding of the continental-scale crustal model due to its inaccessibility and the resulting logistical difficulties when executing geophysical field work, such as seismograph deployment. Resolving a high spatial resolution crustal model for Antarctica where seismographs are sparsely distributed stimulates scientific interest in this relatively less studied continent. In this study, we utilize satellite gravity observations from the global gravity model EIGEN-6C4 to create an alternative crustal thickness model of Antarctica. The gravity data were corrected for sediments, topography and ice cover. Furthermore, we considered the gravity effect due to vertical deformation of the lithosphere caused by ice load besides the earth's curvature in the modelling. We inverted the corrected gravity data using the regularized Bott's inversion method in spherical approximation and constrained the results by seismic observations. This crustal thickness model shows a thicker average crust in East Antarctica and a thinner one in West Antarctica. The thickest crust is in the Gamburtsev Subglacial Mountains with a Moho depth of over 40 km. The thicker crust is particularly evident along the Transantarctic Mountains and the Dronning Maud lands. Comparisons with existing models show a good correlation in gravity-constrained areas. Differences appear in the sedimentary basins and crust with thickness closer to seismic point observations. Overall, our crustal model is relatively improved than the existing gravity derived models.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The formation of fold-thrust belts at convergent margins is a dynamic process. Accretion of weak sediments to the front of the overriding plate results in crustal thickening and continued flexural subsidence of the underthrusting plate. Fold-thrust belts are often treated as a Coulomb wedge having self-similar geometries with a critical taper, and either a rigid or isostatically compensated base. In this paper, we build upon this work by developing a new dynamic model to investigate both the role of the thickness and material properties of the incoming sediment, and the flexure in the underthrusting plate in controlling the behaviour and evolution of fold-thrust belts. Our analysis shows that the evolution of fold-thrust belts can be dominated by either gravitational spreading or vertical thickening, depending on the relative importance of sediment flux, material properties and flexure. We apply our model to the Makran accretionary prism and the Indo–Burman Ranges, and show that for the Makran flexure must be considered in order to explain the dip of the sediment–basement interface from seismic reflection profiles. In the Indo–Burman Ranges, we show that incoming sediment thickness has a first-order control on the variations in the characteristics of the topography from north to south of the Shillong Plateau.〈/span〉