Method for obtaining the Knudsen diffusion coefficient

Graphical abstract


Theoretical basis
In this section, we first present the basic theory of the DGM as it applies to a binary gas system in first section. Then we explain in second section how we derived advection-diffusion equations for calculating the effective compound diffusion coefficient and the effective compound velocity of each component gas from the DGM equations. In final section, we explain how the Knudsen diffusion coefficient of each component gas is calculated using the derived advection-diffusion equations.
Basic theory of dusty gas model for a binary gas system The dusty gas model (DGM), which was developed by Mason et al. [1], can take account of the molecular diffusion flux, the Knudsen diffusion flux, and nonequimolar fluxes, and it can solve the migration of gases in multicomponent gas mixtures in a porous medium. If a gas flowing horizontally in a column filled with dry soil and a different gas is assumed, then the effect of gravity can be omitted from the equations for the migration of the first gas in the porous medium. In a binary gas system, the DGM equations are as follows: where X A and X B are the molar fractions (dimensionless), N D A and N D B are the molar fluxes (mol/L 2 T), and P A and P B are the partial pressures (M/LT 2 ) of component gases A and B, respectively; R is the gas constant (ML 2 /molKT 2 ); T is absolute temperature (K); D AB is the molecular diffusion coefficient between components A and B (L 2 /T); D A and D B are the Knudsen diffusion coefficients (L 2 /T) of components A and B, respectively; t m and t p are the tortuosity for molecular diffusion and Knudsen diffusion, respectively; and u g is the gas-filled porosity [2][3][4]. t m (0 < t m 1.0) is the ratio of the molecular diffusion coefficient of the porous medium to the molecular diffusion coefficient in the absence of any obstruction by solid particles. t p (0 < t p 1.0) is the ratio of the Knudsen diffusion coefficient in the DGM modified for a real porous medium to the Knudsen diffusion coefficient in the original DGM, which was developed by assuming an ideal parallel-pore arrangement and a uniform pore size [1].

Transformation of DGM equations to advection-dispersion equations
The molar concentrations of component gases A and B (C A and C B , respectively) (mol/L 3 ) are derived from the equation of state for an ideal gas as follows: The following equation is derived by summing Eqs. (1a) and (1b).
where P g is the total gas pressure (M/LT 2 ). Eq. (3) is rearranged as follows: Then both sides of Eq. (4) are multiplied by X A /t m u g D AB to obtain Eq. (5): Next, Eq. (6) is derived by subtracting Eq. (5) from Eq. (1a) and substituting Eq. (2a) into the result: Then, because C A = X A P g / RT, Eq. (6) can be transformed as follows: Thus, we define the compound diffusion coefficient of component A, D * A (L 2 /T), which includes both molecular and Knudsen diffusion, as in Eq. (7).
Therefore, by substituting Eq. (8) back into Eq. (7), and rearranging the result with respect to N D A , we can obtain: The total flux of component A, N T A , is calculated from N D A (Eq. 9) and the advection flux expression -C A k a rP g / m gmix , where k a is the apparent permeability coefficient (L 2 ) and m gmix is the viscosity of the gas mixture consisting of components A and B (M/LT), as follows: Here, we define V * gA (compound velocity) (L/T) as follows: We obtain Eq. (12) for N T A by combining Eqs. (10) and (11).
The mass conservation equation for component A in a porous medium is derived from N T A and u g as follows: @(u g C A ) /@t + r N T Next, Eq. (12) is substituted into Eq. (13) to obtain Eq. (14): The first and second terms on the left side of Eq. (14) can be differentiated as follows: Finally, by substituting the mass conservation equation for the total gas flux V * gA in a porous medium, @u g /@t + r V * gA = 0, into Eq. (15), the DGM advection-diffusion equation for component A can be derived from Eq. (15) as follows: where D' * A and V' * gA are the effective compound diffusion coefficient and the effective compound velocity of component A, respectively, and D' * A = D * A /u g and V' * gA = V * gA /u g . By referring to Eqs. (8) and (11), These coefficients can be expressed as follows: The advection-diffusion equation of component B is similarly derived: where D' * B = D * B / u g and V' * gB = V * gB /u g are the effective compound diffusion coefficient and the effective compound velocity of component B, respectively, which similar to those of component A can be expressed as follows: Calculation of the Knudsen diffusion coefficients To calculate the Knudsen diffusion coefficients, we first calculate the reciprocals of the effective compound diffusion coefficients (see Eqs. 17 and 20) and then substitute X B = 1 À X A and α = D B / D A into the result: Eqs. (22) and (23) explicitly indicate that X A and X B are proportional to 1 / D' * A and 1 / D' * B , respectively. If the slope of the proportional relationship between X A and 1/D' * A is m A and that between X B and 1/D' * B is m B, then m A and m B must be as follows: By substituting Eq. (27) into Eq. (24a) and rearranging, t m D AB can be obtained from the slopes of the respective proportional relationships as follows: The intercepts of the regression lines between X A and 1 / D' * A and between X B and 1 / D' * B are defined as respectively. Eq. (30b) is then subtracted from Eq. (30a) to obtain Eq. (31): Next, D B = αD A is substituted into Eq. (31), and the result is rearranged to obtain Eq. (32): Then Eq. (32) is rearranged: Finally, t p D B is derived from α = t p D B / t p D A as follows: Experimental procedure and method protocols

Tracer experiments in a column
In Section 2.1, we describe the tracer gas experiments conducted previously [5]. Then, in Section 2.2, we explain how an inversion simulation is used to fit the advection-diffusion equation to the experimentally obtained molar fraction distribution of a tracer gas for obtaining the effective compound diffusion coefficient and the effective compound velocity of that gas. In Section 2.3, the protocol for obtaining the Knudsen diffusion coefficient of tracer gases A and B from D'* A and D'* B is described Tracer experiments were conducted using columns of two different lengths made of acrylic resin [5]. For experiments with a rapidly diffusing tracer gas, a column with an outer diameter of 7 cm, an inner diameter of 5 cm, and a length of 150 cm was used (Fig. 1). This column was equipped with seven electrical pressure gauges (Tokyo Sokki Kenkyujo Co., Ltd., PW-100 kPa; range, 0-100 kPa; accuracy, 0.005 kPa) and 43 gas sample ports, each of which consisted of a connector (Koyo Co., Ltd., RGB 05,818) and a septum (Hamilton 9 mm DIA 12/PK) ( Fig. 1). For experiments with a tracer gas with a slowly diffusion tracer gas, a column with a length of 90 cm long and the same outer and inner diameters was used. This column was equipped with five electrical pressure gauges, installed at 20-cm intervals between 5 and 85 cm from one end of the column, and 17 gas sample ports, installed at 5-cm intervals from one end of the column. The pressure gauges of each column were connected to a data logger (Tokyo Sokki Kenkyujo Co., Ltd., TDS-303), and the gas pressure in the column was measured automatically. In each experiment, gases were injected via the pressure regulators (Fairchild Model 10; range, 0-15 kPa) into the column from gas tanks (Fig. 1, right side). The gases flowed through the column (from right to left in Fig. 1) and were released to the atmosphere (Fig. 1, left side) after first passing through a flow meter (Horiba Stc Co., Ltd., SEC-N100; maximum discharge, 300 cc/min; accuracy, 0.1 cc/min) and a backpressure regulator (Fairchild Model 10BP; maximum pressure, 15 kPa). The tracer experiments were conducted at a temperature of 20 C and a humidity of 50%. Each column was filled with a porous medium by placing it vertically, with cover plates and flow meter removed, on the laboratory floor. The filled column was then tightly closed with the cover plates and placed horizontally on a table.
In each experiment, a gas was first injected into the column until it filled the void volume in the porous medium. It was confirmed that this gas, called the background gas, had filled the column by using gas chromatography (GL-Sciences Inc., GC-3200 with a thermal conductivity detector; carrier gas, helium; oven temperature, 50 C; temperature in the injector, 50 C; temperature in the detector, 70 C) to analyze samples of gas from the column. Gas pressures at the inlet and outlet were adjusted to predetermined values by means of a gas pressure regulator installed at the inlet ( Fig. 1, right side) and a gas backpressure regulator installed at the outlet ( Fig. 1, left side). After it was confirmed that the column was filled with the background gas, a second gas, called the tracer gas, was injected into the column via the three-way valve. As the tracer gas was being injected into the column, gas pressure, barometric pressure, and the discharge of gas from the column were measured by the electrical pressure gauges, a barometer (Isuzu Seisakusho CO., LTD. Aneroid Barometer B-180-NO; range, 930-1070 hPa; accuracy, 1 hPa), and the flow meter, respectively. Then 0.15-mL gas samples were withdrawn from the column via the gas sample ports with syringes (SGE Analytical Science; maximum volume, 1.0 mL; minimum volume, 0.02 mL) at a specific time. The gas sampling did not affect the flow of the tracer gas in the column because only 0.15 mL of gas was sampled from each gas sample port and because the samples were collected immediately before the end of each experiment. These gas samples were analyzed by gas chromatography to determine the molar fraction of the tracer gas at each sample port at the time of sampling. Then the tracer gas and the background gas were exchanged, and another tracer experiment was conducted in the same way. For more accurate determination of the Knudsen diffusion coefficient, two to five tracer experiments were conducted with each tracer gas, one for each predetermined pressure difference between the inlet and outlet of the column. The distributions of the molar fractions of the two tracer gases were obtained from the results of these two sets of tracer experiments. For instance, the experimentally obtained distributions of the molar fractions of nitrogen and carbon dioxide (used as tracer gases) when the porous medium was oven-dried field soil (collected at Inazawa City, Japan; particle density, 2.61 g/cm 3 ; porosity, 0.533; particle diameter, 0.009 to 0.7 mm) are shown in Fig. 2.

Inversion simulation
The effective compound diffusion coefficients and the effective compound velocities of two tracer gases (gases A and B) are obtained by using an inversion simulation to fit the advection-diffusion equation to the experimentally determined distribution of the molar fraction of each tracer gas. Although for this fitting we used the ISCCFEM-adq inversion simulator developed by Hibi and colleagues [6,7], which employs a conjugate direction method [8] and a characteristic finite element scheme [9], any general inversion simulator for the advection-diffusion equation can be employed to obtain the effective compound diffusion coefficients and the effective compound velocities of the tracer gases from tracer experiment results.
The length of the one-dimensional domain for the inversion simulation, which is the length of the column used (i.e., 90 and 150 cm), is first For instance, in the case of the tracer experiments conducted with field soil, carbon dioxide gas, and nitrogen gas described in Section 2.1, the reciprocals of the effective compound diffusion coefficients of carbon dioxide and nitrogen, 1/D' * CO2 and 1/D' * N2 , respectively, were plotted against the molar fractions of carbon dioxide and nitrogen, X CO2 and X N2 , respectively (Fig. 3). When the pressure difference between the inlet and the outlet of the column was 3.0 kPa, the slope and the intercept of regression line of CO 2, m CO2 , and Y CO2 , were 11.286 s/m 2 and 21.57 s/m 2 , respectively, and those of N 2 , m N2 , and Y N2 , were -5.3805 s/m 2 and 12.419 s/m 2 , respectively (Fig. 3a). Therefore, the molecular diffusion coefficient between carbon dioxide and nitrogen, t m D CO2-N2 , was calculated to be 0.0973 by Eq. (28) from m CO2 and m N2 , and then α was calculated to be 2.1 by Eq. (29) from m CO2 and t m D CO2-N2 . The Knudsen diffusion coefficient of CO 2 , t p D CO2 , was calculated to be 0.0572 m 2 /s by Eq. (33) from α, Y CO2 , and Y N2 . Finally, the Knudsen diffusion coefficient of N 2 , t p D N2 , was calculated to be 0.120 m 2 /s by Eq. (33) from α and t p D CO2 .

Estimation of dispersion
The correlation coefficients of the regression lines, which were calculated from the coefficients of determination shown in Fig. 3, ranged from 0.68 to 0.72 for carbon dioxide and from -0.71 to -0.89 for nitrogen; thus, the reciprocal of the effective compound diffusion coefficient was correlated with the molar fraction of gas in each case. Therefore, the experimentally determined effective compound diffusion coefficient was directly proportional to the molar fraction of the tracer gas, as indicated by Eqs. (22) and (23). The effective compound diffusion coefficients obtained from the results of the tracer experiments, however, might include the effect of dispersion, which is directly proportional to the gas velocity. The effective dispersion coefficient of tracer gas A, D' dispÀA , includes the molecular diffusion coefficient t m D AB , the Knudsen diffusion coefficient t p D A , and the mechanical dispersion coefficient D mech . Thus, the effective dispersion coefficient D' dispÀA of gas A can be expressed as D' dispÀA = D' * A + D mech , if the assumption of Sleep [10] that such dispersion is similar to the dispersion of chemical components dissolved in groundwater [11] is accepted. As a result, Eq. (22) becomes as follows: The reciprocal of Eq. (35) is D' dispÀCO2 of carbon dioxide was calculated by using Eq. (36), the values of t m D CO2-N2 , α, and t p D CO2 obtained in Section 2.3, where gas A is CO 2 and gas B is N 2 , and each of four different values for D mech (0.0, 0.0005, 0.005, and 0.001 cm 2 /s). Then, 1/D' dispÀCO2 was plotted against the molar fraction of CO 2 , X CO2 , as shown in Fig. 4. However, although D' dispÀCO2 appears to be directly proportional to X CO2 in the figure, the mathematical relationship between D' dispÀCO2 and X CO2 described by Eq. (36) is not rigorously proportional. In fact, the coefficient of determination of the regression lines between D' dispÀCO2 and X CO2 , R 2 , shown in Fig. 4, decreased slightly, from 1.0000 to 0.9995, as the mechanical dispersion coefficient D mech increased; therefore, Eq. (36) is not obviously consistent with any particular line. Furthermore, the slope of the regression line also decreased as D mech increased, instead of remaining constant. In contrast, the relationships between 1/D'* CO2 and X CO2 and between 1/D'* N2 and X N2 obtained from the tracer experiment results by the protocol described in Section 2.3 could be expressed as linear equations. Further, when we plotted the slope m of the regression line for the relationship between 1/D' dispÀCO2 and X CO2 and the slope m of the regression line for the relationship between 1/ D' disp-N2 and X N2 against the pressure difference DP (Fig. 5), the R 2 values of the m-DP relationships for CO 2 and N 2 were very small, 0.0051 and 0.2016, respectively. Further, the slopes of the regression lines of the m-DP relationships were very small, -0.083 and -0.6115, respectively; thus, m did not change with DP but was approximately constant and independent of DP. Because in general D mech increases as DP increases, the slopes of the regression lines between 1/D' disp-CO2 and X CO2 and between 1/D' dispÀN2 and X N2 were constant, regardless of the value of D mech . Therefore, the relationship between D' dispÀCO2 and X CO2 obtained by Eq.  where D mech is a parameter that we call the "difference coefficient with mechanical dispersion", the value of which depends on the properties of the porous medium. We can similarly obtain Eqs. (24) to (34) (see Section 1.3) from Eqs. (37) and (38), except for Eqs. (30a) and (30b), which become as follows: Because when Eq. (39b) is subtracted from Eq. (39a), D mech is eliminated, Eq. (31) remains as shown in Section 1.3.

Accuracy of the experimentally determined Knudsen diffusion coefficients
We estimated the accuracy of the method used to obtain the Knudsen diffusion coefficients from the results of tracer experiment results. Ogata and Banks [12] mathematically solved a onedimensional advection-diffusion equation, with the boundary conditions that the concentration was C 0 at x = 0 and t = 0 and 0.0 at x = 1 and t = 0 and the initial condition that the concentration was 0.0 at x > 0 and t = 0, where x is the distance from the gas inlet and t is the elapsed time from the injection of gas as follows: where C is the concentration of the gas, D is the diffusion coefficient, and V is the velocity of the gas. The concentration calculated with Eq. (42) should have no numerical noise, unlike solutions obtained by numerical simulation. We use the ISCCFEM-adq inversion simulator (see Section 2.2) to obtain the diffusion coefficient from the distributions of concentrations calculated by Eq. (40) with C 0 = 1.0. The diffusion coefficient obtained by the inversion simulation should thus be equal to the diffusion coefficient D used in Eq. (40) for calculating the distribution of the concentrations. Therefore, the difference between the dispersion coefficient obtained by the inversion simulation and that used for calculating the distribution of concentrations is the error in the method we use to obtain the Knudsen diffusion coefficient from the tracer experiment results.
The length of the one-dimensional domain, the range used to divide the molar fraction into blocks, the cell size of the grid, and the number of time steps used in the inversion simulation were the same as those in the real tracer experiment. The length of the simulation domain was set equal to 90 cm, the length of the column used in the real tracer experiments. The inversion simulation was conducted for each block, and then the diffusion coefficient of each block was obtained from the concentrations calculated with Eq. (40). The observation points for the inversion simulations were arranged at x = 10, 20, and 25 cm and at intervals of 2.5 cm over the x = 30 cm, consistent with the locations of the gas sampling points in the real tracer experiments.
The distribution of concentrations was calculated for three cases: a velocity of 0.05 cm/s, a diffusion coefficient of 0.1 cm 2 /s, and an elapsed time from the injection of 1200s (Case 1); a velocity of 0.02 cm/s, a diffusion coefficient of 0.02 cm 2 /s, and an elapsed time from the injection of 2000s (Case 2); and a velocity of 0.005 cm/s, a diffusion coefficient 0.01 cm 2 /s, and an elapsed time from the injection of 8000 s and 11,000 s (Case 3). Concentrations obtained by the numerical inversion simulation were then fitted to the concentrations calculated with Eq. (40) at the indicated observation points. The results show that the diffusion coefficients obtained by the inversion simulation differed from those used to compute Eq. (40) by less than 3%, with the exception of Block 4 at the elapsed time 11,000 in Case 3 (Error 8%). The F statistic value, which is equal to the sum of the squares of the differences between the concentrations calculated by Eq. (42) and those obtained by the inversion simulation at the observation points, ranged from 4.62 Â 10 -9 to 8.38 Â 10 -5 (Table 1). Therefore, the ISCCFEM-adq inversion simulator used in this study was able to accurately simulate the diffusion coefficient.
The slope of the regression line between the molar fraction of the tracer gas and the reciprocal of the compound dispersion coefficient is an important parameter in this method for obtaining the Knudsen diffusion coefficient. In Eq. (40), the diffusion coefficient D is constant, independent of the molar fraction of the tracer gas. In Cases 1 and 2, the reciprocal of the diffusion coefficient was not correlated with the molar fraction of the gas but was constant everywhere, and the slopes of the regression lines were 0.0168 and -0.338 s/cm 2 , respectively. Though the slope of the regression line for all data at the elapsed time of 8000 s in Case 3 was 2.51 s/cm 2 , the slope of the regression line for the data of Blocks 1-3 only at the elapsed time of 8000 s in Case 3, 1.13 s/cm 2 , was less than half that for all data. At the elapsed time 11,000 in Case 3, the absolute slope of the regression line was 6.65 s/cm 2 for all data and 1.13 s/cm 2 for Blocks 1-3 only. Thus, similar to the situation at the elapsed time 8000 s in Case 3, at the elapsed time of 11,000 s in Case 3 the absolute slope of the regression line for all data was greater than that for data of Blocks 1-3 only. Furthermore, the absolute slope of the regression for data of Blocks 2-3 at elapsed times of both 8000 s and 11,000 s in Case 3 was 0.839, the smallest slope for Case 3 data (Fig. 6). Accordingly, we confirmed that this method can accurately reproduce the slope of the regression line between the reciprocal of the effective compound dispersion coefficient and the molar fraction of tracer gas when the reciprocal of the effective compound dispersion coefficient is less than 50 s/cm 2 . However, when the reciprocal of the effective compound dispersion coefficient is more than 100 s/cm 2 , the data in Block 1 (molar fraction 0.0-0.25) and Block 4 (molar fraction 0.75-1.0) must be excluded to accurately obtain the slope of the regression line. Further, many reciprocals of the compound dispersion coefficient must be obtained by fitting the advection-diffusion equation to the data at several different elapsed times with the inversion simulator to more accurately obtain the slope of the regression line.

Summary
To apply this method, the distributions of the molar fractions should first be obtained by conducting tracer experiments with a porous medium and a binary gas system consisting of components A and B. The tracer experiments are first conducted with component A as the tracer gas and component B as the background gas, and subsequently with component B as the tracer gas and component A as the background gas. Next, the advection-dispersion equation is fitted to the distribution of the molar fractions of tracer gas A or B obtained by the tracer experiments by means of an inversion simulation in order to obtain the effective compound dispersion coefficient and the effective compound velocity. The inversion simulation is conducted separately for four blocks (Block 1, molar fraction 0-0.25; Block 2, molar fraction 0.25-0.50; Block 3, molar fraction 0.50-0.75; and Block 4, molar fraction 0.75-1.0); thus, the effective compound dispersion coefficient and the effective compound velocity are obtained for each block. Further, the average molar fraction of each block is obtained by averaging the maximum and minimum molar fractions in the block. Then, the effective compound dispersion coefficients are plotted against the averages of the molar fractions, and two regression lines between the averages of the molar fraction and the reciprocal of the effective compound dispersion coefficient are calculated, one for tracer gas A and the other for tracer gas B, by the least-squares method. t m D AB is calculated from the slopes of the regression lines for tracer gas A and B, m A and m B , respectively, by using Eq. (28), and then α is calculated from the slope of the regression line for component A m A and t m D AB by using Eq. (29). t p D A can then be determined from α and the intercepts of the two regression lines between the average of the molar fractions and the reciprocal of the effective compound dispersion for tracer gas A and tracer gas B, Y A and Y B , respectively, by using Eq. (33). Finally, t p D B is calculated from α and t m D AB by using Eq. (34).

Additional information
The Knudsen diffusion coefficient can be obtained from the results of permeability experiments. Reinecke and Sleep [13] and Abu-El-Sha'r and Abriola [14] actually performed permeability experiments with the aim of obtaining the Knudsen diffusion coefficient.
In the case of a compression fluid in a porous medium, apparent gas permeability k a (L 2 ) in the porous medium can be obtained as follows: k a = q / [(P 1 2 -P 2 2 ) / 2Lm g P 2 ] (41) where P 1 and P 2 are the gas pressures at the inlet and the outlet, respectively, of a column filled with gas and a porous medium; q is the volumetric flux (L/T) discharged from the outlet; L is the length of the column packed with porous medium (L); and m g is viscosity (M/LT).
Klinkenberg [15] derived the following equation, which relates effective gas permeability k e (L 2 ) to apparent gas permeability k a by considering the slip effect, that gas velocity at the surface of a solid is not zero: k a = k e + bk e / P gave (42) where b is the Klinkenberg parameter (M/LT 2 ) and P gave is the average gas pressure. It is clear from Eq. (42) that a plot of k a against 1/ P gave yields a straight line with slope bk e and intercept k e . The Klinkenberg parameter can be obtained by dividing the slope by k e . If chemical component A composes a single-gas system in a porous medium, then the average gas pressure can be related to the total molar flux of component A, N T A (mol/L 2 T), by employing the DGM equation for a single-component gas system as follows:

N T
A RT / rP g = D A + k e P gave / m g If the gas in the porous medium is assumed to obey Darcy's law and the ideal gas law, the total molar flux of component A can be expressed as follows: N T A = -(P gave / RT)(k a /m g )rP g .
By rearranging Eq. (44), N T A RT / rP g can be calculated with Eq. (45) as follows: A plot of N T A RT / rP g against P gave in Eq. (43) yields a straight line with slope k e / m g , and the Knudsen diffusion coefficient of component A is the intercept. The effective gas permeability can be calculated by multiplying this slope by the viscosity.