To better understand the relationship between P-wave velocities and ice content in saturated, unconsolidated saline permafrost, we constructed an effective-medium model based upon ultrasonic P-wave data that were obtained from earlier laboratory studies. The model uses a two-end-member mixing approach in which an ice-filled, fully frozen end member and a water-filled, fully unfrozen end member are mixed together to form the effective medium of partially frozen sediments. This mixing approach has two key advantages: (1) It does not require parameter tuning of the mixing ratios, and (2) it inherently assumes mixed pore-scale distributions of ice that consist of frame-strengthening (i.e., cementing and/or load-bearing) ice and pore-filling ice. The model-predicted P-wave velocities agree well with our laboratory data, demonstrating the effectiveness of the model for quantitatively inferring ice content from P-wave velocities. The modeling workflow is simple and is largely free of calibration parameters — attributes that ease its application in interpreting field data sets.

Saline permafrost is widespread in subsea and coastal areas of the Arctic and Antarctic (Osterkamp, 1989; Hivon and Sego, 1993; Brouchkov, 2003; Ingeman-Nielsen et al., 2008). It typically contains soluble salts originating from seawater of the present day or geologic past. Although permafrost of all types is expected to contain some salts, saline permafrost is only recognized when its salinity is sufficiently high. That is, it should contain soluble salts that take up at least 0.05 wt% of its dry weight (Brouchkov, 2002). For fully saturated permafrost, this roughly translates to a pore-water salinity of approximately 1–2 ppt (0.1–0.2 wt%).

The most striking characteristic that distinguishes saline permafrost from its nonsaline counterpart is its heightened temperature sensitivity. Unlike nonsaline permafrost that is generally stable unless its temperature approaches 0°C, saline permafrost is highly responsive to warming even at temperatures that are well below 0°C (Ruffell et al., 1990). This is attributed to two effects of the dissolved salts: (1) freezing-point depression and (2) progressive salination of the residual pore water. In the latter process, the residual pore water become more saline because salts are excluded from ice (known as the freeze desalination process that results from the inability of ice crystals to accommodate salts). Such salination further lowers the freezing point of the residual water and consequently slows down freezing. As a result, saline permafrost usually takes the form of a partially frozen “slush,” in which ice and brine coexist in a delicate equilibrium that can be easily perturbed.

Being widely distributed, saline permafrost has profound influences on the future climate and on engineering practices of the present day. Given its temperature sensitivity, saline permafrost is likely to thaw at an accelerated pace in a warming climate. Being a partially frozen slush, saline permafrost deforms easily under external loads or self-weight, which makes it an undesirable candidate as foundation ground and a frequent culprit behind landslides and infrastructure damage (Ogata et al., 1983; Hivon and Sego, 1993; Brouchkov, 2003). For example, Nixon (1988) concludes that the foundation-bearing capacities of permafrost can be reduced by factors of two to three if the pore-water salinity exceeds 10–20 ppt. Characterized by these unusual thermal and mechanical properties, saline permafrost is becoming increasingly important to the changing Arctic landscape (Portnov et al., 2013, 2014). To predict its likely responses in a warming climate, as well as to prevent and mitigate the associated hazards, we need to be able to detect and characterize saline permafrost.

Geophysical methods are effective at mapping and characterizing permafrost primarily because of the marked contrasts between the geophysical properties of ice and water (e.g., seismic refraction and surface-wave methods, electrical resistivity imaging, electromagnetic induction, and ground-penetrating radar; see reviews by Hauck and Kneisel [2008], Kneisel et al. [2008], and Hauck [2013]). Among the commonly used geophysical parameters for permafrost investigations, seismic velocities are the most relevant to engineering practices because they are often indicative of the mechanical strengths of geomaterials (Schön, 2011). Prior permafrost studies have shown that increases in ice content generally coincide with increases in mechanical strengths and seismic velocities (Timur, 1968; Nakano and Froula, 1973; Zimmerman and King, 1986). Containing less ice, saline permafrost can be identified as low-velocity anomalies that stand out against a background of nonsaline permafrost (Collett and Bird, 1988, 1993; Schmitt et al., 2005; Hubbard et al., 2013; Dou and Ajo-Franklin, 2014). Such a contrast enables seismic mapping of saline permafrost units on the scales of geotechnical interests (i.e., tens to hundreds of meters).

Besides its spatial distributions, the ice content of saline permafrost is a key parameter that we hope to infer from seismic measurements. This is because ice content plays a central role in assessing present-day strength of saline permafrost as well as its likely risk of thaw settlement in the future (Nelson et al., 2001, 2002). However, quantitatively relating seismic velocities to ice content has proven challenging because the mapping between the two depends heavily on microstructural distributions of ice that are rarely known a priori. For example, if ice acts like cement that bonds sediment grains together, even a small amount of ice is enough to cause marked increases in seismic velocities. In contrast, ice could be suspended in the pore space without contacting sediment grains (i.e., the pore-filling model); in this case, ice becomes part of the pore fluid and influences seismic velocities by changing density and compressibility of the pore fluids (ice is less dense and less compressible than water). Because seismic waves transmit energy more efficiently through the solid frame than through the pore fluids, the same amount of pore-filling ice will only yield moderate increases in seismic velocities. Examples of the cementing and pore-filling models are well-represented in the gas hydrate literature (Ecker et al., 1998; Helgerud, 2001; Waite et al., 2004).

Fortunately, it is possible to infer ice distributions through rock-physics-based models, in which an assumed ice distribution can be verified by comparing the associated model predictions against measurements. Measurements used for this purpose must be of high quality, and enough prior knowledge must already exist so that no other unknowns can bias the inferred ice distributions. One suitable measurement approach is time-lapse seismic monitoring of permafrost that undergoes controlled cooling and/or warming. Although the time-lapse approach is amenable to field and laboratory experiments, development of rock-physics data sets is preferably done at the laboratory scale in which tight control of temperature and fluid composition is attainable.

In Dou et al. (2016), we report a laboratory experiment in which we acquire ultrasonic P-wave measurements from saturated, unconsolidated saline permafrost samples during freeze/thaw cycles. In this paper, we present a set of rock-physics modeling efforts that aim at quantitatively relating the P-wave velocity (VP) measurements of Dou et al. (2016) to the corresponding ice content. We construct our model using a two-end-member mixing approach that improved upon the original work of Minshull et al. (1994). The workflow of the modeling is simple and is largely free of calibration parameters, and the model predictions agree well with our laboratory measurements.

The direct measurements of Dou et al. (2016) are sets of (VP, T) points (i.e., P-wave velocities VP and the associated temperatures T). To analyze these data sets using established rock-physics relationships between VP and ice saturation si (si=Volice/Volpore: volume fraction of the total pore space occupied by ice), we need to convert the temperature measurements into the corresponding ice saturations (hereinafter referred to as T-to-si conversion).

For controlled freezing and thawing of saline permafrost in which NaCl is the predominant salt, the T-to-si conversion is straightforward for two reasons: First, NaCl-H2O brine has phase-diagram expressions that are well-established and readily usable for the T-to-si conversion. Second, because the water-to-ice phase transitions are progressive in the presence of dissolved salts, the associated changes in VP occur over a broad temperature range. This allows the VP versus T mapping of the partially frozen sediments to be measured accurately within temperature resolutions that are achievable in typical laboratory settings (approximately 0.1°C in our case). In this section, we present the underlying principles, assumptions, and procedures for the T-to-si conversion.

Brief overview of the NaCl-H2O phase diagram

We first review the nomenclatures and expressions of the NaCl-H2O phase diagram that are essential to this study, namely, the liquidus, the eutectic point, and the expressions for determining the freezing point and the equilibrium salinity.

As shown in Figure 1, the liquidus on the NaCl-H2O phase diagram is the temperature-versus-salinity (T-Sn) curve above which the solution is entirely liquid. The most basic information it conveys is the positive correlation between the initial salinity of the solution (Sn0) and its freezing point (Tfp) lying on the liquidus (we denote it as “the Tfp versus Sn0 interpretation”): the higher the initial salinity, the lower the freezing point. In this sense, the liquidus quantitatively depicts the freezing-point depression effect of the dissolved salts. Note that the lowering of the freezing point is limited by a critical temperature known as the eutectic point (Teutectic; approximately 21°C for NaCl-H2O brine). Once the temperature drops below Teutectic, all dissolved salts, regardless of the initial salinity of the solution, rapidly precipitate as solid crystals, leaving behind desalinated residual water that quickly freezes at subeutectic temperatures.

An alternative to the Tfp versus Sn0 interpretation of the liquidus is the equilibrium salinity interpretation. At a given temperature within the partially frozen regime (T[pf];Teutectic<T[pf]<Tfp; i.e., between the freezing point and the eutectic point), the residual solution maintains an equilibrium salinity (SnEQ) that equals the minimum salinity needed for the solution to remain liquid at this temperature. Under this interpretation, we view the temperature axis as the partially frozen temperatures, and the corresponding salinities on the liquidus as the equilibrium salinities.

In short, the liquidus has twofold utility (the expressions shown below come from the polynomial regressions of experiment data by Potter et al., 1978):

  • Determining freezing point Tfp (in °C) for a given initial salinity Sn0 (in wt%; Sn0=Sn|TT[af], where T[af] denotes the above-freezing temperature T[af]>Tfp):
  • Determining equilibrium salinity SnEQ (in wt%; SnEQ=Sn|TT[pf]) at a given temperature in the partially frozen regime (in °C; TT[pf], where Teutectic<T[pf]<Tfp):

Note that although NaCl brines of different initial salinities have different freezing points, once they become partially frozen, the residual brines will have the same equilibrium salinity at any given partially frozen temperature. If we only measure the salinities of the residual brines, it would look as if the initial salinities had been “forgotten” upon freezing. However, as we will soon illustrate, the initial salinity is “remembered” by the ice saturation (lower initial salinity yields higher ice saturation, and vice versa). This leads to a simple way of experimentally sampling various levels of ice saturations by subjecting samples of different initial salinities to the same freeze/thaw processes. For our laboratory measurements (as summarized in the upcoming section), this is a key factor that has ensured good measurement repeatability and dense sampling at various levels of ice saturations.

Brief overview of the laboratory measurements

As mentioned earlier, our rock-physics modeling is based upon ultrasonic P-wave measurements of Dou et al. (2016). Here, we briefly summarize this earlier experimental study.

Experiment objectives and materials

To better understand how seismic properties of saline permafrost vary with respect to temperatures and pore-water salinities, Dou et al. (2016) conduct ultrasonic P-wave measurements for saturated, unconsolidated sediment samples over a broad range of salinities (0.0–2.5 M [0.0–130 ppt]) and temperatures (10°C to 30°C). Two types of samples were used in the experiment: six coarse-grained Ottawa sand samples saturated with distilled water or brine of different salinities and a fine-grained saline permafrost core sample at in situ salinity (extracted from the Barrow Peninsula on the Alaskan Arctic Coastal Plain) (Table 1 shows the texture and compositions of the samples). For the coarse sand samples, the measurements are mainly affected by pore-water salinities. For the fine-grained core sample, the measurements are affected by salinities and surface effects (i.e., unfrozen water retention that is caused by capillary effects of narrow pores and adsorption effects of charged mineral surfaces).

Key experiment observations

We condense the key results from Dou et al. (2016) into Figure 2. The major observations are summarized as follows:
  • The influences of water-to-ice phase transitions: Regardless of salinity and sediment texture, the onset of freezing is signified by rapid increases of P-wave velocities and sharp declines of P-wave amplitudes.

  • The influences of pore-water salinity: For the nonsaline sample, the freezing-induced velocity increases and amplitude decreases starts near 0°C. Then, as the temperatures decrease further, the velocities and amplitudes quickly plateau at high values. For the saline samples, the higher the initial pore-water salinity, the lower the freezing point. Neither velocities nor amplitudes show signs of plateauing until the temperature drops below the eutectic point. Between the freezing point and the eutectic point, P-wave velocities of the saline samples are markedly lower than the nonsaline sample, manifesting reductions in stiffness because less ice can form in the presence of dissolved salts. Interestingly, the amplitudes remain low over the course of the partially frozen stage, whereas the velocities increase monotonically with decreasing temperatures, contradicting the commonly expected positive correlation between amplitudes and velocities.

  • The influence of fine-grained particles: For the fine-grained core sample, surface effects, though expected to cause freezing-point depression in nonsaline samples (Anderson et al., 1973), do not appear to deepen the extent of freezing-point depression that is expected from salinity alone. It is only until the temperature has fallen below the eutectic point when the influences of surface effects become noticeable, and eventually result in an apparent shift of the eutectic point from the expected 21°C down to 31.8°C (Figure 2). We interpret this shift not as a true decrease of the eutectic point, but as a manifestation of the surface-effect-induced freeze inhibition exerted on the largely desalinated residual water. In other words, due to surface effects, the desalinated water does not immediately freeze near the eutectic point (unlike in the coarse sand samples); instead, they remain liquid until the temperature has become so low that surface effects can no longer inhibit freezing.

Figure 2.

Key experimental results presented in Dou et al. (2016). (a) The P-wave velocities and (b) relative peak-to-peak amplitudes (relative to amplitude values at the highest above-freezing temperature Tmax; ΔAmplitude(%)=(A|TTmaxA|T=Tmax)/(A|T=Tmax)×100%, where A denotes the peak-to-peak amplitude). Ppt, parts per thousand (=10×wt%); M, molar (1  M=1  mol/L).

Figure 2.

Key experimental results presented in Dou et al. (2016). (a) The P-wave velocities and (b) relative peak-to-peak amplitudes (relative to amplitude values at the highest above-freezing temperature Tmax; ΔAmplitude(%)=(A|TTmaxA|T=Tmax)/(A|T=Tmax)×100%, where A denotes the peak-to-peak amplitude). Ppt, parts per thousand (=10×wt%); M, molar (1  M=1  mol/L).

Observation-based assumptions and simplifications for modeling

At above-eutectic temperatures, dissolved salts are considered to be the sole factor that dictates the freezing point and the subfreezing variations in water/ice content, whereas surface effects are assumed to be negligible. Because the surface effects mainly affect the freezing process of bound water (water that is held in pores or adsorbed on particle surface and does not move under the pull of gravity), and bound water only starts to freeze when most of the free water (water that moves through the sediments under the pull of gravity) has frozen, this assumption implies a persistent presence of free water that is only possible if the pore-water salinity of the sediments is sufficiently high (pinpointing the lower-limit salinity is beyond the scope of this study).

At the eutectic point, all of the salts are assumed to precipitate out of the solution. In other words, although freezing-point depression and eutectic-point depression can occur in saline permafrost (Toner et al., 2014), we do not allow eutectic-point depression in our modeling, and thus all the residual pore water at the subeutectic temperature is assumed to be nonsaline. This simplification is necessary for three reasons: (1) In our experiment, salt precipitation was already visible to the naked eye when the temperature of the sample reached the eutectic point, (2) unlike the deterministic processes of freezing-point depression caused by salts, eutectic-point depression is the result of subeutectic supercooling, a stochastic process for which universal equations are unavailable, and (3) given that the extreme conditions of near-eutectic and subeutectic temperatures are rarely encountered in terrestrial permafrost systems, simplifications can be applied without degrading the usefulness of the model.

One key outcome of these assumptions is a clear-cut division between the salinity-dominated regime and the surface-effects-dominated regime (Figure 3). This division simplifies the procedure required for inferring ice saturations from the temperature measurements, the initial pore-water salinity, and the sediment texture information. For the salinity-dominated regime, we only need to use phase-diagram expressions that relate ice saturations to temperatures and initial pore-water salinities. For the surface-effects-dominated regime, we only use empirical relationships that relate ice saturations of nonsaline permafrost to temperatures and the particle sizes of the sediments. In the next section, we will present the relevant equations in detail.

In this section, we first present the parameters and assumptions involved in estimating ice saturations (the volume ratios of the ice and the total pore space). We then derive the governing equations that are specific to the salinity-dominated and the surface-effects-dominated regimes, respectively. Note that regardless of the specific regime, the porosity of the sediment frame is one of the parameters required for estimating ice saturations.

Parameter requirements

For the salinity-dominated regime, phase-diagram expressions serve as the foundation for estimating ice saturations, for which the initial salinity of the pore water must be known a priori. In addition, because the direct output of the phase diagram expressions is the mass ratio between the residual unfrozen brine and the initial brine content (prior to freezing), brine densities are required for converting these mass ratios into volume ratios.

For the surface-effects-dominated regime, because no expressions are universally applicable, we use the empirical expressions of Anderson et al. (1973) to relate a specific area of the sediment grains to the mass ratio of the unfrozen water relative to the total mass of the sediments. An underlying requirement is that the grain sizes of the sediments should be sufficiently uniform such that a representative value of the specific surface area can be found. The conversion from mass ratio to volume ratio requires density values for the sediment grains and the nonsaline water.

Assumption of constant pore volume

In the equation for obtaining ice saturations, si=Volice/Volpore, we assume that the total volume of the pore space Volpore remains constant, despite the fact that approximately 9% of the volume expansion is expected as water turns into ice. The rationales behind the constant pore volume assumption are as follows:

  • Unconfined sample: Our data were acquired from samples that were not confined at the top. Although ice takes up more space than water and thus could cause a pore-space deficit that forces some portions of the water to leave its original position, the sediments appear to have adequate vertical conductivity that the extra water could migrate toward the unconfined top without having to “break” grain contacts apart to create extra pore space.

  • Avoid introducing extra fitting parameters: In order for ice to push the grains apart, it needs to be in contact with the surrounding grains. However, some fractions of the ice (particularly at the early stage of freezing) may not be in contact with the grains and thus are unlikely to change the pore volume. Although one could specify a certain amount of ice to be in contact with the grains, this quantity itself is not known a priori and would become a fitting parameter.

Estimating ice saturations

Equations for the salinity-dominated regime

For the salinity-dominated regime, the governing equations are derived from the NaCl-H2O phase-diagram expressions described earlier. Key assumptions in this step are (1) phase diagrams of unconfined brines are directly applicable to brines that are confined in porous media and (2) during the freeze desalination process, all the dissolved salts are always fully rejected into the residual pore water. Prior to estimating the ice saturations, we first need to determine the freezing point (Tfp) that is specific to a given initial salinity (Sn0) using equation 1.

Within the salinity-dominated regime, we use equation 2 to determine the salinity of the residual brine. Because we assume that all the dissolved salts are always fully rejected into the residual pore water, the mass conservation of salts yields the following relationship:
where msalt is the total mass of the dissolved salts, mw0 is the initial mass of the pore water prior to freezing, mw is the mass of the residual pore water at partially frozen temperatures, and Sn is the equilibrium salinity that is determined by equation 2.
Finally, by replacing mass m with the product of density ρ and volume Vol (m=ρVol) in equation 3, we obtain
where ρw0 and ρw denote the densities of above-freezing brine and subfreezing residual brine. Considering that the initial volume of the pore water equals the total pore volume of the saturated sediment (Volpore=Volw0), we rearrange equation 4 and obtain the saturation of the residual brine as
Because the saturation values of water and ice are related by sw+si=1, we obtain the equation for determining ice saturation after replacing Sn with equation 2:

In equation 6, brine densities ρw0 and ρw are obtained by using equation 27a and 27b in Batzle and Wang (1992). The constant pore volume assumption described at the beginning of this section is implicitly used to arrive at equation 6.

Equations for the surface-effects-dominated regime

For the surface-effects-dominated regime that is associated with subeutectic temperatures, the residual pore water is assumed to be nonsaline under the assumption that all the salts precipitate at the eutectic point. We use the empirical equations of Anderson et al. (1973) as a starting point to relate the mass fraction of unfrozen water (relative to the total mass of the sediments mtotal) to the specific surface area of the sediment grains (Sa; in m2/kg):
where specific surface area Sa is related to the diameter (Dg; in m) and density (ρg; in kg/m3) of sediment grains via Sa=6/ρgDg; θ is relative temperature (in °C) in reference to the eutectic point in saline permafrost (θ=|TTeutectic|).
In equation 7, the total mass of the sediments mtotal can be expressed as the above-freezing mass sum of the sediment grains (mg) and pore water (mw0): mtotal=mg+mw0. By substituting in mg=ρgVolg=ρgVoltotal(1ϕgf) and mw0=ρw0Volw0=ρgVtotalϕgf (where Voltotal is the total volume of the saturated sediments and ϕgf is the grain frame porosity), we can rewrite the total mass mtotal as
Returning to equation 7, by replacing mtotal with equation 8 and the total mass of residual water mw with mw=ρwVolw, we obtain
After rearranging equation 9, we obtain the volume fraction of residual water relative to the total volume of the frozen sediments:
Note that in equation 10, the water density in the numerator (ρw0) corresponds to the brine density at the initial salinity prior to freezing, whereas the water density in the denominator (ρw) corresponds to the nonsaline water density at subeutectic temperature. By replacing Voltotal with Voltotal=Volpore/ϕgf and rearranging equation 10, we arrive at the expression of unfrozen water saturation:
Finally, because the water and ice saturations are related by sw+si=1, we arrive at the expression of ice saturation after replacing ww with equation 9:

In this section, we focus on modeling the velocities of saturated, unconsolidated saline permafrost as a function of temperature and pore-water salinity. The previously mentioned P-wave measurements from Dou et al. (2016) are used as the reference data set for model testing and calibration. Our procedure is based upon effective-medium approximation (EMA), which represents the physical properties of a composite by the average properties of its constituents. For EMA to be applicable, the composite must be statistically homogeneous. Simply put, the composite is microscopically heterogeneous but it appears homogeneous on a larger scale (Guéguen and Palciauskas, 1994). For this study, the size of the “larger scale” should be comparable with seismic wavelengths.

We organize this section into three parts according to the three groups of required inputs for EMA: (1) the elastic properties of the constituents, (2) the volume fractions of the constituents, and (3) the geometric details depicting the arrangements of the constituents.

Elastic properties of the constituents

Our EMA model requires the elastic properties of the three constituent phases, ice, brine, and the solid mineral grains. The ice phase is assumed to be isotropic and homogeneous. Its density ρi is determined using the expression from Pounder (1965) (equation A-1). Its bulk and shear moduli (Ki and Gi) are derived from the ρi expression and the P- and S-wave velocity measurements of polycrystalline ice (VPi and VSi) by Vogt et al. (2008) (equation A-2) (Ki=VPi2ρi4/3Gi and Gi=VSi2ρi). The parameters ρi, Ki, and Gi are all temperature-dependent (i.e., ρi(T), Ki(T), and Gi(T), where T is the temperature).

Relevant properties of the pore water are expressed as a function of temperature and salinity using the model of Batzle and Wang (1992). The bulk modulus Kw is derived from ρw and P-wave velocity expressions (VPw) via Kw=VPw2ρw. The ρw and Kw are temperature and salinity dependent (i.e., Kw(T,Sn) and ρw(T,Sn), where T is the temperature and Sn is the pore-water salinity), and ambient pressure (1 atm) is used in expressions that contain pressure terms.

Properties of the solid mineral grains are approximated by those of a single mineral phase or by the average properties of a multi-mineral mixture. Because the coarse sand samples are composed almost entirely of quartz, the elastic properties of the grains are approximated as those of pure quartz. The fine-grained saline permafrost core sample is constituted mainly of quartz, plagioclase, and clay minerals. Mixing laws are applied to obtain a single-phase equivalence of the mineral mixture, including the use of Hashin-Shtrikman (HS) average for elastic moduli (Hashin and Shtrikman, 1963) and arithmetic average for density. Table 2 lists the elastic properties of all the mineral phases that are used in the modeling.

Volume fractions of the constituents

The second set of key components of our permafrost EMA are the volume fractions among mineral grains, ice, and water. The volume fractions of mineral grains fgf are fgf=1ϕgf, where ϕgf is the porosity of the sediment grain frame. Here, we reemphasize the simplification of assuming constant ϕgf (and as a result, constant fgf) throughout the freeze/thaw processes. For the volume fractions of pore ice and pore water, we treat salinity-dominated and surface-effects-dominated regimes separately as described previously. Equations 6 and 12 are used to estimate ice saturations for these two regimes.

Arrangements of the constituents

For velocity modeling of unconsolidated permafrost, the pore-scale distribution of ice is the most crucial and the least constrained among all the elements involved in depicting the arrangements of the constituents. Because direct measurements of representative scales are rarely feasible, ice distributions need to be inferred indirectly through model evaluation. That is, we use rock-physics models that assume certain scenarios of ice distributions to predict the values of observables (e.g., seismic velocities). We then judge the validity of the assumptions based upon how well these models can predict the data.

Most of the existing models use ice distributions that are either of single or mixed types. The former assumes that during the initial growth of ice, only one type of ice distribution is present (Dvorkin et al., 1994, 1999; Dvorkin and Nur, 1996; Ecker et al., 1998; Helgerud et al., 1999; Helgerud, 2001). The latter assumes that more than one types of ice distributions can coexist, but the ratio among the various types often is a free parameter that needs to be calibrated (Chand et al., 2006; Minshull and Chand, 2009; Chand, 2013; Schindler et al., 2016). Although certain choices of the mixing ratio can achieve good data fits, it is difficult to verify the validity of the mixing ratio itself. Moreover, the calibrated mixing ratios often differ from data set to data set; and even within one data set, different levels of ice saturations often require different mixing ratios (e.g., Figure 1 in Chand, 2013). This ultimately weakens the applicability of these models.

Our model, as will be presented later in this section, assumes a mixed ice distribution. However, in contrast with many existing models that also assume mixed ice distributions (e.g., Carcione and Tinivella, 2000; Chand, 2013; Schindler et al., 2016), the mixing ratio in our model is not a free parameter. This is a desirable feature that makes the modeling procedure generalizable and simple to apply in scenarios lacking a calibration data set.

Brief overview of models with single ice distribution

In commonly used models with a single ice distribution, each portion of added ice is assumed to be only cementing, pore filling, or load bearing. Although in some of these models, ice distributions are allowed to transition from one type to another as the amount of ice increases, the assumed ice type for the freezing onset often dictates the overall characteristics of the velocity versus ice saturation trends. A good example is the extended contact cement theory (CCT) proposed by Dvorkin et al. (1999), in which an ice saturation threshold of approximately 30% acts as a transition point between cementing ice and pore-filling ice. As a result, this model predicts rapid increase of seismic velocities at the early stage of freezing and slower velocity increase for ice saturations of 30% and above. Despite the slowed velocity increase, however, velocities predicted by this model are generally high due to the stiffening effects of cementing ice assumed at the freezing onset.

This class of models, with their unambiguous assumptions about ice distributions, can provide useful insights into the pore-scale structure. Figure 4 shows two examples based upon the laboratory data of Dou et al. (2016), with the shaded areas denoting the partially frozen temperatures. Figure 4a corresponds to a high-salinity sample in which the slope of the P-wave velocity versus temperature trend remains moderate throughout the partially frozen stage. Figure 4b corresponds to a low-salinity sample in which the trend is steep near the onset of freezing and becomes more gradual as the temperature decreases further. Regardless of the pore-water salinity, none of the models assume that a single type of ice distribution can accurately fit this data set. In the partially frozen regime, the cementation models (models 1 and 2 in Figure 4) overestimate the velocities by up to 48%–80% and the pore-filling model (model 4 in Figure 4) underestimates the velocities by up to 15%–18%. The low velocities predicted by the load-bearing model (lower than the predictions of the pore-filling model for the high-salinity sample) appears surprising at first glance (model 3 in Figure 4) because the strengthening effect of the load-bearing ice is expected to be intermediate between the cementing ice and the pore-filling ice. The low velocities result from the very low confining pressure (approximately 1.16 kPa) that the samples were subjected to during the experiment. Because the strengthening effect of load-bearing ice has a strong stress dependency, the low effective stress reduces the influence of the load-bearing ice.

Noticeable velocity misfits are also present at subeutectic temperatures when the samples are approaching 100% ice saturation. Specifically, the extended CCT model underestimates the velocities by 71142  m/s and the load-bearing model overestimates the velocities by 356427  m/s.

The presence of velocity misfits shown in Figure 4 indicates that the class of models assuming a single type of ice distributions cannot accurately represent the seismic properties of saturated, unconsolidated saline permafrost. Given that the observed velocities are intermediate between the predictions of the cementing and pore-filling models, a natural step is to devise a model that includes both ice types.

Brief overview of models with mixed ice distribution

Several existing studies have explored various ways of mixing the cementing and pore-filling models to improve the representation of permafrost or hydrate-bearing sediments (Carcione and Tinivella, 2000; Chand, 2013; Schindler et al., 2016). Two good examples of such models are the partial cementation model (Chand et al., 2006; Minshull and Chand, 2009; Chand, 2013) and the isoframe model (Fabricius, 2003; Schindler et al., 2016). Although these two models differ in the procedure used to mix the cementing and pore-filling end members, both models assume some portion of the ice to be cementing and the remainder of the ice to be pore filling. The ratio between the ice types is a free parameter that requires fitting.

Besides requiring fitting to determine a mixing ratio, another source of uncertainty in such models lies in the difficulties in assessing the accuracy of the cementing and pore-filling end members themselves. Understandably, the performance of the end-member models dictates the performance of the mixture model. To assess the performance of the end-member models, model predictions should be compared against their observational counterparts. However, as we have demonstrated via examples in Figure 4, none of the velocity measurements can be directly related to cementing or pore-filling models. Without a sound observational basis, the accuracy of the end-member predictions cannot be evaluated, and consequently, the reliability of the mixture model appears unclear.

Minshull et al.’s (1994) two-end-member mixing model

A different model that assumes mixed ice distributions is proposed by Minshull et al. (1994). Instead of using cementing and pore-filling end members, this model uses fully unfrozen, water-filled sediments and fully frozen, ice-filled sediment as the soft and stiff end members. Such a setting has obvious benefits: Observational counterparts can be found for both end members (particularly for laboratory studies in which above-freezing and subeutectic temperatures can be reached), and thus one can check if the velocity predictions of the end members are sufficiently good before using them to model the velocities related with intermediate ice fractions.

Another advantage of this model is that the mixing ratio is not a free parameter. It directly equals the volume ratio between ice and unfrozen water, a quantity that can be determined via phase-diagram expressions of salt solutions for saline permafrost (as illustrated in the previous section “Estimating ice saturation”).

Despite having the above-mentioned merits, however, the original form of Minshull et al.’s (1994) model relies heavily on Wyllie’s slowness average (Wyllie et al., 1956) in the construction of the fully frozen end member and in the mixing strategy of the two end members. Unfortunately, slowness averaging is generally ineffective in representing unconsolidated sediments (Dvorkin and Nur, 1998).

However, the merits of Minshull et al.’s (1994) model outweigh its drawbacks. As we will soon illustrate, with modifications, it can become highly effective in modeling P-wave velocities of saturated, unconsolidated permafrost.

Before elaborating the modification steps, it is worthwhile to take a closer look at the key equations and procedures of the original model. The primary equation takes a similar form to Wyllie’s slowness average, but differently than the usual average between the fluid and mineral phases (1/V=sfluid/Vfluid+smineral/Vmineral), the average is taken between the water-filled, fully unfrozen end member and the ice-filled, fully frozen end member as follows:
where VWF and VIF are the P-wave velocities of fully unfrozen, water-filled end member and fully unfrozen, ice-filled end member, respectively, and sw and si are the relative volume fraction of water and ice (sw+si=1 for saturated permafrost).

Besides equation 13, the modeling of the fully frozen end member (VIF) invokes another use of the slowness average between mineral grains and ice: That is, 1/VIF=ϕgf/Vi+(1ϕgf)/Vm, where Vi and Vm are the velocities of ice and mineral grains and ϕgf is the porosity of the grain frame. At first glance, the second use of the slowness average seems justifiable given that the slowness average remains a frequently used model for consolidated sediments, and the ice-filled end member is expected to be as stiff as consolidated sediments. But a simple calculation quickly dismisses this notion: If we are to model the velocity observed at approximately 30°C from our coarse sand samples (VP approximately 4454  m/s; porosity ϕ=36%) by applying the slowness average to the velocities of ice and quartz (3922 and 6008  m/s, respectively), we obtain an estimated velocity of 5043  m/s, which amounts to an overestimation of 13%. Therefore, we conclude that the use of the slowness average should be avoided entirely.

An improved two-end-member mixing model

Once we identify the use of slowness average equations as a key limitation of Minshull et al.’s (1994) original model, we replace it with EMA procedures that are suitable for unconsolidated sediments. While keeping the choices of the end members (fully frozen and fully unfrozen sediments) and the values of the mixing ratios (equals the volume ratios between ice and unfrozen water), we change the construction of the end-member models and the EMA relationship used in mixing the two end members.

We approximate the fully frozen, ice-filled end member as a binary composite constituted of sediment grains and ice. Note that salt crystals are not considered as a component due to their negligible effects on the P-wave velocities. We seek solutions that accommodate both of the following two scenarios: (1) Some portions of the pore ice can be interconnected, in which case ice can be treated as a host medium and grains are the embedded inclusions and (2) some portions of the pore ice can be disconnected, in which case, the grains forming the host and ice are the embedded inclusions. The self-consistent approximation (SCA) (Berryman, 1995) suits this purpose because it yields an effective-medium equivalence of the composite that can self-consistently act as a host medium for the grain and ice inclusions. For the SCA of the fully frozen end member, we treat mineral grains as spherical inclusions and ice as penny-shaped inclusions. The aspect ratio of the ice inclusions is the only free parameter in this model, with smaller aspect ratios corresponding to lower velocities, and vice versa. Even though we cannot completely eliminate this free parameter, our sensitivity test shows that the plausible range of the aspect ratio is fairly narrow (approximately 0.001 and 0.06, as shown in Figure 5), and thus the required calibration is much less arbitrary. When the fully frozen state of the sample is attainable (e.g., in laboratory settings in which subeutectic temperatures are feasible), the aspect ratio could be obtained by calibrating the model predictions against the data. Otherwise, the Voigt upper bound (computed using the grain and ice properties) can serve as a useful reference for practitioners to choose a value that yields velocity predictions that are close to but lower than the Voigt velocities.

We model the unfrozen, water-filled end member as a random dense pack of sediment grains saturated with water (saline or nonsaline). We first estimate the elastic moduli of the dry granular pack with the Hertz-Mindlin contact theory (Mindlin, 1949). We then include the effects of pore water by using Biot’s fluid substitution equations for the high-frequency limiting velocities (Biot, 1956; Bouzidi and Schmitt, 2009). The validity test for the use of Biot’s high-frequency limiting velocity and for the choice of tortuosity value is detailed in Appendices  C and  D.

We model the partially frozen permafrost by mixing the above-mentioned two end members. As in Minshull et al.’s (1994) model, we equate the volume ratio of ice and water to the volume ratio of the ice-filled and water-filled end members. To compute the properties of the two-end-member mixture, we replace the original use of the slowness average with heuristic modifications of the HS average. That is, we first compute the modified upper and lower HS bounds by treating the two end members as two “phases.” We then compute the properties of the mixture by taking arithmetic average of the modified HS bounds. We summarize the workflow of the two-end-member mixing model in Figure 6. All of the associated equations are detailed in Appendix  E.

Microstructural implications

We can draw an intuitive understanding of the microstructure based upon heuristic interpretation of the modified HS bounds (Figure 6b): The upper bound is realized when the stiff end-member envelopes the soft end member, in which case the ice-filled, stiff end member is largely interconnected and thus the pore ice mostly acts as a frame-strengthening material (i.e., cementing and/or load bearing). The lower bound is realized when the soft end member surrounds the stiff end member, in which case the ice-filled, stiff end member is largely disconnected and thus the pore ice mostly acts as a pore-filling material. By taking the arithmetic average of the modified HS bounds, the coexistence of frame-strengthening ice and pore-filling ice is taken into account.

We can now compare the results of our velocity modeling with the previously described laboratory measurements of Dou et al. (2016). Our modeling was conducted for the coarse sand samples and the fine-grained permafrost core sample described in the paper. The modeling for the coarse sand samples can be viewed as a validity test for the underlying EMA because the simple mineral composition and texture of the samples (quartz sand with nearly uniform grain sizes) minimize uncertainties in estimating end-member properties. In contrast, the fine-grained core sample is more complicated in its mineral makeup and grain size distribution; thus, larger uncertainties are inevitable.

Modeling results: Coarse sand samples

For the coarse sand samples, the model-predicted P-wave velocities match well with the measured velocities, particularly for temperature ranges that are above freezing, subeutectic, and immediately below the freezing point (Figure 7). Table 3 shows the chosen value of the ice inclusion aspect ratio for each sample and the extent of the maximum velocity misfit. Overall, the misfit, even in the worst case, is less than 8%, and the aspect ratio values only need to be adjusted slightly within a narrow range that is consistent with our sensitivity test shown in Figure 5. This demonstrates the effectiveness of the underlying EMA, as well as confirming the coexistence of frame-strengthening ice and pore-filling ice in saturated, unconsolidated saline permafrost.

Modeling results: Fine-grained saline permafrost core sample

Due to the surface effects of the fine-grained particles, the saline permafrost core sample did not reach its fully frozen state during the experiment of Dou et al. (2016), despite being subjected to a subeutectic temperature of 35°C. This makes it impossible to directly calibrate the aspect ratio of the ice inclusions from the measurement of a fully frozen sample. Instead, the aspect ratio needs to be chosen empirically from a plausible range.

Another source of uncertainty lies in the selection of clay properties because published values vary broadly and lack consensus (Mondol et al., 2008). For instance, among the published bulk moduli of kaolinite, the lowest and highest values differ from each other by more than 30 fold (Table 4). As a result, choices of clay properties are inevitably subjective and the associated uncertainties are difficult to quantify.

We conduct the modeling in the form of a parameter sensitivity test targeted at clay properties and ice-inclusion aspect ratios. Starting from our prior knowledge concerning the reasonable range of these two parameters, we use the high and low extrema in the modeling. Namely, we compare the effects of using soft and stiff clay (Table 4), and soft and stiff ice inclusions (aspect ratio of 0.001 and 0.1, respectively). The resulting velocity predictions form bounds that both demonstrate how well the model can fit the data as well as how sensitive the model is to the corresponding parameters.

As shown in Figure 8, with the exception of “stiff clay” in Figure 8a, the rest of the model-predicted velocities show a good fit to the data in the partially frozen regime (with overall velocity mismatch no larger than 170  m/s). But for the subeutectic regime in which the temperature-dependent variations of ice saturations are mainly controlled by surface effects, all of the velocity predictions exhibit a noticeable deviation from the data (up to approximately 700  m/s). Such a mismatch is caused by the nonideal performance of the empirical relationship that is being used to estimate ice saturations for the surface-effects-dominated regime. However, because subeutectic temperatures are rarely encountered in nature, this mismatch does not impact the utility of the model.

Next, we examine the impact of the clay properties and inclusion aspect ratio to the velocity predictions of the partially frozen regime. In the temperature range of 2.5°C to 6°C (i.e., at the immediate vicinity of the freezing point), all of the model-predicted velocities can fit the data with a misfit no larger than 8%, demonstrating the robustness of the model when applied to the early stage of freezing. As temperatures decrease further, whereas the “soft clay” and “soft inclusion” cases yield a mismatch no larger than 8%, the stiff clay and “stiff inclusion” cases both produce a larger misfit (approximately 11% and 13%, respectively). Despite the uncertainties, however, the quality of the data fits remains superior to the examples shown in Figure 4, demonstrating the robustness of the model in spite of the larger uncertainties associated with the fine-grained core sample.

Up to this point, our modeling and analyses have been based upon velocity-versus-temperature (VP-T) trends. In this section, we move on to analyzing the relationship between velocities and ice content (VP-si, where si is the ice saturation), for which we replace the temperature measurements with the associated ice-saturation estimates. We focus only on the coarse sand samples (Figure 9) because the modeling results are free of uncertainties that are related to surface effects and clay properties.

In addition to our own measurements, we also use the ultrasonic measurements of brine-saturated quartz sand presented in Matsushima et al. (2016) as independent references. Besides seismic velocities, it is worth comparing ice saturation values between these two studies. Whereas our ice saturation values are estimated using the phase diagram of the NaCl-H2O solution, the ice saturation values in Matsushima et al. (2016) are derived from unfrozen water content measured with pulsed nuclear magnetic resonance (NMR). As we will soon illustrate, consistent discrepancies exist between these two sets of ice saturation values. Whereas we present model-predicted S-wave velocities VS and VP/VS ratios as functions of ice saturation, the experiments in Dou et al. (2016) did not measure VS. However, we do test VS and VP/VS predictions using the data set of Matsushima et al. (2016).

Analysis of VP-si trends

Once VP measurements are presented as functions of ice saturations, all data points from Dou et al. (2016), regardless of their associated pore-water salinities, collapse along a tight trend (Figure 9a). This has at least two implications: First, it is the variations in ice content that dictate the observed changes in VP, and, second, despite that it takes a different amount of time for samples of different salinities to reach the same level of ice saturations (i.e., faster for low-salinity samples and slower for high-salinity samples), the corresponding pore-scale distributions of ice are statistically similar among all the samples.

Comparing the observed VP-si trends against predictions derived from models that assume a single ice distribution (models 1, 2, 3, and 4 in Figure 9a), we once again see the inadequacy of such models in representing seismic properties of saturated, unconsolidated saline permafrost. In contrast, the VP-si trends predicted by the improved two-end-member mixing model provide a good fit to the data points over the entire range of the ice saturations that are covered by the measurements. Variations in the slope of the VP-si trends reflect the evolution of the pore ice distributions as the ice saturation increases:

  • For ice saturations between zero and approximately 0.6: The VP-si trends are nearly linear, indicating a relatively balanced partitioning between the pore-filling and frame-strengthening ice.

  • For ice saturations greater than 0.6: The VP-si trends show an increasing slope, indicating that the influences of the frame-strengthening ice are becoming more dominant. This is because the additional ice has fewer residual pore spaces available to grow freely. Instead, they often must “squeeze” into residual pore spaces and thus are more likely to be in contact with sediment grains and/or existing ice.

Analysis of model-predicted VP-si and VP/VS-si trends

Given that the S-wave cannot propagate through liquids, the S-wave properties are unaffected by pore fluids. They change mainly in response to changes occurring to the solid frame. For permafrost, how VS changes with respect to ice saturations is dictated by the presence and quantity of frame-strengthening ice.

We first examine the predictions derived from models that assume a single type of ice distribution. Because ice is assumed to be either frame strengthening or pore filling, the differences in VS and VP/VS predictions are more pronounced than those in VP predictions. This is particularly the case for cementing and pore-filling models. For cementing models (models 1 and 2 in Figure 9b), the predicted VS increases rapidly near the onset of freezing due to the frame-strengthening effects of the assumed cementing ice. In contrast, the pore-filling model (model 4 in Figure 9b) predicts VS to remain unchanged throughout the ice saturation range within which the pore-filling assumption is applicable (i.e., si=00.9). This contradicts existing laboratory observations (Zimmerman and King, 1986; Li, 2009; Matsushima et al., 2016), hence rendering the pore-filling model unsuitable for predicting S-wave velocities.

Situated in between the above-mentioned extreme cases are the VS-si predictions of the load-bearing model (model 3 in Figure 9b) and our two-end-member mixing model. Because cementing ice is completely absent in the load-bearing model, the predicted VS increase slowly as the ice saturation increases from zero to 0.7. On the other hand, the two-end-member mixing model predicts a VS-si trend that is generally steeper than that of the load-bearing model. This further confirms that the frame-strengthening ice assumed in the two-end-member mixing model inherently encompasses the cementing and load-bearing types.

Finally, we examine the model-predicted VP/VS-si trends (Figure 9c). The VP/VS eliminates the influence of density, and for its positive correlation to Poisson’s ratio (ν=1/2((VP/VS)22)/((VP/VS)21). The higher the VP/VS, the higher the Poisson’s ratio.), it can be viewed as an indicator of material deformability: the higher the VP/VS ratio, the more deformable the material. Given that liquids are more deformable than solids, a high VP/VS often also indicates high liquid saturations. In this way, VP/VS should be useful in distinguishing saline permafrost of higher unfrozen water content from a background of nonsaline permafrost that contains less water. However, this is not the case for cementing models (models 1 and 2 in Figure 9c) because the predicted VP/VS almost immediately falls to a low value of approximately 1.7 as soon as ice starts to form. The pore-filling model (model 4 in Figure 9c), on the other hand, predicts VP/VS to increase with increasing ice saturation (i.e., decreasing unfrozen water content), countering the expected positive correlation between water content and VP/VS. The two-end-member mixing model and the load-bearing model (model 3 in Figure 9c) predict VP/VS of the partially frozen sediment to decrease with increasing ice saturations while being consistently higher than that of the fully frozen permafrost (VP/VS=1.7). For the load-bearing model, VP/VS is less effective in distinguishing low-ice-content saline permafrost (si<0.15) from unfrozen sediments. For the two-end-member mixing model, VP/VS is less effective in distinguishing high-ice-content saline permafrost (si>0.7) from full frozen sediments.

Comparison with measurements of Matsushima et al. (2016)

In addition to the data from Dou et al. (2016), we also consider the ultrasonic measurements of Matsushima et al. (2016). Although the quartz sand samples used in Matsushima et al. (2016) have a smaller grain diameter (100–300 mm; medium sand) and higher initial porosity (41%), the overall sample properties and experiment methods are adequately similar for conducting meaningful comparisons.

In Figure 9, we present measurements of Matsushima et al. (2016) as functions of ice saturations. We use two sets of ice saturation values: One is the original values of Matsushima et al. (2016) that were derived from unfrozen water content measured with pulsed NMR. The other is our phase-diagram-based values that are estimated using the temperature measurements of Matsushima et al. (2016). Curiously, the two sets of ice saturation values persistently deviate by approximately 0.10–0.16 from each other. Despite such deviation, however, most of the (VP, si) points fall close to our measurements and the predictions of the two-end-member mixing model (Figure 9a).

One exception is the (VP, si) point at the ice saturation of 0.42. This corresponds to –2°C in Matsushima et al. (2016), a temperature that is below the phase-diagram predicted freezing point of –1.2°C. At this point, whereas the two-end-member mixing model considers the sample to be partially frozen, the VP measurement of Matsushima et al. (2016) indicates that the sample is fully unfrozen. This mismatch could have two possible causes: (1) The sample could be experiencing supercooling when the ultrasonic measurements were taken and (2) the freezing point of brine that is confined in porous media could be different from the freezing point of unconfined brine.

For VS and VP/VS values, the measurements of Matsushima et al. (2016) fall between load-bearing and two-end-member mixing model. Because the range of ice saturations sampled by these measurements is narrow (0.6–0.8 according to NMR-based measurements; 0.75–0.9 according to phase-diagram-based estimates), it is not possible to directly choose a better model between the two. Note that both models involve the Hertz-Mindlin contact theory, which is known to overestimate VS of unconsolidated sediments (Zimmer et al., 2006; Bachrach and Avseth, 2008). If a more accurate alternative to Hertz-Mindlin theory could be incorporated, we would expect both sets of model-predicted VS to be lowered. This would improve the fit of the two-end-member mixing while causing the VS predictions of the load-bearing model to deviate further from the data.

In our procedure for estimating ice saturation in saline permafrost, phase-diagram expressions of unconfined brines are assumed to be applicable to brines that are confined in porous media. For the partially frozen stage, all of the dissolved salts are always assumed to be fully rejected into the residual pore water as a result of freeze desalination and the total volume of the pore space is assumed to be unchanged. At above-eutectic temperatures, dissolved salts are assumed to be the only factor that controls changes in ice saturation. This assumption inherently requires that the initial pore-water salinity be sufficiently high so that free water is present at partially frozen temperatures. At the eutectic point, all the dissolved salts are assumed to precipitate out of the residual pore water. At subeutectic temperatures, the effects of the precipitated salt crystals are considered negligible and all the residual pore water is assumed to be nonsaline, rendering all the changes in water/ice fractions to be controlled entirely by surface effects. Because no surface-effects equations are universally applicable, empirical equations are used and the resulting ice saturation estimates are subjected to larger uncertainties.

For the two-end-member mixing model, the key limitations are the highly idealized grain/pore geometries that are assumed when constructing the two end members, the heuristic means of determining the mixing ratio, and the use of the modified HS average for mixing the end members. For the fully unfrozen, water-filled end member, the use of Hertz-Mindlin contact theory approximates the sediment grain frame as a random dense pack of identical spheres, and the best tortuosity value used in computing Biot’s high-frequency limiting velocity assumes highly idealized pore shapes (uniform cylindrical pores with axes parallel to the pore pressure gradient). Other assumptions required by the use of Biot’s velocity predictions are the sediment is isotropic, the pore fluid is Newtonian, the wavelength is much larger than the grain or pore scale, only global flow (i.e., fluid flow resulting from wavelength-scale pressure gradients) is considered, and local or pore-scale squirt flow is neglected. For the fully frozen, ice-filled end member, the mineral grains, and pore ice are given idealized shapes: Mineral grains are approximated as spherical inclusions, and pore ice is approximated as penny-shaped inclusions. The aspect ratio of the ice inclusions is assumed uniform, and it is the sole free parameter (and thus an empirical component) of the model that requires calibration. For the step of mixing the two end members, equating the mixing ratio of the fully frozen and the fully unfrozen end members to the volumetric ratio of ice and unfrozen water is a heuristic approach that has no theoretical basis. Another heuristic element is the modified HS average. Although the HS bounds are rigorous when the mixing ingredients are mineral phases, they are no longer rigorous when the mixing ingredients are effective-medium composites. Although similar microstructural interpretations of the HS upper and lower bounds are intuitively applicable (i.e., the upper bound is realized when the stiffer ingredient forms shells that embrace the softer ingredient, and vice versa), the modified HS bounds cannot be derived from first principles.

Finally, due to the large variations in elastic properties of clay minerals and possibly stronger influences of surface effects, the two-end-member mixing model is not applicable to clay-rich sediments. In addition, note that by far the performance of the model is mainly examined against P-wave data that were acquired under ambient pressures. The pressure-dependent properties of the model and its effectiveness in predicting S-wave velocities are yet to be further investigated.

In this study, we conduct rock-physics modeling to quantitatively relate P-wave velocities to ice content in saturated, unconsolidated saline permafrost. Reference data used to verify the models are the laboratory measurements obtained from our earlier studies. Poor data fits of commonly used models, which assume a single type of ice distribution, indicate that neither frame-strengthening nor pore-filling ice alone is adequate in depicting realistic pore-scale distributions of ice. Instead, coexistence of the frame-strengthening and pore-filling ice must be taken into account via a generalizable mixing strategy, in which ad hoc tuning of the mixing ratio should be avoided.

To model mixed ice distributions, we use a two-end-member mixing workflow that does not require ad hoc fitting of the mixing ratio. Using a modified HS average to mix a water-filled, fully unfrozen end member with an ice-filled, fully frozen end member, the resulting EMA of the partially frozen sediments inherently include frame-strengthening and pore-filling ice. The mixing ratio of the two end members is equated to the volumetric ratio of unfrozen water and ice, hence eliminating the need for tuning the mixing ratio.

The model still contains one free parameter, which is the aspect ratio of the penny-shaped ice inclusions in the EMA of the fully frozen end member. However, our sensitivity tests show that the aspect ratio was typically in a narrow range for the samples considered (e.g., approximately 0.001–0.06 for coarse sand samples); hence, the model can be used with a reasonable estimate for this parameter. The quality of the fits of the two-end-member mixing model confirms the validity of the mixed ice distributions as well as the effectiveness of the mixing strategy.

As part of the Next-Generation Ecosystem Experiments (NGEE-Arctic) project sponsored by the Office of Biological and Environmental Research in the DOE Office of Science, this study is supported through contract DEAC0205CH11231 to the Lawrence Berkeley National Laboratory and through contract DE-AC05-00OR22725 to the Oak Ridge National Laboratory. We thank S. Hubbard (Lawrence Berkeley National Laboratory) and S. Wullschleger (Oak Ridge National Laboratory) for their leadership within the NGEE-Arctic program. We also thank our three reviewers — M. Schindler, D. Schmitt, and one anonymous reviewer, assistant editor L. Valentina Socco, and associate editor M. Asten for their constructive comments.


Velocities, density, and elastic moduli of ice

  1. 1)
    The P- and S-wave velocities of ice VPi and VSi (in m/s) as functions of temperature T (in °C) are given by Vogt et al. (2008) as
  2. 2)
    The density of ice ρi (in kg/m3) as a function of temperature T (in °C) is given by Pounder (1965)
  3. 3)
    The P-wave moduli of ice (Ki, Gi, and Mi; in Pa) as functions of temperature are obtained using the relationship among velocities, densities, and moduli:

P-wave velocity, density, and bulk modulus of water (saline or nonsaline)
We compute the P-wave velocity VPw (in m/s) and density ρw (in kg/m3) of water based upon the equations of Batzle and Wang (1992). We then obtain the bulk modulus Kw (in Pa) of water as a function of temperature T (in °C) and salinity Sn (in wt%):


To verify that precipitated solid salts can be neglected in velocity modeling of the fully frozen end member, we compile subeutectic measurements (near 27°C) of the coarse sand samples in Figure B-1 and examine if any systematic trends exist as a function of salt concentrations. The absence of correlations between the observed P-wave velocities of the fully frozen sample and the initial salinities confirms that the effects of solid salts are negligible.


To examine the validity of using Biot’s high-frequency limiting velocity, we start from the formula of the reference frequency fc that determines the low-frequency range (ffc) and the high-frequency range (ffc):
where ϕ is the grain-frame porosity, κ is the grain-frame permeability, η is the fluid viscosity, and ρfl is the fluid density.

The viscosity and density of brine increase with increasing salinity and decreasing temperature. But because the viscosity increases exponentially with decreasing temperatures, its influence outweighs that of density at low temperatures, causing fc to shift toward higher frequencies. In extreme cases, fc can become so high that it falls into typical ultrasonic frequencies, hence rendering high-frequency limiting velocities inapplicable.

We use the upper and lower limits of the brine viscosity to test if the high-frequency limiting assumption remains valid for the dominant frequencies of our laboratory data (150–750 kHz). The room temperature viscosity of pure water (0.001 Pa·s at 20°C) is used as the lower limit reference. The eutectic-point viscosity of saturated brine estimated by Gaidos and Marion (2003) (0.014 Pa·s) is used as the upper limit. Figure C-1 shows the corresponding Biot dispersion relation in the complete frequency range. The dominant frequencies of Dou et al. (2016) (150–750 kHz) are markedly higher than the highest reference frequency (fc=2.12  kHz), demonstrating the validity of using Biot’s high-frequency limiting velocity for the fully unfrozen end member.


Biot’s high-frequency limiting velocity has a strong dependency on tortuosity (denoted as τ—defined as the ratio between the total flow-path length and the sample length). To select an appropriate tortuosity value, we conduct a sensitivity test using Berryman’s tortuosity equation (Berryman, 1981):
where ϕ is the porosity; r is a pore-shape parameter that varies between zero and one (e.g., for spherical pores, r=1/2; for uniform cylindrical pores with axes parallel to the pore pressure gradient, r=0). In response to changes in r, the value of tortuosity τ varies between one and three. Figure D-1 shows that the r value of zero (and thus the lowest possible tortuosity of τ=1) provides the best data fits for the unfrozen coarse sand samples.


The fully frozen, ice-filled end member
The fully frozen end member is approximated as a composite consisting of spherical sediment grains and penny-shaped pore ice. The effective moduli of the fully frozen end member are calculated based upon the SCA (Berryman, 1995). Note that the SCA is indiscriminate toward the constituents of the composite. Hence, the self-consistent effective moduli are derived from the combination of two terms and must satisfy
where {K1, G1} and {K2, G2} represent the following two scenarios:
  • Ice as penny-shaped inclusions embedded in the self-consistent effective medium
    where fi is the volume fraction of ice inclusions, ϕgf is the porosity of the grain frame, Ki is the bulk modulus of ice, Gi is the shear modulus of ice, K* is the self-consistent effective bulk modulus, G* is the self-consistent effective shear modulus, and the P*ipenny-inclusion and Q*ipenny-inclusion coefficients are calculated with ice as penny-shaped inclusions and the effective medium as the background material:
    with αi being the aspect ratio of the penny-shaped ice inclusions and β*=G*((3K*+G*)/(3K*+4G*)).
  • Sediment grains as spherical inclusions embedded in the self-consistent effective medium
    where fg is the volume fraction of sediment grains, Kg is the bulk modulus of the sediment-grain mineral, Gg is the shear modulus of the sediment-grain mineral, and P*gsphere and Q*gsphere coefficients are calculated with sediment grains as spherical inclusions and the effective medium as the host material:
    with ζ*=G*/6((9K*+8G*)/(K*+2G*)).

Then, a combination of equations E-1, E-2, and E-4 forms a set of coupled equations that are used to iteratively solve for the self-consistent effective moduli K* and G*, starting from an initial guess for K*0 and G*0.

In the end, the effective moduli of the fully frozen end member (KIF and GIF; as the effective moduli of the ice-filled, stiff end member) equal the optimal solution of K* and G*:
The fully unfrozen, water-filled end member

The fully unfrozen end member is approximated as a random dense pack of sediment grains saturated with water (saline or nonsaline).

  • The effective moduli of the dry granular frame Kgf and Ggf are first calculated with Hertz-Mindlin contact theory (Mindlin, 1949):
    where νg=(3Kg2Gg)/(2(3Kg+Gg)) is the Poisson’s ratio of the sediment grain; P is the effective pressure (P=(ρbulkρw)gh, where ρbulk is the bulk density of the water-saturated sample, ρw is the water density, h is the height of the sample above the measurement position, and g is the gravitational acceleration); and C is the coordination number (the average number of contacts that each grain has with surrounding grains), and it is related to the grain frame porosity ϕgf based upon a cubic spline interpolation of empirical coefficients that were compiled by Murphy (1982) (it can be found in Table 5.1.5 in Mavko et al., 2009). As an example, for our coarse-grained Ottawa sand sample with a pore-water salinity of 35 ppt (0.6 M), an h value of 7 cm yields effective pressure P=711  Pa and dry grain-frame moduli Kgf=0.07GPa and Ggf=0.10  GPa.
  • The effective moduli of the water-filled granular pack KWF and GWF are then calculated based upon Biot’s fluid substitution equations for the high-frequency limiting velocities:
    where the average density of the water-filled granular pack is ρWF=ρg(1ϕgf)+ρwϕw and the high-frequency limiting velocities are given by (Biot, 1956)
    where τ is the tortuosity, and the other coefficients are expressed as follows:

The two-end-member mixing approach
With the elastic moduli of the stiff, ice-filled end member {KIF, GIF} and the soft, water-filled end member {KWF, GWF}, the effective moduli of saline permafrost are calculated based upon an arithmetic average of the modified HS bounds (Hashin and Shtrikman, 1963) of the two end members for the entire range of water/ice saturations
where the modified HS upper bounds are given by
and the modified HS lower bounds are given by
and the volume fractions of the two end members are given by fIF=si and fWF=sw, where si and sw are the water and ice saturation that are estimated using equations 11 and 12.
Then, P- and S-wave velocities of saline permafrost for the entire range of water/ice saturations are expressed as
with the average density ρ=ρg(1ϕgf)+ρiSiϕgf+ρwSwϕgf.
Freely available online through the SEG open-access option.