\documentclass[12pt]{article}
% !Rnw weave = knitr
%
%\makeatletter
% \VignetteIndexEntry{Comparison of response functions}
% \VignetteEngine{knitr::knitr}
%\makeatother
%

%\usepackage{url}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{fancyvrb}
\usepackage[pdfborder={0 0 0}]{hyperref}
\usepackage{bookmark}
\usepackage{url}
\usepackage{upquote}
\usepackage{graphicx}
\usepackage{grffile}
\usepackage{float}
\usepackage{natbib}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{hyperref}
\usepackage{geometry}
\geometry{verbose,tmargin=3cm,bmargin=5cm,lmargin=2.5cm,rmargin=2.5cm}
\usepackage[font=sf, labelfont={sf,bf}, margin=2cm]{caption}
\usepackage{color}
\raggedright
%
\author{Andrew J. Barbour$^1$ $\cdot$ Jonathan Kennel$^2$}
\title{Comparison of response functions in \kit{}}
\date{%
    \footnotesize{$1$ - U.S. Geological Survey $\cdot$ $2$ - University of Guelph}\\[2ex]%
    \today
}
%
\begin{document}

<<include=FALSE>>=
library(knitr)
opts_chunk$set(
concordance=TRUE
)
@

%\SweaveOpts{concordance=TRUE}
%
\newcommand{\SC}[1]{\textsc{#1}}
\newcommand{\Rcmd}[1]{\texttt{#1}}
\newcommand{\kit}[0]{\href{https://github.com/abarbour/kitagawa/}{\color{blue}\Rcmd{kitagawa}}}
\newcommand{\bidxa}[1]{\index{#1}{\textbf{#1}}} 
\newcommand{\bidxb}[2]{\index{#2}{\textbf{#1}}} 
\newcommand{\idxa}[1]{\index{#1}{#1}} 
\newcommand{\idxb}[2]{\index{#2}{#1}} 
%
\maketitle
%
\begin{abstract}
In this vignette we demonstrate the response functions found in 
the package
\kit{}, which are appropriate for modeling the effect of
harmonic volumetric strain or pressure-head fluctuations in
sealed and open water wells.
For sealed-wells there is only one response function, from
\citet{kitagawa2011}, and this gives the complex frequency response of
virtual water height $Z$ or pressure $P$ as a function of
areal %volumetric?
strain $\epsilon$.
For open wells there is
a suite of open-well response functions, from 
\citet{cooper1965, hsieh1987, rojstaczer1988, liu1989, wang2018}; and
these give the complex frequency response of
water height as a function of aquifer head $H$ or pressure. \citet{wang2018}
allows for leakage from the aquifer.
\end{abstract}
%
\tableofcontents

%\clearpage
\section{Introduction}

 The underlying physical model of these response functions
 is based upon the assumption that fluid flows radially
 through an homogeneous, isotropic, confined aquifer.
%
 The underlying principle is as follows.  When a harmonic wave induces
 strain in a confined aquifer (one having aquitards above and below it), 
 fluid flows radially into, and out of a well penetrating the aquifer.
 The flow-induced drawdown, $s$, is governed by the following 
 partial differential equation, expressed in radial coordinates($r$):
\begin{equation}
 \frac{\partial^2 s}{\partial r^2} + \frac{1}{r} 
 \frac{\partial s}{ \partial r} - \frac{S}{T}\frac{\partial s}{\partial t} = 0
\end{equation}
 where $S$ and $T$ are the aquifer storativity and transmissivity respectively.
 
 The solution to this PDE, with periodic discharge boundary conditions,
 gives the amplitude and phase response we wish to calculate.
 The solution  for an open well was first presented by
 \citet{cooper1965}, and subsequently modified by \citet{rojstaczer1988, liu1989}.
 \citet{kitagawa2011} adapted the solution
 of \citet{hsieh1987} for the case of a sealed well. 
 \citet{wang2018} provides the leaky aquifer response for an open well.
 
 These models are applicable to any quasi-static process involving harmonic, 
 volumetric strain of an aquifer 
 (e.g., passing Rayleigh waves, or changes in the Earth's tidal potential). 
 In practice, however, the presence of permeable fractures can violate the
 assumption of isotropic permeability, which may substantially
 alter the response by introducing shear-strain coupling. Such
 complications are beyond the scope of these models.

\section{Preliminaries}

<<eval=TRUE, echo=FALSE>>=
library(knitr)
options(width=69)
knitr::opts_chunk$set(tidy = TRUE, size="small")
@
%opts_knit$set(verbose = TRUE)

Load the necessary packages:

<<eval=TRUE, echo=TRUE, label=PRELIM>>=
library(RColorBrewer)
Set1 <- brewer.pal(8, "Set1")
library(signal, warn.conflicts=FALSE)
library(kitagawa)
@

Setup some constants:

<<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]
@

And create the dimensionless frequencies, 
defined by $z^2 \omega / 2 D$,
where $D$ is the hydraulic diffusivity:

<<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
@

And onto the response functions...

\clearpage
\section{Sealed well response}

\subsection{Strain: Kitagawa et al. (2011)}
%
<<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
@

\begin{figure}[htb!]
\begin{center}
<<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)
@
\caption{The response of a sealed well to harmonic areal strain using
the Kitagawa model. The amplitude is normalized by Skempton's coefficient $B$
and the undrained bulk modulus $\kappa_u$.
Frequency is dimensionless, based on the well-depth $z$ and the diffusivity $D$.
}
\label{fig:wrsp}
\end{center}
\end{figure}

\clearpage
\section{Open well response}

\subsection{Pressure head: Cooper et al. (1965)}

<<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)
@

\begin{figure}[htb!]
\begin{center}
<<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)
@
\caption{The response of an open well to harmonic areal strain using
the Cooper model. 
Frequency is dimensionless, based on the well-depth $z$ and the diffusivity $D$.
}
\label{fig:owrsp-coop}
\end{center}
\end{figure}

\clearpage
\subsection{Pressure head: Hsieh et al. (1987)}

<<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)
@

\begin{figure}[htb!]
\begin{center}
<<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)
@
\caption{The response of an open well to harmonic areal strain using
the Hsieh model. 
Frequency is dimensionless, based on the well-depth $z$ and the diffusivity $D$.
}
\label{fig:owrsp-hsi}
\end{center}
\end{figure}

\clearpage
\subsection{Pressure head: Liu et al. (1989)}

<<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)
@

\begin{figure}[htb!]
\begin{center}
<<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)
@
\caption{The response of an open well to harmonic areal strain using
the Liu model. 
Frequency is dimensionless, based on the well-depth $z$ and the diffusivity $D$.
}
\label{fig:owrsp-liu}
\end{center}
\end{figure}



\clearpage
\subsection{Pressure head (with leakage): Wang et al. (2018)}

<<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)
@

\clearpage
\subsubsection{Figure 2 from Wang et al. (2018)}

<<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)
@

\begin{figure}[htb!]
\begin{center}
<<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)

@
\caption{Amplitude and phase shift as a function of the specific leakage (K'/b') using the Wang 2018 model for the M2 tide.}
\label{fig:owrsp-wang}
\end{center}
\end{figure}



\clearpage
\subsection{Strain: Rojstaczer (1988)}

<<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)
@

\begin{figure}[htb!]
\begin{center}
<<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)
@
\caption{The response of an open well to harmonic areal strain using
the Rojstaczer model. In the phase curve, phase wrapping has been removed.
Modified from \citet[][Fig.~3]{rojstaczer1988}.
Frequency is dimensionless, based on the well-depth $z$ and the diffusivity
$D$.}
\label{fig:owrsp-roj}
\end{center}
\end{figure}


\clearpage
\section{Model Comparisons}

\subsection{Responses to strain}
\begin{figure}[htb!]
\begin{center}
<<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)
@
\caption{A comparison of well responses to harmonic strain. 
The phase of the water level is relative to $-180^\circ$ the phase of strain.}
\label{fig:ewrsp-all}
\end{center}
\end{figure}

\clearpage
\subsection{Responses to pressure head (all open)}
\begin{figure}[htb!]
\begin{center}
<<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)
@
\caption{A comparison of well responses to harmonic pressure-head, 
from \citet{cooper1965, hsieh1987, liu1989} (all for unsealed).}
\label{fig:owrsp-all}
\end{center}
\end{figure}

\bibliographystyle{apalike}
\bibliography{REFS}

\end{document}