Ex.2 Exercise Notebook — Building a Surrogate Model with Volve Field Data —¶

Practical Reservoir Engineering — Data-Driven Simulator Series, Session 2 Exercise (Ex.2)

This notebook follows the structure of DataDrivenSimulator_S2_Exercise.pptx (Part F): data reload → input/output design → RSM/regression → random forest → (optional) NN → accuracy validation → time-series-split validation. It records the actual building and validation of surrogate models using the field-wide monthly production and pressure data from the Volve field data loaded in the Session 1 exercise (DataDrivenSimulator_S1_Exercise_notebook.ipynb).

Data used: Volve_Production_Data_Processed.csv (112 monthly rows, September 2007 – December 2016)

Step 1: Environment Setup Check / Connection from Session 1¶

In [1]:
import pandas as pd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import sklearn
import warnings

from sklearn.model_selection import train_test_split, KFold, cross_val_score
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from sklearn.pipeline import make_pipeline

plt.rcParams["font.family"] = "Yu Gothic"
plt.rcParams["axes.unicode_minus"] = False

print("pandas", pd.__version__, "/ numpy", np.__version__,
      "/ matplotlib", matplotlib.__version__, "/ scikit-learn", sklearn.__version__)

DATA_DIR = r"Research/References/Volve"  # local path to the Volve exercise CSVs

pandas 3.0.3 / numpy 2.5.0 / matplotlib 3.11.1 / scikit-learn 1.9.0

Step 2: Data Reload and Feature Design¶

In Session 1's Ex.1 we handled per-well bottomhole-pressure data (Volve_BHFP_data.csv) and a per-well feature table (7 wells). A 7-well aggregate table has too few rows (samples) to serve as training data for a surrogate model, so in Session 2's Ex.2 we switch the main material to the field-wide monthly production and pressure data (Volve_Production_Data_Processed.csv, 112 rows).

Preprocessing note (the same pitfall as Session 1): this file also uses DD/MM/YYYY date notation, so we must continue to specify dayfirst=True, as noticed in Session 1.

In [2]:
DATA_DIR = r"Research/References/Volve"  # local path to the Volve exercise CSVs
df["Date"] = pd.to_datetime(df["Date"], dayfirst=True)  # Session 1 insight: dayfirst=True is required
df["elapsed_days"] = (df["Date"] - df["Date"].min()).dt.days
print(df.shape)
print(df["Date"].min(), "to", df["Date"].max())
df.head()

(112, 9)
2007-09-01 00:00:00 to 2016-12-01 00:00:00
Out[2]:
Date p (psia) Np (STB) Gp (SCF) Wp (STB) Gi (SCF) Wi (STB) Rp (SCF/STB) elapsed_days
0 2007-09-01 4780.59 0 0 0.0 0 0 0.0 0
1 2007-10-01 4780.59 0 0 0.0 0 0 0.0 30
2 2007-11-01 4780.59 0 0 0.0 0 0 0.0 61
3 2007-12-01 4780.59 0 0 0.0 0 0 0.0 91
4 2008-01-01 4780.59 0 0 0.0 0 0 0.0 122

Step 3: Designing Input/Output Variables and Checking Multicollinearity¶

  • Output variable (y): p (psia) — reservoir pressure
  • Input variables (X): elapsed_days, Np (STB) / Gp (SCF) / Wp (STB) (cumulative production), Gi (SCF) / Wi (STB) (cumulative injection)
In [3]:
feature_cols = ["elapsed_days", "Np (STB)", "Gp (SCF)", "Wp (STB)", "Gi (SCF)", "Wi (STB)"]
X = df[feature_cols].values
y = df["p (psia)"].values

print("Rows where Gi (cumulative gas injection) is non-zero:", (df["Gi (SCF)"] != 0).sum(), "/", len(df))
corr = df[["Np (STB)", "Gp (SCF)", "Wp (STB)", "Wi (STB)"]].corr()
print(corr.round(4))

Rows where Gi (cumulative gas injection) is non-zero: 0 / 112
          Np (STB)  Gp (SCF)  Wp (STB)  Wi (STB)
Np (STB)    1.0000    0.9998    0.8814    0.9664
Gp (SCF)    0.9998    1.0000    0.8902    0.9711
Wp (STB)    0.8814    0.8902    1.0000    0.9731
Wi (STB)    0.9664    0.9711    0.9731    1.0000

Observation (multicollinearity): the correlation between Np (STB) and Gp (SCF) is 0.9998, an almost perfectly collinear relationship (because Volve is dominated by dissolved gas, once the oil volume is fixed the gas volume is almost uniquely determined). Gi (SCF) (cumulative gas injection) is always 0 across all 112 rows — because Volve is mainly water-flood (water injection) and does no gas injection — so this variable is meaningless as a surrogate-model input (it also appears with importance 0 in the feature importance of Step 5). We confirmed with real data the point from Part B6 that "centering/standardizing input variables mitigates multicollinearity."

Step 4: Building RSM (Polynomial Regression) and Linear Regression Models (Hold-out Validation)¶

In [4]:
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
print("Training data:", len(X_train), "rows / Test data:", len(X_test), "rows")

def evaluate(name, model, X_tr, y_tr, X_te, y_te):
    model.fit(X_tr, y_tr)
    pred = model.predict(X_te)
    rmse = mean_squared_error(y_te, pred) ** 0.5
    mae = mean_absolute_error(y_te, pred)
    r2 = r2_score(y_te, pred)
    print(f"{name:12s}  RMSE={rmse:7.2f} psia  MAE={mae:7.2f} psia  R2={r2:6.3f}")
    return dict(name=name, model=model, rmse=rmse, mae=mae, r2=r2, pred=pred)

results = {}
results["linear"] = evaluate("Linear", make_pipeline(StandardScaler(), LinearRegression()),
                               X_train, y_train, X_test, y_test)
results["poly2"] = evaluate("Poly(deg2)", make_pipeline(StandardScaler(), PolynomialFeatures(degree=2, include_bias=False), LinearRegression()),
                              X_train, y_train, X_test, y_test)

Training data: 89 rows / Test data: 23 rows
Linear        RMSE=  89.74 psia  MAE=  51.21 psia  R2= 0.831
Poly(deg2)    RMSE=  75.51 psia  MAE=  48.77 psia  R2= 0.881

Step 5: Building a Random Forest Model (Hold-out Validation)¶

In [5]:
rf = RandomForestRegressor(n_estimators=300, random_state=42, oob_score=True)
results["rf"] = evaluate("RF", rf, X_train, y_train, X_test, y_test)
print("OOB score:", round(rf.oob_score_, 3))

importances = pd.Series(rf.feature_importances_, index=feature_cols).sort_values(ascending=False)
print("\nFeature importance:")
print(importances.round(3))

RF            RMSE=  49.32 psia  MAE=  31.11 psia  R2= 0.949
OOB score: 0.91

Feature importance:
elapsed_days    0.231
Np (STB)        0.216
Gp (SCF)        0.191
Wi (STB)        0.182
Wp (STB)        0.180
Gi (SCF)        0.000
dtype: float64

Observation: Random forest is the most accurate (R2=0.949, RMSE≈49 psia). In the feature importance, elapsed_days, Np, Gp, Wi, and Wp are close together, while Gi (SCF) has importance 0.000 — the real-data structure of a "variable that is zero throughout," confirmed in Step 3, is reflected directly in the model's importance.

Step 6: Building a Neural Network Model (Optional Task, Hold-out Validation)¶

In [6]:
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    mlp = make_pipeline(StandardScaler(), MLPRegressor(hidden_layer_sizes=(32, 16), max_iter=5000, random_state=42))
    results["mlp"] = evaluate("MLP(NN)", mlp, X_train, y_train, X_test, y_test)

MLP(NN)       RMSE= 137.20 psia  MAE= 109.06 psia  R2= 0.606

Observation: The neural network (MLP) gave RMSE≈137 psia, R2=0.606, below RSM (R2=0.881) and random forest (R2=0.949). With only 89 training rows, the result bears out with real data the general point from Part C10 that "neural networks require large data, and for small-to-medium data random forest is more robust."

Step 7: 5-fold Cross-Validation (Random Forest)¶

In [7]:
kf = KFold(n_splits=5, shuffle=True, random_state=42)
scores = cross_val_score(RandomForestRegressor(n_estimators=300, random_state=42), X, y, cv=kf, scoring="r2")
print("5-fold R2:", np.round(scores, 3))
print("Mean R2:", round(scores.mean(), 3), " / Std:", round(scores.std(), 3))

5-fold R2: [0.95  0.901 0.768 0.897 0.972]
Mean R2: 0.897  / Std: 0.071

Observation: With the hold-out method alone R2=0.949, but in 5-fold cross-validation R2 varied by fold from 0.768 to 0.972, with a mean of 0.897. We confirmed with real data that checking across multiple splits, rather than relying on a single split, increases the robustness of the accuracy evaluation (a concrete example of Part E2).

Step 8: Validation with a Time-Series Split (An Important Insight)¶

So far we validated with a random 80/20 split, but Volve's data is time-series data. As a validation more appropriate for practical prediction (forecasting future pressure), we also try a split that follows the chronological order: the first 80% for training and the remaining 20% (from February 2015 onward) for testing.

In [8]:
n = len(df)
cut = int(n * 0.8)
X_train_t, X_test_t = X[:cut], X[cut:]
y_train_t, y_test_t = y[:cut], y[cut:]
print("Training period:", df["Date"].iloc[0].date(), "to", df["Date"].iloc[cut-1].date())
print("Test period:", df["Date"].iloc[cut].date(), "to", df["Date"].iloc[-1].date())
print("Training pressure range: {:.2f} to {:.2f} psia".format(y_train_t.min(), y_train_t.max()))
print("Test pressure range: {:.2f} to {:.2f} psia".format(y_test_t.min(), y_test_t.max()))

rf_t = RandomForestRegressor(n_estimators=300, random_state=42)
poly_t = make_pipeline(StandardScaler(), PolynomialFeatures(degree=2, include_bias=False), LinearRegression())

time_results = {}
for name, model in [("RF", rf_t), ("Poly(deg2)", poly_t)]:
    model.fit(X_train_t, y_train_t)
    pred = model.predict(X_test_t)
    rmse = mean_squared_error(y_test_t, pred) ** 0.5
    r2 = r2_score(y_test_t, pred)
    time_results[name] = dict(pred=pred, rmse=rmse, r2=r2)
    print(f"[time-series split] {name:12s}  RMSE={rmse:7.2f} psia  R2={r2:7.3f}")

Training period: 2007-09-01 to 2015-01-01
Test period: 2015-02-01 to 2016-12-01
Training pressure range: 4085.44 to 5114.90 psia
Test pressure range: 4815.37 to 5273.87 psia

[time-series split] RF            RMSE= 134.26 psia  R2= -0.419
[time-series split] Poly(deg2)    RMSE= 408.79 psia  R2=-12.153

Observation (important): With a random split, RF R2=0.949 was good, but with a time-series split RF R2 worsened to -0.42 and polynomial regression to R2=-12.15. The cause is that the test period (February 2015 – December 2016) reaches a maximum pressure of 5273.87 psia, whereas the training period (September 2007 – January 2015) only experienced up to 5114.90 psia, forcing the model to predict in the region outside the training data (the extrapolation region). We quantitatively confirmed, with Volve's real data, what Parts E5 and B10 explained conceptually as "be careful about applying a model outside the training-data range." Evaluating accuracy with a random split alone risks overlooking this kind of extrapolation risk — the biggest lesson of this exercise.

Step 9: Visualizing the Results¶

In [9]:
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax = axes[0]
colors = {"linear": "#4A3AA7", "poly2": "#2A78D6", "rf": "#1BAF7A", "mlp": "#EB6834"}
labels = {"linear": "Linear", "poly2": "Poly(deg2)", "rf": "Random forest", "mlp": "NN(MLP)"}
for key, res in results.items():
    ax.scatter(y_test, res["pred"], s=28, alpha=0.75, color=colors[key], label=f"{labels[key]} (R2={res['r2']:.2f})")
lims = [min(y_test.min(), y_test.min()) - 50, y_test.max() + 50]
ax.plot(lims, lims, color="#89877E", linestyle="--", linewidth=1)
ax.set_xlabel("Actual p (psia)"); ax.set_ylabel("Predicted p (psia)")
ax.set_title("Random split: Predicted vs. Actual")
ax.legend(fontsize=8); ax.grid(alpha=0.3)

ax2 = axes[1]
names = list(results.keys())
rmses = [results[k]["rmse"] for k in names]
ax2.bar([labels[k] for k in names], rmses, color=[colors[k] for k in names])
ax2.set_ylabel("RMSE (psia)")
ax2.set_title("RMSE by Method (random split)")
for i, v in enumerate(rmses):
    ax2.text(i, v + 2, f"{v:.1f}", ha="center", fontsize=9)
plt.setp(ax2.get_xticklabels(), rotation=15, ha="right")

fig.tight_layout()
plt.show()
No description has been provided for this image
In [10]:
fig, ax = plt.subplots(figsize=(11, 5))
ax.plot(df["Date"], df["p (psia)"], color="#89877E", linewidth=1.2, label="Actual p (psia)")
ax.plot(df["Date"].iloc[cut:], time_results["RF"]["pred"], color="#1BAF7A", linewidth=1.6,
        linestyle="--", label="RF prediction (time-series split test interval)")
ax.axvline(df["Date"].iloc[cut], color="#EB6834", linestyle=":", linewidth=1.5)
ax.text(df["Date"].iloc[cut], ax.get_ylim()[1], " <- train/test boundary", color="#EB6834", fontsize=9, va="top")
ax.set_xlabel("Date"); ax.set_ylabel("Reservoir pressure p (psia)")
ax.set_title("Time-Series Split: Prediction Diverges by Extrapolation in the Test Period")
ax.legend(fontsize=9); ax.grid(alpha=0.3)
fig.tight_layout()
plt.show()
No description has been provided for this image

Summary / Handover to Session 3¶

  • Using Volve's monthly production and pressure data (112 rows), we built surrogate models of reservoir pressure p with three methods — RSM (linear and quadratic polynomial regression), random forest, and neural network (optional task) — and validated accuracy with the hold-out method and 5-fold cross-validation, running the full workflow along the flow of the pptx material (F1–F10).
  • Three points learned from the real-data validation, to be reflected in the Session 2 material and in cautions to participants:
    1. The per-well feature table (7 wells) has too few samples, so we used the field-wide monthly data of 112 rows as the main material for surrogate-model building.
    2. Np (STB) and Gp (SCF) are almost collinear (correlation 0.9998), and Gi (SCF) is zero throughout (Volve is mainly water-flood and does no gas injection) — usable as real examples of multicollinearity and an uninformative variable.
    3. With a random split RF R2=0.949 is good, but with a time-series split (test from February 2015 onward) RF R2 worsens to -0.42 and polynomial regression to R2=-12.15. The cause is extrapolation, where the test-period pressure exceeds the training-period range. The biggest lesson of this exercise: the conclusion changes completely depending on the validation method (random split vs. time-series split).
  • These insights about "the range in which a surrogate model can be trusted" and "results change depending on the validation design" connect directly to the discussion of uncertainty assessment and prediction range in Session 3 (data-driven history matching).