IREACT=100

[jtwhite79]

Hey guys - I'm having some trouble getting a decent mass balance from a model if use RCT and the zeroth order reaction (IREACT=100). First order good, no reaction good, production ( negative rc1) good, but decay (positive rc1) bad. I created a simple 3-cell model that demonstrates the issue. Please let me know if this is some sort of user/input error.
rct_trouble.zip

[emorway-usgs]

Hedeff Essaid (Menlo office) reported an issue with zero order as well (and sent a test model demonstrating the problem). Vivek (cc'd) is looking in to it, we'll get a fix out soon. I'll append her email to your issue.

Eric Morway Research Hydrologist Nevada Water Science Center U.S. Geological Survey 2730 N. Deer Run Rd. Carson City, NV 89701 (775) 887-7668 orcid: 0000-0002-8553-6140 http://orcid.org/0000-0002-8553-6140

On Thu, Jan 4, 2018 at 5:53 PM, jtwhite79 notifications@github.com wrote:

Hey guys - I'm having some trouble getting a decent mass balance from a model if use RCT and the zeroth order reaction (IREACT=100). First order good, no reaction good, production ( negative rc1) good, but decay (positive rc1) bad. I created a simple 3-cell model that demonstrates the issue. Please let me know if this is some sort of user/input error. rct_trouble.zip https://github.com/langevin-usgs/mt3d-usgs/files/1605293/rct_trouble.zip

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/langevin-usgs/mt3d-usgs/issues/26, or mute the thread https://github.com/notifications/unsubscribe-auth/ADFi4PzVPPHvxeIqLfQcFGzHuNJEpb2Jks5tHYCJgaJpZM4RT8Dl .

[emorway-usgs]

From Hedeff on 10/17/2017: "In one set of simulations (Model#3 in the poster) I am fitting zero-order reaction rates for multiple solutes. One of the solutes is nitrate which is generally completely consumed in the lakebed sediments which are mainly anaerobic. I noticed that zero-order decay led to negative nitrate concentrations. Earlier simulations with first-order decay did not because decay decreased to zero as the concentrations decreased to zero. However, with zero-order decay it seems that maybe the code is not shutting off the reaction term when nitrate disappears. My next simulations will probably incorporate Monod Kinetics which should take care of this problem, but I thought you might want to know about this issue."

And on 10/19/2017: "The model set-up is very simple. I am simulating transport through a 1D vertical column that is 10 cm in length. It represents the upper 10 cm of lakebed sediment below Upper Klamath Lake. I am not using the LAK package. I have been trying out several different conceptual models for transport through the lakebed. The GSA poster outlines three of the models I have run so far. The zero-order rate case is shown as Model #3 with results plotted at the bottom of the central panel of the poster. There is no groundwater flow through the lakebed in this case. All transport is driven by molecular diffusion and production/decay of solutes in the lakebed. This is the common conceptual model that limnologists use for lakebeds. The point of my study is to include the effects of groundwater advection on lakebed exchange with the overlaying lake water. Nitrate is removed by zero-order decay in this scenario and because nitrate concentrations are very low to start with, the concentrations at depths greater than 3 cm are negative during some periods of the simulation. You can see this in the output file NO3MTTVD4n.obs by scrolling to the right and then scrolling down to the middle of the simulation period.

I have put a zip file (ZOinputfiles.zip) of all of the input files, the poster describing the work and the NO3MTTVD4n.obs output file on google drive and will share the file with you." ZOinputfiles.zip

[jtwhite79]

Thanks Eric - I tried to debug, but there's a lot going on in there (at least for my small brain).

[emorway-usgs]

Just tested Hedeff's model, and all results appear to be matching with MT3DMS, except where MT3DMS gives negative concentrations for NO3 and DO. @jtwhite79 mentioned that a cursory review of results from the NZ model also look good. Fixed with commits bcf13e4 and c5bbfc2

From Hedeff's model: DO, cell 67:

DO, cell 111:

NO3, cell 67:

NO3, cell 111:

[emorway-usgs]

R script that was used to generate the plots:

Fe_mtusgs  <- read.table('c:/temp/zoinputfiles/mt3d-usgs/femttvd7n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
DO_mtusgs  <- read.table('c:/temp/zoinputfiles/mt3d-usgs/DOmttvd9n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
DOC_mtusgs <- read.table('c:/temp/zoinputfiles/mt3d-usgs/DOCmttvd5n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
As_mtusgs  <- read.table('c:/temp/zoinputfiles/mt3d-usgs/Asmttvd8n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
Mn_mtusgs  <- read.table('c:/temp/zoinputfiles/mt3d-usgs/Mnmttvd6n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
NH4_mtusgs <- read.table('c:/temp/zoinputfiles/mt3d-usgs/NH4mttvd3n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
NO3_mtusgs <- read.table('c:/temp/zoinputfiles/mt3d-usgs/NO3mttvd4n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
P_mtusgs   <- read.table('c:/temp/zoinputfiles/mt3d-usgs/Pmttvd2n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))

Fe_mtms  <- read.table('c:/temp/zoinputfiles/mt3dms/femttvd7n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
DO_mtms  <- read.table('c:/temp/zoinputfiles/mt3dms/DOmttvd9n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
DOC_mtms <- read.table('c:/temp/zoinputfiles/mt3dms/DOCmttvd5n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
As_mtms  <- read.table('c:/temp/zoinputfiles/mt3dms/Asmttvd8n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
Mn_mtms  <- read.table('c:/temp/zoinputfiles/mt3dms/Mnmttvd6n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
NH4_mtms <- read.table('c:/temp/zoinputfiles/mt3dms/NH4mttvd3n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
NO3_mtms <- read.table('c:/temp/zoinputfiles/mt3dms/NO3mttvd4n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))
P_mtms   <- read.table('c:/temp/zoinputfiles/mt3dms/Pmttvd2n.obs', skip=2, header=FALSE, col.names = c('step','time','c1','c2','c3','c21','c41','c67','c111','c200'))

# Cell 1
# -------

# Fe
png('c:/temp/zoinputfiles/fe_c1.png', height=500, width=800, res=130)
plot(Fe_mtms$time[seq(1,nrow(Fe_mtms),by=10)], Fe_mtms$c1[seq(1,nrow(Fe_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='Fe Concentration', ylim=c(0,0.30))
points(Fe_mtusgs$time[seq(1,nrow(Fe_mtusgs),by=50)], Fe_mtusgs$c1[seq(1,nrow(Fe_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# DO
png('c:/temp/zoinputfiles/do_c1.png', height=500, width=800, res=130)
plot(DO_mtms$time[seq(1,nrow(DO_mtms),by=10)], DO_mtms$c1[seq(1,nrow(DO_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='DO Concentration', ylim=c(0,11.5))
points(DO_mtusgs$time[seq(1,nrow(DO_mtusgs),by=50)], DO_mtusgs$c1[seq(1,nrow(DO_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# DOC
png('c:/temp/zoinputfiles/doc_c1.png', height=500, width=800, res=130)
plot(DOC_mtms$time[seq(1,nrow(DOC_mtms),by=10)], DOC_mtms$c1[seq(1,nrow(DOC_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='DOC Concentration', ylim=c(0,6.75))
points(DOC_mtusgs$time[seq(1,nrow(DOC_mtusgs),by=50)], DOC_mtusgs$c1[seq(1,nrow(DOC_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topleft', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# As
png('c:/temp/zoinputfiles/as_c1.png', height=500, width=800, res=130)
plot(As_mtms$time[seq(1,nrow(As_mtms),by=10)], As_mtms$c1[seq(1,nrow(As_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='As Concentration', ylim=c(0,0.0125))
points(As_mtusgs$time[seq(1,nrow(As_mtusgs),by=50)], As_mtusgs$c1[seq(1,nrow(As_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topleft', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# Mn
png('c:/temp/zoinputfiles/mn_c1.png', height=500, width=800, res=130)
plot(Mn_mtms$time[seq(1,nrow(Mn_mtms),by=10)], Mn_mtms$c1[seq(1,nrow(Mn_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='Mn Concentration', ylim=c(0,0.11))
points(Mn_mtusgs$time[seq(1,nrow(Mn_mtusgs),by=50)], Mn_mtusgs$c1[seq(1,nrow(Mn_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# NH4
png('c:/temp/zoinputfiles/nh4_c1.png', height=500, width=800, res=130)
plot(NH4_mtms$time[seq(1,nrow(NH4_mtms),by=10)], NH4_mtms$c1[seq(1,nrow(NH4_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='NH4 Concentration', ylim=c(0,2))
points(NH4_mtusgs$time[seq(1,nrow(NH4_mtusgs),by=50)], NH4_mtusgs$c1[seq(1,nrow(NH4_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# NO3
png('c:/temp/zoinputfiles/no3_c1.png', height=500, width=800, res=130)
plot(NO3_mtms$time[seq(1,nrow(NO3_mtms),by=10)], NO3_mtms$c1[seq(1,nrow(NO3_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='NO3 Concentration', ylim=c(0,0.1))
points(NO3_mtusgs$time[seq(1,nrow(NO3_mtusgs),by=50)], NO3_mtusgs$c1[seq(1,nrow(NO3_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# P
png('c:/temp/zoinputfiles/p_c1.png', height=500, width=800, res=130)
plot(P_mtms$time[seq(1,nrow(P_mtms),by=10)], P_mtms$c1[seq(1,nrow(P_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='P Concentration', ylim=c(0,0.45))
points(P_mtusgs$time[seq(1,nrow(P_mtusgs),by=50)], P_mtusgs$c1[seq(1,nrow(P_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# Cell 67
# -------

# Fe
png('c:/temp/zoinputfiles/fe_c67.png', height=500, width=800, res=130)
plot(Fe_mtms$time[seq(1,nrow(Fe_mtms),by=10)], Fe_mtms$c67[seq(1,nrow(Fe_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='Fe Concentration', ylim=c(0,1.30))
points(Fe_mtusgs$time[seq(1,nrow(Fe_mtusgs),by=50)], Fe_mtusgs$c67[seq(1,nrow(Fe_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# DO
png('c:/temp/zoinputfiles/do_c67.png', height=500, width=800, res=130)
plot(DO_mtms$time[seq(1,nrow(DO_mtms),by=10)], DO_mtms$c67[seq(1,nrow(DO_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='DO Concentration', ylim=c(-3,4))
points(DO_mtusgs$time[seq(1,nrow(DO_mtusgs),by=50)], DO_mtusgs$c67[seq(1,nrow(DO_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'),  bty='n', bg='white')
abline(h=0)
dev.off()

# DOC
png('c:/temp/zoinputfiles/doc_c67.png', height=500, width=800, res=130)
plot(DOC_mtms$time[seq(1,nrow(DOC_mtms),by=10)], DOC_mtms$c67[seq(1,nrow(DOC_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='DOC Concentration', ylim=c(0,6.75))
points(DOC_mtusgs$time[seq(1,nrow(DOC_mtusgs),by=50)], DOC_mtusgs$c67[seq(1,nrow(DOC_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topleft', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'),  bty='n', bg='white')
dev.off()

# As
png('c:/temp/zoinputfiles/as_c67.png', height=500, width=800, res=130)
plot(As_mtms$time[seq(1,nrow(As_mtms),by=10)], As_mtms$c67[seq(1,nrow(As_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='As Concentration', ylim=c(0,0.0125))
points(As_mtusgs$time[seq(1,nrow(As_mtusgs),by=50)], As_mtusgs$c67[seq(1,nrow(As_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'),  bty='n', bg='white')
dev.off()

# Mn
png('c:/temp/zoinputfiles/mn_c67.png', height=500, width=800, res=130)
plot(Mn_mtms$time[seq(1,nrow(Mn_mtms),by=10)], Mn_mtms$c67[seq(1,nrow(Mn_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='Mn Concentration', ylim=c(0,0.25))
points(Mn_mtusgs$time[seq(1,nrow(Mn_mtusgs),by=50)], Mn_mtusgs$c67[seq(1,nrow(Mn_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'),  bty='n', bg='white')
dev.off()

# NH4
png('c:/temp/zoinputfiles/nh4_c67.png', height=500, width=800, res=130)
plot(NH4_mtms$time[seq(1,nrow(NH4_mtms),by=10)], NH4_mtms$c67[seq(1,nrow(NH4_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='NH4 Concentration', ylim=c(0,5))
points(NH4_mtusgs$time[seq(1,nrow(NH4_mtusgs),by=50)], NH4_mtusgs$c67[seq(1,nrow(NH4_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('bottomright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# NO3
png('c:/temp/zoinputfiles/no3_c67.png', height=500, width=800, res=130)
plot(NO3_mtms$time[seq(1,nrow(NO3_mtms),by=10)], NO3_mtms$c67[seq(1,nrow(NO3_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='NO3 Concentration', ylim=c(-0.05,0.075))
points(NO3_mtusgs$time[seq(1,nrow(NO3_mtusgs),by=50)], NO3_mtusgs$c67[seq(1,nrow(NO3_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
abline(h=0)
dev.off()

# P
png('c:/temp/zoinputfiles/p_c67.png', height=500, width=800, res=130)
plot(P_mtms$time[seq(1,nrow(P_mtms),by=10)], P_mtms$c67[seq(1,nrow(P_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='P Concentration', ylim=c(0,1.3))
points(P_mtusgs$time[seq(1,nrow(P_mtusgs),by=50)], P_mtusgs$c67[seq(1,nrow(P_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('bottomright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# Cell 111
# --------

# Fe
png('c:/temp/zoinputfiles/fe_c111.png', height=500, width=800, res=130)
plot(Fe_mtms$time[seq(1,nrow(Fe_mtms),by=10)], Fe_mtms$c111[seq(1,nrow(Fe_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='Fe Concentration', ylim=c(0,1.50))
points(Fe_mtusgs$time[seq(1,nrow(Fe_mtusgs),by=50)], Fe_mtusgs$c111[seq(1,nrow(Fe_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# DO
png('c:/temp/zoinputfiles/do_c111.png', height=500, width=800, res=130)
plot(DO_mtms$time[seq(1,nrow(DO_mtms),by=10)], DO_mtms$c111[seq(1,nrow(DO_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='DO Concentration', ylim=c(-4,1))
points(DO_mtusgs$time[seq(1,nrow(DO_mtusgs),by=50)], DO_mtusgs$c111[seq(1,nrow(DO_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('bottomleft', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'),  bty='n', bg='white')
abline(h=0)
dev.off()

# DOC
png('c:/temp/zoinputfiles/doc_c111.png', height=500, width=800, res=130)
plot(DOC_mtms$time[seq(1,nrow(DOC_mtms),by=10)], DOC_mtms$c111[seq(1,nrow(DOC_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='DOC Concentration', ylim=c(0,6.75))
points(DOC_mtusgs$time[seq(1,nrow(DOC_mtusgs),by=50)], DOC_mtusgs$c111[seq(1,nrow(DOC_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topleft', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'),  bty='n', bg='white')
dev.off()

# As
png('c:/temp/zoinputfiles/as_c111.png', height=500, width=800, res=130)
plot(As_mtms$time[seq(1,nrow(As_mtms),by=10)], As_mtms$c111[seq(1,nrow(As_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='As Concentration', ylim=c(0,0.0125))
points(As_mtusgs$time[seq(1,nrow(As_mtusgs),by=50)], As_mtusgs$c111[seq(1,nrow(As_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'),  bty='n', bg='white')
dev.off()

# Mn
png('c:/temp/zoinputfiles/mn_c111.png', height=500, width=800, res=130)
plot(Mn_mtms$time[seq(1,nrow(Mn_mtms),by=10)], Mn_mtms$c111[seq(1,nrow(Mn_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='Mn Concentration', ylim=c(0,0.35))
points(Mn_mtusgs$time[seq(1,nrow(Mn_mtusgs),by=50)], Mn_mtusgs$c111[seq(1,nrow(Mn_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'),  bty='n', bg='white')
dev.off()

# NH4
png('c:/temp/zoinputfiles/nh4_c111.png', height=500, width=800, res=130)
plot(NH4_mtms$time[seq(1,nrow(NH4_mtms),by=10)], NH4_mtms$c111[seq(1,nrow(NH4_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='NH4 Concentration', ylim=c(0,7))
points(NH4_mtusgs$time[seq(1,nrow(NH4_mtusgs),by=50)], NH4_mtusgs$c111[seq(1,nrow(NH4_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# NO3
png('c:/temp/zoinputfiles/no3_c111.png', height=500, width=800, res=130)
plot(NO3_mtms$time[seq(1,nrow(NO3_mtms),by=10)], NO3_mtms$c111[seq(1,nrow(NO3_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='NO3 Concentration', ylim=c(-0.1,0.1))
points(NO3_mtusgs$time[seq(1,nrow(NO3_mtusgs),by=50)], NO3_mtusgs$c111[seq(1,nrow(NO3_mtusgs),by=50)], typ='p', pch=16, col='red')
abline(h=0)
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# P
png('c:/temp/zoinputfiles/p_c111.png', height=500, width=800, res=130)
plot(P_mtms$time[seq(1,nrow(P_mtms),by=10)], P_mtms$c111[seq(1,nrow(P_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='P Concentration', ylim=c(0,1.6))
points(P_mtusgs$time[seq(1,nrow(P_mtusgs),by=50)], P_mtusgs$c111[seq(1,nrow(P_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('bottomleft', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# Cell 200
# --------

# Fe
png('c:/temp/zoinputfiles/fe_c200.png', height=500, width=800, res=130)
plot(Fe_mtms$time[seq(1,nrow(Fe_mtms),by=10)], Fe_mtms$c200[seq(1,nrow(Fe_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='Fe Concentration', ylim=c(0,1.50))
points(Fe_mtusgs$time[seq(1,nrow(Fe_mtusgs),by=50)], Fe_mtusgs$c200[seq(1,nrow(Fe_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# DO
png('c:/temp/zoinputfiles/do_c200.png', height=500, width=800, res=130)
plot(DO_mtms$time[seq(1,nrow(DO_mtms),by=10)], DO_mtms$c200[seq(1,nrow(DO_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='DO Concentration', ylim=c(-1,1))
points(DO_mtusgs$time[seq(1,nrow(DO_mtusgs),by=50)], DO_mtusgs$c200[seq(1,nrow(DO_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'),  bty='n', bg='white')
abline(h=0)
dev.off()

# DOC
png('c:/temp/zoinputfiles/doc_c200.png', height=500, width=800, res=130)
plot(DOC_mtms$time[seq(1,nrow(DOC_mtms),by=10)], DOC_mtms$c200[seq(1,nrow(DOC_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='DOC Concentration', ylim=c(0,6.75))
points(DOC_mtusgs$time[seq(1,nrow(DOC_mtusgs),by=50)], DOC_mtusgs$c200[seq(1,nrow(DOC_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topleft', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'),  bty='n', bg='white')
dev.off()

# As
png('c:/temp/zoinputfiles/as_c200.png', height=500, width=800, res=130)
plot(As_mtms$time[seq(1,nrow(As_mtms),by=10)], As_mtms$c200[seq(1,nrow(As_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='As Concentration', ylim=c(0,0.0125))
points(As_mtusgs$time[seq(1,nrow(As_mtusgs),by=50)], As_mtusgs$c200[seq(1,nrow(As_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'),  bty='n', bg='white')
dev.off()

# Mn
png('c:/temp/zoinputfiles/mn_c200.png', height=500, width=800, res=130)
plot(Mn_mtms$time[seq(1,nrow(Mn_mtms),by=10)], Mn_mtms$c200[seq(1,nrow(Mn_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='Mn Concentration', ylim=c(0,0.35))
points(Mn_mtusgs$time[seq(1,nrow(Mn_mtusgs),by=50)], Mn_mtusgs$c200[seq(1,nrow(Mn_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'),  bty='n', bg='white')
dev.off()

# NH4
png('c:/temp/zoinputfiles/nh4_c200.png', height=500, width=800, res=130)
plot(NH4_mtms$time[seq(1,nrow(NH4_mtms),by=10)], NH4_mtms$c200[seq(1,nrow(NH4_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='NH4 Concentration', ylim=c(0,9))
points(NH4_mtusgs$time[seq(1,nrow(NH4_mtusgs),by=50)], NH4_mtusgs$c200[seq(1,nrow(NH4_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('bottomright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# NO3
png('c:/temp/zoinputfiles/no3_c200.png', height=500, width=800, res=130)
plot(NO3_mtms$time[seq(1,nrow(NO3_mtms),by=10)], NO3_mtms$c200[seq(1,nrow(NO3_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='NO3 Concentration', ylim=c(-0.05,0.05))
points(NO3_mtusgs$time[seq(1,nrow(NO3_mtusgs),by=50)], NO3_mtusgs$c200[seq(1,nrow(NO3_mtusgs),by=50)], typ='p', pch=16, col='red')
abline(h=0)
legend('topright', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()

# P
png('c:/temp/zoinputfiles/p_c200.png', height=500, width=800, res=130)
plot(P_mtms$time[seq(1,nrow(P_mtms),by=10)], P_mtms$c200[seq(1,nrow(P_mtms),by=10)], typ='l', col='black', lwd=2, las=1, xaxs='i', yaxs='i', xlab='time', ylab='P Concentration', ylim=c(0,2))
points(P_mtusgs$time[seq(1,nrow(P_mtusgs),by=50)], P_mtusgs$c200[seq(1,nrow(P_mtusgs),by=50)], typ='p', pch=16, col='red')
legend('bottomleft', c('MT3DMS','MT3D-USGS'), pch=c(NA,16), lwd=c(2,NA), col=c('black','red'), bty='n', bg='white')
dev.off()