Growth rate dependence of the permeability and percolation threshold of young sea ice

Abstract

The permeability of sea ice is difficult to observe, and physically based permeability models are lacking so far. Here a model for the permeability of sea ice is presented that combines extensive microstructure observations and modelling with directed percolation theory. The model predicts the dependence of sea ice permeability on brine porosity and growth rate, as well as a percolation transition to impermeable sea ice due to necking of the pores. It is validated by numerical simulations of sea ice permeability on 3D images from X-ray microtomographic imaging and by other existing permeability data. A fundamental model result is that the percolation threshold of sea ice scales as ϕ_c ∝ a₀⁻¹ where a₀ is the plate or brine layer spacing. As the plate spacing decreases with growth velocity V, this implies that the percolation threshold increases as ϕ_c ∝ V^1/3, with the cubic root of the growth rate. For growth rates of natural sea ice the percolation threshold is expected to be in the range of 1 to 4 percent volume fraction of brine. While developed for columnar sea ice, a simple modification for granular surface ice also agrees with observations. The model is valid for sea ice during the growth phase, prior to warming and melting. Permeability modelling of spring and summer sea ice, with wider secondary brine channels present, requires 3D pore space observations in warming sea ice that currently are sparse.

---

1 Introduction

The hydraulic permeability reflects the ability of sea ice to transport fluid through its pore space. It is an important property that controls several sea ice processes and their role in the coupled atmosphere–ice–ocean system. During the melt season, the permeability controls the drainage rate of melt ponds and evolution of its surface properties, like albedo and roughness. During the growth season, permeability controls the salt loss from sea ice and thus is critical for the evolution of sea ice salinity. The latter, in turn, determines many sea ice properties, such as its mechanical, thermal and radiative properties.¹ To predict the salinity and desalination of sea ice, models of different complexity have been suggested, of which the most recent approaches involve a parametrisation of the permeability.^2–8 The permeability is also important for the biogeochemistry and ecosystem of sea ice,^9–11 for air–ice–ocean gas exchange^11–14 and for brine transport to the sea ice surface, from where brine freezing and chemistry may create interactions with atmospheric chemistry as far up as the stratosphere.^15–18

The sea ice permeability is related to its brine volume fraction ϕ. As brine in disconnected inclusions does not contribute to the permeability, knowing the critical brine volume fraction ϕ_c, at which all brine is disconnected, is important. However, understanding this transition and its dependence on growth conditions and microstructure is a challenging problem. Several authors have proposed percolation theory to address it, and have termed the minimum the percolation threshold of sea ice. The first studies on this problem proposed a threshold brine volume fraction of 5 percent.^3,19,20 However, in a recent analysis based on 3D X-ray microtomography of young sea ice, a much lower threshold of 2–3 percent has been found.²¹ In the present study, I extend the latter²¹ results by formulating a new model for the percolation threshold and permeability of sea ice that accounts for microstructure and growth conditions. Based on this model, the differences in permeability observations are explained and most earlier work on sea ice permeability is revised.

2 Permeability – definition and previous work

The permeability K (with unit m²) of a porous medium is defined via Darcy’s equation:^22,23with discharge per unit area W̄, dynamic viscosity μ and pressure gradient dP/dz. Eqn (1) is valid for flow in one direction, while the permeability in 3D coordinates is a tensor. The present study focuses on the vertical permeability of sea ice, for which dP/dz may be created by meltwater on top of the ice or the upward brine salinity increase. The latter is the primary driver of gravity drainage and desalination of sea ice.

The most widely adopted parametrisation of permeability K of sea ice has been proposed by Freitag,²⁴ on the basis of laboratory experiments. It may be written as:

K = 2.00 × 10⁻⁸ϕ^3.1 m², (2)

where ϕ is given as a fraction. According to the data range from ref. 24, this approximation is valid in the regime 0.1 < ϕ < 0.3 and thus reflects younger and warmer ice.

To address lower brine volume fractions, the permeability of sea ice has been discussed on the basis of percolation theory.^3,5,19,20 In these theories, the sea ice permeability becomes zero at some critical brine porosity, for which most studies indicated a threshold of ϕ_c ≈ 0.05. Maus et al.²¹ have revisited these studies and concluded that there is little observational basis for a threshold of ϕ_c ≈ 0.05. Based on numerical simulations on X-ray microtomographic images of young sea ice, they proposed the following equation for the permeability:

K = 1.49 × 10⁻⁸(ϕ − ϕ_c)^2.55 m², (3)

as well as the result ϕ_c = 0.024 ± 0.04, which is roughly half that of earlier estimates.

In Fig. 1, the results from ref. 21 and 24 are compared. In view of the wide range of observed permeability in other studies,^24–26 they appear surprisingly close. However, both represent young sea ice of similar thickness (20–35 cm) that had grown for a few weeks. The concept to obtain the permeability was also similar in the studies, as both are based on centrifuged sea ice samples. In the first study,²⁴ this provided the basis to measure the permeability in laboratory experiments with decane. In the second study,²¹ the centrifuged brine provided enough imaging contrast to obtain 3D X-ray microtomographic images for numerical simulations. For other growth conditions, no similar permeability observations are currently available. Also, little data exist at high brine volume fractions above ϕ = 0.20. This includes, for growing sea ice, the skeletal layer near the ice–water interface. In that regime, typically a few centimeters thick, the convective exchange with the ocean is largest. For columnar sea ice, it is characterised by a lamellar structure of vertically oriented plates that are parallel within each grain. For a brine layer of thickness d between the ice lamellae, the permeability is then given by Hele-Shaw flow as d²/12, such that the bulk permeability becomes K_d = ϕ/12d². With distance from the interface, the ice lamellae thicken due to heat flow and two processes. First, the temperature decreases, and the brine adjusts to a new thermodynamic equilibrium of higher solute concentration. Second, convective exchange with seawater underneath progressively decreases the solute concentration, also increasing the ice solid fraction. Disregarding what exactly happens, this change in brine layer thickness is given as:

which results in a simple analytic expression:for the permeability. It has a similar brine porosity exponent as eqn (2), yet there is an important difference. The conductivity in eqn (5) depends on the plate spacing a₀, which is known to depend on ice growth conditions, to be discussed below. The three relationships between porosity and permeability for young sea ice are compared in Fig. 1, drawing the curves for the respective validity ranges. Eqn (5) is drawn for the plate spacing range of 0.4 to 0.7 mm typical for young sea ice;^27–29 see below. The figure distinguishes three porosity regimes separated by two critical porosities:

(1) ϕ > ϕ₀ – the brine layer regime. In the lamellar regime above ϕ₀, the brine layer width is given by the product of brine porosity ϕ and plate spacing a₀.

(2) ϕ_c < ϕ < ϕ₀ – the percolation regime. The upper critical porosity ϕ₀ marks the transition at which the morphology of sea ice changes from a simple lamellar system of parallel vertical brine layers to a percolation-controlled pore network. In the percolation regime below ϕ₀, pores are shrinking and disconnecting.

(3) ϕ < ϕ_c – the percolation threshold ϕ_c, below which sea ice is impermeable. In the following, a mathematical model will be formulated that quantifies the transitions and the permeability–porosity relationship in dependence on the growth rate of the ice.

Fig. 1 Vertical permeability simulations versus brine porosity with best columnar ice fit from Maus et al.²¹ compared to the relationship obtained by Freitag²⁴ for young ice (red curve), and the lamella model at high porosity (eqn (5) for a typical range in sea ice plate spacing). The regime boundary for ϕ₀ described in the text is chosen tentatively.

3 Fundamental microstructure scales – spacing and thickness of brine layers

Sea ice is henceforth idealized (see eqn (2) above) as an ensemble of parallel vertically oriented plates – ice lamellae. This structure defines the noted plate or brine layer spacing a₀, between which brine layers of thickness d are sandwiched (see Fig. 2). While there are deviations from this structure (e.g. at the ice surface), it is a reasonable simplification of columnar sea ice and has been used in many other studies.¹ Brine porosity ϕ, plate spacing a₀ and brine layer thickness d are simply related as ϕ = d₀/a₀.

Fig. 2 3D (a) and 2D (b) micro-CT-based images of open pores in young sea ice, illustrating the plate (or brine layer) spacing a₀, the width d and the locations where the layer has split/bridged.

The key hypothesis proposed here is that there is a transition from a system of parallel brine lamellae to more complex percolating pore networks that takes place when the thickness of brine layers between the plates becomes less than a critical d₀. The corresponding critical brine porosity ϕ₀ is given as:

This concept of critical microstructure scales has been applied earlier. The first studies were focusing on changes in mechanical properties that are expected when the sea ice brine layers start to neck.^30,31 Other authors^3,5,32 later discussed the relevant length scales d₀ and a₀ in connection with percolation-based permeability models. However, in these earlier studies, detailed observations of 3D pore networks in sea ice were lacking. During recent decades, X-ray computed microtomography (XRT) has emerged as the method of choice to study the 3D microstructure and physics of porous ice and snow media in the environment, including snow,^33,34 polar firn and ice cores,^35,36 and sea ice.^{20,21,37–39} In the following, I show that the recent 3D microstructure data from ref. 21 allow for a consistent quantification of the above basic length scales. Fig. 2 illustrates d₀ and a₀ in 3D and 2D microstructure images from that data.

3.1 Plate spacing versus growth velocity

The fact that most growing sea ice has a lamellar ice–water interface, with concentrated seawater sandwiched between vertical plates, has long been known, with reports of this microstructure dating back up to two centuries.^40–42 The first quantitative investigations of this plate spacing ranged from 0.2 to 0.5 mm for thin sea ice,^27,30,43 via 0.5 to 1 mm for thick Arctic sea ice,^44,45 to 1.3–1.5 mm for a very thick ice island,⁴⁶ suggesting an increasing plate spacing with ice thickness. In later studies, it was shown that the plate spacing decreases with increasing growth velocity.^28,47 Other investigators failed to demonstrate the growth velocity dependence in laboratory experiments.⁴⁸ The present author^29,32,49 has studied the problem in detail and proposed a model for the plate spacing based on morphological stability theory developed by W. W. Mullins and R. F. Sekerka.⁵⁰ In ref. 29, the following relationship between the plate spacing a₀ of sea ice and its growth velocity V was derived:

a₀ = 0.72V^−1/3, (7)

where a₀ is in millimetres and V in cm per day. This equation is valid for sea ice that grows under the presence of solute rejection and convection at the freezing interface, and was validated by a wide range of observations for ice growth rates below ≈15 cm per day.

For the present study, I extended the earlier microstructure analysis that accompanied the above permeability measurements²¹ with plate spacing measurements and an ice growth model. An overview of the field conditions during ice growth is given in Fig. 3, showing air temperature and precipitation from freeze-up in late March to ice core sampling in mid April, when the ice was 3–4 weeks old. The microstructure data analysed in ref. 21 and the present study stem from 15 ice cores of this 35-cm-thick young ice. Each core (of 7.2 cm diameter) was sawed into 3–4-cm-thick segments, which were centrifuged at their in situ temperatures, as well as at lower temperatures. Finally, X-ray microtomographic imaging was performed on samples of 3 cm diameter and thickness, to obtain 3D images with a voxel size of 18 μm. More information was provided in ref. 21.

Fig. 3 (a) Field conditions (air temperature, precipitation) during young ice growth in Advent Bay, Svalbard, likely freeze-up date, sampling dates and thickness, and modelled ice thickness; (b) predicted and measured (on micro-CT images) plate spacing profile (5 ice cores average).

The plate spacing a₀ for the sampled young ice was obtained by model and observation. The modelled plate spacing is based on eqn (7), with growth velocity V obtained with a growth model for young ice.⁵¹ Model inputs are (i) the meteorological observations (air temperature, wind velocity, precipitation, cloudiness, humidity) from Longyearbyen airport (5 km from the field site), (ii) a simple approximation of the bulk ice salinity of young ice,⁵² (iii) starting date and (iv) a guess for the oceanic heat flux. The model than computes the heat budget based on radiative (longwave and shortwave), sensible and latent heat fluxes, the conductive heat flux through the ice and the oceanic heat flux. The model was run with different starting dates and oceanic heat flux within plausible bounds for winter conditions (1 to 10 W m⁻²). The ice thickness 3–4 weeks later was best matched by the starting date shown in Fig. 3a and an oceanic heat flux of 3 W m⁻². The corresponding time series and profile in the growth velocity (not shown) results in the modelled profile of the plate spacing a₀ shown as black circles in Fig. 3b. Other parameter pairs (starting date and oceanic heat flux) and model settings are possible. However, that would change the predicted plate spacing only by a few percent, and have little impact on the general results presented below.

Measurements of plate spacings were made on horizontal microtomographic images. This will overestimate the plate spacing, if the brine layers are vertically inclined. A correction factor cos(α) was applied, based on the measured vertical inclination angle α, which implied maximum corrections of 0.04 mm. Second, at low temperatures, the brine layers and pores become, at our spatial resolution, invisible in parts of a sample. Fig. 4 shows this by comparing 3 horizontal slices that are 1 mm apart. To obtain the plate spacing from such imagery, a₀ was determined from the images with highest porosity. In addition, non-segmented images were used that often allow for detection of thinner air and brine pores than in segmented images. The measured plate spacing profile are shown in Fig. 3b as red square symbols.

Fig. 4 Horizontal sections of columnar sea ice spaced by 1 mm with decreasing porosity from left to right, illustrating the tendency for some pores to disappear, which increases the apparent plate spacing a₀.

Predicted and observed plate spacings agree reasonably, both in terms of vertical trend and several local extrema. On average, the modelled plate spacing is 0.03 mm smaller than measured, and there appears to be a trend of increasing difference towards the surface. As surface temperatures were lower in the experiment, this difference is likely related to the described pore close-off with decreasing temperature. One may also anticipate an imperfect prediction by eqn (7), an overestimate of ice growth velocities due to uncertainties in the freeze-up date and oceanic heat flux. However, most important is that the average plate spacing for this dataset is consistent for model predictions and observations . In the following, a value of will be used as a reference average for the dataset.

3.2 Critical brine layer thickness d₀

The second length scale d₀ marks the condition when brine layers bridge. It is more difficult to estimate. The problem has been first discussed in connection with studies of mechanical properties of sea ice. D. L. Anderson and W. Weeks³⁰ proposed, from the analysis of 2D horizontal thin section photographs, that brine layers would split at 0.07 mm thickness (the thinnest layers observed) and, driven by surface energy minimisation, evolve into circular channels with about two times the diameter. As on the samples in question, a plate spacing a₀ ≈ 0.46 mm was measured, and these authors proposed a critical brine volume of ϕ₀ = d₀/a₀ = 0.07/0.46 = 0.152, below which many mechanical properties of sea ice should change fundamentally, and presented observations of the elastic modulus to support the idea. Assur⁴⁴ investigated the idea in connection with a microstructure-based model of the tensile strength of sea ice. However, the value for d₀, and the surface energy hypothesis, have not been confirmed by later studies, and the number of measurements obtained by Anderson and Weeks³⁰ were rather limited. Weeks and Assur³¹ later estimated d₀ indirectly on the basis of strength tests on ice with known spacing a₀; extrapolating the strength-porosity dependence to zero strength at ϕ₀ = d₀/a₀, they obtained d₀ = 0.112 ± 0.01 mm. As this was much larger than d₀ = 0.07 mm proposed earlier,³⁰ Weeks and Assur noted that “It is hoped that new measurements of d₀…will soon be available”. Microstructure and strength studies performed since then indicate that d₀ most likely lies at the higher end of this range.^32,53

In our recent study of 3D microstructure data,²¹ the author and his colleagues have obtained a more precise estimate of d₀. The relevant data of that study are summarised in Fig. 5. Two length scales were obtained over a wide porosity regime – the median diameter d (or shortest dimension, for non-circular pores) of open brine pores and the throat diameter d_t or thinnest part of a pore. For both length scales, a robust power law dependency on the total brine volume ϕ was found, given as:

d_t = 0.389ϕ^0.46 mm (8)

for the pore throat, and:

d = 0.417ϕ^0.34 mm (9)

for the median open pore diameter. The physical significance of these relationships for the present problem becomes clear when evaluating the equations at the percolation threshold ϕ_c = 0.24 obtained in that study. In Fig. 5, this condition is given by the point where the threshold (vertical line) and the length-scale-porosity relationship cross and thus:

d₀ = d(ϕ_c = 0.024) ≈ 0.117 ± 0.008 mm. (10)

The corresponding throat size at the threshold porosity ϕ_c = 0.024 is d_t0 ≈ 0.070 ± 0.004 mm. This result is interpreted in the way that the pores start disconnecting towards an impermeable pore space when reaching a thickness of d₀ ≈ 0.12 mm, while the characteristic throat size d_t0 is somewhat smaller.

Fig. 5 (a) Median throat width d_t versus brine volume of young ice from ref. 21; (b) median pore diameter d. The critical d_t0 and d₀ values are obtained at the intersection of the relationships with the percolation threshold ϕ_c.

These values, obtained from 3D microstructure imagery, are the first robust results for the critical length scales for ice bridging across brine layers. They improve on earlier, mostly indirect estimates. Interestingly, the thinnest observable thickness of brine layers reported by Anderson and Weeks⁵⁴ was 0.07 mm, and compares to the characteristic throat diameter d_t0 for pore disconnection found by us. The most important scale for the present theory is the average brine layer thickness at that stage, for which d₀ = 0.12 mm is considered as the best empirical estimate.

4 Percolation theory – directed versus isotropic

“Percolation” (from latin “percolate” = to filter) refers to the movement of fluids through a porous material like soil, rocks and porous ice media as firn, snow and sea ice. During the past half century “percolation theory”, first described by Broadbent and Hammersley⁵⁵ has been developed in mathematics and physics. In addition to the structure and properties of pore networks, it deals with broader connectivity problems in other systems, such as forest fires and pandemics.^55–57 Key aspects of percolation theory are that it (i) is a probabilistic model of the system connectivity and (ii) exhibits critical behavior at a percolation threshold. The point of criticality is often termed phase transition.

The first sea ice studies on this problem concluded that the sea ice pore pace undergoes a percolation phase transition and becomes impermeable at a brine volume fraction of 0.05.^3,19,20,26 However, more recent work based on 3D X-ray microtomographic imaging of sea ice has found a threshold of 0.02–0.03 percent for young sea ice.^21,58 While the former studies^3,5,19,20,26 were based on isotropic percolation, the latter work by the present authors^21,58 indicated that directionally growing sea ice should rather be described in terms of directed percolation theory. In directed percolation, fluid is restricted to flow in one spatial direction, and the critical behavior differs from isotropic systems;^59,60 e.g., the important differences with respect to sea ice are summarised as follows.

4.1 Percolation properties and critical exponents

An essential aspect of percolation theory is that many properties P of a medium, like connectivity and permeability, near the percolation threshold may be described by a power law of the form:

P ∼ (ϕ − ϕ_c)^q, (11)

where the critical exponent q only depends on the dimension of the system, and not the microscopic details. Many exponents have been determined numerically for both isotropic and directed percolation. Some of the currently accepted 3D values are given in Table 1. These are the exponents β for the strength of the connected cluster (for sea ice: the connected porosity), μ for the electrical conductivity of the pores space, and t for the permeability. While some exponents are known with higher precision, the table suffices for this study and compares those q values that have been determined by the present analysis of sea ice.

Table 1 Selected critical exponents for d = 3 isotropic and d = 3 + 1 directed percolation systems

Property	Exponent	Isotropic	Directed	This study
Cluster strength	β	0.41 (ref. 61)	0.82 (ref. 62)	0.83 ± 0.03
Conductivity	μ	2.0 (ref. 56)	1.7 (ref. 63)	1.8 ± 0.1
Permeability	t	> 2.0 (ref. 64 and 65)	>1.7	2.55 ± 0.25

---

The critical exponents for sea ice shown in Table 1 have been obtained as follows. As described in ref. 21, their approach was to centrifuge brine from young sea ice core segments at different temperatures. The result obtained from centrifuging 15 young ice cores is shown in Fig. 6a. The fitted relationship between centrifuged (effective, connected) porosity ϕ_eff and total brine porosity ϕ may be written in the form

ϕ_eff = c_ϕ(ϕ − ϕ_c)^β, (12)

where ϕ_c is the critical porosity at which all pores are disconnected. By linear regression the best estimates of ϕ_c and the exponent β were obtained as ϕ_c = 0.024 (95% confidence bounds [0.20, 0.29]) and β = 0.832 (95% confidence bounds [0.803, 0.861]), with c_ϕ = 0.569. In another recent study of young ice by Salomon et al.,⁶⁶ assuming the same ϕ_c, a slightly higher exponent has been obtained (β = 0.869 with confidence bounds [0.803, 0.936]). Next, the critical exponent of 2.55 for the permeability obtained by Maus et al.²¹ has already been mentioned in the introduction, eqn (3). Here I further present numerical simulations of electrical conductivity with ref. 67 on the same CT images on which the permeability was determined (Fig. 6b). The best percolation-based fit to the electrical conductivity data is the equation

σ/σ₀ = c_σ(ϕ−ϕ_c)^μ, (13)

where the exponent μ = 1.8 ± 0.1 and c_σ = 1.194.

Fig. 6 (a) Centrifuged (effective, open) brine porosity of young ice with optimum percolation fit from ref. 21; (b) normalised electrical conductivity with optimum percolation fit against (ϕ − ϕ_c) (new result).

An essential result in Table 1 is that the deduced critical exponent β = 0.83 ± 0.03 for sea ice is very close to the presently accepted β ≈ 0.82 for 3(+1)D directed percolation,^59,60 and far from β ≈ 0.41 in 3D isotropic percolation.^56,57 Also for the electrical conductivity the sea ice exponent μ = 1.8 ± 0.1 is slightly closer to the directed percolation exponent of 1.7 from ref. 63, when compared to the isotropic μ = 2.0.⁵⁶ These results suggest that transport through the pore space of columnar sea ice belongs to the universality class of directed percolation. Note that the critical exponent t for the permeability is not universal, yet depends on the details of the pore space evolution. Its lower bound is the conductivity exponent μ,^64,65 e.g., the present t = 2.55 should thus be viewed as an empirical value for young sea ice.

Note that there is another important exponent ν that describes the behaviour of the correlation length, which may be understood as the average distance of sites (or pores) in a network.⁵⁶ In 3D directed percolation it is direction-dependent and has been determined as ν_∥ ≈ 1.11 in the percolation direction and ν_⊥ ≈ 0.58 normal to it.⁶⁰ We have omitted this property as the sea ice samples analysed here have been too small to retrieve the correlation length and its critical exponent.

4.2 Percolation and porosity threshold

While the critical exponents μ and β are considered universal and micro-structure independent, this is not the case for the percolation threshold ϕ_c, at which the connectivity of a network breaks down. As an example one may consider the percolation of a system of thin rods, for which it is known that the critical solid fraction, at which a random assemblage of these rods is interconnected, decreases with increasing aspect ratio. However, there are some interesting general bounds that have been reported for classical lattices. Here the percolation threshold depends on the coordination number Z, and if one defines it in terms of sites or bonds, e.g., on a 3D simple cubic lattice with Z = 6, the site and bond percolation thresholds are ≈0.31 and ≈0.25, respectively, being smaller (≈0.25 and 0.18) for a body-centered cubic lattice with Z = 8 (see Table 2). However, half a century ago, Scher and Zallen⁶⁸ detected an important aspect for both 2- and 3-dimensional lattices: the higher the coordination number, as well as packing or filling factor F of a lattice, the lower the critical percolation probability P_c, such that their product f_c = FP_c is approximately constant. The property f_c resembles the critical density or filling fraction of a network. For selected 3D lattices a rather constant critical fraction f_c in the range 0.15 to 0.16 was found. In Table 2 we compare this result to the corresponding thresholds for directed percolation that have become available during recent years. An interesting fact is that one can make a similar conjecture for directed percolation, where f_c ≈ 0.11 appears to be the critical density. Furthermore, it appears that the ratio of isotropic to directed threshold densities is close to 1.5, i.e., that the directed percolation threshold is 2/3 of the isotropic value. This aspect will be important when comparing results for columnar and granular sea ice below.

Table 2 Site percolation thresholds for the simple cubic (SC), body-centered cubic (BCC) and diamond (Ice Ic) lattices in three dimensions for isotropic and directed percolation. Z is the coordination number or number of bonds per site, and F is the filling factor. The critical filling fraction f_c is the product of F and the critical percolation probability P_c

Lattice	Z	F	Pc, iso	Pc, dir	fc, iso	fc, dir	fci/fcd
BCC	8	0.680 (ref. 68)	0.246 (ref. 69)	0.161 (ref. 70)	0.167	0.110	1.52
SC	6	0.524 (ref. 68)	0.312 (ref. 69)	0.208 (ref. 71)	0.163	0.109	1.50
Ice (Ic)	4	0.340 (ref. 68)	0.430 (ref. 69)	0.303 (ref. 72)	0.146	0.103	1.42

---

5 Synthesis of a permeability model based on directed percolation

Combining now the general percolation scalings with the sea ice microstructure results from the previous sections one can write

ϕ_c = f_cϕ₀ (14)

for the upper and lower thresholds of the percolation regime. ϕ₀ marks the onset of directed percolation when ice bridges start forming within the vertical brine layers, while the percolation threshold ϕ_c is given by the product of ϕ₀ and critical percolation density f_c. In terms of microstructure scales this reads

These upper and lower thresholds separate the regimes described in Section 3 that now can be defined in terms of the length scales d₀ and a₀, as well as f_c:

(1) ϕ > ϕ₀ – the brine layer regime: the permeability of an ensemble of shrinking brine layers is given by eqn (5) and thus is proportional to a₀² and ϕ³. With a₀ being related to the ice growth velocity through eqn (7), the permeability equations will be growth-velocity dependent.

(2) ϕ_c < ϕ < ϕ₀ – directed percolation regime: at a critical thickness of d₀ ≈ 0.12 mm (and porosity d₀/a₀) ice bridges form across the brine layers and a more complex network pores with variable diameter forms. With decreasing porosity (due to cooling of the ice) this percolation process continues to close more and more pores. The permeability in this regime is given by eqn (3) as K ≈ (ϕ − ϕ_c)^t, with empirically determined exponent t = 2.55.

(3) ϕ_c – directed percolation threshold: below ϕ_c sea ice is vertically impermeable. This happens at a critical pore fraction of f_c ≈ 0.11 in the original brine layers, corresponding to the porosity threshold f_cd₀/a₀ given by eqn (14).

The permeability for these porosity regimes is explicitly formulated as a function of a₀, d₀, f_c, and the percolation critical exponent t for the permeability:

K = c_k(ϕ − f_cd₀/a₀)^t [f_cd₀/a₀ < ϕ < d₀/a₀] (17)

K = 0 [ϕ < f_cd₀/a₀] (18)

The parameter c_k follows from matching K from eqn (16) and (17) at the upper critical porosity ϕ = d₀/a₀:

A first test of the consistency of this model with observations is made by comparing c_k from eqn (19) to the factor from the numerical simulations, given as 1.48 × 10⁻⁸ m² in eqn (3). Inserting the best estimates of d₀ = 0.12 mm, a modelled plate spacing , t = 2.55 and f_c = 0.11 into eqn (19) one obtains c_k = 1.66 × 10⁻⁸ m², while using the observed gives c_k = 1.81 × 10⁻⁸ m². The standard settings thus give predictions that are 12–22% higher than the permeability simulations. The predicted percolation threshold based on eqn (14) is 0.023 or 0.024, depending if one uses the measured or modelled average plate spacing of , respectively. This result also agrees remarkably with the value ϕ_c = 0.024 ± 0.04 from the best centrifuge-based fit (eqn (12)). Hence, the microstructure-based model agrees quantitatively with observations of macroscopic system properties (permeability, connected porosity, percolation threshold), showing the consistency of the macroscopic and microscopic data.

6 Discussion

6.1 Permeability versus growth rate and plate spacing

A central result of the present model is that the percolation threshold (eqn (15)) may be written aswhere the numerator 0.0123 = f_cd₀ mm is the product of the critical bridging length scale d₀ and a network directed percolation threshold f_c, assuming that these are constant. For the product f_cd₀, one can, based on eqn (10) and the variation of f_c in Table 2, estimate an error of 0.0009 mm (when both d₀ and a₀ are given in mm). The percolation threshold is thus inversely proportional to the plate spacing a₀. While such a conjecture has been proposed in previous studies,^3,21,32 the present analysis puts it on a concise theoretical and observational ground. With the growth velocity dependence of a₀ (eqn (7)), one may write this in the form:

ϕ_c ≈ 0.0183V^1/3, (21)

when V is given in cm per day. The percolation threshold thus increases with the cubic root of the ice growth rate. For V in the range 0.5–1 cm per day, typical for thicker first-year ice,^1,73 ϕ_c becomes 0.015–0.018. For rapidly growing young ice, say 10 cm per day, the threshold will be 0.039. The expected range of ϕ_c for naturally growing sea ice is thus likely 0.01–0.04.

The model was further tested by sorting the permeability data from Fig. 1 into different groups of plate spacing and growth velocity (according to Fig. 3). The result is presented in Fig. 7 for the modelled plate spacing (being very similar for the measured plate spacing). In the figure three curves are drawn. The black intermediate curve shows the model result (based on eqn (16) to (19)) for an average plate spacing of a₀ = 0.55 mm. This curve is almost indistinguishable from the best fit to all data that is also shown as a hatched curve. The two other curves show the model predictions for the upper (5 cm per day) and lower (1 cm per day) bounds during the growth of the ice. For comparison, two subgroups of the data have been selected that correspond to a high (>3.5 cm per day) and low (<1.5 cm per day) regime, and are highlighted with red and green dots. While the scatter is large, the distribution clearly supports the hypothesis of a growth-velocity-dependent permeability. In particular, the low growth-velocity results may be much better understood by a model with a lower permeability threshold (upper red curve). Note that Maus et al.²¹ did not include these values in the least squares fit of Fig. 1, arguing that they are from a mixed regime of 0.024 < ϕ < 0.031 where both zero and non-zero permeabilities were observed. Fig. 7 now explains this in terms of an ice growth rate difference.

Fig. 7 Vertical permeability simulations of columnar ice versus model results (eqn (16) to (19)). The high and low growth rate simulations are depicted with green and red dots, respectively, and compared to the predictions for upper (green curve) and lower (red curve) ice growth rates in the model. The best fit and the model prediction for the average plate spacing are also shown very close to each other.

The analysis was extended by another parametric fit. It is seen that each point in Fig. 7 can be interpreted in terms of a unique growth velocity (and/or plate spacing) and percolation threshold. One can thus, based on the model curves, compute ϕ_c for each simulation point, and compare it with the observed plate spacing a₀ to derive the corresponding critical filling fraction f_c of the brine layers through eqn (15). This procedure has been performed for the modelled and observed plate spacings from Fig. 3, and by dividing the data into low- and high-growth-velocity subgroups (simply divided by the average plate spacing). The results are shown in Fig. 8, where the inverse of the plate spacing is plotted against the average percolation threshold obtained from the parametric model fit. For the whole dataset, the average f_c results are shown as full and open black circles, being 0.112 ± 0.007 for the modelled a₀ and 0.119 ± 0.007 for the observed a₀. These values are close to f_c ≈ 0.11 obtained for directed percolation theory (Table 2). While not surprising, as the estimates of ϕ_c from theory and observations agreed well, this comparison reconfirms that the model settings and data are consistent.

Fig. 8 Inverse measured and observed plate spacings a₀⁻¹ versus the average percolation threshold ϕ_c obtained by drawing model curves through each permeability simulation point for all data and for subgroups of high and low growth rate (low and high a₀). The slope of the lines is d₀f_c.

Fig. 8 also shows the percolation thresholds for the low and high growth velocity subgroups. For perfect agreement these should fall on the stippled and full lines for the observed and modelled plate spacings. Due to the reduced number of data points for the subgroups, the standard deviation is large, and the confidence of these results is not very high. However, the consistence is very reasonable, supporting the hypothesis that the percolation threshold is inversely proportional to the plate spacing a₀.

6.2 Granular surface ice

The above data analysis and fitting procedures have focused on columnar sea ice, omitting the results for granular surface ice (not shown so far). In the earlier paper²¹ it was shown that the results for granular surface ice samples indicate a higher percolation threshold. The database was, however, too limited for a quantitative analysis. According to the present analysis, the dependence of permeability on the growth rate of sea ice might serve as an explanation, as ice is initially growing faster (leading to the low plate spacing shown in Fig. 3).

In Fig. 9 only the granular surface ice data (typically the upper 3 centimeters) are compared with the model predictions. The red and green curves span, as shown in Fig. 7 for the columnar ice data, the growth velocity regime of the ice. Although the scatter in this data is large, it is seen that the majority of data points fall below this regime, and that there are also two observations with zero permeability above ϕ ≈ 0.05. As a possible explanation one may consider that the very surface ice is growing faster than in the present ice growth model, e.g., assuming a growth rate of 10 cm per day would raise the percolation threshold to ϕ_c = 0.04. However, there is another aspect that needs to be considered. Surface sea ice is not only growing faster, it also is very often characterised by a random orientation of crystals, and termed granular. For such ice one would expect that the percolation is isotropic. An ad hoc approach to account for this in the present model involves two steps:

(1) The directed f_c ≈ 0.11 is replaced by the isotropic f_c ≈ 0.16 (Table 2) in eqn (17) to (19).

(2) As the high porosity state of flow through granular ice, in contrast to the columnar lamella model, is tortuous, a factor T⁻² is introduced in the permeability equations (eqn (16) and (17)). Results for the permeability of granular ice and the present samples (not shown) indicate that T⁻² ≈ 1/2 is a reasonable value.^21,39

Fig. 9 Vertical permeability simulations of granular surface ice versus model results (eqn (16) to (19)) for high and low growth rates as in Fig. 7, plus a modified high-growth-velocity prediction assuming isotropic percolation (blue dashed curve), explained in the text.

The results of such a granular ice modification are shown in Fig. 9 as a blue dashed curve. Assuming the highest growth velocity at the surface, this predicts a percolation threshold of ϕ_c ≈ 0.046. While the number of permeability simulations is limited, the higher threshold is consistent with the results. More data would be needed for a validation. It also should be noted that the structure of granular ice may differ from the assumed simple plate assemblage with random orientation, and that its critical exponent t may also be different from the columnar value. However, the basic idea introduced here is that granular sea ice has a 50% higher percolation threshold than columnar sea ice. It will be further discussed in connection with previous studies.

6.3 Revisiting previous studies

How do the results compare to earlier studies that had largely agreed on a percolation threshold of ≈0.05?^3,19,20,26 In a recent study²¹ the author and his colleagues qualitatively discussed the following factors as possible reasons: (i) differences in ice growth velocity, (ii) the difference between granular and columnar ice, (iv) aspects of full-thickness permeability and aging of sea ice, (iv) spatial resolution (when analysing CT images²⁰), and (v) limitations in indirect methods, like desalination measurements of sea ice. The present permeability model does now allow a more quantitative assessment of these factors.

6.3.1 Surface ice and full-depth permeability. To my knowledge, Fig. 9 presents the only results so far that contrast the permeability of columnar and granular sea ice in a concise quantitative manner. However, several authors have investigated the full-depth permeability of sea ice, and for these studies one would expect that the granular surface layer constrains the percolation threshold. In the first account on sea ice percolation, K. Golden and his colleagues¹⁹ referred to several observations that all involved permeation through the full ice depth. One observation from the East Antarctic Sea Ice Zone involved the flooding of 60-cm-thick ice from below that happened when the ice warmed due to increasing air temperature. Before the flooding event the authors reported a brine volume of 0.04 to 0.05 at the ice surface, which during flooding increased rapidly to 0.1. However, from the data presented it is difficult to constrain at which porosity between 0.04 and 0.1 the ice surface became permeable. A second observation reported by these authors refers to the cutoff of algae growth in thin granular Weddell Sea ice, interpreted as a transition to an impermeable state. It was reported at a temperature of ≈−4 °C and a salinity of ≈5, which corresponds to a brine volume fraction of 0.062. However, also here the salinity and temperature were only approximately known and one may anticipate that the uncertainty of the estimate is at least ± 0.01.

The third dataset of full depth permeability discussed by Golden et al.¹⁹ is the experimental data from Ono and Kasai,⁷⁴ who measured the upward and downward permeability through 6-cm-thick laboratory grown sea ice at different temperatures. Golden et al.¹⁹ also associated these experimental results with their proposed percolation threshold (≈0.05). However, the author and his colleagues²¹ pointed out that this (very thin) ice must have had a high salinity and, if at all, the data would indicate a percolation transition at a brine volume above 0.1. On the basis of the present model it is now possible to revise the interpretation of these data. Fig. 10 compares the data from Ono and Kasai⁷⁴ to the present permeability simulations and model predictions. The low- and high-growth-rate model results are again shown with red and green curves, respectively. However, the velocity bounds have been extended with respect to Fig. 7. The lower bound reflects now a plate spacing of 1.0 mm (V ≈ 0.4 cm per day), while the upper bound corresponds to the average growth rate in the experiments from ref. 74, which was V ≈ 8.6 cm per day (giving a₀ ≈ 0.35 mm). The green and red curves thus span the typical growth velocity regime from rapidly growing lead ice to thick ice, as well as the observed range in plate spacing.²⁹ It is clear that the experimental results from ref. 74 fall several orders of magnitude below the columnar sea ice model predictions. However, as the experiments from ref. 74 refer to the full-depth permeability, they very likely imply the flow through a granular surface ice regime. To account for this the modifications as described in Section 6.2 were made, and the high-growth-velocity solution was also drawn for this granular ice model, shown as the hatched blue curve denoted as ‘isotropic’. This curve now shows much better agreement with the upward permeability tests from ref. 74. The modelled (granular) percolation threshold in that case is ϕ_c ≈ 0.055. One may further argue that the very surface ice likely has grown faster than the average of V ≈ 8.6 cm per day reported by Ono and Kasai;⁷⁴ e.g., for a 50% larger growth rate (≈13 cm per day) the model would closely fit the upward permeability data, with a porosity threshold ϕ_c ≈ 0.055 (not shown).

Fig. 10 Vertical permeability simulations versus model results (eqn (16) to (19)). Compared to Fig. 7 the model bounds for high and low growth rate are chosen somewhat higher and lower, to resemble the range in naturally growing sea ice (with percolation threshold range 0.013 < ϕ_c < 0.038). The upper growth rate is set to V = 8.6 cm per day for comparison to Ono and Kasai’s⁷⁴ experiments (blue symbols). A modified high-growth-rate prediction, assuming isotropic percolation, is shown as blue dashed curve, and the corresponding ϕ_c = 0.055 is also highlighted.

Note that Golden et al.¹⁹ only showed the much smaller downward permeability observations from ref. 74, and that their hypothesis of a percolation threshold of 0.05 for these data implies a salinity of ≈5. The data points in Fig. 10 are based on a two-times higher salinity that appears much more likely for this thin, young ice.²¹ The comparison to the present model indicates that the downward permeability rather reflects a porosity threshold of 0.1 or higher. It is not clear why the downward permeability should be so much smaller than the upward permeability (and the threshold so large). The author and his colleagues²¹ suggested that, if the downward experiments were performed after the upward experiments, the ice salinity could have been reduced (which would shift the data points to lower porosity to the left). Another explanation could be that pouring cold brine on the surface, as done by Ono and Kasai⁷⁴ in downward tests, led to additional crystal growth affecting the permeability.

In summary, when modifying the present model for granular surface ice (using an isotropic instead of directed percolation approach), one would expect a 50% higher percolation threshold than for columnar sea ice. For an initial growth rate of 5 to 10 cm per day for surface ice, this implies a percolation threshold in the range 0.04 < ϕ_c < 0.06. A review of permeability measurements involving flow through the ice surface appears to be consistent with such a hypothesis.

6.3.2 Effects of sample dimensions and centrifuging. As shown in Fig. 1, the permeability relationship from Freitag²⁴ is close to the field values from our study,²¹ with Freitag’s relationship being ≈30% below the present best fit. Freitag’s equation was obtained through experiments with young sea ice grown in a tank to a thickness of ≈16 cm. The average growth velocity of 2.3 cm per day, as well as the velocity range of 1 to 5 cm per day, were close to the field conditions in our study.²¹ Both investigations were based on centrifuging of ice samples for which three potential biases have been pointed out:^21,24 (i) vertical channels are connected for short (in the vertical) samples, but closed in taller samples; (ii) inclination of the channels makes them leave a sample laterally, not contributing to vertical permeability; (iii) when cooled after centrifugation, the brine in dead-end pores, that could not be centrifuged out, may be expelled into open pores and lock them. While (i) leads to an overestimation, (ii) and (iii) tend to underestimate the permeability. Freitag²⁴ investigated aspect (i) by comparing the permeability of core segments of different lengths 1.5, 3.5, 7 and 16 cm (full core), and proposed that a length of 6 cm was sufficient to reduce this effect. However, much of the vertical permeability resolution was lost in this way, and the results were dominated by the minimum permeability in core segments. In my opinion, effect (i) is not a primary error source for young ice when 3–4 cm segments are studied. The author and his colleagues²¹ rather suggested (ii) to be a major issue. To estimate this, consider a cylinder with a quadratic cross section, and assume that the channels are randomly oriented by the angle α against the vertical. The fraction of channels that will not reach from top to bottom and run out laterally may then be estimated as tan(α)ε, where ε is the ratio of sample height to length. For the present ice we found α to decrease from 20–25° at the surface to less than 10° at the ice–water interface, which is quantitatively similar to the c-axis alignment described by others.^1,75 Using 10 < α < 20° as a characteristic inclination angle, and ε = 1/4 for our simulations,²¹ we get tan(α)ε in the range 0.04 to 0.09 as the underestimate of the permeability. For the experiments by Freitag,²⁴ with larger ε = 2/3, the underestimate would be 0.12 to 0.24. The effect thus may also explain some difference between the numerical simulations and the experiments by Freitag.²⁴ For (iii) where the closure of open pores is due to brine expulsion from dead-end pores, there are currently no quantitative estimates. I suspect that it may play some role near the percolation threshold, when the connected porosity is small. In a study with low (10×g) centrifuge acceleration, considerable amounts of brine remained in the samples.⁷⁶

6.3.3 Spatial resolution. Only a few other authors have investigated the permeability or percolation threshold of sea ice by means of micro-CT image analysis. Pringle et al.²⁰ studied columnar sea ice grown in a tank by seeding with a layer of natural sea ice on top. This procedure allowed a contrast agent to be added to the water from which the ice was grown, which provided reasonable micro-CT data quality without centrifuging. The authors used finite size scaling methods of isotropic percolation theory⁵⁶ to estimate the percolation threshold. Different approaches gave results for ϕ_c in the range from 0.039 ± 0.003 to 0.067 ± 0.007. From two of their published images (Fig. 1a and b in ref. 20), one can estimate a plate spacing a₀ of 0.5 to 0.6 mm, and thus would expect a similar percolation threshold as in the present study. Why are the values reported a factor of 2 larger? As the fitting method by Pringle et al.²⁰ is based on β = 0.41 for isotropic percolation (and a correlation length exponent ν = 0.88), I have recalculated their results with directed percolation exponents β = 0.82 and ν_∥ = 1.11. This indeed changes the range to slightly lower values for ϕ_c (0.037 to 0.056), yet they are still at least 50% larger than the present data and predictions. As discussed in our related study,²¹ another aspect is important. Considering that the detectability of pores is limited by the Nyquist criterion (2 times the voxel size), in the study by Pringle et al.²⁰ with voxel size 41.5 μm this detection limit was 83 μm. This is larger than the critical pore necking scale shown in Fig. 5, d_t0 ≈ 70 ± 0.04 micrometer that we found on the basis of 2.3-times-better spatial resolution (voxel size 18 μm). It may then be anticipated that the percolation threshold found by Pringle et al.²⁰ was limited by spatial resolution. As an example one may insert 83 μm on the left-hand side of eqn (9) to estimate the porosity at which the ice from ref. 20 would appear impermeable – which gives a value of ϕ_c = 0.035. When the same exercise is done for the maximum diameter of a flow path (not shown, see Fig. 11d in ref. 21) a limit of 0.042 is obtained. Spatial resolution thus may reasonably explain the discrepancy between the present (and ref. 21) results and those from Pringle et al.²⁰

Among the few other micro-CT studies, Salomon et al.⁶⁶ reported, for a voxel size of 25 μm, vertical connectivity down to a porosity of 0.05, yet the number of samples below that value was limited (at 0.03 two were permeable, one not). The authors of ref.39 have analysed granular and columnar Arctic sea ice samples with the same spatial resolution as in ref. 20, a voxel size of 41.5 μm, and reported a percolation threshold of ϕ_c ≈ 0.08 for columnar ice. As this ice was from intermediate depths (0.3–1 m) of 1.4-m-thick sea ice, this is a large value compared to the present predictions. Other factors like raw data quality, noise and filtering, as well as aspects of centrifuging and sample dimensions (see above) may have influenced the results. In that context another study of young sea ice⁷⁶ is interesting: here the voxel size was 11.8 μm and the percolation threshold of columnar samples appeared to be close to 0.05. However, in that study considerable amounts of brine could not be centrifuged out, which was likely due to a low (10×g) centrifuge acceleration.

6.3.4 Desalination as an indirect permeability indicator. Cox and Weeks² were the first to propose a critical brine porosity ϕ_c = 0.05 when deriving an empirical relation for the desalination of sea ice. The data has been later analysed by several authors,^3,4,32 e.g. Petrich et al.³ proposed, based on their analysis, a porosity threshold of ϕ_c = 0.054 below which desalination ceases. They also proposed an equation for the percolation threshold similar to 15, yet with different values for d₀ and a₀, and by employing isotropic percolation. While that study appeared to support a threshold of 0.05, a closer look is needed. First, sea ice desalination is known to depend on other factors than the permeability, like the brine salinity gradient and the critical Rayleigh number.^5–8 Second, desalination involves convection patterns that are constrained by horizontal permeability, and thus the horizontal percolation threshold. The latter appears to be 2–3 times higher (depending on horizontal direction) than in the vertical.²⁰ While more observations are needed, one can argue as follows. As the data from Cox and Weeks² are for young ice with similar growth velocities as in our work, the present model would suggest ϕ_c of 0.02–0.03. Assuming a ratio of 2 between the horizontal and vertical percolation threshold would give a horizontal threshold range of 0.04–0.06. This in turn agrees reasonably with the noted desalination threshold.

6.3.5 Electrical conductivity. For the electrical conductivity a similar model may be formulated based on the fits and experimental data in Fig. 6. While such an analysis is beyond the scope of the present paper, one can say that electrical conductivity observations are not in conflict with the present theory. The resistivity measurements in ref. 77 indicated a change in the mode of the electrical conduction mechanism below a brine porosity of 0.03. Also, more recent observations have shown that electrical conduction does not show anomalous behavior above a porosity of 0.04–0.05.⁷⁸ However, there is little data on electrical conduction at lower porosity.

6.4 Relevance for sea ice modelling

The essential results of the present work for sea ice modelling summarised in Fig. 10 are highlighted as follows:

(1) The percolation threshold of columnar sea ice depends on its growth rate and can be expected to vary in the range 0.01 < ϕ_c < 0.04 for typical natural growth conditions.

(2) The corresponding permeability, when corrected for brine porosity, may vary by almost 2 orders of magnitude.

(3) The permeability can be parametrised in terms of brine volume and ice growth velocity alone.

(4) The estimate of a 50% larger percolation threshold for granular surface ice (compared to columnar ice grown at the same rate) is another important conjecture supported by observations.

The present work provides the relationship of permeability and brine volume ϕ for different growth conditions. The brine volume dependence on ice temperature and salinity is approximately ϕ ≈ c_mS_i/T_i, where T_i is given in °C and c_m is typically in the range −0.05 to −0.07 °C/ppt for natural sea ice conditions (for more exact formulations see ref. 44 and 79). The sea ice salinity S_i is thus of particular importance to obtain the permeability. Starting with the work by Cox and Weeks² there have been considerable efforts to numerically predict the sea ice salinity evolution based on ‘mushy layer’ theories of gravity-driven brine drainage.^5–8,80,81 All these models require as input a permeability–porosity function that in most cases has been taken from the discussed study by Freitag.²⁴ The present model improves this considerably and makes it possible to account for microstructure and growth rate effects, as well as physically consistent percolation thresholds. However, some caution is in order: the mentioned large-scale models have so far only used the vertical permeability in their model frameworks, while horizontal permeability and percolation thresholds are likely to play a role. Hence one needs more observations and modelling effort on the anisotropy of permeability to properly predict salinity. The present author is currently working on a salinity prediction model based on the microstructure scales d₀ and a₀, in line with the present work and some earlier approaches.^5,32,82

The permeability equations (eqn (16) to (18)) may, using the relationship in eqn (7) between a₀ and growth velocity V, be written as:

K = 0.0432V^−2/3ϕ³ [ϕV^−1/3 > 0.167], (22)

K = c_k(ϕ−0.0183V^1/3)^2.55 [0.0183 < ϕV^−1/3 < 0.167], (23)

K = 0 [ϕV^−1/3 < 0.0183], (24)

where

c_k = 2.60 × 10⁻⁸ V^(−1.55/3) m² (25)

was obtained from eqn (19) by inserting the best estimates for d₀, f₀ and t. Note that V is given in cm per day, which is a convenient unit for sea ice growth, while the permeability K is given in m². In this form the permeability is a function of brine porosity ϕ and growth velocity V that now effectively parametrise the microstructure information. While current sea ice models do not compute or track microstructure characteristics, the ice growth velocity V is a property that is routinely computed and could be tracked to implement the present permeability model. Eqn (22) to (25) may then improve the representation of sea ice processes that depend on permeability and precise porosity thresholds in models. The approach likely also bears potential to improve our understanding of mechanical properties, for which the question of a critical d₀ was once raised.^30,44 A recent study by the author on the sea ice tensile strength supports this.⁵³ With fluid transport and mechanical properties controlled by the same porosity thresholds also their interaction may be better understood, e.g., the permeability of sea ice under a compressive load.⁸³

6.5 Ocean–ice–atmosphere fluxes of gas and matter

The vertical permeability of sea ice is also essential for understanding sea ice biogeochemistry and air–ice gas fluxes.^11–14 Sea ice plays a different role than open water or cracks in the ice. Due to its low temperature it contains solutes and chemical species at higher concentration, and also allows for bio-geochemical processes in its pore space that would not take place in the oceanic surface layer. This implies interactions between sea ice and atmospheric chemistry, for which the sea ice permeability is important.

An example for the effect of sea ice permeability and chemistry on gas exchange with the ocean and the atmosphere is the CO₂ budget.^84,85 During wintertime and sea ice growth it is the near bottom permeability of sea ice that controls gravity-driven brine convection, and a transport of CO₂ to deeper ocean layers.⁸⁶ When the ice melts, the surface is undersaturated with CO₂, leading to enhanced CO₂ uptake by the ocean. This process is complicated by the carbonate chemistry and thermodynamics of cooling sea ice, where calcium may precipitate as calcite or as ikaite. This implies differences in the alkalinity and buffering capacity of brine, and in the temperature (and permeability) of calcium precipitation. While brine convection transfers the CO₂ to deeper ocean layers, the carbonate precipitates may then remain in the ice, with consequences for the carbon budget.^17,18,85 The permeability may play a critical role in this puzzle.

A related impact on atmospheric chemistry is the enhanced bromine release by a reduced buffering capacity of sea ice brine due to calcite precipitation,¹⁶ with bromine acting as a catalyst for the destruction of stratospheric ozone.^15,87 The process depends on the ability of liquid brine to migrate to the surface, and hence on ice microstructure and permeability. Others⁸⁸ have discussed sea ice as the source of iodine in the atmosphere, with production by micro-algae in the ice linked to brine transport to the surface. Also this process will depend on the detailed permeability and connectivity of the sea ice pore space, and the details of sea ice thermodynamics and growth conditions.

The ocean–ice–atmosphere exchange of dimethyl sulfide (DMS), as a major source of polar atmospheric aerosols, was recently studied by Gourdal et al.⁸⁹ The fluxes were linked to algae growth (during spring and summer) and DMS production at the sea ice bottom in connection with brine circulation and gas diffusion in the ice (molecular diffusion of gas dissolved in brine as well as upward movement of gas bubbles). Also the exchange of methane between sea ice and the atmosphere during the winter season has been studied.⁹⁰ In these and many other studies on gas fluxes (ref. 11 and 13, for example), the permeability has been suggested as an important property. However, most investigators have interpreted their results relying on a constant percolation threshold of ϕ_c ≈ 0.05 and the permeability equation (eqn (2)) from Freitag.^13,89,90 The present model for the permeability and its threshold provides the basis for a revised analysis of such studies.

Two processes are of particular relevance for the noted air–ice gas fluxes. One is the transport of brine to the sea ice surface and the overlying snow layer, triggering air–ice exchange processes as mentioned for the halogen chemistry.^15,87,88 This transport is still poorly understood. During wintertime it may relate to internal freezing and upward expulsion of brine towards the surface. While the colder the ice, the more brine will be expelled, colder ice will also approach the permeability threshold, which in turn depends on growth conditions. This interaction may lead to complex bounds in the brine fluxes. However, what role the micro-scales d₀ and a₀ play at the ice surface, and how permeability evolves there, needs to be further studied. Not only is the initial ice growth not columnar, but the microstructure is also affected by other factors like large temperature fluctuations, drainage processes in the freeboard, interaction with snow and weathering. The second is the transport of gas by rising gas bubbles that form in warming ice during the melt season.^14,90 For such bubble movement in old ice, the percolation thresholds may differ from those of thin brine layer networks during ice growth.

6.6 Relevance for directed percolation in general

Directed percolation has been used to theoretically describe a wide range of processes in statistical physics:^59,60 catalytic reactions, epidemics, calcium dynamics, forest fires, directed polymers, porous media, and even turbulence. However, only very few experimental realisations have been reported.^59,60 While gravity-driven directed percolation of natural porous media had already been discussed by Broadbent and Hammersley,⁵⁵ so far it has not been realised experimentally. To understand why sea ice may exhibit directed percolation behavior, the following brief overview is given following ref. 59.

Directed percolation may be described as a spreading process restricted to a given direction, where activity may either spread over the entire lattice or die out. The inactive state is called absorbing, as it can be reached and not left again, and directed percolation is thus a non-equilibrium process. When interpreting the direction as time, it may be viewed as a dynamic reaction–diffusion process in d + 1 dimensions, governed by the following operations of active (A) and inactive (I) sites:

(i) Diffusion: I + A → A + I

(ii) Self-destruction: A → I

(iii) Offspring production: A → A + A

(iv) Coagulation: A + A → A

Depending on the ratio between (ii) self destruction and (iii) offspring production, the spreading may remain active or reach the absorbing state from where it cannot escape.

The pore space evolution of sea ice may be interpreted as such a reaction–diffusion process. The direction is given by heat flow and gravity resulting in a columnar anisotropic microstructure and pore networks. The non-conservative character is related to thermo- and fluid dynamics: cooling leads to decreasing porosity but also to vertical movement of brine, while gravity drainage exchanges fluid vertically with the ocean underneath. As the ocean salinity is less than that of brine inside pores, this process also decreases (increases) the pore size and volume for upflow (downflow). The change in brine porosity, pore sizes and pore connectivity thus spreads vertically. The redistribution of solute and brine corresponds to diffusion (i), the porosity decrease due to thermodynamics and upward flow corresponds to self destruction (ii), the downward flow and brine transport may create (iii) offsprings and coagulating channels (iv). Sea ice reaches an inactive state below the percolation limit, once salt fluxes between ice and ocean, in concert with thermodynamic transitions, have locked the pore space. It cannot escape from this new absorbing state before it warms again. While the absorbing state is not strictly met in sea ice, which intermittently warms and cools, one expects a fluctuation between active and inactive states. Whereas the present study (as well as ref. 21) cannot be considered as a proof that DP manifests itself in sea ice, the above considerations and the present data analysis support the hypothesis.

6.7 Limitations and future needs

The present study helps to better understand the scatter in observations and simulations of permeability, Fig. 7 and 9 relating it to differences in ice type and growth velocity. How large are the expected uncertainties in the model results and observations? As discussed in Section 5, the equations based on the standard settings for d₀, f_c and t predict a permeability that is only slightly higher (12–22%) than the best fit to simulations. The underestimate of permeability by simulation and experiment, related to finite sample sizes and pore inclination/tortuosity, has a similar magnitude (Section 6.3.2). In granular ice with higher tortuosity the underestimate will be larger. Such underestimates may be reduced by simulations on samples with larger diameter. To keep the spatial resolution this requires micro-CT optics with more pixels.

The model results depend mainly on parametrisations for the plate spacing a₀, the value of d₀, and the critical exponent t. The question if d₀ is indeed constant, and the mechanism that leads to bridging at d₀, needs to be investigated. The growth rate dependence of a₀ is reasonably understood,²⁹ and taking the plate spacing error as 10% results in a 20% permeability error. The effect of t increases with decreasing porosity, e.g., varying t within its confidence range 2.55 ± 0.25 at ϕ = 0.07 and a₀ = 0.55 mm changes the permeability by up to 50%. More data, and also here simulations on larger samples, would improve the confidence of t. However, t is not a true universal exponent, yet is expected to vary with pore size distribution, and progress may be made by studies that constrain t based on certain pore space metrics.^64,65 There is another refinement that one may introduce in the model: the formulation of the columnar permeability model assumed no initial (at ϕ₀) tortuosity. However, as discussed above, typical brine layer inclination angles in young ice are 10–20°. The corresponding permeability reduction, given as cos(α)², is 3–12%.

To model the permeability of granular ice the same critical exponent t as obtained empirically for the columnar ice data has been used. This ad hoc approach was necessary as there is insufficient data for granular surface ice. As the pore space evolution with porosity may be different for granular compared to columnar ice, and as granular ice may be better described by isotropic percolation (due to which argument a 50% higher ϕ_c was proposed), a different exponent t could expected. Also adopting the plate spacing and bridging concept for granular ice may need some modification. However, the simple approach reasonably explains the difference in percolation threshold of granular and columnar ice as indicated by simulations and observations.

Older ice has often experienced severe temperature fluctuations that will change the microstructure from their initial scales (d₀ and a₀) to coarser brine channels. The permeability may then reach values one to two orders of magnitude above the range of young ice.²⁴ Important questions are: Is this transition smooth and does the permeability evolution remain related to the original growth conditions? Or does the permeability of old summer ice become independent of the young ice permeability during its formation, and thus the growth conditions? Some information on pore size evolution during warming and aging of sea ice is available but limited,^38,39,66,91 and more detailed studies are needed.

It should also be noted that percolation theory based on critical exponents is strictly valid only close to the percolation threshold ϕ_c. In practice it often holds at higher porosity, where the effect of ϕ_c due to the (ϕ − ϕ_c)^q dependence is weak. With more exact fits closer to ϕ_c, however, the simplicity of the present model would be lost. Also, the constant f_c approach from ref. 68 is semi-empirical and does not work for all isotropic lattices.⁹² The agreement we find for sea ice might be coincidence and deserves further studies.

7 Conclusions

In this study I have derived equations that predict the permeability of sea ice, and its percolation threshold, in dependence on the growth rate of ice. While eqn (16) to (19) appear rather simple, their derivation is the result of half a century of studies by many investigators, starting with an approach to relate the mechanical behavior of sea ice to its microstructure and a critical porosity ϕ₀ ≈ d₀/a₀.^30,44 A value of ϕ₀ ≈ 0.2 has since then been discussed in many investigations of the mechanical properties of sea ice mechanics.^31,53,93,94 Studies on sea ice thermodynamics and desalination then indicated another critical porosity of about ϕ_c ≈ 0.05², and in later years this transition was interpreted in terms of permeability and percolation theory.^3,19,20,26 Only recently, X-ray microtomographic imaging has provided sufficient data to study these microstructure details, derive a more accurate ϕ_c, and draw conclusions about the relevant (directed) percolation model.²¹ The present study has extended that work by a model for the percolation threshold and permeability that is supported by most observations. It puts the above critical porosities into a unified form ϕ_c = f_cϕ₀ = f_cd₀/a₀, where f_c is based on directed percolation theory, d₀ is the critical brine layer width, and the plate spacing a₀ is related to the growth rate of ice.²⁹

It is anticipated that d₀ and a₀ are relevant for other sea ice properties and processes, some of which were briefly discussed in the present paper. Increasing availability and quality of 3D sea ice microstructure data in recent years^{21,38,39,66,91,95} imply an increasing potential to study the relationships of physical properties, microstructure and growth conditions. Such studies are, in the author’s opinion, inevitable for obtaining a good understanding of the role of sea ice in the environment.

Data availability

The microCT data (170 GB) are available online.²¹

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This project was partly funded through the Research Council of Norway (RCN) programme PETROMAKS2, grant no. 308786.

References

1. W. F. Weeks, On Sea Ice, Univ. of Alaska Press, 2010, p. 600.
2. G. F. N. Cox and W. F. Weeks, J. Geophys. Res., 1988, 93, 12449–12460.
3. C. Petrich, P. J. Langhorne and Z. F. Sun, Cold Reg. Sci. Technol., 2006, 44, 131–144.
4. M. Vancoppenolle, C. M. Bitz and T. Fichefet, J. Geophys. Res.:Oceans, 2007, 112, C04022.
5. S. Maus, Report Series of Geophysics, University of Helsinki, 2008, pp. 99–112.
6. P. J. Griewank and D. Notz, J. Geophys. Res.:Oceans, 2013, 118, 3370–3386.
7. A. K. Turner and E. C. Hunke, J. Geophys. Res.:Oceans, 2015, 120, 1253–1275.
8. D. Rees Jones and M. Grae Worster, J. Geophys. Res.:Oceans, 2014, 119, 5599–5621.
9. I. Melnikov, J. Geophys. Res.:Oceans, 1995, 100, 4673–4680.
10. D. Thomas and G. S. Dieckmann, Sea Ice: an Introduction to its Physics, Chemistry, Biology and Geology, Blackwell, 2003.
11. M. Vancoppenolle, et al., Quat. Sci. Rev., 2013, 79, 207–230.
12. I. Semiletov, A. Makshtas, S. Akasofu and E. L. Andreas, Geophys. Res. Lett., 2004, 31, L05121.
13. M. Angelopoulos, et al., Front. Earth Sci., 2022, 10, 864523.
14. N. Geilfus, B. Delille, J. Tison, M. Lemes and S. Rysgaard, Elem. Sci. Anth., 2023, 11, 00056.
15. T. Koop, A. Kapilashrami, L. T. Molina and M. J. Molina, J. Geophys. Res., 2000, 105, 26393–26402.
16. R. Sander, J. Burrows and L. Kaleschke, Atmos. Chem. Phys., 2006, 6, 4653–4658.
17. G. S. Dieckmann, G. Nehrke, S. Papadimitriou, J. Göttlicher, R. Steininger, H. Kennedy, D. Wolf-Gladrow and D. N. Thomas, Geophys. Res. Lett., 2008, 35, L08501.
18. S. Morin, G. M. Marion, R. von Glasow, D. Voisin, J. Bouchez and J. Savarino, Atmos. Chem. Phys., 2008, 8, 7317–7324.
19. K. M. Golden, S. F. Ackley and V. I. Lytle, Science, 1998, 282, 2238–2241.
20. D. J. Pringle, J. E. Miner, H. Eicken and K. M. Golden, J. Geophys. Res., 2009, 114, C12017.
21. S. Maus, M. Schneebeli and A. Wiegmann, The Cryosphere, 2021, 15, 4047–4072.
22. F. A. L. Dullien, Porous Media: Fluid Transport and Pore Structure, Academic Press, 2nd edn, 1991.
23. D. A. Nield and A. Bejan, Convection in Porous Media, Springer, 2nd edn, 1999.
24. J. Freitag, PhD thesis, Universität Bremen, 1999.
25. T. Maksym and M. O. Jeffries, J. Geophys. Res., 2000, 105, 26313–26331.
26. K. M. Golden, H. E. A. Heaton, J. Miner, D. J. Pringle and J. Zhu, Geophys. Res. Lett., 2007, 34, L16501.
27. W. F. Weeks and W. L. Hamilton, Am. Mineral., 1962, 47, 945–961.
28. M. Nakawo and N. K. Sinha, Atmos.-Ocean, 1984, 22, 193–206.
29. S. Maus, Ann. Glaciol., 2020, 61, 408–425.
30. D. L. Anderson and W. F. Weeks, Trans., Am. Geophys. Union, 1958, 39, 632–640.
31. W. F. Weeks and A. Assur, Ice and Snow, 1963, pp. 258–276.
32. S. Maus, On Brine Entrapment in Sea Ice: Morphological Stability, Microstructure and Convection, Logos, Berlin, 2007.
33. F. Flin, J. B. Brzoska, B. Lesaffre, C. Coleou and R. A. Pieritz, Ann. Glaciol., 2004, 38, 39–44.
34. M. Schneebeli and S. A. Sokratov, Hydrol. Processes, 2004, 18, 3655–3665.
35. J. Freitag, F. Wilhelms and S. Kipfstuhl, J. Glaciol., 2004, 50, 243–259.
36. M. W. Hörhold, M. R. Albert and J. Freitag, J. Glaciol., 2009, 55, 625–630.
37. R. Obbard, G. Trodermann and I. Baker, J. Glaciol. 2009, 55, 1113–1115.
38. S. Maus, J. Becker, S. Leisinger, M. Matzl, M. Schneebeli and A. Wiegmann, 23rd Proc. Port and Ocean Engineering under Arctic Conditions, Trondheim, Norway, 2015, https://www.poac.com.
39. M. Oggier and H. Eicken, J. Glaciol., 2022, 68, 833–848.
40. W. Scoresby, Mem. Wernerian Soc., 1815,(2), 328–336.
41. D. Walker, Proc. R. Soc. London, 1859, 9, 609–611.
42. F. Ruedorff, Ann. Phys., 1861, 190, 63–81.
43. T. Fukutomi, M. Saito and Y. Kudo, Low Temp. Sci., 1952, 9, 113–123.
44. A. Assur, Arctic Sea Ice, 1958, pp. 106–138.
45. W. Schwarzacher, J. Geophys. Res., 1959, 64, 2357–2367.
46. N. V. Cherepanov, Tr. Arkt. Antarkt. Nauchno-Issled. Inst., 1964, 267, 13–18.
47. N. K. Sinha and C. Zhan, J. Mater. Sci. Lett., 1996, 15, 2118–2121.
48. W. F. Weeks and G. Lofgren, Physics of Snow and Ice, Sapporo, Japan, 1967, pp. 579–597.
49. S. Maus, in Physics and Chemistry of Ice, ed. W. F. Kuhs, Royal Society of Chemistry, 2007, pp. 371–382.
50. W. W. Mullins and R. F. Sekerka, J. Appl. Phys., 1964, 35, 444–451.
51. G. A. Maykut, J. Geophys. Res., 1978, 83, 3646–3658.
52. G. F. N. Cox and W. F. Weeks, J. Glaciol., 1974, 13, 109–120.
53. S. Maus, Proceedings - 27th Conference on Port and Ocean Engineering under Arctic Conditions, Glasgow, Scotland, 2023, p. 41, https://www.poac.com.
54. W. F. Weeks and D. L. Anderson, Trans., Am. Geophys. Union, 1958, 39, 641–647.
55. S. R. Broadbent and J. M. Hammersley, Math. Proc. Cambridge Philos. Soc., 1957, 53, 629–641.
56. D. Stauffer and A. Aharony, Introduction to Percolation Theory, Taylor & Francis, 2nd edn, 1992.
57. M. Sahimi, Rev. Mod. Phys., 1993, 65, 1393–1534.
58. S. Maus, S. Bahafid, M. Hendriks, S. Jacobsen and M. Geiker, Cold Reg. Sci. Technol., 2023, 208, 103780.
59. H. Hinrichsen, Adv. Phys., 2000, 49, 815–958.
60. M. Henkel, H. Hinrichsen and S. Lübeck, Non-equilibrium Phase Transitions. Volume I: Absorbing Phase Transitions, Springer, 2008.
61. A. Sur, J. Lebowitz, J. Marro, M. Kalos and S. Kirkpatrick, J. Stat. Phys., 1976, 15, 345–353.
62. H. Janssen and U. Täuber, Ann. Phys., 2005, 315, 147–192.
63. H. Janssen and O. Stenull, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2001, 63, 025103(R).
64. A. J. Katz and A. H. Thompson, Phys. Rev. B:Condens. Matter Mater. Phys., 1986, 34, 8179–8181.
65. P. Le Doussal, Phys. Rev. B:Condens. Matter Mater. Phys., 1989, 39, 4816–4819.
66. M. Salomon, S. Maus and C. Pertrich, J. Glaciol., 2022, 68, 571–590.
67. GeoDict, Geometric material Models and Computational PreDictions of Material Properties, 2012–2024, https://www.geodict.com.
68. H. Scher and R. Zallen, J. Chem. Phys., 1970, 53, 3759–3761.
69. D. S. Gaunt and M. F. Sykes, J. Phys. A: Math. Gen., 1983, 16, 783–799.
70. S. Lübeck and R. Willmann, J. Stat. Phys., 2004, 115, 1231–1250.
71. P. Grassberger, J. Stat. Mech.:Theory Exp., 2009, 2009, P08021.
72. J. Wang, Z. Zhou, Q. Liu, T. Garoni and Y. Deng, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2013, 88, 042102.
73. M. Shokr and N. Sinha, Sea Ice Physics and Remote Sensing, John Wiley and Sons, 2015.
74. N. Ono and T. Kasai, Ann. Glaciol., 1985, 6, 298–299.
75. T. Kawamura, Contrib. Inst. Low Temp. Sci., Hokkaido Univ., Ser. A, 1988, 3120, 1–29.
76. J. Buettner, MSc thesis, University Bergen, 2011.
77. R. Morey, A. Kovacs and G. F. N. Cox, Cold Reg. Sci. Technol., 1984, 9, 53–75.
78. K. Jones, M. Ingham and H. Eicken, J. Geophys. Res., 2012, 117, C02005.
79. G. F. N. Cox and W. F. Weeks, J. Glaciol., 1986, 32, 371–375.
80. N. Jeffery, E. Hunke and S. Elliott, J. Geophys. Res., 2011, 116, C07020.
81. J. R. Urrego-Blanco, N. Urban, E. C. Hunke, A. Turner and N. Jeffery, J. Geophys. Res.:Oceans, 2016, 121, 2709–2732.
82. V. L. Tsurikov, Okeanologiia, 1965, 5, 463–472.
83. C. Renshaw, O. Marchenko, E. Schulson and E. Karulin, J. Glaciol., 2018, 64, 443–449.
84. E. P. Jones and A. R. Coote, J. Geophys. Res., 1981, 86, 11041–11043.
85. S. Rysgaard, R. N. Glud, M. K. Sejr, J. Bendtsen and P. B. Christensen, J. Geophys. Res., 2007, 112, C03016.
86. L. G. Anderson, E. Falck, E. P. Jones, S. Jutterström and J. H. Swift, J. Geophys. Res., 2004, 109, C06004.
87. W. R. Simpson, R. von Glasow and K. Riedel, et. al., Atmos. Chem. Phys., 2007, 7, 4375–4418.
88. A. Saiz-Lopez., C. S. Blaszczak-Boxe and L. Carpenter, Atmos. Chem. Phys., 2015, 15, 9731–9746.
89. M. Gourdal, O. Crabeck, M. Lizotte, V. Galindo, M. Gosselin, M. Babin, M. Scarratt and M. Levasseur, Elem. Sci. Anth., 2019, 7, 33.
90. E. Damm, S. Thoms, M. Aneloppulos, L. von Albedyll, A. Rinke and C. Haas, Front. Earth Sci., 2019, 12, 1338246.
91. C. M. Frantz, B. Light, S. Farley, S. Carpenter, R. Lieblappen, Z. Courville, M. V. Orellana and K. Junge, The Cryosphere, 2019, 13, 775–793.
92. I. Balberg and N. Binenbaum, Phys. Rev. B:Condens. Matter Mater. Phys., 1987, 35, 8749–8752.
93. B. Michel, Ice Mechanics, Les Presses De L’Université Laval, Québec, 1978.
94. J. A. Richter-Menge and K. F. Jones, J. Glaciol., 1993, 39, 609–618.
95. R. Lieb-Lappen, E. Golden and R. Obbard, Cold Reg. Sci. Technol., 2017, 138, 24–35.

---