Nothing
inst/doc/ResponseModels-knitr.R
In kitagawa: Spectral Response of Water Wells to Harmonic Strain and Pressure Signals
## ----include=FALSE------------------------------------------------------------
library(knitr)
opts_chunk$set(
concordance=TRUE
)
## ----eval=TRUE, echo=FALSE-----------------------------------------
library(knitr)
options(width=69)
knitr::opts_chunk$set(tidy = TRUE, size="small")
## ----eval=TRUE, echo=TRUE, label=PRELIM----------------------------
library(RColorBrewer)
Set1 <- brewer.pal(8, "Set1")
library(signal, warn.conflicts=FALSE)
library(kitagawa)
## ----eval=TRUE, echo=TRUE, label=PARAMS1---------------------------
S. <- 1e-5 # Storativity [nondimensional]
T. <- 1e-4 # Transmissivity [m**2 / s]
D. <- T./S. # Diffusivity [m**2 / s]
Ta <- 50 # Aquifer thickness [m] #100
Hw <- z <- 50 # Depth to water table [m] #10
# Using ANO1 stats from Kit Tbl 1
Rc. <- 0.075 # Radius of cased portion of well [m]
Lc. <- 570 # Length of cased portion of well [m]
Rs. <- 0.135 # Radius of screened portion of well [m]
Ls. <- 15 # Length of screened portion of well [m]
Vw. <- sensing_volume(Rc., Lc., Rs., Ls.) # volume of fluid [m**3]
#
# parameters assumed by well_response:
# rho=1000 # density of rock [kg/m**3]
# Kf=2.2e9 # Bulk modulus of fluid [Pascals]
# grav=9.81 # gravitational acceleration [m/s**2]
rhog <- 9.81*1000
# Kitagawa Fig 7: Ku B / Kw Aw = 3 => Aw==4.8 at 40GPa
Ku. <- 40e9 # Bulk modulus [Pascals]
B. <- 0.5 # Skemptons ratio [nondimensional]
## ----eval=TRUE, echo=TRUE, label=PARAMS2---------------------------
# Frequencies
Q <- 10**seq(-5,2,by=0.05) # [nondimensional]
lQ <- log10(Q)
omega <- omega_norm(Q, z, D., invert=TRUE) # [Hz]
Phase <- function(Z){
Phs. <- Arg(Z) # will wrap to -pi/pi
uPhs. <- signal::unwrap(Phs., tol=pi/30)
return(data.frame(Phs=Phs., uPhs=uPhs.))
}
# Responses converted to pressure if TRUE
asP <- FALSE
ZasP <- FALSE
## ----eval=TRUE, echo=TRUE, fig.height=4.7, fig.width=4, label=KITRESP----
wrsp <- well_response(omega, T.=T., S.=S., Vw.=Vw., Rs.=Rs., Ku.=Ku., B.=B., Avs=1, Aw=1, as.pressure=asP)
plot(wrsp) # uses plot.wrsp method
crsp <- wrsp[["Response"]][,2] # Complex response
kGain <- Mod(crsp)/Ku./B. # Amplitude (or Gain)
kP <- Phase(crsp) # Phase
## ----eval=TRUE, echo=FALSE, fig.height=7, fig.width=5.5, label=KITRESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
plot(lQ, (kGain), type="l", #ylim=c(0, 1.15),
xaxt="n", #yaxs="i",
lwd=2,
ylab=expression( Z / Epsilon * B * kappa[u] ),
xlab="",
main="")
log10_ticks()
mtext("Sealed Well Response (KITAGAWA): Harmonic Strain", font=2, line=1.5)
mtext("(a) Gain", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, kP$Phs*180/pi, type="l", lty=3, #ylim=c(-190, -130),
xaxt="n",
ylab=expression(Z ~ "rel." ~ Epsilon),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
abline(h=-180, col="grey")
lines(lQ, kP$uPhs*180/pi, type="l", lwd=2)
log10_ticks()
mtext("(b) Phase", adj=0)
## ----eval=TRUE, echo=TRUE, fig.height=4.7, fig.width=4, label=COOPERRESP----
wrsp <- open_well_response(omega, T.=T., S.=S., Ta=Ta, Hw=Hw,
model = "cooper", as.pressure=ZasP)
plot(wrsp)
crsp <- wrsp[["Response"]][,2]
cGain <- Mod(crsp)
cP <- Phase(crsp)
## ----eval=TRUE, echo=FALSE, fig.height=7, fig.width=5.5, label=COOPERRESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
plot(lQ, cGain, type="l", #ylim=c(0, 1.15),
xaxt="n", #yaxs="i",
lwd=2,
ylab=expression( Z / H ),
xlab="",
main="")
log10_ticks()
mtext("Open Well Response (COOPER): Harmonic Strain", font=2, line=1.5)
mtext("(a) Gain", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, cP$Phs*180/pi, type="l", lwd=2, #ylim=c(-190, -130),
xaxt="n",
ylab=expression(Z ~ "rel." ~ H),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
abline(h=-180, col="grey")
lines(lQ, cP$uPhs*180/pi, type="l", lty=3)
log10_ticks()
mtext("(b) Phase", adj=0)
## ----eval=TRUE, echo=TRUE, fig.height=4.7, fig.width=4, label=HSIEHRESP----
wrsp <- open_well_response(omega, T.=T., S.=S., Ta=Ta, Hw=Hw, model = "hsieh", as.pressure=ZasP)
plot(wrsp)
crsp <- wrsp[["Response"]][,2]
hGain <- Mod(crsp)
hP <- Phase(crsp)
## ----eval=TRUE, echo=FALSE, fig.height=7, fig.width=5.5, label=HSIEHRESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
plot(lQ, hGain, type="l", #ylim=c(0, 1.15),
xaxt="n", #yaxs="i",
lwd=2,
ylab=expression( Z / H ),
xlab="",
main="")
log10_ticks()
mtext("Open Well Response (HSIEH): Harmonic Strain", font=2, line=1.5)
mtext("(a) Gain", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, hP$Phs*180/pi, type="l", lty=3, #ylim=c(-190, -130),
xaxt="n",
ylab=expression(Z ~ "rel." ~ H),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
abline(h=-180, col="grey")
lines(lQ, hP$uPhs*180/pi, type="l", lwd=2)
log10_ticks()
mtext("(b) Phase", adj=0)
## ----eval=TRUE, echo=TRUE, fig.height=4.7, fig.width=4, label=LIURESP----
wrsp <- open_well_response(omega, T.=T., S.=S., Ta=Ta, Hw=Hw, model = "liu", as.pressure=ZasP)
plot(wrsp)
crsp <- wrsp[["Response"]][,2]
lGain <- Mod(crsp)
lP <- Phase(crsp)
## ----eval=TRUE, echo=FALSE, fig.height=7, fig.width=5.5, label=LIURESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
plot(lQ, lGain, type="l", #ylim=c(0, 1.15),
xaxt="n", #yaxs="i",
lwd=2,
ylab=expression( Z / H ),
xlab="",
main="")
log10_ticks()
mtext("Open Well Response (LIU): Harmonic Strain", font=2, line=1.5)
mtext("(a) Gain", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, lP$Phs*180/pi, type="l", lwd = 2, #ylim=c(-190, -130),
xaxt="n",
ylab=expression(Z ~ "rel." ~ H),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
abline(h=-180, col="grey")
lines(lQ, lP$uPhs*180/pi, type="l", lty=3)
log10_ticks()
mtext("(b) Phase", adj=0)
## ----eval=TRUE, echo=TRUE, fig.height=4.7, fig.width=4, label=WANGRESP----
wrsp <- open_well_response(omega, T.=T., S.=S., leak = 1e-8,
model = "wang", as.pressure=asP)
plot(wrsp)
crsp <- wrsp[["Response"]][,2]
rGain <- Mod(crsp)
rP <- Phase(crsp)
## ----eval=TRUE, echo=TRUE, label=WANGRESPDAT-----------------------
Transmiss <- c(1e0, 1e-2, 1e-4, 1e-6, 1e-8)
Storativ <- c(1e-2, 1e-4, 1e-6, 1e-8)
omeg <- 1.9322736 / 86400 # M2 in Hz
leak <- 10^seq(-11, -3, 0.2)
## ----eval=TRUE, echo=FALSE, fig.height=4.7, fig.width=4, label=WANGRESPFIG----
pt.cex = 0.35
layout(matrix(c(1,2), ncol=1), heights=c(0.5,0.5))
par(mar=c(2,4,1,1),
oma=c(2,0.1,1,0.1),
tcl=-0.3,
mgp=c(2.5, 0.5, 0), las=1)
plot(c(), c(), log = 'x',
xlim = range(leak), type='l', ylim = c(-80, 100),
ylab = 'Phase shift (deg)', xaxt='n', frame=FALSE)
lats <- seq(-11,-3,by=2)
axis(1, at=10^lats, labels=parse(text=sprintf("10^%s",lats)))
axis(3, at=10^lats, labels=FALSE)
axis(4, labels=FALSE)
box()
for (i in seq_along(Transmiss)) {
T. <- Transmiss[i]
for (j in seq_along(Storativ)) {
S. <- Storativ[j]
wang <- open_well_response(omeg,
T.,
S.,
Rs. = 0.1,
model = 'wang',
leak = leak,
freq.units = 'Hz',
as.pressure = FALSE)
wang <- wang[["Response"]][,2]
col <- ifelse(i <=3, 'firebrick4', i)
if(i == 4) col <- 'dodgerblue4'
if(i == 5) col <- 'forestgreen'
points(x = leak, Arg(wang) * 180/pi, type = 'o', col = col,
pch = j, cex = pt.cex)
}
}
legend('bottomright',
col = c('firebrick4', 'dodgerblue4', 'forestgreen'),
lty = 1,
legend = c('1.0 >= T >= 1e-4',
'T = 1e-6',
'T = 1e-8'),
cex = 0.5)
plot(c(), c(), log = 'xy',
xlim = range(leak), type='l', ylim = c(1e-10, 1),
ylab = 'Amplitude ratio',
xlab = 'Specific Leakage (1/s)', axes=FALSE)
axis(1, at=10^lats, labels=parse(text=sprintf("10^%s",lats)))
lats.y <- seq(-10,-1,by=3)
axis(3, at=10^lats, labels=FALSE)
axis(2, at=10^lats.y, labels=parse(text=sprintf("10^%s",lats.y)))
axis(4, at=10^lats.y, labels=FALSE)
box()
for (i in seq_along(Transmiss)) {
T. <- Transmiss[i]
for (j in seq_along(Storativ)) {
S. <- Storativ[j]
wang <- open_well_response(omeg,
T.,
S.,
Rs. = 0.1,
model = 'wang',
leak = leak,
freq.units = 'Hz',
as.pressure = FALSE)
wang <- wang[["Response"]][,2]
col <- ifelse(i <=3, 'firebrick4', i)
if(i == 4) col <- 'dodgerblue4'
if(i == 5) col <- 'forestgreen'
points(x = leak, Mod(wang), type = 'o', col = col,
pch = j, cex = pt.cex)
}
}
legend('bottomleft',
pch = 1:4,
legend = c('S = 1e-2',
'S = 1e-4',
'S = 1e-6',
'S = 1e-8'),
pt.cex = pt.cex,
cex = 0.5)
mtext("Leakage Coefficient (1/s)", side=1, line=2)
## ----eval=TRUE, echo=TRUE, fig.height=4.7, fig.width=4, label=ROJRESP----
wrsp <- open_well_response(omega, T.=T., S.=S., z=z, model = "rojstaczer", as.pressure=asP)
plot(wrsp)
crsp <- wrsp[["Response"]][,2]
rGain <- Mod(crsp)
rP <- Phase(crsp)
## ----eval=TRUE, echo=FALSE, fig.height=6, fig.width=5.5, label=ROJRESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
plot(lQ, (rGain), type="l",
#ylim=c(0, 1.15),
xaxt="n",
lwd=2,
ylab=expression( Z / Epsilon ),
xlab="",
main="")
log10_ticks()
mtext("Open Well Response (ROJSTACZER): Harmonic Strain", font=2, line=1.5)
mtext("(a) Gain", adj=0)
# relative to static-confined areal strain response", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, rP$Phs*180/pi, type="l", lty=3, ylim=c(130, 180),
xaxt="n",
ylab=expression(Z ~ "rel." ~ Epsilon),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
abline(h=-180, col="grey")
lines(lQ, rP$uPhs*180/pi, type="l", lwd=2)
log10_ticks()
mtext("(b) Phase", adj=0)
## ----eval=TRUE, echo=FALSE, fig.height=6, fig.width=5.5, label=ALLRESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
# kitagawa and rojstaczer
alim <- range(pretty(log10(c(kGain,rGain))))
plot(lQ, log10(kGain), col=Set1[1],
type="l", ylim=alim,
#ylim=c(-2.5, 0.2),
yaxt="n", xaxt="n",
lwd=2,
ylab=expression(log[10] ~ Z / Epsilon),
xlab="",
main="")
lines(lQ, log10(rGain), lwd=2, col=Set1[2])
legend("bottomright",
c("Kitagawa et al (2011) -- sealed","Rojstaczer et al (1988) -- open"),
lwd=3, col=Set1[1:2], bty="n")
log10_ticks(2, major.ticks=-5:5)
log10_ticks()
mtext("Harmonic Strain Well Responses", font=2, line=1.5)
mtext("(a) Gain", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, kP$Phs*180/pi -180, lty=3, lwd=1.5, col=Set1[1],
type="l", ylim=90*c(0,-1),
yaxt="n",
xaxt="n",
ylab=expression(Z ~ "rel. -180" ~ Epsilon),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
lbls <- ats <- seq(-90,90,by=15)
lbls[seq_along(lbls)%%2==0] <- ""
axis(2, at=ats, labels=lbls)
lines(lQ, rP$Phs*180/pi-180, lty=3, lwd=1.5, col=Set1[2])
# unwrapped phase
lines(lQ, kP$uPhs*180/pi-180, lwd=2, col=Set1[1])
lines(lQ, rP$uPhs*180/pi-180, lwd=2, col=Set1[2])
log10_ticks()
mtext("(b) Anti-Phase", adj=0)
## ----eval=TRUE, echo=FALSE, fig.height=6, fig.width=5.5, label=ALLORESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
# cooper, liu, hsieh
alim <- range(pretty(log10(c(hGain,cGain,lGain))))
plot(lQ, log10(hGain), col=Set1[3],
type="l", ylim=alim,
xaxt="n", yaxt="n",
lwd=2,
ylab=expression(log[10] ~ Z/H),
xlab="",
main="")
lines(lQ, log10(cGain), lwd=2, col=Set1[1])
lines(lQ, log10(lGain), lwd=2, col=Set1[2])
legend("bottomleft",
c("Cooper et al (1965)","Liu et al (1989)","Hsieh et al (1987)"),
lwd=3, col=Set1[1:3], bty="n")
log10_ticks(2, major.ticks=-9:2)
log10_ticks()
mtext("Harmonic Pressure-head Well Responses (Open)", font=2, line=1.5)
mtext("(a) Gain", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, cP$uPhs*180/pi, lty=3, lwd=1.5, col=Set1[1],
type="l", ylim=185*c(0,-1),
yaxt="n",
xaxt="n",
ylab=expression(Z ~ "rel." ~ H),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
lbls <- ats <- seq(-180,180,by=30)
lbls[seq_along(lbls)%%2==0] <- ""
axis(2, at=ats, labels=lbls)
lines(lQ, lP$uPhs*180/pi, lty=3, lwd=1.5, col=Set1[2])
lines(lQ, hP$uPhs*180/pi, lty=3, lwd=1.5, col=Set1[3])
# unwrapped phase
lines(lQ, cP$Phs*180/pi, lwd=2, col=Set1[1])
lines(lQ, lP$Phs*180/pi, lwd=2, col=Set1[2])
lines(lQ, hP$Phs*180/pi, lwd=2, col=Set1[3])
log10_ticks()
mtext("(b) Phase", adj=0)
Try the kitagawa package in your browser
Any scripts or data that you put into this service are public.
kitagawa documentation built on July 2, 2020, 1:47 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----include=FALSE------------------------------------------------------------
library(knitr)
opts_chunk$set(
concordance=TRUE
)
## ----eval=TRUE, echo=FALSE-----------------------------------------
library(knitr)
options(width=69)
knitr::opts_chunk$set(tidy = TRUE, size="small")
## ----eval=TRUE, echo=TRUE, label=PRELIM----------------------------
library(RColorBrewer)
Set1 <- brewer.pal(8, "Set1")
library(signal, warn.conflicts=FALSE)
library(kitagawa)
## ----eval=TRUE, echo=TRUE, label=PARAMS1---------------------------
S. <- 1e-5 # Storativity [nondimensional]
T. <- 1e-4 # Transmissivity [m**2 / s]
D. <- T./S. # Diffusivity [m**2 / s]
Ta <- 50 # Aquifer thickness [m] #100
Hw <- z <- 50 # Depth to water table [m] #10
# Using ANO1 stats from Kit Tbl 1
Rc. <- 0.075 # Radius of cased portion of well [m]
Lc. <- 570 # Length of cased portion of well [m]
Rs. <- 0.135 # Radius of screened portion of well [m]
Ls. <- 15 # Length of screened portion of well [m]
Vw. <- sensing_volume(Rc., Lc., Rs., Ls.) # volume of fluid [m**3]
#
# parameters assumed by well_response:
# rho=1000 # density of rock [kg/m**3]
# Kf=2.2e9 # Bulk modulus of fluid [Pascals]
# grav=9.81 # gravitational acceleration [m/s**2]
rhog <- 9.81*1000
# Kitagawa Fig 7: Ku B / Kw Aw = 3 => Aw==4.8 at 40GPa
Ku. <- 40e9 # Bulk modulus [Pascals]
B. <- 0.5 # Skemptons ratio [nondimensional]
## ----eval=TRUE, echo=TRUE, label=PARAMS2---------------------------
# Frequencies
Q <- 10**seq(-5,2,by=0.05) # [nondimensional]
lQ <- log10(Q)
omega <- omega_norm(Q, z, D., invert=TRUE) # [Hz]
Phase <- function(Z){
Phs. <- Arg(Z) # will wrap to -pi/pi
uPhs. <- signal::unwrap(Phs., tol=pi/30)
return(data.frame(Phs=Phs., uPhs=uPhs.))
}
# Responses converted to pressure if TRUE
asP <- FALSE
ZasP <- FALSE
## ----eval=TRUE, echo=TRUE, fig.height=4.7, fig.width=4, label=KITRESP----
wrsp <- well_response(omega, T.=T., S.=S., Vw.=Vw., Rs.=Rs., Ku.=Ku., B.=B., Avs=1, Aw=1, as.pressure=asP)
plot(wrsp) # uses plot.wrsp method
crsp <- wrsp[["Response"]][,2] # Complex response
kGain <- Mod(crsp)/Ku./B. # Amplitude (or Gain)
kP <- Phase(crsp) # Phase
## ----eval=TRUE, echo=FALSE, fig.height=7, fig.width=5.5, label=KITRESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
plot(lQ, (kGain), type="l", #ylim=c(0, 1.15),
xaxt="n", #yaxs="i",
lwd=2,
ylab=expression( Z / Epsilon * B * kappa[u] ),
xlab="",
main="")
log10_ticks()
mtext("Sealed Well Response (KITAGAWA): Harmonic Strain", font=2, line=1.5)
mtext("(a) Gain", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, kP$Phs*180/pi, type="l", lty=3, #ylim=c(-190, -130),
xaxt="n",
ylab=expression(Z ~ "rel." ~ Epsilon),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
abline(h=-180, col="grey")
lines(lQ, kP$uPhs*180/pi, type="l", lwd=2)
log10_ticks()
mtext("(b) Phase", adj=0)
## ----eval=TRUE, echo=TRUE, fig.height=4.7, fig.width=4, label=COOPERRESP----
wrsp <- open_well_response(omega, T.=T., S.=S., Ta=Ta, Hw=Hw,
model = "cooper", as.pressure=ZasP)
plot(wrsp)
crsp <- wrsp[["Response"]][,2]
cGain <- Mod(crsp)
cP <- Phase(crsp)
## ----eval=TRUE, echo=FALSE, fig.height=7, fig.width=5.5, label=COOPERRESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
plot(lQ, cGain, type="l", #ylim=c(0, 1.15),
xaxt="n", #yaxs="i",
lwd=2,
ylab=expression( Z / H ),
xlab="",
main="")
log10_ticks()
mtext("Open Well Response (COOPER): Harmonic Strain", font=2, line=1.5)
mtext("(a) Gain", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, cP$Phs*180/pi, type="l", lwd=2, #ylim=c(-190, -130),
xaxt="n",
ylab=expression(Z ~ "rel." ~ H),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
abline(h=-180, col="grey")
lines(lQ, cP$uPhs*180/pi, type="l", lty=3)
log10_ticks()
mtext("(b) Phase", adj=0)
## ----eval=TRUE, echo=TRUE, fig.height=4.7, fig.width=4, label=HSIEHRESP----
wrsp <- open_well_response(omega, T.=T., S.=S., Ta=Ta, Hw=Hw, model = "hsieh", as.pressure=ZasP)
plot(wrsp)
crsp <- wrsp[["Response"]][,2]
hGain <- Mod(crsp)
hP <- Phase(crsp)
## ----eval=TRUE, echo=FALSE, fig.height=7, fig.width=5.5, label=HSIEHRESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
plot(lQ, hGain, type="l", #ylim=c(0, 1.15),
xaxt="n", #yaxs="i",
lwd=2,
ylab=expression( Z / H ),
xlab="",
main="")
log10_ticks()
mtext("Open Well Response (HSIEH): Harmonic Strain", font=2, line=1.5)
mtext("(a) Gain", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, hP$Phs*180/pi, type="l", lty=3, #ylim=c(-190, -130),
xaxt="n",
ylab=expression(Z ~ "rel." ~ H),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
abline(h=-180, col="grey")
lines(lQ, hP$uPhs*180/pi, type="l", lwd=2)
log10_ticks()
mtext("(b) Phase", adj=0)
## ----eval=TRUE, echo=TRUE, fig.height=4.7, fig.width=4, label=LIURESP----
wrsp <- open_well_response(omega, T.=T., S.=S., Ta=Ta, Hw=Hw, model = "liu", as.pressure=ZasP)
plot(wrsp)
crsp <- wrsp[["Response"]][,2]
lGain <- Mod(crsp)
lP <- Phase(crsp)
## ----eval=TRUE, echo=FALSE, fig.height=7, fig.width=5.5, label=LIURESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
plot(lQ, lGain, type="l", #ylim=c(0, 1.15),
xaxt="n", #yaxs="i",
lwd=2,
ylab=expression( Z / H ),
xlab="",
main="")
log10_ticks()
mtext("Open Well Response (LIU): Harmonic Strain", font=2, line=1.5)
mtext("(a) Gain", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, lP$Phs*180/pi, type="l", lwd = 2, #ylim=c(-190, -130),
xaxt="n",
ylab=expression(Z ~ "rel." ~ H),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
abline(h=-180, col="grey")
lines(lQ, lP$uPhs*180/pi, type="l", lty=3)
log10_ticks()
mtext("(b) Phase", adj=0)
## ----eval=TRUE, echo=TRUE, fig.height=4.7, fig.width=4, label=WANGRESP----
wrsp <- open_well_response(omega, T.=T., S.=S., leak = 1e-8,
model = "wang", as.pressure=asP)
plot(wrsp)
crsp <- wrsp[["Response"]][,2]
rGain <- Mod(crsp)
rP <- Phase(crsp)
## ----eval=TRUE, echo=TRUE, label=WANGRESPDAT-----------------------
Transmiss <- c(1e0, 1e-2, 1e-4, 1e-6, 1e-8)
Storativ <- c(1e-2, 1e-4, 1e-6, 1e-8)
omeg <- 1.9322736 / 86400 # M2 in Hz
leak <- 10^seq(-11, -3, 0.2)
## ----eval=TRUE, echo=FALSE, fig.height=4.7, fig.width=4, label=WANGRESPFIG----
pt.cex = 0.35
layout(matrix(c(1,2), ncol=1), heights=c(0.5,0.5))
par(mar=c(2,4,1,1),
oma=c(2,0.1,1,0.1),
tcl=-0.3,
mgp=c(2.5, 0.5, 0), las=1)
plot(c(), c(), log = 'x',
xlim = range(leak), type='l', ylim = c(-80, 100),
ylab = 'Phase shift (deg)', xaxt='n', frame=FALSE)
lats <- seq(-11,-3,by=2)
axis(1, at=10^lats, labels=parse(text=sprintf("10^%s",lats)))
axis(3, at=10^lats, labels=FALSE)
axis(4, labels=FALSE)
box()
for (i in seq_along(Transmiss)) {
T. <- Transmiss[i]
for (j in seq_along(Storativ)) {
S. <- Storativ[j]
wang <- open_well_response(omeg,
T.,
S.,
Rs. = 0.1,
model = 'wang',
leak = leak,
freq.units = 'Hz',
as.pressure = FALSE)
wang <- wang[["Response"]][,2]
col <- ifelse(i <=3, 'firebrick4', i)
if(i == 4) col <- 'dodgerblue4'
if(i == 5) col <- 'forestgreen'
points(x = leak, Arg(wang) * 180/pi, type = 'o', col = col,
pch = j, cex = pt.cex)
}
}
legend('bottomright',
col = c('firebrick4', 'dodgerblue4', 'forestgreen'),
lty = 1,
legend = c('1.0 >= T >= 1e-4',
'T = 1e-6',
'T = 1e-8'),
cex = 0.5)
plot(c(), c(), log = 'xy',
xlim = range(leak), type='l', ylim = c(1e-10, 1),
ylab = 'Amplitude ratio',
xlab = 'Specific Leakage (1/s)', axes=FALSE)
axis(1, at=10^lats, labels=parse(text=sprintf("10^%s",lats)))
lats.y <- seq(-10,-1,by=3)
axis(3, at=10^lats, labels=FALSE)
axis(2, at=10^lats.y, labels=parse(text=sprintf("10^%s",lats.y)))
axis(4, at=10^lats.y, labels=FALSE)
box()
for (i in seq_along(Transmiss)) {
T. <- Transmiss[i]
for (j in seq_along(Storativ)) {
S. <- Storativ[j]
wang <- open_well_response(omeg,
T.,
S.,
Rs. = 0.1,
model = 'wang',
leak = leak,
freq.units = 'Hz',
as.pressure = FALSE)
wang <- wang[["Response"]][,2]
col <- ifelse(i <=3, 'firebrick4', i)
if(i == 4) col <- 'dodgerblue4'
if(i == 5) col <- 'forestgreen'
points(x = leak, Mod(wang), type = 'o', col = col,
pch = j, cex = pt.cex)
}
}
legend('bottomleft',
pch = 1:4,
legend = c('S = 1e-2',
'S = 1e-4',
'S = 1e-6',
'S = 1e-8'),
pt.cex = pt.cex,
cex = 0.5)
mtext("Leakage Coefficient (1/s)", side=1, line=2)
## ----eval=TRUE, echo=TRUE, fig.height=4.7, fig.width=4, label=ROJRESP----
wrsp <- open_well_response(omega, T.=T., S.=S., z=z, model = "rojstaczer", as.pressure=asP)
plot(wrsp)
crsp <- wrsp[["Response"]][,2]
rGain <- Mod(crsp)
rP <- Phase(crsp)
## ----eval=TRUE, echo=FALSE, fig.height=6, fig.width=5.5, label=ROJRESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
plot(lQ, (rGain), type="l",
#ylim=c(0, 1.15),
xaxt="n",
lwd=2,
ylab=expression( Z / Epsilon ),
xlab="",
main="")
log10_ticks()
mtext("Open Well Response (ROJSTACZER): Harmonic Strain", font=2, line=1.5)
mtext("(a) Gain", adj=0)
# relative to static-confined areal strain response", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, rP$Phs*180/pi, type="l", lty=3, ylim=c(130, 180),
xaxt="n",
ylab=expression(Z ~ "rel." ~ Epsilon),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
abline(h=-180, col="grey")
lines(lQ, rP$uPhs*180/pi, type="l", lwd=2)
log10_ticks()
mtext("(b) Phase", adj=0)
## ----eval=TRUE, echo=FALSE, fig.height=6, fig.width=5.5, label=ALLRESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
# kitagawa and rojstaczer
alim <- range(pretty(log10(c(kGain,rGain))))
plot(lQ, log10(kGain), col=Set1[1],
type="l", ylim=alim,
#ylim=c(-2.5, 0.2),
yaxt="n", xaxt="n",
lwd=2,
ylab=expression(log[10] ~ Z / Epsilon),
xlab="",
main="")
lines(lQ, log10(rGain), lwd=2, col=Set1[2])
legend("bottomright",
c("Kitagawa et al (2011) -- sealed","Rojstaczer et al (1988) -- open"),
lwd=3, col=Set1[1:2], bty="n")
log10_ticks(2, major.ticks=-5:5)
log10_ticks()
mtext("Harmonic Strain Well Responses", font=2, line=1.5)
mtext("(a) Gain", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, kP$Phs*180/pi -180, lty=3, lwd=1.5, col=Set1[1],
type="l", ylim=90*c(0,-1),
yaxt="n",
xaxt="n",
ylab=expression(Z ~ "rel. -180" ~ Epsilon),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
lbls <- ats <- seq(-90,90,by=15)
lbls[seq_along(lbls)%%2==0] <- ""
axis(2, at=ats, labels=lbls)
lines(lQ, rP$Phs*180/pi-180, lty=3, lwd=1.5, col=Set1[2])
# unwrapped phase
lines(lQ, kP$uPhs*180/pi-180, lwd=2, col=Set1[1])
lines(lQ, rP$uPhs*180/pi-180, lwd=2, col=Set1[2])
log10_ticks()
mtext("(b) Anti-Phase", adj=0)
## ----eval=TRUE, echo=FALSE, fig.height=6, fig.width=5.5, label=ALLORESPFIG----
par(mfrow=c(2,1),
mar=c(2.5,4,3.1,2),
oma=rep(0.1,4),
omi=rep(0.1,4),
las=1)
# cooper, liu, hsieh
alim <- range(pretty(log10(c(hGain,cGain,lGain))))
plot(lQ, log10(hGain), col=Set1[3],
type="l", ylim=alim,
xaxt="n", yaxt="n",
lwd=2,
ylab=expression(log[10] ~ Z/H),
xlab="",
main="")
lines(lQ, log10(cGain), lwd=2, col=Set1[1])
lines(lQ, log10(lGain), lwd=2, col=Set1[2])
legend("bottomleft",
c("Cooper et al (1965)","Liu et al (1989)","Hsieh et al (1987)"),
lwd=3, col=Set1[1:3], bty="n")
log10_ticks(2, major.ticks=-9:2)
log10_ticks()
mtext("Harmonic Pressure-head Well Responses (Open)", font=2, line=1.5)
mtext("(a) Gain", adj=0)
par(mar=c(4,4,1.1,2))
plot(lQ, cP$uPhs*180/pi, lty=3, lwd=1.5, col=Set1[1],
type="l", ylim=185*c(0,-1),
yaxt="n",
xaxt="n",
ylab=expression(Z ~ "rel." ~ H),
xlab=expression("Dimensionless frequency," ~ Q ==z^2 * omega / 2 * D ))
lbls <- ats <- seq(-180,180,by=30)
lbls[seq_along(lbls)%%2==0] <- ""
axis(2, at=ats, labels=lbls)
lines(lQ, lP$uPhs*180/pi, lty=3, lwd=1.5, col=Set1[2])
lines(lQ, hP$uPhs*180/pi, lty=3, lwd=1.5, col=Set1[3])
# unwrapped phase
lines(lQ, cP$Phs*180/pi, lwd=2, col=Set1[1])
lines(lQ, lP$Phs*180/pi, lwd=2, col=Set1[2])
lines(lQ, hP$Phs*180/pi, lwd=2, col=Set1[3])
log10_ticks()
mtext("(b) Phase", adj=0)
Try the kitagawa package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.