Joint Distributions of Local Pore Space Properties Quantitatively Explain Simulated Air Flow Variations in Paper

The gas flow through sheet-like porous materials such as paper can show marked lateral variations due to a heterogeneous, locally varying microstructure. Hence, reliable predictions of such lateral flux variations require an appropriate consideration of local variations in the microstructure. The flow through such sheet-like materials is commonly described with Darcy’s law in which permeances are formulated in terms of microstructure properties, such as porosities, tortuosities, or hydraulic radii. This work proposes an extension of existing permeance models that directly considers the variation and the cross-dependence between local microstructure properties. The extended model is applied to local air fluxes through a paper sheet to exemplarily reveal the joint impact of local porosities and local tortuosities on the air flux. The key extension is to consider a joint distribution of porosity and tortuosity. The latter is constructed from the univariate property distributions using a copula approach and yields local tortuosities including their variation for any encountered local porosity. These values jointly enter any permeance model that qualitatively captures the dependence of the air flux on the porosity. To assess the merit of the model, variations in the air flux and in the pore space properties are independently determined from the same measured microstructure of paper. Air flux variations are provided by computational fluid dynamics simulations on multiple, nonoverlapping segments taken from the microstructure. A statistical analysis of the entire microstructure provides the distribution of local porosity, tortuosity, and thicknesses. Our model quantitatively explains that porosity-dependent variations in the tortuosity, in particular the ones associated with high-volume pathways, decisively determine air flux variations. Computational fluid dynamics simulations determine local variation of air flow through an actual microstructure of paper. Construction and use of the copula approach to model joint distribution of local porosities and tortuosities is explained. Existing permeance models equipped with variations in porosity and porosity dependent tortuosity capture most variations in air flow. Computational fluid dynamics simulations determine local variation of air flow through an actual microstructure of paper. Construction and use of the copula approach to model joint distribution of local porosities and tortuosities is explained. Existing permeance models equipped with variations in porosity and porosity dependent tortuosity capture most variations in air flow.


Introduction
It is highly desirable to quantify how the macroscopic transport through porous materials is governed by material properties that are either directly related to pores or rely on the available pore space. With such an understanding, the relative impact of pore space and external conditions on the observed fluxes can be discriminated and fluxes can be even controlled via tailored pore space properties. For gas transport through sheet-like materials, not only porosity, but also pathway properties such as tortuosity, hydraulic radii, or sheet thickness could considerably influence the permeance. Besides sheet-averaged pore properties, also the property variations due to marked heterogeneities of the in-plane composition may decisively affect the through-thickness transport. Such heterogeneities occur, for example, in non-woven fiber-based sheets such as paper due to local fiber accumulations. Local variations within the sheet can be decisive to distinguish paper types and to permit interpretations in terms of the properties of the fibers between which the pore space is formed (Sampson 2001b).
Prevalent attempts to associate permeabilities with pore space properties are predominantly lumped single pipe models (Bernabé 2018). These attempts aim at casting the macroscopic permeability into expressions reminiscent to Poiseuille's equation for viscous flow in a cylindrical pipe such that the complexity of the porous materials enters as "effective pipe" properties. Since early model equations have been introduced, back then formulated to capture the permeability for liquid flow through parallel cylinders (Carman 1939;Kozeny 1927), many different formulations of such permeability expressions have been introduced to account for the structure of highly distinct materials. Nevertheless, all proposed expressions seek to associate the permeability with a product of individual factors, each of which is governed by a pore space property.
Such lumped expressions capture the macroscopic permeance for averaged pore space properties, but are not necessarily suitable to quantitatively account for local variations (Bernabé 2018). In the present article we seek to determine, whether local variations in air permeance can be quantitatively explained by local variations in the porosity and tortuosity as the leading influences. For this purpose, we provide a reference data set of air fluxes and local pore space properties that are determined for a given microstructure of a paper sheet. This microstructure is available from X-ray computed microtomography measurements (Machado Charry et al. 2018). The through-thickness air flow and its local variations in this structure are simulated with computational fluid dynamics (CFD) simulations.
To link the air flux to local variations in porosity, tortuosity, and thickness, we assume that lumped single pipe models predict the local air flux for locally given pore space properties correctly. Then, in a central step, we extend lumped single pipe models to estimate flow variations by accounting for the dependence of tortuosity on porosities. The key extension is to incorporate this interdependence with a copula-based parametric model for the joint distribution of local porosity and local tortuosity as introduced by Neumann et al. (2021b). Note that this approach is, at the first glance, much less complex than other approaches to account for local variations. Other approaches explicitly consider the parallel transport through a distribution of pore sizes (Dodson and Sampson 2000), a pore network with variations in connectivity (Bernabé et al. 2011), and networks with varying pore sizes and connectivities (Bernabé 2018;Gostick et al. 2007). However, it appears to be reasonable to assume that a (correctly chosen) lumped single pipe model is locally valid, since we consider variations in large-area, sheet-like materials such as paper. Most of the flux variations are expected to originate from different numbers of locally accumulated fibers, such that we can expect lateral variations in the number of inlet/outlet pores as well as in the pore size distribution (Dodson and Sampson 2000). Though such fiber accumulations may also cause variations in thickness direction, these variations essentially comprise a few fiber diameters. However, in-plane variations in the pore space may occur on the mm to cm-scale (Neumann et al. 2021a), so that we can assume the local validity of a permeance model provided that the considered locations are far enough apart.
This article describes the implementation of the proposed extension in a step-by-step manner. A first step provides the data set for the analysis. Here we provide the details associated with the CFD simulations and determination of the local pore space properties. In second step, we explain the nature of lumped single pipe-inspired expressions for the permeability and show that such models readily rationalize the dependence of the simulated fluxes on local porosities. Before using this formulation for estimates, a third step is required in which the interdependence between porosity and tortuosity is quantified. In a final step, the contributions to permeability variations are collected and compared.

Material and Imaging
The sample material was a paper from industrial production, made from unbleached softwood Kraft pulp. The paper had a basis weight of 70 g/m 2 . This value, which corresponds to the supplier specifications, was confirmed by a test in accordance to the standard DIN EN ISO 536 (Paper and board: Determination of grammage).
In order to obtain the 3D microstructure of the volume, an X-ray microcomputed tomography (μ-CT) study was made using an Xradia 500 Versa 3D X-ray microscope (Zeiss, Germany). The isotropic voxel has an edge length of 1.5 µm and the final binary volume has a size of 2.001 × 2.802 mm × thickness. Detailed information about the experimental setup as well as the preprocessing and binarization of the data can be found in Machado Charry et al. (2018).
The air permeance was assessed by determining the air resistance in the Gurley method (ISO 5636-5). A reference volume 100 mL of air is pressed through the sample of an area 6.452 cm 2 at an excess pressure is p = 1.22 kPa against ambient pressure of an atmosphere p atm =101.325 kPa.
The test is recording the time necessary to push the reference volume through the paper. According to standard the result is calculated as the mean of ten measurements performed on different paper specimen. For easier comparison with simulation data we are transferring the time measured in the Gurley test into volume fluxes and, eventually, into mean air velocities.

Database for Analyzing Local Variations
The database necessary for an analysis of local variations has been built in two consecutive steps starting out from the microstructure extracted from a -CT scan of a paper sheet. In a first step, local variations in the air permeance have been determined using computational fluid dynamics simulations. We will show below that these simulations readily provide a guide how the mean air velocity, a direct measure for the permeance, is governed by the local porosity. In a second step, a statistical analysis of selected, geometry-related pore space properties has been conducted. This analysis determines local variations of these properties, finds mutual dependencies, and estimates their impact on volume fluxes of air.

Computational Fluid Dynamics Simulations
The CFD simulations based on the stationary, compressible Navier-Stokes equation were conducted in the pore space of the microstructure using the ANSYS code. The volume data was triangulized to obtain the required surface mesh of the pore space. The large aspect ratio of the sample poses a challenge for triangulation. To obtain a good comprise between a high resolution and a uniform coverage of surface triangles, we stretched the volume data in thickness direction by a factor of five, loaded this stretched volume data as images in Fiji (Schindelin et al. 2012) to perform the triangulation, and stretched the resulting surface mesh by a factor of 0.2 in thickness direction. This surface mesh encloses the pore volume in which the CFD simulations take place. As ANSYS relies on a finite volume approach, this pore volume was meshed using Numeca Hexpress Hybride. Filtering was applied. Further details can be found in Leitl (2020). The surface meshing closes all pores that are not connected to the exterior, i.e., either to an inlet or outlet. The air permeance was obtained for conditions that correspond to the standardized method test to assess the porosity of paper (ISO 5636-5:2013(ISO 5636-5: 2013.
The total cross sectional area available for simulations was A = 2.8125 × 1.666 mm 2 . However, to also obtain information on local variations, we divided the total volume into twelve equally-sized segments of A seg = 700 × 667 μm 2 . The boundary conditions for the pressure were chosen to match the Gurley measurements. The pressure at the top surface was p = 1.22 kPa, and p = 0 at the bottom surface. Also the external conditions were chosen to match the standardized conditions demanded by Gurley measurements. The temperature was T =298 K, so that the density of air was air =1.184 g cm −3 the specific heat capacity was c =1006.43 kJkg −1 K −1 , and the dynamic viscosity was = 1.838 ⋅ 10 −5 Pas.
For each segment, the CFD simulations yield the spatially resolved air pressure and air velocity, and at the outlet side the mean volume flux V∕ t and the mean air velocity v, which is given by The porosity responsible for the observed air flow is directly extracted from CFD simulation results using the pressure profile. The apparent extra effort to involve the pressure distribution (rather than the plain mesh) spares us the ambiguity to define which portion of the air-filled space contributes to the fraction of pore space and which portion is to be considered exterior air. The totally simulated volume can be viewed to consist of an inlet region, the sheet region, and the outlet region, cf. Fig. 1. These regions can be conveniently distinguished by introducing iso-pressure volumes. Each iso-pressure volume comprises all portions in the sample in which the pressure varies only in a small interval with respect to a given pressure value. After dividing the total range of observed pressures into evenly sized small intervals, we obtain a series of stacked, irregularly-shaped iso-pressure volumes, as indicated in Fig. 1. While the pressure at the inlet region adopts the value imposed by the external excess pressure, the outlet region is at ambient pressure. This choice readily distinguishes between pores inside and outside the paper sheet, i.e., it defines the surfaces of the paper sheet experienced by the air flux. The volume associated to the paper sheet is, hence, the total volume less the inlet and outlet regions. All regions inside the paper volumes that reside on a pressure level between excess and ambient pressure belong to the volume fraction of paper pores; volume regions that cannot be associated to a pressure trivially correspond to solid material. Then, the porosity is obtained straight forward from dividing the pore volume by the volume of the paper sheet.

Statistical Analysis of the 3D Microstructure
We morphologically analyze the pore space of the 3D microstructure to get more detailed insight into local variations of associated properties. In contrast to the CFD simulations shown above, such a pure-geometry-based analysis allows us to inspect local properties with different spatial resolutions, i.e., to consider a large range of segment sizes with segment lengths between 30 and 150 μm . The intended variation of segment sizes precludes the reuse of CFD-related segments and its CFD-consistent pore space. Rather, we determine the top and bottom surfaces, i.e. the boundaries of model paper in z-direction, with the rolling ball algorithm (Machado Charry et al. 2018). These surfaces are needed to distinguish air-filled pores inside the paper sheets apart from exterior air surrounding the paper.
To obtain the distributions of local properties, we partition the microstructure into sets of non-overlapping segments in the xy-plane (Neumann et al. 2021b). These segments are square-shaped in xy-plane and contain the complete microstructure in transversal z direction. The centers of the segments are arranged on a square grid in xy-plane at a distance of 150 µm to the nearest neighboring segment center. This setting leads to a maximum total number of 204 segments. The side length of the square-shaped segments is varied in the range from 30 to 150 µm in 30 µm steps; this maximum size prevents an overlap of segments.
We determined the local porosity , the mean geodesic tortuosity (0) , and the mean geodesic tortuosity (r) for pathways with a minimum diameter of 2r ( r > 0 ), i.e., of pathways with different minimum local volumes. The local porosity of a segment is defined as the ratio of the number of pore voxels to that of all voxels, which are at the same time contained in the segment and are between the top and bottom surface. The tortuosity is quantified with the mean geodesic tortuosity, because it can be directly and unambigously obtained from the 3D microstructure (Neumann et al. 2019a) and describes the connectivity between the pores, i.e., it incorporates the topology of the underlying pore network Fig. 1 Excerpt of a two-dimensional pressure profile of the CFD simulated volume including the paper sample. Regions residing at the imposed excess pressure are indicated in red (inlet volume), regions at ambient pressure are shown in blue (outlet volume). All regions at pressure between imposed and ambient pressure indicate the volume of connected pores, white regions refer to fiber volume (Clennell 1997). The mean geodesic tortuosity (0) assesses how strongly paths deviate from a straight connection of top and bottom surface. The more a path meanders the more the associated (0) would exceed the minimum value of one. The mean geodesic tortuosity describes the length of the shortest connection that can possibly be formed between two points through a succession of pores. As such connections are permitted to find any shortcut, geodesic paths are shorter than those considered by geometric tortuosities even if the same pores are visited, as in the latter case paths are constructed on the pore skeleton. Tortuosities derived from actual flow patterns, being obtained from diffusion or hydraulic tortuosities, give even larger values as the flow does not always follow the geometrically shortest option (Clennell 1997).
Within each segment, (0) is taken as the average of shortest path lengths from the bottom surface to the top surface divided by the local thickness of the paper sheet (Machado Charry et al. 2018;Neumann et al. 2019a). For a formal definition of mean geodesic tortuosity in the framework of random closed sets, we refer to Neumann et al. (2019a). From all pathways present in a segment, we only consider paths, whose starting points are located in the considered segment. These paths, however, are allowed to leave the segment, because otherwise boundary effects would have a stronger impact on (0) of smaller segments. The lengths of the shortest pathways associated with the starting points in a segment are determined by means of the Dijkstra algorithm (Thulasiraman and Swamy 1992) on the voxel grid.
Pathways differ not only in length, but also in their volume. To discriminate between pathways in terms of their volume, we consider the mean geodesic tortuosity (r) of pathways with a minimum radius r, where (r) characterizes the connected pore space that can be filled by spheres of radius r. Hence, radius r serves as an indicator for the path volume. The larger r, the weaker the bottlenecks along the path, i.e., the larger the volume of the path. For the selection of pathways with a minimum radius equal or larger than r, the value of (r) is determined in the same way as described above for (0) . We consider the cases r > 1.5 and r > 3.0 μm , so that particularly the case r > 3.0 μm will provide insight into the behavior of high volume pathways later on. Figure 2 shows that the mean local air velocity, determined for each segment, tends to increase with local porosity. However, the velocities exhibit a marked scattering beyond their porosity dependence. This is seen best for the segments that share a common porosity of ca. 0.32, where air velocities scatter by ±50% (highlighted region in Fig. 2).

Air Velocity Depends on Porosity
The strong relation between mean local air velocities and local porosities reproduces a key observation for porous materials. In general, this relation triggers ongoing efforts to quantify how permeabilities depend on porosities and other properties of the microstructure. Such quantifications are commonly performed in the spirit of lumped single pipe models (Kozeny 1927;Carman 1939;Bernabé 2018). To get a rapid overview, how the microstructure can influence the air velocities, we collect the most important steps to represent the velocities with lumped single pipe models. As an instructive starting point, we assume that Darcys law holds, so that the volume flux V t relates to the pressure gradient, ∇p , and the dynamic viscosity, , as where is the permeability. Darcy flow is appropriate because it holds for laminar pore flow. Finding a maximum air speed of about 5 ms −1 in the paper pores during the Gurley test (compare Fig. 4 below) and a maximum pore size in the paper of 10 μm , we arrive at a Reynolds number Re of circa three, which is way below the critical value for turbulent flow, Re = 2300.
We further assume a practically constant pressure gradient ∇p that can be approximated by the global pressure difference Δp and the thickness, h, of the penetrated volume, i.e., we assume that Then, a mean velocity v is readily obtained by dividing both sides of Eq. (2) by the nominal area A, i.e., This mean velocity is not only closely associated with permeability , but also depends on the thickness h of the sheet. Since sheet thicknesses h are much smaller than the lateral sheet extensions, the actual value of h could decisively influence the value of v. For example, flows transmitted to a thin sheet experience much less of variation in through-thickness direction compared to the variations in lateral directions. For this reason, we prefer to characterize the flow by a segment resolved mean velocity v rather than by local permeabilities .
In a next step it is desirable to refine the permeability in terms of porosity and other microstructure-derived properties of the sample. The widely used Carman-Kozeny equation illustrates such a refinement for a laminar liquid flow through packed beds: where C CK is a material-specific factor that contains information on the structure of the pore space (Kruczek 2015). When considering, for example, packed beds as an assembly of parallel, cylindrically-shaped pores formed between identical particles, C CK reads as where d p is the diameter of the packed particles, A s the specific surface area per unit volume of the particle, and K an empirical parameter whose value depends on the tortuosity (Kruczek 2015).
A possible way to explicitly account for path lengths which significantly exceed the thickness h of a sheet or a membrane is to express the pressure gradient with the average path length h , i.e., −∇p = Δp∕( h) (Kruczek 2015). This leads to Even when going beyond packed-beds, expressions for v via Darcy's law essentially keep the same structure. The permeability is a product containing a factor accounting for the characteristic size of the involved solid, a factor accounting for the porosity, and a factor related to tortuosity. Models for v and the permeability of such a structure were also proposed for the case of fibrous materials (Jackson and James 1986;Gebart 1992;Koponen et al. 2017Koponen et al. , 1998Clague et al. 2000;Nabovati et al. 2009). Often the fiber diameter, a, serves as characteristic size so that typically scaled permeabilities ∕a 2 are formulated. The model of Gebart may serve here as an example: where C G is a geometric factor, and c is the critical value of porosity below which there is no permeation, i.e., c is a percolation threshold (Nabovati et al. 2009). Both factors depend on the specific packing of fibers and do not account for fiber diameter, , or . Values of an ideal square and hexagonal array of parallelly oriented fibers ( = 1 ) are C G = 16∕(9 , respectively (Gebart 1992). This model has been later refined by Nabovati et al. (2009) and Clague et al. (2000) to cover an even larger porosity range.
Knowing the velocity v, the porosity , and the thickness h of each segment, models of the type given by Eq. (7) or Eq. (8) readily permit to check the consistency of the simulated v CFD with the predicted v( , , h) dependence and to fit associated material-dependent parameters (including the characteristic size). Regardless whether Eq. (7) or Eq. (8) is used to model the v( ) dependence, it is possible to obtain an equal-quality fit in the given range of local porosities, cf. Fig. 2. Even though the trend in is well reproduced, our data set is not suited to verify one or the other model for v( ) . The true value of lumped single-pipe models is that their typical structure readily suggests how to search for other relevant properties of the pore space. This will be explained in the next section.

How to Assess the Roles of Thickness and Tortuosity
In our case, we are interested in the impact of the thickness h and the tortuousity . The thickness is an obvious candidate to provoke local variations in v, because its value depends on the number and shape of the locally stacked fibers and, hence, position in the paper sheet.
But also the tortuosity appears to be an intriguing candidate, as it parameterizes the nature of pathways.
The streamlines of air flow inside a simulated segment, shown in Fig. 3, reveal that the air flux follows extended path segments in lateral direction, rather than straightly penetrating the paper in thickness direction. In this work, we chose to represent the pathways by means of mean local geodesic tortuosity (cf. Sect. 3.2). This particular definition of tortuosity solely accounts for the morphology of the connected pore space (Clennell 1997).
Given that there is a formula for v model ( , , h) , like in the Carman-Kozeny model, we seek to estimate the variation that contains all dependencies in addition to the influence of the porosity seen in Fig. 2. In particular, our velocity model suggests to which extent v varies with thickness and tortuosity. Then, provided that we know the local thickness h, remaining variations in v can only be caused by variations in the local tortuosity . Hence, we need to provide information how strongly tortuosities could vary for a given porosity. The required relation between tortuosity and porosity along with the entailing variation in the tortuosity will be provided in the next section. Before doing that, however, it is useful to have a closer look on the pathways involved in transport.
The CFD simulations strongly suggest that air flows in heavily preferred pathways. Figure 4 illustrates this finding with a cross-section through a segment that is centered between the sheet surfaces and oriented perpendicular to the imposed pressure difference, cf. Fig. 4a. All gray-shaded regions correspond to pores, white regions are fibers. For each pore area, the absolute value of the local velocity, being a measure for the local volume flux V = V∕ t , has been recorded. In Fig. 4a, all pores with a volume flux between the maximum flux, V max , and 0.5V max are highlighted with a color corresponding to the local velocity. Clearly, only a small fraction of pores governs the overall flow. To quantify this fraction, we (i) arranged the pore areas in descending order with respect to the permitted volume fluxes V and (ii) show in Fig. 4b how the volume flux increases by incorporating more and more of such pore areas. Note that the flux accumulated this Fig. 3 a Location of the streamlines seen in top view on the lateral xy-plane for one of the simulated segments. The color bar indicates the location of a streamline segment in transversal direction. b Locations of large lateral velocity components in an xy-cross-section at the center of the segment way appears to exceed the total flux (100% in Fig. 4b); this excess originates from the fact that the total net flux also contains contributions flowing against the external pressure gradient. Figure 4b suggests a marked disparity between relevant and less relevant pore crosssections. Half of the flux is transmitted through only 6% of the total pore area. In fact, an equivalent of less than 20% of the total pore area is needed to account for 90% of the flux.
A quick and convenient way to account for the disparity between offered and actually used pathways in our analysis is explore the role of high-volume pathways using (3.0) . Hence, it is worthwhile to clarify both the relation between the porosity and all pathways, (0) ( ) , and the relation between porosity and high-volume pathways, (3.0) ( ).

Capturing the Relation Between Porosity and Tortuosity
In our effort to estimate Δv due to variations Δ in the tortuosity, we need to know the extent to which the tortuosity varies with porosity. Formally speaking, we require the change of with , ∕ , and the -dependent variation in for each local porosity. In the absence of a known analytical formula for ( ) , we inspect our microstructure to obtain such a relationship empirically.
The key idea is to consider a bivariate porosity-tortuosity probability distribution function K( , ) to model the joint evolution of and and their local variations across the material. From this model of K( , ) , a conditional probability density function k | = � ( ) can be derived to encounter a certain tortuosity for a given porosity ′ . The conditional Fig. 4 a Positions at which velocities at least as large as half of the peak velocity are encountered (highlighted in color) in the cross-section shown in Fig. 3b. The gray regions mark pores, white regions correspond to fibers. b Volume flux vs pore cross-section area. Starting from the pores carrying the largest volume flux, the pore area successively increases by incorporating the pores carrying lesser volume flux (in descending order). Only a small fraction of the pore area (ca. 5%) is responsible for 50% of the total flux. c Derivative of the volume flux with respect to pore area. The region highlighted in gray refers to pores showing backflow, i.e., with a flow against the externally imposed pressure gradient probability density k | = � ( ) provides us two important quantities: The first is the (conditional) distribution of tortuosity values for a given porosity ′ from which the variation in directly follows. The second is a relation between an averaged tortuosity, ̄ , and , from which we can estimate ∕ .
The joint model distribution function K( , ) is directly constructed from available geometry data using the copula approach as formally described in Neumann et al. (2021b). In a step-by-step manner, we now briefly motivate, sketch, and illustrate the actions that are necessary to provide the joint model distribution function K( , ) and to derive k | = � ( ) . For a complete description of all technical details the reader is referred to Neumann et al. (2021b).

Joint Distributions and Copulas
The possibly best known example for joint distributions are bivariate normal distributions that relate a pair of normally distributed random variables to each other. In this case, a variable X with mean X and standard deviation X and a variable Y with mean Y and standard deviation Y possess a joint probability density function with x, y ∈ ℝ , and ∈ (−1, 1) being the correlation coefficient. Expression (9) illustrates the two types of necessary ingredients: There are the two individual (marginal) distributions, each parameterized by its mean and its standard deviation . Then, solely one additional parameter is necessary to judge to which extent variable Y (co-) varies with X. Such bivariate normal distributions have already been successfully used to distinguish paper sheets in terms of characteristic relations between thickness and basis weight (Dodson et al. 2001).
In principle, the bivariate distribution of local porosities and local tortuosities of our material is readily sampled (i) by collecting the value pairs ( , ) in non-overlapping inspection regions of a given size and (ii) by determining how often certain value pairs ( , ) occur across the sample. This yields a sampling-dependent point cloud in the − space. To interpolate to situations located between the probed , values, a continuous model function is called for. Thus, we seek for a bivariate probability density k( , ) , alike the bivariate normal density g(x, y) given in Eq. (9), that reproduces the empirically obtained value pairs ( , ) and contains the relation between and as a parameter (alike in Eq. (9)). This is where copulas come into play.
Copulas combine three aspects in a highly convenient way. First, the use of copulas yields a parametric model for the bivariate distribution function K( , ) . Second, a copulabased K uses the univariate distributions of the individual quantities and (i.e., the marginal distributions) by construction, regardless whether each of them is normally distributed or not. Third, these univariate distributions are joined into K using a single additional parameter that informs on the correlation of the two quantities. Formally, a two-dimensional copula C ∶ [0, 1] 2 → [0, 1] is a particular kind of a bivariate probability distribution function. It starts out from an arbitrary pair of random variables X (modeling, e.g., sampled values of ) and Y (modeling, e.g., sampled values of ) with their joint (bivariate) distribution function K(x, y), where K(x, y) = ℙ(X ≤ x, Y ≤ y) , and their univariate distribution , respectively. Each of these cumulative distribution functions, F(x) and G(y), is associated to a probability density function, f(x) and g(y), that describes how often a given value x or y, respectively, is encountered in our material. Sklar's theorem (Nelsen 2007), the main result of copula calculus, states that there exists a copula C which admits the representation of the joint distribution function K ∶ ℝ 2 → [0, 1] of X and Y, for all x, y ∈ ℝ. Therein, the copula C is the joint (bivariate) distribution function a random vector (U, V) containing two uniformly distributed variables U and V (becoming F(X) and G(Y) in our context) with values in the unit interval [0, 1]. Thus, for all 0 i.e., both marginals U = F(X) and V = G(Y) are uniformly distributed on the unit interval (Schweizer 1991).
The copula function C to be employed in Eq. (10) can be chosen from different function families. Within such a family, the functions C are distinguished with a parameter ∈ ℝ whose value measures the degree of correlation in analogy to the correlation coefficient in the bivariate normal distribution (Eq. (9)). In other words, the cumulative joint distribution function K(x, y) of the random vector (X, Y) is expressed by a copula function C that takes the values of the (random) cumulative distribution functions F(X) and G(Y) as arguments and contains a parameter to quantify the correlation between X and Y.
In our case, X and Y denote the random local porosity and the random local tortuosity, respectively, i.e., If the value pairs ( , ) have a joint probability density function k ∶ ℝ 2 → [0, ∞) , we directly obtain a formula for k from Eq. (11), namely where f ∶ ℝ → [0, ∞) and g ∶ ℝ → [0, ∞) are the univariate probability density functions corresponding to F and G, respectively.
To arrive at expressions of the form (11) and (12), three steps are required: (i) determine the univariate probability density functions f ( ) , g( ) and the corresponding cumulative probability distribution functions F( ) , G( ) , (ii) select a candidate copula function C , and (iii) fit the bivariate density function k( , ) , by determining the parameter of C derived from the " − point cloud".

Local Variations in Porosity and Tortuosity
We estimate the individual (univariate) distributions of local porosity and local tortuosity, respectively, directly from the 3D microstructure for each size of the segments. In Fig. 5a, b, the resulting probability densities are shown for the segment size of 150 µm.
In a second step, we will cast these univariate distributions in parametric models. Such parametric models provide a lean and feasible way to approximate the probability densities estimated from the 3D microstructure by analytical expressions. The few shape parameters of the analytical model function suffice to reconstruct the entire distribution, including the mean value and standard deviation.
The distributions associated to the local porosity, shown for the segment size of 150 µm in Fig. 5a, can be well modeled by beta distributions (Fig. 5c). Note that the probability density f s,q ∶ ℝ → [0, ∞) of the beta distribution is given by for each ∈ ℝ , where s, q > 0 are some shape parameters, 1 [0,1] is the indicator function of the unit interval [0, 1], and Γ denotes the gamma function (Johnson et al. 1995). Though the beta distribution appears to closely resemble a normal distribution (when s = q), the beta distribution is more appropriate to model the distribution of local porosities, as it requires the sample values to vary within the interval [0, 1]. The shape parameters s and q are fitted by the maximum likelihood method. In Fig. 5c, the probability density function of the fitted beta distribution is shown for the segment size of 150 µm. Comparing Fig. 5a, c, it becomes clear that the model density approximates the corresponding empirical density estimated from image data very well. Moreover, the fitted parametric model densities nicely resemble the estimated densities of local porosities across all inspected segment sizes (Neumann et al. 2021b).
The distributions of (0) and (3.0) , estimated from image data, are illustrated in Fig. 5b for the segment size of 150 µm. The average high-volume pathway is considerably longer than the overall average pathway (Neumann et al. 2019b). Both distributions are markedly skewed. Regardless of the segment size, we found that the estimated distributions are well modeled by a shifted gamma distribution (Neumann et al. 2021b) with probability density function g a,b ∶ ℝ → [0, ∞) defined via for all ∈ ℝ , where a, b > 0 are the parameters of the distribution (Johnson et al. 1995). For the segment size of 150 µm, the probability density function of the fitted shifted gamma distribution is shown in Fig. 5d.
Note that also other morphological attributes of the pore space turn out to be non-normally distributed. For example, the pore size distribution being another important pore space property follows a markedly asymmetric gamma distribution (Johnston 1998;Sampson 2001a;Dodson et al. 2003).

Correlation Between Porosity and Tortuosity
The segment-wise collected pairs of and (0) (resp. and (3.0) ) are shown as scatter plots in Fig. 6a (resp. Fig. 6b) for the segment size of 150 µm. Asking how often specific porosity -tortuosity pairs occur transforms the scattered data points into an empirical joint probability density function associated with our model paper. Note that for each segment size considered in this study, another (i.e. different) empirical bivariate probability distribution is obtained. However, for the segment size of 150 µm, it turns out that the values of the parameter of the fitted Gumbel copula are almost identical for (0) and (3.0) . In the following we describe how the empirical bivariate probability distribution function, estimated from the value pairs ( , ) , can be approximated by a parametric bivariate distribution function K ∶ ℝ 2 → [0, 1] , picking an appropriately chosen copula function C ∶ [0, 1] 2 → [0, 1] . In general, one is free to choose the copula function that accounts best for the data at hand. As a large number of possible copulas C exists (Nelsen 2007), it is by no means straight-forward to pick the most suitable one. Rather, the identification of robust criteria to appropriately choose copulas appears to be an active area of current research. Yet, it is possible to line out aspects that can potentially narrow possible choices. Copulas differ in the way data in the tail regions of the univariate distributions is considered. Tail regions comprise data values that are located most remote from the mean value, i.e., at which the cumulative distribution function takes values either very close to zero or unity. Compared to the possibly most employed copula, the Gauss copula related to bivariate normal distributions, the Gumbel copula emphasizes more arguments close to zero, while a Clayton copula provides a larger emphasis on arguments close to one. In our context we Fig. 6 Estimation and parameterization of the joint distribution of local porosity and local mean geodesic tortuosity of a all pathways (0) and b high-volume pathways (3.0) for a segment size of 150 µm. Joint distribution of the two properties as directly obtained from the segments of the microstructure (scatter plots). Parametric model of the joint probability density based on the Gumbel copula with parameter = 1.73 for (0) and = 1.74 for (3.0) (heatmaps). The triangles at the porosity axes indicate specific values of local porosity for which the conditional distribution of the local mean geodesic tortuosity is investigated later on in Sect. 5.4 had no reason to pay particular attention to this tail behavior. Thus, we compared the fits obtained with one-parametric families of so-called Archimedean copulas (Nelsen 2007), i.e., the fits with the families of Gumbel and Clayton copulas to fits with copulas of the Frank type. This comparison revealed that Gumbel copulas achieve a description that is superior to Clayton and Frank copulas, even though the Frank copula type was already shown to appropriately model the correlation between local porosity and local mean geodesic tortuosity in one of our previous papers (Neumann et al. 2019b).
The Gumbel copula C ∶ [0, 1] 2 → [0, 1] with parameter ∈ ℝ ⧵ {0} is given by for all 0 ≤ u, v ≤ 1 . A detailed account on the properties of the Gumbel copula is provided e.g. in Genest (1987). Important is the meaning of the parameter . In the case of a Gumbel copula, the parameter is closely related to Kendall's rank coefficient K via the expression i.e., the value of is closely related with the result of a statistical correlation test. Namely, if the parameter tends to unity, meaning that K tends to zero, then Eq. (15) implies that the Gumbel copula C tends to the product copula C product given by C product (u, v) = uv for all 0 ≤ u, v ≤ 1 , which is equivalent with the independence of the random variables X and Y introduced above in Sect. 5.1. On the other hand, the larger the value of is the stronger is the correlation between X and Y.
The parameter of the Gumbel copula is usually determined via pseudo-likelihood estimation (Genest et al. 1995). To assess the goodness-of-fit, a formal validation can be performed by means of a multivariate two-sample test (Baringhaus and Franz 2004;Neumann et al. 2021b). The values of obtained in this way are similar to those resulting from Eq. (16) via the estimation of Kendall's tau K (cf. Supporting information).
For the joint distribution of local porosity and local mean geodesic tortuosity (0) of all pathways, the value of = 1.73 was obtained for the segment size of 150 μm . The corresponding probability density is shown as heatmap in Fig. 6a. Note that for increasing local porosities, the local mean geodesic tortuosities shift towards smaller values and are more narrowly distributed. For comparison, we visualize the joint probability density of and (3.0) in Fig. 6b for the same segment size of 150 μm . As (3.0) refers to the tortuosity of only such pathways, whose minimum radius exceeds the threshold of 3 μm , the value of (3.0) represents the particular behavior of pathways that allow to transmit high volumes of air. Both the empirical and parameterized joint distributions of and (3.0) were obtained following the same steps as lined out above for and (0) . Despite a much larger variability of (3.0) , the estimated value of = 1.74 indicates that the interdependence between (3.0) and is very similar to the one between (0) and , where we obtained that = 1.73.

Conditional Distribution of Local Tortuosity
With the copula model for the joint probability density k( , ) of local porosity and local mean geodesic tortuosity established, we can easily determine the conditional probability density function k | = � ( ) of local mean geodesic tortuosity for a given value = � ∈ (0, 1) of local porosity , and vice versa. Namely, from the definition of k | = � ( ) , which is given by the quotient k | = � ( ) = k( � , )∕f ( � ) for all > 1 , using Eq. (12) we obtain that In Fig. 7, the conditional probability density functions k | = � ( ) are shown for the tortuosities (0) of all pathways and (3,0) of high-volume pathways, and for the selected porosity values ′ that were indicated by triangles in Fig. 6.
The (conditional) probability densities of the tortuosity (0) of all pathways, shown in Fig. 7a for porosity values ′ between 0.2 and 0.5 markedly differ in shape. For a porosity of 0.5, the density of (0) (red curve) is sharp and possesses a maximum close to unity. At paper locations with such an elevated porosity, there are ample opportunities to form paths that straightly connect top and bottom surfaces. Moreover, the pore space is large enough that very likely pathways form that are almost as short as the local film thickness. However, at locations with porosities as small as 0.2 (blue curve) only little pore space is available, be that due to scarcely available, small, or ill-connected pores. Correspondingly lesser and longer pathways are formed in comparison to a large porosity. Thus, the actual local situation decides whether a formed pathway is short or long in comparison to the film thickness. These associated large local variations are directly reflected in the wide distribution of (0) (blue curve).
The conditional probability densities of the tortuosity (3.0) of high-volume pathways, shown in Fig. 7b for porosities ′ between 0.3 and 0.45, illustrate the above-discussed ′ -dependence even more extremely. With the constraint of a minimum path radius of 3 μm , much less pathways are formed. Chances to find suitable pores to connect in the thickness direction rapidly decrease with increasing porosity. Hence, the mean geodesic tortuosity is massively shifted towards higher values and shows a wide distribution (note the different scales of the horizontal axes in Fig. 7a, b).
(17) k | = � ( ) = g( ) 2 C (F( � ), G( )).  Mean values ̄( 0) ( ) and ̄( 3.0) ( ) of the conditional distribution of (0) and (3.0) , respectively, as function of the local porosity (curves in the middle of shaded areas). The vertical extent of the shading at each porosity value indicates twice the size of the standard deviation of the corresponding distribution of local tortuosity performed on segments large enough to suppress boundary effects. However, the extent of variations always depends on the segment size, i.e., the size of the inspection region. The larger the inspection region, the smaller the variation in pore space properties (Sampson 2001b;Kanit et al. 2003;Rolland du Roscoat et al. 2007;Neumann et al. 2021a, b) and associated fluxes (Dodson and Sampson 2000). Considering that, is our model going to predict substantially smaller variations in v when we assume smaller local variations in and consistent with the segment size? For a similar paper type (Neumann et al. 2021a) we already know that porosity variations become independent of the inspection region size beyond a segment length of 120 μm (cf. also Rolland du Roscoat et al. (2007)), while tortuosity variations (both in (0) and (3.0) ) continue to shrink even for segment sizes well above 300 μm . To quantitatively grasp the consequence of reduced variations in and , we performed the analysis for v( , (0) , h) as described in Sect. 5 for a smaller segment size of 90 μm with correspondingly wider porosity and a (0) distributions. For this segment size, both univariate distributions, their standard deviations (Neumann et al. 2021b), and the Gumbel coefficient = 1.59 of the joint distribution are available in the Supporting information, Fig. 3. Fitting the correspondingly adjusted model to v CFD yielded the same Gebart coefficients C G a 2 Δp∕ as found for segment size 150 μm . Hence, the predicted mean velocities coincide with values shown in Fig. 9 and the velocity spread, Δv , can be directly compared between the segment sizes. However, the standard deviation v of v G , serving as a first, conservative estimate for Δv , reduces by one third when going from a segment size of 90 μm to 150 μm . Taking instead Δv as the range of predicted velocity values v G between the 0.95 and 0.05 quantiles, as was done in Fig. 9, leads to a reduction of Δv by 25% ( cf. Fig. 4 in the Supporting information for details). Regardless which estimate is used for Δv , the variations in v G found for the two segment sizes remain in the same order of magnitude. Hence, local variations of the tortuosity remain a plausible cause for the variations seen among the CFD-simulated velocities.

Summary
We quantified the impact of variations of two correlated, local microstructure properties on locally varying mean velocities of transmitted air through a sheet of paper. As a model system, we considered the simulated air flux of air through the microstructure of a paper sheet as obtained by -CT. The flow simulations readily corroborate the strong dependence of the air velocity on the porosity, but show additional variations. To also relate these additional variations to the microstructure, we analyzed the impact of thickness and the tortuosity representing pathway-related properties, because the stream lines obtained in the CFD simulations suggest that fairly long pathways are formed and that most of the air is transported in high-volume pathways. Our quantification approach combines lumped single pipe-inspired permeability models (here, the Gebart model for fibrous materials), which explicitly contain the microstructure parameters, with an empirical, parametric model to account for the correlation between porosity and tortuosity. This parametric model represents the joint porosity-tortuosity distribution that relates the porosity, the tortuosity, and their local variations obtained from the microstructure. Key to construct such a joint distribution is the copula approach, where the copula joins the univariate distributions of porosity and tortuosity, both obtained from the local microstructure, with a correlation parameter. The resulting parametric model successfully approximates the distribution of local porosity and tortuosity value pairs found in our microstructure and permits the extrapolation to continuously varying property values. While employed here for two property pairs, i.e., the tortuosities related to two classes of pathways and the porosity, the copula approach is suitable for any conceivable pair of microstructure properties. With extending permeability models of the lumped single pipe type with the porosity dependent local variations in tortuosity, we quantitatively predict local variations in mean velocities that can plausibly explain the simulated air fluxes.
Author Contributions KZ, UH, MN, and PL contributed to the study conception and design. Collection, postprocessing of the microstructure data was performed by EM, EB, and MN. MN and VSt performed the statistical microstructure analysis and conceived the copula model. PL performed the computational fluid dynamic simulations. KZ analyzed the simulated fluxes using permeance models combined with copulabased predictions of local variations. The first draft of the manuscript was written by KZ and all authors commented each version of the manuscript. All authors read and approved the final manuscript.
Funding Open access funding provided by Graz University of Technology. The authors gratefully acknowledge financial support from the Christian Doppler Research Association, Federal Ministry for Digital and Economic Affairs and the National Foundation for Research, Technology, and Development, Austria. MN acknowledges funding by the German Research Foundation (DFG) under Project ID 390874152 (POLiS Cluster of Excellence, EXC 2154). KZ thanks the "TU Graz Lead Project LP-03: Porous Materials @ Work for Sustainability" for inspiration.
Data Availability Additional data is available as supporting material. Additional datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Competing interests
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.