An effective-medium model for P-wave velocities of saturated, unconsolidated saline

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. 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 ( V P 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. correlation between the initial salinity of the solution ( S n0 Sn0) and its freezing point ( T fp Tfp) lying on the liquidus (we denote it as “the T fp Tfp versus S n0 Sn0 interpretation”): the higher the initial salinity, the the freezing point. In this sense, the liquidus quantitatively depicts the freezing-point depression effect of the dissolved salts. the lowering of the freezing point is limited by a critical temperature known the point ( approximately O desalinated 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. 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. The upper is 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. 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 −2.5°C to −6°C −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. V P /V S 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 V S 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 V S VS to be lowered. This would improve the fit of the two-end-member mixing while causing the V S VSpredictions of the load-bearing model to deviate further from the data. 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 −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. 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 ( f c =2.12kHz fc=2.12 kHz), demonstrating the validity of using Biot’s high-frequency limiting velocity for the fully unfrozen end member.

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(Nelson et al., , 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 porefilling 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 ( V P 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 ( V P VP, T T) points (i.e., P-wave velocities V P VP and the associated temperatures T T). To analyze these data sets using established rock-physics relationships between V P VPand ice saturation s i si ( s i =Vol ice /Vol pore 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 T-to-s i siconversion). For controlled freezing and thawing of saline permafrost in which NaCl is the predominant salt, the T T-tos i siconversion is straightforward for two reasons: First, NaCl-H 2 O H2O brine has phase-diagram expressions that are well-established and readily usable for the T T-to-s i si conversion. Second, because the water-to-ice phase transitions are progressive in the presence of dissolved salts, the associated changes in V P VP occur over a broad temperature range. This allows the V P VP versus T 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 T-tos i si conversion. Brief overview of the NaCl-H 2 O H2O phase diagram We first review the nomenclatures and expressions of the NaCl-H 2 O 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.

BACKGROU ND
As shown in Figure 1, the liquidus on the NaCl-H 2 O H2O phase diagram is the temperature-versus-salinity ( T T-S n 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 ( S n0 Sn0) and its freezing point ( T fp Tfp) lying on the liquidus (we denote it as "the T fp Tfp versus S n0 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 ( T eutectic Teutectic; approximately −21°C −21°C for NaCl-H 2 O H2O brine). Once the temperature drops below T eutectic 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.
View larger version (55K) Figure 1. Phase diagram of the NaCl-H2OH2O solution.
An alternative to the T fp Tfp versus S n0 Sn0 interpretation of the liquidus is the equilibrium salinity interpretation. At a given temperature within the partially frozen regime ( T [pf] ;T eutectic <T [ (2 ) 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  • 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−21°C down to −31.8°C−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.
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 freezingpoint 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.
View larger version (78K) Figure 3. Schematic illustration of the two-regime division assumed in estimating ice saturations. Sand, silt, and clay are approximated using quartz grains of descending diameters (700, 40, and 3μm3 μm, respectively). Pore water is an NaCl-H2OH2O solution with salinity of 0.7 M (=42ppt=42 ppt): TeutecticTeutectic, eutectic point (−21°C−21°C); SnSn, porewater salinity. 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 salinitydominated 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.

ESTIMATIN G ICE SATURATIO NS
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, s i =Vol ice /Vol pore si=Volice/Volpore, we assume that the total volume of the pore space Vol pore 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-H 2 O 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 ( T fp Tfp) that is specific to a given initial salinity ( S n0 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 m salt msalt is the total mass of the dissolved salts, m w0 mw0 is the initial mass of the pore water prior to freezing, m w mw is the mass of the residual pore water at partially frozen temperatures, and S n Sn is the equilibrium salinity that is determined by equation 2.
Finally, by replacing mass m m with the product of density ρ ρ and volume Vol ( m=ρVol m=ρVol) in equation 3, we obtain Sn0ρw0Volw0=SnρwVolw,Sn0ρw0Volw0=SnρwVolw, where ρ w0 ρw0 and ρ w ρ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 ( Vol pore =Vol w0 Volpore=Volw0), we rearrange equation 4 and obtain the saturation of the residual brine as sw=VolwVolpore=VolwVolw0=ρw0Sn0ρwSn.sw=VolwVolpore=VolwVolw0=ρw0Sn0ρwSn.

(5)
Because the saturation values of water and ice are related by s w +s i =1 sw+si=1, we obtain the equation for determining ice saturation after replacing S n Sn with equation 2: In equation 6, brine densities ρ w0 ρw0 and ρ w ρ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 m total mtotal) to the specific surface area of the sediment grains ( S a Sa; in m 2 /kg m2/kg): where specific surface area S a Sa is related to the diameter ( D g Dg; in m) and density ( ρ g ρg; in kg/m 3 kg/m3) of sediment grains via S a =6/ρ g D g Sa=6/ρgDg; θ θ is relative temperature (in °C) in reference to the eutectic point in saline permafrost ( θ=|T−T eutectic | θ=|T−Teutectic|). In equation 7, the total mass of the sediments m total mtotal can be expressed as the above-freezing mass sum of the sediment grains ( m g mg) and pore water ( m w0 mw0): m total =m g +m w0 mtotal=mg+mw0. By substituting in m g =ρ g Vol g =ρ g Vol total (1−ϕ gf ) mg=ρgVolg=ρgVoltotal(1−ϕgf) and m w0 =ρ w0 Vol w0 =ρ g V total ϕ gf mw0=ρ w0Volw0=ρgVtotalϕgf (where Vol total Voltotal is the total volume of the saturated sediments and ϕ gf ϕgf is the grain frame porosity), we can rewrite the total mass m total mtotal as mtotal=(ρg(1−ϕgf)+ρw0ϕgf)Voltotal.mtotal=(ρg(1−ϕgf)+ρw0ϕgf)Voltotal.

(10)
Note that in equation 10, the water density in the numerator ( ρ w0 ρw0) corresponds to the brine density at the initial salinity prior to freezing, whereas the water density in the denominator ( ρ w ρw) corresponds to the nonsaline water density at subeutectic temperature. By replacing Vol total Voltotal with Vol total =Vol pore /ϕ gf 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 s w +s i =1 sw+si=1, we arrive at the expression of ice saturation after replacing w w 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 ρi is determined using the expression from Pounder (1965) (equation A-1). Its bulk and shear moduli ( K i Ki and G i Gi) are derived from the ρ i ρiexpression and the P-and S-wave velocity measurements of polycrystalline ice ( V Pi VPi and V Si VSi) by Vogt 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 K w Kw is derived from ρ w ρw and P-wave velocity expressions ( V Pw VPw) via K w =V 2Pw ρ w Kw=VPw2ρw. The ρ w ρw and K w Kw are temperature and salinity dependent (i.e., K w (T,S n ) Kw(T,Sn) and ρ w (T,S n ) ρw(T,Sn), where T Tis the temperature and S n 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 singlephase 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.
View Larger Version Table 2. Density and elastic moduli of sediment-grain minerals. 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 f gf fgf are f gf =1−ϕ gf fgf=1−ϕgf, where ϕ gf ϕgf is the porosity of the sediment grain frame. Here, we reemphasize the simplification of assuming constant ϕ gf ϕgf (and as a result, constant f gf 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., 1994Dvorkin 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.
View larger version (42K) . 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, icefilled 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 . 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=s fluid /V fluid +s mineral /V mineral 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 V WF VWF and V IF VIF are the P-wave velocities of fully unfrozen, water-filled end member and fully unfrozen, ice-filled end member, respectively, and s w sw and s i si are the relative volume fraction of water and ice ( s w +s i =1 sw+si=1 for saturated permafrost). Besides equation 13, the modeling of the fully frozen end member ( V IF VIF) invokes another use of the slowness average between mineral grains and ice: That is, 1/V IF =ϕ gf /V i +(1−ϕ gf )/V m 1/VIF=ϕgf/Vi+(1−ϕgf)/Vm, where V i Vi and V m Vm are the velocities of ice and mineral grains and ϕ gf ϕ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 −30°Cfrom our coarse sand samples ( V P VP approximately 4454m/s 4454 m/s; porosity ϕ=36% ϕ=36%) by applying the slowness average to the velocities of ice and quartz (3922 and 6008m/s 6008 m/s, respectively), we obtain an estimated velocity of 5043m/s 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 highfrequency 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.
View larger version (65K) Figure 6. Schematic illustration of the improved two-end-member mixing model. 1, Effective-medium approximation of the fully frozen, ice-filled end member as a self-consistent mixture of spherical grains and penny-shaped pore ice; 2, EMA of the fully unfrozen, water-filled end members as a random dense pack of identical spherical grains saturated with water (saline or nonsaline); and 3, mixing of the two end members via a modified HS average. TfpTfp, freezing point; TeutecticTeutectic, eutectic point.

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 framestrengthening 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.

.CITING ARTICLES
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.   Table 2.

View Larger Version
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  (Table 4). As a result, choices of clay properties are inevitably subjective and the associated uncertainties are difficult to quantify.
View Larger Version Table 4. Elastic properties of clay minerals used in velocity modeling of the fine-grained saline core sample. 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 170m/s 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 700m/s 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.
View larger version (32K) Figure 8. Data fits of the improved two-end-member mixing model for P-wave velocities of fine-grained saline permafrost core sample. (a) Effects of uncertainties in clay properties (with the ice inclusion aspect ratio fixed at 0.01) and (b) effects of uncertainties in ice inclusion aspect ratios (using "soft" clay parameters). Soft and "stiff" clay properties are listed in Table 4. Soft and stiff ice inclusions have aspect ratios of 0.001 and 0.1, respectively. The unfilled circles denote the laboratory data of Dou et al. (2016). The solid and dashed lines denote the model predictions.
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 −2.5°C to −6°C −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.

.CITING ARTICLES
Up to this point, our modeling and analyses have been based upon velocity-versus-temperature ( V P VP-T T) trends.
In this section, we move on to analyzing the relationship between velocities and ice content ( V P VP-s i si, where s i 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. Comparing the observed V P VP-s i 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 V P VP-s i 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 V P VP-s i 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 VPVP-sisi trends are nearly linear, indicating a relatively balanced partitioning between the pore-filling and framestrengthening ice. • For ice saturations greater than 0.6: The VPVP-sisi 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 V P VP-s i si and V P /V S VP/VS-s i 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 V S 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 V S VS and V P /V S VP/VS predictions are more pronounced than those in V P 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 V S 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 V S VS to remain unchanged throughout the ice saturation range within which the pore-filling assumption is applicable (i.e., s i =0-0.9 si=0-0.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 Swave velocities.
Situated in between the above-mentioned extreme cases are the V S VS-s i 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 V S 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 V S VS-s i 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 V P /V S VP/VS-s i si trends (Figure 9c). The V P /V S VP/VS eliminates the influence of density, and for its positive correlation to Poisson's ratio ( ν=1/2((V P /V S ) 2 −2)/ ((V P /V S ) 2 −1) ν=1/2((VP/VS)2−2)/((VP/VS)2−1). The higher the V P /V S VP/VS, the higher the Poisson's ratio.), it can be viewed as an indicator of material deformability: the higher the V P /V S VP/VS ratio, the more deformable the material. Given that liquids are more deformable than solids, a high V P /V S VP/VSoften also indicates high liquid saturations. In this way, V P /V S 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 V P /V S 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 V P /V S VP/VS to increase with increasing ice saturation (i.e., decreasing unfrozen water content), countering the expected positive correlation between water content and V P /V S VP/VS. The two-end-member mixing model and the load-bearing model (model 3 in Figure 9c) predict V P /V S VP/VS of the partially frozen sediment to decrease with increasing ice saturations while being consistently higher than that of the fully frozen permafrost  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 porescale 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.

CONCLUSIONS
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.  2 )

ACKNOWLEDGMENT S
The density of ice ρiρi (in kg/m3kg/m3) as a function of temperature TT (in °C) is given by Pounder (1965) ρi (  P-wave velocity, density, and bulk modulus of water (saline or nonsaline) We compute the P-wave velocity V Pw VPw (in m/s m/s) and density ρ w ρw (in kg/m 3 kg/m3) of water based upon the equations of Batzle and Wang (1992). We then obtain the bulk modulus K w Kw (in Pa) of water as a function of temperature T T (in °C) and salinity S n Sn (in wt%): Kw(T,Sn)=VPw(T,Sn)2ρw(T,Sn).Kw(T,Sn)=VPw(T,Sn)2ρw(T,Sn).

APPENDIX BNEGLIGIB LE EFFECTS OF SOLID SALTS ON P-WAVE
where ϕ ϕ is the grain-frame porosity, κ κ is the grain-frame permeability, η η is the fluid viscosity, and ρ fl ρ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 f c fc to shift toward higher frequencies.
In extreme cases, f c 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  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 r is a pore-shape parameter that varies between zero and one (e.g., for spherical pores, r=1/2 r=1/2; for uniform cylindrical pores with axes parallel to the pore pressure gradient, r=0 r=0). In response to changes in r r, the value of tortuosity τ τ varies between one and three. Figure D-1 shows that the r r value of zero (and thus the lowest possible tortuosity of τ=1 τ=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 pennyshaped 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 K1+K2=0,G1+G2=0,K1+K2=0,G1+G2=0, (E-1)
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 * K* and G * G*, starting from an initial guess for K *0 K*0 and G *0 G*0. In the end, the effective moduli of the fully frozen end member ( K IF KIF and G IF GIF; as the effective moduli of the icefilled, stiff end member) equal the optimal solution of K * K* and G * G*: KIF=K * ,GIF=G * .KIF=K*,GIF=G*.

(E-6)
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 KgfKgf and GgfGgf are first calculated with Hertz-Mindlin contact theory (Mindlin, 1949): where νg=(3Kg−2Gg)/(2(3Kg+Gg))νg=(3Kg−2Gg)/(2(3Kg+Gg)) is the Poisson's ratio of the sediment grain; PP is the effective pressure (P=(ρbulk−ρw)ghP=(ρbulk−ρw)gh, where ρbulkρbulk is the bulk density of the water-saturated sample, ρwρw is the water density, hh is the height of the sample above the measurement position, and gg is the gravitational acceleration); and CC 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ϕgf based upon a cubic spline interpolation of empirical coefficients that were compiled by Murphy (1982) (it can be found in Table 5 where the average density of the water-filled granular pack is ρWF=ρg(1−ϕgf) +ρwϕwρWF=ρg(1−ϕgf)+ρwϕw and the high-frequency limiting velocities are given by (Biot, 1956)