A numerical solution for the diffusion equation in hydrogeologic systems

Abstract

This report documents a computer code for the numerical solution of the linear diffusion equation in one or two dimensions in Cartesian or cylindrical coordinates. Applications of the program include molecular diffusion, heat conduction, and fluid flow in confined ground-water systems. The flow media may be anisotropic and heterogenous. The model is formulated by replacing the continuous linear diffusion equation by discrete finite-difference approximations at each node in a blockcentered grid. The resulting matrix equation is solved by the method of preconditioned conjugate gradients. The conjugate gradient method does not require the estimation of iteration parameters and is guaranteed convergent in the absence of rounding error. The matrices are preconditioned to decrease the steps to convergence. The model allows the specification of various boundary conditions for any number of stress periods and the output of a summary table for selected nodes showing flux and the concentration of the flux quantity for each time step. The model is written in. a modular format for ease of modification. The model is verified by comparison of numerical and analytical solutions for cases of molecular diffusion, two-dimensional heat transfer, and axisymmetric radial saturated fluid flow. Application of the model to a hypothetical two-dimensional field situation of gas diffusion in the unsaturated zone is demonstrated. The input and output files are included as a check on program installation. The definition of variables, input requirements, flow chart, and program listing are included in the attachments. INTRODUCTION Interest in the unsaturated zone as a pathway for the movement of liquids, gases, and solutes has increased in recent years. A recognition of the importance of the process of diffusion in both aeration and the transport of contaminant gases through the unsaturated zone has arisen. The diffusion equation describes this process in mathematical terms. An efficient algorithm, the method of preconditioned conjugate gradients, which is particularly suited for the solution of the diffusion equation, has also received increased attention. The method is not widely available as a matrix-solving algorithm in numerical models. This relatively new method is attractive in that it does not require the estimation of iteration parameters, as do other iterative methods, and has a comparatively good rate of convergence. The report documents a numerical model that was developed to utilize this algorithm in a simple and efficient model that may be used in any process that is described by the diffusion equation. The model is a Fortran-77 1 com-puter program for the numerical solution of the linear diffusion equation in one or two dimensions in Cartesian or cylindrical coordinates. The model is formulated by replacing the continuous equation for a region with a set of finitedifference equations written for each node in a block-centered grid. The time derivatives are discretized by a backwards implicit approximation; the spatial derivative by a central difference. The resulting matrix equation is solved iteratively by the method of conjugate gradients. The processes of heat transport and saturated fluid flow also are described by the diffusion equation. Because the fundamental laws are of identical form, it is possible to model heat conduction and saturated fluid flow using the same solution algorithms. This is accomplished by changing the constant of proportionality appropriately. The constant of proportionality is the same in all directions only if the medium is isotropic. In the case of an anisotropic medium in two dimensions, there are two principal directions of diffusivity. The principal directions may be aligned with the xand z-axis to simplify the analytical and numerical solutions. In the model, the process modeled and the medium anisotropy are indicated by the values of the constant of proportionality and the ratio of anisotropy. The derivation of the diffusion equation is explained in terms of molecular diffusion, and the analogous equations for heat conduction and saturated fluid flow are described. The model is verified by comparison with analytical solutions for simple cases of molecular diffusion, heat transfer, and axisymmetric fluid flow. The application of the model to a hypothetical field situation is demonstrated. The definition of variables, input requirements, flow chart, and program listing are included in the attachments. THEORETICAL DEVELOPMENT The diffusion equation is a parabolic partial differential equation that describes transient flow under concentration, temperature, heat, and hydraulic head gradients. This equation may be derived from the fundamental laws that describe the steady-state flux of each quantity by using the law of conservation of mass or energy to write the equations of continuity. For simplicity, derivation of the equation is presented for the case of molecular diffusion in one dimension. The extension to the analogous transport equations will be described in the next section. of brand and trade names in this report is for identification purposes only and does not constitute endorsement by the U.S. Geological Survey. Molecular-Diffusion Equation Ordinary gaseous diffusion is the movement of gas resulting solely from the concentration gradient of the gas in an isobaric system. Pick's first law states that the mass flux of the diffusing gas is proportional to its change in concentration over space. This may be expressed for a binary system in one dimension as dC ^ = _ n _ FA DAB dx , where A is the mass flux of the diffusing gas A, D A TJ is the constant of proportionality commonly known as the binary 2 -1 molecular diffusivity for gas A into gas B, L T ; CA is the concentration of the diffusing gas A, ML , and x is distance in the direction of decreasing CA / L. The value of DAB depends on physical constants of the gas species as well as the temperature and pressure of the system. Values of DAB for various gases may be found in a standard reference such as the Perry's Chemical Engineers' Handbook (Perry and others, 1984). The diffusion equation is written by considering conservation of mass. Figure 1 shows an element in a one-dimensional closed system. Since there is no mass flux into or out of the system, the net mass flux into any element must equal the time rate of change of storage (concentration) within the element. Thus, the equation of continuity for gas A is _ _ DAB ax ax Since the total density, pressure, and temperature of the system are constant in this analysis, DAB is constant and the equation of continuity may be expressed as 3 2 CA 3CA , which is termed Pick's second law. For two and three dimensions, CA is differentiated twice over each additional direction. For cylindrical coordinates, Fick's second law is JAB 1 JL r 3r aca 3r a^c. at The derivation of the diffusion equation is discussed in depth in Bird and others (1960, p. 554-560).