//********************************************************************************* // Standard ORCHESTRA object definitions // // Version 2014 // // J.C.L. (Hans) Meeussen // Nuclear Research and Consultancy Group, Petten, The Netherlands // www.nrg.eu/re // meeussen@nrg.eu // // This file contains the standard ORCHESTRA objects for implementing // chemical reaction and mass transport models. // // The ENTITY, PHASE, and REACTION classes are the three basic building blocks // for all chemical models. // // Orchestra Objects (classes) are defined in terms of: // 1) Variables: // @Var: // @Var: pH 7 // // 2) Mathematical expressions: // @Calc:(1, " = ") // @Calc:(1, "pH = -{H+.logact}") // // 3) Unknown - equation pairs: // @Uneq2: unknown:(name:, pH, type:, lin) equation:(name:, H+.tot, tol: , 1e-15) // // The definition of the phase and entity object classes can be found within the orchestra2.jar file. // // Changes in the 2014 version // - mineral object with variable tolerance that is adapted during iteration // - updated solid solution object // // Changes in the 2013 version // - added a set of @logKreaction() objects which require // a logK value as input rather than a k value. The logK reactions can be mixed with normal @reactions() in a single // input file to ensure bakcward compatibility. // The graphical user interface understands both formats and will keep existing format by default. // - The donnan model and nica sites were made electroneutral by adding dummy proton entities. // In this way the overall system remains electroneutral when surface area is changed. // The oxide surfaces are are expressed in terms of neutral surface sites, but here the composition of // the DDL is not explicitly taken into consideration. --> // - A mineral entity now has a .con variable to make it consistent with other entities, and to make it possible // to use it as a primary entity. // - A disperse() class was added, that calculates dispersion when combined with the convection object // - a balance_charge(cation,anion) object was added to calculate // // changes in the 2011 version // - half life and decay objects added // - minerals are now phases as well as entities (made optional for efficiciency reasons) // // // Changes in the 2008 version // - conversion to log activities // - changed all .act variables to .logact // - removed all .sum variables (obsolete) // - changed activity correction factors to log equivalents (f-1.logact etc.) // - defined all "unkown" variables as global variables with default values // - changed solid solution iteration from log to linear // - added surface reaction with 6 components // - changed primary entity object for 7 parameters, (used for pH pe) so that it uses given default value for unknown // - updated and tested ddl objects (new version can handle charge unbalance in solution) //************************************************************************************** @logactivities: // indicates that the chemistry calculator uses log activities @GlobalVar: T 298.15 // Temperature with default value @Class: calc_conc(entity, factor){ @Calc:(2,".con = (10^{.logact}) * ") } @Class: calc_conc_with_davies(name, charge, phase){ @Calc:(2,".con = (10^({.logact} + {f.log}))") @Calc:(2,"f.sum = {f.sum} + {.con}") @Calc:(2,". = {.con}") } @Class: reaction(entity, k){ @Var: .k @Calc:(1,".logact = log10({.k})") } // The reaction classes with 1 - 10 reactants @Class: reaction(entity, k, c1, n1){ @reaction(, ) @link(, , ) } @Class: reaction(entity, k, c1, n1, c2, n2){ @reaction(, ) @link(, , ) @link(, , ) } @Class: reaction(entity, k, c1, n1, c2, n2, c3, n3){ @reaction(, ) @link(, , ) @link(, , ) @link(, , ) } @Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4){ @reaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } @Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5){ @reaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } @Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6){ @reaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } @Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7){ @reaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } @Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7, c8, n8){ @reaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } @Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7, c8, n8, c9, n9){ @reaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } @Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7, c8, n8, c9, n9, c10, n10){ @reaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } // The logk reaction allows to use logK values as input for reaction definitions // these can be mixed with ordinary reactions in input files @Class: logKreaction(entity, logk){ @Var: .logk @Calc:(1,".logact = {.logk}") } // The reaction classes with 1 - 8 reactants @Class: logKreaction(entity, logk, c1, n1){ @logKreaction(, ) @link(, , ) } @Class: logKreaction(entity, logk, c1, n1, c2, n2){ @logKreaction(, ) @link(, , ) @link(, , ) } @Class: logKreaction(entity, logk, c1, n1, c2, n2, c3, n3){ @logKreaction(, ) @link(, , ) @link(, , ) @link(, , ) } @Class: logKreaction(entity, logk, c1, n1, c2, n2, c3, n3, c4, n4){ @logKreaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } @Class: logKreaction(entity, logk, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5){ @logKreaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } @Class: logKreaction(entity, logk, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6){ @logKreaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } @Class: logKreaction(entity, logk, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7){ @logKreaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } @Class: logKreaction(entity, logk, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7, c8, n8){ @logKreaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } @Class: logKreaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7, c8, n8, c9, n9){ @logKreaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } @Class: logKreaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7, c8, n8, c9, n9, c10, n10){ @logKreaction(, ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) @link(, , ) } @Class: davies() This object creates a chemical system with a calculated ionic strength. It defines the ionic strength I as an unknown, and the difference between the calculated and estimated ionic strength I.sum as the equation value which has to become zero during the iteration. { @davies(0.1) @Uneq2: unknown:(name:, I, delta:, 1e-1, type:, lin, step: , .1, min:, 1e-9, max:, 10, iia:, true) equation:(name:, I.sum, tol: , 1e-4) } @Class: f(charge) This object represents the ion activity correction factor for ions, according to the Davies model. The "log activity" (f.log) of this object represents the activity factor, the "sum" f.sum the total amount of ions with this charge.
Hans Meeussen { //***** act = activity coefficient, sum = sum * Z^2 of ions with charge z @Var: f.log 0 @Var: f.sum 0 @Calc:(1,"f.sum = 0") @Calc:(1,"f.log = log10(I.act ^ (*) )") @Calc:(3,"I.sum = I.sum + ({f.sum} * (*))") @Calc:(3,"chargebalance = chargebalance + ({f.sum} * ())") } @Class: davies (i) This object calculates the activity correction factors for ions charged 1-9 according to the Davies equation. It also calculates the ionic strength I, and the total sum concentration of ions with a certain charge.
Hans Meeussen, June 2002. { @globalvar: I // estimated ionic strength @Var: I.act 1 @Var: I.sum 0 // sum = sum conc * z^2 for all species @Var: I.calc 0 // I as calculated from actual speciation @Var: EC 0 // Electrical conductivity millimhos/cm @var: chargebalance 0 // chargebalance of dissolved species @Calc: (1, "chargebalance = 0") // initialises chargebalance //@Calc:(1, "I.act = 10^(0.51* ( (sqrt(I)/(sqrt(I)+1))-(0.23*I)) )" ) // In this version I is maximized at 1 //@Calc:(1, "I.act = 10^(0.51* ( (sqrt(min(I,1))/(sqrt(min(I,1))+1))-(0.23*min(I,1))) )" ) @Calc:(1, "I.act = 10^(0.51* ( (sqrt(I)/(sqrt(I)+1))-(0.23*I)) )" ) // Modified version according to E. Samson et al Comput. Mat. Science 15 (1999) 285-294 // This modification linearly changes the 0.2 factor from 0.2 at I=0 to 1.5 at I = 1.5. // which gives a better estimation of ion activity at higher ionic strengths. // @Calc:(1, "I.act = exp( (0.0360*sqrt(I*1000)/(1+0.0309*sqrt(I*1000)) - (-4.17e-5*I*1000+.02)*0.0360*I*1000/sqrt(1000)) )") @Calc:(5, "EC = I.calc/0.013") // Calculate electrical conductivity from calculated I @Calc:(1, "I.sum = 0") @Calc:(5, "I.sum = (I.sum -2*I)") @Calc:(5, "I.calc = (I.sum/2)") @Var: f0.log 0 @Var: f0.sum 0 @f(1) @f(2) @f(3) @f(4) @f(5) @f(6) @f(7) @f(8) @f(9) @f(-1) @f(-2) @f(-3) @f(-4) @f(-5) @f(-6) @f(-7) @f(-8) @f(-9) } @Class: aqphase(){diss} @Class: species(name, charge){ @entity(, "@aqphase()", 1) @calc_conc_with_davies(, , "@aqphase()") } @Class: mineral(name) This object represents an automatically precipitating/dissolving mineral. A mineral is an entity in the "min" phase, which can be combined with a standard (formation)reaction object. The ".logact" variable of a mineral entity is its log IAP + log K value. Because, in contrast with ordinary entities, the amount of this mineral entity cannot be calculated directly from its activity we introduce an extra unknown and equation that is solved by the standard iteration procedure. The amount of mineral can now be calculated from the estimated unknown. If the unknown is negative, the mineral is not present and the unknown represents its log undersaturation. If the unknown is positive it directly represens the amount of mineral present. The equation that has to become zero to reach convergence: unknown < 0 equation = log (IAP*K) - unknown // difference between unknown and log undersaturation (should go to zero) unknown >= 0 equation = log (IAP*K) // log saturation index (should go to zero) May 2014 Hans Meeussen. { // @phase(, min, 1) @entity(, min) @Var: .un -1e-3 @Var: .eq 0 @Var: .si 0 @Var: minTol 1e-3 @Calc:(2,".si = {.logact}") @Calc:(2,".con = if({.un} >= 0, ({.un}), 0) ") @Calc:(2,".min = {.con}") //@Calc:(2,".eq = if(0 > {.un}, ({.si} - {.un}), {.si})") @Calc:(2,".eq = if(0 > {.un}, ({.si} - {.un}), {.si}- {.un} * minTol )") @Uneq3: unknown:(name:, .un, delta:, 1e-6, type:, lin, step: , .1) equation:(name:, .eq, tol: , 1e-9, si:, .si ) } @Class: solid_solution(name, parent_phase) Solid solutions in are solved in ORCHESTRA by introducing an extra unknown and equation. When negative the unknown represents the undersaturation, and indicates that the amount of solid solutio is zero. The accompanying equation then is the difference between unknown and (log) undersaturation. When positive the unknown represents the amount of solid solution. The accomponying equation then is the difference between 1 and the sum of fractions { @GlobalVar: .un -1 // the unknown that is used in the iteration @GlobalVar: .eq 0 // the calculated equation that should become zero during iteration @Var: .est_sum 0 // the estimated sum (is calculated from unknown) // a solid sulution is a phase in which entities can be present @link_phase(, , .est_sum) // a solid solution is also an entity (in its own phase) //(so has an activity and mass balance and can be used in reactions) // logact = 0 @Globalvar: .logact 0 @entity(, , 0, 0) // calculate estimated amount of solid solution (est_sum) directly from the unknown @Calc:(1,".est_sum = if({.un} >= 0, {.un}, 0) ") // calculate equation (error that should become zero) from unknown and calculated sum of fractions // this implies that the sum of fractions needs to become 1 when solid solution exists. // if solid solution does not exist (unknown<0), the equation is undersaturion. // if solid solution exists, (unknown>0) estimated sum should be equal to calculated sum // if solid solution is completely dissolved, (unknown<=0) calculated sum should be zero // the variable . contains the sum of fractions that should go to 1. @Calc:(5, ".eq = if(0> {.un}, ( log({.})/1000 - {.un}), 1-{.})") // iterate for a value of the unknown that results in an equation value of zero. @Uneq2: unknown:(name:, .un, delta:, 1e-3, type:, lin, step: , .1 ) equation:(name:, .eq, tol: , 1e-6) } @Class: e_layer (name, fase) An entity that represents an electrostatic layer with .logact = log boltzmann factor . = charge balance. By placing this entity in an e fase, its sum is directly converted to the units Coulomb/m2 { @entity (, , 0) @Var: .psi .1 @GlobalVar: T 298.15 // Temperature with default value //* calculate psi from boltzman factor name.act @Calc:(1,".psi = (-{.logact}* log(10))/(96484.56/(8.31441*T))") } @Class: cdmusicsurf(name, parent_phase, surface_area) { @cdmusicsurf(, , , 1.1, 5) } @Class: cdmusicsurf(name, parent_phase, surface_area, c1, c2) { @globalvar: .c1 // default values for capacitance @globalvar: .c2 // @phase(, , ) @phase(_e) @link_phase(, _e, 96484.56) // The electrostatic layers are child phases of the _e fase, and linked to this phase // via the conversion factor 96484.56. // This makes the charge balances (_e0._e) available in coulomb/m2 after phase 4. @e_layer (_e0, _e) @e_layer (_e1, _e) @e_layer (_e2, _e) // a diffuse double layer that uses the potential of the e2 layer //@ddl(, _e2) @ddl_old(, _e2) @Var: .eq1 0.0 // eq1= ddl.sigma + e0.sigma + e1.sigma + e2.sigma @Calc:(5,".eq1 = {.ddl} + {_e0._e} + {_e1._e} + {_e2._e} ") @Uneq2: unknown:(name:, _e0.logact, delta:, 1e-6, type:, lin, step: , 1) equation:(name:, .eq1, tol: , 1e-6) @Var: .eq2 0.0 // eq2 = (psi0-psi1)* c1 - e0.sigma @Calc:(5,".eq2 = ({.c1}*({_e0.psi}-{_e1.psi}))-{_e0._e}") @Uneq2: unknown:(name:, _e1.logact, delta:, 1e-6, type:, lin, step: , 1) equation:(name:, .eq2, tol: , 1e-6) @Var: .eq3 0.0 // eq3 = (psi1-psi2)* c2 - e0.sigma - e1.sigma @Calc:(5,".eq3 = {.c2}*({_e1.psi}-{_e2.psi})-{_e0._e} - {_e1._e}") @Uneq2: unknown:(name:, _e2.logact, delta:, 1e-6, type:, lin, step: , 1) equation:(name:, .eq3, tol: , 1e-6) } @Class: bstern(name, phase, surface_area) A Stern model has two electrostatic layers (e0 and e1) and a diffuse double layer that uses the potential of the e1 layer. { @Var: .c 3.9 // capacitance in F/m2 @phase(, , ) @phase(_e) @link_phase(, _e, 96484.56) // The electrostatic layers are positioned in the _e fase. This phase has // the sums (_e0._e)are in the dimension coulomb/m2 after stage 4 @e_layer (_e0, _e) @e_layer (_e1, _e) @ddl(, _e1) @Var: .eq1 0.0 @Calc:(5, ".eq1 = .ddl + _e0._e + _e1._e") @Uneq2: unknown:(name:, _e0.logact, delta:, 1e-6, type:, lin, step: , 1) equation:(name:, .eq1, tol: , 1e-6) @Var: .eq2 0.0 // eq2 = (psi0-psi1)* c - e0.sigma @Calc:(5, ".eq2 = (_e0.psi-_e1.psi) * .c - _e0._e") @Uneq2: unknown:(name:, _e1.logact, delta:, 1e-6, type:, lin, step: , 1) equation:(name:, .eq2, tol: , 1e-6) } //********************************************************************************** // The transport object classes // ********************************************************************************* //------------------------------------------------------------------------------------ // Transport objects translated to new expressions (9 September 2002) @Class: update_mass (name) Updates the mass of a component in dynamic systems (transport or kinetics). It is normally used within the file update_mass.inp. It assumes that the phases "tot" and "d" are used to represent total and delta (per time unit) amounts of component masses.
Hans Meeussen { @globalvar: .tot 0 @globalvar: .d 0 @Var: dt 1 // Add the mass change .d (mol/s) multiplied by the time step dt(s) // to the total amount in the cell and reset name.d.sum to zero. @Calc:(1, ".tot = {.tot} + {.d} * dt") @Calc:(1, ".d = 0") } @Class: convection () This object contains the general part of convection and defines the amount of water transported between cell 1 and cell 2. It should be used in convection.inp before the component specific objects. { @Var: dwater 0 @Var: 1.J 0 //The amount of water moving from cell 1 to cell 2 equals the //flow rate (l/s). @Calc:(1, "dwater=1.J") } @Class: convec (name) This object contains the component specific part of convection and defines the amount of component (mol/s) transported between cell 1 and cell 2. { @globalvar: 1..diss 0 @globalvar: 2..diss 0 @Var: 1..d 0 @Var: 2..d 0 @Var: dmass 0 //Calculate the transported mass between cell 1 and cell 2 //depending on the direction of the flow @Calc:(1, "dmass = if (dwater>0 , {1..diss} * dwater, {2..diss} * dwater)") @Calc:(1,"1..d = {1..d} - dmass") @Calc:(1,"2..d = {2..d} + dmass") } @Class: diffuse (name, D) { @Class: diffuse (, , diss) } @Class: diffuse (name, D, phase) This object is used to calculate diffusion of a solute driven by its concentration gradient according to Fick's law. It uses the concentrations in both cells, the distance between the cells (dx) in meter, and the surface area between the cells (A) in m2.
Wendy van Beinum, Hans Meeussen { @Var: 1.. 0 @Var: 2.. 0 @Var: 1..d 0 @Var: 2..d 0 @Var: .D //Calculate concentration gradient in (mol/m3)/m //dc=1000(C2-C1)/dx @Var: .dc 0 @Var: 1.dx 1 @Calc:(1, ".dc = 1000*({2..} - {1..})/1.dx") //Calculate the transported amount in mol/s //J= -D*dC*A @Var: .J 0 @Var: 2.A 0 @Calc:(1, ".J=(-{.D} * {.dc})*2.A") //Add and substract the transported amount from the mass changes in both cells @Calc:(1, "1..d = {1..d} - {.J}") @Calc:(1, "2..d = {2..d} + {.J}") } @Class: disperse (name, dispersivity, phase) This object is used to calculate dispersion of a solute driven by its concentration gradient, (material) dispersivity and pore-water velocity. It uses the concentrations in both cells, the distance between the cells (dx)in meter, the surface area between the cells (A) in m2, and the pore-water velocity (A*(J/porosity)) as input. J = in liters/s J/A = liter per s per m2 A is total contact area D = m porosity = water filled porosity //
Gijsbert Cirkel, Hans Meeussen 25-09-2012 { @Var: 1.. 0 @Var: 2.. 0 @Var: 1..d 0 @Var: 2..d 0 @Var: .D @Var: porosity 0.38 @Var: 1.porosity 1 @Var: 2.porosity 1 //**** The dx and porosity calculations are generic for all substances, and could bedefined outside this substance specific object. However, defining them here will not affect performance as the ORCHESTRA optimizer will remove the redundant calculations. *******// @Var: 1.dx 0.001 //distance between cell 1 and 2 (m) @Var: dx 1 // the local dx variable with default value @Calc:(1, "dx = 1.dx") // copy dx from cell @Calc:(1, "porosity = (1.porosity+2.porosity)/2") //**** //Calculate concentration gradient in (mol/m3)/m //dc=1000(C2-C1)/dx @Var: .dc 0 @Calc:(1, ".dc = 1000*({2..} - {1..})") //Calculate the transported amount in mol/s //J= -D*dC*A*pore-water velocity @Var: 1.A 1 @Var: .J 0 // .J = delta mass for this component, // dwater = amount of water in liter/sec // dwater/(A*porosity) = pore-water velocity @Calc:(1, ".J=(-porosity * {.D} * {.dc}/dx) * 1.A * (abs(dwater)/(1.A * porosity))") //porosity added GC 20-09-2012 // @Calc:(1, ".J=(-{.D} * {.dc}/dx) * abs(dwater)") //simplified version //Add and substract the transported amount from the mass changes per time in both cells @Calc:(1, "1..d = {1..d} - {.J}") @Calc:(1, "2..d = {2..d} + {.J}") } @Class: neutral_diffusion () This object class initialises the neutral diffusion calculations, and performs part of the calculations that are not component specific. It is used in combination with the n_diffus() object class in the diffusion.inp file.
Wendy van Beinum, Hans Meeussen { @Var: nd.sum1 0 //sum denominator @Var: nd.sum2 0 //sum nominator //**** Stage 1 **** @Calc:(1,"nd.sum1 = 0") // initialise sums with zero @Calc:(1,"nd.sum2 = 0") //**** Stage 2 **** @Var: nd.factor 0 @Calc:(2,"nd.factor = nd.sum1/nd.sum2") } @Class: n_diffus (name, charge, D) This object calculates diffusion for a single component and makes sure that there is no net transport of charge. It can be used after a neutral_diffusion() object is defined.
Wendy van Beinum, Hans Meeussen { @Var: 1..diss 0.0 @Var: 2..diss 0.0 @Var: 1..d 0.0 @Var: 2..d 0.0 @Var: .z @Var: 1..D //**** Stage 1 ***** // dC = C2-C1 @Var: .dc 0 @Calc:(1,".dc = 2..diss.sum - 1..diss") //Calculate contribution to sum D.z.dC @Calc:(1,"nd.sum1 = nd.sum1 + (1..D * .z * .dc) ") //Calculate contribution to sum D.z^2.C @Calc:(1,"nd.sum2 = nd.sum2 + (1..D * .z * .z * (1..diss+2..diss)/2)") //**** Stage 2 **** // Calculate J: J = A * D (-dC +z.C. nd.factor)/dx @Var: .J 0 @Var: 1.dx 1 @Var: 2.A 0 @Calc:(2,".J = 2.A * 1..D*(-.dc + (.z * 1..diss * nd.factor))/1.dx)") //Calculate contribution of J to delta mass of component in both cells @Calc:(2, "1..d = 1..d - .J") @Calc:(2, "2..d = 2..d + .J") } // The user interface objects (used by the graphical chemistry editor) @Class: primary_entity(name, unknownvalue, phase, equationvalue){ @Var: tolerance 1e-13 @globalvar: . // define this variable as global (node) variables + default value @Uneq2: unknown:(name:, .logact, delta:, 1e-6, type:, lin, step: , 1, default:, -12) equation:(name:, ., tol: , tolerance, default:, ) ) } @Class: primary_entity(name, unknownvalue){ @globalvar: .logact } @Class: primary_entity(name, unknown, unknownvalue, type, step, phase, equationvalue){ @Var: tolerance 1e-13 @GlobalVar: . // define this variable as global (node) variables + default value @Uneq2: unknown:(name:, , delta:, 1e-6, type:, , step:, , default:, ) equation:(name:, ., default:, , tol:, tolerance) } @Class: primary_entity(name, unknown, unknownvalue){ @globalvar: } @class: primary_entity_synonym:(synonym, name){ @GlobalVar: .tot 0 @GlobalVar: .diss 0 @GlobalVar: .ads 0 @GlobalVar: .min 0 @GlobalVar: .gas 0 @Synonym: .tot .tot @Synonym: .diss .diss @Synonym: .ads .ads @Synonym: .min .min @Synonym: .gas .gas } // An imineral is just a mineral that is selected as a primary entity // this creates a fixed activity @class: imineral(name){ @entity(, min, 1) } @Class: gas(name){ @entity(, gas, 1) } @Class: nicamodel(name, parentphase, concentration, donnanVolume){ // define the nica surface phase @phase(, , ) //-------------------------------------------------------------------------------------- // Add a donnan phase to this surface, which consist of a phase, linked to the // surface via volume/kg. The donnan phase is also an entity (act = boltzman factor, sum = charge balance) // unknown is activity, equation is charge balance at surface = 0 //-------------------------------------------------------------------------------------- @Var: donvol 1.5 // The Donnan volume + default value @Calc:(1,"donvol = ") @phase(_don, , donvol) @Var: _don.logact -1 @entity(_don, _don, 0) @Uneq2: unknown:(name:, _don.logact, delta:, 1e-6, type:, lin, step: , 1, default:, 0, iia:, true) equation:(name:, _don., tol: , 1e-8) //-------------------------------------------------------------------------------------------- } @Class: nicamodel(name, parentphase, concentration, donnanVolume, charge){ @nicamodel(, , , ) @Stage:(3,"_don. = {_don.} + ") // by adding a dummy phase and dummy proton entity to this surface e make that the net surface charge zero. // otherwise surface would be negative and consumes protons, or cations upon addition to solution // @Calc:(1,"donvol = ") // @phase(_dummy, , donvol) // @entity(dummyproton, _dummy, 1) // @link(dummyproton, H+, 0, ) } //---------------------------------------------------------------------------------------------- // A NICA site can be represented by a standard phase plus a standard entity. // The phase is linked to the total particle surface phase, by the site density sd (moles/kg ) //---------------------------------------------------------------------------------------------- @Class: nicasite(name, surface, donnan, p, sd) { // A nica site is a phase that is linked to its parent surface via its site densitiy @phase(, , ) // it is also an entity in its own phase @entity(, , 0) // by adding a dummy proton entity to this site we make that the net surface charge of nica site = zero. // otherwise adding nica-donnan surface is negative and consumes protons, or cations upon addition to solution @entity(dummyproton, , 1) @link(dummyproton, H+, 0, -1 ) // The amount of charge that is carried by the empty sites is added to the donnan mass (charge) balance @Stage:(3,". = {.} - ") // The value of the unknown C is used to calculate C^(p-1)/(1+C^p) which is the .act of the NICA site. // The following lines are the following formula: act = (unknown^(p-1))/(1+unknown^p) @Var: .unknown .1 @Var: .equation 0 @Stage:(1,".logact = log10(({.unknown}^(

-1))/(1+{.unknown}^

))") // Subtract the calculated sum(K^n* C^n) from the estimated sum, result should become zero // equation = calculated sum ((nic.act*C)/nic.act = C) - estimated sum (unknown) // equation is calculated logC - estimated logC, should become zero @Stage:(3,".equation = log10({.} / (10^{.logact})) - log10({.unknown})") @Uneq2: unknown:(name:, .unknown, delta:, 1e-3, type:, log, step:, 1, iia:, true) equation:(name:, .equation, tol: , 1e-6) } @Class: nicaspecies(name, site, ion, n, nH, logK) { @entity(, , "/") @Calc:(1, ".logact = *") @link(, , 1, "(/)") @link(, , , 1 ) } @Class: donnanspecies(name, phase, ion, charge) { @entity(, , 1) @calc:(1,".logact = {f.log})") @link(, , 1) @link(, , ) } @class: adsmodel(name, parent_phase ,concentration, type){ @(, ,) } @class: surfsite(model, name, density, coef){ //@Calc:(1, "_.logact = 0") @phase(_, , ) @entity(_, _, ) @Uneq2: unknown:(name:, _.logact, delta:, 1e-6, type:, lin, step: , 1, default:, -1, max:, 0) equation:(name:, _._, tol: , 1e-4, default:, 1) } @class: surfsite(model, name, density, coef, charge, plane){ //@Calc:(1, "_.logact = 0") @phase(_, , ) @entity(_, _, ) @link(_, _, 0, ) @Uneq2: unknown:(name:, _.logact, delta:, 1e-6, type:, lin, step: , 1, default:, -1, max:, 0) equation:(name:, _._, tol: , 1e-4, default:, 1) } @class: surfsite(model, name, density, coef, charge, plane, charge2, plane2){ //@Calc:(1, "_.logact = 0") @phase(_, , ) @entity(_, _, ) @link(_, _, 0, ) @link(_, _, 0, ) @Uneq2: unknown:(name:, _.logact, delta:, 1e-6, type:, lin, step: , 1, default:, -1, max:, 0) equation:(name:, _._, tol: , 1e-4, default:, 1) } @class: surfspecies(model, site, name, coef){ @entity(_, _, ) } @class: surfreaction(name, k, c1, n1){ @reaction(, , , ) } @class: surfreaction(name, k, c1, n1, c2, n2, ){ @reaction(, , , , , ) } @class: surfreaction(name, k, c1, n1, c2, n2, c3, n3){ @reaction(, , , , , , , ) } @class: surfreaction(name, k, c1, n1, c2, n2, c3, n3, c4, n4 ){ @reaction(, , , , , , , , , ) } @class: surfreaction(name, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5 ){ @reaction(, , , , , , , , , , , ) } @class: surfreaction(name, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6 ){ @reaction(, , , , , , , , , , , , , ) } @class: ddl_surface(name, parent_phase, surface_area){ @phase(, , ) @phase(_e) @link_phase(,_e, 96484.56) // The electrostatic layer is positioned in the _ee fase. This phase has // the sums (_e._e) is in the dimension coulomb/m2 after stage 4 @e_layer(_e, _e) @ddl(, _e) // @ddl_dzombak(, _e) @GlobalVar: .eq 0 // this defines eq as global variables with initial value 0. @Stage:(5, ".eq = {.ddl} + {_e._e}") @Uneq2: unknown:(name:, _e.logact, delta:, 1e-6, type:, lin, step: , 1, default:, 0, iia:, true) equation:(name:, .eq, tol: , 1e-6) } @Class: ddl(name, e_layer) { // Object to calculate diffuse double layer charge from given potential of electrostatic layer @Var: .ddl 0.0 @Var: .PF_RT 1 @GlobalVar: T 298.15 // Temperature with default value @Var: R 8.31441 @Var: F 96484.56 @Var: eps0 8.85419e-12 // dielectric constant vacuum Sposito p 229 @Var: eps1 78 // dielectric constant water Hiemstra p 46 @Calc:(1,".PF_RT = {.psi} * (F/(R*T)) ") @Var: .sum_of_charge 0 @Calc:(4,".ddl = if({.psi} > 0, {.ddl}, -{.ddl})") @Calc:(4,".ddl = if({.sum_of_charge} > 0, -sqrt(2000*eps0*eps1*R*T) * sqrt({.sum_of_charge}), 0 ))") // if the charge balance of the solution is not zero (which can happen during iterations, or partial system definition) // we assume that the charge unbalance is due to a missing -1 or +1 species and add its contribution to the sum @Calc: (4, ".sum_of_charge = .sum_of_charge + if (chargebalance>0, chargebalance * (exp(.PF_RT) -1), -chargebalance * (exp(-.PF_RT) -1) )" // we use the total concentrations of dissolved species for each distinct charge, as calculated by the ionic strength object @Calc:(4,".sum_of_charge = {f1.sum} * (exp(-1 * .PF_RT) -1) + {f-1.sum} * (exp(.PF_RT) -1) + {f2.sum} * (exp(-2 * .PF_RT) -1) + {f-2.sum} * (exp(2 * .PF_RT) -1) + {f3.sum} * (exp(-3 * .PF_RT) -1) + {f-3.sum} * (exp(3 * .PF_RT) -1)") } @Class: ddl_old (name, e_layer){ // This version of ddl object was used in first Orchestra versions @Var: T 298.15 @Var: R 8.31441 @Var: F 96484.56 @Var: .ddl 0 @Calc:(1,".ddl = -0.0587 * sqrt( I * (exp(-1 * {.psi} * (F/(R*T)) ) -1) + I * ( 1 * exp({.psi} * (F/(R*T)) ) -1)) )") @Calc:(1,".ddl = if({.psi} > 0, {.ddl}, -{.ddl})") } @Class: ddl_dzombak (name, e_layer){ @Var: F 96484.56 // This version of ddl model is used by the Generalized two layer model of Dzombak and Morel. @Var: .ddl 0.0 @Calc:(4,".ddl = -2.5*sqrt(I)*{.psi} ") //dzombak & morel p12 } //** // This object creates a kd entity and phase for a component that contains // an amount of mass equal to total dissolved concentration * kd // // usage in chemical input file: // @kd(K+, 4.8) // @class: kd(name, kd){ @phase(_kd, ads, "(*.diss)") // concentration in ads phase is dissolved concentration * kd @entity(_kd,_kd, 1) @link(_kd, , 0, 1) } // This object calculates the amount of adsorbed protons necessary to agree with measured ANC // usage in chemical input file: // @ANC("4e-3*pH") // @class: ANC(formula){ // here we need a formula that calculates the total amount of adsorbed protons from the pH. // subtract the actual adsorbed amount of protons from this, so the effective total equals ANC @phase(ANC, ads, - H+.ads.sum") @entity(ANC, ANC, 1) @link(_kd, , 0, 1) } @Class: balance_charge(cation,anion) // usage: @balance_charge(Na+, Cl-) // This object balances charge by adding either the amount of a cation or an anion, whichever is appropriate. // The charge_min_act variable represents the lower limit value for the ion that is not corrected. // Make sure that for each of these ions a constant activity is chosen in the primary entities settings // so no iteration on their individual mass balances is carried out. // Hans Meeussen. July 2012 { @Var: chargebalance 0 @Var: charge_unknown -1e-3 // the unknown that is varied until chargebalance is zero. @Var: charge_min_act 1e-3 // the concentration of the cat/anion that is not corrected. @Stage: (1, ".logact = if(charge_unknown>0, log10(charge_min_act + charge_unknown), log10(charge_min_act))") @Stage: (1, ".logact = if(charge_unknown<0, log10(charge_min_act - charge_unknown), log10(charge_min_act))") @Uneq2: unknown:(name:, charge_unknown, delta:, 1e-3, type:, lin, step: , .1) equation:(name:, chargebalance, tol: , 1e-6) } @Class: balance_charge(ion) // usage: @balance_charge(Cl-) // This object balances charge by adding the amount of an ion. // Make sure that for this ion a constant activity is chosen in the primary entities settings // Hans Meeussen. July 2012 { @GlobalVar: chargebalance 0 @Uneq2: unknown:(name:, .logact, delta:, 1e-3, type:, lin, step: , .1) equation:(name:, chargebalance, tol: , 1e-6) } // the radioactive decay objects // Hans Meeussen, NRG, September, 2010 @class: radionuclide(name, value){ @var: .halflife // half life in seconds @globalvar: .tot 1 // The total amount of mother isotope @globalvar: .d 0 // The delta amount of mother isotope @calc:(1, ".d = {.d} - (log(2)/{.halflife} * if({.tot} > 1e-20, {.tot}, 0))") } @class: halflife(name, value){ @radionuclide(, ) } @class: decay(mother, frac, daughter) // mother component, fraction daughter component { @globalvar: .d 0 // The delta amount of daughter product // here we limit the the decay at a minimum total concentration of 1e-20 M to prevent zero concentrations @calc:(1, ".d = {.d} + * (log(2)/{.halflife} * if({.tot} > 1e-20, {.tot}, 0))") } // Merge isotope concentrations into total element concentrations @class: prechem_perphase(phase, isotope, component){ @globalVar: . 0 @Calc: (1, ". = {.} + {.}") } // Split element concentrations back into individual isotope concentrations @class: postchem_perphase(phase, isotope, component){ @globalVar: . 0 @globalVar: . 0 @Calc: (1, ". = {.} * {.frac}") }