data = read.csv('karpur.csv', header = TRUE)
head(data)
##    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   k.core Facies
## 1       -0.033   2.205  33.9000 2442.590     F1
## 2       -0.067   2.040  33.4131 3006.989     F1
## 3       -0.064   1.888  33.1000 3370.000     F1
## 4       -0.053   1.794  34.9000 2270.000     F1
## 5       -0.054   1.758  35.0644 2530.758     F1
## 6       -0.058   1.759  35.3152 2928.314     F1
summary(data)
##      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         
##  Min.   :    0.42   Length:819        
##  1st Qu.:  657.33   Class :character  
##  Median : 1591.22   Mode  :character  
##  Mean   : 2251.91                     
##  3rd Qu.: 3046.82                     
##  Max.   :15600.00

The Lorenz coefficient of permeability variation is obtained by plotting a graph of cumulative kh versus cumulative phi*h, sometimes called a flow capacity plot.

first, let’s calculate the thickness from the depths.

the differences in depths gives the thickness:

d = data$depth
n = length(d)
h = d[2:n]-d[1:n-1] #differences in depths
head(h)
## [1] 0.5 0.5 0.5 0.5 0.5 0.5

we need one more value to have the same number of data samples, since most of the values are 0.5, this what I’m going to add.

h = append(h, 0.5)
head(h)
## [1] 0.5 0.5 0.5 0.5 0.5 0.5

second, lets calculate the cumulative kh, and cumulative phi*h

df        = data.frame(h) # we will use data frames to store the data

df['h']   = h
k         = data$k.core
df['k']   = k
df['phi'] = data$phi.core
df        = df[order(k, decreasing = TRUE),] #sort data such that k is in descending order

df['kh']  = df['k']*df['h'] #multiply k and h
df['sum_kh'] = cumsum(df['kh']) #calculate the cumulative sum
df['sum_kh/total'] = df['sum_kh']/df[n, 'sum_kh'] #divid by the total

# same for phi and h 

df['phih'] = df['phi']*df['h']
df['sum_phih'] = cumsum(df['phih'])
df['sum_phih/total'] = df['sum_phih']/df[n, 'sum_phih']

head(df)
##       h        k     phi       kh   sum_kh sum_kh/total     phih sum_phih
## 809 0.5 15600.00 27.9000 7800.000  7800.00  0.008237623 13.95000 13.95000
## 808 0.5 14225.31 28.5071 7112.657 14912.66  0.015749339 14.25355 28.20355
## 810 0.5 13544.98 27.7563 6772.489 21685.15  0.022901802 13.87815 42.08170
## 807 0.5 13033.53 29.0333 6516.764 28201.91  0.029784192 14.51665 56.59835
## 806 0.5 11841.74 29.5596 5920.872 34122.78  0.036037257 14.77980 71.37815
## 811 0.5 11117.40 27.5865 5558.701 39681.48  0.041907832 13.79325 85.17140
##     sum_phih/total
## 809    0.001241781
## 808    0.002510583
## 810    0.003745967
## 807    0.005038189
## 806    0.006353836
## 811    0.007581664
scatter.smooth(df['sum_phih/total'][1:n,1] ,df['sum_kh/total'][1:n,1],
xlab = "fraction of total volume (phih)",
ylab = "fraction of total flow capacity (kh)")

we need to calculate the area under the above graph, I wrote the fallowing code to find the area, basically this is an implementation of the trapezoid method to fined AUC.

note: I’m familiar with python, and still new to R, writing my own code actully took less time than searching for some libraries that would calculate the area for me.

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

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

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

area = (lower_area + upper_area)/2

print(area)
## [1] 0.7247318

finally we calculate Lorenz coefficient as fallows:

Lorenz coefficient = (Area - 0.5)/0.5

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

The Lorenz coefficient of permeability variation also varies from zero to one. Unfortunately, the Lorenz coefficient is not a unique measure of reservoir heterogeneity. Several different permeability distributions can give the same value of Lorenz coefficient. For log-normal permeability distributions, the Lorenz coefficient is very similar to the Dykstra-Parsons coefficient of permeability variation.