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 (REV) 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 specific surface area . Hydraulic conductivity fluctuates by two orders of

'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 on the hydrology of high latitude continental surfaces. Therefore, some recent efforts have been put to emphasize 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 (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 & Price, 2020).
However, the mechanistic modeling of water and heat fluxes in ground vegetation layers remains difficult, as their porous media 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 3 70 75 80 85 90 95 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) 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 enhance reproducibility and inter-comparisons 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, contrary 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 patch work of biological porous media (Price et al., 2008, Voortman et al., 2014, Hamamoto et al., 2015. [ Table 2 location] 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., 1997, Koponen et al., 1996. This enables the classification of samples according to their homogeneity. Indeed, for homogeneous samples, a smaller sample region can be considered as 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 media modeling methodologies are used throughout this study: Direct Numerical Simulations on computed Representative Elementary Volumes (DNS-REV) and Pore Network Modeling with built-in solver (PNM). The impossibility to collect a substantial number of samples is compensated by a statistical quantification of a REV for each sample. This implies that the REV is smaller than the sample, hence sampling size is chosen to match sizes that were used in previous literature, such as Weber et al. (2017).

Sample collection and digital reconstruction
Samples are 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 Sampling is thoroughly conducted, to minimize structural perturbations. In order to achieve this, each sample's surroundings is 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 location] The samples are then dried at 40 °C for 48 hours 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 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 an EasyTom® XL (RX Solutions, France) with a maximal X-ray emission source set to 90 kV. The obtained resolutions, after tridimensional reconstruction, is 94 µm·voxel -1 , except for 'Lichen2.1' at 88 µm·voxel -1 (due to the scanning settings used specifically for this sample). The voxel number ranges from 4.8·10 8 voxels to 9.6·10 8 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 eight-bit black and white image stack. Supplement S1 contains some of the technical data, such as Usable Volumes and digital reconstructions, for each sample.

Drying impacts assessment on 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 hours at 40°C , at atmospheric pressure, to halt biological degradation. As expected, Sphagnum mosses began to whiten and become papery, as described by Hayward & Clymo (1982). The samples are scanned at the Toulouse Institute of Fluid Mechanics in dry conditions, two months after their primary extraction at Khanymei Research Station and drying.
To ensure the dry samples' representativity, we used an analogous drying experimental protocol than the one carried out by Kämäräinen et al. (2018). This experiment is conducted on similar Sphagnum species 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.

Morphological analysis: Total porosity (εtotal), open porosity (εopen), Specific Surface Area SSA and Pore size distribution
Global porosity (εTotal) is calculated for each sample using built-in ImageJ-Fiji's tools and macro scripting. Porosity is considered to be a ratio between the number of voxels representing the void phase ivoxel (void phase's internal volume, including closed porosity) over the total number of voxels representing a sample Ntotal (void and matrix volume). This Porosity is computed on bi-dimensional horizontal slices along the z axis to evaluate porosity variations along the samples.
Image stacks are then reconstructed along the x and y axis, to create two other image stacks. Finally, porosity is computed along the x, y and z axis using a voxel counting algorithm, shown in 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 axis, sample classification is solely based on vertical porosity z:  Type I: Constant high porosity along z axis excluding border effects;  Type II: Low basal porosity, linearly increasing to the top of the sample;  Type III: No specific trend observed on vertical porosity.
For each sample nature as well as each classified type, a dedicated color palette is choosen and kept consistent throughout the study for the sake of clarity.
Open and connected porosity (ε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): 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 SSA is obtained by dividing this surface with the sample's bounding box volume, expressed in m 2 .m -3 , as shown in Eq. (3): Specific surface area is conventionally expressed in relation to density (in m 2 .g -1 ). 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 with dry density, as shown in Eq. (4): Pore size distribution is calculated using ImageJ-Fiji's implemented image segmentation tools on the binarized image stacks.
On 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.

Darcy scale morphological and hydrological properties' definition: Representative Elementary Volume (REV)
In this study, the collected samples are assumed to form a complex fibrous porous media. 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 Representative Elementary Volume (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 as 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) 8 220 225 230 235 240 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 location] 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 as an effective property.
Finding the Representative Elementary Volumes of some key properties (e.g., porosity and intrinsic permeability) is a routine workflow in porous media 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 if compared to characteristic lengths of heterogeneities at a microscopic scale but small if 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 media property. Porosity and hydraulic conductivity related REVs are characterized throughout this study, leading to two different sizes, one for each property.

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 location] 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 deviation 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 one percent threshold can be assumed to be a high-grade REV, whereas the five-percent threshold corresponds to lower-grade REVs. For 10 of the 12 studied samples, a REV of porosity is found.
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 vi could then be injected into a regular Darcy's law, as shown in Eq. (6), where kii is a tensorial component of intrinsic permeability (m²) and µw is the dynamic viscosity: To avoid artifacts related to the physics of a specific fluid, the simpleFoam solver uses kinematic pressure (expressed in m 2 .s -2 ) and kinematic viscosity ν (in m 2 .s -1 ) to solve Navier-Stokes equations. These equations are based on intrinsic permeability k expressed in m 2 . However, in the field of hydrology the hydraulic conductivity in m.s -1 , abbreviated to Kw, is generally used. One can relate hydraulic conductivity Kw to intrinsic permeability k by using Eq. (7) described by Claisse (2016).
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.  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 showed in Table   2.

No REV for hydraulic conductivity: use of Pore Network Modeling (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 modelling 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 & Blunt, 2009, among others). Various algorithms are used to create the internal pore network structure, such as the maximal ball algorithm (Silin & 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 & 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 Direct Numerical Simulation (DNS) or Pore Network Modeling (PNM) is suitable for a given sample.

Morphological analysis
The global porosity and open porosity proportion (popen) for each sample is shown in Table 3, ranging from less than 40 % (Mound1.1) to more than 95% (Hollow2.7).
[ Table 3  [ Figure 4 location] No specific trend can be accessed from the x and y axes 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:  Open and connected porosity (popen, 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 is not playing 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 mm 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).

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 the twelve studied samples, although two 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.
[ using high performance Python image processing libraries (IPSDK™).Graphical synthesis of the digital porosity assessment is presented in Fig. 7.

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 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 a wider statistical spread around the mean value. Moreover, the higher values for the largest studied sizes can also be correlated to 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. 9-Right.
[ Table 5 location] For three of the Type I samples, REVK 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 REVK. It can be seen that the smallest REVε was also described for Lichen1.3. Size differences can be seen between REVε and REVK, up to five 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 is 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. 9-Left. 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.

[Figures 8 & 9 location]
The values obtained, vary from 1.1·10 -1 m.s -1 to 9.5·10 -1 m.s -1 for Type I samples and from 7.8·10 -3 m.s -1 to 4.8·10 -1 m.s -1 for Type II and 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 -5 m.s -1 , which is in the same order of magnitude as previously published field measurements (Crockett et al., 2016) and computed values (McCarter & Price, 2012). Analogous values for vertical hydraulic conductivity have been found in the literature at kzz ≈10 -2 m.s -1 (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 -4 m.s -1 (Hamamoto et al., 2015). 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 & Whittington, 2018;Howie & Hebda, 2018), which is going to be discussed in the next part.

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 ( Fig.7 and Fig. 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 Direct Numerical Simulation (DNS) and Pore Network Modeling (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).

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 µm.voxel -1 , 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 (macro-pores) predominate over enclosed pore dynamics. This might not be the case for unsaturated hydraulic property assessments.

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, 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 (1.4·10 -1 m 2 .g -1 for Hypnum cupressiforme (Hedw., 1801) moss and 2.4·10 -2 m 2 .g -1 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 2 adsorption isotherms (1.1·10 1 m 2 .g -1 for artificially grown Sphagnum denticulatum (Brid., 1926) in Gonzalez et al., 2016). As discussed in Section 4.1, a lack of micropores could explain the observed discrepancies (one to two orders of magnitude) between calculated geometric and B.E.T.

Numerical results vs. field experiments: hydraulic conductivity
In Supplement S5, a comparison between Direct Numerical Simulations and Pore Network Modelling is made showing that Pore Network Modelling is suitable to bypass the heterogeneity issues observed in our samples. Indeed, the obtained porosity values with PNM are in a five-percent threshold, compared to voxel counting results (Eq. (1)). The hydraulic conductivities computed by PNM and DNS are more contrasted, with one to two 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 fifteen orders of magnitude between coarse gravel (10 -1 m.s -1 ) and unweathered shale (10 -15 m.s -1 ). 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 tenth of days for DNS and hours for PNM). In some cases, (e.g., samples of Type II and III) DNS is simply not possible with the current regional scale supercomputing means.
17 505 510 515 520 525 The obtained numerical hydraulic conductivities tend to show high, and relatively isotropic, hydraulic conductivity tensor values. Hydraulic conductivities found using Direct Numerical Simulations (DNS) are sizably higher than previous values found in the literature using field percolation (Table 2), often by up to one to three 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 procedure (discussed in 4.1); the latter being moss compression during field experiments.
[ Table 6 location] Our digital, constant-head permeameter experiments were conducted in a fully saturated media. 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 & Whittington, 2018;Howie & Hebda, 2018;Price & Whittington, 2010). Field percolation experiments induce a sizeable and rapid mass imbalance on this bryophytic cover, compacting the pore space more than would occur in natural rainfall conditions. This might notably affect flow patterns in macro-pores 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.

Conclusions and perspectives
A numerical assessment of morphological and hydraulic properties was carried out on digital X-CT reconstructions of samples The methods developed for this application show that a numerical workscheme based on image processing allows retrieving the morphological properties of any variety of sample. Using such method permits nearly unlimited number of properties assessment on the same sample whereas an experimental workscheme requires many samples. Numerical methods enables 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 heteregeneous samples, image processing is not optimal. However, in the absence of another method, pore network modelling allows 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 lead 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 hydro-systems 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 on 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. Bernard-Jannin, L., Binet, S., Gogo, S., Leroy, F., Défarge, C., Jozja, N., Zocatelli, R., Perdereau, L., and Laggoun-Défarge, F.: Hydrological control of dissolved organic carbon dynamics in a rehabilitated Sphagnum-dominated peatland: a water-