transfer(nsteps = 500, dmode = "coupled", devmax = 0.1, plot = NULL,
ibalance = 1, fmode = "one", buffers = NULL, alphamax = -2,
alphastart = -10, T = 25, P = "Psat", do.title = TRUE, beta = 0)
draw.transfer(t, ylim = c(-10, 1), ylimbasis = c(-12, -2),
logprogress = FALSE)
feldspar(which = "closed", plot.it = FALSE)
apc(which = "open", basis = "CO2", plot.it = FALSE)
transfer
.transfer
returns a list containing information about the conditions at each step: basis
, data frame of the logarithms of activities of basis species, species
, data frame of the logarithms of activities (moles for solids) of species, alphas
, numeric vector of the values of the destruction exponent, dmodes
, character vector of the destruction mode, istables
, numeric vector of the index of the most stable product, myaffs
, list of the affinities of the formation reactions of species, didwork
, logical vector indicating whether the steps succeeded or failed.species
. At each step, a small amount ($10^{\alpha}$) of the starting composition is provisionally reacted and a relatively more stable product may be formed. The amount of product formed is such that the activity of the primary conservant (the basis species given in ibalance
) is not changed. The changes in the activities of the other basis species are calculated, and the process is iterated until nsteps
is reached or the value of $\alpha$ is driven to a very low value (logpresent
, which is a constant set in the code to $-50$).
If at a given step the most stable product is different from the one before, either the previous products are ignored (for dmode
equal to none, i.e. an open system) or the reaction of the starting material is coupled to that of the existing products (for dmode
equal to coupled, i.e. a closed system) through a secondary conservation constraint. The basis species that are candidates for the secondary conservation are identified in iplot
.
The initial value of alpha is given by alphastart
. After successful steps, the function increases the value of alpha by $1$, and after failed steps decreases the value of alpha by $1$. One condition that can lead to a failed step is that the logarithm of activity of any basis species changes by more than devmax
. Therefore, throughout the simulation the value of $\alpha$ dynamically adjusts based on the devmax
set by the user.
buffers
is a list with elements basis
indicating the basis species to be buffered and buffer
naming the buffers to use for that basis species. If this argument is given, at each step the activity of the basis species in the buffer is calculated. The difference between this activity and the current activity of the basis species in the system is then multiplied by 10 raised to the ( alpha
+ beta
) and this quantity added to the current activity of the basis species in the system. As a result, the value of beta
modifies the strength of the buffer relative to the incremental reaction progress.
draw.transfer
is used to plot the logarithms of activities of basis species, and logarithms of activities (moles for solid species, molalities for aqueous species) of the minerals or proteins as a function of reaction progress, or logarithm of reaction progress if logprogress
is set to TRUE
. The y-limits of the plots can be set using ylim
and ylimbasis
.
feldspar
and apc
encode examples for feldspar weathering and reactions among proteins in the anaphase-promoting complex of yeast.
Steinmann, P., Lichtner, P. C. and Shotyk, W. (1994) Reaction path approach to mineral weathering reactions. Clay Clay Min. 42, 197--206.
data(thermo)
## react potassium feldspar in a closed system
## after Steinmann et al., 1994 and Helgeson et al., 1969
basis(c("Al+3", "H4SiO4", "K+", "H2O", "H+", "O2"), c(0, -6, -6, 0, 0, 0))
species(c("k-feldspar", "muscovite", "pyrophyllite", "kaolinite", "gibbsite"))
a <- affinity(H4SiO4=c(-6, -2), "K+"=c(-3, 8))
diagram(a, fill="heat")
basis("pH", 4)
species(1:5, c(-4, rep(-999, 4)))
tr <- transfer(550, dmode="coupled", plot=c(2, 3), devmax=0.2)
# plot the output from transfer
draw.transfer(tr)
# reset the plot layout
layout(matrix(1))
## can also run the calculation above with
#feldspar("closed")
## or an example for proteins with
#apc("closed")
Run the code above in your browser using DataLab