Numerical assessment of morphological and hydraulic properties of moss, lichen and peat from a permafrost peatland

Due to its insulating and draining role, assessing ground vegetation cover properties is important for high-resolution hydrological modeling of permafrost regions. In this study, morphological and effective hydraulic properties of Western Siberian Lowland ground vegetation samples (lichens, Sphagnum mosses, peat) are numerically studied based on tomography scans. Porosity is estimated through a void voxels counting algorithm, showing the existence of representative elementary volumes (REVs) of porosity for most samples. Then, two methods are used to estimate hydraulic conductivity depending on the sample's homogeneity. For homogeneous samples, direct numerical simulations of a single-phase flow are performed, leading to a definition of hydraulic conductivity related to a REV, which is larger than those obtained for porosity. For heterogeneous samples, no adequate REV may be defined. To bypass this issue, a pore network representation is created from computerized scans. Morphological and hydraulic properties are then estimated through this simplified representation. Both methods converged on similar results for porosity. Some discrepancies are observed for a specific surface area. Hydraulic conductivity fluctuates by 2 orders of magnitude, depending on the method used.

Porosity values are in line with previous values found in the literature, showing that arctic cryptogamic cover can be considered an open and well-connected porous medium (over 99 % of overall porosity is open porosity). Meanwhile, digitally estimated hydraulic conductivity is higher compared to previously obtained results based on field and laboratory experiments. However, the uncertainty is less than in experimental studies available in the literature. Therefore, biological and sampling artifacts are predominant over numerical biases. This could be related to compressibility effects occurring during field or laboratory measurements. These numerical methods lay a solid foundation for interpreting the homogeneity of any type of sample and processing some quantitative properties' assessment, either with image processing or with a pore network model. The main observed limitation is the input data quality (e.g., the tomographic scans' resolution) and its pre-processing scheme. Thus, some supplementary studies are compulsory for assessing syn-sampling and syn-measurement perturbations in experimentally estimated, effective hydraulic properties of such a biological porous medium.

1 Introduction

Covering a quarter of the Northern Hemisphere's land surface (Brown et al., 1997), permafrost soils are the most representative soil types in arctic and subarctic regions. Permafrost is a soil layer in which temperature remains below 0 ^∘C for at least 2 consecutive years, thus holding ice in its porous structure. Frozen layers make permafrost hydrology peculiar, resulting in complex couplings between heat and water fluxes (Grenier et al., 2018; Tananev et al., 2020). Seasonal structural variations occur in permafrost soils, as surface thawing forms an “active layer”. Most permafrost-related biogeochemical processes (especially organic matter degradation) take place in this layer. The active layer is at its maximum thickness in the early fall and is generally meters in scale (Clayton et al., 2021; Aalto et al., 2018; Guo and Wang, 2017). Active layer thickness is, nonetheless, spatially variable due to climatic conditions, land cover, and the micro- and macro-topography. The impact of hydrological climate change is particularly drastic in permafrost-dominated environments because of deepening thaw fronts (Hinzman and Kane, 1992). Indeed, between 2008 and 2016, the average annual temperature of arctic permafrost soil increased by 0.4 (±0.25) ^∘C (Biskaborn et al., 2019; Fox-Kemper et al., 2023). This causes positive feedback on average atmospheric temperatures (Meredith et al., 2019), reduces latent heat effects (Walwoord and Kurylyk, 2016) and increases water drainage in arctic watersheds (Liljedahl et al., 2016). Hence, quantifying heat and water transfer properties in permafrost-affected regions is compulsory.

Previous studies have addressed this quantification through field observations (Olefeldt and Roulet, 2014; Streletskiy et al., 2015; Throckmorton et al., 2016; O'Connor et al., 2020) or field and laboratory experiments (Vedie et al., 2011; Roux et al., 2017; Wagner et al., 2018). Some recent studies have also dealt with this question using a modeling approach (Bense et al., 2012; Genxu et al., 2017; Burke et al., 2020; Du et al., 2020; Fabre et al., 2017).

Bryophytes (mosses) and lichens are widely present in tundra and taiga environments. The dominant ground cover consists of Sphagnum mosses and lichens in permafrost peatlands (Soudzilovskaia et al., 2013; Volkova et al., 2018). Sphagnum mosses are part of the Bryophyta plant division, which represents non-vascular plants (without xylem or phloem). Sphagnum colonies grow indeterminately from their apical structure, named the capitula. Their water content mainly relies on capillary forces maintained by each individual's density (Hayward and Clymo, 1982; Howie and Hebda, 2018). Lichens are not “vegetation” but consist of a symbiotic association between heterotrophic Fungi and autotrophic Algae. Both Sphagnum and lichens can be gathered into the Cryptogamae subkingdom. This cryptogamic layer has an important impact on permafrost dynamics because it influences heat and water exchanges between the soil and atmosphere (Soudzilovskaia et al., 2013; Launiainen et al., 2015; Porada et al., 2016; Park et al., 2018; Loranty et al., 2018). Boreal vegetation is assumed to be a major nutrient and inorganic solute exchange medium at a watershed scale (Shirokova et al., 2021). Boreal vegetation is likely to accumulate in lowlands at a low degradation rate, resulting in the formation of Sphagnum peatlands, such as the Western Siberian Lowlands.

Ground vegetation transfer properties are key information for high-resolution hydrological modeling of permafrost-related catchments. Thus, reliable estimates of them are necessary for water flux studies for boreal soils and for climate change impact assessment of the hydrology of high-latitude continental surfaces. Therefore, some recent efforts have been put into emphasizing the role of the cryptogamic layer in Earth system models (Stepanenko et al., 2020; Shi et al., 2021). Devoted modeling tools have also been created to predict Sphagnum dynamics (the Peatland Moss Simulator by Gong et al., 2020). Furthermore, specific modeling work has been conducted on restored Sphagnum peatlands to link hydrological properties with dissolved organic carbon dynamics (Bernard-Jannin et al., 2018) or soil moisture dynamics (Elliott and Price, 2020). However, the mechanistic modeling of water and heat fluxes in ground vegetation layers remains difficult, as their porous-medium transfer properties are not straightforward to evaluate (Orgogozo et al., 2019).

Many studies are available for the decayed Sphagnum layer: peat. The hydrological and thermal properties of peat are well documented. Extensive reviews of the relation between hydrogeochemical processes in peatlands and peat's porous-medium structure were conducted by McCarter et al. (2020). A study of peatland's hydraulic properties was initiated during the 1920s for peatland drainage (Malmström, 1925). Then, some introductory field experiments were conducted on Finnish peatlands (Virta, 1962; Heikurainen, 1963; Sarasto, 1963) as well as in the United States (Boelter, 1968) and Ireland (Galvin, 1976).

Only a few studies were conducted on the living part of this upper permafrost layer. Hence, quantitative assessments of some key hydrological properties of ground vegetation layers are needed, such as total, open and enclosed porosity, hydraulic conductivity and specific surface area. In terms of hydraulic properties, hydraulic conductivity has been assessed in the laboratory using constant or falling-head permeameters (Quinton et al., 2000; Price et al., 2008; Hamamoto et al., 2016; Weber et al., 2017) or via field measurements (Päivänen, 1973; Crockett et al., 2016; this study). The results are presented in Table 2, with some peat results for comparison. Otherwise, arctic lichens have received little attention to date. To our knowledge, only one study has estimated lichens' hydraulic properties, considering unsaturated hydraulic conductivity without taking into account macropores (Voortman et al., 2014). However, the specific surface areas of some other lichen species are documented in the literature (Adamo et al., 2007). Some studies quantified arctic lichen properties in response to acid rain (Tarhanen et al., 1999) to clarify their interaction with the rhizosphere (Banfield et al., 1999) or in relation to their albedo properties (Bernier et al., 2011). Some transmembrane transfer properties are also available in the literature (Potkay et al., 2020)

Table 1 Notation glossary.

Table 2 Synthesis of saturated hydraulic conductivity (m s⁻¹) of peat and Sphagnum found in the literature. This study's values refer to field experiments conducted during sample collection (CHP: constant head permeameter; FHP: falling head permeameter; FP: field percolation; IM: inverse modeling; NM: numerical model, Hydrus-1D; see McCarter and Price, 2012).

Thus, many field and experimental studies are available throughout the literature. However, field work and experimental studies are known to bring their own difficulties, due to the local conditions' variabilities, sampling biases, disturbances and measurement uncertainties. The aim of this study is to assess some of the transfer properties that are well documented with field and experimental studies with an innovative numerical scheme. Indeed, the use of numerical workflows enhances reproducibility and intercomparison capability between samples. Numerical workflows are often used when experimental studies are complicated to implement (reservoir engineering, aerodynamics, micro-fluidics).

In this work, such numerical workflow is intended to be used for evaluating the hydrological transfer properties of representative vegetation types of the Western Siberian Lowlands. To this end, natural samples collected from the Western Siberian Lowlands are digitally analyzed to characterize some morphological and hydraulic transfer properties. Thus, in contrast to previous works compiled in Table 2, this study aims to assess the hydraulic properties of lichens and Sphagnum mosses by numerical methods rather than experimental measurements. The arctic cryptogamic layer is assumed, hereafter, to represent a complex patchwork of biological porous media (Price et al., 2008, Voortman et al., 2014, Hamamoto et al., 2015).

To validate this hypothesis, a thorough analysis of sample homogeneity is carried out, based on porosity, as it is the main driver of flow dynamics in porous media (Koponen et al., 1996, 1997). This enables the classification of samples according to their homogeneity. Indeed, for homogeneous samples, a smaller sample region can be considered an effective medium sharing the same properties as the whole sample. Multidimensional porosity description leads to a statistical study of the existence of a representative elementary volume (REV). Two standard porous-medium modeling methodologies are used throughout this study: direct numerical simulations on computed representative elementary volumes (DNS-REVs) and pore network modeling with a built-in solver (PNM). The impossibility of collecting a substantial number of samples is compensated for by a statistical quantification of a REV for each sample. This implies that the REV is smaller than the sample, and hence sampling size is chosen to match sizes that were used in the previous literature, such as Weber et al. (2017).

2 Material and methods

2.1 Sample collection and digital reconstruction

Samples were collected at Khanymei Research Station (N63^∘43^′19,73^′′ E75^∘57^′47,91^′′) in Western Siberia (Autonomous District of Yamal-Nenets, Russian Federation) in July 2018. The mean annual temperature is −5.6 ^∘C, and average precipitation is 540 mm (Payandi-Rolland et al., 2020; Raudina et al., 2018). Eight moss samples (Sphagnum sp.) were collected, either on moss mounds or in thermokarstic hollows. Additionally, two lichen samples (Cladonia sp., named Lichen1.3 and Lichen2.1) and two peat samples (named Peat2.2 and Peat2.3) were collected. Of the eight Sphagnum samples, three are S. angustifolium (C.E.O. Jensen ex Russow, named Mound2.4, Mound2.5 and Mound2.6), two are S. lindbergii Schimp (named Hollow1.2 and Hollow1.4) and two are S. majus (Russow) (C.E.O. Jensen, named Hollow2.7 and Hollow2.8). The last moss sample is S. lenense H. Lindb. ex L.I. Savicz (named Mound1.1).

Sampling is thoroughly conducted to minimize structural perturbations. In order to achieve this, each sample's surroundings are cleared with special care prior to extraction. Then, the sample is extracted using a ceramic knife directly at the right dimensions to fit in a high-density polyethylene box, where it remains from the moment of sampling and drying to the tomographic examination. Additionally, four in situ hydraulic conductivity measurements are performed on various Sphagnum plots, using a double-ring infiltrometer (Table 2). An overview of the sample collection method is shown in Fig. 1 as well as 3D tomographical visualizations of each sample type.

Figure 1 Sampling collection method overview.

The samples are then dried at 40 ^∘C for 48 h after sampling. Thereafter, each sample is scanned and digitally reconstructed using high-resolution X-ray imagery.

X-ray computed tomography (X-CT) has been widely studied and is extensively used for medical purposes and geoscientific applications (Christe et al., 2011). Tomography is a non-destructive technique which enables the observation of pore structure data at micron scale, especially for pore space assessment in sedimentary rocks. X-CT scanning has been acknowledged as being an efficient method for accessing morphological information, such as the pore structure of peat soils (Turberg et al., 2014). Cnudde and Boone (2013) published an exhaustive review of X-ray tomography applications for the Earth sciences. Rezanezhad et al. (2016) demonstrated that X-CT peat scanning showed a satisfactory spatial resolution for the study of peat's pore morphology. Since bryophytes can be assumed to represent a cluster of individuals, X-CT permits the segmenting of each plant structure, which cannot otherwise be achieved without destructive techniques. Tomographical scans of studied samples are produced using EasyTom^®XL (RX Solutions, France) with a maximal X-ray emission source set to 90 kV. The obtained resolution, after tridimensional reconstruction, is 94 µm voxel⁻¹, except for “Lichen2.1” at 88 µm voxel⁻¹ (due to the scanning settings used specifically for this sample). The voxel number ranges from 4.8×10⁸ to 9.6×10⁸ voxels. Virtual samples are cropped and reduced to form a usable volume (UV) to avoid sampling border effects. Using ImageJ-Fiji (Schindelin et al., 2012), computed samples are then binarized using an intra-class variance-reducing algorithm (Otsu, 1979). This resulted in a sample consisting of an 8-bit black and white image stack. Supplement S1 contains some of the technical data, such as usable volumes and digital reconstructions, for each sample.

2.2 Drying impact assessment of sample representativity

The sampling locations and processing facilities were far away from each other. To ensure structural preservation, special care is taken throughout the sampling, transportation and scanning operations. The samples are oven-dried for 48 h at 40 ^∘C, at atmospheric pressure, to halt biological degradation. As expected, Sphagnum mosses began to whiten and become papery, as described by Hayward and Clymo (1982). The samples are scanned at the Toulouse Institute of Fluid Mechanics under dry conditions, 2 months after their primary extraction at Khanymei Research Station and drying.

To ensure the dry samples' representativity, we used an analogous drying experimental protocol to the one carried out by Kämäräinen et al. (2018). This experiment is conducted on similar Sphagnum species (S. fuscum (Schimp.) H.Klinggr., S. majus (Russow) C.E.O. Jensen) and sampled according to the same method at Clarens mire (southwestern France, N43^∘08^′41.3", E0^∘25^′12.9") in March 2021, notwithstanding that Sphagnum acrotelm (growing section of a Sphagnum individual) is much thinner than Siberian samples. Four samples are collected, two being dried at 40 ^∘C for 48 h and the other two being left untouched, as control samples. Each of the four samples is scanned 2 d after their extraction and then again 14 d after extraction. Additionally, one Sphagnum majus individual is extracted and left to dry under ambient conditions.

A comparative study between each of the two sample lots and the lone individual show that drying does not affect structural preservation. Our validation experiment converges with the results found by Kämäräinen et al. (2018). This also confirms hyaline cells' structural durability: the early work of Puustjärvi (1977) showed that hyaline cells were well preserved during biological decay. Drying impacts aside, Sphagnum's continuous growth on non-dried control samples seems the most impactful structural factor, as each individual was striving to adapt itself to the sampling box hydric conditions. Fast drying before tomographic examination can be a reasonable solution for preserving the morphological structure, in conjunction with careful on-site sampling.

2.3 Morphological analysis: total porosity (ε_total), open porosity (ε_open), specific surface area (S_SA) and pore size distribution

Global porosity (ε_Total) is calculated for each sample using built-in ImageJ-Fiji tools and macro-scripting. Porosity is considered to be a ratio between the number of voxels representing the void phase i_voxel (void phase's internal volume, including closed porosity) over the total number of voxels representing a sample N_total (void and matrix volume). This relation is shown in Eq. (1):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathit{ϵ}}_{\mathrm{Total}}& ={\displaystyle \frac{V_{i=\mathrm{0}}}{V_{\mathrm{Total}}}}\times \mathrm{100}={\displaystyle \frac{\int \int \underset{\mathrm{0}}{\overset{N_{i=\mathrm{0}}}{\int }}{i}_{\mathrm{voxel}}\left(x,y,z\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z}{\int \int \underset{\mathrm{0}}{\overset{N_{\mathrm{Total}}}{\int }}{i}_{\mathrm{voxel}}\left(x,y,z\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z}}\\ & \times \mathrm{100}={\displaystyle \frac{\sum _{\mathrm{0}}^{N_{i=\mathrm{0}}}{i}_{\mathrm{voxel}}}{\sum _{\mathrm{0}}^{N_{\mathrm{Total}}}{i}_{\mathrm{voxel}}}}\times \mathrm{100}.\end{array}\end{array}$$

Porosity is computed on bidimensional horizontal slices along the z axis to evaluate porosity variations along the samples. Image stacks are then reconstructed along the x and y axes to create two other image stacks. Finally, porosity is computed along the x, y and z axes using a voxel-counting algorithm shown in Sect. 3.1.

The samples could then be classified into three types according to the porosity profile along the vertical axis (Table 3 and Supplement S2). As porosity appears to be almost constant over the x and y axes, sample classification is solely based on vertical porosity z.

• Type I: constant high porosity along the z axis, excluding border effects

• Type II: low basal porosity, linearly increasing to the top of the sample

• Type III: no specific trend observed in vertical porosity

For each sample nature as well as each classified type, a dedicated color palette is chosen and kept consistent throughout the study for the sake of clarity.

Table 3 Computed global porosity (ε_total), ratio of open porosity (p_open), classification (I, II, II) and average planar porosity (%) for each sample obtained using a voxel-counting algorithm.

Open and connected porosities (ε_open) are retrieved using dedicated shape analysis and labeling tools provided in the IPSDK™image-processing toolkit (a Reactiv'IP product used in Goubet et al., 2021). This enables a precise segmentation to associate each connected void space into a unique identifier. Here, since the samples have more void than matter, this first label is assumed to be connected void space, which plays a major role in the flow and transfers (porosity), the latter being a closed or non-communicating element. From the raw dataset, voxel intensity is integrated to get the first label's voxel sum divided by the overall voxel number, as shown in Eq. (2):

$$\begin{array}{}\text{(2)}& p_{\mathrm{open}}={\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{open}}}{{\mathit{ϵ}}_{\mathrm{Total}}}}={\displaystyle \frac{\sum _{\mathrm{0}}^{N_{i\in \text{label\hspace{0.17em}0}}}{i}_{\mathrm{voxel}}}{\sum _{\mathrm{0}}^{N_{\mathrm{Total}}}{i}_{\mathrm{voxel}}}}\times \mathit{ϵ}{{}_{\mathrm{Total}}}^{-\mathrm{1}}.\end{array}$$

The specific surface area is deduced using the same shape analysis and labeling tools included in IPSDK™. Integrating the surface between both phases (void and solid) yields the total surface S. Thus, volumetric-specific surface area S_SA is obtained by dividing this surface by the sample's bounding box volume, expressed in m² m⁻³, as shown in Eq. (3):

$$\begin{array}{}\text{(3)}& S_{{\mathrm{SA}}_V}={\displaystyle \frac{S_{\mathrm{solid}}}{V_{\mathrm{BBox}}}}={\displaystyle \frac{\underset{\mathrm{0}}{\overset{N-\mathrm{1}}{\iint }}S\left(\left(i,j\right)\left(i+\mathrm{1},j+\mathrm{1}\right)\right)\mathrm{d}i\mathrm{d}j}{L_X{L}_Y{L}_Z}}.\end{array}$$

Specific surface area is conventionally expressed in relation to density (m² g⁻¹). For this purpose, each dried sample mass is obtained using an analytical balance and the sample's dry bulk density ρ_dry. Then, volumetric-specific surface values are converted into a mass-related specific surface by dividing volumetric-specific surface area by dry density, as shown in Eq. (4):

$$\begin{array}{}\text{(4)}& S_{{\mathrm{SA}}_M}={\displaystyle \frac{S_{{\mathrm{SA}}_V}}{{\mathit{\rho }}_{\mathrm{dry}}}}.\end{array}$$

Pore size distribution is calculated using ImageJ-Fiji's implemented image segmentation tools on the binarized image stacks. In each stack's image, a Euclidean distance transformation of the matrix phase from the void phase is first applied. Then, for each isolated void patch, the Feret diameter is computed.

2.4 Darcy-scale morphological and hydrological property definition: REV

In this study, the collected samples are assumed to form a complex fibrous porous medium. Resolving mechanistic equations in such large domains is not straightforward due to the extensive computational resources required. Conversely, resolving such equations on an arbitrary cropped sample would not aid the hydraulic property assessment. To make the link between microscale and macroscale phenomena, a reproducible pattern is required to avoid microscale heterogeneities and lack of information due to a diminutive sample size. To do this, finding a representative region that validates scale separation assumptions with both microscale and macroscale heterogeneities is compulsory, thus defining the volumetric average of a microscale property that is continuous and informative at a macroscale. One of the first volume-averaging methods consists of finding a statistical REV for the given studied property.

Indeed, REV is a theoretical concept clarifying the definition of the macroscopic scale (Darcy scale) and the microscopic scale (pore scale) and characterizing a given porous medium. This REV can be assumed to be a specific sample volume in which transfer-governing equations (single-phase flow, for example) may be defined along with the associated effective properties. A proper mathematical definition of a REV is given in Bachmat and Bear (1987), Quintard and Whitaker (1989) and Whitaker (1999) along with a thorough definition of volume-averaging methods. A generic profile for a given property φ_β is shown in Fig. 2.

Figure 2 Schematic representation of fluctuations of a generic property φ_β in conjunction with volume (adapted following Brown and Hsieh, 2000).

The fluctuation profile shows three main domains. Here, the REV is defined as the smallest volume for which statistical fluctuations of a given property in a given space are sufficiently low to consider its average value an effective property. Finding the representative elementary volumes of some key properties (e.g., porosity and intrinsic permeability) is a routine workflow in porous-medium sciences. It is often used for fractured oil reservoirs (Durlofsky, 1991) or artificially packed glass bead media (Leroy et al., 2008). A REV is, by definition, large when compared to characteristic lengths of heterogeneities at a microscopic scale but small when compared to characteristic lengths of heterogeneities at the macroscopic scale. Thus, the properties computed for a REV of a porous medium may be defined and computed as continuous functions of space and even constant, in the case of a homogeneous porous medium, as defined by Bear (1972). In general, REVs are described on the basis of morphological characteristics such as porosity, although a distinct REV can be found for any given porous-medium property. Porosity and hydraulic-conductivity-related REVs are characterized throughout this study, leading to two different sizes, one for each property.

2.4.1 Porosity: binarization and voxel counting

From previously binarized image stacks, a statistical REV analysis is conducted using dedicated high-performance image processing Python libraries (IPSDK™), encapsulated in a specifically designed batch process for which the flowchart is shown in Fig. 3.

Figure 3 Flowchart of the representative elementary volume of porosity (REV_ε) from a binarized image. For each sample, three thresholds are tested (5 %, 3 % and 1 % and standard deviation variation).

First, porosity (Eq. 1) is computed for a given sub-sampling volume within the whole sample. Then, the sub-sampling volume location is incrementally reduced and moved in every spatial direction. For each sub-volume, intermediate porosities are computed. The average and standard deviations are stored for each chosen sub-sampling volume. Then, an algorithmic routine is used to find the maximal size that satisfies a given threshold (1 %, 3 % or 5 % of porosity fluctuation). These thresholds define the statistical representativity of these REVs. Thus, a REV satisfying a 1 % threshold can be assumed to be a high-grade REV, whereas the 5 % threshold corresponds to lower-grade REVs. For 10 of the 12 studied samples, a REV of porosity is found. The two remaining samples, Hollow1.2 and Peat2.2, do not exhibit a REV for the chosen thresholds. A collection of tridimensional reconstructions of the samples as well as some examples of REVs for each sample are shown in Supplement S1. The graphical plots for each studied sample are available in Supplement S3. The black dot represents the smallest reached representative elementary volume for each sample.

2.4.2 Hydraulic conductivity: direct numerical simulation

Hydraulic conductivity is estimated through single-phase flow computations performed by solving a Navier–Stokes equation in the pore space of the considered sample. The concept is to carry out the numerical simulation of fluid flow, reproducing the conditions occurring in a constant-head (CHP) permeameter. Then, a sample's hydraulic conductivity is computed from the obtained velocity field. A virtual CHP is created by imposing a constant pressure on two opposite faces to one direction (inlet and outlet). Water-tight wall boundary conditions are applied to other faces, as shown in a conceptual representation of the initial and boundary conditions available in Supplement S2.

Due to computation time limitations, the biggest studied sub-volume with this approach corresponds to a quarter of the total sample. In Sect. 2.4, we stated that the REVs of effective physical properties were valid for that particular physical property. Thus, a hydraulic conductivity REV is required to statistically assess hydraulic conductivity. For that purpose, instead of counting voxel value algorithms (as made for porosity), retrieving a representative elementary volume for hydraulic conductivity requires extensive fluid mechanics simulations. Here, a laminar single-phase flow induced by a pressure gradient is computed for each sub-sample, being consistent with the idea of reducing and moving a defined sub-volume inside the overall sample. From open-porosity data in Table 3, we can show that most porosity is connected to each one. This leads to the assumption that considering all pores can be considered effective in permeability-driven phenomena. As these simulations are resource-costly, Type-I samples (constant porosity) are selected as they are sufficiently homogeneous for the establishment of REVs. Other types are treated by another method presented in Sect. 2.4.3. The implemented method relies on Mohammadmoradi and Kantzas (2016) in conjunction with automatic mesh manipulation tools (trimesh Python library, Dawson-Haggerty et al., 2019). For each Type-I sample, single-phase flow simulation through a fraction of the solid volume (representing a sample) is conducted. The simpleFoam algorithm is used to solve a velocity field from an initial pressure gradient. The simpleFoam solver is an enhanced version of the original SIMPLE algorithm (Semi-Implicit Method for Pressure-Linked Equations – Patankar, 1980) nested in the open-source Computational Fluid Dynamics toolkit OpenFOAM (Weller et al., 1998, https://openfoam.org/, last access: 17 January 2023; https://www.openfoam.com/, last access: 17 January 2023).

For each sample, four potential REV sizes are computed (23.5, 15.7, 11.8 and 9.4 mm), consisting of 8, 27, 64 and 125 simulations on the x, y and z axes, respectively, representing 672 simulations per sample. These sizes are chosen from the beginning to limit computation times, as using the same scanning method than what was developed for porosity is computationally prohibitive. This is run on the tier-2 supercomputer Olympe (CALMIP computational mesocenter, Toulouse, France). These calculations are run simultaneously, each occupying one node (36 physical cores), representing 10 500 h CPU (about 12 d of physical time) per sample. For each simulation, the velocity field u_i is integrated with the overall outlet surface S_outlet (including the surface occupied by the solid matrix) to get an averaged outlet flux value v_i, according to Eq. (5):

$$\begin{array}{}\text{(5)}& v_i={\displaystyle \frac{\mathrm{1}}{S_{\mathrm{outlet}}}}\underset{\mathrm{0}}{\overset{S}{\int }}{u}_i\mathrm{d}S.\end{array}$$

A careful convergence study is also conducted so that numerical errors, associated with discretization resolutions and iterative procedures for the approximated inversions of the linear systems involved, are low enough to be neglected in the analysis of the results. Inlet pressures are chosen to avoid turbulent flows (Re <<1). The computed Darcy velocity v_i could then be injected into a regular Darcy law, as shown in Eq. (6), where k_ii is a tensorial component of intrinsic permeability (m²) and μ_w is the dynamic viscosity:

$$\begin{array}{}\text{(6)}& k_{ii}=v_i{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{w}}}{\mathrm{\nabla }P}}\mathrm{with}\mathrm{\nabla }P={\displaystyle \frac{P_{\mathrm{inlet}}-P_{\mathrm{outlet}}}{L_i}}.\end{array}$$

To avoid artifacts related to the physics of a specific fluid, the simpleFoam solver uses kinematic pressure (expressed in m² s⁻²) and kinematic viscosity ν (m² s⁻¹) to solve Navier–Stokes equations. These equations are based on intrinsic permeability k expressed in m². However, in the field of hydrology the hydraulic conductivity (m s⁻¹), abbreviated to K_w, is generally used. One can relate hydraulic conductivity K_w to intrinsic permeability k by using Eq. (7) described by Claisse (2016).

$$\begin{array}{}\text{(7)}& k={\displaystyle \frac{K_{\mathrm{w}}{\mathit{\mu }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{w}}g}}\end{array}$$

In continental surface hydrology, liquid water's physical property variations (e.g., volumetric mass ρ_w and dynamic viscosity μ_w) are generally neglected. Thus, intrinsic permeability values obtained from the numerical computations were converted using water's thermodynamic properties at 293.15 K and 1.013 kPa (Chemical Rubber Company and Lide, 2004), considering the following conversion equation (Eq. 8):

$$\begin{array}{}\text{(8)}& k=\mathrm{1.0217}\times {\mathrm{10}}^{-\mathrm{7}}{K}_{\mathrm{w}}.\end{array}$$

This method is suitable for samples meeting porosity homogeneity requirements, classified into Type-I samples. However, another method is needed to compensate for Type-II and Type-III sample heterogeneity, as using direct numerical simulations on a complete usable volume is prohibitive in terms of computational resources. The results obtained for hydraulic conductivity REV computations are available in Supplement S4.

A double-ring infiltrometry test was also conducted during the sampling campaign. For the sake of comparison, hydraulic conductivity values obtained using this method are also shown in Table 2.

Table 4 Obtained representative elementary volume based on porosity (REV_ε). $L_{\mathit{\epsilon }}^{\mathrm{REV}}$ is the side length of a cubic REV of porosity. ${\mathit{\epsilon }}_{\mathrm{mean}}^{\mathrm{REV}}$ is the average porosity of a given cubic REV. The ratio represents the volumetric percentage of the sample included in the REV.

2.4.3 No REV for hydraulic conductivity: use of PNM

For the samples that do not exhibit a REV for hydraulic conductivity (Type-II and Type-III samples), the hydraulic conductivity is then studied using a pore network model, generated from the binarized image stacks. Pore network models are based on the structural simplification of a complex pore structure (rocks or reactive porous industrial media, for example) into a two-state model: spheres and throats. This method often uses various image processing and segmentation tools to generate a network of spheres and linking throats, based on an initial tridimensional volume. Introduced by Fatt (1956), pore network modeling was first studied in conjunction with predefined network properties. Then, pore network generation was adapted to model some porous media, scanned with X-ray tomography using image-processing algorithms as accurately as possible (Dong and Blunt, 2009). Various algorithms are used to create the internal pore network structure, such as the maximal ball algorithm (Silin and Patzek, 2006). More recently, other algorithms based on the morphological properties of the studied porous media have emerged, such as the Sub-Network of the Oversegmented Watershed (SNOW) algorithm (Gostick, 2017). This alternative algorithm is considered to be computationally efficient, allowing a porous medium to be accurately modeled by numerical imagery (Khan et al., 2020). The SNOW algorithm showed a good fit with the standard maximal ball algorithm. Generating a pore network and simulating a flow in it is often cheaper, in terms of computational resources, when compared to direct numerical simulation. However, more complex transfer mechanisms, such as imbibition and drainage, are still in the study phase, and some extensive work on computational optimization has yet to be conducted, specifically on non-user-generated porous media (Maalal et al., 2021).

For each binarized Type-II and Type-III image stack, a direct pore network extraction is conducted using the SNOW algorithm implemented in the OpenPNM and PoreSpy open-source Python libraries (Gostick et al., 2019). Then, a synthetic porosity and a synthetic specific surface area may be computed for the obtained simplified representation of the pore space of the porous medium. Using the implemented Stokes equation solver, a diagonal permeability tensor is retrieved from the generated pore networks, applying the identical boundary conditions based on the method given by Sadeghi and Gostick (2020). Once again, intrinsic permeability tensors are converted into a hydraulic conductivity tensor using the relation in Eq. (8).

In Supplement S5, a comparative study is described based on Type-I samples between both developed workflows. Then, some clues are given as to whether DNS or PNM is suitable for a given sample.

3 Results

3.1 Morphological analysis

The global-porosity and open-porosity proportion (p_open) for each sample is shown in Table 3, ranging from less than 40 % (Mound1.1) to more than 95 % (Hollow2.7).

On average, lichens are the most porous of the collection, and peat is the least porous. Porosity values are in line with previously obtained data from the literature for the highest porous media of the collection (Yi et al., 2009). However, important variability can be observed for Sphagnum samples, gathering minimal and maximal porosity values. Mound mosses have an average porosity of 65.9±22.3 %, whereas the average hollow moss sample porosity is 79.6±20.2 %. Porosity profiles for each sample are presented in Fig. 4.

Figure 4 Planar porosity plot along the x, y and z axes for moss, lichen and peat samples. An averaged value is computed for each sample type, each color nuance representing each type. Type I: stable high-porosity profile samples, excluding border effects. Type II: low basal porosity, linearly increasing to the top of the sample. Type III: no specific trend observed in vertical porosity.

No specific trend can be accessed from the x- and y-axis porosity profiles, and yet variations can be observed on the z axis. Again, three trends can be observed, clustering samples into three groups according to their respective porosity profile trends.

• Type I: stable high-porosity profile samples, excluding border effects (ε_total>80 %): Mound2.6, Hollow2.7, Hollow2.8, Lichen2.1, Lichen1.3

• Type II: medium- to high-porosity profile samples associated with a progressive increase from the bottom to the top: Hollow1.2, Hollow1.4, Peat2.2, Peat2.3

• Type III: medium to low porosity associated with no specific trend porosity profiles: Mound1.1, Mound2.4, Mound2.5.

The Type-I class contains both lichen samples (Lichen2.1, Lichen1.3), whereas Type III only consists of mound Sphagnum (Mound1.1, Mound2.4, Mound2.5). The Type-II class contains half of the hollow Sphagnum samples as well as both peat samples (Peat2.2, Peat2.3).

Open and connected porosity (p_open, Table 3) represents nearly all the void space volume in each sample. Open porosity ratio values range from 0.99 to 0.9999. Thus, we can assume that, due to the fibrous nature of the studied material, enclosed porosity does not play a major role in the flow dynamics of the studied samples.

Pore size distribution (Fig. 5) is heterogeneous in each sample, and the sizes are concentrated between 0.01 and 1.00 mm of pore radii. The median pore size varies from 0.23 mm (for peat samples) up to 0.88 mm (for lichen samples).

Figure 5 Inscribed pore size distribution by classified type using particles' Feret diameter measurement. An averaged value is computed for each sample type, each color nuance representing each type.

Intermediate median pore size values can be found for mound Sphagnum samples at average values between 0.34 and 0.70 mm (for hollow Sphagnum samples). According to the previous classification, the median pore size for each sample type (I, II and III) is 0.66, 0.42 and 0.33 mm, respectively. While Type-II and Type-III curves differ in bidimensional porosity along z, they share similar global pore size distributions. Type-I samples are distinct from Type-II and Type-III curves. Specific surface area (S_SA) values for each sample are shown in Fig. 6.

Figure 6 Specific surface area (m² g⁻¹) plots for each sample. Colors refer to each sample type.

Specific surface area values seem to be uneven between each sample type. For instance, low specific surface areas can be observed for some hollow Sphagnum samples ($\mathrm{2.6}\times {\mathrm{10}}^{-\mathrm{2}}$ and $\mathrm{2.9}\times {\mathrm{10}}^{-\mathrm{2}}$ m² g⁻¹ for Hollow2.7 and Hollow2.8, respectively). Higher specific area values can be found for one mound of Sphagnum ($\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{2}}$ m² g⁻¹ for Mound2.4) and for one hollow Sphagnum sample ($\mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{1}}$ m² g⁻¹ for Hollow1.2).

3.2 Porosity

Representative elementary volumes for porosity have been computed when possible. For samples exhibiting a REV, porosity has been computed using Eq. (1) applied to the REV. For samples admitting no REV, porosity has still been computed using Eq. (1) but applied to the whole usable volume of the sample. A REV retrieval algorithm was applied to all 12 studied samples, although 2 of them (Hollow1.2 and Peat2.2) did not admit a REV. Obtained REV sizes are shown in Table 4. Some examples of tridimensional visualizations of REVs of porosity are shown in Supplement S1. Due to the numerous graphs obtained during REV computation, tridimensional porosity plots are available in Supplement S3.

REV_ε sizes vary from 2 mm to 2 cm, representing 8.0×10³ to 2.19×10⁶ voxels. Substantial morphological variations are visible, spanning from simple tubular structures (visible in the REV of Hollow2.7) to a complex and fibrous medium (for the Mound2.6 sample, Supplement S1). The average porosity obtained for these REVs (shown in Supplement S3) varies from 83.4 % to 96.0 %, which confirms the high porosity factor of these biological media. Computation times for porosity–REV retrieval range from 3 to 6 h, using two Intel^® Xeon^® E5-2680 v2 (2.80 GHz) processors and 128 GB of RAM, using high-performance Python image-processing libraries (IPSDK™). Graphical synthesis of the digital porosity assessment is presented in Fig. 7.

Figure 7 Numerical porosity estimations (%) for each sample. An averaged value is computed for each identified sample type (I, II, III) with corresponding color nuances. Peat 2.2 and Hollow1.2 did not admit any REV.

3.3 Hydraulic conductivity

Due to the time and computational resources needed to achieve a careful study of a representative elementary volume of hydraulic conductivity, only Type-I samples were studied by DNS, as they represent the most homogeneous samples of the collection. Computed REVs of the hydraulic conductivity sizes are given in Table 5. Diagonal hydraulic conductivity tensor components are shown in Fig. 8, and box plots are available in Supplement S4. Computations for the largest sub-sample size (on a 23.5 mm edge) showed higher component hydraulic conductivity values than for the three smaller sizes. This discrepancy can be related to an insufficient computation number for obtaining a good average value, hence the wider statistical spread around the mean value. Moreover, the higher values for the largest studied sizes can also be correlated with heterogeneous hydraulic conductivity behavior, as theoretically shown in Fig. 2, such as effects related to the existence of macropores. An example of a pressure field obtained on a sub-sample of Hollow2.8 through DNS is shown in Fig. 9b.

Figure 8 Diagonal components of the hydraulic conductivity tensor (m s⁻¹) on the x (K_xx), y (K_yy) and z axes (K_zz) based on DNS on a representative elementary volume of hydraulic conductivity (REV_K) for Type-I samples and with a PNM for Type-II and Type-III samples.

Figure 9 (a) Pressure field (Pa) after a single-phase flow simulation thorough a pore network model based on a Mound2.5 sample. Spherical pore sizes are represented according to their respective size in the network. (b) Pressure field lines (Pa) after a single-phase flow simulation through a sub-sample of the Hollow2.8 sample. The gray mesh corresponds to the isolated biological phase.

Table 5 Diagonal components of the hydraulic conductivity tensor (m s⁻¹) for the studied REV_K for Type-I samples using direct numerical simulations.

For three of the Type-I samples, REV_K length is computed as 15.7 mm, which is the second-largest computed size. Variations in hydraulic conductivity, with respect to study volume reduction, are smaller than those found for porosity, although study points were scarcer in the case of hydraulic conductivity assessment. Lichen1.3 shows the smallest REV_K. It can be seen that the smallest REV_ε was also described for Lichen1.3. Size differences can be seen between REV_ε and REV_K, up to 5 times larger for Lichen1.3 and half the REV_ε for Mound2.6. This seems to show that the representative elementary volume found for porosity cannot accurately describe properties such as hydraulic conductivity. This is often the case, as porosity REV is smaller than the REVs defined for other properties (Zhang et al., 2000; Costanza-Robinson et al., 2011).

Numerical estimations of hydraulic conductivity are presented in Fig. 8. For each sample of Type I, the axial components of the hydraulic conductivity tensor are given, based on the representative elementary volume of hydraulic conductivity. For Type-II and Type-III samples, hydraulic conductivity estimates are given based on pore network modeling. An example of pressure field computation is shown on sample Mound2.5 in Fig. 9a. Using a pore network allows the estimation of properties in a model based on the whole sample. The use of a pore network is an affordable alternative to direct numerical simulations at the cost of accuracy.

The values obtained vary from $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{1}}$ to $\mathrm{9.5}\times {\mathrm{10}}^{-\mathrm{1}}$ m s⁻¹ for Type-I samples and from $\mathrm{7.8}\times {\mathrm{10}}^{-\mathrm{3}}$ to $\mathrm{4.8}\times {\mathrm{10}}^{-\mathrm{1}}$ m s⁻¹ for Type-II and Type-III samples. Type-I samples can be assumed to be highly water-conductive biological media. Mean hydraulic conductivity decreases when the computed region size becomes smaller, for each direction and each sample (Supplement S4). In situ measurements, conducted by infiltration (Table 2), give an average of 10⁻⁵ m s⁻¹, which is on the same order of magnitude as previously published field measurements (Crockett et al., 2016) and computed values (McCarter and Price, 2012). Analogous values for vertical hydraulic conductivity have been found in the literature at $k_{zz}\approx {\mathrm{10}}^{-\mathrm{2}}$ m s⁻¹ (Päivänen, 1973; Crockett et al., 2016; Golubev et al., 2021). However, other studies showed results of a different order of magnitude for Sphagnum samples, with values under 10⁻⁴ m s⁻¹ (Hamamoto et al., 2016). These differences could be explained by the experimental method used to retrieve hydraulic conductivity as well as Sphagnum bog oscillation occurring during sampling (mire breathing) (Strack et al., 2009; Golubev and Whittington, 2018; Howie and Hebda, 2018), which is going to be discussed in the next part.

4 Discussion

Digital assessments of the morphological and hydraulic properties of Sphagnum and lichens of the Western Siberian Lowlands presented in this work suggest extremely porous, connected media with high specific surfaces and high hydraulic conductivities. These results are in line with the biogeochemical observations of Shirokova et al. (2021), demonstrating the overwhelming role of Sphagnum mosses in organic carbon, nutrient and inorganic solute fluxes in the Western Siberian Lowlands. Nonetheless, discrepancies between the numerical results presented in this work (Figs. 7 and 8) and previously published measurements of the hydraulic properties of Sphagnum are noteworthy (Table 2). Weaker but still sizeable differences can be seen between the results given by both of the numerical methods used here for the estimation of hydraulic conductivity on the same sample, namely, DNS and PNM. This last methodological point is discussed in Supplement S5, where a comparative validation is performed between DNS and PNM on homogeneous samples (Type I).

4.1 Numerical reconstruction after scanning

Due to technical limitations, scanning devices have a minimal resolution that causes a loss of information, acting as a threshold. In this study, minimal resolution fluctuated between 88 and 94 mum voxel⁻¹, meaning that two elements of this size could not be distinguished. In our study, technically unreachable porosity (porosity that is smaller than the minimal scanning resolution) is assumed to play a negligible role in transfers through a saturated medium, reacting as an enclosed porosity. Pre-processing algorithms (especially binarization) can cause information loss due to the arbitrary categorization of each voxel. This erroneous description can be seen for small elements (such as Sphagnum leaves) which shrink them. Mesh generation may also bring some additional “over-erosion” that helps flows inside a sample. These impacts could be studied by reducing scanning resolution, albeit not available at the time of the scans. However, hydraulic conductivity overestimation in DNS that could be related to these pre-processing effects is likely to be negligible. Indeed, the high porosities encountered and the preferential flow paths that occur in the largest pores (macropores) predominate over enclosed pore dynamics. This might not be the case for unsaturated hydraulic property assessments.

4.2 Numerical results vs. field experiments: porosity and specific surface

As described in previous sections of this study, the samples collected are considerably porous. Porosity values are in line with past results found in the literature (Yi et al., 2009; Kämäräinen et al., 2018), with porosities above 90 % for some of the samples. Interestingly, a volumetric digital-specific surface can be well linked with the porosity of complete samples as well as the average porosities found for representative elementary volumes.

A clustering can be seen for the three studied sample types (Fig. 10), although mathematical relations between specific surface and porosity are not well defined for such porous media. The specific surface values obtained are of the same magnitude as previous values obtained for other natural moss and lichen species using geometrical calculations ($\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{1}}$ m² g⁻¹ for Hypnum cupressiforme (Hedw., 1801) moss and $\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{2}}$ m² g⁻¹ for Pseudevernia furfuracea ((L.) Zopf, 1903) lichen in Adamo et al., 2007). These values are still notably lower than the values obtained using the B.E.T. method of N₂ adsorption isotherms (1.1×10¹m² g⁻¹ for artificially grown Sphagnum denticulatum (Brid., 1926) in Gonzalez et al., 2016). As discussed in Sect. 4.1, a lack of micropores could explain the observed discrepancies (1 to 2 orders of magnitude) between calculated geometry and B.E.T.

Figure 10 Specific surface (m⁻¹) as a function of total porosity ε_Total (%) and representative elementary averaged porosity ${\mathit{\epsilon }}_{\mathrm{mean}}^{\mathrm{REV}}$ (%).

4.3 Numerical results vs. field experiments: hydraulic conductivity

In Supplement S5, a comparison between direct numerical simulations and pore network modeling is made showing that pore network modeling is suitable for bypassing the heterogeneity issues observed in our samples. Indeed, the obtained porosity values with PNM are in a 5 % threshold compared to voxel-counting results (Eq. 1). The hydraulic conductivities computed by PNM and DNS are more contrasted, with 1 to 2 orders of magnitude of difference. One should bear in mind that the range of hydraulic conductivity of natural porous media is huge, with up to 15 orders of magnitude between coarse gravel (10⁻¹ m s⁻¹) and unweathered shale (10⁻¹⁵ m s⁻¹). Besides this, it is logical that the simplifications involved in the PNM method result in information loss compared to the DNS method. On the other hand, computational time savings (by using the PNM method) are huge (counted in tens of days for DNS and hours for PNM). In some cases (e.g., samples of Types II and III), DNS is simply not possible with the current regional-scale supercomputing means.

The obtained numerical hydraulic conductivities tend to show high and relatively isotropic hydraulic conductivity tensor values. Hydraulic conductivities found using DNS are sizably higher than previous values found in the literature using field percolation (Table 2), often by up to 1 to 3 orders of magnitude. The hydraulic conductivities found using pore network modeling seem to be more in line with the values in Table 2 as well as the field experiment results shown in Table 6. Nevertheless, it should be kept in mind that the results obtained by this method are less structurally accurate that those obtained from DNS, since they rely on a simplified description of the pore structure. Some clues can be advanced to explain this discrepancy, the first being the impact of numerical reconstruction routines and mesh generation procedures (discussed in Sect. 4.1), the latter being moss compression during field experiments.

Table 6 Hydraulic conductivity values (m s⁻¹) obtained using a double-ring infiltrometer during the sampling campaign.

Our digital, constant-head permeameter experiments were conducted in a fully saturated medium. Technically unreachable porosity (porosity that is smaller than the minimal scanning resolution) is assumed to play a negligible role in transfers through a saturated medium, reacting as enclosed porosity. In the case of low-permeability porous media, such sub-resolution porosity may affect flow (Soulaine et al., 2016). However, in the case of highly porous and connected media like mosses and lichens, the effects related to sub-resolution porosity are assumed to be low when compared to the effects of the large macropores, which has been shown by Baird (1997). It should also be noted that most of the porosity is opened and connected in our case.

However, moss and lichen samples are compressible (Golubev and Whittington, 2018; Howie and Hebda, 2018; Price and Whittington, 2010). Field percolation experiments induce a sizeable and rapid mass imbalance on this bryophytic cover, compacting the pore space more than would occur under natural rainfall conditions. This might notably affect flow patterns in macropores and explain the lower hydraulic conductivities found in field experiments. Indeed, some clues are given with the results of Weber et al. (2017) on hydraulic conductivity variations according to water saturation. Therefore, the numerical hydraulic conductivity assessments carried out in this study enable property quantification of the medium without perturbation, such as compression of the biological pore structure, which is not possible in field experiments.

5 Conclusions and perspectives

A numerical assessment of morphological and hydraulic properties was carried out on digital X-CT reconstructions of samples of Sphagnum moss, lichen and peat from the Western Siberian Lowlands' bryophytic cover. This porous-medium-centered approach confirmed the high porosities (from 70 % to 95 % for most samples) already found in previous studies involving experimental measurements. Hydraulic conductivity estimation was conducted using direct numerical simulations for Type-I samples and pore network modeling for Type-II and Type-III samples, both fluctuating around 10⁻¹ m s⁻¹. Indeed, both methods used in this study converge to classify macroscopic lichen, Sphagnum moss and peat as being considerably porous and pervious biological media. The hydraulic conductivity tensor shows isotropic horizontal components; however, some differences can be seen, particularly in the vertical component. Both methods reach higher values than seen before in the literature. This may have been caused by interfering phenomena, such as moss compressibility, occurring during field experiments.

The methods developed for this application show that a numerical work scheme based on image processing allows retrieval of the morphological properties of any variety of sample. Using such a method permits a nearly unlimited number of property assessments of the same sample, whereas an experimental work scheme requires many samples. Numerical methods enable a qualitative classification of the overall homogeneity of a sample, which is not easily doable using solely experimental methods. Image processing seems to be a satisfactory method, provided that the studied sample is sufficiently homogeneous for the studied property. For heterogeneous samples, image processing is not optimal. However, in the absence of another method, pore network modeling allows us to obtain some information on the studied property which is close to the one found for the homogeneous samples using image processing.

These results provide firm ground for quantitative hydrological modeling of the bryophytic cover in permafrost-dominated peatland catchments, which is crucially important for a better understanding of the global climate change impacts on arctic areas. Using numerical methods potentially enables the assessment of moss and lichen's structural hydraulic conductivity without disturbance by any biological or physical phenomena. Therefore, the porous-medium approaches developed throughout this study led to unprecedented qualitative and quantitative descriptions of such peculiar, highly porous biological media.

These physical properties can then be used as input parameters to describe ground vegetation layers in high-resolution hydrological models of arctic hydrosystems and extensively refine simulations of this critical compartment of boreal continental surfaces. For example, they will be used in further modeling studies of permafrost under climate change at the Khanymei Research Station in the framework of the HiPerBorea project (hiperborea.omp.eu). Further studies are needed to assess variable water content consequences for peat and vegetation pore structure. Indeed, water content is one of the main drivers controlling effective transport properties, such as unsaturated flow, volume change and thermal conductivity.