Material Balance - Gas Condensate Reservoirs"

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Examples

Example 1: Lean Gas Condensate Reservoir [@Walsh2003]

Part I: History Match

library(Rmbal)
library(Rrelperm)
library(pracma)
library(minpack.lm)
library(ggplot2)
library(dplyr)
library(magrittr)

p_pvt <- c(3700, 3650, 3400, 3100, 2800, 2500, 2200, 1900, 1600, 1300, 1000, 700, 
           600, 400)

Bo <- c(10.057, 2.417, 2.192, 1.916, 1.736, 1.617, 1.504, 1.416, 1.326, 1.268, 1.205, 
        1.149, 1.131, 1.093)

Rv <- c(84.11765, 84.11765, 70.5, 56.2, 46.5, 39.5, 33.8, 29.9, 27.3, 25.5, 25.9, 
        28.3, 29.8, 33.5) / 1e6

Rs <- c(11566, 2378, 2010, 1569, 1272, 1067, 873, 719, 565, 461, 349, 249, 218, 141)

Bg <- c(0.87, 0.88, 0.92, 0.99, 1.08, 1.20, 1.35, 1.56, 1.85, 2.28, 2.95, 4.09, 
        4.68, 6.53) / 1000

cw <- 3e-6

Bwi <- 10.05

Bw <- Bwi * exp(cw * (p_pvt[1] - p_pvt))

muo <- c(0.0612, 0.062, 0.1338, 0.1826, 0.2354, 0.3001, 0.3764, 0.4781, 0.6041, 
         0.7746, 1.0295, 1.358, 1.855, 2.500)

mug <- c(0.0612, 0.062, 0.0554, 0.0436, 0.0368, 0.0308, 0.0261, 0.0222, 0.0191, 
         0.0166, 0.0148, 0.0135, 0.0125, 0.0115)

muw <- rep(0.25, length(p_pvt))

pvt_table <- data.frame(p = p_pvt, Bo = Bo, Rs = Rs, Rv = Rv, Bg = Bg, 
                               Bw = Bw, muo = muo, mug = mug, muw = muw)

p <- c(3700, 3650, 3400, 3100, 2800, 2500, 2200, 1900, 1600, 1300, 1000, 700, 
           600)

We <- rep(0, length.out = length(p))

Np <- c(0, 28.6, 93, 231, 270, 379, 481, 517.2, 549, 580, 675, 755, 803) *1e3

Gp <- c(0, 0.34, 1.2, 3.3, 4.3, 6.6, 9.1, 10.5, 12, 12.8, 16.4, 19.1, 20.5) * 1e9

Wp <- rep(0, length.out = length(p))

Wi <- rep(0, length.out = length(p))

wf <- rep(1, length.out = length(p))

mbal_optim_gas_lst <- mbal_optim_param_gas(input_unit = "Field", output_unit = "Field",  
                                          unknown_param = "G", aquifer_model = NULL, 
                                          phi = 0.1, swi = 0.2, Np = Np, 
                                          Gp = Gp, Wp = Wp, Wi = Wi, We = We, 
                                          pd = 3650, p = p, pvt = pvt_table, M = 0, 
                                          cf = 2e-6, wf = wf, sgrw = 0.15)

time_lst <- mbal_time(c(1:length(p)), "year")

# a number of plots will be automatically generated for quality check

optim_results <- mbal_optim_gas(mbal_optim_gas_lst, time_lst)

glimpse(optim_results)

Part II: Reservoir Performance

mbal_results <- mbal_perform_gas(optim_results, time_lst)

mbal_results

Part III: Reservoir Forecast

# Step I: generating a set of pseudo relative permeability curves using 
# laboratory 'Kr' values

sg_lab <- c(0.05, 0.152, 0.248, 0.352, 0.448, 0.552, 0.65)

krg_lab <- c(0, 0.05, 0.09, 0.18, 0.3, 0.5, 1)

kro_lab <- c(1, 0.6, 0.35, 0.13, 0.04, 0.01, 0)

swcrit_lab <- 0.2

sgcrit_lab <- 0.05

sorgr_lab <- 0.15

fun_kr <- function(x, swcrit, sgcrit, sorg, sg, krg, kro) {

  kr_est <- Rrelperm::kr2p_gl(SWCON = swcrit, SOIRG = sorg, SORG = sorg,
                               SGCON = sgcrit, SGCRIT = sgcrit, KRGCL = 1,
                               KROGCG = 1, NG = x[1], NOG = x[2], NP = 101)

   krg_est_sub <- approx(x = kr_est[,1], y = kr_est[,3], xout = sg, rule = 2)$y

   kro_est_sub <- approx(x = kr_est[,1], y = kr_est[,4], xout = sg, rule = 2)$y

   error <- (krg - krg_est_sub) ^ 2 + (kro - kro_est_sub) ^ 2

   return(error)
}

par <- c(2, 2)

opt_results <- minpack.lm::nls.lm(par = par, fn = fun_kr, swcrit = swcrit_lab,
                                  sgcrit = sgcrit_lab, sorg = sorgr_lab, sg = sg_lab,
                                  krg = krg_lab, kro = kro_lab,
                                  lower = c(0.1,0.1), upper = c(10,10))

sol <- opt_results$par

sol

rel_perm <- as.data.frame(Rrelperm::kr2p_gl(SWCON = swcrit_lab, SOIRG = sorgr_lab,
                                            SORG = sorgr_lab, SGCON = sgcrit_lab,
                                            SGCRIT = sgcrit_lab, KRGCL = 1,
                                            KROGCG = 1, NG = sol[1], NOG = sol[2],
                                            NP = 101))

colnames(rel_perm) <- c("Sg", "Sl", "Krg", "Krog")

p_forecast <- p

wf_forecast <- wf

time_lst_forecast <- mbal_time(c(1:length(p_forecast)), "year")

forecast_lst <- mbal_forecast_param_gas(input_unit = "Field", output_unit = "Field",
                                        G = 2.41e10, phi = 0.1, swi = 0.2,
                                        pd = 3650, p = p_forecast, pvt = pvt_table,
                                        M = 0, cf = 2e-6, wf = wf_forecast,
                                        rel_perm = rel_perm)

glimpse(forecast_lst)

forecast_results <- mbal_forecast_gas(forecast_lst, time_lst_forecast)

forecast_results

p1 <- forecast_results %>% ggplot(aes(`P (psia)`, SOg, color = "Forecast")) +
  geom_point(size = 3) +
  geom_point(data = mbal_results, aes(`P (psia)`, SOg, color = "Field"))+
  scale_color_manual(name="Data", values=c("Forecast" = "red", "Field" = "black")) +
  ggtitle("Condensate Saturation Plot") +
  theme_bw()

p1


p2 <- forecast_results %>% ggplot(aes(`P (psia)`, `GOR (SCF/STB)`, color = "Forecast")) +
  geom_point(size = 3) +
  geom_point(data = mbal_results, aes(`P (psia)`, `GOR (SCF/STB)`, color = "Field")) +
  scale_color_manual(name="Data", values=c("Forecast" = "red", "Field" = "black")) +
  ggtitle("GOR Plot") +
  theme_bw()

p2


p3 <- forecast_results %>% ggplot(aes(`P (psia)`, `RF_oil`, color = "Forecast")) +
  geom_point(size = 3) +
  geom_point(data = mbal_results, aes(`P (psia)`, `RF_oil`, color = "Field")) +
  scale_color_manual(name="Data", values=c("Forecast" = "red", "Field" = "black")) +
  ggtitle("Condensate Recovery Plot") +
  theme_bw()

p3


p4 <- forecast_results %>% ggplot(aes(`P (psia)`, `Liq_volume`, color = "Forecast")) +
  geom_point(size = 3) +
  scale_color_manual(name="Data", values=c("Forecast" = "red")) +
  ggtitle("CCE Liquid Volume Plot") +
  theme_bw()

p4


p5 <- forecast_results  %>%
  tidyr::pivot_longer(cols = c("Igd", "Ifwd"), names_to = "Drive Mechanism",
                      values_to = "Fraction", values_drop_na = TRUE) %>%
  ggplot(aes(`P (psia)`, Fraction, fill = `Drive Mechanism`)) +
  geom_area() +
  ggtitle("Energy Plot") +
  theme_bw()

p5

Example 2: Lean Gas Condensate Reservoir [@Walsh1994a]

Part I: Reservoir Forecast

library(Rmbal)
library(Rrelperm)
library(pracma)
library(ggplot2)
library(dplyr)
library(magrittr)

p_pvt <- c(8000, 7500, 7280, 7250, 7000, 6500, 6000, 5500, 5000, 4500, 4000, 3500,  
           3000,    2500,   2000,   1500,   1000) # psia

Bo <- c(12.732, 13.044, 13.192, 1.054, 1.041,   1.018,  1.002,  0.983,  0.965, 0.947,
        0.93,   0.913,  0.896,  0.877,  0.858,  0.839,  0.819) # RB/STB

Rs <- c(22527,  22527,  22527,  860,    819,    754,    704,    648,    593,    541, 490, 440,
        389,    336,    280,    228,    169) #SCF/STB

Bg <- c(0.565,  0.579,  0.586,  0.587,  0.595,  0.613,  0.634,  0.661, 0.694, 0.737,
        0.795,  0.817,  0.997,  1.178,  1.466,  1.963, 2.912) / 1000 # RB/SCF

Rv <- c(44.4,   44.4,   44.4,   44.3,   43.9,   40.3,   36.5,   32.9,   29.2,   25.4,   21.4,   17.6, 
        13.7, 10.5, 7.9, 5.8,   4.4) / 1e6 # STB/SCF

cw <- 2e-6

Bwi <- 1.0

Bw <- Bwi * exp(cw * (p_pvt[1] - p_pvt))


muo <- c(0.049, 0.047,  0.046,  19.541, 20.965, 23.958, 26.338, 29.633, 33.319, 
         37.401, 42.161,    47.465, 53.765, 61.887, 72.143, 83.478, 99.049)  # cp

muw <- rep(0.25, length(p_pvt))

mug <- c(0.049, 0.047,  0.046,  0.046,  0.0449, 0.042,  0.0393, 0.0366, 0.0339, 
         0.0312,    0.0283, 0.0254, 0.0225, 0.0198, 0.0174, 0.0154, 0.014) # cp

pvt_table <- data.frame(p = p_pvt, Bo = Bo, Rs = Rs, Rv = Rv, Bg = Bg, 
                               Bw = Bw, muo = muo, mug = mug, muw = muw)

p <- c(8000,    7500,   7280,   7250,   7000,   6500,   6000,   5500,   5000,   4500,   4000,   3500,   
           3000,    2500,   2000,   1500,   1000) # psia

Gi <- rep(0, length.out = length(p))

wf <- rep(1, length.out = length(p))

# generating a set of pseudo relative permeability curves using 
# laboratory 'Kr' values

sg_lab <- c(0.05, 0.152, 0.248, 0.352, 0.448, 0.552, 0.65)

krg_lab <- c(0, 0.05, 0.09, 0.18, 0.3, 0.5, 1)

kro_lab <- c(1, 0.6, 0.35, 0.13, 0.04, 0.01, 0)

swcrit_lab <- 0.2

sgcrit_lab <- 0.05

sorgr_lab <- 0.15

fun_kr <- function(x, swcrit, sgcrit, sorg, sg, krg, kro) {

  kr_est <- Rrelperm::kr2p_gl(SWCON = swcrit, SOIRG = sorg, SORG = sorg, 
                               SGCON = sgcrit, SGCRIT = sgcrit, KRGCL = 1, 
                               KROGCG = 1, NG = x[1], NOG = x[2], NP = 101)

   krg_est_sub <- approx(x = kr_est[,1], y = kr_est[,3], xout = sg, rule = 2)$y

   kro_est_sub <- approx(x = kr_est[,1], y = kr_est[,4], xout = sg, rule = 2)$y

   error <- (krg - krg_est_sub) ^ 2 + (kro - kro_est_sub) ^ 2

   return(error)
}

par <- c(2, 2)

opt_results <- minpack.lm::nls.lm(par = par, fn = fun_kr, swcrit = swcrit_lab, 
                                  sgcrit = sgcrit_lab, sorg = sorgr_lab, sg = sg_lab, 
                                  krg = krg_lab, kro = kro_lab, 
                                  lower = c(0.1,0.1), upper = c(10,10))

sol <- opt_results$par

sol

rel_perm <- as.data.frame(Rrelperm::kr2p_gl(SWCON = swcrit_lab, SOIRG = sorgr_lab, 
                                            SORG = sorgr_lab, SGCON = sgcrit_lab, 
                                            SGCRIT = sgcrit_lab, KRGCL = 1, 
                                            KROGCG = 1, NG = sol[1], NOG = sol[2], 
                                            NP = 101))

colnames(rel_perm) <- c("Sg", "Sl", "Krg", "Krog")

p_forecast <- p

Gi_forecast <- Gi

wf_forecast <- wf

time_lst_forecast <- mbal_time(c(1:length(p_forecast)), "year")

forecast_lst <- mbal_forecast_param_gas(input_unit = "Field", output_unit = "Field",
                                        G = 1e6, phi = 0.1, swi = 0.2,
                                        pd = 7255, p = p_forecast, pvt = pvt_table,
                                        M = 0, cf = 2e-6, wf = wf_forecast,
                                        rel_perm = rel_perm)

forecast_results <- mbal_forecast_gas(forecast_lst, time_lst_forecast)

forecast_results


reservoir_performance_table <- data.frame(p = p)

reservoir_performance_table$`RF_oil` <- c(0,    2.4,    3.5,    4.2,    5.8, 7.6,   10.4,
                                                13.8,   16.6,   20, 23.1,   25.9,   28.5,   30.8,
                                                32.7,   34.1,   35.2) / 100

reservoir_performance_table$`Sg` <- c(100,  100,    100,    99.99,  99.91,  99.29,  98.68,  
                                      98.15,    97.66,  97.21,  96.79,  96.45,  96.15,  
                                      95.99,    95.9,   95.9, 95.95) * 0.8 / 100


reservoir_performance_table$`GOR` <- c(22.53,   22.53, 22.53, 22.55, 22.78, 24.81,
                                       27.4, 30.4, 34.25,   39.37, 46.73, 56.82, 72.99, 
                                       95.24, 126.58, 172.41, 227.27) * 1000 # SCF/STB


p1 <- forecast_results %>% ggplot(aes(`P (psia)`, SGg, color = "Forecast")) +
  geom_point(size = 3) +
  geom_point(data = reservoir_performance_table, aes(`p`, Sg, color = "Field"))+
  scale_color_manual(name="Data", values=c("Forecast" = "red", "Field" = "black")) +
  ggtitle("Gas Saturation Plot") +
  theme_bw()

p1


p2 <- forecast_results %>% ggplot(aes(`P (psia)`, `GOR (SCF/STB)`, color = "Forecast")) +
  geom_point(size = 3) +
  geom_point(data = reservoir_performance_table, aes(`p`, GOR, color = "Field"))+
  scale_color_manual(name="Data", values=c("Forecast" = "red", "Field" = "black")) +
  ggtitle("GOR Plot") +
  theme_bw()

p2


p3 <- forecast_results %>% ggplot(aes(`P (psia)`, `RF_oil`, color = "Forecast")) +
  geom_point(size = 3) +
  geom_point(data = reservoir_performance_table, aes(`p`, RF_oil, color = "Field"))+
  scale_color_manual(name="Data", values=c("Forecast" = "red", "Field" = "black")) +
  ggtitle("Oil Recovery Plot") +
  theme_bw()

p3


p4 <- forecast_results %>% ggplot(aes(`P (psia)`, `Liq_volume`, color = "Forecast")) +
  geom_point(size = 3) +
  scale_color_manual(name="Data", values=c("Forecast" = "red")) +
  ggtitle("CCE Liquid Volume Plot") +
  theme_bw()

p4


p5 <- forecast_results  %>%
  tidyr::pivot_longer(cols = c("Igd", "Ifwd"), names_to = "Drive Mechanism",
                      values_to = "Fraction", values_drop_na = TRUE) %>%
  ggplot(aes(`P (psia)`, Fraction, fill = `Drive Mechanism`)) +
  geom_area() +
  ggtitle("Energy Plot") +
  theme_bw()

p5

Example 3: Rich Gas Condensate Reservoir [@Walsh1994a]

Part I: Reservoir Forecast

library(Rmbal)
library(Rrelperm)
library(pracma)
library(ggplot2)
library(dplyr)
library(magrittr)

p_pvt <- c(5800,    5550,   5450,   5420,   5300,   4800,   4300,   3800,   3300,   2800,   2300,   
           1800,    1300,   800) # psia

Bo <- c(4.382,  4.441,  4.468,  2.378,  2.366,  2.032,  1.828,  1.674,  1.554,  
        1.448,  1.36,   1.279,  1.2,    1.131) # RB/STB

Rs <- c(6042,   6042,   6042,   2795,   2750,   2128,   1730,   1422,   1177,   960,    776,    607,    
        443,    293) #SCF/STB

Bg <- c(0.725,  0.735,  0.739,  0.74,   0.743,  0.758,  0.794,  0.854,  0.947,  
        1.09,   1.313,  1.677,  2.316,  3.695) / 1000 # RB/SCF

Rv <- c(165.5,  165.5,  165.5,  164.2,  156.6,  114,    89, 65.2,   48.3,   35, 25, 
        19, 15, 13.5) / 1e6 # STB/SCF

cw <- 2e-6

Bwi <- 1.0

Bw <- Bwi * exp(cw * (p_pvt[1] - p_pvt))


muo <- c(0.0612,    0.062,  0.0587, 0.135,  0.1338, 0.1826, 0.2354, 0.3001, 0.3764, 
         0.4781,    0.6041, 0.7746, 1.0295, 1.358)  # cp

muw <- rep(0.25, length(p_pvt))

mug <- c(0.0612,    0.062,  0.0587, 0.0581, 0.0554, 0.0436, 0.0368, 0.0308, 0.0261, 
         0.0222,    0.0191, 0.0166, 0.0148, 0.0135) # cp

pvt_table <- data.frame(p = p_pvt, Bo = Bo, Rs = Rs, Rv = Rv, Bg = Bg, 
                               Bw = Bw, muo = muo, mug = mug, muw = muw)

p <- c(5800,    5550,   5450,   5420,   5300,   4800,   4300,   3800,   3300,   2800,   2300,   
           1800,    1300,   800) # psia

Gi <- rep(0, length.out = length(p))

wf <- rep(1, length.out = length(p))

# generating a set of pseudo relative permeability curves using 
# laboratory 'Kr' values

sg_lab <- c(0.05, 0.152, 0.248, 0.352, 0.448, 0.552, 0.65)

krg_lab <- c(0, 0.05, 0.09, 0.18, 0.3, 0.5, 1)

kro_lab <- c(1, 0.6, 0.35, 0.13, 0.04, 0.01, 0)

swcrit_lab <- 0.2

sgcrit_lab <- 0.05

sorgr_lab <- 0.15

fun_kr <- function(x, swcrit, sgcrit, sorg, sg, krg, kro) {

  kr_est <- Rrelperm::kr2p_gl(SWCON = swcrit, SOIRG = sorg, SORG = sorg, 
                               SGCON = sgcrit, SGCRIT = sgcrit, KRGCL = 1, 
                               KROGCG = 1, NG = x[1], NOG = x[2], NP = 101)

   krg_est_sub <- approx(x = kr_est[,1], y = kr_est[,3], xout = sg, rule = 2)$y

   kro_est_sub <- approx(x = kr_est[,1], y = kr_est[,4], xout = sg, rule = 2)$y

   error <- (krg - krg_est_sub) ^ 2 + (kro - kro_est_sub) ^ 2

   return(error)
}

par <- c(2, 2)

opt_results <- minpack.lm::nls.lm(par = par, fn = fun_kr, swcrit = swcrit_lab, 
                                  sgcrit = sgcrit_lab, sorg = sorgr_lab, sg = sg_lab, 
                                  krg = krg_lab, kro = kro_lab, 
                                  lower = c(0.1,0.1), upper = c(10,10))

sol <- opt_results$par

sol

rel_perm <- as.data.frame(Rrelperm::kr2p_gl(SWCON = swcrit_lab, SOIRG = sorgr_lab, 
                                            SORG = sorgr_lab, SGCON = sgcrit_lab, 
                                            SGCRIT = sgcrit_lab, KRGCL = 1, 
                                            KROGCG = 1, NG = sol[1], NOG = sol[2], 
                                            NP = 101))

colnames(rel_perm) <- c("Sg", "Sl", "Krg", "Krog")

p_forecast <- p

Gi_forecast <- Gi

wf_forecast <- wf

time_lst_forecast <- mbal_time(c(1:length(p_forecast)), "year")

forecast_lst <- mbal_forecast_param_gas(input_unit = "Field", output_unit = "Field",
                                        G = 1e6, phi = 0.1, swi = 0.2,
                                        pd = 5430, p = p_forecast, pvt = pvt_table,
                                        M = 0, cf = 2e-6, wf = wf_forecast,
                                        rel_perm = rel_perm)

forecast_results <- mbal_forecast_gas(forecast_lst, time_lst_forecast)

forecast_results

reservoir_performance_table <- data.frame(p = p)

reservoir_performance_table$`RF_oil` <- c(0,    1.3,    1.9,    2.1,    2.6,    7,  10.1,   13.3,
                                          16.2, 18.4,   20.2,   21.6,   22.8,   23.7) / 100

reservoir_performance_table$`Sg` <- c(100,  100,    100,    99.2,   95, 81.8,   78.7,   76.7,   76.3,
                                      76.6, 77.2,   78.3,   79.5,   80.6) * 0.8 / 100


reservoir_performance_table$`GOR` <- c(6.04,    6.04,   6.04,   6.09,   6.39,   8.77,   11.22,  
                                       15.29,   20.62, 28.45, 39.82, 52.42, 66.48,  
                                       73.99) * 1000 # SCF/STB


p1 <- forecast_results %>% ggplot(aes(`P (psia)`, SGg, color = "Forecast")) +
  geom_point(size = 3) +
  geom_point(data = reservoir_performance_table, aes(`p`, Sg, color = "Field"))+
  scale_color_manual(name="Data", values=c("Forecast" = "red", "Field" = "black")) +
  ggtitle("Gas Saturation Plot") +
  theme_bw()

p1


p2 <- forecast_results %>% ggplot(aes(`P (psia)`, `GOR (SCF/STB)`, color = "Forecast")) +
  geom_point(size = 3) +
  geom_point(data = reservoir_performance_table, aes(`p`, GOR, color = "Field"))+
  scale_color_manual(name="Data", values=c("Forecast" = "red", "Field" = "black")) +
  ggtitle("GOR Plot") +
  theme_bw()

p2


p3 <- forecast_results %>% ggplot(aes(`P (psia)`, `RF_oil`, color = "Forecast")) +
  geom_point(size = 3) +
  geom_point(data = reservoir_performance_table, aes(`p`, RF_oil, color = "Field"))+
  scale_color_manual(name="Data", values=c("Forecast" = "red", "Field" = "black")) +
  ggtitle("Oil Recovery Plot") +
  theme_bw()

p3


p4 <- forecast_results %>% ggplot(aes(`P (psia)`, `Liq_volume`, color = "Forecast")) +
  geom_point(size = 3) +
  scale_color_manual(name="Data", values=c("Forecast" = "red")) +
  ggtitle("CCE Liquid Volume Plot") +
  theme_bw()

p4


p5 <- forecast_results  %>%
  tidyr::pivot_longer(cols = c("Igd", "Ifwd"), names_to = "Drive Mechanism",
                      values_to = "Fraction", values_drop_na = TRUE) %>%
  ggplot(aes(`P (psia)`, Fraction, fill = `Drive Mechanism`)) +
  geom_area() +
  ggtitle("Energy Plot") +
  theme_bw()

p5

References



Try the Rmbal package in your browser

Any scripts or data that you put into this service are public.

Rmbal documentation built on July 8, 2020, 7:16 p.m.