Heterogeneity is the spatial variability in permeability and can be calculated by tow methods:

1- Dykstra-parson coefficient

2- Lorentz Coefficient

1-Dykstra-parson coefficient

This methods implies simple calculations on permeability values in a table form. It starts with calculating the number of data points with greater or equal permeability (Number of samples >=k).This can be used in further calculation of the samples portion (%)>=k).Then, when plotting the latter against permeability in log-normal distribution,a straight trend line should be observed.Finally, we can model the relationship and calculate the heterogeneity index (HI) through mathematical formula.

First, reading the whole data set

data = read.csv("karpur.csv")
head(data)

For simplicity, sorting permeability values in a descending order

data = data[order(data$k.core,decreasing = TRUE), ]
K = data$k.core

Since values are sorted, Number of sample >=k is actually the row indices

sample = c(1: length(K))

Calculating %>=k

k_percent = (sample * 100) / length(K)

when plotting the output, a straight line is observed.

xlab = "portion of Total Samples Having larger or Equal K "
ylab = "permeability (md)"
plot(k_percent, K, log = 'y', xlab = xlab, ylab = ylab, pch =10, cex = 0.5, col = "#001c49")

A linear model is fitted to the data.

log_k = log(K)
model = lm(log_k ~ k_percent)
plot(k_percent,log_k, xlab = xlab, ylab = ylab, pch = 10, cex = 0.5, col = "#001c49")
abline(model, col = 'blue', lwd = 2)

Calculating the HI mathematically.

new_data = data.frame(k_percent = c(50, 84.1))
predicted_values= predict (model, new_data)
heterogenity_index = (predicted_values [1]- predicted_values [2]) / predicted_values [1]
heterogenity_index
##         1 
## 0.2035464

2-Lorentz Method

The Lorenz coefficient of permeability variation is obtained by plotting graph of cumulative k*h versus cumulative phi*h

Reading data set again.

data = read.csv("karpur.csv")

Calculating thickness values for each layer from their tops.

#thickness calculation 
thickness = c(0.5)
for (i in 2:819) {
   h = data$depth[i] - data$depth[i-1]
   thickness = append(thickness, h)
}

data = data[order(data$k.core,decreasing = TRUE), ]

Then,calculate the flow capacity and storage capacity which is easily the multiplication between permeability and porosity with corresponding thickness.Then calculate the cumulative sum of each one using cumsum() function.

KH = thickness * data$k.core
PH = thickness * data$phi.core
cum_sum_kh = cumsum(KH)
cum_sum_ph = cumsum(PH)

Finally we need to calculate the following terms:

∑cum.k.coretotal & ∑cum.phi.coretotal

frac_tot_vol = (cum_sum_ph/100) / max(cum_sum_ph/100) #to convert phi between (0,1)
frac_tot_flow = cum_sum_kh / max(cum_sum_kh)

Now, plotting the percent storage capacity on X-axis and the percent of flow capacity on Y-axis:

plot(frac_tot_vol,frac_tot_flow,  pch = 10, cex = 0.2, col = '#001c49', text(0.4, 0.8, "D"))
text(0,0, "A", pos = 2)
text(1,1, "B", pos = 1)
text(1,0, "C", pos = 2)
abline(0,1, col = 'red' , lwd = 2)

require("DescTools")
## Loading required package: DescTools
library(DescTools)
area = AUC(frac_tot_vol, frac_tot_flow, method = "trapezoid")

The formula of Lorenz Coefficient:

\[ Lorenz\ Coefficient =\frac{area\ ADBA}{area\ ABCA} \]

Lcoff = (area - 0.5) / 0.5
Lcoff
## [1] 0.4524741