STELLA AND CHEMICAL REACTIONS

John C. Drake

Stella provides an ideal platform for representing and evaluating reactions of various types. It can be advantageously used to facilitate analysis of chemical equilibria in complex systems and to examine evolutionary changes in dynamic systems. The following models will examine some of these approaches; they are not intended to provide an introduction to underlying chemical or geological principles, but rather to illustrate the ease with which object oriented simulation modeling can facilitate problem definition and problem solving in complex systems.

  Reaction Kinetics

            In order to truly represent reactions kinetically it is necessary to know all of the steps in the reaction path as well as the rate laws and constants for each.  In some instances it may be possible by measuring reactants and products to formulate an empirical rate law and determine an empirical rate constant. For a discussion of reaction kinetics see, for example:

Langmuir, D., 1997, Aqueous Environmental Chemistry; Chapter 2: Chemical Kinetics; Upper Saddle River, NJ, Prentice Hall

Lasaga, A. and Kirkpatrick, R.J. (Eds.), 1981, Kinetics of Geochemical Processes; Reviews in Mineralogy, Vol. 8,  Min. Soc. Amer., Washington, D.C.

Stumm, W. and Morgan, J., 1996, Aquatic Chemistry 3rd edition, Chapter 2: Chemical Thermodynamics and Kinetics; New York, John Wiley

Rate laws are frequently classified by the order of reaction corresponding to the exponential of concentration affecting the reaction rate.  Thus a non-reversible reaction consuming reactant A could be:

            Zeroth order    dA/dt = -k(A)0 = -k

            First order        dA/dt = -k(A)1 = -k(A)

            Second order    dA/dt = -k (A)2

Reaction order can often be evaluated by monitoring reactant (or product) concentration over time. For reversible reactions it is necessary to know rate constants and order for both forward and reverse reaction.  Thus a first order reversible reaction with A as the reactant and B as the product would be first order with respect to the forward reaction and first order with respect to the reverse reaction. A reaction can be considered elementary (Langmuir, pp. 56 - 60) if the reaction order can be determined from the stoichiometry of the equation.  For the non-reversible elementary reaction:

                                     aA + bB products

the rate of the forward reaction would be  vf = kf(A)a(B)b

For a reversible elementary reaction:                            

aA + bB = cC + dD

the rate law for the forward reaction would be vf = kf(A)a(B)b  and the rate of the reverse reaction would be vr = kr(C)c(D)d.  At equilibrium the rates of the forward and reverse reactions are equal. If vf =vr, then

kf(A)a(B)b = kr((C)c(D)d    and Keq = kf/kr = (C)c(D)d/(A)a(B)b

where Keq is the equilibrium constant for the reaction.

  Although reaction rate constants (or actual reaction paths) are not known for many geological processes it is illustrative to develop kinetic based models The following exercises will introduce you to two different reaction scenarios: a) consecutive reactions; b) concurrent reactions. You should construct models of each reaction mechanism and can experiment by using reaction rates and starting compositions different from the ones given. Consider all reactions to be elementary in nature. Note that although stocks are best considered in units of mass, they will correspond to units of concentration if all volumes are equal and constant. You should make this assumption for the present problem. Converters can be added in subsequent model development to accommodate such factors as variable volumes, non-ideality of systems (i.e. concentration is different from activity), various units of concentration, etc.

CONSECUTIVE REVERSIBLE REACTIONS:

                                    1)         A + B   =  C      

                                    2)         C  =  D

Using starting compositions of 100 for all reactants and products run this model when

                                    k1f = 0.005;  k1r = 0.002;  k2f =  0.05;  k2r = 0.005

Using the same reaction rates run the model when A = B =100 and C =D = 0; when A = B = 0 and C = D = 100        

CONCURRENT REVERSIBLE REACTIONS:

                        1)         A + B = C

                        2)         A = D

Using starting compositions of 100 for all reactants and products run this model under the following conditions:

                                    k1f = 0.005;  k1r = 0.002;  k2f = 0.05;  k2r = 0.005

Run the same model using the same reaction rates when A = B = 100 and C = D = 0; when A = B = 0 and C = D = 100

For each of the above models experiment with different starting compositions and different relative reaction rates to see the effect on final composition. Are final compositions in agreement with equilibrium calculations? 

Chemical Fluxes:

  As another example of a Stella application to a dynamic simulation, and also as an example of importing/exporting data to another application, I have constructed the following simplified model to illustrate diffusion within a sediment column. This is a 7 layer model in which there is an initial quantity of a pollutant in layer 4 but not in any of the adjacent layers. I have assumed that Fick's first law of diffusion is applicable, that the bulk diffusion coefficient takes into consideration tortuosity and porosity, and that each layer is 1 unit thick. All assigned values are arbitrary. This model has been initially developed such that the composition of both boundary layers (layer 1 and layer 7) are constant (i.e. flux in =  flux out). It should also be noted that stocks represent quantity.  By making each stock identical in volume (e.g. area and thickness) that quantity can be used as a proxy for concentration. Otherwise values of quantity can be changed to concentrations by using converters.  It is also convenient to use the output of this model to demonstrate transfer of data to spreadsheets to expand data analysis and graphing. We will review the model GSAFlux has been constructed using Fick’s first law of diffusion:

                                                JA = -D dc/dx

Where J is the rate of diffusive transfer per unit area, D is the diffusion coefficient and dc/dx is the concentration gradient.  This model is an example of “main chain” construction” in which a sequence of stocks is linked together using a similar architecture. After analyzing the construct of this model, try changing boundary conditions and creating new graphs. Three possibilities are: a) close the system at both ends so that there is no material leaving layers 1 and 7; b) provide a source of material to layer 4 that replenishes the initial source as rapidly as it is depleted; c) vary the value of the diffusion coefficient.

Radioactive Decay:

Stella also provides a convenient platform with which to model radioactive decay processes. This is really nothing more than another application of first order reaction kinetics

We will use the following terminology:

            N =  amount of parent;             NO =  amount of parent at time zero

            D =  amount of daughter           DO = amount of daughter at time zero

            D* = radiogenically produced daughter

            l = decay constant                   T1/2 = half life

 

-dN/dt =  lN   where “-“ indicates a decrease;

 

-ln N = lT + C   where C = -ln NO

-ln N = lT – ln NO   or  ln N – ln NO =  -lT

ln N/NO = e -lT             N = NOe-lT

Providing no daughters are present other than those produced by radioactive decay the number of radiogenic daughters at any time can be calculated:

D* = NO – N = NO – NOe-lT   = NO (1- e-lT)

Half life can be incorporated into these equations as a substitute for the decay constant because of the following relationship:

T1/2 is the time required for  NO to decrease to 1/2 NO

1/2 NO = NO e-lT1/2  so that T =  T1/2 when N = 1/2 NO

N = NOe-lT   and 1/2 NO = NOe-lT1/2

ln 1/2 + lnNO/NO =  ln1/2  =  -lT1/2 or ln 2 = lT1/2  or  T1/2 = (ln 2)/l

It is sometimes convenient to relate D* to N using the following relationships

D* = NO – N where N = NOe-lT  or  NO = NelT

D* = NelT  - N = N(elT – 1)

D = DO + D* = DO + N(elT –1) = DO + NelT - N

These relationships can be used to calculate a radiometric age (T)

T = (1/l)*ln((D – DO + N)/N

BACK TO AGENDA