1- Dykstra-Parsons Method
setwd("/Users/khalilahmed/Desktop/User1/Stage 4/RC/Rstudio files")
kdata<- read.csv("karpur.csv",header=TRUE)
##phi core in karpur.csv file is normal number so in excel conveted to fractin under title (phi.core.)
head(kdata)
## depth caliper ind.deep ind.med gamma phi.N R.deep R.med SP
## 1 5667.0 8.685 618.005 569.781 98.823 0.410 1.618 1.755 -56.587
## 2 5667.5 8.686 497.547 419.494 90.640 0.307 2.010 2.384 -61.916
## 3 5668.0 8.686 384.935 300.155 78.087 0.203 2.598 3.332 -55.861
## 4 5668.5 8.686 278.324 205.224 66.232 0.119 3.593 4.873 -41.860
## 5 5669.0 8.686 183.743 131.155 59.807 0.069 5.442 7.625 -34.934
## 6 5669.5 8.686 109.512 75.633 57.109 0.048 9.131 13.222 -39.769
## density.corr density phi.core phi.core. k.core Facies
## 1 -0.033 2.205 33.9000 0.339000 2442.590 F1
## 2 -0.067 2.040 33.4131 0.334131 3006.989 F1
## 3 -0.064 1.888 33.1000 0.331000 3370.000 F1
## 4 -0.053 1.794 34.9000 0.349000 2270.000 F1
## 5 -0.054 1.758 35.0644 0.350644 2530.758 F1
## 6 -0.058 1.759 35.3152 0.353152 2928.314 F1
### Dykstra-Parsons Method
##Sorting base on core permeability
kdata = kdata[order(kdata$k.core, decreasing = TRUE), ]
K = kdata$k.core
##creating another vector containing the number of samples ≥ k
sample = c(1:819)
##percent of samples that have permeability greater or equal to the current K value
kp = (sample * 100) / 820
## logarithmic permeability(y) as function of probability(x)
xlab = "% Higher K "
ylab = "Permeability (md)"
plot(kp, K, log = 'y', xlab = xlab, ylab = ylab, pch = 15, cex = 0.5, col = 'blue')
##use simpe Linear regression Model
kl = log(K)
slr1 = lm(kl ~ kp)
plot(kp,kl, xlab = xlab, ylab = ylab, pch = 15, cex = 0.5, col ='blue')
abline(slr1, col = 'red', lwd = 2)
summary(slr1)
##
## Call:
## lm(formula = kl ~ kp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.8697 -0.2047 0.1235 0.3150 0.4280
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.2584172 0.0377994 244.94 <2e-16 ***
## kp -0.0426137 0.0006549 -65.07 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5404 on 817 degrees of freedom
## Multiple R-squared: 0.8382, Adjusted R-squared: 0.838
## F-statistic: 4234 on 1 and 817 DF, p-value: < 2.2e-16
##find K values at K50 & K84.1
ndata = data.frame(kp = c(50, 84.1))
predictk = predict(slr1, ndata)
## find Heterogeneity index
heterogenity_index = (predictk[1] - predictk[2]) / predictk[1]
heterogenity_index
## 1
## 0.2038692
##heterogeneity index = 0.204 , which indicate slightly heterogeneous, so we can use arithmetic mean for permeability distribution
2- Lorenz Coefficient Method
###Lorenz Coefficient Method
#thickness calculation
thickness = c(0.5)
for (i in 2:819) {
h = kdata$depth[i] - kdata$depth[i-1]
thickness = append(thickness, h)
}
kdata = kdata[order(kdata$k.core, decreasing = TRUE), ]
#calculate storage capacity(cumulative(porosity multiply thickness))
#calculate the flow capacity(cumulative(permeability multiply thickness))
Qh= thickness * kdata$phi.core.
cQh= cumsum(Qh)
kh= thickness * kdata$k.core
ckh= cumsum(kh)
##plotting the percent of storage capacity on X-axis and the percent of flow capacity on Y-axis
pQh = (cQh) / sum(Qh)
pkh = (ckh) / sum(kh)
plot(pQh,pkh, pch = 15, cex = 0.3, col = 'gold', text(0.4, 0.8, "D"), xlim = c(0,1), ylim = c(0,1))
text(0,0, "A", pos = 2)
text(1,1, "B", pos = 1)
text(1,0, "C", pos = 2)
abline(0,1, col = 'blue', lwd = 2)
require(DescTools)
library(DescTools)
area = AUC(pQh, pkh, method="trapezoid")
Lcoff = (area - 0.5) / 0.5
Lcoff