Arps Decline Analysis in R

Here’s a quick example to show how to use R to fit a best-fit line to the most linear portion of a data series. This example fits two common Arps models to well production data in order to estimate ultimate recovery.

If you haven’t already done so, download R from http://www.r-project.org/ and install it.

Once you’ve declared a variable for rate (“q” in the attached code, with units of Mcf/d) and a variable for cumulative production (“Q” in the attached code, with units of MMcf), it only takes a couple more lines of code to calculate exponential and harmonic declines. A final rate of 75 Mcf/d was used.


require(MASS)
abandon_rate <- 75 URR <- c(0,0) q <- c(721.2, 2143.4, 1114.6, 823.7, 620.1, 554.9, 692.9, 597, 592.7, 584.5, 737, 609.8, 355.5, 318.4, 567, 488, 455.8, 148.8, 711.5, 76, 170.3, 385.7, 445.6, 500.3, 385.1, 341.7, 349.1, 331, 244.1, 198.8, 315.8, 344.8, 308, 310.9, 330.2, 307.3, 290.6, 273.3, 270.1, 275.7, 275, 264.2, 262.8, 260.3, 246.5, 235.6, 254.1, 234.4, 252.7, 235.7, 230.9, 234.3, 237.1, 203, 202.9, 258.2, 219.6, 213, 208, 200.7, 209.1, 214.3, 209.5, 207.6, 213.3, 203, 92.4, 99.5, 131.7, 204.5, 198.6, 211.4, 195.8, 166, 180.1, 200.1, 123.8, 161.1, 173.2, 132.3, 10.4, 162.8) Q <- c(21.6, 88.1, 121.5, 147.1, 166.3, 182.9, 204.4, 222.3, 240.7, 258.8, 279.4, 298.3, 309, 318.9, 335.9, 351, 365.2, 369.6, 391.7, 394, 399.2, 411.2, 423.7, 439.2, 450.7, 461.3, 471.8, 482.1, 489.6, 495.6, 505.4, 515.7, 525.3, 534.9, 544.1, 553.7, 562.4, 570.9, 579, 587.5, 596, 604, 612.1, 619.9, 627.6, 634.9, 642.2, 649.5, 657.1, 664.4, 671.3, 678.6, 685.9, 692, 698.3, 706.1, 712.9, 719.5, 725.3, 731.5, 737.8, 744.4, 750.7, 757.1, 763.8, 769.9, 772.7, 775.7, 779.8, 786.1, 791.7, 798.2, 804.1, 809.3, 814.7, 820.9, 824.7, 829.5, 834.9, 838.9, 839.2, 839.2) decline_exp <- lmsreg(Q[q>0], q[q>0])
URR[1] <- (abandon_rate-decline_exp$coefficients[1])/decline_exp$coefficients[2] decline_hyp <- lmsreg(Q[q>0], log(q[q>0]))
URR[2] <- (log(abandon_rate)-decline_hyp$coefficients[1])/decline_hyp$coefficients[2] png(file="arps_decline_exp.png",width=400,height=350) if(URR[1]>0) xmax<-1000*ceiling(URR[1]/1000) else xmax <- 1000*ceiling(max(Q/1000000)) plot(Q[q>0],q[q>0],xlim=c(0,xmax),xlab="Cumulative Gas Production (MMcf)",ylab="Calendar Day Gas Rate (Mcf/d)",col="00000016",cex=1.0,main=paste("Exponential URR = ",formatC(50*round(URR[1]/50), big.mark = ",", format = "d")," MMcf",sep=""),type = "p")
abline(decline_exp, col="blue")
abline(h=abandon_rate,col="gray")
dev.off()

png(file="arps_decline_harm.png",width=400,height=350)
if(URR[2]>0) xmax<-1000*ceiling(URR[2]/1000) else xmax <- 1000*ceiling(max(Q/1000000)) plot(Q[q>0],log(q[q>0]),ylim=log(c(10,10000)),xlim=c(0,xmax),xlab="Cumulative Gas Production (MMcf)",ylab="Calendar Day Gas Rate (Mcf/d)", yaxt="n",col="00000016",cex=1.0,main=paste("Harmonic URR = ",formatC(50*round(URR[2]/50), big.mark = ",", format = "d")," MMcf",sep=""))
axis(2,at=log(10^(1:4)),labels=formatC(c(10^(1:4)), big.mark = ",", format = "d"))
axis.at <- 10^(1:4) axis(2, at = log(1:10 * rep(axis.at[-1] / 10, each = 10)),tcl = -0.5, labels = FALSE,tck=-0.01) abline(decline_hyp, col="blue") abline(h=log(abandon_rate),col="gray") dev.off()

The function lmsreg (from the MASS package) is used to fit a linear relationship (ie a best-fit line) through the good data (the data that is best represented by a straight line). Data that does not fit a straight-line relationship, or is dramatically different than the main trend, is given very little weight. The end results are very similar to what a human would assign - the code given above produces the following plots for exponential and harmonic declines.