Systems and methods for in-situ characterization of permafrost sites

Aspects of this disclosure relate to the characterization of permafrost.

The adverse effects of climate change and global warming on the built environment are indisputable, undeniable and accelerating at an alarming pace. According to Canada's Changing Climate Report (2019), in the Arctic regions, temperatures have been warming at approximately twice the rate of the rest of the world. This drastic trend in climate warming will no doubt affect permafrost temperatures and conditions, continued growth in greenhouse gas emissions, and further adding to the high cost of northern development. Planning and designing climate-resilient northern infrastructure as well as projecting associated reductions in permafrost thickness from climate model simulations requires an estimation of permafrost properties. Important properties include the current thickness of active layer, porosity, saturation degree of water, saturation degree of ice, soil type, as well as mechanical properties such as shear modulus and bulk modulus.

Permafrost is defined as the ground that remains at or below 0° C. for at least two consecutive years. The upper layer of the ground in permafrost areas, termed as the active layer, seasonally thaws and freezes. The thickness of the active layer greatly depends on the local geological and climate conditions such as vegetation, soil composition, air temperature, solar radiation, and wind speed.

Within the permafrost, the distribution of the ground ice is highly variable. Also, ground ice can be present under distinctive forms including (1) pore ice, (2) segregated ice, and (3) ice-wedge [7, 22]. Pore water, which fills or partially fills the pore space of the soil, freezes in place if the temperature drops below the freezing point [27]. On the other hand, segregated ice is formed when water migrates to the freezing front and it can cause excessive deformations in frost-susceptible soils. Frost-susceptible soils, e.g. silty or silty clay soils, have relatively high capillary potential and moderate relative and intrinsic permeability. The ice wedge is referred to large masses of ice formed in thermal contraction cracks. During the winter months, the ice expands as the ground freezes to form cracks in the subsurface [18]. Ice wedges are large masses of ice formed over many centuries by repeated frost cracking and ice vein growth.

Design and construction of engineering works on permafrost normally follow one of two broad principles which are based on whether the frozen foundation soil is thaw-stable or thaw-unstable (ice-rich permafrost). This distinction is determined by the amount of ice content within the permafrost. Ice-rich permafrost contains ice in excess of the water content at saturation. The construction on thaw-unstable permafrost is challenging and requires remedial measures since upon thawing, permafrost will experience significant thaw-settlement and suffer loss of strength to a value significantly lower than that for similar material in an unfrozen condition. Consequently, remedial measures for excessive soil settlements or new design of infrastructure in permafrost zones affected by climate warming require a reasonable estimation of the amount of ice content within the permafrost (frozen soil). The rate of settlement relies also on the mechanical properties of the foundation permafrost at the construction site. Furthermore, a warming climate can accelerate the microbial breakdown of organic carbon stored in permafrost and can increase the release of greenhouse gas emissions, which in return would accelerate climate change [30].

Several in-situ techniques have been employed to characterize or monitor permafrost conditions. For example, techniques such as remote sensing [34, 1, 37], and ground penetrating radar (GPR) [25, 6, 33] have been used to detect ice-wedge formations within the permafrost layers. Also, electrical resistivity tomography (ERT) has been extensively used to qualitatively detect pore-ice or segregated ice in permafrost based on the correlation between the electrical conductivity and the physical properties of permafrost (e.g., unfrozen water content and ice content) [9, 10, 29, 35]. The apparent resistivity measurement by ERT is higher in areas having high ice contents [35]; however, at high resistivity gradients, the inversion results become less reliable, especially for the investigation of permafrost base [12, 23]. Furthermore, in ERT investigations, the differentiation between ice and certain geomaterials can be highly uncertain due to their similar electrical resistivity properties [15]. GPR has been also used for mapping the thickness of the active layer; however, its application is limited to a shallow penetration depth in conductive layers due to the signal attenuation and high electromagnetic noise in ice and water [15]. It is worth mentioning that none of the above-mentioned methods characterizes the mechanical properties of permafrost layers.

Non-destructive seismic testing, including multi-channel analysis surface waves (MASW) test [8, 9], passive seismic test with ambient seismic noise [14, 26], seismic reflection [2], and seismic refraction method have been previously employed to map the permafrost layer based on the measurement of the shear wave velocity. In the current seismic testing practice, it is commonly considered that the permafrost layer (frozen soil) is associated with a higher shear wave velocity due to the presence of ice in comparison to unfrozen ground. However, the porosity and soil type can also significantly affect the shear wave velocity [20]. In other words, a relatively higher shear wave velocity could be associated to an unfrozen soil layer with a relatively lower porosity or stiffer solid skeletal frame, and not necessarily related to the presence of a frozen soil layer. Therefore, the detection of permafrost layer and permafrost base from only the shear wave velocity may lead to inaccurate and even misleading interpretations.

In one example, a method of quantifying a plurality of parameters of a subsurface structure comprising a region of interest, the method comprising the steps of:

In another example, a system for quantifying a plurality of parameters of a subsurface structure comprising a region of interest, the system comprising:

In another example, a computer-readable medium storing instructions which, when executed by a processor, cause a computing device to perform operations comprising:

Advantageously, an in-situ seismic surface wave technique estimates the physical and mechanical properties of permafrost sites, based on the existence of two Rayleigh waves, that is, a first Rayleigh wave (R1) signal which is dependent on at least one of the plurality of parameters of the region of interest, and a second Rayleigh wave (R2) signal. Generally, the R2 wave is highly sensitive to the physical properties (e.g., unfrozen water content, ice content, and porosity) of permafrost while it is relatively insensitive to the mechanical properties (e.g., shear modulus and bulk modulus). On the other hand, the R1 wave velocity depends strongly on the mechanical properties of foundation permafrost. The in-situ surface wave measurements revealed the experimental dispersion relations of both types of Rayleigh waves from which relevant properties of a permafrost site can be derived by means of a hybrid inverse and multi-phase poromechanical approach.

The identification of two distinctive types of Rayleigh waves in the seismic measurement at permafrost sites is useful in characterizing permafrost. The identification of the R2 wave also allows the quantitative characterization of permafrost properties (e.g., unfrozen water content, ice content, and porosity) independently without making any assumptions on the mechanical properties of the permafrost. This significantly reduces the difficulties in the inversion of the multi-layered three-phase poromechanical model since the dependent optimization variables are largely reduced. The inversion results from the R2 wave dispersion relation can be further used in the characterization of the mechanical properties of permafrost based on the R1 wave. This also increases the stability and convergence rate of the inversion model.

Accordingly, the physical properties of the permafrost, such as the amount of unfrozen water content, ice content, and porosity as well as the mechanical properties such as the shear modulus and bulk modulus of permafrost layers, may be quantified. In addition, the location of the permafrost table and permafrost base may also be determined.

shows a system for determining characteristics of a subsurface structure, general designated by numeral, in one implementation. Subsurface structurecomprises active layerwith surface, permafrost layerand unfrozen ground. Seismic sourcemay be selectively and periodically actuated to impart energy in the form of an acoustic or pressure wave through subsurface structureand to generate seismic waves within subsurface structure. In one example, seismic sourcemay be ambient noise or an object impacting surface, or other means. Various types of controlled seismic sources may be used. Seismic detectoris preferably secured in essentially fixed relation to subsurface structure. In one example, seismic detectoris a geophone which measures the vertical particle velocity induced by seismic waves in subsurface structure. Seismic sourcemay be in electrical communication with computing deviceto control the actuation of seismic source, and seismic detectoris in electrical communication with computing devicefor transmission of detected seismic signals thereto.

Seismic detectorsmay be arranged in a line, or in an array, or positioned in specific distances (radius) from a center (with an equally spaced radius Δr, but an arbitrary angle θ) for performing a seismic survey operation with respect to a subsurface structure. The subsurface structuremay have at least one subsurface region of interest, such as permafrost layer.

Computing devicecomprises processor, memoryfor storing data associated with the seismic signals received from seismic detectors, and for storing computer readable instructions executable by processorto at least, via input/output moduleand receive the detected seismic signals from seismic detectors, and determine the characteristics of a region of interest of subsurface structure. Accordingly, through analysis of these detected seismic signals, the characteristics of permafrost layer, such as the amount of ice content, unfrozen water content, and porosity, as well as the shear modulus and bulk modulus of permafrost layers, including the location of the permafrost table and permafrost base, may be determined.

shows flow diagramdepicting a computer-implemented method for determining the characteristics of a region of interest of subsurface structureby processing seismic data in accordance with various implementations described herein, andshows a general schematic of a method for determining the characteristics of a region of interest of subsurface structureby processing seismic data from a multi-channel analysis surface waves (MASW) testing method. It should be understood that while the operational flow diagram indicates a particular order of execution of the operations, in other implementations, the operations might be executed in a different order. Further, in some implementations, additional operations or blocks may be added to the method. Likewise, some operations or blocks may be omitted.

Looking at, at block, seismic sourceis actuated source to selectively and/or periodically impart energy in the form of an acoustic or pressure wave through subsurface structure. For example, the imparted energy may comprise a known amplitude, frequency and duration.

At block, geophoneacquires seismic signals from subsurface structure, including signals for a region of interest (i.e., “the received seismic data”). The region of interest may include an area of the subsurface in the earth that may be of particular interest, such as permafrost layer. Surface wave measurements are then derived from the acquired seismic signals.

At block, computing devicegenerates an initial model to obtain a dispersion relation for the R1 and R2 wave modes.

At block, computing devicereceives the acquired signals from geophones, uses the surface wave measurements and computes a dispersion relation of R1 and R2 wave modes. The dispersion of R2 waves is then used to characterize the physical properties (e.g., unfrozen water content, ice content, and porosity) of permafrost layerof subsurface structure.

At block, an initial estimate of the physical properties of active layer, permafrost layerand unfrozen groundis postulated.

At block, a forward three-phase poromechnical dispersion prediction model is used to compute the theoretical dispersion relation of the R2 mode. The dispersion relation prediction model offers theoretical dispersion relation functions for propagation of stress waves.

At block, samples within a parameter space of the poromechanical dispersion prediction model are ranked based on the L2 norm between the experimental and theoretical dispersion relations between the R1 wave and the R2 wave.

At block, a Neighborhood sampling for the reduction of L2 norm is performed.

At block, the best samples with the minimum loss function are selected to obtain the most likely physical properties of active layer, permafrost layer, and unfrozen ground.

At block, the steps-for dispersion inversion for the mechanical properties based on the R1 dispersion relation are repeated.

At block, the best samples are selected and the mechanical properties are obtained. The mechanical properties can be derived based on the dispersion relation of the R1 wave mode in a similar manner.

Looking back at step, in one implementation surface wave measurements are employed to obtain a dispersion relation of R1 and R2 wave modes from the experimental measurement. A frozen soil from subsurface structuremay be composed of three phases: solid skeletal frame, pore-water, and pore-ice, and the Green-Lagrange strain tensor (ϵ) for infinitesimal deformations expressed as displacement vector

for solid skeleton, pore water and pore ice may be expressed as:

The constitutive models may be defined as the relation between the stress and strain tensors for solid skeleton, pore water and pore ice may be expressed as:

(m, ranging from 1 to 3, represents the different phases) are defined as follows:

The momentum conservation considers the acceleration of each component and the existing relative motion of the pore ice and pore water phases with respect to the solid skeleton. The momentum conservation for the three phases is given by Equation D.3.

Through the infinitesimal kinematic assumption (Equation D.1), the stress-strain constitutive model [3] (Equation D.2), and the conversation of momentum (Equation D.3), the field equation can be written in the matrix form

The matrix,,andare given in Appendix E. The field equations can also be written in the frequency domain by performing convolution with eiωt. The field equations in the Laplace domain are obtained by replacing ω with i·s (i=−1 and s the Laplace variable).

To obtain the analytical solution, the Helmholtz decomposition is used to decouple the P waves (P1, P2, and P3) and S waves (S1 and S2). The displacement vector (ū) is composed of the P wave scalar potentials q and S wave vector potentials=(ψ, ψ, ψ). Since P waves exist in the solid skeleton, pore-ice and pore-water phases, three P wave potentials are used, including φ, φand φmay be expressed as follows:

The detailed steps for obtaining the analytical solutions for P waves and S waves using the Eigen decomposition are summarized in Appendix D. After obtaining the stiffness matrix for each layer, the global stiffness matrix can be obtained by applying the continuity conditions between layer interfaces. The stiffness assembling method is shown in. The global stiffness is denoted as H matrix for simplicity.

The dispersion relation is obtained by setting a zero stress condition at the surface (z=0). To obtain the non-trivial solution, the determinant of the global stiffness matrix has to be zero, as expressed in Equation 1 [38].

The global stiffness matrix, H(ω, k), is a function of angular frequency ω and wavenumber k. For a constant frequency, the value of the wavenumber can be determined when the determinant of the global stiffness matrix is zero. The dispersion curve is also commonly displayed as frequency versus phase velocity,

The different wavenumbers determined at a given frequency correspond to dispersion curves of different modes. To extract the fundamental mode of the R1 wave, the velocities of P1 wave and S1 wave are calculated first for the given physical properties and mechanical properties of each layer. The global stiffness matrix for the R1 wave can be decomposed into the components related only to the P1 and S1 wave velocities. This is viable since we have proved that the R1 wave is generated by the interaction between the P1 and S1 waves. This approach avoids the difficulties in differentiating the higher modes of R2 wave from the fundamental mode of the R1 wave. The detailed root search method has been documented in [18].

Next, inversion is performed. The aim function is defined as the Euclidean norm between the experimental and numerical dispersion relations. The problem is formulated in Equation 2:

Here, we used the neighborhood algorithm that benefits from the Voronoi cells to search the high-dimensional parameter space and reduce overall cost function [28]. The algorithm contains only two tuning parameters. The neighborhood sampling algorithm includes the following steps: a random sample is initially generated to ensure the soil parameters are not affected by the local minima. Based on the ranking of each sample, the Voronoi polygons are used to generate better samples with a smaller objective function. The optimization parameters are scaled between 0 and 1 to properly evaluate the Voronoi polygon limit, which is calculated as shown in paragraph below. The detailed description of the neighborhood algorithm is described by [28].