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

1- Dykstra-Parson Coefficient

2- Lorentz Coefficient

1-Dykstra-Parson Coefficient

Let’s load the data

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

Sorting permeability values in a descending order and finding samples >= k portions

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

#Calculating Number of Samples >= k
sample = c(1: length(K))

# Calculating % >= k
k_percent = (sample * 100) / length(K)

Straight line can be shown

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 = 'red', lwd = 2)

Calculating the Heterogeneity Index

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 

So, the heterogeneity index = 0.2 , which means slightly heterogeneous reservoir, which can be approximated by homogeneous numerical model with minimal error

2-Lorentz Method

The Lorenz coefficient of permeability variation is obtained by plotting graph of cumulative kh 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 of 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")
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

is 0.45, in Lorenz method if the value is closer to zero, it considers more heterogeneous. So the current value considers as heterogeneous reservoir.

LS0tDQp0aXRsZTogIkhldGVyb2dlbmVpdHkgUXVhbnRpZmljYXRpb24gYnkgRHlrc3RyYS1QYXJzb24gYW5kIExvcmVudHogQ29lZmZpY2llbnRzIE1ldGhvZHMiDQphdXRob3I6ICJJYnJhaGltIEtoYWxlZWwgSWJyYWhpbSINCmRhdGU6ICJEZWNlbWJlciAxNiwgMjAyMyINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNCkhldGVyb2dlbmVpdHkgaXMgdGhlIHNwYXRpYWwgdmFyaWFiaWxpdHkgaW4gcGVybWVhYmlsaXR5IGFuZCBjYW4gYmUgY2FsY3VsYXRlZCBieSB0d28gbWV0aG9kczoNCg0KMS0gRHlrc3RyYS1QYXJzb24gQ29lZmZpY2llbnQNCg0KMi0gTG9yZW50eiBDb2VmZmljaWVudA0KDQojIDEtRHlrc3RyYS1QYXJzb24gQ29lZmZpY2llbnQNCg0KTGV0J3MgbG9hZCB0aGUgZGF0YQ0KDQpgYGB7cn0NCmRhdGEgPSByZWFkLmNzdigna2FycHVyLmNzdicpDQpoZWFkKGRhdGEpDQpgYGANCg0KU29ydGluZyBwZXJtZWFiaWxpdHkgdmFsdWVzIGluIGEgZGVzY2VuZGluZyBvcmRlciBhbmQgZmluZGluZyBzYW1wbGVzIFw+PSBrIHBvcnRpb25zDQoNCmBgYHtyfQ0KZGF0YSA9IGRhdGFbb3JkZXIoZGF0YSRrLmNvcmUsIGRlY3JlYXNpbmc9VFJVRSksIF0NCksgPSBkYXRhJGsuY29yZQ0KDQojQ2FsY3VsYXRpbmcgTnVtYmVyIG9mIFNhbXBsZXMgPj0gaw0Kc2FtcGxlID0gYygxOiBsZW5ndGgoSykpDQoNCiMgQ2FsY3VsYXRpbmcgJSA+PSBrDQprX3BlcmNlbnQgPSAoc2FtcGxlICogMTAwKSAvIGxlbmd0aChLKQ0KYGBgDQoNClN0cmFpZ2h0IGxpbmUgY2FuIGJlIHNob3duDQoNCmBgYHtyfQ0KeGxhYiA9ICJQb3J0aW9uIG9mIFRvdGFsIFNhbXBsZXMgSGF2aW5nIExhcmdlciBvciBFcXVhbCBLICINCnlsYWIgPSAiUGVybWVhYmlsaXR5IChtZCkiDQpwbG90KGtfcGVyY2VudCwgSywgbG9nID0gICd5JywgeGxhYiA9IHhsYWIsIHlsYWIgPSB5bGFiLCBwY2ggPSAxMCwgY2V4ID0gMC41LCBjb2wgPSAiIzAwMWM0OSIpDQpgYGANCg0KQSBsaW5lYXIgbW9kZWwgaXMgZml0dGVkIHRvIHRoZSBkYXRhLg0KDQpgYGB7cn0NCmxvZ19rID0gbG9nKEspDQptb2RlbCA9IGxtKGxvZ19rIH4ga19wZXJjZW50KQ0KcGxvdChrX3BlcmNlbnQsbG9nX2ssIHhsYWIgPSB4bGFiLCB5bGFiID0geWxhYiwgcGNoID0gMTAsIGNleCA9IDAuNSwgY29sID0gIiMwMDFjNDkiKQ0KYWJsaW5lKG1vZGVsLCBjb2wgPSAncmVkJywgbHdkID0gMikNCmBgYA0KDQpDYWxjdWxhdGluZyB0aGUgSGV0ZXJvZ2VuZWl0eSBJbmRleA0KDQpgYGB7cn0NCm5ld19kYXRhID0gZGF0YS5mcmFtZShrX3BlcmNlbnQgID0gYyg1MCwgODQuMSkpDQpwcmVkaWN0ZWRfdmFsdWVzID0gcHJlZGljdChtb2RlbCwgbmV3X2RhdGEpDQpoZXRlcm9nZW5pdHlfaW5kZXggPSAocHJlZGljdGVkX3ZhbHVlc1sxXSAtIHByZWRpY3RlZF92YWx1ZXNbMl0pIC8gcHJlZGljdGVkX3ZhbHVlc1sxXQ0KaGV0ZXJvZ2VuaXR5X2luZGV4DQpgYGANCg0KU28sIHRoZSBoZXRlcm9nZW5laXR5IGluZGV4ID0gMC4yICwgd2hpY2ggbWVhbnMgc2xpZ2h0bHkgaGV0ZXJvZ2VuZW91cyByZXNlcnZvaXIsIHdoaWNoIGNhbiBiZSBhcHByb3hpbWF0ZWQgYnkgaG9tb2dlbmVvdXMgbnVtZXJpY2FsIG1vZGVsIHdpdGggbWluaW1hbCBlcnJvcg0KDQojIDItTG9yZW50eiBNZXRob2QNCg0KVGhlIExvcmVueiBjb2VmZmljaWVudCBvZiBwZXJtZWFiaWxpdHkgdmFyaWF0aW9uIGlzIG9idGFpbmVkIGJ5IHBsb3R0aW5nIGdyYXBoIG9mIGN1bXVsYXRpdmUgaypoIHZlcnN1cyBjdW11bGF0aXZlIHBoaSogaA0KDQpSZWFkaW5nIGRhdGEgc2V0IGFnYWluLg0KDQpgYGB7cn0NCmRhdGEgPSByZWFkLmNzdigia2FycHVyLmNzdiIpDQpgYGANCg0KQ2FsY3VsYXRpbmcgdGhpY2tuZXNzIHZhbHVlcyBmb3IgZWFjaCBsYXllciBmcm9tIHRoZWlyIHRvcHMuDQoNCmBgYHtyfQ0KI3RoaWNrbmVzcyBjYWxjdWxhdGlvbg0KdGhpY2tuZXNzID0gYygwLjUpDQpmb3IgKGkgaW4gMjo4MTkpIHsNCiAgaCA9IGRhdGEkZGVwdGhbaV0gLSBkYXRhJGRlcHRoW2ktMV0NCiAgdGhpY2tuZXNzID0gYXBwZW5kKHRoaWNrbmVzcywgaCkNCn0NCg0KZGF0YSA9IGRhdGFbb3JkZXIoZGF0YSRrLmNvcmUsIGRlY3JlYXNpbmcgPSBUUlVFKSwgXQ0KYGBgDQoNClRoZW4sIGNhbGN1bGF0ZSB0aGUgZmxvdyBjYXBhY2l0eSBhbmQgc3RvcmFnZSBjYXBhY2l0eSB3aGljaCBpcyBlYXNpbHkgdGhlIG11bHRpcGxpY2F0aW9uIGJldHdlZW4gcGVybWVhYmlsaXR5IGFuZCBwb3Jvc2l0eSB3aXRoIGNvcnJlc3BvbmRpbmcgdGhpY2tuZXNzLiBUaGVuIGNhbGN1bGF0ZSB0aGUgY3VtdWxhdGl2ZSBzdW0gb2YgZWFjaCBvbmUgdXNpbmcgY3Vtc3VtKCkgZnVuY3Rpb24uDQoNCmBgYHtyfQ0KS0ggPSB0aGlja25lc3MgKiBkYXRhJGsuY29yZQ0KUEggPSB0aGlja25lc3MgKiBkYXRhJHBoaS5jb3JlDQpjdW1fc3VtX2toID0gY3Vtc3VtKEtIKQ0KY3VtX3N1bV9waCA9IGN1bXN1bShQSCkNCmBgYA0KDQpGaW5hbGx5IHdlIG5lZWQgdG8gY2FsY3VsYXRlIHRoZSBmb2xsb3dpbmcgdGVybXM6DQoNCuKIkWN1bS5rLmNvcmV0b3RhbCAmIOKIkWN1bS5waGkuY29yZXRvdGFsDQoNCmBgYHtyfQ0KZnJhY190b3Rfdm9sID0gKGN1bV9zdW1fcGgvMTAwKSAvIG1heChjdW1fc3VtX3BoLzEwMCkgI3RvIGNvbnZlcnQgcGhpIGJldHdlZW4gKDAsMSkNCmZyYWNfdG90X2Zsb3cgPSBjdW1fc3VtX2toIC8gbWF4KGN1bV9zdW1fa2gpDQpgYGANCg0KTm93LCBwbG90dGluZyB0aGUgcGVyY2VudCBvZiBzdG9yYWdlIGNhcGFjaXR5IG9uIFgtYXhpcyBhbmQgdGhlIHBlcmNlbnQgb2YgZmxvdyBjYXBhY2l0eSBvbiBZLWF4aXM6DQoNCmBgYHtyfQ0KcGxvdChmcmFjX3RvdF92b2wsZnJhY190b3RfZmxvdywgcGNoID0gMTAsIGNleCA9IDAuMiwgY29sID0gJyMwMDFjNDknLCB0ZXh0KDAuNCwgMC44LCAiRCIpKQ0KdGV4dCgwLDAsICJBIiwgcG9zID0gMikNCnRleHQoMSwxLCAiQiIsIHBvcyA9IDEpDQp0ZXh0KDEsMCwgIkMiLCBwb3MgPSAyKQ0KYWJsaW5lKDAsMSwgY29sID0gJ3JlZCcsIGx3ZCA9IDIpDQpgYGANCg0KYGBge3J9DQpyZXF1aXJlKCJEZXNjVG9vbHMiKQ0KYGBgDQoNCmBgYHtyfQ0KbGlicmFyeShEZXNjVG9vbHMpDQoNCmFyZWEgPSBBVUMoZnJhY190b3Rfdm9sLCBmcmFjX3RvdF9mbG93LCBtZXRob2Q9InRyYXBlem9pZCIpDQpgYGANCg0KVGhlIGZvcm11bGEgb2YgTG9yZW56IENvZWZmaWNpZW50Og0KDQokJA0KTG9yZW56XCBDb2VmZmljaWVudCA9IFxmcmFje2FyZWFcIEFEQkF9e2FyZWFcIEFCQ0F9DQokJA0KDQpgYGB7cn0NCkxjb2ZmICA9IChhcmVhIC0gMC41KSAvIDAuNQ0KYGBgDQoNCmBgYHtyfQ0KTGNvZmYNCmBgYA0KDQppcyAwLjQ1LCBpbiBMb3JlbnogbWV0aG9kIGlmIHRoZSB2YWx1ZSBpcyBjbG9zZXIgdG8gemVybywgaXQgY29uc2lkZXJzIG1vcmUgaGV0ZXJvZ2VuZW91cy4gU28gdGhlIGN1cnJlbnQgdmFsdWUgY29uc2lkZXJzIGFzIGhldGVyb2dlbmVvdXMgcmVzZXJ2b2lyLg0K