The numerical model developed in this thesis uses a method similar in nature to operator-splitting (OS) called integrated operator-splitting (IOS). Integrated operator-splitting was chosen because it allows different reactions to be easily included. This section describes the IOS method used in this model and how it relates to standard OS.

Standard operator-splitting divides equation (3.1) into two parts.

	(3.2)
	(3.3)

The method integrates these equations separately and then recombines them to form equation (3.4) below:

(3.4)

Where:	D C = change in control volume concentration over a timestep
	Cf = concentration at end of timestep
	Co = concentration at beginning of timestep
	D CA&D = change in concentration due to advection and dispersion
	D CR = change in concentration due to reactions

Operator-splitting solves the non-reactive solute transport term (D C_A&D) first. The result is an intermediate solution (C*) that includes only the effects of advection and dispersion over the time interval D t. It then uses this intermediate solution as the initial concentration to solve the reaction term (D C_R) for the concentration at the end of the timestep. Figure 3.2.1 graphically illustrates this process.

Figure 3.2.1 Standard Operator-Splitting Procedure

The dashed arrows in Figure 3.2.1 illustrate dependencies. In OS D C_A&D is a function of the initial concentration, C_o. The reaction term, D C_R, is a function of the intermediate solution, C*.

Other researchers have used methods similar to the Operator-Splitting method. Kinzelbach and Schäfer (1991) modified the standard OS method to partially re-couple equations (3.2) and (3.3). The value of D C_R enters into the solution of the D C_A&D term as an added sink/source term in equation (3.2). In turn, D C_A&D enters into the solution of D C_R as an extra constant sink/source reaction in equation (3.3). Figure 3.2.2 illustrates the dependencies associated with this method of solution.

Figure 3.2.2 Operator-Splitting Procedure modified by Kinzelbach and Schäfer.

Kinzelbach and Schäfer iterate steps 1 and 2 until the final concentration, C_f, converges.

Integrated Operator-Splitting borrows from each of these ideas. The difference between this method and the others is that IOS only re-couples equation (3.3). (See Figure 3.2.3.) The calculation of D C_A&D is independent of kinetic reactions. However, D C_A&D still enters into the calculation of D C_R. In addition, this method makes a distinction between kinetic and equilibrium reactions. D C_R is split into two components, kinetic reactions (D C_K) and equilibrium reactions (D C_E). Also, the model simulates simultaneous movement and reaction of multiple species. Equation (3.5) rewrites equation (3.4) in a form used by IOS.

(3.5)

Where:	D C = change array of concentration values over a time interval
	Cf = array of concentration values at end of timestep
	Co = array of concentration values at beginning of timestep
	D CA&D = change in concentration array due to solute flux or advection-dispersion
	D CK = change in array due to kinetic reactions
	D CE = change in array due to equilibrium reactions

Vector notation is appropriate due to the simulation of multiple species. These species react with each other via kinetic or equilibrium equations. It would be inappropriate to write a separate version of equation (3.5) for each solute since these solutes can be interdependent.

Similar to Kinzelbach and Schäfer’s method, IOS uses an iterative procedure to calculate the concentration at the end of a timestep. Table 3.2.1 and Figure 3.2.3 present this iterative procedure.

Table 3.2.1 IOS Order of Calculation

Figure 3.2.2 Integrated Operator-Splitting Procedure

The solution to each term on the right hand side of equation (3.5) uses a different set of underlying assumptions and methods of solution. The change in concentration due to advection and dispersion (D C_A&D) is a function of the length of the time step and the concentration at the beginning of the timestep. The model uses a finite difference method corrected for numerical dispersion developed by Poulsen (1994) to solve this term. The solution to the D C_A&D term requires the calculation of fluid velocities. The model estimates these velocities using a finite difference technique described later in this section. The change in concentration due to kinetic reactions (D C_K) is a function of the time interval length and the concentration at both the beginning and end of the time interval. The linear integrated method described below is used to calculate this term. The change in concentration due to equilibrium reactions is solely a function of the initial concentrations and changes in concentration due to kinetic and flux components (C’).

The following sections describe how the IOS method calculates the components of equation (3.5). The first section describes how the aquifer is discretized. The following section discusses how the model estimates groundwater velocities. A discussion of how the model calculates advective and dispersive fluxes follows in the next section. The subsequent section presents the Linear Integrated method and how it works. Equilibrium calculations appear in the following section. The final section discusses stability criteria. For additional information on the IOS method, see Appendix B.