The Dyktra-Parsons coefficient of variation is a widely employed method for evaluating permeability variation. It relies on determining the coefficient value, V, using the formula: V = (k50 - k84.1)/k50. Here, k50 represents the permeability value at %k>=50 (log mean), and k84.1 is the permeability value at %k>=84.1 (one standard deviation). To compute this coefficient, knowledge of k50 and k84.1 is essential. The process involves using Rstudio to plot the frequency distribution of permeability on log-normal probability graph paper. This requires organizing permeability values in descending order and calculating the percentage of samples with permeability greater than each value. To avoid issues with values of zero or 100%, the percentage greater than or equal to a value is normalized by n+1, where n is the number of samples. Before proceeding, the data must be loaded into Rstudio.

data2 = read.csv('D:/karpur.csv')
 k_core = data2$k.core
 depth=data2$depth
 summary(k_core)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##     0.42   657.33  1591.22  2251.91  3046.82 15600.00

Then,we take the permeability and arrange them in descending order.

sorted_k = sort(k_core,decreasing = TRUE)
head(sorted_k)
## [1] 15600.00 14225.31 13544.98 13033.53 11841.74 11117.40

Following the computation of permeability values, the subsequent task is to ascertain the percentage of samples exhibiting permeability greater than each specific value. To circumvent complications associated with values hitting extremes of zero or 100%, the calculated percentage, for values greater than or equal to a particular threshold, is adjusted through normalization using the formula n+1, where n denotes the number of samples.

n = length(sorted_k) 
k_percent = c(1:n)/(n+1) 
head(k_percent)
## [1] 0.001219512 0.002439024 0.003658537 0.004878049 0.006097561 0.007317073

Create a log-normal probability graph by plotting k percent on the x-axis and sorted k on the y-axis.

plot(
  x = k_percent, 
  y = sorted_k,
  type="l",
  xlab = "% >= k",
  ylab = "k",
  log='y',
  lwd=2
  )

As the curve is largely within the 0.1-0.9 interval, we can approximate it effectively using a straight line. This allows for the possibility of performing a straight line fitting on the curve.

nearly_linear_k_values = sorted_k[floor(0.1*n): floor(0.9*n)]
nearly_linear_k_percent = k_percent[floor(0.1*n): floor(0.9*n)]
model = lm(log10(nearly_linear_k_values) ~ nearly_linear_k_percent)
summary(model)
## 
## Call:
## lm(formula = log10(nearly_linear_k_values) ~ nearly_linear_k_percent)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.19128 -0.03880  0.02148  0.04110  0.07176 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              3.885895   0.005010   775.7   <2e-16 ***
## nearly_linear_k_percent -1.461904   0.009112  -160.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05402 on 655 degrees of freedom
## Multiple R-squared:  0.9752, Adjusted R-squared:  0.9751 
## F-statistic: 2.574e+04 on 1 and 655 DF,  p-value: < 2.2e-16
plot(
  x = k_percent, 
  y = sorted_k,
  type="l",
  xlab = "% >= k",
  ylab = "k", 
  log='y',
  lwd=2
  )
abline(model, col="red", lwd=2)

The permeability variation coefficient, V, is derived from the comparison of permeability values at specific percentiles on a graph. To be specific, V is calculated as the difference between the permeability value at %k>=50 (log mean) and the permeability value at %k>=84.1 (one standard deviation), divided by the permeability value at %k>=50, represented as V=(k50−k84.1)/k50.

log_k50_and_84.1 = predict(model, 
  data.frame(nearly_linear_k_percent=c(0.50, 0.841))
  )
print(10^log_k50_and_84.1)
##         1         2 
## 1428.7056  453.3498
k50   = 10^log_k50_and_84.1[1]
k84.1 = 10^log_k50_and_84.1[2]
V = (k50 - k84.1)/k50
print(V)
##        1 
## 0.682685

The Dykstra-Parsons coefficient reveals a permeability coefficient of variation exceeding 0.5, indicating pronounced heterogeneity in the reservoir.

Lorenz coefficient

The Lorenz coefficient for permeability variation is calculated by constructing a flow capacity plot, which entails graphing cumulative kh against cumulative phi*h.

first, let’s calculate the thickness from the depths.the differences in depths gives the thickness:

thickness = depth[2:n]-depth[1:n-1]
summary(thickness)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.5000  0.5000  0.5000  0.5086  0.5000  1.0000
thickness=append(thickness, 0.5)

ta=data.frame(data2)
ta["thickness"]=thickness
j=ta[order(k_core, decreasing = TRUE),]
summary(j)
##      depth         caliper         ind.deep          ind.med       
##  Min.   :5667   Min.   :8.487   Min.   :  6.532   Min.   :  9.386  
##  1st Qu.:5769   1st Qu.:8.556   1st Qu.: 28.799   1st Qu.: 27.892  
##  Median :5872   Median :8.588   Median :217.849   Median :254.383  
##  Mean   :5873   Mean   :8.622   Mean   :275.357   Mean   :273.357  
##  3rd Qu.:5977   3rd Qu.:8.686   3rd Qu.:566.793   3rd Qu.:544.232  
##  Max.   :6083   Max.   :8.886   Max.   :769.484   Max.   :746.028  
##      gamma            phi.N            R.deep            R.med        
##  Min.   : 16.74   Min.   :0.0150   Min.   :  1.300   Min.   :  1.340  
##  1st Qu.: 40.89   1st Qu.:0.2030   1st Qu.:  1.764   1st Qu.:  1.837  
##  Median : 51.37   Median :0.2450   Median :  4.590   Median :  3.931  
##  Mean   : 53.42   Mean   :0.2213   Mean   : 24.501   Mean   : 21.196  
##  3rd Qu.: 62.37   3rd Qu.:0.2640   3rd Qu.: 34.724   3rd Qu.: 35.853  
##  Max.   :112.40   Max.   :0.4100   Max.   :153.085   Max.   :106.542  
##        SP          density.corr          density         phi.core    
##  Min.   :-73.95   Min.   :-0.067000   Min.   :1.758   Min.   :15.70  
##  1st Qu.:-42.01   1st Qu.:-0.016000   1st Qu.:2.023   1st Qu.:23.90  
##  Median :-32.25   Median :-0.007000   Median :2.099   Median :27.60  
##  Mean   :-30.98   Mean   :-0.008883   Mean   :2.102   Mean   :26.93  
##  3rd Qu.:-19.48   3rd Qu.: 0.002000   3rd Qu.:2.181   3rd Qu.:30.70  
##  Max.   : 25.13   Max.   : 0.089000   Max.   :2.387   Max.   :36.30  
##      k.core            Facies            thickness     
##  Min.   :    0.42   Length:819         Min.   :0.5000  
##  1st Qu.:  657.33   Class :character   1st Qu.:0.5000  
##  Median : 1591.22   Mode  :character   Median :0.5000  
##  Mean   : 2251.91                      Mean   :0.5085  
##  3rd Qu.: 3046.82                      3rd Qu.:0.5000  
##  Max.   :15600.00                      Max.   :1.0000

Now we will make table contain (thickness,phi.core,k.core):

ta['kh']= j$k.core * j$thickness
ta['sum_kh']=cumsum(ta['kh'])
ta['sum_kh/total']=ta['sum_kh']/ta[n, 'sum_kh']

ta['phih']           = j$phi.core * j$thickness
ta['sum_phih']       = cumsum(ta['phih'])
ta['sum_phih/total'] = ta['sum_phih']/ta[n, 'sum_phih']

plot(ta['sum_phih/total'][1:n,1], ta['sum_kh/total'][1:n,1], xlab = "fraction of total volume (phih)", 
  ylab = "fraction of total flow capacity (kh)")

x  = c(ta['sum_phih/total'][1:n,1])
dx = x[2:n] - x[1:n-1]

y1 = c(ta['sum_kh/total'][1:n,1])[1:n-1]
y2 = c(ta['sum_kh/total'][1:n,1])[2:n]

lower_area = sum(y1*dx) #Reduce("+", y1*dx)
upper_area = sum(y2*dx) #Reduce("+", y2*dx)

AUC = (lower_area + upper_area)/2

print(AUC)
## [1] 0.7247318

Lorenz coefficient=(AUC−0.5)/0.5

Lorenz_coefficient = (AUC - 0.5)/0.5
print(Lorenz_coefficient)
## [1] 0.4494636