Added time-varying pore pressure field

DESCRIPTION

@@ -17,7 +17,8 @@ Depends:
     sp,
     maptools,
     plyr,
-    reshape2
+    reshape2,
+    abind
 Suggests:
     maps
 LinkingTo: Rcpp

NAMESPACE

@@ -18,6 +18,7 @@ export(rigid_hydraulic_diffusivity)
 export(segall85)
 export(surface_displacement)
 export(timevarying_fluidmass)
+export(timevarying_porepressure)
 export(timevarying_surface_displacement)
 exportPattern("^[[:alpha:]]+")
 importFrom(Rcpp,evalCpp)

NOBUILD/Sandbox/expandgrid.R

@@ -0,0 +1,11 @@
+x <- 1:3
+z <- c(10,20,30,40)
+grd <- as.matrix(expand.grid(VarX=x,VarZ=z))
+
+testfun <- function(XY){
+  sum(XY)
+}
+
+Res <- apply(grd, 1, testfun)
+
+matrix(Res, ncol=length(z))

NOBUILD/Sandbox/segall.R

@@ -1,4 +1,5 @@
 library(deform)
+library(kook)
 library(TeachingDemos)
 
 .x. <- sort(unique(c(-7:-3, seq(-3.90,3.90,by=0.15), 3:7)))
@@ -33,6 +34,7 @@ F2 <- function(){
 # Figs 7,8
 mxx <- 50
 .x.km. <- sort(unique(c((-1*mxx):-3, seq(-2.90,2.90,by=0.1), 3:mxx)))
+.z.km. <- sort(unique(c(seq(0,3,by=0.25),seq(3,12,by=0.75))))
 yr <- 365*86400
 .time. <- seq(2,10,by=2)*10*yr
 .Vdot. <- 2e6/yr # volume rate m^3/yr to m^3/s
@@ -41,9 +43,12 @@ yr <- 365*86400
 .B. <- 0.6 # Skemptons coeff
 .c. <- 0.1 # hydraulic diffusivity m^2/s
 .Sources.x. <- 1e3*c(0)
+.TwoSources.x. <- 1e3*c(0,20)
 # for mass computation
 .t. <- 100 # thickness
 .phi. <- 0.2 #porosity
+# for pressure computation
+.mu. <- 5.6 #GPa -- shear modulus
 
 zz2 <- timevarying_surface_displacement(.x.km.*1e3, .time., .Vdot., .B., .L., .D., .c., Pt.Sources.x=.Sources.x.)
 zz2t <- apply(zz2, 2, function(.z.) matrix(Tilt(.x.km.*1e3, z=.z.)$ztilt))
@@ -53,6 +58,15 @@ F3 <- function(){
   matplot(.x.km., zz2*1e3, type="l", col="black", main="Subsidence, mm, Segall 1985, Fig 8B", sub=Sys.time())
 }
 
+F3t <- function(){
+  #matplot(.time./yr, t(zz2t), type="l")
+  matplot(.x.km., zz2t, type="l", col="black")
+}
+
+#try(F3())
+#try(F3t())
+
+
 zz3 <- timevarying_fluidmass(.x.km.*1e3, .time., .Vdot., .L., .t., .c., phi.=.phi.)
 
 F4 <- function(){
@@ -62,12 +76,33 @@ F4 <- function(){
 
 #try(F4())
 
+zzp <- timevarying_porepressure(.x.km.*1e3, .z.km.*1e3, .time., .Vdot., .B., .L., .D., .c., .t., .mu., Pt.Sources.x=.TwoSources.x.)
 
-F3t <- function(){
-  #matplot(.time./yr, t(zz2t), type="l")
-  matplot(.x.km., zz2t, type="l", col="black")
+F5 <- function(do.log=FALSE){
+  #matplot(.time./yr, t(zz3)*1e2, type="l")
+  X<- zzp[,,length(.time.)]
+  if (do.log) X <- log10(abs(X))
+  matplot(x=.x.km., X, col=NA)
+  aaply(zzp, 3, .fun = function(X) {
+    if (do.log) X <- log10(abs(X))
+    matplot(x=.x.km., X, type="l", add=TRUE)
+    return("x")
+    })
+  invisible()
 }
 
-try(F3())
-#try(F3t())
+try(F5())
+
+
+F5c <- function(){
+  #matplot(.time./yr, t(zz3)*1e2, type="l")
+  #contour(x=.x.km., y=.z.km., zzp[,,length(.time.)], col=NA)
+  aaply(zzp, 3, .fun = function(X) {
+    image(x=.x.km., y=.z.km., X, ylim=c(12,0), col = brewerRamp())
+    contour(x=.x.km., y=.z.km., X, ylim=c(12,0), add = TRUE)
+    return("x")
+    })
+  invisible()
+}
 
+#try(F5c())

R/segall85.R

@@ -88,6 +88,68 @@ timevarying_surface_displacement <- function(x, Time, Vdot., B., L., D., HD., nu
   apply(X=matrix(Time), MARGIN=1, FUN=.FUN, Xi.sources=Sources)
 }
 
+#' @rdname segall85
+#' @export
+timevarying_porepressure <- function(x, z, Time, Vdot., B., L., D., HD., t., mu.gpa., nu.=1/4, nuu.=1/3, Pt.Sources.x=0, x.lim=round(max(abs(x),na.rm=TRUE))/1){
+  # Time varying p.p. associated with fluid extraction
+  # segall85 eq 28
+  #
+  sc <- 1e9
+  mu. <- sc * mu.gpa.
+  #
+  Sources <- matrix(Pt.Sources.x, ncol=1)
+  #
+  .FUN <- function(ti, zi=0, Xi.sources){
+    message(paste("Time:", ti/365/86400, "years"))
+    # time varying component (scaling)
+    tvar1 <- -2 * mu. * (1 + nuu.)^2 * B.^2 * Vdot. * sqrt(ti / HD.) / (9 * L. * t.)
+    #
+    # Function that returns the integral value at
+    # a position Xi away from the source at dXi
+    .fun <- function(Xi., dXi.=0, X.=x, Z.=0, ti.=ti){
+      # from the spatio-temporal component
+      ti.d <- Xi.^2 / (4 * HD. * ti.)
+      C1 <- ierfc2(sqrt(ti.d))
+      #
+      xmxi <- (X. - Xi. - dXi.)^2
+      zpd <- (Z. + D.)
+      C2A <- zpd / (zpd^2 + xmxi^2)
+      zpdt <- zpd + t.
+      C2B <- -1 * zpdt / (zpdt^2 + xmxi^2)
+      sc * (C2A + C2B) * C1
+    }
+    #
+    # function used to integrate over Xi at different sources
+    .xi.integ2D <- function(dxi, ti.=ti){
+      message(paste("\t2D integration of source at:", dxi/1e3, "km"))
+      grd <- as.matrix(expand.grid(VarX=x, VarZ=z))
+      ires <- apply(grd, 1, FUN=function(xxzz){
+        #message(paste(xxzz, collapse="//"))
+        zi <- try(integrate(.fun, -1*x.lim, x.lim, X.=xxzz[1], Z.=xxzz[2], dXi.=dxi, subdivisions = 150L, rel.tol=.Machine$double.eps^0.5))
+        if (inherits(zi,"try-error")){
+          #message(paste("\t==== ==== ====>  shrink integration limits to below",x.lim))
+          message(paste(xxzz, collapse=" x // z "))
+          NA
+        } else {
+          #message(paste("integration",zi$message))
+          ti.d <- xxzz[1]^2 / (4 * HD. * ti.)
+          C1 <- (1 - 2*nu.) * ierfc2(sqrt(ti.d)) / ((nuu. - nu.) * (1 - 2*nuu.))
+          C2 <- -2/(pi * (1 - nuu.))
+          C1 + C2 * zi$value
+        }
+      })
+      matrix(ires, ncol=length(z), byrow = FALSE)
+    }
+    # superpostion of all sources
+    xvar.t2s <- apply(X = abind(lapply(X = Sources, FUN = .xi.integ2D), along=3), MARGIN = c(1,2), FUN = sum, na.rm = TRUE)
+    # combine and return
+    res <- tvar1 * xvar.t2s
+    return(-1*res/sc)
+  }
+  # apply FUN through time
+  abind(lapply(X = Time, FUN = .FUN, Xi.sources=Sources), along=3)
+}
+
 #' Simple numerical deformation estimates
 #' @name Simple-deformation
 #' @param x numeric; the spatial vector in real units (not an index!)

man/segall85.Rd

@@ -4,6 +4,7 @@
 \alias{segall85}
 \alias{surface_displacement}
 \alias{timevarying_fluidmass}
+\alias{timevarying_porepressure}
 \alias{timevarying_surface_displacement}
 \title{Surface deformation associated with fluid withdrawl}
 \usage{
@@ -17,6 +18,10 @@ timevarying_fluidmass(x, Time, Vdot., L., t., HD., phi.)
 
 timevarying_surface_displacement(x, Time, Vdot., B., L., D., HD., nuu. = 1/3,
   Pt.Sources.x = 0, x.lim = 2 * round(max(abs(x), na.rm = TRUE)))
+
+timevarying_porepressure(x, z, Time, Vdot., B., L., D., HD., t., mu.gpa.,
+  nu. = 1/4, nuu. = 1/3, Pt.Sources.x = 0, x.lim = round(max(abs(x),
+  na.rm = TRUE))/1)
 }
 \description{
 Surface deformation associated with fluid withdrawl