ANALYTICAL SOLUTIONS FOR SEQUENTIALLY COUPLED MULTI-SPECIES REACTIVE TRANSPORT PROBLEMS Except where reference is made to the work of others, the work described in this thesis is my own or was done in collaboration with my advisory committee. This thesis does not include proprietary or classified information. Venkatraman Srinivasan Certificate of Approval: Mark O. Barnett Associate Professor Civil Engineering T. Prabhakar Clement, Chair Professor Civil Engineering Dongye Zhao Associate Professor Civil Engineering George T. Flowers Interim Dean Graduate School ANALYTICAL SOLUTIONS FOR SEQUENTIALLY COUPLED MULTI-SPECIES REACTIVE TRANSPORT PROBLEMS Venkatraman Srinivasan A Thesis Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Master of Science Auburn, Alabama December 17, 2007 iii ANALYTICAL SOLUTIONS FOR SEQUENTIALLY COUPLED MULTI-SPECIES REACTIVE TRANSPORT PROBLEMS Venkatraman Srinivasan Permission is granted to Auburn University to make copies of this thesis at its discretion, upon request of individuals or institutions and at their expense. The author reserves all publication rights. Signature of Author Date of Graduation iv THESIS ABSTRACT ANALYTICAL SOLUTIONS FOR SEQUENTIALLY COUPLED MULTI-SPECIES REACTIVE TRANSPORT PROBLEMS Venkatraman Srinivasan Master of Science, December 17, 2007 (B.Tech., Indian Institute of Technology Madras, Chennai, India, 2004) 216 Typed Pages Directed by T. Prabhakar Clement Multi-species reactive transport equations coupled through sorption and sequential first-order reactions are commonly used to model sites contaminated with radioactive wastes, chlorinated solvents and nitrogenous species. Although researchers have been attempting to solve various forms of this reactive transport problem for over fifty years, a general closed-form analytical solution to this problem is not available in the published literature. In the first part of this research work, a closed-form analytical solution to this problem is deduced involving a generic spatially-varying initial condition. Two distinct solutions are derived for Dirichlet and Cauchy boundary conditions each with Bateman-type source terms. The proposed solution procedure employs a combination of Laplace and linear transform methods to uncouple and solve the system of partial differential equations. The final solution is organized and presented in a v general format that represents the solutions to both the boundary conditions. In addition, the mathematical concepts for deriving the solution within a generic framework that can be used for solving similar transport problems are also presented. In the second part of this research work, the computational techniques for implementing the new solutions are discussed. These techniques are then adopted to develop a general computer code which is used to verify the solutions. In addition, several special-case solutions for simpler transport problems involving zero initial condition, identical retardation factors, zero advection, zero dispersion and steady-state condition are also derived. Where ever possible, these special-case solutions are compared against previously published analytical solutions to establish the validity of the new solution. The performance of the new solution is tested against other published analytical and semi-analytical solutions using a set of example problems. Finally, an investigation into extending the general solution to multiple dimensions using the approximate Domenico solution is also presented. This thesis has produced the following three journal publications: 1) V. Srinivasan, T.P. Clement, and K.K. Lee. ?Domenico Solution ? Is it Valid??, Ground Water, 25(2): 136-146, May 2007. 2) V. Srinivasan, and T.P. Clement. ?Analytical Solutions for Sequentially Coupled Reactive Transport Problems ? Part I: Mathematical Derivations?, Submitted May 2007, Advances in Water Resources. 3) V. Srinivasan, and T.P. Clement. ?Analytical Solutions for Sequentially Coupled Reactive Transport Problems ? Part II: Special Cases, Implementation and Testing?, Submitted May 2007, Advances in Water Resources. vi ACKNOWLEDGMENTS I would like to express my sincere gratitude to my advisor Dr. T. Prabhakar Clement for his valuable guidance and constant encouragement during my graduate studies here at Auburn University. I would like to thank my committee members Dr. Mark O. Barnett and Dr. Dongye Zhao for their support in the preparation of my thesis. I would also like to thank Dr. Kang-Kun Lee, Dr. Rein van Genuchten, Dr. Joel Melville, Dr. Jacob Dane and Dr. Willie Harper and for their guidance throughout my research work. Funding for this research work was, in part, supported by the Korean Research Council?s Frontier project, funded through the Seoul National University, Seoul, Korea, and by the Office of Science (BER), U.S. Department of Energy Grant no. DE-FG02- 06ER64213. vii Style manual of journal used Advances in Water Resources Computer software used Microsoft Office Word 2003 viii TABLE OF CONTENTS LIST OF TABLES............................................................................................................. xi LIST OF FIGURES .......................................................................................................... xii I. INTRODUCTION........................................................................................................... 1 1.1 Background......................................................................................................... 1 1.2 Literature Review................................................................................................ 2 1.3 Scope and Objective ........................................................................................... 7 1.4 Organization........................................................................................................ 7 II. GOVERNING EQUATIONS AND SOLUTION DERIVATION................................ 9 2.1 Governing Equations .......................................................................................... 9 2.2 Derivation of the Solution for the Dirichlet Boundary condition..................... 10 2.3 Derivation of the Analytical Solution for the Cauchy Boundary Condition .... 27 2.4 Discussion......................................................................................................... 30 III. TESTING AND VALIDATION OF THE ANALYTICAL SOLUTIONS ............... 33 3.1 Introduction....................................................................................................... 33 3.2 Implementation of the Solution ........................................................................ 34 3.3 Special Cases .................................................................................................... 38 3.3.1 Zero Initial Condition ............................................................................... 38 3.3.2 Identical Retardation Factors .................................................................... 39 3.3.3 Zero Advection Velocity........................................................................... 40 ix 3.3.4 Steady State............................................................................................... 42 3.3.5 Zero Dispersion Coefficient...................................................................... 44 3.4 Example Problems ............................................................................................ 46 3.5 Salient Features of the General Solution .......................................................... 55 3.6 Limitations of the General Solution.................................................................. 56 IV. INVESTIGATION OF ERROR ASSOCIATED WITH THE DOMENICO APPROXIMATION ......................................................................................................... 59 4.1 Introduction....................................................................................................... 59 4.2 Governing Equations ........................................................................................ 61 4.3 Review of the Domenico Solution.................................................................... 62 4.4 A Rigorous Approach to Derive the Domenico Solution ................................. 67 4.5 Mathematical Analysis of the Validity of the Approximation Involved in the Domenico Solution ...................................................................................................... 69 4.6 Analysis of the Error Associated With the Domenico Solution ....................... 71 4.7 Example Problem.............................................................................................. 73 4.7.1 Plume Comparison Analysis Behind the Advective Front ....................... 73 4.7.2 Plume Comparison Analysis Beyond the Advective Front ...................... 78 4.8 Discussions and Recommendations.................................................................. 84 V. SUMMARY AND CONCLUSIONS .......................................................................... 88 REFERENCES ................................................................................................................. 91 APPENDIX A................................................................................................................... 99 APPENDIX B ................................................................................................................. 113 APPENDIX C ................................................................................................................. 118 x APPENDIX D................................................................................................................. 127 APPENDIX E ................................................................................................................. 140 APPENDIX F.................................................................................................................. 154 APPENDIX G................................................................................................................. 162 APPENDIX H................................................................................................................. 182 APPENDIX I .................................................................................................................. 184 APPENDIX J .................................................................................................................. 198 APPENDIX K................................................................................................................. 200 xi LIST OF TABLES Table 1: Parameters used in the Cho [9] example problem.............................................. 47 Table 2: Parameters used in the modified Higashi and Pigford [24] example problem as given by van Genuchten [45]............................................................................................ 50 Table 3: Parameters used in the new 10 species example problem .................................. 52 Table 4: Summary of parameter limitations ..................................................................... 58 Table5: Parameters used in the Domenico and Robbins [17] example problem.............. 74 xii LIST OF FIGURES Fig. 1: Comparison of the concentration profiles of the general solution with the Cho [9] solution (3 Species, Dirichlet boundary condition) .......................................................... 48 Fig. 2: Comparison of the concentration profiles of the general solution with the van Genuchten [45] solution (4 Species, Cauchy boundary condition) .................................. 51 Fig. 3: Comparison of the concentration profiles of the general solution with the Quezada et al. [38] semi-analytical solution (10 Species, Dirichlet boundary condition) .............. 53 Fig. 4: Comparison of the concentration profiles of the general solution with the Quezada et al. [38] semi-analytical solution (10 Species, Cauchy boundary condition) ................ 54 Fig. 5: Concentration contours predicted by the Domenico and Wexler solutions for the base case: solutions behind the advective front for (a) X-Y plane (b) X-Z plane ............ 75 Fig. 6: Sensitivity results for variations in the longitudinal dispersivity value; solutions behind the advective front for (a) ?x*10 (b) ?x/10 .......................................................... 77 Fig. 7: Sensitivity results for variations in the transport velocity; solutions behind the advective front for (a) v*10 (b) v/10............................................................................... 79 Fig. 8: Concentration contours predicted by the Domenico and Wexler solutions; solutions include concentration contours beyond the advective front.............................. 80 Fig. 9: Sensitivity results for variations in the longitudinal dispersivity value; solutions include concentration contours beyond the advective front for (a) ?x*10 (b) ?x/10....... 82 xiii Fig. 10: Sensitivity results for variations in the transport velocity; solutions include concentration contours beyond the advective front for (a) v*10 (b) v/10 ....................... 83 Fig. 11: Concentration profiles predicted by the Domenico solution compared with the Wexler solution at x = 1000m and 1500m (a) along Y axis (b) along Z axis................... 86 1 CHAPTER I INTRODUCTION 1.1 Background Transport problems involving sequentially decaying contaminants are frequently analyzed by groundwater hydrologists to assess water quality issues associated with environmental and health hazards. Examples of sequentially decaying contaminants include radioactive waste materials, chlorinated solvents, and nitrogenous species [4, 11, 45]. Several types of models, using both analytical and numerical procedures, have been formulated for solving these sequentially coupled reactive transport problems [12, 26]. Although numerical models are capable of solving complex and heterogeneous problems, their performance often needs to be tested against experimental datasets or analytical models. Experimental simulations of complex reactive transport problems are not only time consuming but can also be expensive. Therefore, analytical models provide a convenient, cost-effective alternative to test and validate numerical formulations [18, 30, 44]. Furthermore, analytical models also provide computationally efficient screening tools for simulating the fate and transport of reactive contaminants in groundwater systems [3, 13]. 2 1.2 Literature Review The analytical solution given by McLaren [32] and McLaren [33], which describes the steady-state, one-dimensional transport of a five species nitrogen chain, is one of the first multi-species solutions derived for solving sequentially coupled reactive transport problems. This work assumed that the transport was only governed by advection, and the effects of dispersion and sorption were ignored. Cho [9] developed explicit analytical solutions to a three species transport problem that was subjected to advection, dispersion, linear equilibrium sorption, coupled through sequential first-order reactions. Explicit analytical solutions were obtained using Laplace transform procedures for the Dirichlet boundary condition. One of the limitations of this solution is that only the first species in the chain was subjected to sorption. Later, Misra et al. [34] derived semi-analytical solutions to a problem similar to Cho [9] using a pulse source boundary. Burkholder and Rosinger [7] and Lester et al. [29] developed solutions for the advective dispersive transport of radionuclide chains subjected to linear equilibrium sorption. Explicit analytical solutions were presented for a three species problem involving distinct retardation factors for each species, for both impulse and decaying- band release boundary conditions. In addition, they also provided solutions for the case of no dispersion and for the case of identical retardation factors. Harada et al. [23] published a research report presenting general semi-analytical solutions to sequentially coupled one dimensional reactive transport problems of arbitrary chain lengths subjected to arbitrary release modes. However, one of the major limitations of the solution strategy was that, the semi-analytical solution for a given species in the 3 chain required the computation of its entire predecessor species. This would result in computationally inefficient algorithms especially when analyzing transport problems involving long reactive chains. Harada et al. [23] and Higashi and Pigford [24] also provided explicit closed-form solutions for a set of purely advective (no dispersion) transport problems with various types of boundary conditions. Gurehian and Jansen [21] presented an analytical solution to a transport problem involving a three member, first-order decay chain in a multi-layered system, subjected to advection, dispersion and linear equilibrium sorption processes for both continuous and band source release conditions. Convolution theorems and Laplace transform techniques were used to obtain semi-analytical solutions for the case involving both advective and dispersive transport, and explicit closed-form analytical solutions for the case involving non-dispersive transport. van Genuchten [45] developed explicit analytical solutions to model a sequentially coupled four species transport problem governed by advection, dispersion and linear equilibrium sorption processes involving, first-order reactions. It was assumed that all the species had distinct retardation factors. One of the key contributions of this work is that it considered both Dirichlet and Cauchy boundary conditions. Furthermore, van Genuchten [45] developed a robust computer code (CHAIN) for implementing his analytical solution. Angelakis et al. [2] developed a semi-analytical solution to a sequentially coupled two-species transport problem governed by advection, dispersion and linear equilibrium sorption subjected to Dirichlet boundary condition. The transport problem assumed that the reactions were first-order and each of the species had different dispersion coefficients 4 and distinct retardation factors. The authors also demonstrated that when the dispersion coefficients of both the species were equal, their solution reduced to the closed-form solution similar to the solutions presented by Cho [9] and Misra et al. [34]. Furthermore, the authors also provided solutions for the no dispersion (pure advection) case. Angelakis et al. [1] developed an interesting semi-analytical solution for a problem involving the coupled transport of two solutes and a gaseous product in soils. The solute migration was governed by advection, dispersion, linear equilibrium sorption and sequential first-order reaction, whereas the gas migration was governed by diffusive transport coupled with reversible linear equilibrium dissolution. Lunn et al. [30] solved a three-species transport problem, which was similar to the Cho [9] problem, using the Fourier transform method. The authors demonstrated that the use of Fourier transforms enabled them to solve problems having non-zero initial conditions by solving two special case problems. Khandelwal and Rabideau [27] developed semi-analytical solutions for a three species, sequentially-coupled, first-order reactive transport problem. The key contribution of this work was that they addressed cases involving linear, non-equilibrium sorption mechanisms. Eykholt and Li [18] developed a solution method based on kinetic response functions to solve a linearly coupled non-sequential reactive transport problem having different retardation factors. Although, there was no restriction on the number of species in the system, this method required numerical procedures to evaluate the final solution. Furthermore, for the case of the non-ideal plug flow scenarios (advective dispersive 5 transport), the accuracy of this method appears to decrease with decrease in Peclet number. Sun et al. [42] developed a method that can solve multi-species advective dispersive transport equations coupled with sequential first-order reactions involving arbitrary number of species for different types of initial and boundary conditions. Their method was based on the use of a transformation procedure to uncouple the system of equations, which could then be solved analytically in the transformed domain. The final solutions are obtained by retransforming the solutions to the original domain. Later, Sun et al. [43] extended the transformation format to solve problems involving a combination of serial and parallel reactions. Clement [11] presented a more general and fundamental approach to derive the Sun et al. [43] solution by employing the similarity transformation method. The approach presented by Clement [11] can also be used to solve problems involving serial, parallel, converging, diverging and/or reversible first-order reaction network. However, all of these methods are only applicable for solving problems involving identical retardation factors. Bauer et al. [4] presented a method to solve one, two, and three-dimensional sequentially coupled reactive transport problems with distinct retardation factors. This method was based on transforming the system of equations to a Laplace domain and then obtaining a set of fundamental solutions to each of the equations in the transformed domain. The specific solutions in the Laplace domain can then be obtained through a linear combination of the fundamental solutions provided the fundamental solutions are linearly independent. Finally, the Laplace domain solutions can be transformed back to the time domain using the inverse Laplace transform procedure, which could be 6 accomplished either analytically or numerically. Although this method can be applied to solve different types of boundary conditions, the solution procedure is mathematically tedious; specifically obtaining analytical inverse transform expressions for long chain lengths can be a challenge. Montas [35] developed an analytical procedure to solve a three species, multi- dimensional transport problem coupled by a first-order, non-sequential reaction network subject to a pulse type boundary condition. This procedure involved obtaining a basis solution of a convoluted form for the transport equation and then evaluating the basis solution using Laplace transforms. One of the key advantages of this procedure is that it can model transport problems with distinct retardation factors. However, as mentioned earlier this solution was limited to a three species system. Quezada et al. [38] extended the approach given by Clement [11] and developed a method that can solve multi-species transport equations coupled with a network of first- order reactions involving distinct retardation factors. This method involves transforming the system of governing equations to a Laplace domain and then solving the transformed system of equations using the Clement [11] approach. The solutions in the Laplace domain are then retransformed to the time domain using an inverse Laplace transform procedure. One of the key limitations of this approach is that, except for a simple two species transport problem, the solutions are in general semi-analytical since they require a numerical inverse Laplace transform routine to evaluate them. 7 1.3 Scope and Objective The above literature review indicates that one-dimensional reactive transport equations coupled through sorption and sequential first-order reactions have explicit closed-form analytical solutions only for short chains up to four species. To model transport problems involving longer reaction chains, one has to either use semi-analytical solutions or purely numerical solutions. This objective of this work is to develop a general closed-form analytical solution to the sequential transport problem involving arbitrary number of species subjected to a generic exponentially decaying Bateman-type source boundary, for a spatially varying initial condition. The solutions derived in this study are then implemented in a computational platform and tested and validated against other analytical and semi- analytical solutions published in the literature. 1.4 Organization This thesis is organized into five chapters. Chapter I presents a brief introduction with a detailed literature review of analytical solutions to coupled reactive transport problems. Chapter II details the solution derivation of the sequentially coupled one dimensional reactive transport problem. Chapter III discusses some of the techniques in implementing the new solution on a computational platform and testing the new solution through a set of example problems. Chapters II and III are draft versions of manuscripts that have been submitted to Advances in Water Resources journal. 8 Chapter IV investigates the validity of approximations involved in extending the one-dimensional solutions derived in this study into three-dimensions. The contents of this chapter are extracted from the author?s publication titled ?Domenico Solution ? Is it Valid?? in the Groundwater journal. Chapter V summarizes the key findings of this study, and also discusses some recommendations for future work. 9 CHAPTER II GOVERNING EQUATIONS AND SOLUTION DERIVATION 2.1 Governing Equations Consider, a one-dimensional transport problem involving ?n? sequentially decaying contaminants simultaneously subjected to advection, dispersion and linear adsorption processes. The general governing equation for this transport problem can be expressed as: ( ) ( ) ( ) ( ) ( ) ( ) 2 1 12 , , , , , ; 2, 3,... , 1 i i i i x i i i i i i i c x t c x t c x tR v D y k c x t k c x t i n t x x k c x t ; i ? ? ? ? ?+ ? = ? ? = ? ? ? = ? = 0 0t and x ; ? > < < ? (1) where; ? ic ? is the concentration of species i [ML-3]; ? iR ? is the retardation coefficient of species i [Dimensionless];? iy ? is the effective yield factor that describes the mass of a species i produced from species i-1 [MM-1]; ? ik ? is the first-order decay rate constant of species i [T-1]; ?v ?is the transport velocity [LT-1]; ? xD ? is the dispersion coefficient [L2T-1]; and ?n? is the total number of species in the reaction network. Equation(1) is solved for a generic exponentially distributed initial condition given by: ( ),0 , 1, 2,...ixoi ic x c e x i n??= 0 < < ? ; ? = (2) 10 where; ? oic ? is the initial source concentration of species i [ML-3]; ? i? ? is the first-order decay parameter of the initial distribution of species i [L-1]. The boundary condition at ?? ? is given as: ( ), 0, 0 1, 2,...ic t t i nx? ? = > ; ? = ? (3) Explicit solutions including detailed derivation steps are provided for the following two inlet (source) boundary conditions: Dirichlet (Section-3) and Cauchy (Section-4) boundaries. 2.2 Derivation of the Solution for the Dirichlet Boundary condition For the case of the Dirichlet boundary, the boundary condition at the source is described as follows: ( ) 1 1 1 1 , 00, 1, 2,... 0, i i t i i o ii o B e t tc t i n t t ?? = ? < ?? = ; ? = ? ? >? ? (4) where; ? 1iiB ? is the source boundary concentration of specie i1 that contributes to species i [ML-3]; ? 1i ? ? is the first-order decay of the corresponding ? 1iiB ? term [T-1]. Equation(4) can be conveniently re-written as: ( ) ( ) ( ){ } ( ) 1 1 1 1 0, , 0 1, 2,... ; 0, 1, i i t i i i o i c t B e u t u t t t i n where u istheunit step function givenby if t au t a if t a and ais an arbitrary positiveconstant ?? = = ? ? > ; ? = ? ?? ? ? (8) ( ){ } { } ( ) ( ){ }1 1 1 1 0, , 0 ; , 0 1, 2,...i i t i i i o i c t t where B e u t u t t t i n?? = = ? > ? = ? ? > ; ? = ? (9) The solution procedure used here is adopted from Quezada [38]. Applying Laplace transform to equation(6), we get: [ ] { } [ ]{ } { } { } [ ]{ } 2 2( ,0) x d p d pR s p R c x v D K p dx dx? + ? = (10) where; ? s ? is the Laplace variable and ? p ? is the Laplace transformed concentration. Substituting equation(7) in (10) and rearranging we get: { } { } [ ] [ ]( ){ } [ ]{ } 2 2 1 1 o x x x x d p d pv K R s p R c e dx D dx D D ??? ? ?? + ? =? ? ? ? (11) Now in order to uncouple the system of ordinary differential equations (ODEs) given by equation(11), we apply the linear transform procedure described by Clement [11], by performing the following matrix operation. 12 { } [ ]{ }p A b= (12) where; ?{ }b ? is the concentration in the doubly transformed domain and [ ]A is an arbitrary square matrix of order n. Applying this transformation equation(11) gets modified as: [ ]{ } [ ]{ } [ ] [ ]( )[ ]{ } [ ]{ } 2 2 1 1 o x x x x A b A bv K R s A b R c e x D x D D ??? ?? ? ?? + ? =? ? ? ?? ? (13) Pre-multiplying equation(13) with ?[ ] 1A ? ? we get: { } { } { } { } [ ] [ ] [ ]( )[ ] { } [ ] [ ]{ } 2 ~ ~ 2 ~ ~1 1 1 1 ; x x x o x b bv K b C x D x D D where K A K R s A and C A R c e? ? ?? ? ?? ? ?? ?? + = ? ? ? ?? ? ? ?? ? ? ? = ? = ? ?? ? (14) By forcing the columns of the ?[ ]A ? matrix as the eigenvectors of the combined reaction coefficient matrix ? [ ] [ ]K R s? ? we can make the ? ~ K? ?? ?? ? ? matrix a diagonal matrix and thus uncouple the system of equations; the details of this similarity transformation procedure are illustrated in Clement [11]. The corresponding ?[ ]A ? matrix is: 13 [ ] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 2 1 1 1 2 2 3 31 1 1 1 2 2 1 1 1 1 1 ,0,0,..., 1 , ,0,0,..., 2 . . , ,..., , n i i i i i n n i i i i i ii i i i n n n i i i i n n i i i n i n i ni i i i i i k sR k sR n times y k k sR k sR k sR k sR n times y k y k A k sR k sR k sR k sR k sR k sR y k y k y k = ? = =? ? ? ? = = =? ? ? ? + ? ? ? ? + ? ? ? + ? ? ? = ? + ? ? ? + ? ? ? + ? ? ? ? ? ? ? ? 0 1,1,...,ntimes ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (15) The ?[ ] 1A ? ? matrix is: 14 [ ] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 1 1 1, 1 1 1 2 3 2 2 2 21 1, 2 2, 2 1 1 2 3 1, ,0,0,..., 1 , ,0,0,..., 2 . . , n i i i n i i i i n n i i i i i i n n i i i i i i i i n n i i i i i i n n n i i n n i i i n y k n times k sR k sR y k y k n times k sR k sR k sR k sR A y k y k k sR k sR k sR k ? = = ? ? ? = = ? = ? = ? ? ? = = = ? ? ? + ? ? ? ? + ? ? ? + ? ? = ? + ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ( ) ( ) ( ) ( ) 1 2, 1, ,..., ,1 n i i i n n n i n n i i i i n i n i n y k sR k sR k sR ? = = ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? + ? ? ? ? ? ? ? ?? ? ? ? ? (16) 15 The corresponding ? ~ K? ?? ?? ? ? matrix is: ( ) ( ) ( ) ( ) ( ) 1 1 2 2 ~ ,0,0,... 1 0, ,0,0,... 2 . 0,0,... 1 , ,0,0,... . 0,0,... 1 , i i n n k sR n times k sR n times K i times k sR n i times n times k sR ? ? ?? ? ? ?? ? ? ? ? ? ?? ? ? ?=? ? ? ?? ? ? ?? ? ? ? ? ? ? ?? ? ?? ? (17) The corresponding ?{ } ~ C ? vector is: { } ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 2 1 2 1 1 1 2 1 1 1, 1 1 1 1 2 2 1 2 3 ~ 2 2 2 2 1, 2 2, 2 1 1 1 2 2 1 2 1, . . n xo i i i n i i i i n n x xo o i i i i i i n n i i i i i i i i n x xo o i i i i i i n n n i i i i n Rc e y k k sR k sR Rc e y k R c e y k k sR k sR k sR k sR C Rc e y k R c e y k k sR k sR ? ? ? ? ? ? ? = = ? ? ? ? ? = = = ? = ? ? ? ? ? = = ? ? + ? ? + ? + ? ? ? + ? ? = + ? + ? ? ? ? ? ? ? ? ? ? ( ) ( ) 3 2, ... n n xo n nn n n i i i i n R c e k sR k sR ??= = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?+ + ? ?? + ? ? ? ? ? ? ? ?? ? ? ? (18) The explicit expression for ? ~ iC ? in equation(18) is: 16 ( ) ( ) 1 1 1 2 2 2 1 1 2 2 2 1 2 1~ 1 1 , 1, 2,... i nx o i i i ii i i i n i i i i i i i i i R c e y k C i n k sR k sR ?? ? = + = = ? ? ? ? ? ? ?= ; ? = ? ?? + ? ? ? ?? ? ? ? ? (19) Equation(14) describes a set of ?n ? independent second-order non-homogeneous ODEs the boundary conditions of which are obtained by performing Laplace and linear transforms of the boundary conditions given by equations(8) and (9). Laplace transform of equations(8) and (9) yields: ( ), 0p sx? ?? ? =? ?? ? ? (20) ( ){ } { } ( ){ } ( ) 11 1 11 0, 1 ; 1, 2,... o it si i i i i i p s B e where i ns ? +? = = ? ? ? = ; ? = +?? (21) To transform the boundary conditions from ? p ? domain to the ?b? domain, we apply the linear transform given by equation(12). This yields: ( ), 0b sx? ?? ? =? ?? ? ? (22) ( ){ } [ ] { }10,b s A ?= ? (23) The explicit expression for ? ( )0,ib s ? in equation(23) is given as: ( ) ( ) ( ) ( ){ } ( ) 2 2 2 1 2 2 2 1 2 221 1 1 1 1 2 2 , 1 1 1 0, 1, 2,... n o i i i i i n i i i i i i i i t si ii i i i i i y k k sR k sR B e b s s i n ? = + = ? ? +? = =+ ? ? ? ?? ? ? ? ?= +?? ?? ? ?? ? ; ? = ? ? ? ? (24) 17 Since equation(14) is uncoupled, it can now be written as a set of ?n ? independent equations as: ( ) ( ) ( ) ( ) ( ) ( ) 1 1 1 2 2 2 1 1 2 2 2 1 2 12 1 2 1 , , , 1 1, 1, 2,.. i nx o i i i ii i ii i i i i n ix x x i i i i i i i i R c e y kb x s b x sv k sR b x sx D x D D k sR k sR i ?? ? = + = = ? ? ? ? ?? ?? ? ? ? ?? + ? ? =? ? ? ?? ?? ? ? + ? ? ? ?? ? ; ? = ? ? ? .n (25) The general solution to equation(25) is given as: ( ) ( ) ( ), , , 1, 2,...h pi i ib x s b x s b x s i n= + ; ? = (26) where; ? ( ),hib x s ? is the general solution of the homogeneous part of equation(25) and ? ( ),pib x s ? is a particular solution of equation(25). The general solution ? ( ),hib x s ? can be readily obtained as: ( ) ( ) ( )2 2 2 2 4 4 2 21 2, 1, 2,... i i i i x x x xx x k sR k sRx v v x v v D D D DD Dh i i ib x s e e i n ? ? ? ?? ? ? ?+ +? ? ? ? ? ? ? ?+ + ? +? ? ? ? ? ? ? ?? ? ? ?? ? ? ?? ? ? ?= ? + ? ; ? = (27) where; 1i? and 2i? are constants. The particular solution ? ( ),pib x s ? is obtained by using the method of undetermined coefficients. The general form of the particular solution is given as: ( ) ( ) ( ) 1 1 1 1 2 2 2 1 1 2 2 2 1 2 1 1 1 , 1, i nx o i i i i ii i ip i n ix i i i i i i i i M R c e y k b x s D k sR k sR ?? ? = + = = ? ? ? ? ?? ? ?= ? ?? + ? ? ? ?? ? ? ? ? (28) 18 where; ? 1i M ? is a constant. Substituting equation(28) in the governing equation(25) and simplifying, we evaluate the constants ? 1i M ? as: ( )1 1 12 1, 2,... x i i x i i i DM D v k sR i n = ? +? ? ? ; ? = (29) Substituting the values of ? 1i M ? into equation(28) we get the particular solution ? ( ),pib x s ? to equation(25) as: ( ) ( ) ( ) ( ) 1 1 1 2 2 2 1 1 1 1 2 2 2 1 2 1 1 21 , , 1, 2,... i nx o i i i ii i ip i n i i x i i i i i i i i i i i R c e y k b x s D v k sR k sR k sR i n ? ? ? ? ? = + = = ? ? ? ? ? ? ?= ? ? ?+ ? ? ? + ? ? ? ?? ? ; ? = ? ? ? (30) Substituting equations(27) and (30) in equation(26), we get the general solution to equation(25) as: ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 1 1 1 2 2 2 1 1 1 1 2 2 2 1 2 4 4 2 21 2 1 1 21 , , i i i i x x x xx x i k sR k sRx v v x v v D D D DD D i i i nx o i i i ii i i n i i x i i i i i i i i i i i b x s e e R c e y k D v k sR k sR k sR ? ? ? ? ? ? ?? ? ? ?+ +? ? ? ? ? ? ? ?+ + ? +? ? ? ? ? ? ? ?? ? ? ?? ? ? ?? ? ? ? ? ? = + = = ? = ? + ? ? ? ? ? ? ? ? ? ?+ ? ? ? + ? ? ? ?? ? ? ? ? 1, 2,...i n ; ? = (31) In order to apply the boundary condition given by equation(22) we differentiate the general solution with respect to x. Differentiation of equation(31) with respect to x yields: 19 ( ) ( ) ( ) ( ) ( ) 2 2 2 2 42 21 2 42 22 2 4, 1 2 41 2 i i x xx i i x xx k sRx v v D DDi ii i x x x k sRx v v D DDi i i x x x k sRb x s v v e x D D D k sRv v e D D D ? ?? ?+? ? ? ?+ +? ? ? ?? ?? ?? ? ? ?? ?+? ? ? ?? +? ? ? ?? ?? ?? ? ? ?? ?+? ? ? ? ?= ? + +? ?? ? ?? ?? ?? ? ? ?? ?+? ? ? ? + ? ? +? ? ? ?? ?? ?? ? ( ) ( ) ( ) ( ) 1 1 1 1 2 2 2 1 1 1 1 2 2 2 1 2 1 1 21 , 1, 2,... i nx o i i i i ii i i n i i x i i i i i i i i i i i R c e y k D v k sR k sR k sR i n ?? ? ? ? ? = + = = ? ? ?? ? ? ? ?? ? ?+ ? ? ? + ? ? ? ?? ? ; ? = ? ? ? (32) To satisfy the boundary condition given by equation(22), i.e. when ? x ? tends to ???; the exponential function in the first term tends to ???, hence 1i? must vanish. Equation(32) now reduces to: ( ) ( ) ( ) ( ) ( ) 2 2 1 1 1 2 2 2 1 1 1 1 2 2 2 1 2 4 22 1 1 21 , , i i x xx i k sRx v v D DD i i nx o i i i ii i i n i i x i i i i i i i i i i i b x s e R c e y k D v k sR k sR k sR ? ? ? ? ?? ?+? ? ? ?? +? ? ? ?? ?? ?? ? ? ? = + = = ? = ? ? ? ? ? ? ? ? ? ?+ ? ? ? + ? ? ? ?? ? ? ? ? 1, 2,...i n ; ? = (33) Applying the second boundary condition given by equation(24), we get: 20 ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) 2 2 2 1 2 2 2 1 2 221 1 1 1 1 2 2 , 1 1 2 2 2 1 1 1 1 2 2 2 1 2 2 1 1 1 1 21 , 1 1 n o i i i i i n i i i i i i i i t si ii i i i i i n o i i i ii i i n i i x i i i i i i i i i i i y k k sR k sR B e s R c y k i D v k sR k sR k sR ? = + = ? ? +? = = ? = + = = ? + ? ? ? ?? ? ? ? ?? = +?? ?? ? ?? ? ? ? ? ? ? ? + ; ? = ? ?? +? ? ? ? + ? ? ? ?? ? ? ? ? ? ? ? ? , 2,...n (34) Therefore, the solution in the ?b?domain is: ( ) ( ) ( ) ( ){ } ( ) ( ){ } ( ){ } 2 2 2 1 2 2 2 1 2 22 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 , 42 1 1 42 1 1 1 , o i x i i x x i i ix n i i i i n i i i i i i i i t si x ii v v D k sRi D i i i i xn v v D k sR xDo i i i i i i y k k sR k sR B e b x s es R c y k e e ? ? ? ? = + = ? ? + ? ? ? + +? ? ? ? = = ? ?? + + ? ? ?? ? ? = + + ? ? ? ?? ? ? ? ?= + ? ?? ? ?? ? ? ?? ? ?? ? ?? + ? ? ? ? ? ( ) ( ) ( ) 1 1 1 2 2 2 1 2 21 , 1, 2,... i n i i x i i i i i i i i i i i D v k sR k sR k sR i n ? ?= = ? ? ? ? ? ? ??? ? ? ? ?+ ? ? ? + ? ? ? ? ? ?? ? ; ? = ? ? (35) Inverse linear transform of equation(35) is done to obtain the solution in the Laplace domain (? p ? domain) by using equation(12). The solution given by equation(35) can be split into two parts and represented as: 21 ( ) ( ) ( ) ( ) ( ) ( ) ( ){ } ( ) ( ){ } ( ) 2 2 2 1 2 2 2 1 2 22 2 1 1 1 2 2 2 1 1 2 2 2 1 1 1 , 1 2 421 1 1 42 1 12 , , , ; 1 , , o i x i i x x i x n i i i i n i i i i i i i i i i i t si x ii v v D k sRi D i i i i xn v v D k sR Do i i i i i i i y k k sR k sR b x s b x s b x s where B e b x s es R c y k e b x s ? ? ? = + = ? ? + ? ? ? + +? ? ? ? = = ? + + ? = + + ? ? = + ? ?? ? ? ? ?= + ? ?? ? ?? ? = ? ? ? ? ? ( ){ } ( ) ( ) ( ) 1 1 1 1 2 2 2 1 2 21 , 1, 2,... i i x i n i i x i i i i i i i i i i i e D v k sR k sR k sR i n ? ? ? ? ? ? ? ?? ? = = ? ? ?? ?? ? ? ??? ? ? ?? ?? ? ? ? ? ?+ ? ? ? + ? ? ? ? ? ?? ? ; ? = ? ? (36) Using the distributive property of matrix addition, we can apply the inverse linear transform to each of the individual terms and then sum them to get the solution in the ? p ? domain. This is expressed as: { } [ ]{ } [ ]{ } [ ]{ }1 2p A b A b A b= = + (37) The first term ?[ ]{ }1A b ? can be evaluated as: { } [ ]{ }1 1p A b= (38) The explicit expression for ? ( )1 ,ip x s ? in equation(38) is given as: 22 ( ) ( ){ } ( ) ( ){ } ( ) ( ) 22 1 1 2 2 22 1 2 2 2 2 1 2 1 2 2 3 3 3 1 3 2 1 11 1 42 1 , 1 , o i x i i x t si ii i i i ii i ii xi v v D k sR Dii i i i i i i i i i i i B e y k s p x s e k sR k sR ? ? ? + ? == + ? ?? + + ? ?? ?= = = ? ? ??? ? ? ?? ? ? ?+? ? ? ? ? ?= ? ? ? ? ? ?? + ? ? ? ?? ? ?? ? ? ? 1, 2,...i n ; ? = (39) Using a similar approach the second term ?[ ]{ }2A b ? is evaluated and the explicit expression for ? ( )2 ,ip x s ? is: ( ) ( ){ } ( ) ( ) ( ) 1 1 2 2 2 1 2 2 2 1 1 2 1 1 1 2 2 2 2 3 3 3 1 3 2 1 1 42 2 1 2 , , x i i ix i o i i i i i i xi v v D k sR xD i i i i i i i x i i i i i i i i i i i R c y k p x s e e D v k sR k sR k sR ? ? ? ? = + ? ?? + + ? ? ?? ? = = = ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ?= ?? ? ? ?? ? ? ?? ? ? ?+ ? ? ? + ? ? ? ? ? ?? ? ? ? ? ? 1, 2,...i n ; ? = (40) Substituting equations(39) and (40) into equation(37) we get the solution in the Laplace domain as: 23 ( ) ( ) ( ) ( ){ } ( ) ( ){ } ( ) ( ) 22 1 1 2 2 22 1 2 2 2 2 1 2 1 2 2 3 3 3 1 3 2 1 1 1 2 1 11 42 1 , , , , 1 o i x i i x i i i t si ii i i i ii i ii x v v D k sR Dii i i i i i i i i i i i o i i p x s p x s p x s B e y k s e k sR k sR R c y ? ? ? + ? == + ? ?? + + ? ?? ?= = = ? = + = ? ??? ? ? ?? ? ? ?+? ? ? ? ? ? ? ? ? ? ? ?? + ? ? ? ?? ? + ?? ? ? ? ( ){ } ( ) ( ) ( ) 2 2 2 1 2 2 2 1 1 2 1 1 1 2 2 2 2 3 3 3 1 3 2 1 1 42 1 2 , x i i ix i i i i i xi v v D k sR xD i i i i i i x i i i i i i i i i i i k e e D v k sR k sR k sR ? ? ? ? = + ? ?? + + ? ? ?? ? = = = ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ??? ? ? ?? ? ? ?? ? ? ?+ ? ? ? + ? ? ? ? ? ?? ? ? ? ? ? 1, 2,...i n ; ? = (41) The final solution is obtained by taking an inverse Laplace transform of the solution given by equation(41). Inverse Laplace transform is performed as follows: ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 1 1 2 1 1 1 2 , , , , , , , 1, 2,... i i i i i i i c x t c x t c x t p x s p x s p x s p x s i n ? ? ? = + = + = + ; ? = ? ? ? (42) In Appendix A, the terms ? ( )1 1 ,ip x s?? ? and ? ( )1 2 ,ip x s?? ? are evaluated. Substituting equations (A.7) and (A.16) in equation(42) we obtain the final solution in the time domain as: 24 ( ) ( ){ } ( ){ } 1 2 2 1 2 1 32 1 1 1 2 2 1 2 12 1 1 1 1 1 1 1 2 1 11 2 2 2 1 1 1 2 1 1 , iii i i i i i i i ii i ii i o i i i i i i ii i c x t y k G h G G R c y k G h G G ? = = == + ? = == + ? ?? ?= + ? ?? ?? ? ? ?? ? ? ?? ? + + ? ?? ?? ? ? ?? ? ? ??? ? ?? 1, 2,...i n ; ? = (43) where; the ?G ? terms are defined as (See equations (A.7) and (A.16) in Appendix A). 25 [ ] ( ) ( ) ( ) [ ] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 3 2 3 2 33 1 3 2 2 4 2 2 4 4 2 4 3 2 5 2 4 2 5 2 5 2 4 5 1 5 2 5 1 5 2 5 45 2 5 2 , ,0 , ,0 , , , ,1 1 , , , , , , , , , , , , , , , oi oi i i i i t i i o i i oi i t i i i o i i i o i ii i i i i i i i i i i i i i i i i i i R R i i i i i i R R F x t u t t e F x t t B F x t u t t e F x t tG a k R R a a ? ? ? ? ? = = ? = = ? ? ? ? ?? ? ?? ? ? ?? + ? ?? ?= ? ? ? ?? ? ? ? ?? ? ? ? ? ?( ) [ ] ( ) ( ) ( ) ( ) [ ] ( ) [ ] ( ) ( ) ( ) 1 4 2 4 2 33 1 2 3 2 3 2 4 4 1 4 2 4 2 ,, 1 2 31 1 2 2 1 2 2 2 3 2 3 1 2 2 2 4 4 1 4 2 4 2 , , , ,0 , ,0 1 2 , , , , , , ,2 1 , , , , , , , , , i i oi i i i i ii i i i i i i i R R ti i i i o i i o i i i i i i i R R x a tx a t i i i i i i i i i i i i i i i i i R R B F x t u t t e F x t t G k F x t e F x t e G a a R k ? ?? ? ? ? ? = ? = ? ?? ? ? ? = ? = ? ?? ? ?? ? = ? ? ? + = ? ? ? ? ( ) ( ) ( ) [ ] ( ) ( ) 3 1 3 2 3 2 2 3 2 4 2 4 2 3 4 1 4 2 4 3 4 2 ,1 1 2 2 1 2 2 2 3 3 1 3 2 3 2 , , , , , , , , , , ,2 2 , , , , i i i i i i i i i i i ii i i i R R i i i i i i i i i i i i i i R R x a t i i i i i i i i i i i R R R R a a F x t e G R k ? ? = ? ? = ? ? ? ? ? ? = ? = ? ? ? ? ? ?? ? ? ? ? ? = ? ? ? ? ? (44) 26 where; the term ? 1 2 3, ,i i i F ? is given by: (See equation(B.8) in Appendix B) [ ] , ,1 2 3 1 1 2 3 ,2 3 1 1 2 3 , ,1 2 3 1 1 2 3 1 1 1 2 3 1 2 3 1 , ,2 2 , , , ,2 2 , , , 2 , 2 2 ; 4 i i i x i i x i i i x x i i i iD xv a t D i x i i i x i i i iD i x i i i i i x i i i R x te erfc R D te e F x t R x te erfc R D t kwhere v R D a R ? ? ? ? ? ? ? ? ?? ?? ? ? ? ?? ? ? ?? ?? ? ? ?= ? ?? ?+ ? ? ? ?+ ? ? ? ?? ?? ?? ? ? ? = + ? ? ?? ? ? ? (45) The above solution is valid only for real values of ? 1 2 3, ,i i i ? ?. For problems involving complex values for ? 1 2 3, ,i i i ? ? the ? 1 2 3, ,i i i F ? terms are given as: (See equation(B.14) in Appendix B) [ ] ( ) , 1 2 3 1 2 32 3 1 2 3 1 1 2 3 1 2 3 1 1 1 2 3 1 * * , , , ,2 , , * 2 , , , * , , , 2 2 ; 4 2 i i x xv a t i i i i i iD i i i x x i i i i i x i i i i i i i i x x xF x t e e ACos BSin D D kwhere v R D a R R x i tand A iB erfc R D t ? ? ? ? ? ? ?? ? ? ?= ?? ?? ? ? ? ? ? ? ?? ?? ? ? ?? ? ? ? = + ? ? ?? ? ? ? ? ?+ ? ? + = ? ? ? ?? ? (46) 27 2.3 Derivation of the Analytical Solution for the Cauchy Boundary Condition For the case of the Cauchy boundary, the boundary condition at the source is described as follows: ( ) ( ) 1 1 1 1 , 00, 0, 1, 2,... 0, i i t i i oi ix i o B ve t tc tD vc t i n x t t ?? = ? < ?? ? ? + = ; ? = ?? ? >? ? (47) Equation(47) can be conveniently re-written as: ( ) ( ) ( ) ( ){ } ( ) 1 1 1 1 0, 0, , 0 1, 2,... ; 0, 1, i i t i i x i i o i c tD vc t B ve u t u t t t i n x where u istheunit step function givenby if t au t a if t a and ais an arbitrary positiveconstant ?? = ?? + = ? ? > ; ? = ? < < ? (52) Note the additional parameters in the right side of equation(52). The solution to the above equation can be readily obtained from the previous solution given by equation(43) by substituting the value of ?k ? with ? Rk ?. From Sections 2.2 and 2.3, it can be seen that the solutions for the Dirichlet and the Cauchy boundaries share a common structure. However, careful observation indicates that the ? 21G ? and ? 22G ? terms for the Dirichlet boundary (see equation(44)) and the Cauchy Boundary (see equation(49)) are different. Furthermore, from equations (45) and (46) and equations (50) and (51) it can be observed that the ? 1 2 3, ,i i i F ? terms involved in the ?G ? terms are distinctly different for the two boundary conditions. The solutions are presented in a format wherein they share a very similar structure. One of the key advantages of this type of presentation is that the implementation of the solution for one type of boundary condition becomes relatively easy if we have the algorithm for the other boundary implemented. Furthermore, the general solution is presented in a format which enables us to directly obtain explicit solutions for any of the species in the chain without involving the computations of its parent chain members. This unique feature makes the general solution computationally efficient. The solutions previously published in the literature have either been restricted to small chain lengths or have been semi analytical solutions for longer chain lengths. The solutions derived in this study overcome both these difficulties and making them highly suitable for screening level models. 32 In the following chapter, the general solutions are first implemented in a computational platform and then tested against other analytical and semi-analytical solutions using a set of example problems. 33 CHAPTER III TESTING AND VALIDATION OF THE ANALYTICAL SOLUTIONS 3.1 Introduction To test the general solution, one can either compare it analytically against other closed-form solutions or compare the results of its simulation against other analytical or numerical solutions. While testing the general solution using example simulations requires its implementation in a computational platform, analytical comparison is a viable verification option only when we have suitable solutions that solve an identical transport problem. Since explicit closed-form solutions are available only for simplified transport scenarios, the general solution is first simplified to these special cases and then compared against previously published closed-form solutions. The objective of this chapter is to: 1) discuss the computational implementation of the new general analytical solution; 2) test and validate the solution against published analytical and numerical results; and 3) evaluate the salient features and limitations of the solution. In the following sections the general solutions are first validated analytically against the solutions presented by van Genuchten [45] and Higashi and Pigford [24] for simplified transport scenarios. To further test the general solutions for a more general case, example simulations are performed and the results of the general solution are compare against the results generated from the semi-analytical solution given by Quezada et al. [38]. 34 3.2 Implementation of the Solution To develop a computer code for the general solution, one has to formulate an algorithm based on the looping structure presented in equation(43). In order to evaluate the summation loops, one has to compute all the ?G ? terms in equation(43), which are defined in equations (44) and (49). From these, equations one can observe that the ?G ? terms require the evaluation of ? 1 2 3, ,i i i F ? terms. These ? 1 2 3, ,i i i F ? terms, which are defined in equations (45), (46), (50) and (51), involve computations of the product of complementary error functions and exponential functions. After evaluating the ? 1 2 3, ,i i i F ? terms, one can then compute the ?G ? terms. The general solution can then be implemented by substituting the values for the corresponding ?G ? terms in the summation loop described by equation(43). Although the above procedure appears to be a straightforward task, round-off errors and underflow/overflow errors that occur during the computations of the ?G ? and ? 1 2 3, ,i i i F ? terms, can cause severe stability problems even for short chain lengths involving four species [45]. These computational errors are further aggravated when solving problems involving the Cauchy boundary condition for long chain lengths and/or large simulation times. As pointed out by van Genuchten [45] the round-off errors arise due to the approximations in the computation of the complementary error functions (ERFC) involved in the ? 1 2 3, ,i i i F ? terms, especially when employing the Cauchy boundary condition, for large retardation factors and long simulation distances or simulation times. Typically ERFC is evaluated by approximating it either as a closed-form analytical 35 expression [20] or as an infinite summation series [19]. Van Genuchten [45] suggested that this round-off error could be considerably limited by first substituting the approximate expressions for the ERFCs into the ? 1 2 3, ,i i i F ? terms and then combining and simplifying the resulting terms. The mathematical details of this strategy can be obtained from the computer code CHAIN [45]. Although this strategy significantly improves the round-off errors, it involves extensive book keeping procedures. It must be noted that van Genuchten [45] used the Gautschi [20] and Gautschi [19] approximations for computing the ERFC terms. Using a more accurate approximation for the ERFC function, such as the Cody [14] Chebyshev approximation, would provide additional improvement to the round-off errors. Another important computational challenge encountered in implementing this solution is the handling of very large and/or very small numbers involving the computations of the exponential and ?? ? functions. These numbers arise when solving problems involving long chain lengths and/or large simulation times. For example, consider the computations involved in the ? 21G ? term in equation(44). 36 [ ] ( ) [ ] ( ) ( ) ( ) ( ) ( ) ( ) ,, 1 2 31 1 2 2 1 2 2 2 3 3 1 3 2 3 2 2 3 1 2 2 2 4 2 3 2 4 2 4 2 3 4 1 4 2 4 1 4 2 4 34 2 4 2 , , , ,2 1 , , , , , , , , , , , , , , , , i i ii i i i i i i i i x a tx a t i i i i i i i i ii i i i R R i i i i i i i i i i i i i i i i i i i R R i i i i i i R R F x t e F x t e G a a R k R R a a ?? ? ? ?? ? ? = ? ? ? = ? = = ? ? ? ? ? + = ? ? ? ?? ? ? ?? ? ? ? ? ? ? (53) From equation(53) one can observe that the denominator term involves two independent product loops that have to be evaluated. Under certain parameter combinations, each of the arguments within these product loops can assume very small values (e.g., less than 1E-5) and when this loop is run for long chain lengths (e.g., greater than 10) the product of these arguments results in an extremely small number (lesser than 1E-50). Performing division by a small number of this order, results in very large arguments which cannot be represented accurately by using finite computer precision. Furthermore, mathematical operations (additions and/or subtractions) of these large arguments can result in severe round-off errors. Note that this problem is not only restricted to the ? 21G ? term, but can also apply to the other ?G ? terms. Similarly, one can also observe that the numerator of the ? 21G ? terms in equation(53) involves the computation of exponential functions. While performing simulations for large times and small distances, one may encounter very large arguments within these exponential functions that cannot be represented by finite computer precision. Therefore, it is common to observe overflow errors when solving problems involving large simulation times. 37 A powerful solution to tackle these two problems involving long chain lengths and large simulation times is to use a log based formulation. The log based formulation involves transforming the arguments whose products have to be evaluated in to the log space. Within the log space the products are evaluated by performing a set of additions. Finally, the evaluated product is inverse transformed from the log scale to the linear scale. By employing the proposed log based formulation, the underflow and/or overflow errors related to evaluating the ?? ? terms can be virtually eliminated. However, the proposed log based formulation can only be applied to compute product terms and cannot be used to compute summation terms. To tackle the problem of overflow errors while evaluating summation terms, one can make use of the symmetrical property of the ? 1 2 3, ,i i i F ? terms. This involves combining the ? 1 2 3, ,i i i F ? terms that have identical ? 2i ? and ? 3i ? values and then evaluating the sum of these combination terms separately and then substituting them in the main solution [45]. Although this method seems to have a direct approach, the algebraic manipulation involved in combining the various ? 1 2 3, ,i i i F ? terms requires the reformulation of the final solution in a different format. This reformulation would be mathematically tedious involving extensive book keeping, especially for problems having long chain lengths. One can use the log formulation and other computational techniques suggested above for developing a computer code for the general solution. In addition to this, in the section below, another alternative is presented where solutions for several simpler and more practical cases that are computationally less challenging to implement are deduced. These special-case solutions are particularly useful for developing screening level 38 models. Furthermore, some of these simple solutions can be directly compared against analytical solutions presented in the literature to provide powerful arguments for demonstrating the validity of the general solution. 3.3 Special Cases In this section solutions to several simplified transport scenarios are presented and where ever possible, these solutions are validated against previously published analytical solutions. It should be noted that for some of these special cases, a more rigorous mathematical analysis (rather than simple substitution) is required to derive these simplified solutions. Furthermore, unlike the general solution these special case solutions impose fewer restrictions on the transport parameter values (Appendices A and D) 3.3.1 Zero Initial Condition The general solution can be readily simplified to solve transport scenarios where the initial concentrations of all the contaminants are zero. This is done by substituting the value of ? oc ? as zero. For this case, the general solution simplifies to: ( ) ( ){ } 1 2 2 1 2 1 32 1 1 1 1 1 1 1 2 1 11 , 1, 2,... iii i i i i i i i ii i c x t y k G h G G i n ? = = == + ? ?? ?= + ? ?? ?? ? ? ?? ? ; ? = ? ??? (54) where; ? 1 2,i i a ? is given by: 1 2 1 21 2 1 , 2 ,, 2 ; 0 ; 0 i i i ii i i k when i Ra when i? ? > ?= ? ? =? (55) 39 Note that the ? 11G ? and ? 12G ? terms in equation(54) defined in Appendices A and D remain unchanged. It can be shown that the expressions for the first four species of this special case solution match the solutions given by van Genuchten [45]. Furthermore, from equation(54) one can observe that, for the zero initial condition case, the second term of the general solution given by equation(43) is absent. This not only relaxes some of the restrictions on the transport parameter values but also directly helps in improving the round off errors especially during large simulation times. 3.3.2 Identical Retardation Factors Special cases arise when the retardation factors of all the species in the transport problem are identical (for example, the transport of non-sorbing set of sequentially decaying contaminants, where the retardation factors of all the species will be 1). One can obtain the solution for this problem from the general solution by simple substitution. The modified solution for this special case problem for both types of boundary conditions is given as: ( ) { } { } 1 2 2 1 2 1 32 1 1 1 2 2 1 2 12 1 1 1 2 1 11 2 1 2 1 1 , 1, 2, iii i i i i i i i ii i ii i o i i i i i i ii i c x t y k G R c y k G i ? = = == + ? = == + ? ?? ?= ? ?? ?? ? ? ?? ? ? ?? ? + ? ?? ?? ? ? ?? ? ; ? = ? ??? ? ?? ...n (56) Note that the ? 12G ? and ? 22G ? terms defined in Appendices A and D remain unchanged. The exclusion of the ? 11G ? and ? 21G ? terms helps in obtaining a 40 computationally more stable solution for longer chain lengths. Additionally this solution imposes lesser restriction on the transport parameter values. Sun et al. [42] solved a transport problem similar to the above problem assuming zero initial conditions. However, unlike the solution presented in equation(56), Sun et al. [42] did not provide an explicit closed-form expression, instead they only present a computational algorithm to compute the concentration profiles. 3.3.3 Zero Advection Velocity In some situations the transport of contaminants is governed mostly by dispersion. The solution for this condition is identical to the general solution given by equation(43), except that the ? 1 2 3, ,i i i F ? terms will be modified. For the Dirichlet boundary condition, the ? 1 2 3, ,i i i F ? term is given as: [ ] , ,1 2 3 1 1 2 3 ,2 3 1 1 2 3 , ,1 2 3 1 1 2 3 1 1 1 2 3 1 2 3 1 , ,2 , , , ,2 , , , 2 , 2 2 ; 4 i i i x i i i i i x x i i i iD a t i x i i i x i i i iD i x i i i i i x i i i R x te erfc R D te F x t R x te erfc R D t kwhere R D a R ? ? ? ? ? ? ? ? ?? ?? ? ? ? ?? ? ? ?? ?? ? ? ?= ? ?? ?? ? ? ? ?+ ? ? ? ?? ?? ?? ? ? ? = ? ? ?? ? ? ? (57) The above expression is valid only for real values of ? 1 2 3, ,i i i ? ?. If ? 1 2 3, ,i i i ? ? is complex the ? 1 2 3, ,i i i F ? term is modified as: 41 [ ] ( ) , 1 2 3 1 2 32 3 1 2 3 1 1 2 31 1 2 3 1 2 3 1 1 * * , , , , , , * , ,* , , , , 2 2 ; 4 2 i ia t i i i i i i i i i x x i i i ii i i i i x i i i i x x xF x t e ACos BSin D D R x i tkwhere R D a and A iB erfc R R D t ? ? ?? ? ? ?? ? ? ?= ?? ?? ? ? ? ? ? ? ?? ?? ? ? ?? ? ? ?? ? + ? ? = ? ; + =? ? ? ?? ? ? ?? ? ? ? (58) For the case of the Cauchy boundary condition forcing the advection velocity to be zero would result in a zero boundary condition. To avoid this, it is assumed that the boundary condition at the source does not involve the ?v ? term and modify Cauchy source boundary condition as: ( ) 1 1 1 1 , 00, 1, 2,... 0, i i t i i oi ix o B e t tc tD i n x t t ?? = ? < ?? ? ? = ; ? = ?? ? >? ? (59) Under this condition, the ? 1 2 3, ,i i i F ? term is modified as: [ ] , ,1 2 3 1 1 2 3 1 2 3 1, 1 2 3 1 2 3 2 3, ,1 2 3 1 1 1 2 3 1 2 3 1 2 , , , , , , , 2 , , , , 2 , , 2 i i i x i i i i i x x D i i i i i i i i xa t i i i i i ix iD i i i i i i i i x R x te erfc R D t k F x t e when a andR R x te erfc R D t ? ? ? ? ? ? ? ? ? ?? ?? ? ?? ?+? ? ? ?? ? ? ?? ?= ; ? ? ? ? ?? ?+ ? ? ? ?? ? ?? ?? ?? ?? ? ( )211 1 11 1 2 3 1 11 1 1 2 3 1 2 3 1 4 , , , , 2 2 ; 4 ii i i x R x vtk t R R D t i i i i i x x ii x i i i i i x i i i R x kt xe e erfc when a R D D RR D t kwhere R D a R ? ? ? ?? ? ?? ? ? ?? ? = ? ; = ? ?? ?? ? ? ?? ? ? ? = ? ? ?? ? ? ? (60) 42 The above expression is valid only for real values of ? 1 2 3, ,i i i ? ?. If ? 1 2 3, ,i i i ? ? is complex the ? 1 2 3, ,i i i F ? term is modified as: [ ] ( ) 1 2 3 1 2 3 ,2 3 1 2 3 1 2 3 1 2 3 1 2 3 11 1 2 3 1 2 3 1 * , ,* , , 2 , , * 2 * , , , ,* , , * , , , 22 , 2 4 i i x i i i i i ixv xa t D i i i i i i i i i i i i x ii i i i i x i i i xB Cos D F x t e e xA Sin D R x ikwhere R D a and A iB erfc R ?? ? ? ? ? ? ? ?? ?? ? ?? ?? ? ? ?? ?= ? ?? ? ? ?? ? ? ? ?? ?? ?? ? ? ? + ; = ? + = ? ?? ? ? ? 1 2 3 1 * , , 2 i i i i x t R D t ?? ? ? ? ? ? ? ?? ? (61) The solutions obtained for the zero advection velocity case are specialized solutions that can be directly applied to model reactive transport scenarios involving chemical and nuclear repository leakage through diffusion. 3.3.4 Steady State The behavior of a contaminant plume under steady state conditions is of special interest especially in analyzing monitored natural attenuation (MNA) problems. Steady state solutions avoid the problems related to overflow errors which occur when employing the general solution to solve for large simulation times. Although this problem can be tackled by using the algebraic manipulation technique suggested in Section 3.2, it can be a tedious effort. Furthermore, numerical codes (e.g., RT3D [12]) and semi-analytical solutions (e.g., Quezada et al. [38]), require large amounts of computational resources to accurately model steady state solutions. For these reasons, it is very attractive to obtain explicit solutions for steady state conditions. It must be noted 43 that unlike previous special cases, it is not possible to deduce the steady state solution by direct substitution. Therefore, the steady state solution is derived by solving the steady state governing equations assuming a constant source for the Dirichlet and Cauchy boundaries conditions. Since the time term is absent in the governing equations, Laplace transform techniques are not involved, and one can directly use the linear transform method given by Clement [11] to uncouple the system of equations. The detailed solution procedure is given in Appendix F. The steady state solution for the Dirichlet boundary is given as: ( ) { } ( ) ( ) 2 2 1 2 2 2 1 1 2 2 12 1 2 3 3 1 3 2 42 1 1 11 , 1, 2,... x i x x v v D k Diii i i i i i i i i i i ii i i i i i i i ec x y k B k k i n ? ?? + ? ?? ? ? = = == + = ? ? ? ? ?? ?? ? ? ?= ? ?? ? ? ?? ?? ? ? ? ? ?? ? ; ? = ? ? ?? ? (62) The corresponding steady state solution for the Cauchy boundary is given as: ( ) { } ( ) ( ) 2 2 1 2 2 2 1 1 2 2 12 1 2 2 3 3 1 3 2 42 1 21 11 , 1 4 2 2 x i x x v v D k Diii i i i i i i i i i i ii i x i i i i i i i vec x y k B v v D k k k i ? ?? + ? ?? ? ? = = == + = ? ? ? ? ?? ?? ? ? ?= ? ?? ? ? ?? ?? ?? ? + + ? ?? ? ? ?? ?? ? ; ? ? ? ?? ? 1, 2,...n= (63) From equations(62) and (63) it can be observed that the steady state solutions are less complex compared to the generic transient solutions since they do not involve the ERFC terms. 44 3.3.5 Zero Dispersion Coefficient If we ignore the effects of dispersion, the governing transport simplifies to: ( ) ( ) ( ) ( ) ( ) 1 1 1 , , , , ; 2, 3,... , 1 i i i i i i i i i i c x t c x tR v y k c x t k c x t i n t x k c x t ; i ? ? ? ? ?+ = ? ? = ? ? = ? = 0 0t and x ; ? > < < ? (64) Note that in the absence of dispersion, the Cauchy boundary condition will be identical to the Dirichlet boundary condition. In Appendix G the above problem is solved using the generic exponential initial condition given by equation(2). The solution for the no dispersion condition can also be represented by equation(43). However, the ?G ? terms associated with this solution are defined as: (See equations (G.29) and (G.38) in Appendix G) 45 [ ] ( ) [ ] ( ) ( ) ( ) ( ) ( ) ( ) 3 2 3 1 2 33 1 3 2 2 4 2 2 4 4 1 4 2 2 4 3 2 5 2 4 2 5 2 5 2 4 5 1 5 2 5 1 5 2 5 45 2 5 2 , ,0 , , , , , ,1 1 , , , , , , , , , , , , , , , , , o i o i i i i i t i i i i i oi i t i i i i i i o i ii i i i R i i i i i i i i i i i i i i i i i R R i i i i i i R R F x t e F x t t B F x t e F x t tG a k R R a a ? ? ? ? ? = ? = ? = = ? ? ? ? ?? ?? ? ? ?? + ?? ?= ? ? ? ?? ? ? ? ?? ? ? ? ? ?( ) [ ] ( ) ( ) [ ] ( ) [ ] ( ) ( ) ( ) 4 2 3 3 1 2 3 1 2 3 2 4 4 1 4 2 4 2 ,, 1 2 31 1 2 2 1 2 2 2 3 2 3 1 2 2 2 4 2 4 1 4 2 4 2 , ,0 , ,1 2 , , , , , , ,2 1 , , , , , , , , , i i o i i i i i ii i i i i i R ti i i i i i i o i i i i i i i R R x a tx a t i i i i i i i i i i i i i i i i i i i R R B F x t e F x t t G k F x t e F x t e G a a R k R ? ?? ? ? ? = ? = ? ?? ? ? ? = ? = ? ?? ?? ? = ? ? ? + = ? ? ? ?? ?? ? ? ? ? ? ? ( )( )( ) [ ] ( ) ( ) 3 1 3 2 3 2 3 2 4 2 4 2 3 4 1 4 2 4 3 4 2 ,1 1 2 2 1 2 2 2 3 3 1 3 2 3 2 , , , , , , , , , , ,2 2 , , , , i i i i i i i i i i ii i i i R R i i i i i i i i i i i i i R R x a t i i i i i i i i i i i R R R a a F x t e G R k ? ? = ? ? = ? ? ? ? ? ? = ? = ? ? ? ? = ? ? ? ? (65) 46 where; the term ? 1 2 3, ,i i i F ? is given by: (See equation (H.5) in Appendix H) [ ] 1, 2 31 1 2 3, , , i i i R xa t vi i i i R xF x t u t e v ? ?? ?? ? ? ?? ?= ?? ? ? ? (66) It can be shown that in the absence of initial contaminant distributions the solutions derived in this section are identical to those published in Harada et al. [23] and Higashi and Pigford [24]. Note that, although the solution for the no dispersion case appears to be identical in structure to the general solution, the ? 1 2 3, ,i i i F ? terms involved in these ?G ? terms, defined in equation (H.5) in equation(66), are distinctly different from that of the general solution. These new ? 1 2 3, ,i i i F ? terms do not involve the ERFC terms and hence are easier to compute. 3.4 Example Problems A general computer code (FORTRAN) that solves equation(43) was developed based on the computational algorithm provided in Section 3.2. Appendix I gives the details of the code along with a sample input and output files. Using this code, the general solution is first tested against two published examples problems: one involving the Dirichlet boundary and the other involving the Cauchy boundary. The example problem for the Dirichlet boundary condition is taken from Cho [9]. The parameters used in this problem are summarized in Table 1. Note: it is assumed that the decay occurs in the liquid phase only. The results of the general solution are tested against the results of the analytical solution given by Cho [9]. Fig.1 compares the concentration profiles of the two solutions. It can be observed that the two solutions give identical results. 47 Parameter Value Velocity; v (cm/h) 1 Dispersion coefficient; Dx (cm2/h) 0.18 Simulation Time; t (h) 200 Pulse Time; to (h) 200 Retardation factors; R, 1 to 3 2, 1, 1 First-order decay coefficients k, 1 to 3 (1/h) 0.01, 0.1, 0 Yield coefficients; y, 1 to 2 1,1 First-order source decay coefficients; ?, 1 to 3 (1/h) 0,0,0 Boundary constant value species 1; B, 1 to 1 1 Boundary constant value species 2; B, 1 to 2 0, 0 Boundary constant value species 3; B, 1 to 3 0, 0, 0 Initial condition constant; co, 1 to 3 0, 0, 0 Initial condition Exponent; ?, 1 to 3 0, 0, 0 Table 1: Parameters used in the Cho [9] example problem 48 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 120 140 160 180 200 Distance (cm) Re lati ve C on ce ntr ati on Species1_Cho Species2_Cho Species3_Cho Species1 Species2 Species3 Fig. 1: Comparison of the concentration profiles of the general solution with the Cho [9] solution (3 Species, Dirichlet boundary condition) 49 The Cauchy boundary example problem is taken from van Genuchten [45]. This problem is a modified version of the problem given by Higashi and Pigford [24]. Here the source is modeled as a step source and decay is assumed to occur in both phases. The transport parameters used in this problem are summarized in Table 2. The results of the general solution are compared against the van Genuchten [45] solution and are summarized in Fig. 2. Observation of the concentration profiles given in Fig. 2 indicates that there is an excellent match between the two solutions. To further test the performance of the general solution for long chain lengths, a hypothetical 10 species problem was formulated that included an exponentially decaying boundary condition along with spatially varying initial contaminant distributions. Decay is assumed to occur only in the liquid phase and furthermore, the retardation factors of some of the species were made equal. Table 3 summarizes the model parameters used in this problem. The above problem was solved using a Dirichlet boundary condition assuming a pulse time for the source release ? ot ? as 10 years. The results of the general solution were compared against the semi-analytical solution given by Quezada et al. [38]. Fig. 3 shows the concentration distributions of all the species. These results indicate that the two solutions match well. An identical problem with a source pulse time ? ot ? of 20 yrs was formulated using a Cauchy boundary. The concentration distributions of all the species were compared against the Quezada et al. [38] solution. Fig. 4 gives the comparative profiles of both these solutions. Again, it can be seen that the profiles predicted by the two solutions have a good match. These simulations further validate the performance of the general solution for transport problems involving long chain lengths. 50 Parameter Value Velocity; v (m/yr) 100 Dispersion coefficient; Dx (m2/yr) 10 Simulation Time; t (yr) 1.0E+5 Pulse Time; to (yr) 1.0E+4 Retardation factors; R, 1 to 4 1.0E+4, 1.4E+4, 5.0E+4, 5.0E+2 First-order decay coefficients k, 1 to 4 (1/yr) 7.9E-3, 2.8E-6, 8.71E-6, 4.3E-4 Yield coefficients; y, 1 to 3 1,1,1 First-order source decay coefficients; ?, 1 to 4 (1/yr) 8.9E-3, 1.0028E-3, 1.0087E-3, 1.43E-3 Boundary constant value species 1; B, 1 to 1 1.25 Boundary constant value species 2; B, 1 to 2 -1.2504, 1.2504 Boundary constant value species 3; B, 1 to 3 4.43684E-4, 5.93431E-1, -5.93874E-1 Boundary constant value species 4; B, 1 to 4 -5.1674E-6, 1.20853E-2, -1.22637E-2, 1.78958E-4 Initial condition constant; co, 1 to 4 0, 0, 0, 0 Initial condition Exponent; ?, 1 to 4 0, 0, 0, 0 Table 2: Parameters used in the modified Higashi and Pigford [24] example problem as given by van Genuchten [45] 51 1.E-09 1.E-08 1.E-07 1.E-06 1.E-05 1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1 10 100 1000 10000 100000 Distance (m) Re lat ive C on ce nt ra tio n Species2_VG Species3_VG Species4_VG Species2 Species3 Species4 Fig. 2: Comparison of the concentration profiles of the general solution with the van Genuchten [45] solution (4 Species, Cauchy boundary condition) 52 Parameter (Should be updated) Value Velocity; v (m/yr) 5 Dispersion coefficient; Dx (m2/yr) 50 Simulation Time; t (yr) 20 Pulse Time; to (yr) 10 (For Dirichlet) and 20 (For Cauchy) Retardation factors; R, 1 to 10 1.9, 1, 1.4, 1, 5, 8, 1.4, 3.1, 1,1 First-order decay coefficients k, 1 to 10 (1/yr) 3, 2, 1.5, 1.25, 2.75, 1, 0.75, 0.5, 0.25, 0.1 Yield coefficients; y, 1 to 9 1, 2, 1.5, 0.4, 1, 1, 0.7, 0.9, 1 First-order source decay coefficients; ?, 1 to 10 (1/yr) 0.1, 0.75, 0.5, 0.25, 0, 0, 0.3, 1, 0, 0.65 Boundary constant value species 1; B, 1 to 1 10 Boundary constant value species 2; B, 1 to 2 0, 5 Boundary constant value species 3; B, 1 to 3 0, 0, 2.5 Boundary constant value species 4; B, 1 to 4 0, 0, 0, 0 Boundary constant value species 5; B, 1 to 5 0, 0, 0, 0, 10 Boundary constant value species 6; B, 1 to 6 0, 0, 0, 0, 0, 5 Boundary constant value species 7; B, 1 to 7 0, 0, 0, 0, 0, 0, 2.5 Boundary constant value species 8; B, 1 to 8 0, 0, 0, 0, 0, 0, 0, 0 Boundary constant value species 9; B, 1 to 9 0, 0, 0, 0, 0, 0, 0, 0, 0 Boundary constant value species 10; B, 1 to 10 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 Initial condition constant; co, 1 to 10 0, 0.1, 0.2, 0, 0.25, 0.3, 0.15, 0, 0, 0 Initial condition Exponent; ?, 1 to 10 0, 0.01, 0, 0, 0.02, 0.01, 0.1, 0, 0, 0 Table 3: Parameters used in the new 10 species example problem 53 1.E-08 1.E-07 1.E-06 1.E-05 1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 0 20 40 60 80 100 120 140 160 180 200 Distance (m) Co nc en tra tio n ( Mg /L) Species1_Q Species2_Q Species3_Q Species4_Q Species5_Q Species6_Q Species7_Q Species8_Q Species9_Q Species10_Q Species1 Species2 Species3 Species4 Species5 Species6 Species7 Species8 Species9 Species10 Fig. 3: Comparison of the concentration profiles of the general solution with the Quezada et al. [38] semi-analytical solution (10 Species, Dirichlet boundary condition) 54 1.E-08 1.E-07 1.E-06 1.E-05 1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 0 20 40 60 80 100 120 140 160 180 200 Distance (m) Co nc en tra tio n (M g/L ) Species1_Q Species2_Q Species3_Q Species4_Q Species5_Q Species6_Q Species7_Q Species8_Q Species9_Q Species10_Q Species1 Species2 Species3 Species4 Species5 Species6 Species7 Species8 Species9 Species10 Fig. 4: Comparison of the concentration profiles of the general solution with the Quezada et al. [38] semi-analytical solution (10 Species, Cauchy boundary condition) 55 3.5 Salient Features of the General Solution One of the most important features of the general solution is that it is a generic solution that can solve an arbitrary number of species. Furthermore, the final solution is given in an explicit closed-form format that avoids numerical methods. In addition to this, the general solution is formulated and presented in format such that solutions for the Dirichlet and the Cauchy boundaries share a common structure. This enables the convenient implementation of both these solutions by designing one code that shares several common routines having the capability to easily switch between the two boundary conditions. Another interesting advantage of this solution is that it can be used to directly obtain the solution for a particular species in the chain without involving the computations of other species. This is extremely efficient when modeling scenarios where there the daughter product would be more toxic or more mobile than the parent compound. In this case, one would be interested in obtaining the concentration profiles for the daughter product only. Solving this problem numerically would involve the computation of all the predecessor species and this can be time consuming. A similar situation also arises when employing the recursive type semi-analytical solutions presented by Harada et al. [23] (See section4, p6, eq4.20). Their solution expression for any species ?i>1' involves the concentration term of its predecessor ?i-1?th species. The third feature of this solution is that it can be readily extended to a general diverging reaction network [43]. One can conceptualize a diverging network as a superposition of a series of parallel reaction networks each of which can be solved independently. It must be noted that one diverging species will have no effect on the 56 other diverging species and hence they can be considered independent of each other. Hence in order to solve a diverging reaction network, one has to split the network into a set of parallel sequential reaction networks and solve each sequential chain independently. The one dimensional solutions presented in this work can be readily extended to solve multi-dimensional transport problems involving transverse dispersion terms using the approximate Domenico solution [16]. The mathematical details of this extension strategy are provided in Quezada [37]. As pointed out by Quezada [37] , this strategy yields an approximate solution. In the following chapter this approximation is investigated to identify transport situations where the Domenico solution can be advantageously applied with high accuracy [41]. 3.6 Limitations of the General Solution Despite its advantages, the general solution and its simplified solutions have some key limitations. As with the case of all analytical solutions, these solutions can model idealized transport problems only. Situations involving heterogeneity, variation in the flow field, presence of source/sink terms, etc. cannot be modeled using these solutions. Furthermore, these solutions are restricted to modeling first-order kinetic reactions only. Although this is highly suitable for modeling radioactive waste transport and certain types of simplified chemical reactions, these solutions cannot capture the more generic reaction kinetics such as the Monod kinetics. The general solution places several limitations on the values of the transport parameters which are summarized in Table 4 (Appendices A and D). These limiting 57 conditions arise when factorizing the denominator of the ?G ? terms during the computation of the inverse Laplace transforms. An easy approach to tackle this problem would be to marginally perturb the individual model parameters such that the limiting conditions are averted [45]. This would enable us to tackle the problem using the general code. However, this technique can aggravate the existing computational problems related to overflow/under flow errors. Another more rigorous approach to tackle this problem would be to perform a mathematical analysis incorporating these special cases and evaluating the corresponding alternate expressions, similar to the one performed on the case of identical retardation factors (i.e. when ? 1 2, 1 2 1 2 0 ; , 1,2,...i iR i i n where i i= ? = ; ; ? ?). This analysis would involve modifying the factorization of the ?G ? terms to include the limiting conditions and performing generalized inverse Laplace transform operations on the alternate expressions. Finally, it must be noted that unlike the semi-analytical solutions presented by Clement [11] or Quezada et al. [38] which can model arbitrary reaction networks, the general solution can only model sequential and diverging reaction networks. However, it is possible to extend this solution to solve systems involving complex reaction networks using a similar approach. The key step required to perform this extension would be to obtain generalized expressions for the linear transform matrices used to uncouple the coupled system of equations. 58 1 2, 1 2 1 2 0 ; , 1,2,...i ik i i n where i i ? ? = ; ? 1 2 1 3 1 2 1 3 , , 1 2 3 1 2 1 3 2 3 ; , , 1,2,...i i i i i i i i a a i i i n where i i i i i i R R and R R ? ? = ; ? , ? , ? , ? ? 1 2 3 1 2 , 1 2 3 1 2 ; , ; 1,2,...i i i i i a i i i n i i i n where i i and R R ?? ? ? 1? ? ; ? = ; ? ? 1 2 1 1 2 , , 1 2 1 2 ; , ; 1,2,...i i i i i i a a i i i n i n where i i and R R ?? ? ? ? = ; ? ? Table 4: Summary of parameter limitations 59 CHAPTER IV INVESTIGATION OF ERROR ASSOCIATED WITH THE DOMENICO APPROXIMATION 4.1 Introduction This chapter investigates some of the approximations involved in extending the one-dimensional solutions derived earlier into multiple-dimensions. One of the most popular multi-dimensional analytical solutions used for modeling multi-dimensional ground water contaminant plumes is the Domenico solution [16]. The Domenico solution is an approximate three-dimensional solution that describes the fate and transport of a decaying contaminant plume evolving from a finite patch source. This solution was based on an approach previously published by Domenico and Robbins [17] for modeling a non-decaying contaminant plume. Prior to this work, Cleary and Ungs [10] presented an analytical solution to a similar three dimensional transport problem for a domain finite in y and z directions. Later, Sagar [39] published an exact analytical solution to the transport problem considered by Domenico and Robbins [17]. Wexler [48] extended the Sagar [39] solution to include the effects of reaction and presented an exact analytical solution to the transport problem considered by Domenico [16]. However, these solutions are not closed form expressions since they involve numerical evaluation of a definite integral. This numerical integration step can be computationally demanding and can also introduce numerical errors. The key advantage of the Domenico and Robbins 60 [17] approach is that it provides a closed form solution without involving numerical integration procedures. Due to this computational advantage, the Domenico solution has been widely used in several public domain design tools including the USEPA tools BIOCHLOR and BIOSCREEN [3, 36]. Although the Domenico solution is extensively employed in several ground water transport models, its approximate nature has received mixed reviews over the years. For example, West and Kueper [47] compared the BIOCHLOR model against a more rigorous analytical solution and observed considerable discrepancies. By comparing the near field concentration profiles they concluded that the Domenico solution can produce errors up to 50%. Guyonnet and Neville [22]compared the Domenico solution against the Sagar [39] solution and presented the results in a non-dimensional form. They concluded that for ground water flow regimes dominated by advection and mechanical dispersion the discrepancies between the two solutions can be considered negligible along the plume centerline. They further added that the errors may increase significantly outside the plume centerline. The above review indicates that there are conflicting opinions regarding the performance of the Domenico solution. Furthermore, since the development of the Domenico solution was based on a heuristic approach, researchers have expressed skepticism regarding its validity [47]. Presently, there are several unanswered issues related to the performance of this solution that include: Is there a mathematical basis for deriving the Domenico solution? If so, what are the approximations involved in deriving the solution? What are the errors associated with these approximations? And finally, under what conditions are these approximations valid? To answer these questions, we 61 need a fundamental understanding of the nature of the approximations involved in the Domenico solution. The focus of this chapter is to perform a rigorous mathematical analysis on the origin and development of the Domenico solution. The outcomes of this analysis are used to develop some general guidelines for the appropriate use of the solution. 4.2 Governing Equations The transport problem considered by Domenico [16] assumes a patch source of constant concentration ? oc ? located at the origin in a clean, semi-infinite aquifer. The contaminant is subjected to advection in the x direction and dispersive mixing in all three directions. Furthermore, it is assumed that the contaminant decays through a first order process. The governing transport equation considered by Domenico [16] is: 2 2 2 2 2 2 ? ? ? ? ?= ? + + + ? ? ? ? ? ?x y z c c c c cv D D D kc t x x y z (67) The initial and boundary conditions are: 0 ( , , ,0) 0 0 , , (0, , , ) , , 02 2 2 2 0 , 0 ( , , , )lim 0 ( , , , )lim 0 ( , , , )lim 0 x y z c x y z x y z Z Z Y Yc y z t c z y t otherwise t c x y z t x c x y z t y c x y z t z ?? ??? ??? = ? < < ? ?? < < ? ?? < < ? = ? < < ? < < ? > = ? > ? = ? ? = ? ? = ? (68) 62 where; ?c ? is the concentration of the contaminant [ML-3]; ? oc ? is the concentration at the source [ML-3]; ?Y ? and ?Z ? are the source dimensions in y- and z- directions, respectively [L]; ? xD ?, ? yD ?, and ? zD ? are the dispersion coefficients in x, y and z directions, respectively [L2T-1]; ?v ? is the advection velocity in the x direction [LT- 1]; and ?k ? is the first-order decay coefficient [T-1]. 4.3 Review of the Domenico Solution The Domenico [16] solution was based on an approximate approach given by Domenico and Robbins [17]. Therefore, we first present a detailed review of the development of the Domenico and Robbins solution. Domenico and Robbins [17] began their analysis by presenting the following exact analytical solution that describes the transport of an instantaneous pulse source in a three-dimensional domain [25]: ( ) ( ) ( ) ( ) 1/ 2 1/ 2 1/ 2 1/ 2 ( , , , ) ( , ) ( , ) ( ,8 2 2; ( , ) 2 2 2 2( , ) 2 2 o x y z x x x y y y cc x y z t f x t f y t f z t X Xx vt x vt where f x t erf erf D t D t Y Yy y f y t erf erf D t D t = ? ?? ? ? ?? + ? ? ? ?? ? ? ? = ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ?+ ? ? ?? ? ? ?= ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ? ? ?? ? ? ? ? ? ? ( ) ( )1/ 2 1/ 2 2 2( , ) 2 2z z z Z Zz z and f z t erf erf D t D t ? ?? ? ? ?+ ? ? ?? ? ? ? = ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ? ? ?? ? ? (69) They then present the following one-dimensional analytical solution [15]: 63 ( )1/ 2 ( , ) ( , )2 ; ( , ) 2 o x x x cc x t f x t x vtwhere f x t erfc D t = ? ??? ? = ? ? ? ?? ? (70) Note that the above expression is the solution to the standard one-dimensional advection dispersion equation for an instantaneous source extending from zero to negative infinity [5]. To account for the transverse dispersion due to a finite sized two-dimensional source they employed the following two analytical solutions [15]: ( ) ( ) ( ) ( ) 1/ 2 1/ 2 1/ 2 1/ 2 ( , ) ( , ) ( , ) ( , )2 2 2 2; ( , ) 2 2 2 2( , ) 2 2 o o y z y y y z z z c cc y t f y t and c z t f z t Y Yy y where f y t erf erf D t D t Z Zz z and f z t erf erf D t D t = = ? ?? ? ? ?+ ? ? ?? ? ? ? = ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ?+ ? ? ?? ? ? ? = ? ? ?? ? ? ? ? ? ? ? ?? ? ? ? ?? ? ?? (71) Note that ( , )c y t and ( , )c z t are solutions to two independent one-dimensional transient diffusion equations subjected to an instantaneous line source of widths Y and Z respectively. Further, it can be observed that the terms ( , )yf y t and ( , )zf z t in equation(71) are identical to the terms ( , )yf y t? and ( , )zf z t? in the Hunt [25] solution. Domenico and Robbins [17] multiplied the one-dimensional solution ( , )xf x t with these ?transverse spreading terms? ( , )yf y t and ( , )zf z t and presented the following expression. 64 ( , , , ) ( , ) ( , ) ( , )8o x y zcc x y z t f x t f y t f z t= (72) However, the authors did not justify this superposition step. Note that the Hunt solution [25] was never used in this analysis. At this stage, Domenico and Robbins presented the following arguments: ?The product of these three integral solutions [shown in equation(72)] describes a semi-infinite contaminated parcel which moves in the positive x direction with a one-dimensional velocity but which continuously expands in size in directions transverse to x throughout the whole domain of x, i.e., in the positive and negative regions. This is because the time t in the transverse spreading terms is interpreted as running time. Reinterpreting this time as x/v for a moving coordinate system, as is common in all transverse spreading models (Bruch and Street, 1967; Ogata, 1970; Domenico and Palciauskas, 1982) has the effect of maintaining the original source dimensions at x=0 so that the condition C=Co is maintained at x=0 for t>0.? Using these arguments they reinterpret the time term t in the transverse spreading terms ( , )yf y t and ( , )zf z t as x/v. However, the authors did not provide a mathematical reasoning for this time reinterpretation step. Further, all the references cited in the above text solve fundamentally different problems and we will address this issue in a later section. Using this time reinterpretation step, equation(72) was modified as: 65 ( ) ( ) ( ) ( ) 1/ 2 1/ 2 1/ 2 1/ 2 ( , , , ) ( , ) ( , ) ( , )8 ; ( , ) 2 2 2( , ) 2 / 2 / 2( , ) 2 / o x y z x x y y y z z cc x y z t f x t f y x f z x x vtwhere f x t erfc D t Y Yy y f y x erf erf D x v D x v Zz and f z x erf erf D x v = ? ??? ? = ? ? ? ?? ? ? ?? ? ? ?+ ? ? ?? ? ? ?= ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ? ? ?? ? ? ?+ ? ? = ?? ? ? ?? ? ( )1/ 2 2 2 /z Zz D x v ? ?? ?? ? ?? ? ? ?? ? ? ?? ?? ? ? ?? ? (73) Equation(73) was presented as the final solution to the continuous finite patch source problem considered by Domenico and Robbins [17]. Later, Domenico [16] incorporated the effects of first-order decay by replacing the ( , )xf x t term with an analytical solution for the semi-infinite pulse source problem with a decay term, presented by Bear [5]. The final solution was given as [16]: ( ) ( ) 1/ 2 1/ 2 1/ 2 1/ 2 ( , , , ) ( , ) ( , ) ( , )8 41 4; ( , ) exp 1 1 2 2 2( , ) 2 o x y z x x x x x y y cc x y z t f x t f y x f z x kx vt kx vwhere f x t erfc v vt Yy y f y x erf erf x = ? ?? ??? ? ? ?? +? ?? ?? ?? ??? ? ? ?? ? ? ? ? ? = ? + ? ?? ? ? ?? ? ? ?? ? ? ?? ?? ? ? ?? ?? ? ? ?? ?? ? ? ?? ? ? ?+ ? ?= ?? ? ?? ? ? ? ( ) ( ) ( ) 1/ 2 1/ 2 1/ 2 2 2 2 2( , ) 2 2 y z z z Y x Z Zz z and f z x erf erf x x ? ?? ?? ? ?? ? ? ?? ?? ? ?? ?? ? ? ?? ? ? ?? ? ? ?+ ? ? ?? ? ? ? = ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ? ? ?? ? (74) 66 where; xx Dv? = , yy Dv? = and zz Dv? = are the dispersivities in the x, y and z directions respectively [L]. Martyn-Hayden and Robbins [31] later modified the Domenico [16] solution, referred in this work as modified-Domenico solution, by incorporating the following one- dimensional solution (which describes a constant source boundary) in the ( , )xf x t term [5]: ( ) ( ) 1/ 2 1/ 2 1/ 2 1/ 2 1/ 2 1/ 2 41 4exp 1 1 2 2 41 4exp 1 1 2 2 x x x x x x x x x kx vt kx vf erfc v vt kx vt kx verfc v vt ? ?? ??? ? ? ?? +? ?? ?? ?? ??? ? ? ?? ? ? ? ? ?= ? + ? ?? ? ? ?? ? ? ?? ? ? ?? ?? ? ? ?? ?? ? ? ?? ?? ? ? ?? ? ? ??? ?+ + ? ?? ?? ?? ??? ? ? ?? ? ? ? + + + ? ?? ? ? ?? ?? ? ? ?? ?? ? ? ?? ?? ? ? ?? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? (75) As pointed out by Bear [5] if the value of x x ? is sufficiently large, a condition usually satisfied in practice, the additional term in the above equation can be safely ignored. The above review shows that the development of various forms of the Domenico solution does not have a rigorous mathematical basis. The empirical arguments provided by the authors are vague because the mathematical procedures implied by these arguments are inexplicit and nebulous. In the following section a more rigorous approach to derive the Domenico solution that clearly stated the approximations involved is presented. 67 4.4 A Rigorous Approach to Derive the Domenico Solution The exact semi-analytical solution for the three dimensional transport problem described by equations (67) and (68), considered by Domenico [16], was provided by Wexler [48] as: ( ) ( ) 0 2 2 3/2 1/ 2 1/2 ( , , , ) ' , ) ' , ) ' , )8 exp 4 4 ; ' , ) exp 2 2 2' , ) 2 2 ' t o x y z x x x xx y y y z cc x y z t f x f y f z d v xk D Dx vxwhere f x DD Y Yy y f y erf erf D D and f ? ? ? ? ? ? ? ? ? ? ?? ? ? ? = = = ( ( ( ? ?? ?? + ? ?? ? ? ? ( = ? ? ? ? ? ?? ? ? ?+ ? ? ?? ? ? ?( = ? ? ?? ? ? ? ? ?? ? ? ?? ?? ? ? ? ? ? ( ? ( ) ( )1/ 2 1/2 2 2, ) 2 2z z Z Zz z z erf erf D D ? ? ? ? ?? ? ? ?+ ? ? ?? ? ? ?= ? ? ?? ? ? ? ? ?? ? ? ?? ?? ? ? ? ? ? (76) To obtain the Domenico solution from the above exact solution, we replace the value of ? in the transverse spreading terms ' , )yf y ?( and ' , )zf z ?( with x/v (the validity of this substitution will be discussed later). This yields the following expression: ( ) ( ) ( ) ( ) 0 1/2 1/ 2 1/ 2 1/ 2 ( , , , ) ( , ) ( , ) ' ( , )8 2 2; ( , ) 2 2 2 2( , ) 2 2 t o y z x y y y z z z cc x y z t f y x f z x f x d Y Yy y where f y x erf erf x x Z Zz z and f z x erf erf x x ? ? ? ? ? ? ? ? = = = ? ?? ? ? ?+ ? ? ?? ? ? ? = ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ?+ ? ? ?? ? ? ? = ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ?? ? ? ? ?? (77) 68 where; yy Dv? = and zz Dv? = . Note that by substituting ? = x/v we have made the transverse spreading terms ( , )yf y x and ( , )zf z x independent of time and hence they will not participate in the integration process. Without the transverse terms, the definite integral can be evaluated analytically as shown in Appendix J. Therefore, the above equation can be simplified as: ( ) 1/2 1/ 2 1/2 ( , , , ) ( , ) ( , ) ( , )8 41 4; ( , ) exp 1 1 2 2 exp 1 12 o x y z x x x x x x cc x y z t f x t f y x f z x kx vt kx vwhere f x t erfc v vt x ? ? ? ? ? = ? ?? ?? ? ? +? ?? ?? ?? ?? ?? ? ? ?? ? ? ? ? ? = ? + ? ?? ? ? ?? ? ? ?? ?? ?? ? ? ?? ?? ? ? ?? ?? ? ? ?? ? + + ( ) 1/ 2 1/2 1/ 2 41 4 2 x x x kx vt k verfc v vt ? ? ? ? ?? ?? ? + +? ?? ?? ?? ?? ?? ? ? ?? ? ? ? ? ?+ ? ?? ? ? ?? ? ? ?? ?? ?? ? ? ?? ?? ? ? ?? ?? ? ? ?? ? (78) Equation(78) is identical to the modified-Domenico solution shown in equation(75). If we set the first-order decay coefficient k to zero equation(78) reduces to: ( ) ( )1/2 1/2 ( , , , ) ( , ) ( , ) ( , )8 ; ( , ) exp 2 2 o x y z x xx x cc x y z t f x t f y x f z x x vt x x vtwhere f x t erfc erfc vt vt?? ? = ? ? ? ?? ?? +? ? ? ? = +? ? ? ? ? ? ? ?? ? ? ?? ? ? ? (79) Equation(79) is similar to the Domenico and Robbins (1985) solution given by equation(73). The additional expression in the ( , )xf x t term in equations (78) and (79) is due to the use of the expanded form of the one-dimensional solution that describes a constant concentration boundary condition instead of a semi-infinite pulse source boundary condition. 69 The above analysis shows that the only approximation required for rigorously deriving the Domenico solution is the time reinterpretation step, where ? is replaced by x/v in the transverse dispersion terms. In the following section, a detailed mathematical analysis that investigates the validity of this approximation is presented 4.5 Mathematical Analysis of the Validity of the Approximation Involved in the Domenico Solution Review of transport modeling literature indicates that it is common to replace ? with x/v in the transverse dispersion terms when solving convection dominated problems that have low longitudinal mixing. For example, Bruch and Street [6] used a similar assumption to solve the advection-dispersion problem when the longitudinal mixing was smaller than the transverse mixing. Another example of a convection-dominated problem that employs this approximation is the air pollution model used for predicting the fate and transport of smoke plumes evolving from chimneys [46]. Here, the transport is dominated by convection along the wind direction and dispersive mixing is restricted to the transverse directions only. Neglecting the effects of longitudinal dispersion in such problems simplifies the governing transport equation as: 2 2 2 2y z c c c cv D D kc t x y z ? ? ? ?= ? + + ? ? ? ? ? (80) Consider solution to the above transport problem subject to the following Domenico-type initial and boundary conditions: 70 0 ( , , ,0) 0 0 , , (0, , , ) , , 02 2 2 2 0 , 0 ( , , , )lim 0 ( , , , )lim 0 y z c x y z x y z Z Z Y Yc y z t c z y t otherwise t c x y z t y c x y z t z ??? ??? = ? < < ? ?? < < ? ?? < < ? = ? < < ? < < ? > = ? > ? = ? ? = ? (81) In Appendix K, Laplace transform techniques are used to solve the above problem and the resulting exact analytical solution is: ( , , , ) ( , ) ( , ) ( , )8 ; ( , ) 2exp ; 0 1 oo x y z o x cc x y z t f x t f x y f x z k x xwhere f x t u t v v xand u t isthe step function givenby v xif t x vu t xv if t v = ? ? ? ? = ? ?? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??? ?? = ? ? ?? ? ? > ?? (82) where; ( , )yf x y and ( , )zf x z are identical to the expressions given in equation(77). Since the Domenico approach approximates ? as x/v in the transverse dispersion terms one can hypothesize that the Domenico approximation must be valid when x? is zero. To test this hypothesis, a limiting analysis is performed on the modified-Domenico solution by forcing x? to zero; this is expressed as: 0 ( , , , ) lim , ) , ) ( , )8 x o x y z cc x y z t f x t f y x f z x ? ? = ( ( (83) 71 The mathematical details of this limiting analysis are given in Srinivasan et al. [41]. The analysis shows that when x? approaches zero the modified-Domenico solution relaxes to the exact analytical solution given by equation(82). This proves that the Domenico approximation indeed yields an exact analytical solution when x? is equal to zero. 4.6 Analysis of the Error Associated With the Domenico Solution The mathematical analysis presented in the previous section demonstrates that the time reinterpretation step, where ? is replaced with x/v, is exactly valid when x? = 0. From these results one could also infer that this time reinterpretation process provides a reasonable approximation when longitudinal dispersion plays an insignificant role in the overall transport. Hence, the Domenico solution can be expected to produce reasonable estimates for advection-dominated problems; however, it can introduce significant errors for longitudinal dispersion-dominated problems. Another important feature of the time reinterpretation step is that it forces a quasi steady-state condition along the transverse directions at all times. In other words, the ?conceptual? residence time (x/v value) associated with a point located at the centerline to disperse contaminant mass in the transverse directions is independent of the simulation time. Further, this residence time is also assumed to increase linearly with respect to x. These unrealistic assumptions regarding residence times will lead to erroneous predictions, especially beyond the advective front. For example, consider a problem where v = 50 m/year and we are interested in predicting the concentration distribution of 72 a two year plume (t = 2 years) at a location x = 200 meters. The Domenico solution will estimate the residence time ? for our location of interest x = 200 meters as ? = x/v = 4 years; this in fact is greater than the total simulation time it self! This is an unrealistic assumption since a 2 year old plume simply cannot have the time to disperse for 4 years! For a particle located at the advective front the residence time assumed by the Domenico solution is 2 years (the simulation time), and for all the particles located behind the advective front the residence time assumed by the Domenico solution will be equal to x/v (which will be less than 2 years); these seem to be reasonable estimates. However, for all points beyond the advective front i.e. x > 100 meters, the Domenico solution will assign unrealistic conceptual residence times which will be greater than the simulation time t = 2 years. It must be noted that this incorrect behavior will vanish when x? is zero because, for this case, the plume will abruptly end at the advective front and the residence time for each particle located at or behind the advective front will in fact be equal to x/v. When solving steady state problems, the assumption related to residence time should be a reasonable approximation. This is because, at steady state, the theoretical advective front will be at infinity. Therefore, the time reinterpretation should be reasonable for any finite domain. Hence, the performance of the Domenico solution under steady state conditions can be expected to be better. However, it is important to note that even under steady state conditions the solution will not be exact because it will still ignore the transport due to longitudinal mixing. In general, it can be concluded that the Domenico solution can be expected to perform better behind the advective front. In the following section we use an example problem to illustrate the implication of these theoretical results. 73 4.7 Example Problem The example problem presented by Domenico and Robbins [17] is considered in this analysis. The transport parameters used in the problem are summarized in Table 5. The performance of the modified-Domenico solution was tested by comparing its results against those generated using the exact solution given by Wexler [48]. It has been established in the previous sections that the Domenico approximation makes unreasonable assumptions regarding the residence time beyond the advective front, and reasonable assumptions behind the front. Therefore we analyze the results of this comparison in two parts-- one behind the advective front and the other beyond the advective front (note that for our base case, the front is at x = 1100 meters). 4.7.1 Plume Comparison Analysis Behind the Advective Front Figures 5a and 5b compare the two-dimensional concentration contours of both solutions on the X-Y and X-Z planes (Note: an aspect ratio of ?2.2 : 1? was maintained for the X-Y plots and an aspect ratio of ?55 : 1? was maintained for the X-Z plots). Since the problem is symmetric about the X axis, only half of the plume is presented. It can be observed from Figure 5 that the modified-Domenico solution is reasonably close to the true solution, though there are some noticeable discrepancies. To explore the limits of these discrepancies, we performed a series of sensitivity simulations. 74 Parameter Value Longitudinal dispersivity (?x) 42.58 m Transverse dispersivity (?y) 8.43 m Transverse dispersivity (?z) 0.00642 m Velocity (v) 0.2151 m/day Source width in Y direction (Yo) 240.0 m Source width in Z direction (Zo) 5.0 m Source concentration (Co) 850 mg/L Simulation Time (T) 5110 days Table5: Parameters used in the Domenico and Robbins [17] example problem 75 (a) (b) Fig. 5: Concentration contours predicted by the Domenico and Wexler solutions for the base case: solutions behind the advective front for (a) X-Y plane (b) X-Z plane 76 In the first set of sensitivity simulations, we varied the value of the longitudinal dispersivity (?x) by an order of magnitude above and below the assumed baseline value. Figures 6a and 6b compare the two-dimensional concentration contours of the solutions for both the cases. Comparison of the data shown in Figures 5a and 6 indicate that the discrepancies between the two solutions were large when the value of longitudinal dispersivity was high. Also, as expected, when the longitudinal dispersivity was low there was an excellent match between the solutions. Similar trends were also observed in the concentration contours predicted on the X-Z plane. Since the spreading terms in the y and z directions are identical in structure, the contours in the X-Y and X-Z planes will exhibit identical trends. Therefore, from this point onwards the sensitivity analysis will be restricted to X-Y contours. In the second set of sensitivity simulations, we varied the value of the transverse dispersivity (?y) by an order of magnitude above and below the base line value. These results (refer supplementary material for figures) indicated that the transverse dispersivity in the y-direction does not play a significant role in influencing the error associated with the modified-Domenico solution. Similar sensitivity analysis performed on other transport parameters including the transverse dispersivity ?z and the source dimensions Y and Z also showed minimal sensitivity. A third set of sensitivity simulations were completed for a decaying contaminant plume by assuming various first-order rate coefficients (k). Comparison of the concentration contours for k values of 0.0001 day-1 and 0.001 day-1 indicated that the presence of a decay term does not introduce any significant additional error [41]. 77 (a) (b) Fig. 6: Sensitivity results for variations in the longitudinal dispersivity value; solutions behind the advective front for (a) ?x*10 (b) ?x/10 78 The fourth set of sensitivity simulations involved varying the advection velocity (v) by an order of magnitude above and below the baseline value. Figures 7a and 7b compare the concentration contours of the solutions for the two cases. From these figures it can be concluded that the advection velocity has very little effect in determining the accuracy of the solution. Note that in the absence of first order decay, varying the advection velocity will have the same effect as varying the total simulation time (t). Since decay does not play any significant role in determining the accuracy of the modified-Domenico solution, it can be safely concluded that variations in the total simulation time will have little sensitivity on its accuracy. The above results indicate that within the advective front the longitudinal dispersivity plays a very important role in determining the accuracy of the modified- Domenico solution. All the other transport parameters have negligible effect on the accuracy of the solution. 4.7.2 Plume Comparison Analysis Beyond the Advective Front Figure 8 compares the concentration contours of the two solutions in the X-Y plane for the base case parameters (Note: here an aspect ratio of ?4:1? is maintained for the X:Y plane to capture the plume beyond the advective front; also the location of the advective front is indicated by an arrow on the x axis). It can be observed from Figure 8 that, as we move beyond the advective front, the accuracy of the modified-Domenico solution reduces rapidly. As pointed out in the earlier sections, this is due to the unrealistic assumptions made by the Domenico solution when computing the conceptual residence times beyond the advective front. 79 (a) (b) Fig. 7: Sensitivity results for variations in the transport velocity; solutions behind the advective front for (a) v*10 (b) v/10 80 Fig. 8: Concentration contours predicted by the Domenico and Wexler solutions; solutions include concentration contours beyond the advective front 81 The results of a sensitivity analysis performed on the parameter ?x are summarized in Figures 9a and 9b. These figures indicate a trend similar to those present for regions within the advective front. Higher the value of the longitudinal dispersivity, greater the error associated with the modified-Domenico solution. Further, it can be observed that the error systematically increases when the contaminant is transported beyond the advective front. A similar set of sensitivity simulations were performed on different transport parameters including ?y, ?z, Y, Z and K for regions beyond the advective front as well. As expected, the results indicated that these parameters had negligible contribution in determining the accuracy of the solution. A final set of sensitivity simulations were performed by varying the value of the advection velocity (v) by an order of magnitude above and below the base line value. The results of this sensitivity analysis are summarized in Figures 10a and 10b. Initial observations of these figures may indicate that at higher velocities, the modified- Domenico solution appears to perform better. However, a closer analysis of these figures with respect to their respective advective front locations indicates that at higher velocities a greater portion of the plume is behind the advective front whereas, at lower velocities a relatively lesser portion of the plume is behind the advective front confirming that the advection velocity has little effect in determining the accuracy of the solution. However, it must be noted that the advection velocity itself plays an important role in determining the location of the advective front, which is one of the key parameters that affects the performance of the solution. Variations in the total simulation time (t) will have a similar effect as that of the advection velocity. 82 (a) (b) Fig. 9: Sensitivity results for variations in the longitudinal dispersivity value; solutions include concentration contours beyond the advective front for (a) ?x*10 (b) ?x/10 83 (a) (b) Fig. 10: Sensitivity results for variations in the transport velocity; solutions include concentration contours beyond the advective front for (a) v*10 (b) v/10 84 From the results of these sensitivity simulations, it can be safely concluded that, the two most important factors that affect the accuracy of the modified-Domenico solution are the value of the longitudinal dispersivity (?x) and the position of the advective front (v*t). The solution will have minimal errors when the value of ?x is low and when the advective front is farther away from the source. It must be noted that the conclusions obtained for the modified-Domenico solution apply to the Domenico solution as well, provided the value of x x ? is sufficiently large [5]. 4.8 Discussions and Recommendations Since the original Domenico solution lacked a theoretical basis, several misconceptions regarding its performance have evolved over the years. One of the common misconceptions is that the error will be a minimum along the plume centerline. For example, Guyonnet and Neville [22] compared the Domenico solution against the Sagar [39] solution and concluded that ?The results of the evaluation confirm that along the plume centerline, and for ground water flow regimes dominated by advection and mechanical dispersion rather than by molecular diffusion, discrepancies between the two solutions (namely the Domenico solution and the Sagar solution) can be considered negligible for all practical purposes. However the errors in the Domenico (1987) solution may increase significantly outside the plume centerline.? However, the above simulation results indicate that this conclusion might not be true for all cases. To illustrate this, compare the y and z concentration transects predicted by the two solutions for our base case scenario. Figure 11a compares the concentration profiles along y 85 direction at x = 1000 m and 1500 m, and similar results for the z direction are shown in Figure 11b. It is evident from these figures that the error is not minimal along the centerline, but rather at a point which will always be away from the centerline. This error pattern can also be observed in all the two dimensional contours. Further, it can be observed from Figures 11a and 11b that the absolute error is, in fact, maximum along the plume centerline. Another important issue that should be addressed here is the nature in which the error associated with the Domenico approximation is propagated spatially. The sensitivity results presented in this chapter show that the plumes predicted by the modified-Domenico solution are always wider than the actual plumes. This phenomenon can easily be observed in all the figures presented in this study. This can be attributed to the fact that the Domenico approximation over predicts the conceptual residence times of all particles along the centerline (hence allows more time to disperse in the transverse directions). This over prediction would lead to a decrease in the centerline concentrations; therefore, solutions that employ the Domenico approximation will always under predict the overall extent of the plume in the longitudinal direction. An important transport parameter not addressed so far is the retardation factor (R). Retardation affects the advection velocity and possibly the decay constant (depending on the phase where the decay occurs). Since the presence of a decay term does not introduce any significant additional error to the Domenico solution, its effect can be ignored. However, retardation changes the location of the advective front by changing the advection velocity and hence would influence the overall accuracy of the Domencio solution. 86 0 50 100 150 200 250 0 100 200 300 400 500 Y (m) Co nc en tra tio n ( mg /L) Domenico (1000m) Wexler (1000m) Domenico (1500m) Wexler (1500m) (a) 0 50 100 150 200 250 0 5 10 15 20 Z (m) Co nc en tra tio n ( mg /L) Domenico (1000m) Domenico (1500m) Wexler (1000m) Wexler (1500m) (b) Fig. 11: Concentration profiles predicted by the Domenico solution compared with the Wexler solution at x = 1000m and 1500m (a) along Y axis (b) along Z axis 87 Based on the theoretical results presented in this study, it can be concluded that the key assumption used to derive the Domenico solution is the time reinterpretation step where the time ? in the transverse dispersion terms is replaced with x/v. The derivations presented in section 4.5 prove that this substitution process is valid only when the longitudinal dispersivity is zero. For all non-zero longitudinal dispersivity values the solution will have a finite error. The spatial distribution of this error is highly sensitive to the value of ?x and the position of the advective front (v*t), and is relatively less sensitive to other transport parameters. Based on the results of this study, it can be concluded that the error in the Domenico solution will be low when solving transport problems that have low longitudinal dispersivity values, high advection velocities, and large simulation times. Despite its limitations, the Domenico approximation offers a simple alternative for extending one dimensional analytical solutions to three dimensions. This approach is useful for developing approximate solutions for unsolved, three dimensional, multi- species reactive transport problems that have explicit one dimensional solutions. However, such solutions should be used carefully after understanding the limitations identified in this study. 88 CHAPTER V SUMMARY AND CONCLUSIONS In this research a set of analytical solutions to multi-species reactive transport equations coupled through sorption and sequential first-order reactions is presented. In the first part of this work, the mathematical derivations of the general solution which incorporates a generic Bateman type exponentially decaying Dirichlet and Cauchy source boundary conditions and a spatially varying initial condition are presented. The solution strategy involves uncoupling the system of governing equations through the use of transformations and then solving the uncoupled system of equations individually to obtain independent solutions for each species in the transformed domain. A combination of Laplace and linear transforms were used to uncouple the system of equations. Finally the independent solutions in the transformed domain are analytically re-transformed back to the original time domain to obtain the final solutions. Some of the key challenges in performing these analytical operations include formulation of the generic linear transform matrices and evaluation of the inverse Laplace expressions. The solutions to both the Dirichlet and Cauchy boundaries are presented in a common format to enable simultaneous implementation of both these solutions. Furthermore, the solutions are presented in a closed form format that avoids numerical integrations processes. Solutions for two cases of sorption; one involving decay in liquid 89 phase only and the other involving decay in both liquid and solid phases are also presented. In the second part of this research, the computational techniques for implementing the solution are discussed and the general solution is validated against other published solutions. Some of the key challenges faced when implementing the solutions are round- off errors and underflow/overflow errors. A combination of techniques involving algebraic manipulation and logarithmic transforms are suggested to tackle these errors. A FORTRAN computer code ?SEQUENTIAL? that solved the general solution was developed and is provided in Appendix I. The new code was used to simulate four example problems and the results generated by the code matched the results of previously published analytical and semi-analytical solutions. In addition, several special-case solutions for simpler transport problems involving zero initial condition, identical retardation factors, zero advection, zero dispersion and steady-state condition are also presented. Where ever possible these special case solutions were successfully verified against previously published analytical solutions. Strategies for extending the one dimensional sequential solutions to a generic diverging reaction network and to multiple dimensions (using the approximate Domenico solution) are also presented. A detailed investigation into the errors involved in this approximation is also discussed. The solutions proposed in this research work can be used to develop efficient screening tools for assessing ground water quality issues at sites contaminated with radioactive materials or chlorinated solvents. One of the key limitations of the general solution is that it can model sequential/diverging reaction networks only. More generic reaction networks involving 90 converging and/or reversible reaction networks cannot be modeled using the general solution. The main challenge in incorporating a more generic reaction network is obtaining the analytical linear transform matrices. Simplifying the problem to steady state may perhaps help in reducing the complexity in obtaining the matrices. It must be noted that Laplace transform techniques were used in this study to uncouple the system of equations. The use of Fourier transforms for this purpose may possibly help in easing some of the other limitations on the general solution including the restriction on the transport parameters (See Table 4) and the limitation to modeling first- order kinetics. Future research efforts in this area should focus on implementing more robust mathematical induction techniques and using Fourier transforms to obtain solutions to more general multi-species transport problems. 91 REFERENCES [1] Angelakis AN, Kadir TN, Rolston DE. Analytical solutions for equations describing coupled transport of 2 solutes and a gaseous product in soil. Water Resources Research 1993;29:945-956. [2] Angelakis AN, Kadir TN, Rolston DE. Solutions for transport of 2 sorbed solutes with differing dispersion coefficients in soil. Soil Science Society of America Journal 1987;51:1428-1434. [3] Aziz CE, Newell CJ, Gonzales JR, Hass P, T.P C, Sun Y. BIOCHLOR: Natural attenuation decision support system V.1.0. User's manual. US EPA Report, EPA/600/R- 00/008. Cincinnati, Ohio. 2000; [4] Bauer P, Attinger S, Kinzelbach W. Transport of a decay chain in homogenous porous media: analytical solutions. Journal of Contaminant Hydrology 2001;49:217-239. [5] Bear J. Hydraulic of Groundwater. McGraw-Hill, New York 1979;268-269. 92 [6] Bruch JC, Street RL. Two-dimensional dispersion. Journal of the Sanitary Engineering Division Proceeding of the American Society of Civil Engineers 1967;93:17-39. [7] Burkholder HC, Rosinger ELJ. A model for the transport of radionuclides and their decay products through geologic media. Nuclear Technology 1980;49:150-158. [8] Carslaw HS, Jaeger JC. Conduction of heat in solids. Oxford, Great Britan: Oxford university Press; 1980. [9] Cho CM. Convective transport of ammonium with nitrification in soil. Canadian Journal of Soil Science 1971;51:339-350. [10] Cleary R, Ungs. MJ. Analytical models for groundwater pollution and hydrology, Water Resources Program, Princeton University, 1978. [11] Clement TP. Generalized solution to multispecies transport equations coupled with a first-order reaction network. Water Resources Research 2001;37:157-163. [12] Clement TP, Sun Y, Hooker BS, Petersen JN. Modeling multispecies reactive transport in ground water. Ground Water Monitoring and Remediation 1998;18:79-92. 93 [13] Clement TP, Truex MJ, Lee P. A case study for demonstrating the application of U.S. EPA's monitored natural attenuation screening protocol at a hazardous waste site. Journal of Contaminant Hydrology 2002;59:133-162. [14] Cody WJ. Rational Chebyshev approximations for the error function. Mathematics of Computation 1969;23:631-638. [15] Crank J. The mathematics of diffusion. New York: Oxford university press; 1975. [16] Domenico PA. An analytical model for multidimensional transport of a decaying contaminant species. Journal of Hydrology 1987;91:49-58. [17] Domenico PA, Robbins GA. A new method of contaminant plume analysis. Ground Water 1985;23:476-485. [18] Eykholt GR, Li L. Fate and transport of species in a linear reaction network with different retardation coefficients. Journal of Contaminant Hydrology 2000;46:163-185. [19] Gautschi W. Efficient computation of the complex error function. SIAM Journal on Numerical Analysis 1970;7:187-198. [20] Gautschi W. Error function and fresnel integrals. Washington D.C.: U.S. Government Printing Office; 1964. 94 [21] Gureghian AB, Jansen G. One-dimensional analytical solutions for the migration of a 3-member radionuclide decay chain in a multilayered geologic medium. Water Resources Research 1985;21:733-742. [22] Guyonnet D, Neville C. Dimensionless analysis of two analytical solutions for 3- D solute transport in groundwater. Journal of Contaminant Hydrology 2004;75:141-153. [23] Harada M, Chambre PL, Foglia M, Higashi K, Iwamoto F, Leung D, Pigford TH, Ting D. Migrations of radionuclides through sorbing media, analytical solutions - 1. Report no. LBL-10500 (UC-70, Lawrence Berkely Laboratory, University of California Berkely 1980; [24] Higashi K, Pigford TH. Analytical Models for Migration of Radionuclides in Geologic Sorbing Media. Journal of Nuclear Science and Technology 1980;17:700-709. [25] Hunt B. Dispersive sources in uniform ground-water flow. Journal of the Hydraulic Division 1978;104:75-85. [26] Jones NL, Clement TP, Hansen CM. A three-dimensional analytical tool for modeling reactive transport. Ground Water 2006;44:613-617. 95 [27] Khandelwal A, Rabideau AJ. Transport of sequentially decaying reaction products influenced by linear nonequilibrium sorption. Water Resources Research 1999;35:1939-1945. [28] Kreyszig E. Advanced engineering mathematics. New York, USA: John Wiely & Sons; 2004. [29] Lester DH, HJansen G, Burkholder HC. Migration of radionuclide chains through an adsorbing medium. Adsorption and Ion Exchange, AICHE Symposium Series 152, Ameriican Institute of Chemical Engineers, New York 1975;71: [30] Lunn M, Lunn RJ, Mackay R. Determining analytic solutions of multiple species contaminant transport, with sorption and decay. Journal of Hydrology 1996;180:195-210. [31] Martyn-Hayden J, Robbins GA. Plume distortion and apparent attenuation due to concentration averaging in monitoring wells. Ground Water 1997;35:339-346. [32] McLaren AD. Nitrification in soil - systems approaching a steady state. Soil Science Society of America Proceedings 1969;33:551-&. [33] McLaren AD. Temporal and vectorial reactions of nitrogen in soil - a review. Canadian Journal of Soil Science 1970;50:97-&. 96 [34] Misra C, Nielsen DR, Biggar JW. Nitrogen transformations in soil during leaching .1. Theoretical considerations. Soil Science Society of America Journal 1974;38:289-293. [35] Montas HJ. An analytical solution of the three-component transport equation with application to third-order transport. Water Resources Research 2003;39: [36] Newell CJ, McLeod RK, Gonzales JR. BIOSCREEN: Natural attenuation decision support system user?s manual EPA/600/R-96/087, Robert S. Kerr Environmental Research Center, 1996. [37] Quezada CR. Generalized Solution to Multi-Species Transport Equations Coupled with a First-Order reaction Network with Distinct Retardation factors, Auburn University, Auburn, AL, USA, 2004. [38] Quezada CR, Clement TP, Lee KK. Generalized Solution to Multi-dimensional Multi-species Transport Equations Coupled With a First-order Reaction Network Involving Distinct Retardation Factors. Advances in Water Resources 2004;27:508-521. [39] Sagar B. Dispersion in three dimensions: Approximate analytical solutions. Journal of the Hydraulic Division, Proceeding of the American Society of Civil Engineers. 1982;108:47-62. 97 [40] Selby SM. Standard mathematical tables. Cleveland, Ohio: The Chemical Rubber Co.; 1971. [41] Srinivasan V, Clement TP, Lee KK. Domenico solution - Is it valid? Ground Water 2007;45:136-146. [42] Sun Y, Petersen JN, Clement TP. Analytical solutions for multiple species reactive transport in multiple dimensions. Journal of Contaminant Hydrology 1999;35:429-440. [43] Sun Y, Petersen JN, Clement TP, Skeen RS. Development of analytical solutions for multispecies transport with serial and parallel reactions. Water Resources Research 1999;35:185-190. [44] Tim US, Mostaghimi S. Modeling transport of a degradable chemical and its metabolites in the unsaturated zone. Ground Water 1989;27:672-681. [45] van Genuchten MT. Convective-dispersive transport of solutes involved in sequential 1st-order decay reactions. Computers and Geosciences 1985;11:129-147. [46] Wark K, Warner CF. Air Pollution, Its Origin and Control. New York: Addison Wesley Longman Inc; 1981. 98 [47] West M, Kueper BH. Natural attenuation of solute plumes in bedded fractured rock. USEPA/NGWA Fractured Rock Conference,Portland, Maine, National Ground Water Association 2004;388?401. [48] Wexler EJ. Analytical solutions for one-, two-, and threedimensional solute transport in ground-water systems with uniform flow, USGS?TWRI Book 3, 1992. 99 APPENDIX A Factorization of the Laplace Terms for the Dirichlet Boundary Condition The first term ? ( )1 1 ,ip x s?? ? can be evaluated as follows. From equations (41) and (42) we get: ( ) ( ) ( ){ } ( ) ( ){ } ( ) ( ) 22 1 1 2 2 22 1 2 2 2 2 1 2 1 2 2 3 3 3 1 3 2 1 11 1 1 1 1 42 1 , 1 , , o i x i i x t si ii i i i ii i ii xi i v v D k sR Dii i i i i i i i i i i i B e y k s c x t p x s e k sR k sR ? ? ? + ? == + ? ? ? ? ? + +? ? ? ?= = = ? ? ??? ? ? ?? ? ? ?+? ? ? ? ? ?= = ? ? ? ? ? ?? + ? ? ? ?? ? ?? ? ? ? ? ? 1, 2,...i n ; ? = (A.1) Equation(A.1) can be rewritten as: 100 ( ) ( ){ } ( ){ } ( ) ( ) ( ) 2 2 2 1 2 2 23 3 1 11 2 1 3 3 2 4 2 4 4 1 4 2 1 1 1 1 42 1 1 ,0 , , , , 1 x i io i x i i i i i i x v v D k sRt s Di i ii ii i i i i i i i i i i i i i y k c x t B e e s a R s a ? ? = + ? ?? ? + + ? ?? + ? ? = = = = ? ? ?? ? ? ?? ? ? ?? ? ? ? = ? ?? ? ? ? ? ? ?+ ? + ? ?? ? ? ? ?? ? ? 1 2 1 21 2 1 2 1 2 1 2 1 2 1 , 2 ,, , , 2 1, 2,... ; 0 ; ; 0 i i i ii i i i i i i i i i i i n k when i Rwhere k k k R R R and a when i? ; ? = ? > ? = ? , = ? = ? ? =? (A.2) It must be noted that for equation(A.2) to be valid the condition ? 2 4, 0i iR ? ? must be satisfied. This means that no two species in the transport problem can have identical retardation factors. However, in practice, we do have situations where the retardation factors of some of the species are equal [9]. To overcome this limitation, we reformulate equation(A.2) to accommodate a generic case when the transport problem has any number of sets of species having any number of species with identical retardation factors. Incorporating this special case scenario, equation(A.2) is reformulated as: 10 1 ( ) ( ){ } ( ){ } ( ) ( ) ( ) ( ) 2 2 2 1 2 2 233 1 1 1 2 1 3 3 2 4 2 4 2 4 4 1 4 2 4 1 4 24 2 4 2 1 1 41 1 2 1 1 ,0 , , , , , , , , 1 x i io i x i i i i i i i i i xi v v D k sR t s Di i ii i i i ii i i i i i i i i i i i i i R R i i i i R R y k c x t B e e s a k R s a ? ? = + ? ?? + + ? ?? ? + ? ? = = = = ? = = ? ? ? ?? ? ? ?? ? ? ?? ? ? ? ? ?= ? ? ? ? ?? ? ? ?? ?+ ? ? + ? ?? ?? ? ? ?? ? ? ? ?? ? ? ? 1, 2,...i n ; ? = (A.3) Note that in equation(A.3) the condition ? 2 4, 0i ik ? ? must be satisfied for the solution to be determinate. Factorization of equation(A.3) gives: 10 2 ( ) ( ){ } ( ){ } ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 1 2 2 23 3 1 3 2 4 4 1 4 2 4 2 4 1 4 2 4 2 2 4 2 4 2 5 2 5 2 4 5 1 5 2 5 4 5 2 1 1 42 1 1 ,0 , , , , , , , , , , , , , 1 , 1 x i io i x i i i i i i i i i i i x v v D k sR t s Di i i i i i i i i i i R R i i i i i i R R i i i i i i i i i i i i i i i i R R y k B e e c x t s a k R s a R a a ? ? = + ? ?? + + ? ?? + ? ? ? = ? = = ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? = ? + ? ? ? + ? ? ? ? ? ? ? 11 2 1 3 1 1 1, 2,... i iii i i i i n = = = ? ? ? ? ? ? ? ?? ? ?? ? ?? ? ?? ? ?? ? ? ?? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ?? ? ; ? = ? ?? (A.4) It must be noted that the solution formulation as given by equation(A.4) is valid only when the condition ? 2 5 2 4, ,i i i i a a? ? is satisfied. Equation(A.4) can be further factorized and simplified as: 10 3 ( ) ( ){ } ( ){ } 1 2 2 1 2 1 32 1 2 3 3 1 1 1 1 1 1 1 1 2 1 11 42 2 1 1 1 , 1, 2,... ; 1 x io i x x iii i i i i i i i ii i x v D k t s Dxv i D i c x t y k G h G G i n where e e B e G ? ? = = == + ? + ? + ? ? ?? ?= + ? ?? ?? ? ? ?? ? ; ? = ? = ? ??? ? ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ){ } ( ) ( ) 2 2 3 2 4 4 1 4 2 4 2 2 5 2 4 2 5 2 5 2 4 5 1 5 2 5 1 5 2 5 45 2 5 2 2 2 23 3 1 3 ,0 , , , , , , , , , , , , , 42 2 1 ,0 1 2 1 i i i i i i i x i io i x x sR i i i i i ii i i i R R i i i i i i i i i i i i i i R R i i i i i i R R x v D k sR t s Dxv i D i i s a s a k R R a a e e B e s a G ? ? ?+ ? ?? ? = ? ? = ? = = ? ? ? ? ?? + + ? ?? + ? ? ? + + ? ? ? ?? ? ? ?? ? ? ? ? + = ? ? ? ? ( ) ( ) 2 4 4 1 4 2 4 2 , , , 1 0 i i i i i i i i i R R k if loopis not executedh M if loopis executed = ? = ? ? ?= ? ? ? ? (A.5) 10 4 The term ? 11G ? in equation(A.5) can be further factorized and simplified as: ( ){ } ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) ( ) 2 2 23 3 3 1 2 2 23 2 4 2 4 3 2 5 2 4 2 5 2 5 1 5 2 5 2 42 1 ,0 2 42 1 , 1 1 , , , , , , 1 1 x i io i x x x i io i x i i x v D k sR t s D i xv i D i x v D k sRt s D i i i i i i i i i i i i i i i i i R R e e s a B e e e s a G a k R R a ? ? ? ? ?? + + ? ?? + ? ? ? ? ?? + + ? ?? + ? ? ? = ? = ? ? ? ?? ? ? ? ?+ ? ? ? ? ? ? ? ?? ? ?? ? ?+ ? ? ? ?? ?= ? ? ? ?? ? ? ?? ? ? ? ? ? ? ( ) ( ) ( )4 1 4 2 5 2 5 2 4 5 1 5 2 5 4 5 2 , , , , , , , i i i i i ii i i i R R i i i i i i i i i R R a= ? = = ? ? ? ? ? ? (A.6) Again it must be noted that equation(A.6) is valid only when the condition ? 2 4 3, ,0i i i a a? ? where ? 3 3,0i i a ?= ? is satisfied. Using the results from Appendix B, inverse Laplace expressions for the terms ? 11G ? and ? 12G ? can be evaluated and the solution for ? ( )1 ,ic x t ? is obtained as: 10 5 ( ) ( ){ } [ ] ( ) ( ) ( ) 1 2 2 1 2 1 32 1 3 2 3 2 33 1 1 1 1 1 1 1 1 2 1 11 , ,0 , ,0 1 1 , 1, 2,... , ,oi iii i i i i i i i ii i t i i o i i oi i c x t y k G h G G i n F x t u t t e F x t t B G ? ? = = == + ? ? ?? ?= + ? ?? ?? ? ? ?? ? ; ? = ? ?? ? ?? ? = ? ??? [ ] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ] ( ) ( ) 3 2 2 4 2 2 4 4 1 4 2 4 2 2 4 3 2 5 2 4 2 5 2 5 2 4 5 1 5 2 5 1 5 2 5 45 2 5 2 33 1 2 3 , , , , , , , , , , , , , , , , , , ,0 1 2 , , , oi i i i i i i oi t i i i i o i i i o i ii i i i R R i i i i i i i i i i i i i i i i i R R i i i i i i R R ti i i i o F x t u t t e F x t t a k R R a a B F x t u t t e F G ? ? ? ? = ? ? = ? = = ? ? ? ? ? ?? + ? ?? ? ? ? ? ?? ? ? ? ?? ? ? ? ? ? = ? ? ? ( ) ( ) 2 3 2 4 4 1 4 2 4 2 , ,0 , , , , i i i i o i i i i i i i R R x t t k = ? = ? ??? ? ?? (A.7) where; the term ? 1 2 3, ,i i i F ? is given by equation (B.8) or (B.14) in Appendix B. The second term ? ( )1 2 ,ip x s?? ? is evaluated as follows. From equations (41) and (42) we get: 10 6 ( ) ( ) ( ){ } ( ) ( ) ( ) 1 1 2 2 2 1 2 2 2 1 1 2 1 1 1 2 2 2 2 3 3 3 1 3 2 2 1 2 1 1 41 2 1 2 , , , x i i ix i i i o i i i i i i xi v v D k sR xD i i i i i i x i i i i i i i i i i i c x t p x s R c y k e e D v k sR k sR k sR ? ? ? ? ? = + ? ?? + + ? ?? ?? ? = = = ? = ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? = ?? ? ? ?? ? ? ?? ? ? ? ? ?+ ? ? ? + ? ? ? ?? ? ? ? ? ? ? ? 1, 2,...i n ; ? = (A.8) Equation(A.8) can be written as: ( ) ( ){ } ( ) ( ) ( ) ( ) 1 1 2 2 2 1 2 2 2 1 1 2 1 1 2 2 2 3 2 3 2 3 3 1 3 2 3 1 3 23 2 3 2 1 1 42 2 1 , , , , , , , , , x i i ix i i i i i o i i i i i i x v v D k sR xD i i i i ii i i i i i i i i i i i i i i R R i i i i R R R c y k e ec x t s a R k R s a ? ? = + ? ?? + + ? ? ?? ? ? = ? = ? = = ? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ??= ? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ?? + ? ? + ? ?? ?? ?? ? ? ? ? ? ? 1 1, 2,... i i n = ; ? = ? (A.9) 10 7 Note: the term ? 1 2,i i a ? is modified as: 1 2 1 2 1 2 1 1 1 2 2 , 2 , , 2 2 2 ; 0 ; 0 ; 0 i i i i i i i i x i i i k when i R a when i D v k when i R ? ? ? ? >? ?? = =? ?? ? + ? < ?? (A.10) Note that the condition ? 2 3, 0i ik ? ? must be satisfied in equation(A.9). Factorization of equation(A.9) yields: 10 8 ( ) ( ){ } ( ){ } 1 1 2 2 1 2 12 1 2 2 2 2 2 2 2 1 1 1 2 1 1 42 2 1 2 1 , 1, 2,... ; x i i ix x ii i o i i i i i i i ii i x v v D k sR D xv D c x t R c y k G h G G i n where e e e G ? ? = == + ? ?? + + ? ? ?? ? ? ? ?? ?= + ? ?? ?? ? ? ?? ? ; ? = ? = ? ?? ? ( )( ) ( ) ( ) ( ) ( ) ( ){ } 1 1 2 2 3 3 1 3 2 3 2 2 2 4 2 3 2 4 2 4 2 3 4 1 4 2 4 1 4 2 4 34 2 4 2 2 2 2 1 , , , , , , , , , , , , , , 42 2 1 2 2 i i i i i i x i i ix x x i i i i i i ii i i i R R i i i i i i i i i i i i i i i R R i i i i i i R R x v v D k sR xD xv D s a s a R k R R a a e e e G ? ? = ? ? = ? = = ? ? ? ? ?? + + ? ? ?? ? ? ? ?? ? ? ? ? ?? ? + + ? ? ? ?? ? ?? ? ? ? ? ?? ? ?? ? ? ?? ?? = ? ? ? ? ( ) ( ) 1 2 2 2 3 3 1 3 2 3 2 , , , , i i i i i i i i i i i i R R s a R k ? = ? = + ?? (A.11) 10 9 Note that the condition ? 2 4 2 3, ,i i i i a a? ? must be satisfied in equation(A.11). The term ? 21G ? in equation(A.11) can be further factorized and simplified as: ( ){ } ( ) ( ){ } ( ) ( ) ( ) ( ) 2 2 2 2 2 2 1 1 1 2 2 3 2 3 1 2 2 2 4 2 3 2 4 2 4 2 3 4 1 4 2 44 2 4 42 2 2 1 , , 2 1 , , , , , , , , , x i i x i i i ix x x i i x xv v D k sR v v D k sR x xD D xv D i i i i i i i i i i i i i i i i i i i i i i i i R R i e e e e e s a s a G a a R k R R a a ? ? ? ? ? ?? + + ? + + ? ? ? ?? ?? ? ? ? ? ? ? = ? = ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? ? ?? ? ? ?? + + = ? ? ? ?? ? ? ?? ? ? ? ? ? ( ) ( )3 1 3 2 3 2 1 4 2 4 3 4 2 , , , , , i i i i i ii i i i R R i i i i i R R = ? ? = ? ? ? ? ? (A.12) Note that the condition ? 2 4 1 2, ,i i i i a a ?? ? must be satisfied in equation(A.12). ? 21G ? can be further simplified as: 11 0 ( ){ } ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) 2 2 2 1 1 2 1 2 2 2 2 1 2 3 2 3 2 3 1 2 2 2 4 4 1 4 2 4 2 42 1 1 , , 2 42 1 1 , , 2 1 , , , , , 1 1 x i i x i x x i i x i i i x v v D k sR D x i i i ixv D x v v D k sR D x i i i i i i i i i i i i i i i i R R e e s a s a e e e s a s a G a a R k ? ? ? ?? + + ? ?? ? ?? ? ? ? ? ?? + + ? ?? ? ?? ? ? = ? = ? ? ? ?? ? ?+ + ? ? ? ? ? ? ? ? ? ?? + ? ?+ + ? ?? ?= ? ? ?? ?? ? ? ? ? ? ? ? ( ) ( ) ( )3 1 3 2 3 2 2 3 2 4 2 4 2 3 4 1 4 2 4 3 4 2 , , , , , , , , , i i i i i ii i i i R R i i i i i i i i i i i i i i R R R R a a= ? ? = ? ? ? ? ? ?? ? ? (A.13) ? 22G ? in equation(A.11) can be further simplified as: ( ){ } ( ) ( ) ( ) 2 2 2 1 1 2 1 2 2 2 3 3 1 3 2 3 2 42 2 1 1 , , 2 2 , , , 1x i ix ix i i x v v D k sR xv D xD i i i i i i i i i i i i R R ee e s a s a G R k ? ? ?? + + ? ?? ? ?? ? ? ? = ? = ? ? ? ?? ? ? ?+ + ? ? ? ?= ?? ? ? (A.14) From inverse Laplace transform tables we get: [8] (p494, eq3). 11 1 ( ) ,1 2 1 2 1 , 1 i ia t i i es a ?? =+? (A.15) Using equation(A.15) and Appendix B, inverse Laplace expressions for the terms ? 21G ? and ? 22G ? can be evaluated and the solution for ? ( )2 ,ic x t ? is obtained as: ( ) ( ){ } [ ] ( ) 1 1 2 2 1 2 12 1 ,1 1 2 2 1 2 2 2 3 2 2 2 2 1 1 1 2 1 1 , , , ,2 1 , 1, 2,... ; , ,i i i ii i o i i i i i i i ii i x a t i i i i i i c x t R c y k G h G G i n where F x t e F x G ? ? ? = == + ? ? ? ? ?? ?= + ? ?? ?? ? ? ?? ? ; ? = ? ? = ? ?? [ ] ( ) ( ) ( ) ( ) ( ) ( ) [ ] ( ) ,1 2 3 3 1 3 2 3 2 2 3 1 2 2 2 4 2 3 2 4 2 4 2 3 4 1 4 2 4 1 4 2 4 34 2 4 2 ,1 1 2 2 1 2 2 2 3 3 1 3 2 , , , , , , , , , , , , , , , ,2 2 , , , , i i i i i i i i i i i i x a t i i ii i i i R R i i i i i i i i i i i i i i i i i i i R R i i i i i i R R x a t i i i i i i i i i i R t e a a R k R R a a F x t e G R k ? ? ? ? ? = ? ? ? = ? = = ? ? ? ? ? ? = ? + ? ? ? ?? ? ? ?? ? ? ? ? ? = ? ? ? ? ( )3 2i i i R= ? (A.16) 112 where; the term ? 1 2 3, ,i i i F ? is given by equation (B.8) or (B.14) in Appendix B. 113 APPENDIX B Evaluation of Inverse Laplace Expressions for the Dirichlet Boundary Condition ( ){ } ( ) ( ) 2 1 14 1 2 3 4 2 3 42 2 1 , , , , 1 x i io i x x x v D k sR t s Dxv D i i i i i i e e e s a ? ? ? ?? + + ? ?? + ? ? ? ? = +? (B.1) Equation(B.1) can be simplified as: ( )1 2 4 3 2 1 1 1 1 1 1 2 3 1 1 1 1 2 1 1 1 14 2 , , , 1 2 4 1 1 2 2 , 4 1 2 ; 4 4 x i i x i x i i i x i x io i o xv D i i i i R kvx s D R D R i i i i i x i i x i R kvx s D R D Rt t s e where e k kv vs a R D R R D R e e e s ? ? ? ? ? ? ? ?? ?? ? ?? ?? ? ?? ? ? ?? ?? ? ? ? ?? ?? ? ?? ?? ? ?? ? ? ?? ? ?? ? ? ? = ? = ? ?? ? ?? ?? ? ? ? +? ?? ? ? ?? ?? ?? ? = ? ? 1 1 2 3 1 1 1 1 2 2 ,4 4 i i i i i x i i x i k kv v a R D R R D R ? ?? ??? ?? ? ? ? + ? ?? ?? ? ? ?? ?? ? (B.2) ? 1? ? can be evaluated as follows: Invoking the ?First Shifting Theorem? (see [28], p253) we can simplify ? 1? ? as: 114 12 1 1 1 1 2 3 1 1 4 1 1 2 ,4 i i x i x i Rx s kv Dt R D R i i i i x i ee kvs a R D R ? ? ? ? ??? ?? ? ??? ? ? ? ? ? ?? ?= ? ?? ?? ?? + ? ? ?? ?? ? ? ?? ?? ? ? (B.3) The Laplace inverse of equation(B.3) can be readily obtained from the tables (see [8], p495, eq19). Applying this inversion, ? 1? ? can be evaluated as: 2 1 1 , 2 3 1 1 1 1 2 3 ,2 3 1 1 2 1 1 , 2 3 1 1 1 1 2 3 1 1 24 , 1 24 , 2 4 2 2 4 i i i i x i x i i i i i i i x i x i R kvx a D R D R i i i i a t x i x i R kvx a D R D R i i i i x i x i R kx ve erfc a t D t R D Re R kx ve erfc a t D t R D R ? ? ?? ? + ? ? ? ?? ? ? ? ?? + ? ? ? ?? ? ? ? ?? ?? ? ? ? + ?? ?? ? ? ?? ?? ? ? ?= ? ?? ?? ? + + + ?? ?? ?? ? ? ?? ?? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? (B.4) ? 1? ? can be further rearranged and simplified as: 1 2 1 1 1 2 3 ,1 2 3 11 1 ,2 3 1 2 1 1 1 2 3 ,1 2 3 11 1 2 ,4 2 1 2 ,4 2 4 2 2 4 2 i i x i i x i i i i i x i i x i i kx i i x i iv R D a iD R i x a t i kx i i x i iv R D a iD R i x kR x v R D a t R e erfc R D t e kR x v R D a t R e erfc R D t ? ? ?? + ? ? ? ? ?? ? ? ? ?+ ? ? ? ? ?? ? ? ?? ? ? ?? + ? ? ?? ? ? ?? ?? ? ? ? ? ? ? ?? ? ? ?= ? ? ? + + ? ? ?? ? ? ?+ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ?? ? ? ?? ? ? ?? ?? ? ? ?? ? ? ?? ? ? ?? ?? ? ? ?? ? (B.5) ? 2? ? is evaluated as follows: ? 2? ? can be simplified as; 115 ( )4 12 1o i ot t se e?? ?? ??= ? ? (B.6) Now, we invoke the ?Second Shifting Theorem? (see [28], p265) and evaluate ? 2? ? as: [Note: We have already evaluated ? 1? ?; equation(B.4)] ( ) ( ) ( ) ( ) ,2 34 1 2 1 1 1 2 3 ,1 2 3 11 1 2 1 1 1 ,1 2 3 1 2 2 ,4 2 2 42 4 2 4 i i oi o i i x i i x i i i x i i x i a t tt o i kx i i x i i ov R D a iD R i x o i kx i i xv R D a D R e u t t e kR x v R D a t t R e erfc R D t t kR x v R D e erfc ?? ? ?? ? ?? + ? ? ? ? ?? ? ? ?+ ? ? ? ? ?? ? = ? ? ?? ? ? ?? + ? ?? ?? ? ? ?? ?? ? ? ?? ? ? ? ?? ? ? ? + + + ( ) ( ) 1 2 3 1 1 , 2 i i o i i x o a t tR R D t t ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ?? ?? ? ? ?? ?? ?? ?? ? ? ?? ?? ?? ? ? ??? ? ? ?? ? ? ?? ?? ? ? ?? ? (B.7) The final solution can be compactly expressed as: [ ] [ ] ( ) ( ) [ ] 4 1 2 3 4 1 2 3 1 2 3 , ,1 2 3 1 1 2 3 ,2 3 1 1 2 3 , ,1 2 3 1 1 2 3 1 , , , , , , , , ,2 2 , , , ,2 , , , ; 2 , 2 2 i o i i i x i i x i i i x t i i i i i i i o i i i o x i i i iD xv a t D i x i i i x i i i iD i x x t F x t u t t e F x t t where R x te erfc R D te e F x t R x te erfc R D t ? ? ? ? ? ? ? ? ? ? ?= ? ? ?? ? ? ?? ? ? ? ? ? ?? ?= ? ?+ ? ? + ? ? ?? ? 1 1 2 3 1 2 3 1 2 , , ,4 i i i i i x i i i kv R D a R? ? ? ? ? ? ? ? ? ? ? ? ? ? ??? ? ? ?= + ? ? ?? ? ? ? (B.8) 116 The above solution is valid only for real values of ? 1 2 3, ,i i i ? ?. For problems involving complex values for ? 1 2 3, ,i i i ? ? the ? 1 2 3, ,i i i F ? terms are evaluated as follows: 1 1 2 3 1 2 3 1 2 3 1 1 1 2 3 1 2 3 1 2 * , , , , , * 2 , , , ; 4 ; 4 i i i i i x i i i i i i i i i i i x i i i kLet v R D a i R kwhere v R D a R ? ? ? ? ? = + ? = ? ?? ? ? ? ? ? = + ? ? ?? ? ? ? (B.9) Now the ? 1 2 3, ,i i i F ? terms can be written as: [ ] * , ,1 2 3 1 1 2 3 ,2 3 1 *1 2 3 , ,1 2 3 1 1 2 3 1 * , ,2 2 , , * , ,2 2 , 2 2 i i i x i i x i i i x xi i i i iD xv a t D i x i i i xi i i i iD i x R x i te erfc R D te e F x t R x i te erfc R D t ? ? ? ? ? ? ? ?? ?? ? ? ? ?+? ? ? ?? ?? ? ? ?= ? ?? ?+ ? ? ? ?? ? ? ?? ?? ?? ? (B.10) From the symmetric relations we get: { } { }; ;, , ,if erfc a ib A iB then erfc a ib A iBwhere a b A B R + = + ; ? = ? ; ? (B.11) Also the exponent of a complex number can be expressed as: ( ) ( )( ) ( ) i i e Cos i Sin e Cos i Sin ? ? ? ? ? ? ? = ? = + (B.12) Using equations(B.11) and (B.12) equation(B.10) can be simplified as: 117 [ ] ( ) ( ) ( ) 1 2 3 1 2 3 ,2 3 1 2 3 1 2 3 1 2 3 1 * * , , , , 2 , , * * , , , , 2 2 , 2 2 2 ; i i x i i i i i ixv a t D x x i i i i i i i i i x x i x xCos i Sin A iB D De e F x t x xCos i Sin A iB D D where RA iB erfc ? ? ? ? ? ? ?? ?? ? ? ?? ? ? ?? ? +? ? ? ?? ?? ? ? ? ? ?? ?? ? ? ?? ?= ? ?? ? ? ? ? ?? ?? ?+ + ? ? ? ?? ?? ?? ? ? ? ? ?? ?? ? ? ?? ?? ? ? = ( )1 2 3 1 1 2 3 1 1 * * , , , , 2 2 i i i i i i i i x i x x i t R x i tand A iB erfc R D t R D t ? ?? ? ? ?? + ? ? ? ? ; + = ? ? ? ? ? ? ? ?? ? ? ? (B.13) Further simplification of equation(B.13) yields: [ ] ( ) , 1 2 3 1 2 32 3 1 2 3 1 1 2 3 1 2 3 1 1 1 2 3 1 * * , , , ,2 , , * 2 , , , * , , , 2 2 ; 4 2 i i x xv a t i i i i i iD i i i x x i i i i i x i i i i i i i i x x xF x t e e ACos BSin D D kwhere v R D a R R x i tand A iB erfc R D t ? ? ? ? ? ? ?? ? ? ?= ?? ?? ? ? ? ? ? ? ?? ?? ? ? ?? ? ? ? = + ? ? ?? ? ? ? ? ?+ ? ? + = ? ? ? ?? ? (B.14) Hence in problems where ? 1 2 3, ,i i i ? ? is complex the ? 1 2 3, ,i i i F ? term is given as by equation(B.14). 118 APPENDIX C Derivation of the General Solution for the Cauchy Boundary Condition Since the governing equations, the initial condition and the boundary condition at ?? ? are identical for the Dirichlet and the Cauchy boundaries, the procedures involved in uncoupling the system of equations will be identical and hence we can use the analysis in ?Section 2.2? to obtain the semi-determined general solution and the linear transform matrices ?[ ]A ? and ?[ ] 1A ? ?. The semi-determined general solution [identical to equation(33) is given by: ( ) ( ) ( ) ( ) ( ) 2 2 1 1 1 2 2 2 1 1 1 1 2 2 2 1 2 4 22 1 1 21 , , i i x xx i k sRx v v D DD i i nx o i i i ii i i n i i x i i i i i i i i i i i b x s e R c e y k D v k sR k sR k sR ? ? ? ? ?? ?+? ? ? ?? +? ? ? ?? ?? ?? ? ? ? = + = = ? = ? ? ? ? ? ? ? ? ? ?+ ? ? ? + ? ? ? ?? ? ? ? ? 1, 2,...i n ; ? = (C.1) The linear transform matrices ?[ ]A ? and ?[ ] 1A ? ? are given by equations(15) and (16) respectively. The constant ? 2i? ? in equation(C.1) is evaluated by using the boundary condition given by equation(47) after transforming them into the ?b? domain; this is done as follows. Laplace transform of equation(47) gives: 119 ( ){ } ( ){ } { } ( ){ } ( ) 11 1 11 0, 0, 1 ; 1, 2,... o i x t si i i i i i p sD v p s x B v e where i ns ? +? = ?? + = ? ? ? ? = ; ? = +?? (C.2) In order to transform the boundary condition from the ? p ? domain to the ?b? domain, we apply the linear transform given by equation(12). ( ){ } ( ){ } [ ] { }10, 0,x b sD v b s Ax ??? + = ?? (C.3) The explicit expression for ? [ ] { }1 i A ? ? ? in equation(C.3) is given as: [ ] { } ( ) ( ) ( ){ } ( ) 2 2 2 1 2 2 2 1 2 221 1 1 1 1 2 2 , 1 1 1 1 1, 2,... n o i i i i i n i i i i i i i i t si ii i i i i i y k k sR k sR v B e A s i n ? = + = ? ? +? ? = =+ ? ? ? ?? ? ? ? ?? = + ?? ?? ? ?? ? ; ? = ? ? ? ? (C.4) From equations(C.3) and (C.4) we obtain the boundary condition in the ?b? domain as: ( ) ( ) ( ) ( ) ( ){ } ( ) 2 2 2 1 2 2 2 1 2 221 1 1 1 1 2 2 , 1 1 10, 0, 1, 2,... n o i i i i i n i i i i i i i i t si ii i i x i i i i y k k sR k sR v B eb s D vb sx s i n ? = + = ? ? +? = =+ ? ? ? ?? ? ?? ? ?? + =? +?? ?? ? ?? ? ; ? = ? ? ? ? (C.5) Substituting the expression for ? ( ),ib x s ? from the semi-determined general solution given by equation(C.1) into the transformed boundary condition given by equation(C.5), we can evaluate the constant ? 2i? ? as follows: 120 ( ) ( ) ( ) ( ) ( ) ( ) 1 1 1 2 2 2 1 1 1 1 2 2 2 1 2 1 1 2 2 2 1 1 1 2 2 2 2 2 2 1 1 21 , 1 12 2 41 2 i i x i x xx n o i i i i ii i i x n i i x i i i i i i i i i i i n o i i i i i i i i x i i i i i i i i k sRv vD D DD R c y k D D v k sR k sR k sR R c y k v v D v k sR k sR k sR ? = + = = ? ? = + = ? ?+? ? ? ? ? +? ? ? ?? ? ? ??? ? ? ? ?+ ? ?? +? ? ? ? + ? ? ? ? ? ? + ? ? ? + ? ? ? ? + ? ? ? ? ? ? ( ) ( ) ( ) ( ){ } ( ) 2 2 2 1 2 2 2 1 2 1 1 2 221 1 1 1 1 2 2 , 1 , 1 1 1 1, 2,... n o i i i i i n i i i i i i i i i n i i i i t si ii i i i i y k k sR k sR v B e s i n ? = + = ? = ? ? +? = =+ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ?= ? ?+ ?? ? ?? ? ; ? = ? ? ? ? ? ? (C.6) From equation(C.6) ? 2i? ? is evaluated as: ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) ( ) 2 2 2 1 2 2 2 1 2 221 1 11 1 2 2 , 1 1 1 2 2 2 1 1 1 1 2 2 2 1 2 1 1 1 1 21 ,2 1 n o i i i i i n i i i i i i i i t si ii i i i i n o i x i i i ii i i n i i x i i i i i i i i i i i i y k k sR k sR v B e s v D R c y k D v k sR k sR k sR ? = + = ? ? +? = = ? = + = = ? + ? ? ? ? ? ?? ? ? ? ? ? ? ? ?+ ?? ? ? ?? ?? ? ? ? +?? ? ? ?+ ? ?? + ? ? ? ? + ? ? ? ? ? ??? = ? ? ? ? ? ? ? ( )21 42 2 1, 2,... x i i v v D k sR i n ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ?? ? ?+ + +? ? ? ? ; ? = (C.7) Therefore, the solution in the ?b?domain is: 12 1 ( ) ( ) ( ) ( ){ } ( ) ( ) ( ) 2 2 2 1 2 2 2 1 2 2 222 1 1 1 2 2 1 1 1 2 2 2 1 1 1 , 4 2 21 1 1 1 1 , 1 42 2 1 i i o i x xx n i i i i n i i i i i i i i k sRx v vt s i D DDi i i i i i i x i i i n o i i i i i i y k k sR k sR B e ve b x s vs v D k sR R c y k ? ? ? ? = + = ? ? ?? ?+? ? ? ?? +? + ? ? ? ?? ?? ?? ? = = ? = + + ? ? ? ?? ? ? ? ?= ? ?+ ? ?? + + +? ? ? ? ? ?? ? + + ? ? ? ? ? ( ) ( ) ( ) ( ) ( ) 2 2 1 1 1 1 2 2 2 1 2 4 2 2 21 , 1 4 2 2 i i x xx i k sRx v v D DDx x x i i i n i i x i i i i i i i i i i i D ve v e v v D k sR D v k sR k sR k sR ? ? ? ? ?? ?+? ? ? ?? +? ? ? ?? ?? ?? ? ? = = ? ? ?? ?? ? ? ?? ?? ? ? ?? ?? ?? ? ? ??? ?? ? ? ?? ?+ + +? ? ? ?? ?? ? ? ?? ?? ? ? ? ? ?+ ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? 1, 2,...i n ; ? = (C.8) Inverse linear transform of equation(C.8) is done to obtain the solution in the Laplace domain (? p ? domain) by using equation(12). The solution given by equation(C.8) can be split into two parts and represented as: 12 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) 2 2 2 1 2 2 2 1 2 2 222 1 1 1 2 2 1 1 1 , 1 2 4 2 1 21 1 2 , , , ; 1 , 1 42 2 , i i o i x xx n i i i i n i i i i i i i i i i i k sRx v vt s i D DDi i i i i i i x i i i i y k k sR k sR b x s b x s b x s where B e ve b x s vs v D k sR R c b x s ? ? ? = + = ? ? ?? ?+? ? ? ?? +? + ? ? ? ?? ?? ?? ? = =+ ? ? = + ? ?? ? ? ? ?= ? ?+ ? ?? + + +? ? ? ? ? ?? ? = ? ? ? ? ( ) ( ) ( ) ( ) ( ) 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 2 4 2 1 21 2 , 1 1 4 2 2 i i x xx i k sRx v v D DDi x n x o i i i i i x i i n i x i i i i i i i i i i i D ve vy k e v v D k sR D v k sR k sR k sR ? ? ? ? ? ?? ?+? ? ? ?? +? ? ? ?? ?? ?? ? ? ? = + = ? ? ?? ?? ? ? ?? ?+? ? ? ?? ?? ?? ? ? ??? ?? ? ? ?? ?+ + +? ? ? ?? ?? ? ? ?? ?? ? ? ? ? ?+ ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ?1 1 1, 2,... i i i n = ? ; ? = ? (C.9) 123 Using the distributive property of matrix addition, we can apply the inverse linear transform to each of the individual terms and then sum them to get the solution in the ? p ? domain. This is expressed as: { } [ ]{ } [ ]{ } [ ]{ }1 2p A b A b A b= = + (C.10) The first term ?[ ]{ }1A b ? can be evaluated as: { } [ ]{ }1 1p A b= (C.11) The explicit expression for ? ( )1 ,ip x s ? in equation(C.11) is given as: ( ) ( ){ } ( ) ( ){ } ( ) ( ) 22 1 1 2 2 22 1 2 2 2 2 1 2 1 2 2 2 2 3 3 3 1 3 2 1 11 1 42 1 2 , 1 , 1 4 ( ) 2 2 o i x i i x t si ii i i i ii i ii xi v v D k sR Dii i i i x i i i i i i i i i i B e y k s p x s ve v v D k sR k sR k sR ? ? ? + ? == + ? ?? + + ? ?? ?= = = ? ? ??? ? ? ?? ? ? ?+? ? ? ? ? ?= ? ? ? ? ? ?? ?+ + + ? + ? ?? ? ? ?? ?? ? ?? ? ? ? 1, 2,...i n ; ? = (C.12) Using a similar approach the second term ?[ ]{ }2A b ? is evaluated and the explicit expression for ? ( )2 ,ip x s ? is: 124 ( ) ( ){ } ( ) ( ) ( ) 1 1 2 2 2 1 2 2 21 1 2 2 2 1 1 1 2 2 2 2 3 3 3 1 3 2 1 1 42 2 2 2 , 1 , 1 4 ( ) 2 2 x i i x i i o i i i i i i x v v D k sR Di x x i x i i i i i i i x i i i i i i i i i i i R c y k D ve vp x s e v v D k sR D v k sR k sR k sR ? ? ? ? ? = + ? ?? + + ? ?? ? ? = = ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ?? ? ?+ ? ?? ? ?? ? = ? ?? ? ? ?? ? ?+ + +? ? ? ? ?? ? ? ? ? ? + ? ? ? + ? ? ?? ? ? ? ? 1 1 1, 2,... i i i n = ? ? ? ? ? ? ? ? ? ? ?? ? ; ? = ? (C.13) Substituting equations(C.12) and (C.13) into equation(C.10) we get the solution in the Laplace domain as: 125 ( ) ( ) ( ) ( ){ } ( ) ( ){ } ( ) ( ) 22 1 1 2 2 22 1 2 2 2 2 1 2 1 2 2 2 2 3 3 3 1 3 2 1 2 1 11 42 1 2 , , , , 1 1 4 ( ) 2 2 o i x i i x i i i t si ii i i i ii i ii x v v D k sR Dii i i i x i i i i i i i i i i p x s p x s p x s B e y k s ve v v D k sR k sR k sR ? ? ? + ? == + ? ?? + + ? ?? ?= = = ? = + = ? ??? ? ? ?? ? ? ?+? ? ? ? ? ? ? ? ? ? ? ?? ?+ + + ? + ? ?? ? ? ?? ?? ? ?? ? ? ? ( ){ } ( ) ( ) ( ) 1 1 2 2 2 1 2 2 21 1 2 2 2 1 1 1 2 2 2 2 3 3 3 1 3 2 1 1 42 2 2 , 1 1 4 ( ) 2 2 x i i x i i o i i i i i i x v v D k sR Di x x x i i i i i i i x i i i i i i i i i i i R c y k D ve v e v v D k sR D v k sR k sR k sR ? ? ? ? ? = + ? ?? + + ? ?? ? ? = = ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ?? ?+ ? ?? ? ?? ? + ? ?? ? ? ?? ? ?+ + +? ? ? ? ?? ? ? ? ? ? + ? ? ? + ? ? ?? ? ? ? ? ? 1 1 1, 2,... i i i n = ? ? ? ? ? ? ?? ; ? = ? (C.14) The final solution is obtained by taking an inverse Laplace transform of the solution given by equation(C.14). Inverse Laplace transform is performed as follows: ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 1 1 2 1 1 1 2 , , , , , , , 1, 2,... i i i i i i i c x t c x t c x t p x s p x s p x s p x s i n ? ? ? = + = + = + ; ? = ? ? ? (C.15) In Appendix D, the terms ? ( )1 1 ,ip x s?? ? and ? ( )1 2 ,ip x s?? ? are evaluated. Substituting equations (D.6) and (D.13) in equation(C.15) we obtain the final solution in the time domain as: 126 ( ) ( ){ } ( ){ } 1 2 2 1 2 1 32 1 1 1 2 2 1 2 12 1 1 1 1 1 1 1 2 1 11 2 2 2 1 1 1 2 1 1 , iii i i i i i i i ii i ii i o i i i i i i ii i c x t y k G h G G R c y k G h G G ? = = == + ? = == + ? ?? ?= + ? ?? ?? ? ? ?? ? ? ?? ? + + ? ?? ?? ? ? ?? ? ? ??? ? ?? 1, 2,...i n ; ? = (C.16) where; the ?G ? are defined in equations (D.6) and (D.13) in Appendix D. 127 APPENDIX D Factorization of the Laplace Terms for the Cauchy Boundary Condition The first term ? ( )1 1 ,ip x s?? ? can be evaluated as follows. From equations (C.14) and (C.15) we get: ( ) ( ) ( ){ } ( ) ( ){ } ( ) ( ) 22 1 1 2 2 22 1 2 2 2 2 1 2 1 2 2 2 2 3 3 3 1 3 2 1 1 1 1 11 1 42 1 2 , , , 1 1 4 ( ) 2 2 o i x i i x i i t si ii i i i ii i ii x v v D k sR Dii i i i x i i i i i i i i i i c x t p x s B e y k s ve v v D k sR k sR k sR ? ? ? ? + ? == + ? ? ? ? + +? ? ? ?= = = ? = ? ??? ? ? ?? ? ? ?+? ? ? ? ? ?= ? ? ? ? ? ?? ?+ + + ? + ? ?? ? ? ?? ?? ? ?? ? ? ? ? ? 1, 2,...i n ; ? = (D.1) Equation(D.1) can be rewritten as: 12 8 ( ) ( ){ } ( ){ } ( ) ( ) ( ) ( ) 2 2 2 1 2 2 23 3 1 1 2 1 3 2 2 3 2 4 2 4 2 4 4 1 4 2 4 1 4 24 2 4 2 1 1 41 1 2 1 2 , , , , , , , , 1 1 4 ( ) 2 2 x i io i x i i i i i i i i i x v v D k sR t s Di i ii i i ii i i x i i i i i i i i i i i i i R R i i i i R R y k c x t B e ve v v D k sR s k R s a ? ? ? = + ? ?? + + ? ?? ? + ? ? = = = ? = = ? ? ?? ? ?? ? ?? ? ? = ? ? ?? ? ? ?+ + + + ? ? +? ? ? ? ? ? ? ?? ? ?? ? ? ? 1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 , , , , , 1, 2,... ; i i i i i i i i i i i i i i i i i n kwhere k k k R R R and a R = ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ; ? = = ? , = ? = ? (D.2) Note that in equation(D.2) the condition ? * 2 , 0i ik ? ? must be satisfied for the solution to be determinate. Factorization of equation(D.2) gives: 12 9 ( ) ( ){ } ( ){ } ( ) ( ) ( ) ( ) ( ) 2 2 2 1 2 2 23 3 1 2 2 3 2 4 4 1 4 2 4 2 2 4 2 4 2 5 2 5 2 4 5 1 5 2 5 4 5 2 1 1 42 1 1 2 , , , , , , , , , , , 1 , 1 4 ( ) 2 2 1 x i io i x i i i i i i i i i x v v D k sR t s Di i i i x i i i i i i i i i R R i i i i i i i i i i i i i i i i i R R y k B e ve c x t v v D k sR s k R s a R a a ? ? ? = + ? ?? + + ? ?? + ? ? ? = ? = = ? ? ? ? ? ? ? ? ? ? = ? ?+ + + + ?? ? ? ? ? + ? ? ? ? ? ? ( ) 11 2 1 3 4 1 4 2 4 2 1 1 , , 1, 2,... i i i iii i i i i i i i i R R i n = = = = ? ? ? ? ? ? ? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ?? ? ; ? = ? ?? ? (D.3) It must be noted that the solution formulation as given by equation(D.3) is valid only when the condition ? 2 5 2 4, ,i i i i a a? ? is satisfied. Equation(D.3) can be further factorized and simplified as: 13 0 ( ) ( ){ } ( ){ } 1 2 2 1 2 1 32 1 2 3 3 1 1 1 1 1 1 1 1 2 1 11 42 2 1 1 1 , 1, 2,... ; 1 xo i x x iii i i i i i i i ii i x v D k t s Dxv i D i c x t y k G h G G i n where e ve B e G ? ? = = == + ? + ? + ? ? ?? ?= + ? ?? ?? ? ? ?? ? ; ? = ? = ? ??? ? ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ){ } 2 2 2 2 3 2 4 4 1 4 2 4 2 2 5 2 4 2 5 2 5 2 4 5 1 5 2 5 1 5 2 5 45 2 5 2 3 3 1 2 , , , , , , , , , , , , , 2 2 1 1 2 1 4 ( ) 2 2 1 i i i i i i i i o i x x sR x i i i i ii i ii i i i R R i i i i i i i i i i i i i i R R i i i i i i R R x t s Dxv i D i v v D k sR s s a k R R a a e ve B e G ? ? ? ?+ ? ?? ? = ? ? = ? = = ? ? ? ? ? + ? ? ?+ + + + +? ? ? ? ? ? ? ?? ? ? ?? ? ? ? ? = ? ? ? ? ( ) ( ) ( ) ( ) 2 2 2 2 2 3 2 4 4 1 4 2 4 2 4 2 , , , 1 4 ( ) 2 2 1 0 x i i i i v D k sR x i i i i i i i i i i R R v v D k sR s k if loopis not executedh M if loopis executed ? ? ?+ + ? ?? ? = ? = ? ?+ + + +? ? ? ? ? ? ?= ? ? ? ? (D.4) 13 1 The term ? 11G ? in equation(D.4) can be further factorized and simplified as: ( ){ } ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) 2 2 23 2 2 3 3 1 2 2 23 2 2 2 4 2 4 3 2 5 5 1 5 2 5 2 42 2 2 1 42 2 , 1 1 , , , , 1 1 4 ( ) 2 2 1 1 4 ( ) 2 2 x i io i x x x i io i x i i x v D k sR t s D xv x i i i i D i x v D k sRt s D x i i i i i i i i i i i i i i R R e ve v v D k sR s B e e e v v D k sR s a G a k ? ? ? ? ? ?? + + ? ?? + ? ? ? ? ?? + + ? ?? + ? ? = ? = ? ?? ? + + + +? ?? ? ? ? ?+ + + +? ? ? ?= ? ? ?? ? ( ) ( ) ( ) ( )4 1 4 2 4 2 2 4 2 5 2 5 2 4 5 1 5 2 5 4 5 2 , , , , , , , , , i i i i i ii i i i R R i i i i i i i i i i i i i i R R R R a a= ? ? = ? ? ? ? ? ? ? ? ?? ? ? ? ? ? (D.5) Again it must be noted that equation(D.5) is valid only when the condition ? 2 4 3, ,0i i i a a? ? where ? 3 3,0i i a ?= ? is satisfied. Using the results from Appendix E, inverse Laplace expressions for the terms ? 11G ? and ? 12G ? can be evaluated and the solution for ? ( )1 ,ic x t ? is obtained as: 13 2 ( ) ( ){ } [ ] ( ) ( ) ( ) 1 2 2 1 2 1 32 1 3 2 3 2 33 1 1 1 1 1 1 1 1 2 1 11 , ,0 , ,0 1 1 , 1, 2,... , ,oi iii i i i i i i i ii i t i i o i i oi i c x t y k G h G G i n F x t u t t e F x t t B G ? ? = = == + ? ? ?? ?= + ? ?? ?? ? ? ?? ? ; ? = ? ?? ? ?? ? = ? ??? [ ] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ] ( ) ( ) 3 2 2 4 2 2 4 4 1 4 2 4 2 2 4 3 2 5 2 4 2 5 2 5 2 4 5 1 5 2 5 1 5 2 5 45 2 5 2 33 1 2 3 , , , , , , , , , , , , , , , , , , ,0 1 2 , , , oi i i i i i i oi t i i i i o i i i o i ii i i i R R i i i i i i i i i i i i i i i i i R R i i i i i i R R ti i i i o F x t u t t e F x t t a k R R a a B F x t u t t e F G ? ? ? ? = ? ? = ? = = ? ? ? ? ? ?? + ? ?? ? ? ? ? ?? ? ? ? ?? ? ? ? ? ? = ? ? ? ( ) ( ) 2 3 2 4 4 1 4 2 4 2 , ,0 , , , , i i i i o i i i i i i i R R x t t k = ? = ? ??? ? ?? (D.6) where; the term ? 1 2 3, ,i i i F ? is given by equation (E.20) or (E.26) Appendix E. The second term ? ( )1 2 ,ip x s?? ? is evaluated as follows. From equations (C.14) and (C.15) we get: 13 3 ( ) ( ) ( ){ } ( ) ( ) ( ) 1 1 2 2 2 1 2 2 21 1 2 2 2 1 1 2 2 2 2 3 3 3 1 3 2 2 1 2 1 1 42 1 2 2 , , , 1 1 4 ( ) 2 2 x i i x i i i i o i i i i i i x v v D k sR i x D x x i i i i i i x i i i i i i i i i i i c x t p x s R c y k D ve v e v v D k sR D v k sR k sR k sR ? ? ? ? ? ? = + ? ?? + + ? ?? ? ? ? = = ? = ? ? ? ? ? ? ? ?? ? ? ?+? ? ? ?? ? = ?? ? ? ?? ?+ + +? ? ? ?? ?? ? + ? ? ? + ? ? ? ? ? ? 1 1 1 1, 2,... i i i i n = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ; ? = ? ? (D.7) Equation(D.7) can be written as: 13 4 ( ) ( ){ } ( ) ( ) ( ) 1 1 2 2 2 1 2 2 21 1 2 2 1 2 2 2 3 2 3 2 3 3 1 3 2 3 13 2 1 1 42 2 1 2 , , , , , , , 1 , 1 4 ( )2 2 x i i x i i i i o i i i i i i x v v D k sR Di x x i x i i i i i i i i i i i i i i i i R R i i R c y k D ve v ec x t v v D k sR s a R k R s a ? ? ? = + ? ?? + + ? ?? ? ?? ? = ? = = ? ? ? ? ? ? ? ?? ? ? ?+? ? ? ?? ? ?= ? ? ? ?? ?+ + +? ? ? ?? ?? ? ? ? ? ?? + ? ? +? ? ? ? ? ? ? ( ) 1 2 1 3 2 3 2 1 , 1, 2,... i i i i i ii i i i R R i n = = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ; ? = ? ? ? (D.8) Note that the condition ? 2 3, 0i ik ? ? must be satisfied in equation(D.8). Factorization of equation(D.8) yields: 135 ( ) ( ){ }1 1 2 2 1 2 12 1 2 21 2 2 2 2 1 1 1 2 1 1 42 2 1 2 1 , 1, 2,... ; 1 x ix x ii i o i i i i i i i ii i x v v D k Di x xv D c x t R c y k G h G G i n where D ve v e G ? ? = == + ? + + ? ? ?? ?= + ? ?? ?? ? ? ?? ? ; ? = ? ?+ ? ?? ? = ? ?? ? ( ){ } ( )( ) ( ) ( ) ( ) ( ) 2 1 2 2 1 2 2 3 3 1 3 2 3 2 2 2 4 2 3 2 4 2 4 2 3 4 1 4 2 4 1 4 2 4 34 2 4 2 2 , , , , , , , , , , , , , , 2 2 2 1 4 ( ) 2 2 i i i i i i i i x sR x x i i i i i i i i ii i i i R R i i i i i i i i i i i i i i i R R i i i i i i R R xv D ev v D k sR s a s a R k R R a a e G ? ? ? ? ?? ? ? ? = ? ? = ? = = ? ? ? ? ? ? ? ? ??? ? ? ?? ?+ + +? ? ? ?? ?? ? + + ? ? ? ?? ? ?? ? ? ? ? = ? ? ? ? ( ){ } ( ) ( ) 2 2 21 1 2 2 1 2 2 2 3 3 1 3 2 3 2 42 2 1 , , , , 1 1 4 ( ) 2 2 x i i x i i i x v v D k sR Di x x x i i i i i i i i i i i i R R D ve v e v v D k sR s a R k ? ? ? ?? + +? ?? ? ? ? ? = ? = ? ?? ? ? ?+? ? ? ?? ? ?? ? ? ?? ?+ + +? ? ? ?? ?? ? + ?? (D.9) 13 6 Note that the condition ? 2 5 2 4, ,i i i i a a? ? must be satisfied in equation(D.9). The term ? 21G ? in equation(D.9) can be further factorized and simplified as: ( ){ } ( ) ( ){ } 2 2 21 1 2 2 1 2 2 2 21 1 2 2 42 2 ,2 1 42 2 2 1 1 1 4 ( ) 2 2 1 1 4 ( ) 2 2 x i i x i x x i i x i x v v D k sR Di x x x i i xv i iD x v v D k sR Di x x x i i D ve v e v v D k sR s a e D ve v e v v D k sR G ? ? ? ? ? ?? + + ? ?? ? ? ?? ? ?? + + ? ?? ? ? ? ?? ? ? ?+? ? ? ?? ? ?? ? ? ?? ?+ + +? ? ? ?? ?? ? + ? ?? ? ? ?+? ? ? ?? ? ?? ? ? ?? + + +? ? ? ? ?? ? = ? ( ) ( ) ( ) ( ) ( ) ( ) 2 3 3 1 3 2 3 2 2 3 1 2 2 2 4 2 3 2 4 2 4 2 3 4 1 4 2 4 1 4 2 4 34 2 4 2 , , , , , , , , , , , , , , , i i i i i i i i i i ii i i i R R i i i i i i i i i i i i i i i i i i i R R i i i i i i R R s a a a R k R R a a= ? = ? = ? = = ? ? ? ? ?? + ? ? ? ?? ? ? ?? ? ? ? ? ? ? (D.10) Note that the condition ? 2 4 1 2, ,i i i i a a ?? ? must be satisfied in equation(D.10). ? 21G ? can be further simplified as: 13 7 ( ){ } ( ) ( ) ( ){ } ( ) 2 2 21 1 1 2 1 2 2 2 2 2 21 2 3 2 2 42 1 1 2 , , 2 1 42 1 2 , 2 1 1 1 1 4 ( ) 2 2 1 1 4 ( ) 2 2 x i i x i x x i i x i x v v D k sR Di x x i ii i x i i xv D x v v D k sR Di x i i x i i D ve v e v s as a v D k sR e D ve v e vs a v D k sR G ? ? ? ? ? ?? + + ? ?? ? ?? ? ?? ? ? ?? + + ? ?? ? ?? ? ?+ ? ?? ? ?? ? + + + + +? ?? ? ? ?+ ? ?? ? ? +? ? + + + +? ?? ? = ? ? ? ? ( ) ( ) ( ) ( ) ( ) ( ) 1 2 3 3 1 3 2 3 2 2 3 1 2 2 2 4 2 3 2 4 2 4 2 3 4 1 4 2 4 1 4 2 4 34 2 4 2 1 , , , , , , , , , , , , , , , 1 i i i i i i x i i i i ii i i i R R i i i i i i i i i i i i i i i i i i i R R i i i i i i R R s a a a R k R R a a ? = ? ? ? = ? = = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?+ ? ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ? ? ? ? ? (D.11) ? 22G ? in equation(D.9) can be further simplified as: 13 8 ( ){ } ( ) ( ) ( ) 2 2 21 1 1 2 1 2 2 2 2 2 3 3 1 3 2 3 2 42 2 1 1 2 , , 2 2 , , , 1 1 1 4 ( ) 2 2 x i i x ix i i x v v D k sR Di x xv xD i ii i x i i i i i i i i i i R R D ve ve e v s as a v D k sR G R k ? ? ? ?? + +? ?? ? ?? ? ?? = ? = ? ?? ? ? ?+? ? ? ?? ?? ? ? ?? ? + + + + +? ?? ?? ? ? ?? ? = ?? ? ? (D.12) Using equation(A.15) and Appendix E, inverse Laplace expressions for the terms ? 21G ? and ? 22G ? can be evaluated and the solution for ? ( )2 ,ic x t ? is obtained as: 13 9 ( ) ( ){ } [ ] 1 1 2 2 1 2 12 1 ,1 11 2 1 2 2 2 2 2 1 1 1 2 1 1 , , 2 1 , 1, 2,... ; 1 , i i i ii i o i i i i i i i ii i x ai x i i i c x t R c y k G h G G i n where D F x t e vG ?? ? ? = == + ? ? ? ? ?? ?= + ? ?? ?? ? ? ?? ? ; ? = ? ?+ ? ? ?? ? = ? ?? ( ) [ ] ( ) ( ) ( ) ( ) ( ) ( ) ,1 2 32 1 2 2 3 3 1 3 2 3 2 2 3 1 2 2 2 4 2 3 2 4 2 4 2 3 4 1 4 2 4 1 4 2 4 34 2 4 2 1 2 1 2 , , , , , , , , , , , , , , , , , , 2 2 1 , 1 i i i i i i i i i x a tt i x i i ii i ii i i i R R i i i i i i i i i i i i i i i i i i i R R i i i i i i R R i x i i i D F x t e v a a R k R R a a D F vG ?? ? ? ? = ? ? ? = ? = = ? ? ? ? ? ?? + + ? ?? ? ? ? ? ?? ? ? ?? ? ? ? ? ?? + ? ?? ? = ? ? ? [ ] ( ) ( ) ,1 1 2 2 2 3 3 1 3 2 3 2 , , , , i i i i i x a t i i i i i i i i R R x t e R k ? ?? ? = ? = ? ?? (D.13) where; the term ? 1 2 3, ,i i i F ? is given by equation (E.20) or (E.26) Appendix E. 140 APPENDIX E Evaluation of Inverse Laplace Expressions for the Cauchy Boundary Condition ( ){ } ( ) ( ) ( ) 2 1 14 1 2 3 4 1 1 2 3 42 2 1 , , , 2 , 1 1 4 2 2 x i io i x x x v D k sR t s Dxv D i i i i x i i i i e ve e v v D k R s s a ? ? ? ?? + + ? ?? + ? ? ? ? = ? ? + + + +? ?? ? ? (E.1) Equation(E.1) can be simplified as: 14 1 ( )1 2 3 4 2 1 1 1 1 1 1 1 1 2 3 1 1 1 1 1 11 2 , , , 1 2 4 1 1 2 2 2 , 1 2 ; 4 4 42 x i i x i x i o xv D i i i i R kvx s D R D R i i i i x i i i x i i x i i x ii x t e where ve k k kv v v vR D s s a R D R R D R R D RR D e ? ? ? ? ? ? ? ?? ?? ? ?? ?? ? ?? ? ? ?? ?? ? ? ? ? = ? = ? ? ? ?? ? ? ?? ?? ?? ? + ? ? ? ? ? ? +? ? ? ?? ?? ?? ? ? ? ? ?? ?? ? ? ?? ?? ? = ? ? 2 1 1 1 14 1 1 1 1 2 3 1 1 1 1 1 11 4 2 2 2 ,4 4 42 i i x i x ii o R kvx s D R D Rt s i i i i x i i i x i i x i i x ii x e ve k k kv v v vR D s s a R D R R D R R D RR D ? ?? ?? ? ?? ?? ? ?? ? ? ?? ?? ? ? ? ?? ?? ? ? ?? ?? ?? ? + ? ? ? ? ? ? +? ? ? ?? ?? ?? ? ? ? ? ?? ?? ? ? ?? ?? ? (E.2) ? 1? ? can be evaluated as follows: Invoking the ?First Shifting Theorem? (see [28], p253) we can simplify ? 1? ? as: 14 2 2 1 1 1 1 1 1 2 3 1 11 4 1 1 2 ,42 ii xi x i Rkv x st DR D R i x i i i i x ii x e ve R D kv v s s aR D RR D ? ? ?? ?? ? ???? ? ? ?? ?? ? ? ? ?= ? ?? ?? ?? ?? ? + ? + ?? ?? ?? ?? ? ? ?? ? ? ?? ?? ? ? (E.3) The Laplace inverse of equation(E.3) can be readily obtained from the tables (see [8], p496, eq31). Applying this inversion, ? 1? ? can be evaluated as: 14 3 2 1 1 , 2 3 1 1 1 1 2 3 1 11 2 3 ,2 3 1 1 1 2 1 11 , 2 3 1 1 1 4 2 ,2 , 1 4 2 2 4 42 2 42 i i i i x i x i i i i i i i x i x i R kvx a D R D R i i i i x i x ii i ia t i x ii x R kvi x x a D R D R i x R ke x verfc a t D t R D Rkv v aR D RR D ve R D e v v R D ? ? ?? ? + ? ? ? ?? ? ? ? ?? + ? ? ? ?? ? ? ?? ?? ? ? + ? +? ?? ?? ?? ? ? ?? ? ? ?? ?+ + ?? ? ? ?? ?= ? 1 1 2 3 1 11 2 3 1 1 2 1 1 1 1 1 1 2 3 1 1 1 2 , , 4 2 2 , 2 4 2 4 4 i i x i i i i i x i x ii i i i x i kv t R D R i x i i x i i i x i x i R kx verfc a t D t R D Rk aR D R vve R D kv vR D a R D R D R ? ?? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ?? ?? ?+ + ? ? ?? ?? ?? ?? ? ? ?? ?? ? ? ?? ?+ ?? ? ? ?? ? ? ?? ?? ? ? ? ? ? ? ? + ?? ?? ? ?? ?? 12 1 11 1 42 2 2 i x i xi x iRv v t x D R DR D x i x Rx D v te erfc t R D ? ? ? ?+ ? ?? ? ? ? ? ?? ? +? ?? ? ? ?? ? ?? ?? ?? (E.4) Note: the solution given by equation(E.5) is valid only when ? 1 2 3 1 , i i i i k a R ? ?. ? 1? ? can be further rearranged and simplified as: 14 4 12 1 ,1 2 3 1 1 2 3 1 1 11 1 2 3 1, 2 3 2 1 ,1 2 3 1 1 1 24 ,2 2 , 1 42 2 4 2 4 4 i i x i i x i i i i i x i i x i k ix v R D a i i x i iD R i i xi i x i i ia t kx v R D a D R i i x i kR x v R D a t Rve erfc R D t kv v R D a R e ve kv v R D R ? ? ?? + ? ? ? ? ?? ? ? ? ?+ ? ? ? ? ?? ? ? ?? ? ? ?? + ? ? ?? ? ? ?? ?? ? + ? ?? ? ? ?+ + ?? ? ? ?? ?? ? ? ?? ?= ? + 1 1 1 2 3 1 1 2 3 1 1 11 1 1 2 3 1 2 , , 22 , 4 2 22 i x i i i i x i i i i x i i k txv D R i i i x i i i kR x v R D a t R erfc R D t a R x vtve erfc RkR D a R ? ?? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ?? ?? ? + + ? ? ?? ?? ?? ? ? ?? ?? ?? ? ? ?? ? ? ? ? ??? ? ? ? ? ?? ?? ? ? ?? ? ? ?? ? ? ? ? ? + ?? ? ? ? ? ? ?? ?? ? ?? ?? ? ? ?? ? 1i x D t ? ?? ? ? ? ? ?? ? (E.5) ? 2? ? is evaluated as follows: ? 2? ? can be simplified as; ( )4 12 1i o ot t se e?? ?? ??= ? ? (E.6) 14 5 Now, we invoke the ?Second Shifting Theorem? (see [28], p265) and evaluate ? 2? ? as: [Note: We have already evaluated ? 1? ?; equation(E.5)] ( ) ( ) ( ) ( ) 4 12 1 ,1 2 3 1 1 2 3 1 1 11 1 2 3 1, 2 3 2 1 ,1 2 3 1 2 24 ,2 2 , 42 4 24 i o i i x i i x i i i o i i x i i x i t o k ix v R D a i i x i i oD R i i x oi i x i i ia t t kx v R D a D R e u t t kR x v R D a t t Rve erfc R D t tkv v R D a R e ve ?? ? ? ?? + ? ? ? ? ?? ? ? ? ?+ ? ? ? = ? ? ?? ? ? ?? + ? ?? ?? ? ? ?? ?? ? + ? ??? ? ? ?+ + ?? ? ? ?? ?? ? ? ?? ? ( ) ( ) ( ) 1 1 1 2 3 1 11 1 2 3 1 1 1 1 1 2 3 1 2 , 2 , 22 , 4 24 2 i o x i i i i x i i o i i x oi i x i i i k t txv D R i i x i i i kR x v R D a t t R erfc R D t tkv v R D a R ve kR D a R ?? ? ?? ? ??? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ?? ?? ? + + ? ?? ?? ?? ?? ? ? ?? ?? ?? ? ? ??? ? ? ? ? ?? + ?? ? ? ? ? ?? ?? ? ? ?? ? ? ?? ? ? ? ? ?? ?? ? ( ) ( ) 1 1 2 i o i x o R x v t terfc R D t t ? ? ? ? ? ?+ ? ? ? ? ?? ? ? ? ?? ? ? ?? ?? ? ???? ? (E.7) 146 Now, for the case when ? 1 2 3 1 , i i i i k a R = ?, we evaluate ? 1? ? as follows: From equation(E.3) we get: 2 1 1 1 1 1 11 4 1 1 2 42 ii xi x i Rkv x st DR D R i x i xi x e ve R D v v s s R DR D ? ? ?? ?? ? ???? ? ? ?? ?? ? ? ? ?= ? ?? ?? ?? ? + ?? ?? ? ? ?? ?? ?? ? ? (E.8) ? 1? ? can be factorized as: 2 1 11 1 1 1 1 1 4 1 1 2 3 21 2 2 i ii x i x kv t RR D R x s Di x i x i x i x vs R D e R De v v v s sR D R D ? ? ?? ? ? ? ? ? ? ?? ? ? ?? ? ?? ? ? ? ? ? ?+ ? ?? ? = ?? ?? ? ? ?? ?? ?? ? ?+? ?? ? ? ? ? ?? ?? ?? ? ? (E.9) ? 1? ? can be further simplified as: 1 12 1 1 1 1 1 4 1 1 22 4 2 i i i x x i x i R Rx s x s kv D Dt R D R i x i x e ee v vs sR D R D ? ? ? ? ? ? ? ? ?? ?? ?? ? ? ? ??? ? ? ? ? ? ? ? ?? ?= ? ? ? ? ?? ? ? ? ?+? ? ? ?? ? ? ? ? (E.10) Equation(E.10) can be conveniently rewritten as: 147 ( ) 2 1 1 1 1 1 1 1 4 1 1 2 1 1 1 2 22 ; 4 2 i i x i i i x x kv t R D R R Rx s x s D D i x i x e where e eand v vs sR D R D ? ? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? = ? = =? ? ? ??? ? ? ?+? ? ? ?? ? ? ? ? ? (E.11) Laplace inverse of ? 1? ? can be obtained directly by using the tables (see [8], p495, eq19). Applying this inversion, ? 1? ? can be evaluated as: 1 1 1 2 1 1 1 1 1 1 2 4 1 2 2 2 2 2 2 i x i x i x i x i x iR v x D R D x i xv t R D iR v x D R D x i x Rx D v te erfc t R D e Rx D v te erfc t R D ? ? ? ?? ? ? ?? ?? ? ? ??? ? ? ?? ? ? ?? ? ? ?? ?= ? ?? ? ? ?? ? ? ?? ? ? ?+ +? ? ? ?? ? ? ?? ? ? ?? ?? ? (E.12) ? 1? ? can be further simplified as: 2 1 1 1 1 1 4 2 2 1 2 2 2 i x x x v t xv xvR D i iD D x i x i R x vt R x vte e erfc e erfc D R t D R t? ?? ?? ? ? ?? +? ? ? ? ? ?= +? ? ? ? ? ?? ? ? ?? ? ? ?? ? (E.13) Laplace inverse of ? 2? ? can also be obtained directly by using the tables (see [8], p495, eq17). Applying this inversion, ? 2? ? can be evaluated as: 148 2 1 1 12 1 111 11 1 4 2 2 42 2 2 1 2 42 2 2 i x i x i xi x x R D t i x iRv tv x D R DR D i x x i xi x i x v t e R D Rx R Dv v t v tx e erfc D R DR D t R D ? ? ? ? ? ? ?+ ? ?? ? = ? ? ? ? ?? ? ? ? ? ?+ + + +? ?? ? ? ?? ? ? ?? ? (E.14) ? 2? ? can be further simplified as: 22 1 11 1 1 1 2 2 4 4 2 1 2 2 2 i x i xx xv v tR x D R D iD t i x x i x x i R x vtt xv v tv e e erfc R D D R D D R t? ? ? ?? + ? ? ? ?? ? ? ?? ? +? ?= ? + + +? ? ? ? ? ? ? ?? ? ? ? (E.15) From equations(E.11), (E.13) and (E.15) we can evaluate ? 1? ? as: 1 1 1 1 1 1 2 2 1 1 1 1 1 1 2 2 1 24 4 2 1 1 2 22 2 1 2 2 2 x x i i i x i x x xv xv i iD D k t x i x iR R x v t xv D t R D iD i x x i x x i R x vt R x vte erfc e erfc D R t D R t e R x vtt xv v tv e e erfc R D D R D D R t ? ? ? ? ? ?? ? ?? ? ?? ? ? ?? ? ? ?? ?? ? ? ? ? ?+? ? ? ? ? ?? ? ? ?? ? ? ? ? ?= ? ?? ?? ? +? ? ? ?+ ? + +? ? ? ? ? ?? ?? ?? ? ? ?? ? (E.16) ? 1? ? can be further simplified as: ( )21 1 1 1 111 1 1 1 2 4 2 1 2 1 2 2 1 1 2 2 i x i i i x x R x vt D R ti k t xv i xx iR D xv iD x i x x i R x vt v terfc e R DD R t e e R x vtxv v t e erfc D R D D R t ?? ? ? ? ? ? ?? ??? ? ? ?+? ? ? ?? ?? ? ? ?= ? ?? ?? ? +? ? ? ?? + +? ? ? ? ? ?? ? ? ?? ? ? ?? ? (E.17) 149 In order to evaluate ? 2? ? for the case when ? 1 2 3 1 , i i i i k a R = ?, we use the same method as before. From equation(E.2) ? 2? ? can be written as: ( )4 12 1i o ot t se e?? ?? ??= ? ? (E.18) Now, we invoke the ?Second Shifting Theorem? (see [28], p265) and evaluate ? 2? ? as: [Note: We have already evaluated ? 1? ?; equation(E.16)] ( ) ( ) ( ) ( ) ( ) ( )( )( ) ( ) ( ) ( ) 2 1 1 1 1 114 1 1 1 1 2 4 2 2 2 1 2 2 1 1 2 2 i o x i o i o i o i x x R x v t t D R t ti o o k t t xv i xx i ot R D o xv i oDo x i x x i o R x v t t v t terfc e R DD R t t e u t t e e R x v t tv t txv e erfc D R D D R t t ? ?? ? ? ? ? ? ? ? ? ? ?? ?? ? ?? ? ? ?+? ? ? ??? ? ? ?? ?= ? ? ?? ?? ? + ?? ? ?? ?? + +? ? ? ? ? ?? ? ?? ?? ? ? ?? ?? ? (E.19) The final solution can be compactly expressed as: 150 [ ] [ ] ( ) ( ) ( ) [ ] ( ) ( ) 3 1 2 3 4 1 2 3 1 2 3 , ,1 2 3 1 1 2 3 11 2 3,2 3 1 2 3 , ,1 2 3 1 1 1 2 3 , , , , , , , 2 , , , ,2 , , 2 , , , , , , ; 2 , i o i i i x i i x i i i x t i i i i i i i o i i i o x D i i i i xv i xi i i a t D i i i x D i i i i i i x t F x t u t t e F x t t where R x te erfc R D tv F x t ve e R xe erfc v ? ? ? ? ? ? ? ? ? ? ? ? ?= ? ? ?? ? ? ?? ? ? +? ?+ ? ?? ?= + ? ( ) 2 3 1 1 1 11 2 3 111 2 3 1 1 11 , 2 ,2 2 , , 2 2 2 , 2 1 2 2 i x i i i i i x k txv D R i i i i ii xi i i i k t i xR t R D t R x vt kv e erfc when a and RR D tv R x vt verfc R D te ? ? ?? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ?? ? ? ?? ?? ?? ? ? ?+? ? + ; ? ? ?? ? ?? ? ? ?? ? ?+? ? ? ? = ( )21 1 1 1 2 3 1 1 1 1 1 1 2 3 1 2 3 1 4 , 2 2 , , , 1 1 2 2 4 i i x x R x vt R D t i x i i i xv i iD x i x i x i i i i i x i i i t e R D k when aR R x vtxv v t e erfc D R D R D t kwhere v R D a R ? ? ? ?? ? ? ? ? ? ? ? ; = ? ?? ?? ? + ? ?? ?? + +? ? ? ?? ? ? ?? ? ? ?? ? ? ? ; = + ? ? ?? ? ? ? (E.20) The above solution is valid only for real values of ? 1 2 3, ,i i i ? ?. For problems involving complex values for ? 1 2 3, ,i i i ? ? in the case when ? 1 2 3 1 , i i i i k a R ? ? the ? 1 2 3, ,i i iF ? terms are evaluated as follows: (Note: when ? 1 2 3 1 , i i i i k a R = ? the ? 1 2 3, ,i i iF ? terms are unchanged) 151 1 1 2 3 1 2 3 1 2 3 1 1 1 2 3 1 2 3 1 2 * , , , , , * 2 , , , ; 4 ; 4 i i i i i x i i i i i i i i i i i x i i i kLet v R D a i R kwhere v R D a R ? ? ? ? ? = + ? = ? ?? ? ? ? ? ? = + ? ? ?? ? ? ? (E.21) Now the ? 1 2 3, ,i i i F ? terms can be written as: [ ] ( ) ( ) * , ,1 2 3 1 1 2 3 11 2 3,2 3 *1 2 3 , ,1 2 3 1 1 2 3 11 2 3 1 *2 , , * , ,2 , , *2 , , * , , 2 2 , 2 2 i i i x i i x i i i x xi D i i i i xv i xi i i a t D i i i xi D i i i i i xi i i i R x i te erfc R D tv i F x t ve e R x i te erfc R D tv i v R D ? ? ? ? ? ? ? ? ? ? ?? ? ? ?? ? ?+ ? ?? ?+ ? ?? ?? ?= ? ? ? ?? ?+ ? ? ? ?? ? ? ?? ? ?? ?? ? + 1 11 11 2 3 1 , 2 i x i k txv D R i i xi x i i i R x vte erfc R D tk a R ? ?? ? ? ? ?? ? ? ?+? ?? ? ? ? ? ?? ?? ? ?? ? ? ? (E.22) Using the symmetric relations: { } { }; ;, , ,if erfc a ib A iB then erfc a ib A iBwhere a b A B R + = + ; ? = ? ; ? (E.23) Also the exponent of a complex number can be expressed as: ( ) ( )( ) ( ) i i e Cos i Sin e Cos i Sin ? ? ? ? ? ? ? = ? = + (E.24) Using equations(E.23) and (E.24) equation(E.22) can be simplified as: 15 2 [ ] ( ) ( ) ( ) ( ) 1 2 3 1 2 3 1 2 3 ,2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 * * , , , ,* , , 2 , , 2 * 2 * * , , , , , ,* , , 2 2 , 2 2 i i x i i i i i i i i ixv x xa t D i i i i i i i i i i i i i i i x x x xv i Cos i Sin A iB D Dv F x t e ev x xv i Cos i Sin D D ? ?? ? ? ? ? ? ? ?? ? ? ?? ? ? ? ? +? ? ? ?? ?? ? ? ? ? ?? ? ? ?? ?= + ? ?? ? ? ?? ? + + ? ? ? ?? ?? ? ? ? ? ?? ? ? ?? ? ( ) ( ) ( ) 1 11 11 1 2 3 1 1 1 2 3 1 1 2 3 1 1 2 , * * , , , , 22 ; 2 2 i x i k txv D R i i xi i x i i i i i i i i i i i i x i x A iB R x vtv e erfc R D tkR D a R where R x i t R x i tA iB erfc and A iB erfc R D t R D ? ? ? ?? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ?+ ? ? ? ?? ? ? ?+? ? + ? ?? ? ? ?? ??? ? ? ?? ? ? ?? + ? ? ? = ; + =? ? ? ?? ? t ? ?? ? ? ? ? ?? ? (E.25) 153 Further simplification of equation(E.25) yields: [ ] ( ) ( ) ( ) 1 2 3 1 2 3 ,2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 1 1 1 2 3 1 * , ,* , , 2 , , 2 * 2 * , , , ,* , , 2 , 22 , 2 2 i i x i x i i i i i i ixv xa t D i i i i i i i i i i i i x k txv D R i i i x i i i xAv B Cos Dv F x t e ev xA Bv Sin D Rv e erfc kR D a R ?? ? ?? ? ? ?? ? ? ? ?? ? ? ?? ?? ? ?? ?? ? ? ?? ?= ? ?+ ? ? ? ?? + ? ? ? ?? ?? ?? ? + ? ? ?? ?? ? ? ? ( ) 1 1 1 1 2 3 1 2 3 1 1 1 2 3 1 * 2 , , , * , , 2 4 2 i x i i i i i x i i i i i i i i x x vt R D t kwhere v R D a R R x i tand A iB erfc R D t ? ? ? ?+? ? ? ? ? ?? ? ? ? ; = + ? ? ?? ? ? ? ? ?+ ? ? + = ? ? ? ?? ? (E.26) Hence in problems where ? 1 2 3, ,i i i ? ? is complex; for the case when ? 1 2 3 1 , i i i i k a R ? ? the ? 1 2 3, ,i i i F ? term is given as by equation(E.26). 154 APPENDIX F Derivation of the Steady State Solutions Dirichlet Boundary The system of governing transport equations can be written in a matrix format as [11, 38]: { } { } [ ]{ } 2 2x d c d cv D K c dx dx? = (F.1) where; ?[ ]? denotes a square matrix and ?{ }? denotes a column vector. The corresponding boundary conditions can be written as: ( ) 0dcdx? ?? =? ? ? ? (F.2) ( ){ } { } 1 1 1 0 ; 1, 2,... i i i i i c where B i n = = ? ? = ; ? = ? (F.3) Now in order to uncouple the system of ordinary differential equations (ODEs) given by equation(F.1), we apply a linear transform procedure described by Clement [11]. i.e. we perform the following matrix operation. { } [ ]{ }c A b= (F.4) 155 where; ?{ }b ? is the concentration in the transformed domain. Applying this transformation equation(F.4) gets modified as: [ ]{ } [ ]{ } [ ]( )[ ]{ } 2 2 1 0 x x d A b d A bv K A b dx D dx D ? ?? + = ? ?? ? (F.5) Pre-multiplying equation(F.5) with ?[ ] 1A ? ? we get: { } { } { } [ ] [ ][ ] 2 ~ 2 ~ 1 1 0 ; x x d b d bv K b dx D dx D where K A K A? ? ? ? ?? + = ? ? ? ?? ?? ? ? ? = ? ?? ? (F.6) By using the similarity transformation procedure described in Clement [11], we can make ? ~ K? ?? ?? ? ? a diagonal matrix and thus uncouple the system of equations. The corresponding ?[ ]A ? matrix is given as: [ ] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 1 1 2 3 31 1 11 2 1 1 1 1 ,0,0,..., 1 , ,0,0,..., 2 . . , ,..., ,0 1,1,..., n i i i i n n i i i ii i i i n n n i i n i i n i n i ni i i i i i k k n times y k k k k k n times y k y k A k k k k k k y k y k y k ntimes = ? = =? ? ? = = =? ? ? ? ?? ? ? ? ? ? ? ? ?? ? ? ? ?? ? ? ? ? ?= ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (F.7) The ?[ ] 1A ? ? matrix is: 156 [ ] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 1 1, 1 1 1 2 3 2 21 1, 2 2, 2 1 1 1 2 3 1, 2, 1, ,0,0,..., 1 , ,0,0,..., 2 . . , ,..., n i i i n i i i n n i i i i i i n n i i i i i i n n n i i i i i i i i i n n n n i n i n i i i n i i n i n y k n times k k y k y k n times k k k k A y k y k y k k k k k k k ? = = ? ? ? = = ? = ? = ? ? ? ? = = = = ? = ? = ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ( ) ,1n i n? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? (F.8) The corresponding ? ~ K? ?? ?? ? ? matrix is: ( ) ( ) ( ) ( ) ( ) 1 2 ~ ,0,0,... 1 0, ,0,0,... 2 . 0,0,... 1 , ,0,0,... . 0,0,... 1 , i n k n times k n times K i times k n i times n times k ? ?? ? ? ?? ? ? ? ? ?? ? ? ?=? ? ? ?? ? ?? ? ? ? ? ? ? ?? ?? ? (F.9) Equation(F.5) describes a set of ?n ? independent second-order homogeneous ODEs the boundary conditions of which are obtained by linear transforming equations(F.2) and(F.3). This yields: 157 ( ) 0dbdx? ?? =? ? ? ? (F.10) ( ) ( ) ( ) 2 2 2 1 2 2 1 2 1 1 1 2 1 1 2 , 1 1 0 1, 2,... n i i i i n i i i i i i ii i i i i i y k k k b B i n ? = + = ? = =? ? ? ? ? ? ?= ? ?? ? ?? ? ; ? = ? ? ? ? (F.11) Since equation(F.5) is uncoupled, it can now be written as a set of ?n ? independent equations as: ( ) ( ) ( ) ( )2 2 1 0 1, 2,... i i i i x x d b x db xv k b x dx D dx D i n ? ?? + ? = ? ?? ? ; ? = (F.12) The general solution to equation(F.12) is given as: ( ) 2 2 2 2 4 4 2 21 2 1, 2,... i i x x x xx x k kx v v x v v D D D DD D i i ib x e e i n ? ? ? ?? ? ? ?? ? ? ? ? ? ? ?+ + ? +? ? ? ? ? ? ? ?? ? ? ?? ? ? ?? ? ? ?= ? + ? ; ? = (F.13) where; 1i? and 2i? are constants. In order to apply the boundary condition given by equation(F.10) we differentiate the general solution with respect to ? x ?. Differentiation of equation(F.13) with respect to x yields: ( ) 22 2 2 42 21 2 42 22 2 1 4 2 1 4 2 i x xx i x xx kx v v D DDi i i x x x kx v v D DDi i x x x db x v v k e dx D D D v v k e D D D ? ?? ?? ? ? ?+ +? ? ? ?? ?? ?? ? ? ?? ?? ? ? ?? +? ? ? ?? ?? ?? ? ? ?? ?? ? ? ?= ? + +? ? ? ?? ?? ?? ? ? ?? ?? ? ? ? + ? ? +? ? ? ?? ?? ?? ? 1, 2,...i n ; ? = (F.14) 158 To satisfy the boundary condition given by equation(F.10), i.e. when ? x ? tends to ???; the exponential function in the first term tends to ???, hence 1i? must vanish. Equation(F.13) now reduces to: ( ) 2 2 4 22 1, 2,... i x xx kx v v D DD i ib x e i n ? ?? ?? ? ? ?? +? ? ? ?? ?? ?? ?= ? ; ? = (F.15) Applying the second boundary condition given by equation(F.11), we get: ( ) ( ) 2 2 2 1 2 2 1 2 1 1 1 2 1 1 2 , 2 1 1 1, 2,... n i i i i n i i i i i i ii i i i i i y k k k B i n ? = + = ? = =? ? ? ? ? ? ?? = ; ? = ? ?? ? ?? ? ? ? ? ? (F.16) Note that the condition ? 1 2, 0i ik ? ? where 1 21 ,i i n? ? and 1 2i i? must be satisfied for equation(F.16) to be valid. Therefore, the solution in the ?b?domain is: ( ) ( ) ( ) 2 2 2 1 2 2 1 2 2 1 2 2 1 1 2 1 1 , 4 2 1 1 1, 2,... i x xx n i i i i n i i i i i i kx v vi i D DD i i i i i y k k k b x B e i n ? = + = ? ? ?? ?? ? ? ?? +? ? ? ?? ?? ?? ? = =? ? ? ? ? ? ?= ? ?? ? ?? ? ; ? = ? ? ? ? (F.17) Inverse linear transform of equation(F.17) is done to obtain the solution in the ?c ? domain by using equation(F.4). The final steady state solution for the Dirichlet boundary is: ( ) { } ( ) ( ) 2 2 1 2 2 2 1 1 2 2 12 1 2 3 3 1 3 2 42 1 1 11 , 1, 2,... x i x x v v D k Diii i i i i i i i i i i ii i i i i i i i ec x y k B k k i n ? ?? + ? ?? ? ? = = == + = ? ? ? ? ?? ?? ? ? ?= ? ?? ? ? ?? ?? ? ? ? ? ?? ? ; ? = ? ? ?? ? (F.18) 159 Cauchy Boundary For the case of the Cauchy boundary, the boundary condition at the source is described as follows: ( ) ( ) 1 1 1 0 0 1, 2,...ii i x i i i dcD vc B v i n dx =? + = ; ? = ? (F.19) Since the governing equations and the boundary condition at ?? ? are identical for the Dirichlet and the Cauchy boundaries, the procedures involved in uncoupling the system of equations will be identical and hence we can use the previous section to obtain the semi-determined general solution and the linear transform matrices ?[ ]A ? and ?[ ] 1A ? ?. The semi-determined general solution [identical to equation(F.15)] is given by: ( ) 2 2 4 22 1, 2,... i x xx kx v v D DD i ib x e i n ? ?? ?? ? ? ?? +? ? ? ?? ?? ?? ?= ? ; ? = (F.20) The linear transform matrices ?[ ]A ? and ?[ ] 1A ? ? are given by equations(F.7) and (F.8) respectively. The constant ? 2i? ? in equation(F.20) is evaluated by using the boundary condition given by equation(F.19) after transforming them into the ?b? domain by applying the linear transform given by equation(F.4). ( ) ( ) ( ) ( ) 2 2 2 1 2 2 1 2 1 1 1 2 1 1 2 , 1 1 0 0 1, 2,... n i i i i n i i i i i i ii i i x i i i i y k k k vdb D vb Bdx i n ? = + = ? = =? ? ? ? ? ? ?? + = ? ?? ? ?? ? ; ? = ? ? ? ? (F.21) 160 Note that the condition ? 1 2, 0i ik ? ? where 1 21 ,i i n? ? and 1 2i i? must be satisfied for equation(F.21) to be valid. Substituting the expression for ? ( ),ib x s ? from the semi-determined general solution given by equation(F.20) into the transformed boundary condition given by equation(F.21), we can evaluate the constant ? 2i? ? as follows: ( ) ( ) 2 2 2 1 2 2 1 2 1 1 1 2 1 1 2 , 2 2 2 2 1 1 41 2 1, 2,... n i i i i n i i i i i i ii ii x i i i i ix x x y k k k vkv v D v BD D D i n ? = + = ? = =? ? ? ? ?? ?? ?? ? ? ?? ?? ? ? + + ? =? ? ? ?? ?? ?? ? ?? ? ? ?? ? ; ? = ? ? ? ? (F.22) From equation(F.22) ? 2i? ? is evaluated as: ( ) ( ) 2 2 2 1 2 2 1 2 1 1 1 2 1 1 2 , 1 1 2 2 1, 2,...1 42 2 n i i i i n i i i i i i ii i i i i i x i y k k k v B i nv v D k ? = + = ? = =? ? ? ? ? ? ? ? ?? ? ?? ? ? = ; ? = ? ? + +? ?? ? ? ? ? ? (F.23) Therefore, the solution in the ?b?domain is: ( ) ( ) ( ) 2 2 2 1 2 2 1 2 2 2 1 2 1 1 2 1 1 , 4 2 21 1 1 4 2 2 1, 2,. i x xx n i i i i n i i i i i i kx v v D DDii i i i i i x i y k k k veb x B v v D k i ? = + = ? ? ?? ?? ? ? ?? +? ? ? ?? ?? ?? ? = =? ? ? ? ? ? ?= ? ? ? ?? + +? ? ? ? ? ?? ? ; ? = ? ? ? ? ..n (F.24) Inverse Linear transform of equation(F.24) is done to obtain the solution in the ?c ? domain by using equation(F.4). The final steady state solution for the Dirichlet boundary is: 161 ( ) { } ( ) ( ) 2 2 1 2 2 2 1 1 2 2 12 1 2 2 3 3 1 3 2 42 1 21 11 , 1 4 2 2 x i x x v v D k Diii i i i i i i i i i i ii i x i i i i i i i vec x y k B v v D k k k i ? ?? + ? ?? ? ? = = == + = ? ? ? ? ?? ?? ? ? ?= ? ?? ? ? ?? ?? ?? ? + + ? ?? ? ? ?? ?? ? ; ? ? ? ?? ? 1, 2,...n= (F.25) 162 APPENDIX G Derivation of the Solution for the zero Dispersion Case The system of governing transport equations can be written in a matrix format as [11, 38]: [ ] { } { } [ ]{ }c cR v K ct x? ?+ =? ? (G.1) where; ?[ ]? denotes a square matrix and ?{ }? denotes a column vector. The corresponding initial and boundary conditions are: ( ){ } { },0 ,o xc x c e x??= 0 < < ? (G.2) ( ){ } { } ( ) ( ){ }1 1 1 1 0, , 0 ; , 0 1, 2,...i i t i i i o i c t t where B e u t u t t t i n?? = = ? > ? = ? ? > ; ? = ? (G.3) The solution procedure used here is adopted from Quezada et al. [38]. Applying Laplace transform to equation(G.1), we get: [ ] { } [ ]{ } { } [ ]{ }( ,0) pR s p R c x v K px?? + =? (G.4) where; ? s ? is the Laplace variable and ? p ? is the Laplace transformed concentration. Substituting equation(G.2) in (G.4) and rearranging we get: 163 { } [ ] [ ]( ){ } [ ]{ }o xpv K R s p R c ex ???? + ? = ?? (G.5) Now in order to uncouple the system of ordinary differential equations (ODEs) given by equation(G.5), we apply a linear transform procedure described by Clement [11]. i.e. we perform the following matrix operation. { } [ ]{ }p A b= (G.6) where; ?{ }b ? is the concentration in the doubly transformed domain. Applying this transformation equation(G.5) gets modified as: [ ]{ } [ ] [ ]( )[ ]{ } [ ]{ }o xA bv K R s A b R c ex ???? + ? = ?? (G.7) Pre-multiplying equation(G.7) with ?[ ] 1A ? ? we get: { } { } { } [ ] [ ] [ ]( )[ ] { } [ ] [ ]{ } ~ ~ ~ ~1 1 o x bv K b C x where K A K R s A and C A R c e? ? ?? ? ? ?? + = ? ? ?? ? ? ? ? ; = ? = ? ?? ? (G.8) By using the procedure described in Clement [11], we can make ? ~ K? ?? ?? ? ? a diagonal matrix and thus uncouple the system of equations. The corresponding ?[ ]A ? matrix, ?[ ] 1A ? ? matrix, ? ~ K? ?? ?? ? ? matrix and ?{ } ~ C ? vector is given in ?Section 2.2?. Equation(G.8) describes a set of ?n ? independent first-order non-homogeneous ODEs the boundary conditions of which are obtained by a combination of Laplace and Linear transforms of equation(G.3). Laplace transform of equation(G.3) yields: 164 ( ){ } { } ( ){ } ( ) 11 1 11 0, 1 ; 1, 2,... o it si i i i i i p s B e where i ns ? +? = = ? ? ? = ; ? = +?? (G.9) In order to transform the boundary conditions from ? p ? domain to the ?b? domain, we apply the Linear transform given by equation(G.6). This yields: ( ){ } [ ] { }10,b s A ?= ? (G.10) The explicit expression for ? ( )0,ib s ? in equation(G.10) is given as: ( ) ( ) ( ) ( ){ } ( ) 2 2 2 1 2 2 2 1 2 221 1 1 1 1 2 2 , 1 1 1 0, 1, 2,... n o i i i i i n i i i i i i i i t si ii i i i i i y k k sR k sR B e b s s i n ? = + = ? ? +? = =+ ? ? ? ?? ? ? ? ?= +?? ?? ? ?? ? ; ? = ? ? ? ? (G.11) Since equation(G.8) is uncoupled, it can now be written as a set of ?n ? independent equations as: ( ) ( ) ( ) ( ) ( ) 1 1 1 2 2 2 1 1 2 2 2 1 2 1 1 1 , , , 1, 2,... i nx o i i i ii i ii i i i n i i i i i i i i i R c e y kb x s v k sR b x sx k sR k sR i n ?? ? = + = = ? ? ? ? ?? ? ?? + ? ? = ? ? ?? ? + ? ? ? ?? ? ; ? = ? ? ? (G.12) Equation(G.12) can be re-arranged and expressed as: ( ) ( ) ( ) ( ) ( ) 1 1 1 2 2 2 1 1 2 2 2 1 2 1 1 1 , , 1, 1, 2,... i nx o i i i ii i ii i i i n i i i i i i i i i R c e y kb x s k sR b x sx v v k sR k sR i n ?? ? = + = = ? ? ? ? ?? + ? ?+ = ? ?? ? + ? ? ? ?? ? ; ? = ? ? ? (G.13) The general solution to equation(G.13) is given as: 165 ( ) ( ) ( ) ( ) ( )11 1 2 2 2 1 1 1 2 2 2 1 2 1 1 1 , , 1, 2,... i i i nx o i i i i k sR xi i i v i in i i i i i i i i i i i i R c e y k b x s e v k sR k sR k sR i n ?? ? ? + = + = = ? ? ? ? ? ? ?= ? + ? ? ?? ? ? ? + ? ? ? ?? ? ; ? = ? ? ? (G.14) where; i? is a constant. Applying the boundary condition given by equation(G.11) we get: ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) 2 2 2 1 2 2 2 1 2 1 1 2 2 2 1 1 2 2 2 1 2 221 1 1 1 1 2 2 , 1 1 1 , 1 1 1 1 n o i i i i i n i i i i i i i i n o i i i i i i n i i i i i i i i i i i t si ii i i i i i i i y k k sR k sR R c y k k sR k sR k sR B e s v ? = + = ? ? = + = ? ? +? = = = + ? ? ? + ? ? ? ?? ? ? ? ?? = +?? ?? ? ?? ? ? ? ? ? ? ? + ? ?? ? ? ? ?? ? ? ? ? ? ? ? ? 1, 2,...i n ; ? = (G.15) Therefore, the solution in the ?b?domain is: ( ) ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) ( ) 2 2 2 1 2 2 2 1 2 1 1 2 2 2 1 1 2 2 2 221 1 1 1 1 2 2 , 1 1 1 1 1 1 , n o i i i i i i i n i i i i i i i i i in o i i i i i i i i i i i i i i i t si k sR xii i v i i i i k sR x xv y k k sR k sR R c y k k sR k sR k sR B e b x s es e e v ? = + = ? ? = + ? +? ? + = = ? + ?? + ? ? ? + ? ? ? ?? ? ? ? ?= +?? ?? ? ?? ? ? ?? ?? ? ?? ? ? ? + ? ? ? ? ? ? ? ? ( )1 2 1 , 1 1, 2,... n i i i i i i n = ? = ? ? ? ? ? ? ? ? ? ? ? ?? ? ; ? = ? ? (G.16) Inverse Linear transform of equation(G.16) is done to obtain the solution in the Laplace domain (? p ? domain) by using equation(G.6). The solution given by equation(G.16) can be split into two parts and represented as: 166 ( ) ( ) ( ) ( ) ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) 2 2 2 1 2 2 2 1 2 1 1 2 2 2 1 1 22 1 1 1 2 2 1 1 1 , 1 1 1 2 1 1 1 2 , , , ; 1 , , o i i i i i i n i i i i n i i i i i i i i n o i i i i i i i i i i i t si k sR xii i v i i i i k sR x xv i y k k sR k sR R c y k k b x s b x s b x s where B e b x s es e e b x s v ? ? ? ? ? = + = ? ? = + ? + ? + = = ? + ? + ? ? = + ? ?? ? ? ? ?= + ? ?? ? ?? ? ? ?? ?? ? ?? ? ? ?= ? ? ? ? ? ? ( ) ( ) ( ) 2 22 1 2 1 , 1 1, 2,... n i i i i i i i i i i i sR k sR k sR i n = ? = + ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ?? ? ; ? = ? ? (G.17) Using the distributive property, we can apply the inverse linear transform to each of the individual terms and then sum them to get the solution in the ? p ? domain. This is expressed as: { } [ ]{ } [ ]{ } [ ]{ }1 2p A b A b A b= = + (G.18) The first term ?[ ]{ }1A b ? can be evaluated as: { } [ ]{ }1 1p A b= (G.19) The explicit expression for ? ( )1 ,ip x s ? in equation(G.19) is given as: 167 ( ) ( ){ } ( ) ( ) ( ) ( ) 22 1 1 2 2 22 1 2 2 2 1 2 1 2 2 3 3 3 1 3 2 1 11 1 1 , 1 , o i i i t si ii i i i ii i ii k sR xi i vi i i i i i i i i i i i B e y k s p x s e k sR k sR i ? ? ? + ? == + ? + = = = ? ? ??? ? ? ?? ? ? ?+? ? ? ? ? ?= ? ? ? ? ? ?? + ? ? ? ?? ? ; ? ?? ? ? ? 1, 2,...n= (G.20) Using a similar approach the second term ?[ ]{ }2A b ? is evaluated and the explicit expression for ? ( )2 ,ip x s ? is: ( ) ( ) ( ) ( ) ( ) 1 1 2 2 2 1 2 2 1 1 2 1 1 2 2 2 2 3 3 3 1 3 2 1 1 2 1 , , i i i i o i i i i i i k sR xi xv i i i i i i i i i i i i i i i i i R c y k p x s e e v k sR k sR k sR ? ? ? = + ? + ? = = = ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ?= ?? ? ? ?? ? ? ?? ? ? ?? ? ? + ? ? ? ?? ? ? ? ? ? ? ? 1, 2,...i n ; ? = (G.21) Substituting equations(G.20) and (G.21) into equation(G.18) we get the solution in the Laplace domain as: 168 ( ) ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) 22 1 1 2 2 22 1 2 2 2 1 2 1 2 2 3 3 3 1 3 2 1 1 2 2 2 1 1 2 1 11 1 , 1 1 , , , 1 o i i i i i i t si ii i i i ii i ii k sR x i vi i i i i i i i i i i i i o i i i i i i p x s p x s p x s B e y k s e k sR k sR R c y k ? ? ? + ? == + ? + = = = ? ? = + = + = ? ??? ? ? ?? ? ? ?+? ? ? ? ? ? ? ? ? ? ? ?? + ? ? ? ?? ? + ?? ? ? ? ( ) ( ) ( ) ( ) 2 2 1 1 2 1 1 2 2 2 2 3 3 3 1 3 2 1 , 1, 2,... i i i k sR xi xv i i i i i i i i i i i i i i i i e e v k sR k sR k sR i n ? ? ? + ? = = = ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ??? ? ? ?? ? ? ?? ? ? ?? ? ? + ? ? ? ?? ? ? ? ; ? = ? ? ? ? (G.22) The final solution is obtained by taking an inverse Laplace transform of the solution given by equation(G.22). Inverse Laplace transform is performed as follows: ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 1 1 2 1 1 1 2 , , , , , , , 1, 2,... i i i i i i i c x t c x t c x t p x s p x s p x s p x s i n ? ? ? = + = + = + ; ? = ? ? ? (G.23) The first term ? ( )1 1 ,ip x s?? ? can be evaluated as follows. From equations (G.22) and (G.23) we get: 169 ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) 22 1 1 2 2 22 1 2 2 2 1 2 1 2 2 3 3 3 1 3 2 1 11 1 1 1 1 1 , 1 , , o i i i t si ii i i i ii i ii k sR xi i i vi i i i i i i i i i i i B e y k s c x t p x s e k sR k sR ? ? ? + ? == + ? ? ? + = = = ? ? ??? ? ? ?? ? ? ?+? ? ? ? ? ?= = ? ? ? ? ? ?? + ? ? ? ?? ? ?? ? ? ? ? ? 1, 2,...i n ; ? = (G.24) Equation(G.24) can be rewritten as: 17 0 ( ) ( ){ } ( ) ( ) ( ) ( ) ( ) 2 2 2 1 2 2 33 1 1 1 2 1 3 3 2 4 2 4 2 4 4 1 4 2 4 1 4 24 2 4 2 1 1 1 1 1 1 ,0 , , , , , , , , 1 i i o i i i i i i i i i i k sR xi t si v i ii i i i ii i i i i i i i i i i i i i R R i i i i R R y k c x t B e e s a k R s a ? ? = + ? + ? ? + = = = = ? = = ? ? ? ?? ? ? ?? ? ? ?? ? ? ? ? ?= ? ? ? ? ?? ? ? ?? ?+ ? ? + ? ?? ?? ? ? ?? ? ? ? ?? ? ? ? 1 2 1 21 2 1 2 1 2 1 2 1 2 1 , 2 ,, , , 2 1, 2,... ; 0 ; ; 0 i i i ii i i i i i i i i i i i n k when i Rwhere k k k R R R and a when i? ; ? = ? > ? = ? , = ? = ? ? =? (G.25) Note that in equation(G.25) the condition ? 2 4, 0i ik ? ? must be satisfied for the solution to be determinate. Factorization of equation(G.25) gives: 17 1 ( ) ( ){ } ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 1 2 2 33 1 3 2 4 4 1 4 2 4 2 4 1 4 2 4 2 2 4 2 4 2 5 2 5 2 4 5 1 5 2 5 4 5 2 1 1 1 1 ,0 , , , , , , , , , , , , , 1 , 1 i i o i i i i i i i i i i i i k sR x t si v i i i i i i i i i i R R i i i i i i R R i i i i i i i i i i i i i i i i R R y k B e e c x t s a k R s a R a a ? ? = + ? + ? + ? = ? = = ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? = ? ?+ ? ? ?? ? ? ? ? ? ? ? + ? ? ?? ? ? ? ? ? ? 11 2 1 3 1 1 1, 2,... i iii i i i i n = = = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ?? ? ?? ? ; ? = ? ?? (G.26) It must be noted that the solution formulation as given by equation(G.26) is valid only when the condition ? 2 5 2 4, ,i i i i a a? ? is satisfied. Equation(G.26) can be further factorized and simplified as: 17 2 ( ) ( ){ } ( ){ } ( ) ( ) 1 2 2 1 2 1 32 1 2 2 3 3 1 3 1 1 1 1 1 1 1 2 1 11 1 ,0 1 1 , 1, 2,... ; 1 i i o i iii i i i i i i i ii i k sR x t s v i i i c x t y k G h G G i n where e e B s a s G ? ? = = == + ? + ? + ? ? ?? ?= + ? ?? ?? ? ? ?? ? ; ? = ? + = ? ??? ? ( ) ( ) ( ) ( ) ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) 2 4 4 1 4 2 4 2 2 5 2 4 2 5 2 5 2 4 5 1 5 2 5 1 5 2 5 45 2 5 2 2 2 3 3 1 3 2 4 4 1 4 2 4 2 , , , , , , , , , , , , , 1 ,0 1 2 , , , 1 1 i i i i i i i i o i i i i i i i ii i i i R R i i i i i i i i i i i i i i R R i i i i i i R R k sR x t s v i i i i i i i i i i R R a k R R a a e e B s a G k ifh M ? = ? ? = ? = = ? ? ? ? + ? + ? = ? = + ? ? ? ?? ? ? ?? ? ? ? ? + = ? ? = ? ? ? ? ? 0 loopis not executed if loopis executed ? ? ? ? (G.27) 17 3 The term ? 11G ? in equation(G.27) can be further factorized and simplified as: ( ){ } ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 3 3 3 1 2 43 4 1 4 2 2 4 3 2 5 2 4 2 5 2 5 2 4 5 1 5 2 5 1 5 2 5 45 2 5 2 1 1 ,,0 1 1 , , , , , , , , , , , , 1 1 i i i i o i o i i i i i k sR x k sR x t s t sv v i i i ii i ii i i i i i i i i i i i i i i i i i i i i R R i i i i i i R R e e e e B s as a G a k R R a a ? ? ? ? + ? + ? + ? + ? ? = ? = ? = = ? ? ? ? ? ? ?? ? ? ?? ? ?++ ? ? ? ?= ? ? ? ?? ? ? ? ?? ? ? ? ? ? ? ? ( )4 2, i i i R R? ? (G.28) Again it must be noted that equation(G.28) is valid only when the condition ? 2 4 3, ,0i i i a a? ? where; ? 3 3,0i i a ?= ? is satisfied. Using the Appendix H, inverse Laplace expressions for the terms ? 11G ? and ? 12G ? can be evaluated and the solution for ? ( )1 ,ic x t ? is obtained as: 17 4 ( ) ( ){ } [ ] ( ) 1 2 2 1 2 1 32 1 3 2 3 1 2 33 1 1 1 1 1 1 1 1 2 1 11 , ,0 , , 1 1 , 1, 2,... ; , ,o i iii i i i i i i i ii i t i i i i i oi i c x t y k G h G G i n where F x t e F x t t B G ? ? = = == + ? ? ?? ?= + ? ?? ? ? ?? ?? ? ; ? = ? ? = ? ??? [ ] ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ] 3 2 2 4 2 2 4 4 1 4 2 4 2 2 4 3 2 5 2 4 2 5 2 5 2 4 5 1 5 2 5 1 5 2 5 45 2 5 2 3 3 1 2 3 1 2 , , , , , , , , , , , , , , , , , , ,0 , ,1 2 , , , o i i i i i i i o i ti i i i i i i o i ii i i i R R i i i i i i i i i i i i i i i i i R R i i i i i i R R ti i i i i i F x t e F x t t a k R R a a B F x t e F G ? ? ? ? = ? ? = ? = = ? ? ? ? ? ?? ? ? ?? + ?? ? ? ? ? ?? ? ? ? ?? ? ? ? ? = ? ? ? ( ) ( ) 3 2 4 4 1 4 2 4 2 , , , , i i i o i i i i i i i R R x t t k = ? = ? ??? ? ?? (G.29) where; the term ? 1 2 3, ,i i i F ? is given by equation(H.5) in Appendix H. The second term ? ( )1 2 ,ip x s?? ? is evaluated as follows. From equations (G.22) and (G.23) we get: 17 5 ( ) ( ) ( ) ( ) ( ) ( ) 1 1 2 2 2 1 2 2 1 1 2 1 1 2 2 2 2 3 3 3 1 3 2 1 1 2 1 2 1 1 , , , i i i i o i i i i i i k sR xi xv i i i i i i i i i i i i i i i i i i R c y k c x t p x s e e v k sR k sR k sR ? ? ? = + ? + ? ? ? = = = ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ?= = ?? ? ? ?? ? ? ?? ? ? ?? ? ? + ? ? ? ?? ? ? ? ? ? ? ? ? ? 1, 2,...i n ; ? = (G.30) Equation(G.30) can be written as: ( ) ( ) ( ) ( ) ( ) ( ) 1 1 2 2 2 1 2 2 1 1 2 1 1 2 2 2 3 2 3 2 3 3 1 3 2 3 1 3 23 2 3 2 1 1 2 1 1 , , , , , , , , , i i i i i i i i o i i i i i i k sR xi xv i i i i ii i i i i i i i i i i i i i i R R i i i i R R R c y k e ec x t s a R k R s a ? ? = + ? + ?? = = ? = ? = = ? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ??= ? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ?? + ? ? + ? ?? ?? ?? ? ? ? ? ? ? ? 1, 2,...i n ; ? = (G.31) Note: the term ? 1 2,i i a ? is modified as: 17 6 1 2 1 2 1 2 1 1 2 2 , 2 , , 2 2 ; 0 ; 0 ; 0 i i i i i i i i i i k when i R a when i v k when i R ? ? ? > ? ?? = =? ?? + ? < ?? (G.32) Note that the condition ? 2 3, 0i ik ? ? must be satisfied in equation(G.31). Factorization of equation(G.31) yields: 17 7 ( ) ( ){ } ( ) ( ) 1 1 2 2 1 2 12 1 2 2 1 1 2 2 2 2 2 1 1 1 2 1 1 1 , 2 1 , 1, 2,... ; i i i ii i o i i i i i i i ii i k sR x xv i i c x t R c y k G h G G i n where e e s a G ? ? = == + ? + ? ? ? ? ?? ?= + ? ?? ?? ? ? ?? ? ; ? = ? ?? ? ?? ? ? ?? ? + = ? ?? ? ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 3 3 1 3 2 3 2 2 2 4 2 3 2 4 2 4 2 3 4 1 4 2 4 1 4 2 4 34 2 4 2 2 2 1 1 2 2 2 3 3 1 3 2 3 2 , , , , , , , , , , , , , 1 , 2 2 , , , i i i i i i i i i i i i i i i ii i i i R R i i i i i i i i i i i i i i i R R i i i i i i R R k sR x xv i i i i i i i i i i R R s a R k R R a a e e s a G R k ? = ? ? = ? = = ? ? ? ? + ? ? ? = ? = + ? ? ? ?? ? ?? ? ? ? ? ?? ? ?? ? ? ?? ?? + = ? ? ? ? ? ? (G.33) 17 8 Note that the condition ? 2 4 2 3, ,i i i i a a? ? must be satisfied in equation(G.33). The term ? 21G ? in equation(G.33) can be further factorized and simplified as: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 1 1 1 2 2 3 3 1 3 2 2 3 1 2 2 2 4 2 3 2 4 2 4 2 3 4 1 4 2 4 1 4 2 4 34 2 4 2 1 , , 2 1 , , , , , , , , , , , , , i i i i i i i i i i k sR x k sR x x xv v i i i i i ii i i i i i i i i i i i i i i i i i i i i i i R R i i i i i i R R e e e e s a s a G a a R k R R a a ? ? ? + ? + ? ? ? ? = ? ? = ? = = ? ? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? ? ?? ? ? ?? + + = ? ? ? ?? ? ? ?? ? ? ? ? ? ? ( )3 2, i i i R R? ? (G.34) Note that the condition ? 2 4 1 2, ,i i i i a a ?? ? must be satisfied in equation(G.34). ? 21G ? can be further simplified as: 17 9 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 1 1 2 1 2 2 2 1 2 3 2 3 2 3 1 2 2 2 4 2 3 2 4 2 4 2 3 4 1 4 2 4 1 44 2 1 1 , , 1 1 , , 2 1 , , , , , , , , , , 1 1 i i i i i i i i k sR x v x i i i i k sR x v x i i i i i i i i i i i i i i i i i i i i i i i i R R i i i e e s a s a e e s a s a G a a R k R R a a ? ? ? + ?? ? ? ? ? + ?? ? ? = ? = = ? ? ? ?? ? ?+ + ? ? ? ? ? ? ? ? ? ?? + + +? ? ? ?? ?= ? ? ? ?? ? ? ?? ? ? ? ? ? ? ? ? ( ) ( )3 1 3 2 3 2 2 4 3 4 2 , , , , i i i i i ii i i i R R i i i R R = ? ? ? ? ? ? ? (G.35) ? 22G ? in equation(G.33) can be further simplified as: ( ) ( ) ( ) ( ) 2 2 1 1 2 1 2 2 2 3 3 1 3 2 3 2 1 1 , , 2 2 , , , 1 i i i i i k sR x v x i i i i i i i i i i i i R R e e s a s a G R k ? ? + ?? ? ? ? = ? = ? ? ? ?? ? ? ?+ + ? ?? ? = ?? ? ? (G.36) From inverse Laplace transform tables we get: (see [8], p494, eq3). 18 0 ( ) ,1 2 1 2 1 , 1 i ia t i i es a ?? =+? (G.37) Using equation(G.37) and Appendix H, inverse Laplace expressions for the terms ? 21G ? and ? 22G ? can be evaluated and the solution for ? ( )2 ,ic x t ? is obtained as: ( ) ( ){ } [ ] ( ) 1 1 2 2 1 2 12 1 ,1 1 2 2 1 2 2 2 3 2 2 2 2 1 1 1 2 1 1 , , , ,2 1 , 1, 2,... ; , ,i i i ii i o i i i i i i i ii i x a t i i i i i i c x t R c y k G h G G i n where F x t e F x G ? ? ? = == + ? ? ? ? ?? ?= + ? ?? ? ? ?? ?? ? ; ? = ? ? = ? ?? [ ] ( ) ( ) ( ) ( ) ( ) ( ) [ ] ( ) ,1 2 3 3 1 3 2 3 2 2 3 1 2 2 2 4 2 3 2 4 2 4 2 3 4 1 4 2 4 1 4 2 4 34 2 4 2 ,1 1 2 2 1 2 2 2 3 3 1 3 2 , , , , , , , , , , , , , , , ,2 2 , , , , i i i i i i i i i i i i x a t i i ii i i i R R i i i i i i i i i i i i i i i i i i i R R i i i i i i R R x a t i i i i i i i i i i R t e a a R k R R a a F x t e G R k ? ? ? ? ? = ? ? ? = ? = = ? ? ? ? ? ? = ? + ? ? ? ?? ? ? ?? ? ? ? ? ? = ? ? ? ? ( )3 2i i i R= ? (G.38) 18 1 where; the term ? 1 2 3, ,i i i F ? is given by equation(H.5) in Appendix H. The final solution is obtained by substituting equations(G.29) and (G.38) into equation(G.23). ( ) ( ){ } ( ){ } 1 2 2 1 2 1 32 1 1 1 2 2 1 2 12 1 1 1 1 1 1 1 2 1 11 2 2 2 1 1 1 2 1 1 , iii i i i i i i i ii i ii i o i i i i i i ii i c x t y k G h G G R c y k G h G G ? = = == + ? = == + ? ?? ?= + ? ?? ?? ? ? ?? ? ? ?? ? + + ? ?? ?? ? ? ?? ? ? ??? ? ?? 1, 2,...i n ; ? = (G.39) 182 APPENDIX H Evaluation of Inverse Laplace Expressions for the Zero Dispersion Case [ ] ( ){ } ( ) ( ) 1 1 4 1 2 3 4 2 3 1 , , , , 1 , i i o i k sR x t s v i i i i i i e e x t s a ? ? ? + ? + ? ? = +? (H.1) Equation(H.1) can be simplified as: [ ] ( ) ( ) 1 1 1 4 1 2 3 4 2 3 2 3 1 , , , , , , i i oi o i R x R xs t s xk tv v v i i i i i i i i e e ex t e s a s a ? ? ? ? ? ?? ? +? ? ? ? ? ?? ? ? ? ?= ? + +? (H.2) Invoking the ?Second Shifting Theorem? (see [28], p267) we can simplify ? 1 2 3 4, , ,i i i i ? ? as: 18 3 [ ] ( ) ( ) 1 1 1 1 14 1 2 3 4 2 3 , , , 1 , , 1; i o i xk ti i i iv i i i i o o i i R x R x R x R xx t e u t f t e u t t f t t v v v v where f t s a ?? ? ? ? ? ? ? ? ? ? ? ?= ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? = +? (H.3) ? ( )f t ? can be evaluated using equation(G.37). Using equation(G.37) in equation(H.3) we get: [ ] 1 11 , , 2 3 2 31 1 4 1 2 3 4, , , , i ii i i i i o o i R x R xxk a t a t t tv vi iv i i i i o R x R xx t e u t e e u t t e v v ?? ? ? ? ?? ? ? ? ? ?? ? ? ? ?? ? ? ?? ? ? ?= ? ? ? ?? ? ? ? ? ? ? ? (H.4) Equation(H.4) can be compactly expressed as: [ ] [ ] ( ) [ ] 1 4 1 2 3 4 1 2 3 1 2 3 1, 2 31 1 2 3 , , , , , , , , , , , , ; , i o i i i i xk tv i i i i i i i i i i o R xa t vi i i i x t e F x t e F x t t where R xF x t u t e v ?? ? ? ? ?? ?? ? ? ? ? ?= ? ?? ? ? ?= ? ? ?? ? (H.5) 184 APPENDIX I FORTRAN Code SEQUENTIAL C ANALYTICAL SOLUTION FOR MULTI-SPECIES SOLUTE TRANSPORT C CODED BY: VENKATRAMAN SRINIVASAN C AUBURN UNIVERSITY C JULY 2006 C PROGRAM SEQUENTIAL IMPLICIT NONE C C PROBLEM PARAMETERS C INTEGER N,BC,ADS DOUBLE PRECISION V,AX,AY,AZ,LX,LY,LZ,DELX,DELY,DELZ,T,TA,DELT, $ YS,ZS,X DOUBLE PRECISION, ALLOCATABLE :: R(:),KM(:),Y(:),C(:), $ LAM(:),B(:,:),CO(:),MU(:) C C PROGRAM PARAMETERS C CHARACTER FILENAME*80,TITLE*80,COMMENT*80 INTEGER I,I1,I2,I3,I4,I5,PARA,NEGPARA DOUBLE PRECISION,ALLOCATABLE :: KIJ(:,:),KAIIJKL(:,:,:,:),RIJ(:,:) $ ,AIJ(:,:),FIJK(:,:,:) DOUBLE PRECISION TIME,TEMP,TEMP1,TEMP2,TEMP3,TEMP4 C C READ THE INPUT FILE NAME C WRITE(*,*) "ENTER THE NAME OF THE INPUT FILE" READ(*,*) FILENAME C OPEN(5,FILE=FILENAME,FORM='FORMATTED',STATUS='UNKNOWN') OPEN(10,FILE="OUTPUT.DAT",FORM='FORMATTED',STATUS='UNKNOWN') C READ(5,10) TITLE WRITE(10,10) TITLE READ(5,*) N,COMMENT !NUMBER OF SPECIES ALLOCATE (R(N),KM(N),Y(N),LAM(N),B(N,N),CO(N),MU(N)) C C READING THE PROBLEM GEOMETRY C READ (5,*) LX,DELX,LY,DELY,LZ,DELZ,COMMENT !PROBLEM GEOMETRY READ (5,*) YS,ZS,COMMENT !SOURCE GEOMETRY READ (5,*) T,TA,DELT,COMMENT !SIMULATION TIME, PULSE TIME, TIME STEP C C READING THE PROBLEM PARAMETERS C READ(5,*) V,AX,AY,AZ,COMMENT !ADVECTION DISPERSION PARAMETERS READ(5,*) (R(I),I = 1, N), COMMENT !RETARDATION FACTORS READ(5,*) (KM(I),I = 1, N), COMMENT !REACTION CONSTANTS READ(5,*) (Y(I),I = 2, N), COMMENT !YIELD COEFFICIENTS READ(5,*) BC, COMMENT !BOUNDARY CONDITION C IF BC.EQ.0 THEN CONCENTRATION BOUNDARY C IF BC.NE.0 THEN FLUX BOUNDARY READ(5,*) (LAM(I),I = 1, N), COMMENT !SOURCE DECAY DO I = 1, N READ(5,*) (B(I,I1),I1 = 1, I), COMMENT !BATEMAN BOUNDARY VALUES 185 END DO READ(5,*) ADS !REACTION CONDITION C IF ADS.EQ.0, REACTION IN LIQUID PHASE C IF ADS.NE.0, REACTION IN BOTH SOLID AND LIQUID PHASES READ(5,*) (CO(I), I = 1, N), COMMENT !INITIAL CONDITION CONCENTRATION READ(5,*) (MU(I), I = 1, N), COMMENT !INITIAL CONDITION EXPONENTIAL DECAY CLOSE (5) C IF (ADS .NE. 0) THEN DO I = 1, N KM(I) = R(I)*KM(I) END DO END IF C C ALLOCATING THE PROGRAM PARAMETERS C ALLOCATE (KAIIJKL(N,N,0:N,N),FIJK(N,N,-N:N), $ KIJ(N,N),RIJ(N,N),AIJ(N,-N:N)) C CALL CKIJ(KIJ,KM,N) CALL CRIJ(RIJ,R,N) CALL CAIJ(AIJ,KIJ,RIJ,LAM,R,KM,MU,AX,V,N) C C ALLOCATING THE CONCENTRATION VECTOR MATRIX C ALLOCATE (C(N)) C C CHECKING CONDITIONS OF PARAMETERS FOR DETERMINATE SOLUTION C PARA = 0 C IF ANY OF THE CONDITIONS ARE NOT SATISFIED THEN PARA = 1 C IF PARA = 1, THEN PROGRAM STOPS C C CONDITION 1: KIJ(I2,I4) .NE. 0 C DO I2 = 1, N DO I4 = 1, N IF ((I2 .NE. I4) .AND. (R(I2) .EQ. R(I4)) .AND. $ (KM(I2) .EQ. KM(I4))) THEN PARA = 1 WRITE(10,*) "SOLUTION INDETERMINATE DUE TO PARAMETER $VIOLATION" WRITE(10,*) "CONDITION KIJ(",I2,",",I4,") = 0 ; $WHEN RIJ(",I2,",",I4,") = 0) IS VIOLATED" END IF END DO END DO C C CONDITION 2: AIJ(I2,I5) .NE. AIJ(I2,I4) C DO I2 = 1, N DO I4 = 1, N IF ((I2 .NE. I4) .AND. (R(I2) .NE. R(I4))) THEN DO I5 = 1, N IF ((I5 .NE. I2) .AND. (I5 .NE. I4) .AND. $ (R(I5) .NE. R(I2)) .AND. $ (AIJ(I2,I5) .EQ. AIJ(I2,I4))) THEN PARA = 1 WRITE(10,*) "SOLUTION INDETERMINATE DUE TO $PARAMETER VIOLATION" WRITE(10,*) "CONDITION AIJ(",I2,",",I5,") = $AIJ(",I2,",",I4,") ;WHEN RIJ(",I2,",",I4,") AND RIJ(",I2,",",I5,") $ NE 0 IS VIOLATED" END IF END DO END IF END DO END DO C 186 C CONDITION 3&4: AIJ(I2,I4) .NE. LAM(I3) & AIJ(I2,I4) .NE. AIJ(I1,-I2) C DO I1 = 1, N DO I2 = I1, N DO I3 = 1, I1 DO I4 = I1, N IF ((I4 .NE. I2) .AND. (R(I4) .NE. R(I2)) $ .AND. (AIJ(I2,I4) .EQ. LAM(I3))) THEN PARA = 1 WRITE(10,*) "SOLUTION INDETERMINATE DUE TO $PARAMETER VIOLATION" WRITE(10,*) "CONDITION AIJ(",I2,",",I4,") = $LAM(",I3,") ; RIJ(",I4,",",I2,") NE 0 IS VIOLATED" ELSE IF ((I4 .NE. I2) .AND. (R(I4) .NE. R(I2)) $ .AND.(AIJ(I2,I4) .EQ. AIJ(I1,-I2))) THEN PARA = 1 WRITE(10,*) "SOLUTION INDETERMINATE DUE TO $PARAMETER VIOLATION" WRITE(10,*) "CONDITION AIJ(",I2,",",I4,") = $AIJ(",I1,",",-I2,") ; RIJ(",I4,",",I2,") NE 0 IS VIOLATED" END IF END DO END DO END DO END DO C C C STARTING THE TIME LOOP C IF (PARA .EQ. 0) THEN !CHECKING THE PARAMETER CONDITIONS TIME = 0.0 DO WHILE (TIME .LT. T) TIME = TIME + DELT WRITE(10,*) "TIME = ", TIME WRITE (10,20) (I, I = 1, N) C C STARTING THE DISTANCE LOOP C X = 0.0 DO WHILE (X .LT. LX) X = X + DELX C C STARTING THE LOOP FOR CALCULATING THE CONCENTRATION OF ALL SPECIES C DO I = 1, N TEMP = 0.0 DO I1 = 1, I TEMP1 = 1.0 DO I2 = I1+1, I TEMP1 = TEMP1*Y(I2)*KM(I2-1) END DO !I2 LOOP (1) TEMP2 = 0.0 DO I2 = I1, I C C EXPONENTIALLY DECAYING INITIAL CONDITION C PARA = 0 IF (CO(I1) .NE. 0.0) THEN C WRITE (*,*) "INITIAL LOOP" DO I3 = I1, I IF ((I3 .NE. I2) .AND. $ (R(I2) .NE. R(I3))) THEN PARA = 1 TEMP3 = 0.0 NEGPARA = 1.0 DO I4 = I1, I IF ((I4 .NE. I2) .AND. $ (I4 .NE. I3) .AND. $ (R(I2) .NE. R(I4))) THEN 187 IF ((-KIJ(I2,I4)+RIJ(I2,I4) $ *AIJ(I2,I3)) .LT. 0.0) THEN NEGPARA = NEGPARA*(-1) END IF TEMP3 = TEMP3+LOG(ABS( $ (-KIJ(I2,I4)+RIJ(I2,I4)*AIJ(I2,I3)))) END IF END DO !I4 LOOP (1) DO I4 = I1, I IF ((I4 .NE. I2) .AND.(R(I2) .EQ. R(I4))) THEN IF (-KIJ(I2,I4) .LT. 0.0) THEN NEGPARA = NEGPARA*(-1) END IF TEMP3 = TEMP3+LOG(ABS(-KIJ(I2,I4))) END IF END DO !I4 LOOP (2) IF (((KIJ(I2,I3)-RIJ(I2,I3)* $ AIJ(I1,-I2))*R(I2)) .LT. 0.0) THEN NEGPARA = NEGPARA*(-1) END IF TEMP3 = TEMP3+LOG(ABS((KIJ(I2,I3) $ -RIJ(I2,I3)*AIJ(I1,-I2))*R(I2))) CALL CFIJK(FIJK,AIJ,LAM,TIME,X,V, $ AX,R,KM,N,I2,I1,-I2,BC) CALL CFIJK(FIJK,AIJ,LAM,TIME,X,V, $ AX,R,KM,N,I2,I2,I3,BC) TEMP4 = DEXP(-MU(I1)*X)*(DEXP(-AIJ(I1,-I2)*TIME)- $ DEXP(-AIJ(I2,I3)*TIME)) IF (ABS(TEMP4) .GT. 1E+10) TEMP4 = 0.0 C CONDITION FOR DIRICHLET BOUNDARY IF (BC .EQ. 0) THEN TEMP2 = TEMP2 + R(I1)*CO(I1)* $ (FIJK(I2,I1,-I2)-FIJK(I2,I2,I3) $ -TEMP4)/(DEXP(TEMP3)*NEGPARA) C CONDITION FOR NEWMAN BOUNDARY ELSE TEMP2 = TEMP2 + R(I1)*CO(I1)* $ ((1.0+MU(I1)*AX)* $ (FIJK(I2,I1,-I2)-FIJK(I2,I2,I3)) $ -TEMP4)/(DEXP(TEMP3)*NEGPARA) END IF END IF END DO !I3 LOOP IF (PARA .EQ. 0) THEN TEMP3 = 0.0 NEGPARA = 1 DO I3 = I1, I IF ((I3 .NE. I2) .AND. (R(I2) .EQ. R(I3))) THEN IF ((-KIJ(I2,I3)) .LT. 0.0) THEN NEGPARA = NEGPARA*(-1) END IF TEMP3 = TEMP3+LOG(ABS((-KIJ(I2,I3)))) END IF END DO !I3 LOOP (2) TEMP3 = TEMP3+LOG(R(I2)) CALL CFIJK(FIJK,AIJ,LAM,TIME,X,V, $ AX,R,KM,N,I2,I1,-I2,BC) TEMP4 = DEXP(-MU(I1)*X-AIJ(I1,-I2)*TIME) IF (ABS(TEMP4) .GT. 1E+10) TEMP4 = 0.0 C CONDITION FOR DIRICHLET BOUNDARY IF (BC .EQ. 0) THEN TEMP2 = TEMP2+R(I1)*CO(I1)*(-FIJK(I2,I1,-I2)+ $ TEMP4)/(DEXP(TEMP3)*NEGPARA) C CONDITION FOR NEWMAN BOUNDARY ELSE TEMP2 = TEMP2+R(I1)*CO(I1)* $ (-(1.0+MU(I1)*AX)*FIJK(I2,I1,-I2) $ +TEMP4)/(DEXP(TEMP3)*NEGPARA) END IF 188 END IF END IF C C BOUNDARY CONDITION C DO I3 = 1, I1 IF (B(I1,I3) .NE. 0.0) THEN C WRITE (*,*) "BOUNDARY LOOP" PARA = 0 DO I4 = I1, I IF ((I4 .NE. I2) .AND. (R(I2) .NE. R(I4))) THEN PARA = 1 TEMP3 = 0.0 NEGPARA = 1 DO I5 = I1, I IF ((I5 .NE. I2) .AND. (I5 .NE. I4) .AND. $ (R(I2) .NE. R(I5)))THEN IF ((-KIJ(I2,I5)+RIJ(I2,I5) $ *AIJ(I2,I4)) .LT.0.0) THEN NEGPARA = NEGPARA*(-1) END IF TEMP3 = TEMP3+LOG(ABS((-KIJ(I2,I5)+ $ RIJ(I2,I5)*AIJ(I2,I4)))) END IF END DO !I5 LOOP (1) DO I5 = I1, I IF ((I5 .NE. I2) .AND. (R(I2) .EQ. R(I5))) THEN IF ((-KIJ(I2,I5)) .LT. 0.0) THEN NEGPARA = NEGPARA*(-1) END IF TEMP3 = TEMP3+LOG(ABS((-KIJ(I2,I5)))) END IF END DO !I5 LOOP (2) IF ((-KIJ(I2,I4)+RIJ(I2,I4)*LAM(I3)) .LT.0.0)THEN NEGPARA = NEGPARA*(-1) END IF TEMP3 = TEMP3+LOG(ABS((-KIJ(I2,I4) $ +RIJ(I2,I4)*LAM(I3)))) CALL CKAIIJKL(KAIIJKL,AIJ,LAM,TIME $ ,X,V,AX,R,KM,TA,N,I2,I3,0,I3,BC) CALL CKAIIJKL(KAIIJKL,AIJ,LAM,TIME $ ,X,V,AX,R,KM,TA,N,I2,I2,I4,I3,BC) TEMP2 = TEMP2 + B(I1,I3)*(KAIIJKL(I2,I3,0,I3)- $ KAIIJKL(I2,I2,I4,I3))/ $ (DEXP(TEMP3)*NEGPARA) END IF END DO !I4 LOOP (1) IF (PARA .EQ. 0) THEN TEMP3 = 0.0 NEGPARA = 1 DO I4 = I1, I IF ((I4 .NE. I2) .AND. (R(I2) .EQ. R(I4))) THEN IF ((-KIJ(I2,I4)) .LT.0.0) THEN NEGPARA = NEGPARA*(-1) END IF TEMP3 = TEMP3+LOG(ABS((-KIJ(I2,I4)))) END IF END DO !I4 LOOP (2) CALL CKAIIJKL(KAIIJKL,AIJ,LAM,TIME $ ,X,V,AX,R,KM,TA,N,I2,I3,0,I3,BC) TEMP2 = TEMP2 + B(I1,I3)*KAIIJKL(I2,I3,0,I3)/ $ (DEXP(TEMP3)*NEGPARA) END IF END IF END DO !I3 LOOP END DO !I2 LOOP (2) TEMP = TEMP + TEMP1*TEMP2 END DO !I1 LOOP C(I) = TEMP 189 END DO C C WRITING THE OUTPUT C WRITE (10,30) X, (C(I), I = 1, N) C C END OF DISTANCE LOOP C END DO C C END OF TIME LOOP C END DO END IF 10 FORMAT (A80) 20 FORMAT(' DISTANCE', (" SPECIES",I3.3)) 30 FORMAT (<(N+1)>E15.5) CLOSE(10) END PROGRAM C C SUBROUTINE FOR CALCULATING THE KIJ COEFFICIENTS C SUBROUTINE CKIJ(KIJ,KM,N) IMPLICIT NONE INTEGER N,I1,I2 DOUBLE PRECISION KIJ(N,N),KM(N) DO I1 = 1, N DO I2 = 1, N KIJ(I1,I2) = KM(I1) - KM(I2) END DO END DO RETURN END C C SUBROUTINE FOR CALCULATING THE RIJ COEFFICIENTS C SUBROUTINE CRIJ(RIJ,R,N) IMPLICIT NONE INTEGER N,I1,I2 DOUBLE PRECISION RIJ(N,N),R(N) DO I1 = 1, N DO I2 = 1, N RIJ(I1,I2) = R(I1) - R(I2) END DO END DO RETURN END C C SUBROUTINE FOR CALCULATING THE AIJ COEFFICIENTS C SUBROUTINE CAIJ(AIJ,KIJ,RIJ,LAM,R,KM,MU,AX,V,N) IMPLICIT NONE INTEGER N,I1,I2 DOUBLE PRECISION AIJ(N,-N:N),KIJ(N,N),RIJ(N,N),LAM(N),R(N),KM(N), $ MU(N),AX,V C C INITIALISING THE VALUE OF AIJ C DO I1 = 1, N DO I2 = 0, N AIJ(I1,I2) = 0.0 END DO END DO DO I1 = 1, N DO I2 = -N, N IF (I1 .NE. I2) THEN IF (I2 .LT. 0) THEN AIJ(I1,I2) = (-MU(I1)**2.0*AX*V-MU(I1)*V+KM(-I2))/R(-I2) ELSE IF (I2 .EQ. 0) THEN 190 AIJ(I1,I2) = LAM(I1) ELSE IF (R(I1) .NE. R(I2)) THEN AIJ(I1,I2) = KIJ(I1,I2) / RIJ(I1,I2) END IF END IF END DO END DO RETURN END C C SUBROUTINE FOR CALCULATING KAIIJKL COEFFICIENTS C SUBROUTINE CKAIIJKL(KAIIJKL,AIJ,LAM,TIME,X,V,AX,R,KM,TA,N,I1,I2, $ I3,I4,BC) IMPLICIT NONE INTEGER N,I1,I2,I3,I4,BC DOUBLE PRECISION KAIIJKL(N,N,0:N,N),AIJ(N,-N:N),LAM(N),TIME,X,V, $ AX,R(N),KM(N),TA,FIJK(N,N,-N:N),TEMP CALL CFIJK(FIJK,AIJ,LAM,TIME,X,V,AX,R,KM,N,I1,I2,I3,BC) KAIIJKL(I1,I2,I3,I4) = FIJK(I1,I2,I3) IF (TIME .GT. TA) THEN TEMP = TIME-TA CALL CFIJK(FIJK,AIJ,LAM,TEMP,X,V,AX,R,KM,N,I1,I2,I3,BC) KAIIJKL(I1,I2,I3,I4) = KAIIJKL(I1,I2,I3,I4) - $ DEXP(-LAM(I4)*TA)*FIJK(I1,I2,I3) END IF RETURN END C C SUBROUTINE FOR CALCULATING THE FIJK COEFFICIENTS C SUBROUTINE CFIJK(FIJK,AIJ,LAM,TIME,X,V,AX,R,KM,N,I1,I2,I3,BC) IMPLICIT NONE INTEGER N,I1,I2,I3,BC,PARA DOUBLE PRECISION FIJK(N,N,-N:N),AIJ(N,-N:N),LAM(N),TIME,X,V,AX $ ,R(N),KM(N),WIJK,PI,CALEXF,TEMP1,TEMP2, $ TEMP3,TEMP4,TEMP5,TEMP6,TEMP7 PI = 22.0/7.0 C C CALCULATING THE WIJ TERMS C CALL CWIJK(WIJK,V,R,AX,KM,AIJ,N,I1,I2,I3,PARA) C C CASE FOR CONCENTRATION BOUNDARY CONDITION C IF (BC .EQ. 0) THEN C C CASE FOR NORMAL ERROR FUNCTION C IF (PARA .EQ. 0) THEN TEMP1 = -AIJ(I2,I3)*TIME + (V-WIJK)*X/(2.0*AX*V) TEMP2 = (R(I1)*X-WIJK*TIME)/(4.0*R(I1)*AX*V*TIME)**0.5 FIJK(I1,I2,I3) = 0.5*CALEXF(TEMP1,TEMP2) C C CASE WHEN X << AX C C IF (X .LT. (100.0*AX)) THEN TEMP3 = -AIJ(I2,I3)*TIME + (V+WIJK)*X/(2.0*AX*V) TEMP4 = (R(I1)*X+WIJK*TIME)/(4.0*R(I1)*AX*V*TIME)**0.5 FIJK(I1,I2,I3)=FIJK(I1,I2,I3)+0.5*CALEXF(TEMP3,TEMP4) C END IF C C CASE FOR COMPLEX ERROR FUNCTION C ELSE TEMP1 = R(I1)*X/(4.0*R(I1)*AX*V*TIME)**0.5 TEMP2 = WIJK*TIME/(4.0*R(I1)*AX*V*TIME)**0.5 TEMP3 = -AIJ(I2,I3)*TIME+V*X/(2.0*AX*V) CALL CEXF(TEMP1,TEMP2,TEMP3,TEMP4,TEMP5) 191 FIJK(I1,I2,I3) = DCOS(WIJK*X/(2.0*V*AX))*TEMP4- $ DSIN(WIJK*X/(2.0*V*AX))*TEMP5 END IF C C CASE FOR FLUX BOUNDARY CONDITION C ELSE IF (BC .EQ. 1) THEN C C CASE FOR (KM(I1)/R(I1)) .NE. AIJ(I2,I3) C IF ((KM(I1)/R(I1)) .NE. AIJ(I2,I3)) THEN C C CASE FOR NORMAL ERROR FUNCTION C IF (PARA .EQ. 0) THEN TEMP1 = -AIJ(I2,I3)*TIME + (V-WIJK)*X/(2.0*AX*V) TEMP2 = (R(I1)*X-WIJK*TIME)/(4.0*R(I1)*AX*V*TIME)**0.5 TEMP3 = -AIJ(I2,I3)*TIME + (V+WIJK)*X/(2.0*AX*V) TEMP4 = (R(I1)*X+WIJK*TIME)/(4.0*R(I1)*AX*V*TIME)**0.5 TEMP5 = V*X/(V*AX)-KM(I1)*TIME/R(I1) TEMP6 = (R(I1)*X+V*TIME)/(4.0*R(I1)*V*AX*TIME)**0.5 FIJK(I1,I2,I3) = V*(CALEXF(TEMP1,TEMP2)/(V+WIJK) $ +CALEXF(TEMP3,TEMP4)/(V-WIJK)) $ +2.0*V**2.0/(4.0*R(I1)*V*AX* $ (KM(I1)/R(I1)-AIJ(I2,I3)))* $ CALEXF(TEMP5,TEMP6) C C CASE FOR COMPLEX ERROR FUNCTION C ELSE TEMP1 = R(I1)*X/(4.0*AX*V*R(I1)*TIME)**0.5 TEMP2 = WIJK*TIME/(4.0*AX*V*R(I1)*TIME)**0.5 TEMP3 = -AIJ(I2,I3)*TIME+V*X/(2.0*AX*V) TEMP4 = V*X/(V*AX)-KM(I1)*TIME/R(I1) TEMP5 = (R(I1)*X+V*TIME)/(4.0*R(I1)*AX*V*TIME)**0.5 CALL CEXF(TEMP1,TEMP2,TEMP3,TEMP6,TEMP7) FIJK(I1,I2,I3) = 2.0*V/(V**2.0+WIJK**2.0)*((TEMP6*V-TEMP7* $ WIJK)*DCOS(X*WIJK/(2.0*V*AX))-(TEMP6*WIJK $ +TEMP7*V)*DSIN(X*WIJK/(2.0*V*AX)))+V**2.0 $ /(2.0*R(I1)*V*AX*(KM(I1)/R(I1)-AIJ(I2,I3) $ ))*CALEXF(TEMP4,TEMP5) END IF C C CASE FOR KM(I1) .EQ. AIJ(I2,I3) C ELSE C C ORIGINAL FORMULATION C C TEMP1 = -KM(I1)*TIME/R(I1) C TEMP2 = (R(I1)*X-V*TIME)/(4.0*R(I1)*V*AX*TIME)**0.5 C TEMP3 = -KM(I1)*TIME/R(I1)-(R(I1)*X-V*TIME)**2.0/ C $ (4.0*R(I1)*V*AX*TIME) C TEMP4 = -KM(I1)*TIME/R(I1)+V*X/(V*AX) C TEMP5 = (R(I1)*X+V*TIME)/(4.0*R(I1)*V*AX*TIME)**0.5 C FIJK(I1,I2,I3) = 0.5*CALEXF(TEMP1,TEMP2)+ C $ (V**2.0*TIME/(PI*R(I1)*AX*V))**0.5*DEXP(TEMP3) C $ -0.5*(1.0+X*V/(AX*V)+V**2.0*TIME/(R(I1)*AX*V)) C $ *CALEXF(TEMP4,TEMP5) C C NEW FORMULATION C REFORMULATING TO AVOID THE CONDITION WHEN THE ARGUMENT IN THE C SECOND TERM (EXPONENTIAL FUNCTION) BECOMES ZERO AND THE ARGUMENT C IN THE THIRD TERM (ERFC FUNCTION) BECOMES GREATER THAN 25 C WRITE (*,*) "EQUAL" TEMP1 = -KM(I1)*TIME/R(I1) TEMP2 = (R(I1)*X-V*TIME)/(4.0*R(I1)*V*AX*TIME)**0.5 TEMP3 = -(R(I1)*X-V*TIME)**2.0/(4.0*R(I1)*V*AX*TIME)-X/AX 192 TEMP4 = -KM(I1)*TIME/R(I1)+V*X/(V*AX) TEMP5 = (R(I1)*X+V*TIME)/(4.0*R(I1)*V*AX*TIME)**0.5 TEMP6 = 0.0 IF (TEMP5 .LT. 20.0) THEN FIJK(I1,I2,I3) = 0.5*CALEXF(TEMP1,TEMP2)+DEXP(TEMP4)*( $ (V**2.0*TIME/(PI*R(I1)*AX*V))**0.5*DEXP $ (TEMP3)-0.5*(1.0+X*V/(AX*V)+V**2.0*TIME/ $ (R(I1)*AX*V))*CALEXF(TEMP6,TEMP5)) ELSE FIJK(I1,I2,I3) = 0.5*CALEXF(TEMP1,TEMP2) END IF END IF END IF IF (ABS(FIJK(I1,I2,I3)) .GE. 1E+8) THEN C WRITE (*,*) "FLAG, FIJK" FIJK(I1,I2,I3) = 0.0 END IF RETURN END C C SUBROUTINE FOR CALCULATING THE WIJK COEFFICIENTS C SUBROUTINE CWIJK(WIJK,V,R,AX,KM,AIJ,N,I1,I2,I3,PARA) IMPLICIT NONE INTEGER N,I1,I2,I3,PARA DOUBLE PRECISION WIJK,V,R(N),AX,KM(N),AIJ(N,-N:N),TEMP PARA = 0 TEMP = V**2.0+4.0*R(I1)*V*AX*(KM(I1)/R(I1)-AIJ(I2,I3)) IF (TEMP .LT. 0.0) THEN PARA = 1 C WRITE (*,*) "NEGATIVE" END IF WIJK=(ABS(TEMP))**0.5 RETURN END C C FUNCTION TO CALCULATE EXP(TEMP1)*ERFC(TEMP2) C FUNCTION CALEXF(TEMP1,TEMP2) IMPLICIT NONE C C NEW ERFC ROUTINE C DOUBLE PRECISION TEMP1,TEMP2,TEMP,CALEXF C CALL CALERF(TEMP2,TEMP,1) CALEXF = 0.0 IF ((TEMP .NE. 0.0) .AND. ((TEMP1+LOG(TEMP)) .LT.500.0))THEN CALEXF = DEXP(TEMP1+LOG(TEMP)) END IF C C VANGENUCHTEN ERFC ROUTINE C C DOUBLE PRECISION TEMP1,TEMP2,CALEXF,EXF C IF (TEMP2 .GT. 3.5) THEN C CALEXF = 0.0 C ELSE IF (TEMP2 .LT. -3.5) THEN C CALEXF = 2.0*DEXP(TEMP1) C ELSE C CALEXF = EXF(TEMP1,TEMP2) C END IF C RETURN END C C FUNCTION TO CALCULATE EXP(A) ERFC(B) C FUNCTION EXF(A,B) C 193 C PURPOSE: TO CALCULATE EXP(A) ERFC(B) C IMPLICIT REAL*8 (A-H,O-Z) EXF=0.0 IF((DABS(A).GT.170.).AND.(B.LE.0.)) RETURN IF(B.NE.0.0) GO TO 1 EXF=DEXP(A) RETURN 1 C=A-B*B IF((DABS(C).GT.170.).AND.(B.GT.0.)) RETURN IF(C.LT.-170.) GO TO 4 X=DABS(B) IF(X.GT.3.0) GO TO 2 T=1./(1.+.3275911*X) Y=T*(.2548296-T*(.2844967-T*(1.421414-T*(1.453152-1.061405*T)))) GO TO 3 2 Y=.5641896/(X+.5/(X+1./(X+1.5/(X+2./(X+2.5/(X+1.)))))) 3 EXF=Y*DEXP(C) 4 IF(B.LT.0.0) EXF=2.*DEXP(A)-EXF RETURN END C C SUBROUTINE FOR CALCULATING COMPLEX ERROR FUNCTION C SUBROUTINE CEXF(A,B,Z,U,V) C C COMPLEX ERFC-FUNCTION: U+IV=EXP(Z)ERFC(A+IB) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION W(10),H(10) DATA W/.4622437,.2866755,.1090172,.02481052,.3243773D-2, 1.2283386D-3,.7802556D-5,.1086069D-6,.4399341D-9,.2229394D-12/, 2H/.2453407,.7374737,1.234076,1.738538,2.254974, 32.788806,3.347855,3.944764,4.603682,5.387481/ C C ---------- X=DABS(A) Y=DABS(B) U=0.0 V=0.0 E=0.0 XYZ=Y*Y+Z-X*X IF(DABS(XYZ).GT.170.)RETURN COS=DCOS(2.*X*Y) SIN=DSIN(2.*X*Y) IF((X+Y).LT.6.) GO TO 2 T=0.0 DO 1 K=1,10 T=T+W(K)*((X/((Y-H(K))**2+X*X))+(X/((Y+H(K))**2+X*X))) 1 V=V+W(K)*((Y-H(K))/((Y-H(K))**2+X*X)+(Y+H(K))/((Y+H(K))**2+X*X)) U=.3183099*DEXP(XYZ)*(T*COS-V*SIN) V=.3183099*DEXP(XYZ)*(-T*SIN-V*COS) IF(X.LE.0.)U=DEXP(DMIN1(Z,1.7D2)) GO TO 8 2 IF(X.GT.2.5) GO TO 3 T=1./(1.+.3275911*X) U=T*(.2548296-T*(.2844967-T*(1.421414-T*(1.453152-1.061405*T)))) GO TO 4 3 U=.5641896/(X+.5/(X+1./(X+1.5/(X+2./(X+2.5/(X+3./(X+1.))))))) 4 IF(Y.LE.0.) GO TO 7 IF(X.LE.0.) V=-.3183099*Y IF(X.LE.0.) GO TO 5 U=U-.1591549*(1.D0-COS)/X V=V-.1591549*SIN/X 5 NT=12.+2.*Y DO 6 I=1,NT P=I ARG=P*Y F1=X*(DEXP(ARG)+DEXP(-ARG)) 194 F2=0.5*P*(DEXP(ARG)-DEXP(-ARG)) EX=.6366198*DEXP(-0.25*P*P)/(4.*X*X+P*P) U=U-EX*(2.*X-F1*COS+F2*SIN) 6 V=V-EX*(F1*SIN+F2*COS) V=V*DEXP(Z-X*X) 7 U=U*DEXP(Z-X*X) 8 IF(B.LT.0.) V=-V IF(A.LT.0.) U=2.*EXF(Z,E)-U RETURN END C C SUBROUTINE FOR CALCULATING ERFC C SUBROUTINE CALERF(ARG,RESULT,JINT) C------------------------------------------------------------------ C C This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) C for a real argument x. It contains three FUNCTION type C subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX), C and one SUBROUTINE type subprogram, CALERF. The calling C statements for the primary entries are: C C Function Parameters for CALERF C call ARG Result JINT C C ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0 C ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1 C ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2 C C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C Latest modification: March 19, 1990 C C------------------------------------------------------------------ INTEGER I,JINT CS REAL DOUBLE PRECISION 1 A,ARG,B,C,D,DEL,FOUR,HALF,P,ONE,Q,RESULT,SIXTEN,SQRPI, 2 TWO,THRESH,X,XBIG,XDEN,XHUGE,XINF,XMAX,XNEG,XNUM,XSMALL, 3 Y,YSQ,ZERO DIMENSION A(5),B(4),C(9),D(8),P(6),Q(5) C------------------------------------------------------------------ C Mathematical constants C------------------------------------------------------------------ CS DATA FOUR,ONE,HALF,TWO,ZERO/4.0E0,1.0E0,0.5E0,2.0E0,0.0E0/, CS 1 SQRPI/5.6418958354775628695E-1/,THRESH/0.46875E0/, CS 2 SIXTEN/16.0E0/ DATA FOUR,ONE,HALF,TWO,ZERO/4.0D0,1.0D0,0.5D0,2.0D0,0.0D0/, 1 SQRPI/5.6418958354775628695D-1/,THRESH/0.46875D0/, 2 SIXTEN/16.0D0/ C------------------------------------------------------------------ C Machine-dependent constants C------------------------------------------------------------------ CS DATA XINF,XNEG,XSMALL/3.40E+38,-9.382E0,5.96E-8/, CS 1 XBIG,XHUGE,XMAX/9.194E0,2.90E3,4.79E37/ DATA XINF,XNEG,XSMALL/1.79D308,-26.628D0,1.11D-16/, 1 XBIG,XHUGE,XMAX/26.543D0,6.71D7,2.53D307/ C------------------------------------------------------------------ C Coefficients for approximation to erf in first interval C------------------------------------------------------------------ CS DATA A/3.16112374387056560E00,1.13864154151050156E02, CS 1 3.77485237685302021E02,3.20937758913846947E03, CS 2 1.85777706184603153E-1/ CS DATA B/2.36012909523441209E01,2.44024637934444173E02, CS 1 1.28261652607737228E03,2.84423683343917062E03/ 195 DATA A/3.16112374387056560D00,1.13864154151050156D02, 1 3.77485237685302021D02,3.20937758913846947D03, 2 1.85777706184603153D-1/ DATA B/2.36012909523441209D01,2.44024637934444173D02, 1 1.28261652607737228D03,2.84423683343917062D03/ C------------------------------------------------------------------ C Coefficients for approximation to erfc in second interval C------------------------------------------------------------------ CS DATA C/5.64188496988670089E-1,8.88314979438837594E0, CS 1 6.61191906371416295E01,2.98635138197400131E02, CS 2 8.81952221241769090E02,1.71204761263407058E03, CS 3 2.05107837782607147E03,1.23033935479799725E03, CS 4 2.15311535474403846E-8/ CS DATA D/1.57449261107098347E01,1.17693950891312499E02, CS 1 5.37181101862009858E02,1.62138957456669019E03, CS 2 3.29079923573345963E03,4.36261909014324716E03, CS 3 3.43936767414372164E03,1.23033935480374942E03/ DATA C/5.64188496988670089D-1,8.88314979438837594D0, 1 6.61191906371416295D01,2.98635138197400131D02, 2 8.81952221241769090D02,1.71204761263407058D03, 3 2.05107837782607147D03,1.23033935479799725D03, 4 2.15311535474403846D-8/ DATA D/1.57449261107098347D01,1.17693950891312499D02, 1 5.37181101862009858D02,1.62138957456669019D03, 2 3.29079923573345963D03,4.36261909014324716D03, 3 3.43936767414372164D03,1.23033935480374942D03/ C------------------------------------------------------------------ C Coefficients for approximation to erfc in third interval C------------------------------------------------------------------ CS DATA P/3.05326634961232344E-1,3.60344899949804439E-1, CS 1 1.25781726111229246E-1,1.60837851487422766E-2, CS 2 6.58749161529837803E-4,1.63153871373020978E-2/ CS DATA Q/2.56852019228982242E00,1.87295284992346047E00, CS 1 5.27905102951428412E-1,6.05183413124413191E-2, CS 2 2.33520497626869185E-3/ DATA P/3.05326634961232344D-1,3.60344899949804439D-1, 1 1.25781726111229246D-1,1.60837851487422766D-2, 2 6.58749161529837803D-4,1.63153871373020978D-2/ DATA Q/2.56852019228982242D00,1.87295284992346047D00, 1 5.27905102951428412D-1,6.05183413124413191D-2, 2 2.33520497626869185D-3/ C------------------------------------------------------------------ X = ARG Y = ABS(X) IF (Y .LE. THRESH) THEN C------------------------------------------------------------------ C Evaluate erf for |X| <= 0.46875 C------------------------------------------------------------------ YSQ = ZERO IF (Y .GT. XSMALL) YSQ = Y * Y XNUM = A(5)*YSQ XDEN = YSQ DO 20 I = 1, 3 XNUM = (XNUM + A(I)) * YSQ XDEN = (XDEN + B(I)) * YSQ 20 CONTINUE RESULT = X * (XNUM + A(4)) / (XDEN + B(4)) IF (JINT .NE. 0) RESULT = ONE - RESULT IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT GO TO 800 C------------------------------------------------------------------ C Evaluate erfc for 0.46875 <= |X| <= 4.0 C------------------------------------------------------------------ ELSE IF (Y .LE. FOUR) THEN XNUM = C(9)*Y XDEN = Y DO 120 I = 1, 7 XNUM = (XNUM + C(I)) * Y XDEN = (XDEN + D(I)) * Y 196 120 CONTINUE RESULT = (XNUM + C(8)) / (XDEN + D(8)) IF (JINT .NE. 2) THEN YSQ = AINT(Y*SIXTEN)/SIXTEN DEL = (Y-YSQ)*(Y+YSQ) RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT END IF C------------------------------------------------------------------ C Evaluate erfc for |X| > 4.0 C------------------------------------------------------------------ ELSE RESULT = ZERO IF (Y .GE. XBIG) THEN IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 300 IF (Y .GE. XHUGE) THEN RESULT = SQRPI / Y GO TO 300 END IF END IF YSQ = ONE / (Y * Y) XNUM = P(6)*YSQ XDEN = YSQ DO 240 I = 1, 4 XNUM = (XNUM + P(I)) * YSQ XDEN = (XDEN + Q(I)) * YSQ 240 CONTINUE RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5)) RESULT = (SQRPI - RESULT) / Y IF (JINT .NE. 2) THEN YSQ = AINT(Y*SIXTEN)/SIXTEN DEL = (Y-YSQ)*(Y+YSQ) RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT END IF END IF C------------------------------------------------------------------ C Fix up for negative argument, erf, etc. C------------------------------------------------------------------ 300 IF (JINT .EQ. 0) THEN RESULT = (HALF - RESULT) + HALF IF (X .LT. ZERO) RESULT = -RESULT ELSE IF (JINT .EQ. 1) THEN IF (X .LT. ZERO) RESULT = TWO - RESULT ELSE IF (X .LT. ZERO) THEN IF (X .LT. XNEG) THEN RESULT = XINF ELSE YSQ = AINT(X*SIXTEN)/SIXTEN DEL = (X-YSQ)*(X+YSQ) Y = EXP(YSQ*YSQ) * EXP(DEL) RESULT = (Y+Y) - RESULT END IF END IF END IF 800 RETURN C---------- Last card of CALERF ---------- END 197 Sample INPUT file 3 SPECIES SEQUENTIAL CHO EXAMPLE PROBLEM 3 !NUMBER OF SPECIES 200.0 10.0 1.0 1.0 1.0 1.0 !PROBLEM GEOMETRY, DISTANCE:DISTANCE STEP (X,Y,Z) 1.0 1.0 !SOURCE GEOMETRY 200.0 200.0 200.0 !SIMULATION TIME, PULSE TIME, TIME STEP 1.0 0.18 0.0 0.0 !VELOCITY, ALPHA (X,Y,Z) 2.0 1.0 1.0 !RETARDATION FACTORS 0.01 0.1 0.0 !DECAY CONSTANTS 1.0 1.0 !YIELD COEFFICIENTS 0 !BOUNDARY CONDITION TYPE 0.0 0.0 0.0 !SOURCE DECAY LAMDA 1.0 !SPECIE1 0.0 0.0 !SPECIE2 0.0 0.0 0.0 !SPECIE3 0 !SORPTION CONDITION 0.0 0.0 0.0 !INITIAL CONDITION CONCENTRATION 0.0 0.0 0.0 !INITIAL CONDITION EXPONENTIAL DECAY Sample OUTPUT file 3 SPECIES SEQUENTIAL CHO EXAMPLE PROBLEM TIME = 200.000000000000 DISTANCE SPECIES001 SPECIES002 SPECIES003 0.10000E+02 0.90500E+00 0.58963E-01 0.36037E-01 0.20000E+02 0.81902E+00 0.75434E-01 0.10554E+00 0.30000E+02 0.74122E+00 0.76529E-01 0.18225E+00 0.40000E+02 0.67080E+00 0.72352E-01 0.25685E+00 0.50000E+02 0.60707E+00 0.66636E-01 0.32629E+00 0.60000E+02 0.54940E+00 0.60739E-01 0.38986E+00 0.70000E+02 0.49721E+00 0.55131E-01 0.44766E+00 0.80000E+02 0.44984E+00 0.49951E-01 0.50007E+00 0.90000E+02 0.39123E+00 0.44654E-01 0.54725E+00 0.10000E+03 0.19746E+00 0.31598E-01 0.58224E+00 0.11000E+03 0.19071E-01 0.10525E-01 0.58123E+00 0.12000E+03 0.17579E-03 0.18928E-02 0.54780E+00 0.13000E+03 0.10664E-06 0.30148E-03 0.50195E+00 0.14000E+03 0.48551E-11 0.47855E-04 0.44975E+00 0.15000E+03 0.14559E-16 0.75816E-05 0.39182E+00 0.16000E+03 0.28198E-23 0.11956E-05 0.32774E+00 0.17000E+03 0.34881E-31 0.18650E-06 0.25691E+00 0.18000E+03 0.27368E-40 0.28315E-07 0.17887E+00 0.19000E+03 0.13558E-50 0.40066E-08 0.97222E-01 0.20000E+03 0.42276E-62 0.46934E-09 0.32185E-01 198 APPENDIX J Evaluation of Integral The integral in equation(77) can be evaluated as: 2 2 0 0 exp 4 4 ' exp 2 t t x x x xx v xk D Dx vxI f d d DD ? ? ? ? ? ? ? ? ??? = = 3/2 = = ? ?? ?? +? ? ? ? ? ?= = ? ? ? ?? ? (J1) Applying Laplace transform to equation(J1) we get: [ ] 2 2 exp 4 41 exp 2 x x xx v xk D Dx vxI D sD ? ? ? ?? 3/2 ? ?? ?? ?? + ? ?? ?? ? ? ? ? ?= ? ? ? ?? ? ? ?? ? ? ? ? ? (J2) Where, ? is the Laplace transform operator and ?s? is the Laplace variable. Equation(J2) can be expressed as: [ ] 2 2 exp 41 exp exp2 4 x x xx x Dx vx vI k s D DD ?? ?? 3/2 ? ?? ?? ? ?? ?? ?? ? ? ? ? ? ? ?= ? +? ? ? ? ? ?? ?? ? ? ? ? ?? ? ? ? ? ? (J3) The second term within the Laplace operator can be evaluated by using Selby [40] (See equation(82), p: 497) as: 199 2 exp 4 2 expx x x x D D x s x D ? ? ? 3/2 ? ?? ?? ? ?? ? ? ?? ?? ?? ? ? ? = ?? ?? ? ? ? ? ?? ?? ?? ? ? ?? ? ? ? ? (J4) The entire expression within the Laplace operator can be evaluated by using Selby [40] (See equation(11), p: 491) and equation(J4) as: [ ] 22 exp exp2 4 x xx vx x vI s k s D DD ? ?? ?? ? ? ?? ? ? ?= ? + +? ? ? ?? ? ? ?? ? ? ?? ?? ?? ?? (J5) Inverse Laplace transform of the above equation yields[5]: ( ) ( ) 1/2 1/ 2 1/ 2 1/ 2 1/ 2 1/ 2 41 4exp 1 1 2 2 41 4exp 1 1 2 2 x x x x x x x x kx vt kx vI erfc v vt kx vt kx verfc v vt ? ? ? ? ? ? ? ? ? ?? ?? ? ? +? ?? ?? ?? ?? ?? ? ? ?? ? ? ? ? ?= ? + ? ?? ? ? ?? ? ? ?? ?? ?? ? ? ?? ?? ? ? ?? ?? ? ? ?? ? ? ? ?? ? + +? ?? ?? ?? ?? ? ? ?? ? ? ? + + + ? ?? ? ? ?? ? ? ?? ?? ? ? ?? ?? ? ? ?? ? , xx Dwhere v? ? ? ? ? ? ? ? ? ?? ? ? ? = (J6) 200 APPENDIX K Derivation of the three-dimensional solution to zero longitudinal dispersion case The governing transport equation is: 2 2 2 2y z c c c cv D D kc t x y z ? ? ? ?+ ? ? = ? ? ? ? ? (K1) The initial and boundary conditions are: 0 ( , , ,0) 0 0 , , (0, , , ) , , 02 2 2 2 0 , 0 ( , , , )lim 0 ( , , , )lim 0 y z c x y z x y z Z Z Y Yc y z t c z y t otherwise t c x y z t y c x y z t z ??? ??? = ? < < ? ?? < < ? ?? < < ? = ? < < ? < < ? > = ? > ? = ? ? = ? (K2) Applying Laplace transform to ( , , , )c x y z t in equation(K1) gives: 2 2 2 2 ( ) ; y z y z y z p p p s k p x y z v D Dwhere and v v ? ? ? +?? ?? = ? ? ? ? ? = ? = (K3) where; ?s? is the Laplace variable and ?p? is the concentration in the Laplace domain. The boundary conditions get modified as: 201 0(0, , ) , 2 2 2 2 0 ( , , )lim 0 ( , , )lim 0 y z c Z Z Y Yp y z z y s otherwise p x y z y p x y z z ??? ??? = ? < < ? < < = ? = ? ? = ? (K4) Equation(K3) can be interpreted as a transient two dimensional diffusive reactive transport problem. Its boundary conditions given by equation(K4) represent an instantaneous pulse of a plane source. The solution to this problem without the decay term can be readily deduced from Hunt [25] by ignoring the advection term and reducing the problem to two dimensions. Thus the solution to the above problem without the reaction term is given by: ( ) ( ) ( ) ( ) 1/2 1/ 2 1/ 2 1/ 2 '( , , ) ( , ) ( , )4 2 2; ( , ) 2 2 2 2( , ) 2 2 o y z y y y z z z cp x y z f x y f x z s Y Yy y where f x y erf erf x x Z Zz z and f x z erf erf x x ? ? ? ? = ? ?? ? ? ?+ ? ? ?? ? ? ? = ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ?+ ? ? ?? ? ? ? = ? ? ?? ? ? ? ? ?? ? ? ?? ? ? ? ? ?? ? (K5) Now we make use of a method similar to the Danckwert?s method described by Crank [15] to include the reaction term. If 'p is the solution for the diffusion problem without reaction; the solution for the same problem with a first order reaction (with a rate constant s kv+? ?? ? ? ? ) for the same initial and boundary condition is given as: 202 'exp s kp p xv? + ?? ?= ?? ?? ? ? ?? ? (K6) This can be easily verified by checking if the solution p satisfies the governing equation and the initial and boundary conditions. Since 'p is the solution to equation (K3) without the reaction term, it must satisfy: 2 2 2 2 ' ' ' 0 y z p p p x y z ? ? ??? ?? = ? ? ? (K7) Also, differentiating p with respect to x, y and z to the respective orders yields: 2 2 2 2 2 2 2 2 'exp ' exp 'exp 'exp dp s k s k s k dpx p x dx v v v dx p s k px y v y p s k px z v z + ? + ? ? + ?? ? ? ? ? ?= ? ? + ?? ? ? ?? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? + ??? ?= ?? ?? ? ? ?? ?? ? ? ? + ??? ?= ?? ?? ? ? ?? ?? ? (K8) From equations (K6), (K7) and (K8) we get equation(K3). This proves that solution p satisfies the governing differential equation. To check for the initial condition, we substitute x = 0 in equation(K6). When x = 0, the exponential term becomes unity and hence equation(K6) reduces to p = 'p ; thus the initial condition is satisfied. In order to check for the boundary condition in the y direction, we need to take the derivative of the solution p with respect to y. This is given as: 'expdp s k dpxdy v dy? + ?? ?= ?? ?? ? ? ?? ? (K9) In the limiting case, when y approaches ? ?, the derivative of 'p with respect to y becomes 0. From equation(K9) we can conclude that the derivative of p with respect 203 to y also becomes 0. This proves that the boundary condition in the y direction is satisfied. On similar lines, we see that the boundary condition is also satisfied in the z direction. Hence it is proved that equation(K6) is the solution for the system of equations described by (K3) and (K4). Equation(K6) can be written as: ( , , ) exp ( , ) ( , )4o y zc s kp x y z x f x y f x zs v? + ?? ? = ? ? ?? ? ? ?? ? (K10) Inverse Laplace transform of equation(K10) gives the final solution as [40] (See equation(61), p: 495): ( , , , ) ( , ) ( , ) ( , )8 ; ( , ) 2exp ; 0 1 oo x y z o x cc x y z t f x t f x y f x z k x xwhere f x t u t v v xand u t isthe step function givenby v xif t x vu t xv if t v = ? ? ? ? = ? ?? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??? ?? = ? ? ?? ? ? > ?? (K11)