diff --git a/DESCRIPTION b/DESCRIPTION
index 9bb0d2e..1b8be54 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -33,5 +33,5 @@ Suggests:
     testthat,
     covr
 LinkingTo: Rcpp
-RoxygenNote: 7.1.0
+RoxygenNote: 7.1.1
 Encoding: UTF-8
diff --git a/NAMESPACE b/NAMESPACE
index f2f37d8..bc35ad1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -33,6 +33,7 @@ export(ierfc_re)
 export(lame_first)
 export(mogi.pressure)
 export(mogi.volume)
+export(r96_drainage_pressure)
 export(rigid_hydraulic_diffusivity)
 export(rudnicki86)
 export(segall85)
diff --git a/R/roeloffs96.R b/R/roeloffs96.R
new file mode 100644
index 0000000..14016e2
--- /dev/null
+++ b/R/roeloffs96.R
@@ -0,0 +1,55 @@
+#' Earthquake-related poroelastic solutions in Roeloffs 1996
+#' 
+#' @details \code{r96_drainage_pressure} refers to the solution
+#' given in equation 50 of Roeloffs (1996), the pore pressure
+#' response to a static strain. The temporal change in pore
+#' pressure represents water table drainage as a function of time.
+#' 
+#' @param Times numeric; times to calculate at \code{depth} [s]
+#' @param depth numeric; observation depth [m]
+#' @param depth_water numeric; depth to water table [m]
+#' @param diffusiv numeric; hydraulic diffusivity [m^2/s]
+#' @param strain0 numeric; static strain change [unscaled strain]
+#' @param B numeric; Skempton's coefficient B
+#' @param Ku numeric; undrained bulk modulus
+#' @param ... additional parameters
+#' 
+#' @export
+#' 
+#' @references Roeloffs, E. (1996),
+#' Poroelastic techniques in the study of earthquake-related hydrologic phenomena,
+#' \emph{Advances in Geophysics}, 37, 135-195, \url{https://doi.org/10.1016/S0065-2687(08)60270-8}
+#' 
+#' @examples 
+#' # Remake Fig 16
+#' d2s <- 86400
+#' lti <- 10^seq(-3,0, length.out=21)
+#' ti <- unique(sort(c(0, lti, seq(-10, 50, length.out=201)) * d2s))
+#' zsq4c <- c(1, 24, 24*30) * d2s / 24 # timescales in hours
+#' z <- 1
+#' vert_c <- z^2 / 4 / zsq4c # vertical diffusivity
+#' r1  <- r96_drainage_pressure(Times=ti, depth=z, diffusiv=vert_c[1]) # 1 hour
+#' r2  <- r96_drainage_pressure(Times=ti, depth=z, diffusiv=vert_c[2]) # 1 day
+#' r3  <- r96_drainage_pressure(Times=ti, depth=z, diffusiv=vert_c[3]) # 30 days
+#' plot(c(-10, 50), c(-0.2,1.2), col=NA, xaxs='i', yaxs='i', xlab='TIME, IN DAYS', ylab='P(t)/BKu')
+#' lines(ti/d2s, r1, lty=2)
+#' lines(ti/d2s, r2, lty=3)
+#' lines(ti/d2s, r3)
+#' text(c(20,22,32), c(0.12,0.35,0.9), parse(text=sprintf("z^2/4*c == %s", c('1 ~ "hour"','1 ~ "day"','30 ~ "days"'))))
+#' 
+#' plot(range(10**lti), c(5e-2, 1), col=NA, xaxs='i', log='y', xlab='TIME, IN DAYS', ylab='P(t)/BKu')
+#' lines(ti/d2s, r1, lty=2)
+#' lines(ti/d2s, r2, lty=3)
+#' lines(ti/d2s, r3)
+r96_drainage_pressure <- function(Times, depth, depth_water=0, diffusiv=1, strain0=-1, B=1, Ku=1, ...){
+  # Equation 50
+  neg <- Times < 0
+  Times[neg] <- NA
+  zsq <- (depth - depth_water)^2
+  evol <- sqrt(zsq / 4 / diffusiv / Times)
+  resp <- erf_re(evol)
+  scaling <- -1 * B * Ku * strain0
+  result <- resp * scaling
+  result[neg] <- 0
+  return(result)
+}
\ No newline at end of file
diff --git a/man/r96_drainage_pressure.Rd b/man/r96_drainage_pressure.Rd
new file mode 100644
index 0000000..8c01af2
--- /dev/null
+++ b/man/r96_drainage_pressure.Rd
@@ -0,0 +1,70 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/roeloffs96.R
+\name{r96_drainage_pressure}
+\alias{r96_drainage_pressure}
+\title{Earthquake-related poroelastic solutions in Roeloffs 1996}
+\usage{
+r96_drainage_pressure(
+  Times,
+  depth,
+  depth_water = 0,
+  diffusiv = 1,
+  strain0 = -1,
+  B = 1,
+  Ku = 1,
+  ...
+)
+}
+\arguments{
+\item{Times}{numeric; times to calculate at \code{depth} [s]}
+
+\item{depth}{numeric; observation depth [m]}
+
+\item{depth_water}{numeric; depth to water table [m]}
+
+\item{diffusiv}{numeric; hydraulic diffusivity [m^2/s]}
+
+\item{strain0}{numeric; static strain change [unscaled strain]}
+
+\item{B}{numeric; Skempton's coefficient B}
+
+\item{Ku}{numeric; undrained bulk modulus}
+
+\item{...}{additional parameters}
+}
+\description{
+Earthquake-related poroelastic solutions in Roeloffs 1996
+}
+\details{
+\code{r96_drainage_pressure} refers to the solution
+given in equation 50 of Roeloffs (1996), the pore pressure
+response to a static strain. The temporal change in pore
+pressure represents water table drainage as a function of time.
+}
+\examples{
+# Remake Fig 16
+d2s <- 86400
+lti <- 10^seq(-3,0, length.out=21)
+ti <- unique(sort(c(0, lti, seq(-10, 50, length.out=201)) * d2s))
+zsq4c <- c(1, 24, 24*30) * d2s / 24 # timescales in hours
+z <- 1
+vert_c <- z^2 / 4 / zsq4c # vertical diffusivity
+r1  <- r96_drainage_pressure(Times=ti, depth=z, diffusiv=vert_c[1]) # 1 hour
+r2  <- r96_drainage_pressure(Times=ti, depth=z, diffusiv=vert_c[2]) # 1 day
+r3  <- r96_drainage_pressure(Times=ti, depth=z, diffusiv=vert_c[3]) # 30 days
+plot(c(-10, 50), c(-0.2,1.2), col=NA, xaxs='i', yaxs='i', xlab='TIME, IN DAYS', ylab='P(t)/BKu')
+lines(ti/d2s, r1, lty=2)
+lines(ti/d2s, r2, lty=3)
+lines(ti/d2s, r3)
+text(c(20,22,32), c(0.12,0.35,0.9), parse(text=sprintf("z^2/4*c == \%s", c('1 ~ "hour"','1 ~ "day"','30 ~ "days"'))))
+
+plot(range(10**lti), c(5e-2, 1), col=NA, xaxs='i', log='y', xlab='TIME, IN DAYS', ylab='P(t)/BKu')
+lines(ti/d2s, r1, lty=2)
+lines(ti/d2s, r2, lty=3)
+lines(ti/d2s, r3)
+}
+\references{
+Roeloffs, E. (1996),
+Poroelastic techniques in the study of earthquake-related hydrologic phenomena,
+\emph{Advances in Geophysics}, 37, 135-195, \url{https://doi.org/10.1016/S0065-2687(08)60270-8}
+}