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