plot.aquaenv

Description

PUBLIC function: high level plot function for objects of class aquaenv

Arguments

x

object of class aquaenv

xval

only valid if bjerrum=FALSE: a vector of the (maximal) length of the elements of aquaenv against which they are to be plotted

what

a list of names of the elements of aquaenv that are to be plotted, if not supplied and bjerrum=FALSE and cumulative=FALSE: all elements are plotted, if not supplied and bjerrum=TRUE then what is set to be c("CO2", "HCO3", "CO3", "BOH3", "BOH4", "OH", "H3PO4", "H2PO4", "HPO4", "PO4", "SiOH4", "SiOOH3", "SiO2OH2", "H2S", "HS", "S2min", "NH4", "NH3", "H2SO4", "HSO4", "SO4", "HF", "F", "HNO3", "NO3", "HNO2", "NO2"), needs to be supplied for cumulative=TRUE

bjerrum

flag: TRUE = a bjerrum plot is done (by calling bjerrumplot)

cumulative

flag: TRUE = a cumulative plot is done (by calling cumulativeplot)

newdevice

flag: if TRUE, new plot device is opened

setpar

flag: if TRUE parameters are set with the function par

xlab

x axis label

log

only valif if bjerrum=TRUE: should the plot be on a logarithmic y axis?

total

only valid if cumulative=TRUE: should the sum of all elements specified in what be plotted as well?

device

the device to plot on; default: "x11" (can also be "eps" or "pdf")

filename

filename to be used if "eps" or "pdf" is selected for device

size

the size of the plot device; default: 12 (width) by 10 (height) inches

ylim

standard plot parameter; if not supplied it will be calculated by range() of the elements to plot

lwd

standard plot parameter; width of the lines in the plot

mgp

standard plot parameter; default: axis title on line 1.8, axis labels on line 0.5, axis on line 0

mar

standard plot parameter; default: margin of 3 lines bottom and left and 0.5 lines top and right

oma

standard plot parameter; default: no outer margin

palette

only valid if bjerrum=TRUE or cumulative=TRUE: a vector of colors to use in the plot (either numbers or names given in colors())

legendposition

only valid if bjerrum=TRUE or cumulative=TRUE: position of the legend

legendinset

only valid if bjerrum=TRUE or cumulative=TRUE: standard legend parameter inset

legendlwd

only valid if bjerrum=TRUE or cumulative=TRUE: standard legend parameter lwd: line width of lines in legend

bg

only valid if bjerrum=TRUE or cumulative=TRUE: standard legend parameter: default background color: white

y.intersp

standard legend parameter; if cumulative=TRUE then default: 1.2 lines space between the lines in the legend

...

further arguments are passed on to the plot function

Details

Top level generic usage is

1
2
3
plot.aquaenv(x, xval, what=NULL, bjerrum=FALSE,
             cumulative=FALSE, newdevice=TRUE, setpar=TRUE,
	     device="x11", ...)

Generic usages for standard plotting are

1
plot.aquaenv(x, xval, ...)
1
plot.aquaenv(x, xval, what, mfrow=c(1,1), size=c(7,7), ...)

Generic usage for creating a bjerrum plot is

1
2
3
4
5
6
7
plot.aquaenv(x, what, log=FALSE, palette=NULL,
             device="x11", filename="aquaenv",
             size=c(12,10), ylim=NULL, lwd=2, xlab="free scale pH",
	     mgp=c(1.8, 0.5, 0), mar=c(3,3,0.5,0.5), oma=c(0,0,0,0),
	     legendposition="bottomleft", legendinset=0.05,
	     legendlwd=4, bg="white", newdevice=TRUE, setpar=TRUE,
	     device="x11",...)

Generic usage for creating a cumulative plot is

1
2
3
4
5
6
plot.aquaenv(x, xval, what, total=TRUE, palette=NULL,
             device="x11", filename="aquaenv", size=c(12,10), ylim=NULL,
             lwd=2, mgp=c(1.8, 0.5, 0), mar=c(3,3,0.5,0.5), oma=c(0,0,0,0),
             legendposition="bottomleft", legendinset=0.05, legendlwd=4,
             bg="white", y.intersp=1.2, newdevice=TRUE, setpar=TRUE,
	     device="x11",...)

Author(s)

Andreas F. Hofmann. Maintained by Karline Soetaert (Karline.Soetaert@nioz.nl).

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
## Not run: 
### 0
#####
A <- aquaenv(35, 15, SumCO2=0.003, TA=seq(0.001,0.004, 0.0001))
plot(A, xval=A$TA, xlab="[TA]/(mol/kg-soln)")
plot(A, what=c("CO2", "HCO3", "CO3"), bjerrum=TRUE, log=TRUE)
plot(A, xval=A$TA, xlab="[TA]/(mol/kg-soln)", what=c("CO2", "HCO3", "CO3"),
     cumulative=TRUE, ylab="mol/kg-soln", ylim=c(0,0.0031))


### 1
#####

SumCO2 <- 0.0020
pH     <- 8

S      <- 30
t      <- 1:15
p      <- 10
ae <- aquaenv(S, t, p, SumCO2=SumCO2, pH=pH, revelle=TRUE, dsa=TRUE)
plot(ae, xval=t, xlab="T/(deg C)", newdevice=FALSE)



### 2
#####
S <- 35
t <- 15

SumCO2 <- 0.003500
SumNH4 <- 0.000020

mass_sample  <- 0.01 # the mass of the sample solution in kg
mass_titrant <- 0.02 # the total mass of the added titrant solution in
                     # kg
conc_titrant <- 0.01 # the concentration of the titrant solution in
                     # mol/kg-soln
S_titrant    <- 0.5  # the salinity of the titrant solution (the
                     # salinity of a solution with a ionic strength of
                     # 0.01 according to: I = (19.924 S) / (1000 - 1.005S)
steps        <- 50   # the amount of steps the mass of titrant is added
                     # in
type         <- "HCl"

pHstart <- 11.3


ae <- titration(aquaenv(S=S, t=t, SumCO2=SumCO2, SumNH4=SumNH4,
                pH=pHstart), mass_sample, mass_titrant, conc_titrant,
                S_titrant, steps, type)


# plotting everything
plot(ae, xval=ae$delta_mass_titrant, xlab="HCl solution added [kg]",
mfrow=c(10,10))


# plotting selectively
size  <- c(12,8) #inches
mfrow <- c(4,4)
what  <- c("TA", "pH", "CO2", "HCO3", "CO3", "BOH3", "BOH4", "OH",
           "NH4", "NH3", "H2SO4", "HSO4", "SO4", "HF", "F", "pCO2")


plot(ae, xval=ae$delta_mass_titrant, xlab="HCl solution added [kg]",
     what=what, size=size, mfrow=mfrow)

plot(ae, xval=ae$pH, xlab="free scale pH", what=what, size=size,
     mfrow=mfrow)


# different x values
plot(ae, xval=ae$delta_conc_titrant, xlab="[HCl] offset added
     [mol/kg-soln]", what=what, size=size, mfrow=mfrow)

plot(ae, xval=ae$delta_moles_titrant, xlab="HCl added [mol]", what=what,
     size=size, mfrow=mfrow, newdevice=FALSE)


# bjerrum plots
plot(ae, bjerrum=TRUE)

what  <- c("CO2", "HCO3", "CO3")
plot(ae, what=what, bjerrum=TRUE)
plot(ae, what=what, bjerrum=TRUE, lwd=4, palette=c("cyan", "magenta",
     "yellow"), bg="gray", legendinset=0.1, legendposition="topleft")



what  <- c("CO2", "HCO3", "CO3", "BOH3", "BOH4", "OH", "NH4", "NH3",
           "H2SO4", "HSO4", "SO4", "HF", "F")

plot(ae, what=what, bjerrum=TRUE, log=TRUE, newdevice=FALSE)
plot(ae, what=what, bjerrum=TRUE, log=TRUE, ylim=c(-6,-1),
     legendinset=0, lwd=3, palette=c(1,3,4,5,6,colors()[seq(100,250,6)]))


### 3
#####
parameters <- list(             
    t           = 15        , # degrees C
    S           = 35        , # psu       

    SumCO2_t0   = 0.002     , # mol/kg-soln  (comparable to Wang2005)
    TA_t0       = 0.0022    , # mol/kg-soln  (comparable to Millero1998)

    kc          = 0.5       , # 1/d	         proportionality factor
                              #                  for air-water exchange
    kp          = 0.000001  , # mol/(kg-soln*d)	 max rate of calcium
                              #                  carbonate precipitation
    n           = 2.0       , # -                exponent for kinetic
                              #                  rate law of precipitation
 
    modeltime   = 20        , # d              duration of the model
    outputsteps = 100         #                number of outputsteps
                   )

boxmodel <- function(timestep, currentstate, parameters)
{
  with (
        as.list(c(currentstate,parameters)),
        {        
          ae    <- aquaenv(S=S, t=t, SumCO2=SumCO2, pH=-log10(H), SumSiOH4=0, 
                           SumBOH3=0, SumH2SO4=0, SumHF=0, dsa=TRUE)
                   
          Rc    <- kc * ((ae$CO2_sat) - (ae$CO2)) 
          Rp    <- kp * (1-ae$omega_calcite)^n               

          dSumCO2 <- Rc - Rp

          dHRc    <- (      -(ae$dTAdSumCO2*Rc   ))/ae$dTAdH
          dHRp    <- (-2*Rp -(ae$dTAdSumCO2*(-Rp)))/ae$dTAdH
          dH      <- dHRc + dHRp
          
          ratesofchanges <- c(dSumCO2, dH)
          
          processrates   <- c(Rc=Rc, Rp=Rp)
          outputvars     <- c(dHRc=dHRc, dHRp=dHRp)
          
          return(list(ratesofchanges, list(processrates, outputvars, ae)))
        }
        )
}

with (as.list(parameters),
      {
        aetmp <- aquaenv(S=S, t=t, SumCO2=SumCO2_t0,
                         TA=TA_t0, SumSiOH4=0, SumBOH3=0,
                         SumH2SO4=0, SumHF=0)
        H_t0  <- 10^(-aetmp$pH)
        
        initialstate <<- c(SumCO2=SumCO2_t0, H=H_t0)
        times        <<- seq(0,modeltime,(modeltime/outputsteps))       
        output       <<- as.data.frame(vode(initialstate,times,
                                    boxmodel,parameters, hmax=1))
      })

what   <- c("SumCO2", "TA", "Rc", "Rp",
            "omega_calcite", "pH", "dHRc", "dHRp")
plot(aquaenv(ae=output, from.data.frame=TRUE), xval=output$time,
     xlab="time/d", mfrow=c(3,3), size=c(15,10), what=what) 

what <- c("dHRc", "dHRp")
plot(aquaenv(ae=output, from.data.frame=TRUE), xval=output$time,
     xlab="time/d", what=what, ylab="mol-H/(kg-soln*d)",
     legendposition="topright", cumulative=TRUE) 

## End(Not run)