A simple numerical approach for assessing coupled transport processes in partitioning systems

Applied Mathematical Modelling 25 (2001) 479±498

www.elsevier.nl/locate/apm

A simple numerical approach for assessing coupled transport processes in partitioning systems 1  M.G. Trefry *, J. Ohman , G.B. Davis

Centre for Groundwater Studies, CSIRO Land and Water, Private Bag No. 5, P.O. Wembley, WA 6913, Australia Received 9 August 1999; received in revised form 8 September 2000; accepted 5 October 2000

Abstract Recent concepts in the solution of multidomain equation systems are applied to the problem of distinct transport processes coupled over geometrically disjoint domains. The (time dependent) transport equations for the composite system are solved using a simple domain decomposition approach, with parallel implementations of detailed Schwarz balances for the system subdomain interfaces. An existing numerical partial di€erential equation (PDE) solver is coupled with the interface algorithms to provide a code capable of handling a wide range of dynamical equations within the subdomains. Interface partitioning conditions corresponding to sharply discontinuous Dirichlet constraints, and to (discontinuous) rate-limited Neumann constraints are also incorporated into the code. A variety of transport operators can be handled simply by altering the equation system code block. The code is validated against analytical solutions for representative parabolic transport equations including recent solutions for di€usive transport in partitioning laminates, useful for describing the movement of chemical species in composite materials. The code is then applied to an example problem of coupled multiphase chemical transport in a variably saturated soil column with a low-permeability capping. Ó 2001 Elsevier Science Inc. All rights reserved. Keywords: Partitioning; Domain decomposition; Interface; Transport; Validation

1. Introduction The application of continuum transport models [1] to the migration of chemical species in variably saturated porous media is widespread. This activity has been motivated by the need to assess, monitor and remediate chemical spills in natural systems. Recently, attention has moved to the risk assessment of more complicated systems, involving multiple interacting species moving within structured porous and permeable media with several coexisting physical phases, e.g., multiple gasoline products migrating from soil pro®les into vegetation cover, surface water systems, urban buildings and the atmosphere. For systems where media properties are discontinuous at the macroscale, continuum models must be enhanced to incorporate physico-chemical laws appropriate to the media interfaces, and to observe detailed balance conditions for mass and ¯ux. *

Corresponding author. Tel.: +61-8-9333-6286; fax: +61-8-9333-6211.  E-mail addresses: [email protected] (M.G. Trefry), [email protected] (J. Ohman), greg.davis@per. clw.csiro.au (G.B. Davis). 1 Present address: School of Engineering, Uppsala University, Sweden. 0307-904X/01/$ - see front matter Ó 2001 Elsevier Science Inc. All rights reserved. PII: S 0 3 0 7 - 9 0 4 X ( 0 0 ) 0 0 0 6 3 - 9

480

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

Furthermore, such macroscale discontinuities in media properties can also be co-located with changes in the fundamental transport physics, de®ning critical interfaces in the system. These interfaces act as dynamic mediators between the di€erent transport physics in the system submedia. These multiphysics concepts have previously been applied in ¯uid±solid systems [2,3], with recent algorithmic discussions by Doltsinis [4] and Bailey et al. [5]. In an environmental engineering context, one example of this type of multiphysics system is given by the heat ¯ow dynamics of a porous dam wall system, where the dam wall is heated by solar radiation and cooled by air convection and water evaporation on one side and water convection on the other. A simple single-medium solution may be deduced by representing the convection/radiation terms by effective boundary conditions along the dam wall. However a multiphysics approach may seek to generate coupled convective heat distributions in all three submedia (air, water and dam wall) in order to examine the environmental e€ects induced by the engineered structure, and even the mechanical e€ects induced within the structure itself. This may involve simultaneous solution of distinct sets of coupled ¯ow and heat transport equations and mechanical stress/deformation equations appropriate to the di€erent subdomains. These solutions would then be parameterised by the heat and mass transfer conditions applying at the submedia boundaries, e.g., the interfaces between the dam wall material and the atmosphere or the dam water, etc. Whilst interfacial conditions appropriate for heat ¯ow between static media are well understood, yielding continuous temperature distributions, not all dependent variables obey these types of continuity relations. More problematic interface conditions are presented by systems involving chemical partitioning interactions, as discussed by Trefry and Whyte [6]. Partitioning interactions can lead to sharp discontinuities in concentrations across physical phase interfaces. In these cases the discontinuities are parameterised by partitioning coecients, which are physico-chemical properties of the physical phases involved and are hence problem dependent. For porous media with microscopic structure these diculties may be ameliorated by using volume-averaged formulations (see Ref. [7] and references therein), however for systems with macroscopic structure such formulations are less useful and explicit solutions are required for the coupled transport and interface equations. Discretised numerical codes are routinely used to solve transport equations pertaining to irregular and structured domains, however many code implementations are designed for particular classes of partial di€erential equations (PDEs), making it hard to alter the system dynamics without rewriting the discretised di€erential operators. The ability to solve systems of coupled PDEs is less common, making codes suitable for multispecies transport rare. Furthermore, it is rarer still for discretised codes to incorporate discontinuity conditions at internal boundaries, limiting their ability to simulate partitioned systems. This work presents a ®rst attempt to develop a general solver code suitable for modelling multispecies partitioned transport in structured multiphysics systems. The approach used is to couple an existing non-linear PDE system solver [8,9] with a domain decomposition algorithm implemented with sharply discontinuous macroscale interface conditions. 2. General problem structure We con®ne attention to problems de®ned over 2D rectangular (macroscopic) domains. Neglecting the possibility of statistical distributions of media properties, heterogeneity in the solution domain is accounted for by decomposing the overall domain into non-overlapping, but contiguous, subdomains [10]. The choice of subdomain decomposition is arbitrary, but is usually matched with regions of similar material properties or with similar dynamics.

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

P k

481

P k

Fig. 1. Two-subdomain decomposition of the total domain X ˆ X1 [ X2 , showing the shared interface C12 between nonoverlapping subdomains X1 and X2 with respective material properties, k1 and k2 , and physical dynamics, P1 and P2 .

Fig. 1 shows a two-subdomain decomposition for a rectangular domain X with two zones, X1 and X2 , of di€ering material properties, k1 and k2 , and/or physical dynamics, P1 and P2 . Here X ˆ X1 [ X2 , the boundary of X is denoted by dX, and the shared interface between X1 and X2 is given by C12 ˆ X1 \ X2 . Extensions to cases involving more discontinuities in material properties k are obvious. For the present work, decompositions are further restricted to laminate structure as seen, for example, in strati®ed beds. Within each subdomain j, PDE operators Lj act on the Qj -dimensional solution vectors Q

uj ˆ …u1j ; . . . ; uqj ; . . . ; uj j †T ;

…1†

where T denotes the transposition operator. These PDE operators re¯ect the physical dynamics in the respective subdomains and need not be identical in structure, thereby permitting the description of a diverse range of physical problems in composite systems. Communication between the subdomains is mediated by interface operators Gi;j acting at the shared interfaces Ci;j . The problem statement is completed by specifying boundary conditions along dX. In view of the partitioned chemical transport application at issue, we restrict the PDE operators Lj to parabolic, time-dependent form, although the approach may be applied in a broader multiphysics context. This restriction, in turn, requires us to consider a time-stepping problem within the domain decomposition context. Again within the transport context, the interface operators Gi;j are restricted to linear forms appropriate to general ¯ux balance and continuity constraints. In particular, the chemical partitioning conditions treated by Trefry and Whyte [6] are catered for. These are the Instantaneous Equilibrium Condition (IEC) and the Mass-Limited Condition (MLC). The IEC describes a local equilibrium partitioning condition applying at an interface Cj;k between neighbouring solution vector components umj and unk , i.e., IEC interface operators umj …x; t† qmj …x; t†

aj;k unk …x; t† ˆ 0 qnk …x; t† ˆ 0

8x 2 Cj;k ;

…2†

where the q quantities are gradient ¯uxes (e.g., calculated by Fick laws) corresponding to the solution components, and aj;k is the equilibrium partitioning coecient between the two components. For gas±liquid interfaces, aj;k can be identi®ed with the Henry coecient [11]. The MLC yields an identical equilibrium partitioning relationship between the neighbouring components, but provides a contact resistance at intermediate times. This is achieved by de®ning an interface m;n which measures the departure from partitioning equilibrium of the interacting ``potential'' Pj;k solutions

482

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

MLC interface operators   m;n …x; t† ˆ bj;k umj …x; t† aj;k unk …x; t† Pj;k m;n qmj …x; t† ˆ qnk …x; t† ˆ Pj;k …x; t†

8x 2 Cj;k :

…3†

Here bj;k is a mass transfer or conductivity parameter for the interface Cj;k , mediating the interface equilibration rate. Note, in this case, that the ¯uxes either side of the interface remain equal, but are driven by the disequilibria in the dependent variables across the interface. Other interface operator formulations are possible within the present framework, but are not considered here. Together with the problem boundary conditions, the interface conditions represented by (2) or (3) are sucient to ®x the subdomain solution vectors throughout the problem domain, even though some subdomains may not intersect the system boundary dX. The overall solution is the union of all subdomain solutions, noting that subdomain solution vectors may be of di€erent lengths …Qj †, according to the dimensionalities of the individual Lj operators. 3. Method of solution Automated solution of the above problem can be achieved by coupling a generic PDE solver with a time-stepped interface algorithm within a structural domain decomposition. 3.1. PDE solvers The development of generic PDE solvers is an active ®eld. Apart from solver packages incorporated in major proprietary codes such as Mathematica [12,13] and Maple [14,15] and techniques for importing external objects into their environments, e.g., see Ref. [16], several standalone codes exist for solving arbitrary coupled systems of second-order PDEs. The DIFFPACK library [17] is a set of C++ objects for solving systems of PDEs over complicated geometric domains using ®nite element techniques. The commercial FIDISOL/CADSOL package [18] is a Fortran-77 ®nite di€erence solver for 2D and 3D elliptic and parabolic PDE systems. CADSOL allows the speci®cation of internal boundaries at which interface matching conditions can be applied. The public domain VLUGR package [8,9,19,20] solves arbitrary (second-order in space and ®rst order in time) PDE systems over quasi-rectangular domains in 2D and 3D. Attention is henceforth con®ned to the 2D ®nite di€erence solver VLUGR2 [9]. VLUGR2 solves the non-linear equation L…t; x; y; u; ut ; ux ; uy ; uxx ; uyy ; uxy † ˆ 0

…4†

with boundary conditions and initial conditions given by, respectively, B…t; x; y; u; ut ; ux ; uy † ˆ 0;

u…t0 ; x; y† ˆ u0 …x; y†:

…5†

In Eqs. (4) and (5), u is a vector of solution components, and the problem operators L and B are unrestricted in form, thereby permitting the solution of fully coupled systems of non-linear dynamical equations. Di€erent problems are speci®ed and solved by altering simple statements governing residual terms in two small subroutines. Clearly L and B are general enough to handle the parabolic transport equations outlined above, and related codes have already been applied to 2D density coupled systems [8,20]. A restriction of VLUGR2 is that it does not handle discontinuities in PDE coecients, as would typically arise in describing systems with composite heterogeneities. Special modi®cations to the code are required for such cases [20,21]. However since these discontinuities can be handled conveniently within the domain decomposition method,

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

483

VLUGR2 is suitable as a PDE solver engine for the present purposes. For example, VLUGR2 allows the primary dependent variables of each dynamical transport operator Lj to be associated with separate components of the corresponding solution vector uj , generating fully coupled solutions to the problem and avoiding operator splitting approximations. Accuracy of VLUGR2 solutions is monitored internally, with automated local grid re®nement algorithms activated where curvatures of dependent variables exceed pre-speci®ed tolerances. 3.2. Time-stepping interface algorithms In order to complete the solution method, a technique is required by which solution transients can be transmitted from one subdomain to another via the interface conditions. Iterative domain decomposition methods for steady problems [10,22] provide a starting point. The iterative Schwarz Alternating Method (SAM) and variants are in widespread use, with applications to viscous ¯ow [23], elliptic [24,25] and hyperbolic systems [26]. Domain decomposition methods have been described for local mesh re®nement algorithms [27], and ecient preconditioners have been constructed to solve the Schur complement subproblems [28,29]. The iterative SAM solves subdomain problems alternately, coupling the problems by exchanging function data at the subdomain boundaries [30]. Convergence of the iterative procedure has been proved for elliptic systems and for reduced hyperbolic and parabolic systems. More recently iteration-free time stepping methods have been discussed for parabolic systems [31]. In the present work a simpli®ed domain decomposition approach was used to facilitate the solution of the dynamical equations across discontinuities in the physico-chemical parameters. The partitioning interface operators described above are also ¯exible enough to handle normal discontinuities in media parameters, such as di€usion coecient, whilst still preserving ¯ux balance. For example, choosing a ˆ 1 in Eq. (2) yields a simple continuity condition appropriate for a discontinuity in medium properties, and a ˆ 1 in Eq. (3) yields a standard leakage condition at an interface between two media parameterised by the leakage coecient b. After decomposing the model domain into subdomains corresponding to distinct phases and/or media discontinuities, the resulting sets of partitioning interface operators were solved for each time step using an iterationfree approach. The approach is now sketched for a single Cartesian dimension with ®nite di€erence discretisation over the range ‰x0 ; xM‡N Š. Let the scalar solutions at time t Dt for a twodomain problem be represented by the arrays 

u1t

Dt

Subdomain 1 …x0 †; . . . ; ut1 Dt …xj †; . . . ; ut1

Dt

…xM † ;

 t u2

Dt

Subdomain 2 …xM †; . . . ; u2t Dt …xM‡k †; . . . ; ut2

Dt

…xM‡N † ;

…6†

where subdomains 1 and 2 contain M and N ®nite di€erence points separated by Dx1 and Dx2 , respectively; a partitioning interface is located at x ˆ xM . For the IEC case, the solutions can be advanced to t within each domain by solving L1 and L2 subject to zero-¯ux boundary conditions at x ˆ xM and physical system boundary conditions at x ˆ x0 and xM‡N . At the completion of the time step, the points adjacent to the interface (ut1 …xM Dx1 †, ut1 …xM † and ut2 …xM †, ut2 …xM ‡ Dx2 †) are updated to conform to Eq. (2) using a partitioned mass balance technique to estimate the interface Dirichlet conditions appropriate for the next time step, as follows. Let the volumetric measures for the dependent variables u1 and u2 be h1 and h2 , respectively. These measures are (problem dependent) input parameters, e.g., air-®lled porosity in an unsaturated soil transport application. At the end of time step t the integrated mass, W…t†, local to the interface, xM , is approximated by the four-point sum of dependent variables multiplied by volumetric measures

484

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

W…t† ˆ

 Dx1  h1 …xM Dx1 †ut1 …xM Dx1 † ‡ h1 …xM †ut1 …xM † 2  Dx2  ‡ h2 …xM †ut2 …xM † ‡ h2 …xM ‡ Dx2 †ut2 …xM ‡ Dx2 † : 2

…7†

A four-point scheme is used to overcome indeterminacies in boundary values arising from the zero-¯ux conditions applying at the interfacial subdomain boundaries. Before commencing the next time step, updated local values (denoted by ^ ) for u1 and u2 are determined to satisfying the mass conservation and partitioning conditions: ^ W…t† ˆ W…t†; t u^1 …xM † ˆ a1;2 u^t2 …xM †;

…8† …9†

where

i Dx1 h ^ W…t† ˆ h1 …xM Dx1 †^ ut1 …xM Dx1 † ‡ h1 …xM †^ ut1 …xM † 2 i Dx2 h ‡ ut2 …xM † ‡ h2 …xM ‡ Dx2 †^ ut2 …xM ‡ Dx2 † : h2 …xM †^ 2

…10†

The two inertia conditions u^t1 …xM Dx1 † ˆ u^t1 …xM † u^t2 …xM ‡ Dx1 † ˆ u^t2 …xM †

…11†

are used to supplement Eqs. (8) and (9). Solving Eqs. (8), (9) and (11) simultaneously yields the updated value u^t1 …xM † ˆ

‰a1;2 Dx1 …h1 …xM † ‡ h1 …xM

2a1;2 W…t† : Dx1 †† ‡ Dx2 …h2 …xM † ‡ h2 …xM ‡ Dx2 ††Š

…12†

Solutions for the remaining quantities u^t1 …xM Dx1 †, u^t2 …xM †, u^t2 …xM ‡ Dx2 † follow trivially. These updated values are inserted into the subdomain solution arrays for time t and the next time step is …xM Dx1 †, ut‡Dt …xM † and u2t‡Dt …xM †, then integrated in each subdomain, yielding ut‡Dt 1 1 t‡Dt u2 …xM ‡ Dx2 †. The above mass partitioning process is then applied to these values before the next time step is integrated, and so on. For an MLC interface, the subproblems are integrated for a time step using Neumann ¯uxes at the interfaces constructed from the potentials P (cf. Eq. (3)) determined from the previous time step solution using two-point estimators. That is out out ˆ h2 …xM †D2 …xM † 2 ˆ P1;2 …xM ; t†; …13† h1 …xM †D1 …xM † 1 ox xˆxM ox xˆxM where P1;2 …xM ; t† ˆ b1;2 ut1 …xM †

a1;2 ut2 …xM †




…14†

and D1 …x† is the di€usion coecient for subdomain 1, etc. Again the volumetric measures are problem input parameters. Eq. (13) provides the subdomain ¯ux boundary conditions appropriate for interface x ˆ xM for the time step commencing at time t. At the completion of each time step, the new values of ut1 …xM † and ut2 …xM † are used to set the next interfacial ¯uxes according to Eq. (13). In these interface solution methods, Dirichlet and Neumann information is gathered from the subdomains at the end of every time step and apportioned according to the appropriate

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

485

partitioning operators. For 2D cases where the subdomains do not have matching discretisations at the interface, a simple 6-point Aitken interpolation technique is used to resample the time step solutions so that the 1D interface equations outlined above can be applied on a point-by-point basis along the interface. These interface solution methods are simple and iteration-free, hinging on the solution of linear equations at each time step. However errors involved in this approximation can be signi®cant, especially for systems involving large temporal gradients at the interfaces. In particular, the time stepping errors are of O…Dt†, which are minimised by employing ®ne spatial and temporal discretisations, risking high computational expense. In succeeding sections the errors will be assessed in comparison with analytical solutions. 3.3. Structured Multiphase Flow in 2D (SMF2D)-code origin A simple domain decomposition code, SMF2D, was implemented using the above interface solution algorithms, with VLUGR2 as a subdomain PDE solver engine. SMF2D allows a wide range of non-linear PDE operators in multiple, disjoint, layered domains to be solved simultaneously subject to a choice of the two partitioning interface operators above. Coupled multiphase problems are explicitly available through the vector formulations inherent in VLUGR2. Even though the time-stepping and interface algorithms described above are clearly amenable to parallelisation, in the ®rst implementation stage a serial code was preferred, both for simplicity and benchmarking purposes. Another simpli®cation employed in SMF2D was to de-activate the automatic grid re®nement logic in VLUGR2. This is not a fundamental requirement of the interface algorithms; the grid re®nement logic can be accommodated within SMF2D with a modicum of coding e€ort. 4. Validation for representative transport systems The SMF2D code was recently validated against algebraic solutions for several key 1D transport equations, including advection±di€usion, ®rst-order loss, reversible sorption, and diffusive partitioning [32]. In the present work more scrutiny is placed on the behaviour of the interfacial algorithms under sharp front advection and partitioning with transient forcing. Validation 1 tests the IEC interface algorithm of SMF2D against a standard sharp front problem, while validations 4 and 5 test both SMF2D interface solution algorithms against new transient solutions for partitioning systems. All variables in the validations are expressed in arbitrary units. 4.1. Validation 1: Sharp-front advection Consider a saturated soil column of length L with an advecting front described by the equation oC o2 C ˆ Dxx 2 ot ox

v

oC ox

with boundary conditions oC ˆ 0 and C…0; t† ˆ 1 ox xˆL

…15†

…16†

and initial condition C…x; 0† ˆ 0. If L is suciently large the solution for the semi-in®nite domain may be used, thereby avoiding truncation errors associated with ®nite-domain solutions gained via eigenfunction expansions, e.g., Ref. [33]. The semi-in®nite domain solution is [34]

486

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

C…x; t† ˆ

       1 x vt vx x ‡ vt erfc p : erfc p ‡ exp 2 Dxx 2 Dxx t 2 Dxx t

…17†

Neglecting di€usion yields Dxx ˆ kv, where k is the longitudinal dispersion coecient and v is the steady ¯ow velocity. Smaller values of k lead to steeper advecting fronts. Numerical dispersion e€ects in the gridded solution of such problems are gauged in terms of the grid Peclet and Courant numbers, given for this example by Pˆ

Dx k

and C ˆ

vDt : Dx

…18†

Fig. 2 shows numerical results calculated for various values of P and C, using v ˆ 1 and L ˆ 50 and expressed as functions of dimensionless length x=vt. All solutions are calculated for t ˆ 15. Fig. 2(a) compares results from VLUGR2 with the analytical result of Eq. (17). In these cases VLUGR2 was used with local grid re®nement disabled (®xed spatial grids) but with variable time stepping (variable C). In the spirit of Ref. [34], where spatial errors were examined in isolation

Fig. 2. Code validations against algebraic solutions (solid lines) at t ˆ 15 of Eq. (17) with Dxx ˆ kv and v ˆ 1, unit concentration at x ˆ 0, zero ¯ux at x ˆ L, zero initial condition: (a) shows the performance of VLUGR2 with grid re®nement disabled; (b) shows the nature of the interfacial discretisation of SMF2D for an IEC condition with partitioning coecient a1;2 ˆ 1.

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

487

through the use of low Courant numbers, time stepping in the present case was similarly constrained to low values of C by prescribing ®ne VLUGR2 error tolerances. Values of dispersion coecient are as indicated. For the k ˆ 1 case, a Peclet number of 0.5 …Dx ˆ 0:5† gave reasonable agreement with the analytical result (root-mean-square (RMS) error 6:64  10 4 ), showing that the VLUGR2 solver displays little numerical dispersion for this spatial discretisation. Only at the most advanced regions of the front do the numerical predictions depart appreciably from the exact result. However for k ˆ 0:1 a grid Peclet number of 5 yielded signi®cant numerical dispersion (RMS error 1:06  10 2 ). Reducing P to 2 removes most of the numerical dispersion (RMS error 1:64  10 3 ). This result can be compared with the result of Unger et al. [34] who demonstrated a variable spatial and temporal weighting scheme for advection±dispersion problems with superior performance for higher values of P. The above results indicate that for advecting front problems, in the absence of grid re®nement, VLUGR2's centred in space backwards in time ®nite di€erence scheme is best restricted to reasonably low Peclet numbers, of the order of unity. This is not necessarily an onerous restriction, bearing modern computing facilities in mind. Furthermore, it is only to be expected that a general PDE solver cannot always compete in speci®c situations with purpose-built specialist codes, e.g. Refs. [34,35]. The neglect of grid re®nement improvements in the current example is an unrealistic restriction on the performance of VLUGR2. The major importance of the VLUGR2 approach is its ¯exibility in describing multiple physics processes. Fig. 2(b) examines the interfacial partitioning behaviour of the time-stepping scheme presented above and culminating in Eq. (12). Again, the advecting front example is used. This time an IEC interface is located at x ˆ 30. The IEC partitioning parameter a1;2 is chosen to be unity, along with the phase volumetric measures hi , thereby yielding a simple Fickian continuity condition at the interface. Analytically, the resulting solution is identical to the single domain case, but the numerical implementation of Eq. (12) will lead to interfacial discretisation errors. Fig. 2(b) plots the VLUGR2 numerical solution (triangles) for k ˆ 1 and P ˆ 0:5, as presented in Fig. 2(a). Also shown is a solution from SMF2D (dots) calculated using Dx ˆ 0:1 …P ˆ 0:1†, highlighting the behaviour of the four-point partitioning scheme at the interface x=vt ˆ 2. This solution is not optimal (RMS error 9:62  10 5 ) ± instead a sub-optimal discretisation is used for explanatory purposes to show the nature of the IEC interface discretisation process. Note that the two subdomains share the interfacial boundary, so that each subdomain contributes a concentration point on the interface. The unit value for a1;2 ensures that these points lie over each other for small values of C. The SMF2D solution was calculated with an upper limit of C 6 10 3 ; this mandates the VLUGR2 time integrations in each subdomain to cease every 10 3 time units so that the interfacial partitioning mass balances can be recalculated. The interfacial discretisation error is O…Dt†, leading to slight underestimations of concentration to the left of the interface, and overestimations to the right of the interface. This apparent negative numerical dispersion to the left of the interface is the result of a negative interfacial discretisation contribution added to a smaller positive numerical dispersion contribution. The overall interfacial errors can be minimised further by resorting to lower P and C values, which will also reduce numerical dispersion. Advection-free validations for a1;2 6ˆ 1 are shown in succeeding sections. 4.2. Validation 2: Partitioning interfaces with steady forcing Fig. 3 compares the output of SMF2D against an algebraic result at the indicated times for a twin-domain di€usive system (Eq. (15) with v ˆ 0 and h ˆ 1 in both domains) with a single partitioning interface at x ˆ x1 ˆ 1. The two interface partitioning conditions are used

488

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

Fig. 3. SMF2D code validation (squares) against algebraic solution (solid lines) at indicated times (numbers near

curves) for the twin-domain partitioning system of Eqs. (15) and (19) with (a) D1 ˆ 5, D2 ˆ 1=20, v1 ˆ v2 ˆ 0, zero ¯ux at x ˆ 0 and 2, unit initial concentration in domain 1, zero in domain 2, with IEC parameter a1;2 ˆ 2 at interface x ˆ 1; (b) Eqs. (15) and (19) with D1 ˆ 5, D2 ˆ 1=20, v1 ˆ v2 ˆ 0, zero ¯ux at x ˆ 0 and 2, unit initial concentration in domain 1, zero in domain 2, with MLC parameters a1;2 ˆ 1=2 and b1;2 ˆ 1 at interface x ˆ 1.

Instantaneous Equilibrium …IEC† u1 …x1 ; t† ˆ a1;2 u2 …x1 ; t†;

Mass-Limited …MLC† P1;2 …x1 ; t† ˆ b1;2 …u1 …x1 ; t† a1;2 u2 …x1 ; t†† ˆ D2 ou2 ˆ P1;2 …x1 ; t†: ox

D1 ouox1 x

1

x1

…19† Fig. 3(a) shows results for an interface subject to IEC partitioning, while Fig. 3(b) shows an MLC case. Boundary conditions and parameter values are supplied in the ®gure caption. The algebraic results of Ref. [6] are used in the ®gures. For the IEC case, discretisations of Dx ˆ 0:01 and Dt ˆ 0:001 were used in the simulations, yielding a total RMS numerical error over the two di€erent time solutions of 3:16  10 2 , while for the MLC case discretisations of Dx ˆ 0:005 and Dt ˆ 0:001 yielded a total RMS numerical error of 1:49  10 2 . This validation tests the accuracy of both the internal VLUGR2 solver engine and the subdomain interface time-stepping algorithms. The interface solution method requires that the relaxation of solutions within subdomains during each time step be low, so that the partitioning relationships set up at the beginning of each step are still useful estimates of the true partitioning relationships by the end of the subdomain solution time step. For twin di€usive subdomain problems characterized by di€usion coecients D1 and D2 , and which are stable for spatial discretisations Dx1 and Dx2 , the empirical relationship  2  Dx1 Dx22 ; …20† Dt < min D1 D2 provides a useful gross upper bound on the choice of time step Dt for the overall problem. Note that Eq. (20) is independent of the partitioning coecient.

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

489

4.3. Validation 3: Partitioning interfaces with transient boundary forcing This validation section tests the accuracy of both the internal VLUGR2 solver engine and the subdomain interface time-stepping algorithms for transient boundary conditions. Fig. 4(a) compares the output of SMF2D against algebraic results at the indicated times for a harmonically forced twin-domain di€usive system (Eq. (15) with v ˆ 0 and h ˆ 1 in both domains) with a single partitioning interface at x ˆ x1 ˆ 1. A stationary sinusoidal Dirichlet condition is applied to the boundary x ˆ 2. Algebraic solutions for the harmonically forced partitioning system are presented by Trefry [36]. A numerical solution from SMF2D was gained by de®ning a piecewise

Fig. 4. SMF2D code validation (dots) against independent solutions (solid lines, see text) at indicated times (numbers on curves in (a) and (b)) for a twin-domain partitioning system subject to transient boundary conditions. Subdomains are governed by Eq. (15) with D1 ˆ 1=10, D2 ˆ 1=2, v1 ˆ v2 ˆ 0, zero ¯ux at x ˆ 0: (a) has equilibrated boundary transient 1 cos…2pt† with IEC parameter a1;2 ˆ 3 and MLC parameters a1;2 ˆ b1;2 =10 ˆ 2 at interface x ˆ 0:5; results shown at 3/4 period. (b) and (c) show concentrations in the IEC system of part (a) induced by the Gaussian boundary pulse exp… …t 4†2 †, as functions of location for the four di€erent times indicated (b), and as functions of time for the three di€erent locations (numbers on curves) indicated (c) including the boundary condition (dashed line).

490

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

constant initial concentration condition corresponding to the equilibrium partitioning state, and exposing this condition to a repeating sinusoidal Dirichlet condition at the boundary. The simulation was run for twenty periods of the sinusoid in order to achieve quasi-stationary conditions. The total RMS error for the IEC simulation …a1;2 ˆ 3† was 1:54  10 2 , and for the MLC simulation …a1;2 ˆ 2; b1;2 ˆ 0:2† was 4:39  10 3 , each evaluated after more than 200 000 time steps. Discretisations of Dx ˆ 0:01 and t ˆ 10 4 were used in both simulations. Algebraic solutions for a pulse boundary condition can be derived using the above harmonic solutions [36] in conjunction with a Fourier spectral method. A Dirichlet pulse w…t† is applied at the boundary x ˆ 1. In brief, the boundary pulse function is mapped into the frequency domain by the Fourier transformation, F, Z 1 1 w…x†  F …w…t†† ˆ p w…t† exp… ixt† dt …21† 2p 1 yielding the pulse spectrum w…x†. The dynamics of the pulse inside the model partitioning system are governed by the appropriate system Fourier coecients (which represent the modal dispersion relations for the system). These Fourier coecients have composite structure (according to the partitioned subdomain structure) and have been determined exactly for the relevant partitioning systems [36]. In the following validation example, attention is con®ned to the IEC partitioning case to simplify the Fourier inversions. If the Fourier coecient is represented by nIEC …x; x†, then the full pulsed system solution in x t space is recovered by performing the inverse Fourier transformation Z 1
1 1 w…x†nIEC …x; x† exp…ixt† dx: …22† CIEC …x; t† ˆ F w…x†nIEC …x; x† ˆ p 2p 1 In practice, even though the integrand of Eq. (22) is available in closed form, as is the MLC counterpart, the integrals are often intractable and numerical quadratures are resorted to. For the present work, the numerical inverse Fourier transformation function of the computational package Mathematica [12] is used to perform the quadratures of Eq. (22), with standard parameter settings. Figs. 4(b) and (c) show results for the system in Fig. 4(a) but subject to a Gaussian boundary pulse at x ˆ 1. Boundary conditions and parameter values are supplied in the ®gure caption. Discretisations of Dx ˆ 0:01 and Dt ˆ 10 4 were used in the simulations, yielding a total RMS numerical error with respect to the quadrature solution of 2:24  10 2 , summed over the four snapshots in Fig. 4(b), and of 1:33  10 2 over the two time series in Fig. 4(c). Overall Fig. 4 shows good agreement for the harmonic system (a), but worse agreement for the pulsed system (b), (c) especially in subdomain 1 …x 6 1=2†. This discrepancy can be reduced by choosing ®ner spatial and temporal discretisations, at the cost of greater computational expense. Nevertheless it should be recognised that the numerical pulsed solutions from SMF2D are compared with other (independent) numerical solutions and, arguments on the relative accuracies of the two classes of solutions notwithstanding, the comparisons still show useful agreement between the two. The validation studies of this section show that the SMF2D code provides useful numerical solutions to a broad variety of problems involving chemical transport in mixed media characterized by moderate partitioning coecients. Numerical convergence properties for more stringent partitioning conditions are discussed elsewhere [37]. Many other problems could be devised to test the code, including the addition of non-linear operators, e.g., for Freundlich sorption isotherms, however algebraic solutions are often lacking for such systems. As is emphasized in Validation 1, the core PDE solver cannot always compete with purpose-built solvers in all scenarios, e.g., advecting fronts. This de®ciency can often be compensated for by moving to ®ner grid

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

491

discretisations. Nevertheless, for this reason, the approach outlined here is best viewed as a rapid assessment tool for complex interacting systems, able to indicate key system processes/relationships with little coding development. In the next section a realistic application of the SMF2D code is presented. 5. Application: Volatile transport in variably saturated partitioning composite media As an illustrative example of the utility of the SMF2D code, we consider the problem of unforced 1D transport of Benzene vapour through the vadose zone above a steady water table. In the following text the problem is decomposed into two subdomains with a partitioning interface. In subdomain one (the soil column) the transport physics is summarised in terms of 3 coupled PDEs, in the other subdomain (an arti®cial cap) a single PDE is sucient. Interfacial coupling of the subdomain solutions occurs between a single component phase in the soil column and the sole component in the arti®cial cap. This system corresponds to an experimental arrangement used in a study of in situ vapour degradation kinetics in the vadose zone (see [38] for details). The solution domain extends from the water table at z ˆ 0 to the ground surface at zmax ˆ 3:0015 m, including a permeable polymer capping (1.5 mm thickness) at ground level (Fig. 5). The capping is intended to provide an approximately impermeable upper boundary to the vadose zone; the behaviour of the capping is examined in detail here. Benzene vapour is produced by a smeared NAPL (non-aqueous phase liquid) distribution of gasoline product just above the water table. In the vadose zone soil pro®le, mass transfer of Benzene takes place between the soil gas (g), soil water (w) and soil organic content (oc) phases of the soil continuum. Benzene is mobile in the soil gas and water phases, and is immobile (sorbed) in the soil organic content phase. At the microscopic level, instantaneous equilibration across gas±water interfaces is assumed, whilst at the macroscopic interface between the soil and the polymer capping a partitioning interface condition is assumed for Benzene between the soil gas and polymer (p) phases.

Fig. 5. Schematic of the solution domain for the example application. The domain extends from z ˆ 0 to zmax , and is subdivided into two subdomains at the interface z ˆ z1 . The saturated zone and the air atmosphere provide system boundary conditions. Schematic depth dependent pro®les are also shown for hg ; hw and foc , and for initial values of CNAPL and Cg .

492

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

A dynamic sorption relation is de®ned between the water phase and organic fraction. The relevant transport equations are, for the soil region (subdomain 1)     oCg oCw oS o oCg o oCw …23† ‡ hw ‡ foc qb ‡ ‡ kNAPL CNAPL ; ˆ Dg Dw hg ot ot oz oz ot oz oz Cg ˆ k H Cw ; oS ˆ ot

l…S

oCNAPL ˆ ot

…24† kfoc qb Cw † ˆ

 l S

 kfoc qb Cg ; kH

kNAPL CNAPL

…25† …26†

and for the permeable polymer capping region (subdomain 2)   oCp o oCp Dp ˆ : oz ot oz

…27†

In Eqs. (23)±(27), the quantities h are volumetric measures of the phases, foc is the soil organic carbon fraction, qb the soil bulk density, kH the Henry (partitioning) coecient for Benzene across an air±water interface, k and l are partitioning and mass-transfer parameters, respectively, for Benzene sorption between the water phase and the organic carbon fraction, and ag;p and bg;p are the interfacial parameters for partitioned transport between the soil gas and the polymer capping. The present system is insensitive to changes in l, but more sensitive to changes in k. The quantities Dj …j ˆ g; w† are the e€ective di€usion coecients, which for the gas and water phases are assumed in the respective phases by the empirical to be related to the molecular di€usion coecient Dmol j phase porosity relation [39] h10=3 j Dj ˆ 2 Dmol j ; n

…28†

where n ˆ hg ‡ hw is the total soil porosity. Inserting the local equilibrium relation (24) into Eq. (23) yields   oCg oS o oCg ˆ Dtot ‡ foc qb …hg ‡ hw =kH † ‡ kNAPL CNAPL ot oz ot oz ˆ Dtot

o2 Cg oDtot oCg ‡ ‡ kNAPL CNAPL oz2 oz oz

…29†

with the total e€ective di€usion coecient given by Dtot ˆ Dg ‡ Dw =kH : Using Eq. (28), the derivative of Dtot can be expressed by ! ! 10=3 10=3 mol mol 7=3 oDtot o hg h D 10 oh D h oh g w 7=3 w ˆ ‡ w w ; ˆ 2 Dmol Dmol ‡ w2 g hg oz 3n oz n2 g n kH oz kH oz

…30†

…31†

where the total porosity n is assumed to be independent of depth. The quantity kNAPL is the volatilisation constant for Benzene production from the NAPL pro®le CNAPL , modelled in depthdependent form by

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

kNAPL ˆ kvap

493

…Csat Cg † …Dg ‡ Dw =kH † ; Csat …hg ‡ hw =kH †

…32†

where kvap is an adjustable parameter and Csat is the saturated Benzene vapour concentration. Boundary conditions are a Dirichlet condition corresponding to a saturated vapour phase at the water table, and zero concentration at the upper surface of the capping, corresponding to a well-mixed atmospheric reservoir condition. At the interface …z ˆ z1 ˆ 3 m† between the soil and the polymer capping we have either an IEC condition Cg …z1 † ˆ ag;p Cp …z1 †

…33†

or an MLC condition (mass transfer with contact resistance)
oCg oCp Dtot ˆ Dp ˆ Pg;p …z1 †: Pg;p …z1 † ˆ bg;p Cg …z1 † ag;p Cp …z1 † oz zˆz1 oz zˆz1

…34†

Note that Dtot takes account of the relative contributions to Benzene transport of the soil water and gas phases, and hence the h factor is absent in the left-hand side of Eq. (34) (see Eq. (13)). For the polymer phase in subdomain 2 the h factor is unity. These interfacial conditions mediate partitioning balances between the gas and water phases in the soil and the polymer phase in the capping. Whilst the value of ag;p can be ®xed by equilibrium measurements (and hence the IEC condition is fully speci®ed mathematically), experimental values for the MLC bg;p rate parameter are lacking, as is information on which partitioning condition best represents the true Benzene dynamics at the soil gas ± polymer capping interface. Here, simulations for each interface condition are presented in order to highlight the ¯exibility of SMF2D in describing partitioning processes in a complex ®eld application. In general, the quantities hj ; n; foc ; qb are depth dependent, i.e., functions of z, and in real systems usually display signi®cant variability and structure. For the present purposes some of these quantities will be represented by continuous and smoothly varying functions, the remainder by constant values. Empirical values for these functions and the remaining parameters in Eqs.  (23)±(34), taken from Ohman [38], are summarised in Table 1. Fig. 5 shows depth pro®les for hj ; foc and initial concentration conditions (scaled to have comparable magnitudes), highlighting the non-linearity of the parameters. Polymer transport parameter values ag;p and Dp are typical for polyethylene membranes [40]. In the notation of SMF2D, the solution vector for subdomain 1 is of length 3, with entries corresponding to (1) the soil gas concentration, (2) the sorbed concentration and (3) the remnant NAPL phase concentration of Benzene. The soil water concentration is everywhere identi®ed with the soil gas concentration via Eq. (24). In subdomain 2 only a Table 1

Depth-dependent parameter functions for the Benzene vapour simulationsa Parameter hw =n

foc Cg …t ˆ 0† CNAPL …t ˆ 0† a

Function 0:0943 ‡ 0:0025

0:9057 …1 ‡ 15:06z2:19 †0:5434 0:000479z

0:00041z2 ‡ 0:000267z3

Csat exp‰ 8:0zŠ 0:011103 exp‰ 4:0…z

0:35†2 Š ‡ 0:092523 exp‰ 250:0…z

0:13†2 Š

Constant parameter values are qb ˆ 1:384 g cm 3 , n ˆ 0:477, kvap ˆ 0:2 m 2 , kH ˆ 0:25, Dmol w ˆ 8:8  10 ˆ 0:805 m2 d 1 , Dp ˆ 7:0  10 6 m2 d 1 , l ˆ 10 d 1 , k ˆ 80 cm3 g 1 , ag;p ˆ 0:01, Csat ˆ 3  10 4 M m 3 . Dmol g

5

m2 d 1 ,

494

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

single solution component is required, the polymer phase concentration. Initial conditions correspond to a double-Gaussian NAPL distribution, a Benzene vapour phase declining exponentially with elevation from the water table, and zero sorbed and polymer capping phases. These initial conditions were observed after a nitrogen purge of the soil gas in a vadose zone soil pro®le [38]. A sample code block describing the above transport equations is supplied in Table 2. This code block essentially portrays the work necessary to implement in SMF2D the multidomain coupled equation system used in the example application. According to the value of the subdomain index variable currentdomain (with values 1 or 2 for the example application), di€erent sets of transport equations are de®ned using standard VLUGR2 syntax. The summations are over the grid points in the appropriate domain. Interfacial couplings between the domains are handled elsewhere in the code by generalised mass balancing logic routines, as outlined in an earlier section. In essence, widely di€erent physical models can be studied merely by altering the equation system code block, recompiling, and supplying the necessary physical parameters via input. Fig. 6 shows snapshots of gas phase Benzene concentration calculated from the above transport model. Part (a) shows the concentration assuming a free-air …Cg ˆ 0† boundary at the soil surface, while part (b) shows the Benzene concentration determined using a zero-¯ux (perfect cap) soil±capping interface condition. A common feature is the rapid development of a hump in the concentration pro®le near z ˆ 1. This indicates volatilisation of Benzene from a peak in the assumed NAPL concentration at that location (see Fig. 6 and Table 1). Near the water table the high water content of the capillary zone tends to choke the transport of Benzene in the soil gas, leading to the persisting kink in the concentration pro®le near z ˆ 0. The solutions are grossly di€erent higher in the soil column. Parts (c) and (d) of Fig. 6 show solutions calculated for a leaky cap simulation employing an IEC soil±capping interface condition with vaporisation parameter kvap varying from 0.2 (part (c)) to the default value of 2 (part (d)). Vaporisation from the NAPL phase represents the dominant supply of Benzene in the system. This Benzene production is governed by the rate-limited function in Eq. (32). The leaky cap condition establishes a quasisteady ¯ux from the soil pro®le through the cap to the atmosphere within approximately one to two weeks, although times are longer for highly resistive interface models. Table 2

SMF2D code block for the coupled transport equations of the example applicationa if + + + +

10

20 a

else

endif

(currentdomain.eq.1) then do 10 I ˆ 1, NPTS(currentdomain) RES(I,1) ˆ R1(Z(I))*UT(I,1) + R2(Z(I))*UT(I,2) ) Dtot(Z(I))*UZZ(I,1) ) dDtotdz(Z(I))*UZ(I,1) ) kNAPL(Z(I))*U(I,3) RES(I,2) ˆ UT(I,2) + mu*(U(I,2) ) lambda*R2(Z(I))*U(I,1)/kH) RES(I,3) ˆ UT(I,3) + kNAPL(Z(I))*U(I,3) continue do 20 I ˆ 1, RES(I,1) ˆ RES(I,2) ˆ RES(I,3) ˆ continue

NPTS(currentdomain) UT(I,1) - Dp*UZZ(I,1) U(I,2) U(I,3)

The VLUGR2 solver sets the RES(I,J) residuals to zero in each time step. The above code block can be compared (in order) with Eqs. (29), (25)±(27).

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

495

Fig. 6. Simulated Benzene pro®les: (a) solutions were calculated using a free-air soil boundary condition; (b) solutions

were calculated with zero ¯ux (perfect cap) condition at the soil±capping interface; (c) solutions were calculated with kvap ˆ 0:2 and an IEC soil±capping interface condition; (d) solutions were calculated with an IEC soil±capping interface condition. All parameters are as listed in Table 1 except where otherwise noted. Note the changes in scales for the capping domain in (c) and (d). Numbers on curves indicate elapsed time in days from the cessation of nitrogen ¯ushing.

Fig. 7 plots capping e‚ux rates calculated for MLC interface conditions with di€erent values of interface conductivity bg;p , ranging from 0.1 (resistive interface) to in®nity (equivalent to an IEC interface condition). The ¯ux rates can be compared with the maximum e‚ux in the absence of capping, i.e., a free-air boundary to the soil column. Since the IEC condition does not incorporate a ¯ux-limiting mechanism (as does the MLC condition), at equilibrium the capping e‚ux rate calculated with the IEC interface condition would equal the free-air boundary rate, congruent to the integrated production rate of Benzene in the soil column, assuming the NAPL phase concentration remains constant. However the incorporation in the model of NAPL phase

Fig. 7. Simulated Benzene e‚uxes through the capping to the atmosphere, for di€erent interface conditions and ®xed kvap . For the MLC conditions the value of bg;p is indicated. Note the gradual declines in ¯uxes at large times.

496

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

loss via vaporisation leads to a free-air curve that rises more steeply than the IEC curve, but has lower peak ¯ux. For MLC interfaces, the ¯ux rates across the soil±capping interface are limited. Fig. 7 shows that the quasi-steady state occurs when the capping e‚ux equals the ¯ux across the soil±capping interface, where the magnitude of the latter depends upon the interface conductivity bg;p . Also discernible from Fig. 7 (especially for the IEC and free-air curves) are gradual declines in ¯ux rates for times beyond tens of days. These declines correlate well to a reduction in NAPL phase concentration, showing that cumulative volatilisation losses can impact on net e‚ux rates, although at considerably longer time scales. Overall, the results show that capping ¯uxes stabilise over time scales of the order of weeks after the nitrogen purge, yielding e‚uxes of the order of tens of micromoles of Benzene vapour per square metre of capping each day. Resistive soil±capping interface conditions (i.e., low values for bg;p ) lead to reduced vapour e‚uxes and longer stabilisation times. The capping provides a temporary barrier to Benzene e‚ux, delaying peak ¯uxes for times ranging from several days to several weeks according to the mass transfer conditions employed at the soil±capping interface. Furthermore, ¯ux magnitudes determined with the IEC interface condition were elevated slightly over comparable free-air boundary ¯uxes, but were progressively reduced for more resistive MLC interfacial conditions. Nevertheless, the simulations show that, for times exceeding several weeks after emplacement, the thin polymer capping is largely ine€ective as a barrier to Benzene vapour regardless of the partitioning mechanisms applying at the soil±capping interface, unless extremely resistive MLC conditions are operating. 6. Discussion The use of specialist PDE solver codes provides signi®cant simpli®cation in the development of new computational tools for studying coupled disjoint problems, e.g., chemical transport processes in composite media. Such problems are of particular interest in that the detailed transport processes across the media interfaces are not currently well understood. SMF2D uses one such solver code in a domain decomposition framework to model continuum transport in strati®ed systems with macroscopic phase interfaces. Initial validations of SMF2D with 1D algebraic solutions for key transport systems indicate robust performance and useful accuracy. The ability of the solver code to handle fully coupled multicomponent PDEs raises the prospect of studying in detail mass transfer characteristics in heterogeneous porous media contaminated by many interacting species in multiple physical phases. A salient feature of SMF2D is its implementation of interfacial partitioning laws, which provide ¯exibility in the description of transport between media. The ease of the approach is demonstrated by a soil vapour transport application, where the important code block for the governing equations reduces to a few lines for each subdomain. In this simple implementation of a time-dependent domain decomposition algorithm, subdomains communicate at the end of each time step using iteration-free balance equations applied to the previous time step solutions. Whilst this ``loosely coupled'' technique is simple, robust and easily parallelisable, it often requires ®ne temporal and spatial resolutions to provide acceptable solutions. This may partly be overcome by activating the full grid re®nement functionality of the VLUGR2 subdomain solver within the SMF2D framework. A more fundamental approach would be to incorporate more advanced domain decomposition algorithms appropriate for time dependent problems in order to improve the accuracy of the temporal discretisation scheme. Another drawback of the current method is the relative ineciency of using generic PDE solvers to study speci®c problems. In this sense the ease of application of the generic codes to diverse problems may in some circumstances be counterbalanced by the need for signi®cantly ®ner grids

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

497

and smaller time steps than might otherwise be required for specialist codes. Despite these issues, the general amenity of the approach to modelling multiphase multidomain problems is clear. Acknowledgements This work was part supported by BP Australia, Ltd, and by the Centre for Groundwater Studies. The authors wish to thank Dr P. Binning (University of Newcastle, Australia) for helpful comments. References [1] M.Y. Corapcioglu, A.L. Baehr, A compositional multiphase model for groundwater contamination by petroleum products: 1. Theoretical considerations, Water Resour. Res. 23 (1987) 191±200. [2] C. Bailey, P. Chow, M. Cross, Y. Fryer, K. Pericleous, Multiphysics modelling of the metals casting process, Proc. Roy. Soc. Lond. A Mat. 452 (1996) 459±486. [3] P. Adamidis, A. Beck, U. Becker-Lemgau, Y. Ding, M. Franzke, H. Holtho€, M. Laux, A. Muller, M. Munch, A. Reuter, B. Steckel, R. Tilch, Steel strip production ± a pilot application for coupled simulation with several calculation systems, J. Mater. Process. Tech. 80-1 (1998) 330±336. [4] I. St. Doltsinis, Solution of coupled systems by distinct operators, Eng. Comput. 14 (1997) 829±868. [5] C. Bailey, G.A. Taylor, M. Cross, P. Chow, Discretisation procedures for multi-physics phenomena, J. Comput. Appl. Math. 103 (1999) 3±17. [6] M.G. Trefry, D.S. Whyte, Analytical solutions for partitioned di€usion in laminates: I. Initial value problem with steady Cauchy conditions, Transport in Porous Media 37 (1999) 93±128. [7] W.G. Gray, Thermodynamics and constitutive theory for multiphase porous-media ¯ow considering internal geometric constraints, Adv. Water Resour. 22 (1999) 521±547. [8] R.A. Trompert, J.G. Verwer, J.G. Blom, Computing brine transport in porous media with an adaptive-grid method, Int. J. Numer. Meth. ¯uids 16 (1993) 43±63. [9] J.G. Blom, R.A. Trompert, J.G. Verwer, Algorithm 758: VLUGR2: A vectorizable adaptive grid solver for PDEs in 2D, ACM T. Math. Software 22 (1996) 302±328 (see also URL http://www.netlib.org/toms/758). [10] T.F. Chan, R. Glowinski, J. Periaux, O. Widlund (Eds.), in: Proceedings of the First International Symposium on Domain Decomposition Methods for Partial Di€erential Equations, SIAM, Philadelphia, 1988. [11] C.R. Chiou, D.W. Schmedding, Partitioning of organic compounds in octanol±water systems, Environ. Sci. Technol. 16 (1982) 4±10. [12] S. Wolfram, in: The Mathematica Book, third ed., Wolfram Media/Cambridge University Press, Cambridge, 1996. [13] V.G. Ganzha, E.V. Vorozhtsov, in: Numerical Solutions for Partial Di€erential Equations: Problem Solving Using Mathematica, CRC Press, Boca Raton, 1996. [14] M. Monagan, K. Geddes, K. Heal, G. Labahn, S. Vorkoetter, in: Maple V Programming Guide for Release 5, Springer, Berlin, 1997. [15] G.A. Articolo, in: Partial Di€erential Equations and Boundary Value Problems with Maple V, Academic Press, New York, 1998. [16] M.G. Trefry, P.C. Abbott (Eds.), The InterCall Book, Analytica International, Perth, Western Australia, 1998 URL: http://www.physics.uwa.edu.au/ICB/. [17] A.M. Bruaset, H.P. Langtangen, A comprehensive set of tools for solving partial di€erential equations; Di€pack, in: M. Dhlen, A. Tveito (Eds.), Numerical Methods and Software Tools in Industrial Mathematics, Birkh auser, Boston, Cambridge MA, 1997 (see also Di€pack Online at URL: http://www.oslo.sintef.no/NAM/people/xic/ di€pack_www/). [18] M. Schmauder, R. Weiss, W. Sch onauer, The CADSOL program package, Interner Bericht Nr. 46/92, Universit at Karlsruhe, Rechenzentrum, 1992 URL: http://www.uni-karlsruhe.de/Uni/RZ/Forschung/Numerik/®disol.cadsol/ index.html. [19] J.G. Blom, J.G. Verwer, VLUGR3 ± A vectorizable adaptive grid solver for PDEs in 3D. I. Algorithmic aspects and applications, Appl. Numer. Math. 16 (1994) 129±156. [20] R. Trompert, Local-uniform-grid re®nement and transport in heterogeneous porous media, Adv. Water Resour. 16 (1993) 293±304.

498

M.G. Trefry et al. / Appl. Math. Modelling 25 (2001) 479±498

[21] M. Arad, G. Ben-Dor, A. Yakhot, High-order-accurate discretization of a second-order equation with discontinuous coecients, Appl. Math. Modelling 22 (1998) 69±79. [22] D. Funaro, A. Quarteroni, P. Zanolli, An iterative procedure with interface relaxation for domain decomposition methods, SIAM J. Numer. Anal. 25 (1988) 1213±1236. [23] J.C. Strickwerda, C.D. Scarbnick, A domain decomposition method for incompressible viscous ¯ow, SIAM J. Sci. Comput. 14 (1993) 49±67. [24] D. Yang, A parallel iterative nonoverlapping domain decomposition procedure for elliptic problems, IMA J. Numer. Anal. 16 (1996) 75±91. [25] J.R. Rice, E.A. Vavalis, D. Yang, Analysis of a nonoverlapping domain decomposition method for elliptic partial di€erential equations, J. Comput. Appl. Math. 87 (1997) 11±19. [26] A. Quarteroni, Domain decomposition methods for systems of conservation laws: spectral collocation approximations, SIAM J. Sci. Stat. Comput. 11 (1990) 1029±1052. [27] W.D. Gropp, D.E. Keyes, Domain decomposition with local mesh re®nement, SIAM J. Sci. Stat. Comput. 13 (1992) 967±993. [28] P.E. Bjùrstad, O.B. Widlund, Iterative methods for the solution of elliptic problems on regions partitioned into substructures, SIAM J. Numer. Anal. 23 (1986) 1097±1120. [29] C.H. Tong, T.F. Chan, C.C.J. Kuo, A domain decomposition preconditioner based on a change to a multilevel nodal basis, SIAM J. Sci. Stat. Comput. 12 (1991) 1486±1495. [30] B. Engquist, H.-K. Zhao, Absorbing boundary conditions for domain decomposition, Appl. Numer. Math. 27 (1998) 341±365. [31] T.P. Mathew, P.L. Polyakov, G. Russo, J. Wang, Domain decomposition operator splittings for the solution of parabolic equations, SIAM J. Sci. Comput. 19 (1998) 912±932. [32] M.G. Trefry, SMF2D ± A general continuum code for simulating transport in strati®ed porous media: development and validation, in contaminated site remediation: challenges posed by urban and industrial contaminants, in: C.D. Johnston (Ed.), Proceedings of the 1999 Contaminated Site Remediation Conference, Fremantle, Western Australia, 21±25 March 1999, pp. 275±281. [33] R.W. Cleary, D.D. Adrian, Analytical solution of the convective±dispersive equation for cation adsorption in soils, Soils Sci. Soc. Am. Proc. 37 (1973) 197±199. [34] A.J.A. Unger, P.A. Forsyth, E.A. Sudicky, Variable spatial and temporal weighting schemes for use in multi-phase compositional problems, Adv. Water Resour. 19 (1996) 1±27. [35] A.E. Adenekan, T.W. Patzek, K. Pruess, Modeling of multiphase transport of multicomponent organic contaminants and heat in the subsurface: numerical model formulation, Water Resour. Res. 29 (1993) 3727±3740. [36] M.G. Trefry, Analytical solutions for partitioned di€usion in laminates: II. Harmonic forcing conditions, Transport in Porous Media 37 (1999) 183±212.  [37] M.G. Trefry, J. Ohman, D.S. Whyte, G.B. Davis, A time-synchronous domain decomposition code for multiphysics systems, ANZIAM J. 42E (2000) (to appear) URL: http://anziamj.austms.org.au/V42/CTAC99.  [38] J. Ohman, Modelling transport and degradation dynamics of BTEX vapours in a contaminated vadose zone, M.Sc. Thesis, Uppsala University School of Engineering, Sweden, 1999, 64 pp. [39] R.J. Millington, J.P. Quirk, Permeability of porous solids, Trans. Faraday Soc. 57 (1961) 1200±1207. [40] J.K. Park, M. Nibras, Mass ¯ux of organic chemicals through polyethylene geomembranes, Water Environ. Res. 65 (1993) 227±237.