Solving the inverse Knudsen problem: gas diffusion in random fibrous media

About a century ago, Knudsen derived the groundbreaking theory of gas diffusion through straight pipes and holes, which since then found widespread application in innumerable fields of science and inspired the development of vacuum and related coating technologies, from academic research to numerous industrial sectors. Knudsen's theory can be straightforwardly applied to filter membranes with arrays of extended holes for example, however, for the inverse geometry arrangement, which arises when solid nanowires or fibers are arranged into porous interwoven material (like in carpets or brushes) the derivation of an analytical theory framework was still missing. In this paper, we have identified the specific geometric and thermodynamic parameters that determine the gas diffusion kinetics in arrays of randomly oriented cylinders and provide a set of analytical expressions allowing to comprehensively describe the gas transport in such structures. We confirmed analytical solutions by Monte Carlo simulations. We specify our findings for an atomic layer deposition process, the diffusion of trimethyaluminium molecules into a carbon nanotube array, but highlight the applicability of our derivation for other fields comprising gas diffusion membranes, composite materials, fuel cells and more.


Introduction
A deep understanding of the gas transport in porous media is fundamental in many fields, spanning from the development of fuel cells [1][2][3], molecular sieving [4], direct contact membrane distillation [5][6][7], geologic carbon sequestration to management of nuclear waste, and many more [8]. Moreover, knowledge of the gas transport specifics is particularly important in material science, e.g. for a uniform coating of complex porous structures from gas-phase deposition techniques, such as in atomic layer deposition (ALD) [9][10][11][12][13][14][15] chemical vapor deposition (CVD) [16] or chemical vapor infiltration (CVI), relevant to numerous sectors in industry [17]. The mechanisms of gas transport in the porous media include viscous flow, surface diffusion, molecular diffusion and Knudsen diffusion [18]. Depending on the physicochemical parameters of a system considered, some of the mentioned mechanisms can be neglected, while other ones prevail and govern the process. Particularly, the Knudsen diffusion becomes significant in the so-called molecular regime of gas transport, when the mean free path of gas molecules exceeds the characteristic dimension of the pores. This is why an accurate description of the Knudsen diffusion -and the transition to bulk diffusion -is of great importance for the case of tightly porous media.
Ever since the seminal work of Martin Knudsen [19][20][21], it has been established that modelling the gas transport in the molecular regime requires a detailed analysis of the flight paths of individual molecules. In this regime, the paths are unconstrained by the intermolecular collisions, being hindered only by the solid walls of the medium in which the diffusion occurs. Consequently, the geometry of the porous medium plays a key role in the modelling. Often, the molecular gas transport models for particular geometries are based on probabilistic numerical approaches, such as direct Monte-Carlo methods [22][23][24][25] or Markov chains [26]. Analytical expressions are available for relatively simple geometries such as cylindrical pores [19,27], rectangular cross-section trenches [28], tortuous capillaries [29] or randomly packed hard spheres [30], to name a few examples. These expressions are often used as approximations for gas diffusion in more complex media. To the best of our knowledge however, there has been no analytical model for the Knudsen diffusion in randomly-oriented fibrous structures.
In a recent work Yazdani et al. [10] modelled a case that is close to the one considered here, namely, a gas diffusion in vertically-aligned carbon nanotube (CNT) arrays; with the purpose to uniformly coat the nanostructure with ALD. Therein, a cylindrical pore approximation was employed and an expression for an effective pore diameter was proposed. This approach allowed to grasp the general trends of the diffusion coefficient with respect to the average CNT diameter and areal density. In a subsequent work [31], we introduced a diffusion model derived from basic physical principles and obtained an improved expression for the diffusion coefficient in CNT arrays with ideal verticallyaligned geometry. Here, we culminate in deriving a much more general gas diffusion model for realworld fibrous three-dimensional porous structures. It will allow to greatly expand the range of structures for which the close-form equations for gas transport are available, including electrospun fibers [32], curly CNT arrays [33], gas diffusion layers for polymer electrolyte membranes [1] and more. The previous attempts to describe the gas diffusion within random fibrous media [32,[34][35][36] relied on semi-empirical formulae or probabilistic simulations. Close-form expressions for gas transport parameters in such structures have remained elusive, until now.
The aim of this work was to derive and validate an analytical model for the inverse problem of gas diffusion in random fibrous media. We identify the specific geometrical parameters that determine the gas diffusion kinetics in arrays of randomly oriented cylinders. Expressions for the Knudsen diffusion coefficient and the Knudsen number, which determines the boundary between the Knudsen and bulk diffusion are derived. Moreover, we extend the model to express the effective diffusion coefficient in cases when both Knudsen and bulk diffusion are significant, which is relevant at relatively higher pressures [37]. Additionally, an expression for a gas impingement rate within such structures is given. We validate the analytical model by means of Monte-Carlo simulations.

Assumptions of the model
To specify the applicability of the model presented in this work, we list and discuss the model assumptions. The model presented in this work accounts for both the Knudsen diffusion and a transition to the bulk gas diffusion. While the bulk gas diffusion is already well-described and understood, initially we describe the diffusion in the purely molecular regime. The continuous extension to transition-and bulk diffusion regime is described in section 3.
In the following derivation, the gas transport within the nanostructure is assumed to occur in a molecular regime, i.e. the mean free path in the bulk gas at the gas concentrations considered is much longer than the mean flight path of a molecule in the space confined by the walls of the porous structure . It is established for an ideal gas [38], that the is expressed as where the is the Boltzmann's constant, -temperature, -gas molecule diameter, -gas pressure. The , referred to in the further part of this work as mean flight path, needs to be expressed specifically for the given porous structure geometry and in this work we are establishing it for an array of randomly-oriented cylinders.
Effectively, intermolecular collisions are neglected, only the molecule-wall collisions are accounted for. This assumption comes down to ensure that the Knudsen number Kn, as we defined below, is much higher than 1, while the molecular flow component is stated to become apparent at Kn=1 or greater [39,40]. Notably, in literature, a slip-flow regime is sometimes specified for Kn ∈ (0.001, 0.1) [1], however in this work we are limiting the division into bulk gas diffusion regime, transition regime and molecular regime, for simplicity and clarity. The Kn is conventionally defined as a ratio of to the inner diameter of the given pipe. It has been established so, because for the cylindrical pores, the is in fact equal to the pore diameter. However, in the inverse Knudsen problem that we are considering, there are no pores per se. This is why we are proposing the unambiguous definition of the Knudsen number (2), which is applicable to all kinds of porous structures, including the curly fiber bushes.
We are treating the subsequent flights of a molecule from one wall collision to another as a random walk process, consisting of the following steps:  Collision of a molecule with a cylinder wall,  Release of a molecule from the cylinder surface following the Lambert's law of emission,  Free flight of a molecule from one cylinder wall to another at a speed corresponding to the mean absolute velocity from the Maxwell-Boltzmann distribution,  Collision of a molecule with another cylinder wall.
The mean absolute velocity is expressed as where is temperature, -universal gas constant, -molar mass of the diffusant.
Due to the strictly random orientations and locations of the cylinders, the angle of emission of a molecule is immediately "forgotten" by the molecule upon release -mathematically, it has no bearing on the molecule flight distance. Notably, this is contrary to the anisotropic case of unidirectional cylinders, where the release angle is in fact the key factor determining the path length between subsequent collisions [31].
Following the Einstein-Smoluchowski model [41,42], the diffusion coefficient for a random walk is expressed as where |Δ ⃗| is a single random walk step distance, Δ -its duration, whereas -the number of dimensions in which the diffusion is considered. The triangular brackets represent an expected value of the given random quantity. Given the formula (4) it becomes clear that in order to fundamentally describe the Knudsen gas transport, it is crucial to accurately capture the statistical nature of the molecule scattering, as done for instance by Colson et al. for the case of cylindrical nanopores [43].
In general, if the geometry of the porous structure considered is anisotropic, the diffusion needs to be expressed in each independent dimension. For instance, in our previous work [31] where we considered an anisotropic diffusion within arrays of unidirectional cylinders. Therein, the 1dimensional ( =1) diffusion along the direction of fiber axis and the transverse 2-dimensional diffusion ( =2) needed to be evaluated separately. In this work however, as the considered geometry of tortuous cylinder arrays is isotropic by nature, the 3-dimentional isotropic diffusion model ( =3) is developed.
We define the as the mean time of flight of a molecule in the space confined by the structure and as the mean square flight distance of the molecule in the confined space, or mean square displacement both of which are determined by the geometry of the specific type of the porous structure.

Generalization towards a transition to a viscous regime at high pressures
Although in the derivations in this work a purely molecular regime of gas transport is considered, the results obtained are straightforwardly generalized to account for a continuous transition to a bulk diffusion regime, where Kn takes values of the order of magnitude of unity or smaller. It is important in cases where the gas pressure is relatively high, e.g. in the spatial ALD [37] or CVD [44]. The effective diffusion coefficient can be expressed in terms of the Knudsen diffusivity and the diffusivity through the porous medium in a purely viscous regime , In [37] the has been simply approximated as the bulk gas diffusivity , from the classical gas kinetics [38] being This approximation however does not account for the fact, that even in a purely viscous gas transport regime the diffusivity of gas through the porous medium is not equivalent to the bulk diffusivity [45]. It is so, because the cross section of gas flow is reduced by the factor of porosity and the concentration gradient is applied over a longer path due to tortuosity of the porous structure . Hence, the expression for the diffusion coefficient in porous structures in viscous regime is = .
The geometrical factors and are discussed in more detail in the further part of this work.

Reduced diffusivity
In the context of gas transport in porous media, the diffusivity is often expressed in a reduced form as a ratio between the diffusivity in porous medium to the diffusivity in the bulk gas [46][47][48]. Particularly, the reduced diffusivity in the viscous regime in becomes We can express the effective diffusivity accounting for both Knudsen and viscous diffusion in the dimensionless terms analogously, where the reduced diffusivity is denoted with a and an effective tortuosity is implicitly defined. If we additionally define a Knudsen tortuosity analogously as the in the equation (9), we obtain The form of the equations (11)(12)(13) shows, that the tortuosity contributions originating from different factors add up to an effective tortuosity. The physical meaning of becomes clear with the completed derivation of the in section 4.5.

Transport of gas between randomly-oriented fibers
Examples of locally randomly oriented fibers are shown in Figure 1. The interwoven CNTs ( Figure  1a) were synthesized on an alumina-coated Si wafer by catalytic CVD -the synthesis details available elsewhere [49]. Figure 1b depicts a nanoscale fibrous structure of cellulose aerogel 2 , whereas Figure  1c shows mullite microfibers 3 . Similar arrangements can be found in materials like electrospun fiber mats, recycled carbon fiber, etc. We are deriving a new set of parameters describing the kinetics of the gas transport within such structures accounting for the specifics of their geometry.   In such a case, one can define an average cylinder axes length per volume where Δ is the volume considered, Δ -the substrate surface in the representative volume, Δ is the total length of fibers in this volume, is the areal density of the cylinders on the substrate, whereas the ≔ / , by one of the common simple definitions, represents tortuosity [52]. Consequently, the axes length per volume also has a meaning of an effective areal density of the cylinders on the substrate. For the perfectly isotropic orientations, the equals which is shown in the Appendix A. Due to the random and isotropic nature of the defined system, it is safe to assume, that the tortuosity of the pores between the fibers is the same as the tortuosity of the fibers themselves, therefore from the equation (15) can be substituted into the equation (9) as .
The other key parameters describing the structure are the surface-area-to-volume ratio , = and the volume fraction of the void or porosity , The simplistic expressions (16,17) hold under an assumption that the fibers are not intersecting or that the intersections between the fibers can be neglected. There are cases however, when the intersections do occur and need to be accounted for. One of such cases is the coating of dense carbon nanotube arrays with thin films, especially if the film thickness is of the order of magnitude of the nanotube diameter or greater. The intersections alter the contribution of individual fibers to the surface area and to the void volume fraction. It is so, because the intersecting cylinder volumes count only once, whereas the intersecting cylinder surfaces do not count at all. The following expressions for the surface area to volume ratio and porosity do account for the intersecting cylinders: The derivations of (18,19) can be found in the Appendix B. On note is that, equations (18,19) converge asymptotically to equations (16,17), respectively, for small /4 as shown in the Figure  4a,b. Moreover, they are applicable for every fiber orientation distribution, such as vertically aligned cylinders, as considered in our earlier work [31], and tortuous fibers, as herein. The and are experimentally accessible in a relatively straightforward way, for instance by coupling BET surface area measurements and SEM imaging or directly by means of nano-or microtomography.
Although the expressions (16)(17)(18)(19) are not crucial for the further derivations in this work, we are providing them for the benefit of the scientific community, with a view for the experimental evaluation of the geometrical parameters in the forthcoming studies. Nonetheless, while in the analytical derivations we are assuming perfectly random orientations and positions of the cylinders without any hindrance to cylinders crossing, the equations (18) and (19) reflect the theoretically considered case.

Regular porous structures vs. inverse porous structures
The equation (18) reveals at which fiber diameter the geometry turns from an inverse-porous to a regular porous structure, e.g. in cases where a growing film infills a nanowire array, aerogel or fiber fabric to a compact nanocomposite in the CVI process. We notice, that it is characteristic for the inverse porous structures that their surface area to volume ratio increases as the film thickness increases. Examples of such structures include CNT or nanowire arrays, fiber carpets, fibrous aerogels, et cetera. The regular porous structures exhibit a decrease in with the increasing coating thickness, as schematically illustrated in the Figure 3. Examples of such porous structures include anodic alumina porous membranes (and others geometrically alike), cracks in solids, microchannels and opals, to name a few. This finding allows us to establish a clear geometrical distinction between the inverse and regular porous structures. The expression (18) reflects both the regular-and inverse behaviors and a transition between them, as seen in the Figure 4a. The maximum possible value of for the given is equal 2 / and it is reached at the = 2/( ). This finding is relevant for the applications in which maximizing the surface area to volume ratio of the porous structure is of interest. Interestingly, the = 2/( ) is also an inflection point of the (19), which means that at this critical cylinder diameter the void volume fraction begins to flatten out to converge asymptotically to 0 at high .

Isotropic fiber orientations -angular distribution
In the derivations, we are assuming that the cylinders are randomly oriented and locally straight, i.e. the radius of curvature is much larger than the individual cylinder diameter. To describe the orientations of the cylinders, we are using the isotropic angle distribution expressed in the spherical coordinate system: where and represent the azimuthal and polar angle, respectively, whereas the pdf denotes the probability density function of the quantity given in the subscript s. Due to the isotropy, the orthogonal reference frame can be defined freely at the convenience of particular derivation.

Mean molecule flight path
Analogously, as in our previous work [31], we derive the mean flight path and a mean flight time of a molecule in the given porous geometry. For simplicity, we are assuming that the size of the gas molecules is much smaller than the diameter of the cylinders. To account for the finite molecule diameter, one can use the equations obtained here, with the cylinder diameter increased by 2 , being the molecular diameter.
The molecule on a straight path Δ ⃗ can only encounter a fiber wall if the fiber axis crosses a cylindrical space of a diameter equal axially oriented along the molecule flight path Δ ⃗, which is illustrated in the Figure 5. This cylinder is referred to as molecule path cylinder. The expected length of fiber axes Δ within the volume of molecule path cylinder is expressed following the definition (14), The can be defined as the |Δ ⃗|, for which the expected number of encountered fibers is 1. In other words, the |Δ ⃗| is equal to when the Δ is equal to the average length of a single fiber axis 〈 〉 crossing the cylinder, which is expressed mathematically The geometry that is the basis for the following derivation is illustrated in Figure 6. The distance of the fiber axis to the molecule path cylinder axis is uniformly distributed, It results from the fact that from the perspective along the molecule flight path axis, the fiber presents itself as a two-dimensional stripe-shaped target of a width (see: Figure 5) and there is no preference to any particular distance . The validity of this assumption is further substantiated by a stochastic simulation in the Appendix E. Figure 6. Illustration of the geometry of the fiber axis crossing the molecule path cylinder; projected cross-sectional view along and perpendicular to the molecule path cylinder axis. The is the length of the element of the fiber axis crossing the molecule flight path cylinder projected onto the plane normal to the flight path. The is the actual length of the fiber axis element.
In the derivation of 〈 〉, it is convenient to initially obtain the projection of onto the plane perpendicular to the Δ ⃗, denoted as . Based on the geometry presented in the Figure 6, the is expressed as Given the , the angle and the probability distributions (21) and (24), we obtain the and its expected value, Ultimately, taking the equation (22), the condition (23), and substituting the equation (27), we obtain In our previous work [31], subsection Probability distribution of the transverse penetration distance, we justified that the probability distribution function pdf ( ) of the molecule flight path length follows the Beer-Lambert law, The same holds for the case considered in this work.

Mean flight time and diffusion coefficient
Knowing the distribution of the molecule flight path and taking the mean absolute velocity of a molecule as obtained from Maxwell-Boltzmann distribution, we can define the mean flight time We express the mean square displacement as set =3 in equation (4), as such isotropic diffusion is 3-dimentional, ultimately obtaining the expression for the Knudsen diffusion coefficient of randomly oriented cylinder arrays, We can define the Knudsen tortuosity in analogy to the diffusivity in viscous regime expressed with the equation (9), From the equation (33) we can see, that is proportional to the surface area enhancement of a membrane per thickness equal to and 2/π 2 is a proportionality coefficient. This statement is elaborated on in the Appendix C.

Knudsen number
Knowing the , it is straightforward to establish the Knudsen number Kn, following the definition (2),

Impingement rate
By definition, the impingement rate is the number of molecules Δ that hit the surface area Δ within the time Δ , We set the number of molecules to the average number in the given volume of void within the nanostructure, given the gas concentration (molecules per unit volume), The surface area is The Δ is effectively the mean time between collisions, while in this time each molecule should hit the surface once on average, Putting together the equations (35-38), we obtain which has the same dependency on the gas concentration and mean absolute velocity as the bulk gas impingement rate from the classical gas kinetics , differing only slightly by a numerical multiplicative factor -approx. 0.203 vs. 0.25. The impingement rate onto cylinder walls in the molecular regime needs to be described with Eq. (39), whereas the impingement rate onto a macroscopic object placed within the fiber bushes (e.g. a substrate or large particles), follows the classical expression for , as expanded on in the Appendix D.

Comparison of the analytical model with Monte-Carlo simulation results
The key parameter that fundamentally describes the gas transport kinetics is the mean flight path , while the other parameters (diffusion coefficient , impingement rate and Knudsen number Kn) are intrinsically based on its value. Therefore, we have designed a Monte Carlo (MC) simulation to verify the validity of the derivation of presented in the section 4.4. The simulation is based on generation of randomly-oriented cylinders within the specifically chosen simulation volume and measuring the distance of flight of a molecule from the release from a given cylinder wall to the collision with another cylinder.

Simulation algorithm
We have established, that we can limit the simulation volume Δ to the molecule path cylinder (see: Appendix E). The geometry for the fiber generation is illustrated in the Figure 7. We define the cylinder to be of diameter equal to the diameter of the fibers and of length equal to 10× the expected mean flight path described with the derived equation (28). Subsequently, the fiber axis lines are generated within the cylinder volume the following way. The origin point =( , , ) of each line is randomized, so that where is generated by the distribution (24) and by (20), whereas the is uniformly distributed within The is the same as for generating , and is generated by the distribution (21). A periodic boundary condition is applied to the bases of the simulation volume cylinder, which means that if the line is found to cross any of the bases of a cylinder, it is continued at the other base with the same orientation and origin shifted by ± , the sign depending which base was crossed. Examples of such random fiber generation are visualized in the Figure 8, where the cylinders of =20 nm and =10 11 cm -2 were generated within a cylinder of a larger diameter equal 3× (Figure 8a-c) and 11× ( Figure 8d) and shown only within the diameter 2× (Figure 8a-c) and 10× (Figure 8d). In the simulation we generate the cylinders in ∈(0, 10× ), however the cylinders in the Figure 8 are rendered only within ∈(0, 3× ) for illustrative purposes.
The length of the fiber axis segment crossing the molecule flight cylinder is evaluated for each generated fiber. The generation of fibers continues until a total length of fiber axes Δ = Δ is reached, as required by the equation (14).
We define a 1D grid of 0.5 nm pitch along the flight path coincident with the axis . The evaluation of the flight path length is based on finding the value of , corresponding to the coordinate of release of a molecule from the fiber surface, and finding the corresponding to the collision of a molecule with the other fiber. The flight path length Δ is evaluated as a difference between and . The exact procedure is described in more detail in the Appendix F. The entire process is repeated multiple times for each choice of and in order to obtain enough statistics of the flight path lengths.

MC Simulation results compared with analytical model
We evaluated 7 different values of ranging from 20 to 50 nm at a constant equal 6×10 10 cm -2 and 5 different values of ranging from 2×10 10 to 10×10 10 cm -2 at a constant equal 30 nm. The number of repetitions was chosen to be 10 6 . Histograms of the set of obtained values of Δ have been evaluated at 20 bins between 0 and 5× , normalized to the number of counts in the first bin. The histograms were compared to the analytically obtained probability density functions of the mean flight path (29), where the values of the function (29) were proportionally rescaled in a way for the function to return 1 at the location of the first bin. The results are presented in Figure 9.
The agreement between the analytically obtained flight path distribution and the one obtained by the stochastic simulation is excellent, which allows us to confirm that the formula for the mean flight path derived in the section 4.4 is correct, as well as the , , Kn and , all of which depend on the directly. Figure 9. The histogram data of the flight path length distribution compared to the analytically evaluated probability distribution; a) varying fiber diameter at constant fiber axes length per volume , b) varying at constant . We are emphasizing, that the solid lines are not a fitting result, but they are directly calculated using equation (29) and normalized to the height of the first bin. The number of MC repetitions 10 6 ensures that the relative uncertainties of the bin heights are negligible.

Example calculation of diffusion coefficient for atomic layer deposition on carbon nanotubes
In this section, we are evaluating the diffusion coefficient for a model case of atomic layer deposition of Al2O3 on arrays of CNTs. The diffusing species that we consider here is the trimethylaluminium (TMA) gas precursor-used typically for ALD of Al2O3 together with H2O, O2 or O3 as the oxygen precursor [53]. We set the temperature to 225 °C (498.15 K) and pressure equal to 1 mbar. At these conditions the TMA exists primarily in a monomeric form, although a significant amount of a dimer might be present as well [54]. For simplicity, we are considering only the monomer of a molar mass equal =72.09 g/mol. The diffusion coefficient has been evaluated for a range of values of from 2×10 10 to 10×10 10 cm -2 and from 7 to 30 nm, realistic for the ALD-coated CNTs. The mean absolute velocity is calculated from the Maxwell-Boltzmann distribution. The results are shown in the Figure  10a. The diffusivity exhibits a decreasing behavior both with respect to the increasing and .
The same conditions were used to calculate the effective diffusivity using the equation (7) at a range of pressures from 10 0 to 10 4 mbar, corresponding to the transition from the molecular regime to the viscous regime of gas flow. As an example, let us consider the value of =20 nm. The results of the calculations described previously are shown in Figure 10b, where the bulk diffusivity and Knudsen diffusivity top limits are indicated as well. Notably, the Knudsen diffusivity does not depend on the pressure, because in this regime the intermolecular collisions are neglected. However, at an increasing pressure, the bulk diffusivity becomes a limiting factor for the diffusion, which comes from the fact, that at higher gas concentrations the intermolecular interactions begin to play a role.
Additionally, we evaluated the characteristic time of diffusion-driven infill (alternatively: evacuation) of the structure with gas . It is estimated as a time, for which the diffusion length equals the thickness of the membrane. Consider the CNT mat of a thickness =1 mm. While the diffusion length can be estimated as ( ) 0.5 , we obtain The estimation is shown in the right axis of the graph in Figure 10b. We need to emphasize, that the is not equivalent to the time required for a conformal coating of the CNTs with ALD. For this purpose a scaling law for conformal coating of CNTs needs to be derived accounting for surface reactions, analogously as in [10]. It is however, a good estimate of the time for which a steady state of diffusion is reached within the membrane. Figure 10. Example calculations of the diffusivity of TMA in a molecular regime (a) and a transition regime (b) in arrays of tortuous carbon nanotubes for given values of nanotube diameter and fiber axes length per volume . A bulk diffusivity is plotted for reference. Additionally, a diffusion-driven infill time is indicated on the right axis of (b).

Conclusions and outlook
We have developed a new theoretical framework for the Knudsen diffusion in randomly oriented fibrous media and a transition to bulk gas diffusion regime. The model is derived from the basic physical principles, which gives an advantage over empirical or probabilistic approaches in terms of generality, simplicity and wide applicability. Moreover, the analytical expressions provided allow for straightforward development of scaling laws for various processes involving numerous gas transport processes in various fibrous media. This includes for instance, the prediction of the gas exposure required for a conformal coating of CNTs with ALD, or for a complete infiltration of fibrous structures for composites in CVI, estimation of the growth rate of CNTs synthesized by CVD in case of significant diffusion limitation of the growth, or optimizing the properties of fibrous membranes for the desired gas transport performance in given conditions, among other examples. We have provided a comprehensive set of analytical expressions for the gas diffusion parameters: Knudsen number (34), Knudsen diffusivity (32), effective diffusivity accounting for transition to bulk diffusion (7), mean flight path confined by the structure (28), mean time between subsequent molecule-wall collisions (30) and the impingement rate (39). The set of geometrical parameters determining the Knudsen diffusion is narrowed down to the mean fiber diameter and the fiber length per unit volume, both of which are accessible experimentally in a straightforward way, e.g. by tomography or by coupling microscopy and BET surface area measurements. The present model can serve as a base to generalize the description of gas transport within fibers that are quasi-unidirectional, conceptually between ideally straight (as described in our previous work [31]) and isotropically oriented, like herein. We plan to implement the presented theoretical framework into more complex processes involving e.g. pressure differential, thermal gradient, surface-or gas phase chemical reactions, physisorption, surface diffusion and thermal decomposition. These findings have a substantial impact in various fields of science and technology such as membrane science, composite and nanocomposite technologies, gas phase functionalization, thin film coating by CVD or ALD, CNT synthesis and more.
Consider the projection of the surface area Δ onto the plane inclined at an angle , Δ (see: Figure

B. Porosity and surface area to volume ratio within intersecting fiber mats
In the main part of the work we have presented the equations for and (18,19) for the case, when the fibers are allowed to intersect. The equations can be used in cases when /4 ≪ 1, which means that the average fiber diameter is much smaller than the average spacing between the fibers (effectively: intersections can be neglected), or when films are grown on the fibers, and the film thicknesses are of the same order of magnitude as the fiber diameters or greater, e.g. as grown by means of ALD, CVD or CVI.
It is convenient to begin these derivations from the porosity . Consider the structure of fibers of average diameter and axes length per volume . In the given volume Δ there is a total length Δ = Δ , as per definition (14). The following derivation is constructed as a sequential addition of small fiber lengths d into the structure, while examining the contribution to the void volume fraction (porosity) of each subsequent addition. The ∈ [0, Δ ] is the independent variable.
If there were no intersections, a single fiber fragment d would contribute d ⋅ /4 to the total fiber volume. However, only the fraction of this contribution crosses the empty space, and this fraction is determined by the current porosity ( ). We write where d ≔ d /Δ was defined for convenience. The integration of the above equation should be performed for from 0 up to its desired value, therefore the collision of symbols between the independent variable and the upper integration limit is not an issue. Given the initial condition Where d is the first contribution of the fiber element to the surface area. There is however also a second -negative -contribution. As the fiber element crosses the existing structure, the already present surface area that is crossed by the fiber is excluded. The fraction of the excluded area is equal to the volume fraction occupied by the fiber element d . We write where is the current surface area. Given the initial condition we obtain the solution of (B.7) as introduced in the section 4.1. The above derivation is easily generalized to different shapes of randomly added intersecting objects, such as spherical particles, cubes and other irregular shapes, it is however beyond the scope of this work.

C. Physical meaning of the Knudsen tortuosity
In the section 4.5 we derived the Knudsen diffusivity in arrays of random fibrous media. The expression for the Knudsen tortuosity (33) emerged and its physical meaning was outlined. In this appendix we elaborate on it. Let us consider a fibrous membrane of surface area Δ , thickness and surface area to volume ratio . The dimensionless surface area enhancement of the membrane can be calculated as Hence, can be understood as a surface area enhancement per unit thickness. Consequently, given the equation (33), it becomes clear that has a meaning of the surface area enhancement of the membrane of characteristic thickness equal to 2 /π 2 .
D. Impingement rate onto a macroscopic object within a fiber carpet As stated in the section 4.7, we must differentiate between the macroscopic impingement rate and the impingement rate in molecular regime . Specifically, in this appendix we are showing that the impingement rate within the fiber carpet onto a macroscopic object, such as substrate, is in agreement with the classical gas kinetics. The following derivation is in principle no different than the derivation of the classical bulk gas impingement rate, however we find it necessary to substantiate that the presence of fibers does not change the impingement rate onto macroscopic object. It is of note, that for the simplicity and clarity of the calculations, we treated the molecule movement as if they were always at an average absolute velocity instead of invoking the Maxwell-Boltzmann distribution. It leads however to the same result, while in the following considerations the mean absolute velocity is the determining factor of the kinetics description and the specific profile of the velocity distribution is not relevant.
Consider a very thin slice of the structure adjacent to the macroscopic object surface, illustrated schematically in Figure D  The slice is so thin, that the molecules present within this slice inevitably either hit the underlying surface or drift away from it; the collision with fiber wall is negligibly likely.
Following the definition (35), we need to establish Δ , Δ and Δ in order to obtain . The Δ is The division by 2 is present, because due to the isotropy of the gas flight directions, the movement of only half of all molecules nearby the plane Δ is directed at the plane, whereas the other half is drifting away from it. The Δ is the thickness of the slice described above. Multiplication by reflects the volume fraction available to the gas molecules.
The surface area of the object available to be reached by the gas molecules is where Δ is the element of the surface area of the object, including the parts obscured by the fibers. Taking only the gas molecules moving towards the object plane ( >0), their average drift velocity towards the object is which was calculated using the isotropic distribution of the angle , (21). For all the Δ molecules within the slice Δ to hit the surface, it will take Δ time, Consequently, according to the definition (35), the is which does remain in agreement with the classical gas kinetics. Exactly the same expression for holds at the bulk-to-fibers interface, which ensures that the gas concentration within the void of the fibrous structure is in equilibrium with the bulk gas.

E. Uniform distribution of distance of fiber axes from the flight path
In order to design the MC simulation in an optimal way, it is important to narrow the considered representative volume down to the necessary minimum while avoiding artifacts from boundary effects. In the considered case, the necessary minimum is the molecule path cylinder in which the collisions can occur, as discussed is the section 4.4. To accurately simulate the isotropically random orientations and uniformly random positions of fibers in the narrow molecule path cylinder, we have identified two solutions: a) generating the cylinders in a volume much greater than the volume of interest representing a big part of the porous structure and considering only the fibers that do cross the molecule flight path, b) generating the cylinders directly in the molecule path cylinder with carefully chosen probability distribution of positions and orientations.
The option b) is much less computationally intensive than a). In order to implement the approach b), we validated the probability distribution of the fiber axis distance from the molecule flight path (24) with the approach a). For this purpose, we considered a molecule path cylinder of 70 nm in diameter and of 1 μm length and we generated straight lines corresponding to fiber axes within the cylinder of 10× larger diameter (700 nm) and of 20× larger length (20 μm), oriented co-axially with the molecule path cylinder. Each line was generated with an origin point uniformly distributed within the large cylinder volume and the orientation vector, following the isotropic distribution (20,21). The sequential line generation was continued until we obtained 10 7 axes crossing the small molecule path cylinder. Once the generation was complete, we evaluated the histogram of distances of the lines to the axis of the molecule path cylinder and normalized the histogram so that the histogram bar heights correspond to a probability density. The Figure E.1 presents the result and a comparison to the distribution (24). As we can see in Figure E.1, the distance distribution indeed corresponds to the uniform probability distribution (24), which allows for the optimization of the MC simulation algorithm.

F. Numerical evaluation of the molecule flight path
We define a 1D grid of 0.5 nm pitch along the flight path coincident with the axis . Each element of the grid is assigned its distance to the closest fiber axis -distances to all lines are evaluated and a minimum value is chosen. We obtain dist , where is a natural number indexing the discrete points of the grid and dist is the distance value for the given . The is identified as the first , for which the distance crosses the value of /2 and is increasing with at this point, whereas the is found as the first greater than , for which the distance crosses the value /2 but is decreasing. A linear interpolation of the dist with respect to real values of is implemented for an increase in accuracy, as exemplified in the Figure   The mean value of all collected values of Δ for given and is the evaluated mean flight path .