Formation of bubbly horizon in liquid-saturated porous medium by surface temperature oscillation
We study non-isothermal diffusion transport of a weakly-soluble substance in a liquid-saturated porous medium being in contact with the reservoir of this substance. The surface temperature of the porous medium half-space oscillates in time, which results in a decaying solubility wave propagating deep into the porous medium. In such a system, the zones of saturated solution and non-dissolved phase coexist with the zones of undersaturated solution. The effect is firstly considered for the case of annual oscillation of the surface temperature of water-saturated ground being in contact with atmosphere. We reveal the phenomenon of formation of a near-surface bubbly horizon due to the temperature oscillation. An analytical theory of the phenomenon is developed. Further, the treatment is extended to the case of higher frequency oscillations and case of weakly-soluble solids and liquids.
pacs:47.55.db, 66.10.C-, 92.40.Kf
Diffusion transport in bubbly media Donaldson-etal-1997-1998 (); Haacke-Westbrook-Riley-2008 (); Goldobin-Brilliantov-2011 (); Krauzin-Goldobin-2014 () and media with condensed non-dissolved phase Haacke-Westbrook-Riley-2008 (); Davie-Buffett-2001 (); Goldobin-CRM-2013 (); Goldobin-etal-EPJE-2014 () appears to possess unique features for isothermal conditions and is even more intriguing for non-isothermal ones. For these systems, the non-dissolved phase makes the local solute concentration being equal to the solubility. The solute concentration is not a free variable any more but a function of temperature and pressure. Simultaneously, the non-zero divergence of the solute flux determined by the solubility gradient does not change the enslaved solute concentration but redistributes the mass of non-dissolved phase. Thus, for the dynamics of systems with non-dissolved phase, novel phenomena and mechanisms, which never appear for undersaturated solutions, come into play. These phenomena are especially strongly pronounced for systems where the non-dissolved phase is immobilised (e.g., in porous media) and solubility is small Goldobin-Brilliantov-2011 (). For an immobilised non-dissolved phase, the mass transport operates only through the solution, and when the solubility is small the mass accumulated in the non-dissolved phase can be several orders of magnitude larger than the dissolved mass.
In Ref. Krauzin-Goldobin-2014 (), we reported the effect of surface temperature oscillations, which produce the solubility wave, on the diffusive transport in porous media where the non-dissolved phase is present everywhere and discussed the physical systems for which the effect is relevant. The systems where the zones with non-dissolved phase can coexist with the zones of undersaturated solution can exhibit richer and more sophisticated dynamics. A liquid-saturated porous medium being in contact with the reservoir of weakly-soluble substance (e.g., atmosphere) is a system of that kind. In this paper we consider the effect of the surface temperature oscillation on diffusive transport in the porous medium half-space being in contact with the reservoir of a weakly-soluble substance.
For the sake of definiteness, we will consider transport of gases in detail and comment how the results can be extended to the case of weakly-soluble solids (e.g., limestone) and liquids (e.g., crude oil). The case of solids and liquids is different from the case of gases only by insensitivity of the solubility to pressure. From the view point of mathematics, for the processes under our consideration, the case of gas hydrates is also identical to the case of weakly-soluble solids.
Mathematically, the contact with atmosphere is represented via the boundary condition; at the contact boundary the solute concentration equals the solubility at any instant of time. The same boundary condition and, consequently, the phenomena we report in this paper will hold for the case of contact with reservoir of any substance.
In this paper, we will reveal that the temperature wave leads to formation of a near-surface bubbly horizon and ‘oversaturation’ of the medium with the atmosphere gas compared to the period-mean solubility. In particular, the net molar fraction (non-dissolved phase + solution) of the gas molecules in pores next to the surface equals the maximal-over-period solubility.
The phenomenon we will report is common and can be important for various systems with different origins of the surface temperature oscillations, including technological systems (filters, porous bodies of nuclear and chemical reactors, etc.). However, for the sake of convenience, we first focus on the case characterised by the hydrostatic pressure gradient which is significant for geological systems, where pressure doubles on the depth of meters, leading to a significant change of solubility. The no-pressure-gradient case will be considered as well.
The effect of enhanced filling of water-saturated ground with atmospheric gases creates more favorable conditions for local microflora and fauna and influences conditions for geochemical processes. With weakly-soluble solids, the effect provides opportunities for filling the porous matrix with some weakly-soluble “guest” substance in technology; the spatial mass distribution pattern of the “guest” substance can be controlled by the surface temperature waveform. For natural methane hydrate deposits in seafloor sediments, the effect of temperature waves on the deposit and the gas release from it is of interest in relation to the natural Glacial-Interglacial cycles iceage () and potential global climate change Hunter-etal-2013 ().
The paper is organised as follows. In Sec. II.1, we introduce the physical model for diffusion transport of weakly-soluble substance in the presence of non-dissolved phase and governing equations. With the results of numerical simulation of governing equations for a single-gas-component atmosphere, we demonstrate the phenomenon of formation of a bubbly horizon in Sec. II.2. The analytical theory of the phenomenon is developed in Sec. II.3. In Sec. II.4, we report the dependence of integral quantifiers of the bubbly horizon on the control parameters of the system. In Sec. III, the case of no effect of the hydrostatic pressure gradient is studied (the case corresponds to high-frequency temperature oscillations or solutions of condensed phases). The results are summarised in Conclusion (Sec. IV).
Ii Diffusion in saturated and undersaturated solutions
ii.1 Physical model and governing equations
The diffusion transport of solute in the presence of non-dissolved phase is essentially controlled by solubility. For moderate pressure and far from the solvent boiling point the solubility of gas in liquid reads Pierotti-1976 ()
where molar solubility is the molar amount of solute per 1 mole of solvent, and are reference values, the choice of which is guided merely by convenience reason, and is the solubility at the reference temperature and pressure; the parameter , with being the interaction energy between a solute molecule and the surrounding solvent molecules and being the Boltzmann constant, is provided in table 1 for several typical gases. For condensed matter (solids and liquids) the solubility is nearly independent of pressure and approximately reads
Geological systems are typically much more uniform in the horizontal directions than in the vertical one. Hence, it is reasonable to restrict our consideration to the one-dimensional case; the system is assumed to be homogeneous in the horizontal directions. We assume the -axis to be oriented downwards and its origin to be on the porous medium surface.
Let us consider harmonic oscillation of the surface temperature, , where is the mean temperature, is the oscillation amplitude, is the temperature oscillation cyclic frequency. In particular, annual oscillations of surface temperature only slightly deviate from their harmonic reduction (e.g., see Yershov-1998 ()). The heat diffusion equation with no-heat-flux condition deep below the surface (at infinity) and imposed surface temperature yields
where is the heat diffusivity and is the distance from the surface of porous medium. The hydrostatic pressure field reads
where is the atmospheric pressure, is the liquid density, and is the gravity.
Since the non-dissolved phase is immobilised in pores, the mass transport in the system is contributed solely by the diffusion through the intersticial liquid and governed by equation
where is the molar concentration of the solution, is the net molar fraction of gas molecules in the intersticial fluid, is the molar fraction of the gaseous phase (bubbles) in the intersticial fluid, is the effective molecular diffusion coefficient and is the thermodiffusion constant Bird-Stewart-Lightfoot-2007 (). As compared to the molecular diffusion coefficient in bulk of pure liquid, say , the effective coefficient is influenced by the pore network geometry (tortuosity) and the adsorption of the diffusing agent on porous matrix (on the time scales of our interest the adsorption does not lead to anomalous diffusion; it only changes the effective rate of normal diffusion Gregg-Sing-1982 ()). The importance of thermal diffusion was demonstrated for gases Goldobin-Brilliantov-2011 () and methane hydrate Goldobin-CRM-2013 (); Goldobin-etal-EPJE-2014 () on geological time scales, although for the system of our interest it can be neglected Krauzin-Goldobin-2014 (). Henceforth, we omit the thermal diffusion term. The solute concentration
i.e., it equals the solubility where the net amount of gas molecules exceeds the solubility, and equals , otherwise. In the latter case, .
The formulated mathematical model of the system implies that the dissolution process (as well as opposite process of formation of the non-dissolved phase from solution) occurs much faster than the change of the temperature field and the diffusive redistribution of solute mass on the macroscopic scales. In real systems, the dissolution time scales even for solid non-dissolved phase are assessed as hours (see Buffett-Zatsepina-2000 ()), which is small compared to the reference times of temperature oscillation and diffusive transport on the scale of the system. The hysteresis effects possible for some phase transformations in narrow pore channels Anderson-Tohidi-Webber-2000 () are neglected in our study.
Notice, Eq. (5) is accurate for the case where macroscopic porosity is spatially uniform and the non-dissolved phase occupies a negligible fraction of the pore volume, which holds true for gases and weakly-soluble solids and liquids.
For one-dimensional case, Eq. (5) takes the form
At the upper boundary we assume a contact with the atmosphere, which means
Deep below the surface we assume the no-flux condition and the absence of the non-dissolved phase;
Notice, two boundary conditions are required at ; however, due to specificity of our system, one boundary condition, Eq. (8), is sufficient at . Indeed, since is never less than , the value of at the point does not influence the system dynamics; the condition for it is redundant.
Generally, all material properties of the system depend on temperature and pressure. However, feasible relative variations of the absolute temperature are small. Hence, one can neglect variation of those parameters which depend on temperature polynomially and consider variation of only those parameters which depend on temperature exponentially: the latter parameters are solubility (1) and the molecular diffusion coefficient . Moreover, the only parameter sensitive to pressure is the gas solubility.
We employ the following dependence of molecular diffusion on temperature Bird-Stewart-Lightfoot-2007 ();
where is the dynamic viscosity of the solvent, is the effective radius of the solute molecules with the “coefficient of sliding friction” , . The dependence of dynamic viscosity on temperature can be described by a modified Frenkel formula Frenkel-1955 ()
For water, coefficient , ( is activation energy) and . For the effective diffusion coefficient we assume the same relative variation with temperature as for .
ii.2 Numerical results
Numerical simulation was performed for the atmosphere composed solely by nitrogen—the single-component atmosphere approximation. Technical details on the numerical simulation are provided in Appendix A. The simulation reveals, that for any initial condition, after a transient process, the system reaches a single stable time-periodic regime. Let us consider this regime.
The linear growth of the solubility with depth owned by the hydrostatic pressure gradient is modulated by the temperature wave (3). The oscillating solubility profile (1) for the temperature wave (3) and pressure (4) can be seen in Fig. 1. Oscillations of the solubility profile lead to formation of a nearly constant in time profile of the net amount of gas molecules . The profile of net molar fraction nearly attains the maximal (mid winter) solubility next to the surface, ; there the bubbly fraction exists for nearly entire year except a short coldest period. Further, monotonically decreases with depth, along with the decrease of the time interval when the bubbly phase is present, down to the depth where the bubbly phase never appears. Below the latter depth is nearly uniform and only slightly changes during the year. Non-uniformity of the profile in this zone rapidly decays with depth. The asymptotic value is close to the annual-mean surface solubility of the gas. These features of the regime can be well seen in Figs. 2 and 3.
The net molar fraction profile is nearly constant during the year because the molar diffusivity is 3 orders of magnitude smaller than the heat diffusivity , meaning that diffusive redistribution of mass is a slow process against the background of a rapid temperature (and, hence, solubility) oscillation. This well pronounced separation of time scales provides opportunity for developing an analytical theory of the process, allowing for a better insight into the mechanisms of the formation of the bubbly horizon. The results of numerical simulation can be as well more comprehensively understood in the context of this theory.
ii.3 Analytical theory
The principal assumption of our analytical theory is that the net molar fraction profile is ‘frozen’ for one oscillation period and the period-mean diffusion flux performs a slow diffusive transfer of solute mass. For the analytical treatment we also linearise temperature dependencies of diffusion coefficient
where , , .
Prior to constructing the analytical theory, let us emphasize, that this analytical theory is an approximation but not a limiting case. The lenearisations (12) and (13) require small . Meanwhile, for small the penetration depth of the bubbly zone is small (below in the text it will be shown to be nearly linear function of ) and can become commensurable with the thickness of the diffusion boundary layer . In the latter case the approximation of ‘frozen’ profile is invalid. Thus, the frozen profile approximation is not compatible with the limit of vanishing . Nonetheless, for moderate , both approximations are satisfactory accurate.
For analytical calculations we take in account the order of magnitude of the relevant system parameters; , , , for the annual temperature oscillation in wetlands.
At a steady regime the annual-mean diffusive flux is zero. For the porous medium domain where the bubbly phase appears for some part of the oscillation period, the mean flux reads
where is the oscillation period, are the time instants between which local temperature is high enough so that the local solubility becomes smaller than . Eq. (14) can be rewritten in terms of the temperature oscillation phase ;
with determined by the condition , i.e.,
The latter equation should be treated as an initial value problem; it should be integrated from till attains , the point where is the base of the bubbly zone. For , the bubbly phase never appears and is constant. More convenient is to deal with Eq. (17) in terms of , and ;
The latter equation should be integrated from , , which corresponds to the surface, till the equality is fulfilled, which corresponds to the base of the bubbly zone. In Figs. 2a and b the result of integration of Eq. (18) is plotted with the red dash-dotted lines. One can see that the analytical theory is in a good agreement with the results of numerical simulation.
Eq. (17) can be simplified for the limiting case (where the depth of penetration of the bubbly zone is also small, meaning );
The latter equation can be integrated and yields
For solution (20), coordinate varies from (surface) to (the base of the bubbly zone). The net molar fraction profile is given by . For this limiting solution, one can determine the scaling properties of the solution next to the surface, for ;
This limiting solution is plotted in Fig. 2 with the green dashed lines. Unfortunately, this limiting solution is accurate only for small amplitudes of temperature oscillations, where the diffusive boundary layer is commensurable with the bubbly zone thickness and thus the assumption of the ‘frozen’ profile is not valid for typical diffusivities in liquids.
In order to ensure the correctness of analytical calculations, numerical simulation was performed for a diminished molecular diffusion . The molecular diffusion coefficient was small enough for both assumptions of the analytical theory can be simultaneously accurate. In Figs. 2a and b, one can see that the results of analytical theory (18) match the results of numerical simulation for -fold diminished plotted with blue dotted lines. In Fig. 2c, the limiting analytical solution (20) matches the numerical results for the diffusion coefficient diminished by factor . Considering results of numerical simulation for diminished diffusivity, one can notice disappearance of the kink near the boundary between the bubbly horizon and the zone of undersaturated solution, which can be seen in Figs. 1 and 2a,b for ‘normal’ diffusion strength. Thus, one can conclude that this kink is related to finiteness of the ratio . The -profile is not completely frozen on the time scale of temperature oscillations; at the point of discontinuity of the solute concentration gradient even a small diffusivity can result in an observable distortion of the net profile .
The case of impaired diffusion can be also of physical interest for the porous media where the pore network is not globally connected and diffusive transport necessarily involves diffusion through the solid matrix material separating different connected clusters of pores from each other. In this case the effective diffusion coefficient will be diminished by several orders of magnitude compared to the bulk diffusivity in the pore liquid.
ii.4 Integral quantifiers of the bubbly horizon
The developed analytical theory suggests natural quantifiers of the -profile: the deep solute concentration and the integral strength of the bubbly layer . These characteristics are good quantifiers of the system state even for situations where the -profile is far from being frozen (see Fig. 2). In Fig. 3, quantifiers and are plotted for nitrogen atmosphere and water-saturated porous medium. The deviation of from the annual-mean near-surface solubility is a nonlinear effect, quadratic in amplitude . For non-vanishing the effect is negative, i.e., the temperature oscillation results in ‘ventilation’ of deep areas of the porous medium.
Iii The case of no effect of pressure gradient
For the annual temperature wave the penetration depth is about ; for this depth the hydrostatic pressure increase is comparable to the atmospheric pressure and, thus, the gas solubility is significantly influenced by the hydrostatic pressure gradient. For higher oscillation frequencies in nature (e.g., daily oscillations) and technological systems the hydrostatic pressure increase for the wave penetration depth can be negligible. The case of solutions of solids and liquids is qualitatively similar to the case of no hydrostatic pressure gradient, as the solubility of condensed matter is nearly independent of pressure (see Eq. (2)).
In Fig. 4, one can see that in the absence of pressure gradient the bubbly horizon is not bounded from below. For any depth the solubility profile exceeds the -profile (i.e., the bubbly phase locally disappears) for a small fraction of the oscillation period; is slightly smaller than the envelope of the wave of . For the limit of vanishingly small ratio the -profile coincides with the envelope of . While the envelopes of solubility profiles are different for gaseous and condensed matter, the principle, that -profile tends to coincide with the envelope of the solubility profile wave, remains valid for solids and liquids. Hence, the integral strength of the bubbly horizon approximately is
We have considered the diffusion transport of a weakly soluble substance in a liquid-saturated porous medium half-space being in contact with the reservoir of this substance in the case where the surface temperature harmonically oscillates in time. The surface temperature oscillation creates the temperature wave which propagates deep into porous medium and decays with depth. The solubility wave associated with the temperature wave results in time-dependent intermittency of the zones of non-dissolved phase with saturated solution and the zones of undersaturated solution.
Due to the smallness of the ratio , which is as small as for transport processes in liquids, the diffusion transport in the system is much slower than the temperature (and solubility) variation. As a result, the profile of the net molar fraction of ‘guest’ molecules in pores, (‘net’ means ‘solute + non-dissolved phase’), remains nearly constant over the oscillation period. We have revealed for gases that the -profile nearly attains the maximal-per-period solubility next to the surface and monotonously decays with depth in the zone where the non-dissolved phase can be observed, the bubbly horizon, and reaches a constant level beneath the bubbly horizon (Fig. 2). The bottom boundary of the bubbly horizon is controlled by the hydrostatic pressure gradient.
Without the pressure gradient, the bubbly horizon is not bounded from below and the -profile is slightly smaller than the envelope of the oscillating solubility profile (Fig. 4). At any depth the bubbly phase disappears only for a short part of the oscillation period; in the idealistic limiting case this part of oscillation period tends to . Notice, the solubility profile decays with depth exponentially and, thus, even though the bubbly horizon is not formally bounded from below, it becomes vanishingly weak at the depth. The no-pressure-gradient case is relevant for gases and short-term oscillations, when the penetration depth of the temperature wave is small compered to the depth scale of a noticeable increase of the hydrostatic pressure, and for condensed phases, solubility of which is insensitive to pressure.
For the annual temperature oscillation of amplitude the atmosphere gas bubbly horizon has been found to have penetration depth and near-surface relative increase of by compared to the no-oscillation solubility (Fig. 2a). Such an effect can be considered non-negligible for the local conditions for microflora and fauna. The dependence of integral quantifiers of the bubbly horizon on average surface temperature and temperature oscillation amplitude is plotted in Fig. 3.
The periods of negative temperatures with frozen groundwater are beyond the scope of this study and will be considered elsewhere.
Authors acknowledge financial support by the Government of Perm Region (Contract C-26/0004.3) and the Russian Foundation for Basic Research (project no. 14-01-31380_mol_a).
Appendix A Numerical simulation
The evolution of Eq. (7) was simulated with a finite difference method, explicit scheme with central differences. The time-dependent fields of temperature, molecular diffusion coefficient, and solubility were precalculated for one cycle of the surface temperature oscillation on the space–time grid used for simulation of Eq. (7). The calculation domain was limited from below by the depth , and boundary condition (9) was moved to ; the depth was chosen so that perturbations of the solute concentration profile were not visually resolvable at . The choice of the space stepsize was guided by two requirements. Firstly, the decrease of the stepsize by half should not change the results of simulation more than by . Secondly, for the diffusion coefficient diminished by factor the relative mismatch between the results of numerical simulation and unsimplified analytical solution (17) should not exceed . Practically, the grids with – points per the zone of non-dissolved phase were sufficiently detailed for these requirements to be met. The time stepsize was determined by the requirement of stability of the explicit method.
- (1) J. H. Donaldson, J. D. Istok, M. D. Humphrey, K. T. O’Reilly, et al., Development and Testing of a Kinetic Model for Oxygen Transport in Porous Media in the Presence of Trapped Gas, Ground Water 35, 270 (1997); J. H. Donaldson, J. D. Istok, and K. T. O’Reilly, Dissolved Gas Transport in the Presence of a Trapped Gas Phase: Experimental Evaluation of a Two-Dimensional Kinetic Model, ibid. 36, 133 (1998).
- (2) R. R. Haacke, G. K. Westbrook, and M. S. Riley, Controls on the formation and stability of gas hydrate-related bottom-simulating reflectors (BSRs): A case study from the west Svalbard continental slope, J. Geophys. Res., 113, B05104 (2008).
- (3) D. S. Goldobin and N. V. Brilliantov, Diffusive Counter Dispersion of Mass in Bubbly Media, Phys. Rev. E 84(5), 056328 (2011).
- (4) P. V. Krauzin and D. S. Goldobin, Effect of temperature wave on diffusive transport of weakly soluble substances in liquid-saturated porous media, Eur. Phys. J. Plus 129, 221 (2014).
- (5) M. K. Davie and B. A. Buffett, A numerical model for the formation of gas hydrate below the seafloor, J. Gephys. Res. B 106, 497 (2001).
- (6) D. S. Goldobin, Non-Fickian diffusion affects the relation between the salinity and hydrate capacity profiles in marine sediments, Comptes Rendus Mecanique 341, 386 (2013).
- (7) D. S. Goldobin, N. V. Brilliantov, J. Levesley, M. A. Lovell, et al., Non-Fickian Diffusion and the Accumulation of Methane Bubbles in Deep-Water Sediments, Eur. Phys. J. E 37, 45 (2014).
- (8) J. R. Petit, J. Jouzel, D. Raynaud, N. I. Barkov, et al., Climate and atmospheric history of the past 420000 years from the Vostok ice core, Antarctica, Nature 399, 429 (1999); EPICA community members, Eight glacial cycles from an Antarctic ice core, ibid. 429, 623 (2004).
- (9) S. J. Hunter, D. S. Goldobin, A. M. Haywood, A. Ridgwell, and J. G. Rees, Sensitivity of the global submarine hydrate inventory to scenarios of future climate change, Earth Planet. Sci. Lett. 367, 105 (2013).
- (10) R. A. Pierotti, A scaled particle theory of aqueous and nonaqueous solutions, Chem. Rev. 76(6), 717 (1976).
- (11) E. D. Yershov, General Geocryology (Cambridge University Press, New York, 1998).
- (12) R. B. Bird, W. E. Stewart, and E. N. Lightfoot, Transport Phenomena (Wiley, New York, 2007).
- (13) S. J. Gregg and K. S. W. Sing, Adsorption, Surface Area and Porosity, (Academic Press, New York, 1982).
- (14) B. A. Buffett and O. Y. Zatsepina, Formation of gas hydrate from dissolved gas in natural porous media, Mar. Geol. 164, 69–77 (2000).
- (15) R. Anderson, B. Tohidi, and J. B. W. Webber, Gas hydrate growth and dissociation in narrow pore networks: capillary inhibition and hysteresis phenomena, in Sediment-Hosted Gas Hydrates: New Insights on Natural and Synthetic Systems, edited by D. Long, M. A. Lovell, J. G. Rees, and C. A. Rochelle, Geological Society London Special Publications 319(1), 145–159 (2009).
- (16) J. I. Frekel, Kinetic theory of liquids (Dover Publications, New York, 1955).
- (17) V. I. Baranenko, V. S. Sysoev, L. N. Fal’kovskii, V. S. Kirov, et al., The solubility of nitrogen in water, Atomic Energy 68, 162 (1990); V. I. Baranenko, L. N. Fal’kovskii, V. S. Kirov, L. N. Kurnyk, et al., Solubility of oxygen and carbon dioxide in water, ibid., 342 (1990); S. Yamamoto, J. B. Alcauskas, and T. E. Crozier, Solubility of Methane in Distilled Water and Seawater, J. Chem. Eng. Data 21, 78 (1976).
- (18) P.T.H.M. Verhallen, L.J.P. Oomen, A.J.J.M.v.d. Elsen, A.J. Kruger, and J.M.H. Fortuin, The diffusion coefficients of helium, hydrogen, oxygen and nitrogen in water determined from the permeability of a stagnant liquid layer in the quasi-steady state, Chem. Eng. Sci. 39(11), 1535 (1984); W. Sachs, The diffusional transport of methane in liquid water: method and result of experimental investigation at elevated pressure, J. Petrol. Sci. Eng. 21, 153 (1998); R. E. Zeebe, On the molecular diffusion coefficients of dissloved CO2, HCO3-, and CO32- and their dependence on isotopic mass, Geochimica et Cosmochimica Acta 75, 2483 (2011).