Geostatistical Shale Reservior Modeling

The dataset used in this tutorial consists of two objects. The first object, WELL_PARAMETERS is a data frame that contains all well data such as well locations, completion parameters, oil density etc. The second object, OilRates contains 188 oil production profiles (time series). The dataset was analyzed in great detail in chapter 4 of Grujic (2017).

library(fdagstat)
library(fda)
library(DGSA)
colnames(WELL_PARAMETERS)
FALSE  [1] "Geol3DSpacing"   "GeolAPIGrav"     "GeolX_Rel"      
FALSE  [4] "GeolY_Rel"       "ProdVert200"     "ProdVert300"    
FALSE  [7] "PetroPor"        "PetroSwt"        "PetroTOC"       
FALSE [10] "PetroVClay"      "PetroCarb"       "PetroVQtz"      
FALSE [13] "PetroPyr"        "StagesPumped"    "StimLength"     
FALSE [16] "No..Screenouts"  "TotalFluid.bbl"  "AmtSlickwater"  
FALSE [19] "AmtCrosslink"    "AmtAcid"         "AmtLinear"      
FALSE [22] "CleanFluidTotal" "TotalProppant"   "Amt100mesh"     
FALSE [25] "AmtResin"
par(mfrow=c(1,2))
plot(WELL_PARAMETERS$GeolX_Rel, WELL_PARAMETERS$GeolY_Rel, pch=19, col="red", xlab="X Location", ylab="Y Location")
matplot(OilRates, type="l", col="blue", xlab="Time (days)", ylab="Oil Rate (stb/day)")
Left: Well Locations. Right: Oil Production Profiles

Left: Well Locations. Right: Oil Production Profiles

We will start by performing exploratory data analysis by means of distance based generalized sensitivity analysis (DGSA, Fenwick et al (2015)). The method is appropriate in this application since it quantifies the influence of explanatory variables on the entire production profiles and not some discrete elements of production such as “peak oil” or cumulative production at 3, 6, or 9 months.

The DGSA method starts by clustering the oil responses. For this purpose we will compute Euclidean distances between the production profiles and then apply simple kmeans clustering.

D <- dist(t(OilRates), 'euclidean')
Clustering <- kmeans(D, 2) # two clusters

mds <- cmdscale(D, eig=TRUE)
prVar <- round(100*mds$eig/sum(mds$eig),2)
par(mfrow=c(1,2))
plot(mds$points, main="Low dimensional Distance plot", xlab=paste(prVar[1],"% Variance"), ylab=paste(prVar[2],"% Variance"))
points(mds$points[Clustering$cluster==1, ], col="red", pch=19)
points(mds$points[Clustering$cluster==2, ], col="blue", pch=19)
colours = rep("red", length(Clustering$cluster))
colours[Clustering$cluster==2] <- "blue"
matplot(OilRates, type="l", lty=1, col=colours, xlab="Time (days)", ylab="Oil Rate (stb/day)")
Left: A low dimensional MDS plot. Right: A display of clustering in original data

Left: A low dimensional MDS plot. Right: A display of clustering in original data

To perform EDA and DGSA we will use the DGSA package.

DGSA Scatter Plots

plotDGSAscatter(WELL_PARAMETERS[, -1], mds$points)
Colored MDS Plots

Colored MDS Plots

DGSA CDF Plots

plotCDFS(Clustering$cluster, WELL_PARAMETERS[, -1], .code="all*")
Cluster Specific CDFS of Each Parameter

Cluster Specific CDFS of Each Parameter

DGSA Results

sensitivity <- dgsa(Clustering$cluster, WELL_PARAMETERS[,-1], .interactions = FALSE, .progress = F)
plotParetoDGSA(sensitivity)
DGSA Pareto Plot

DGSA Pareto Plot

Correlation Table

kable(round(cor(WELL_PARAMETERS),2))
Geol3DSpacing GeolAPIGrav GeolX_Rel GeolY_Rel ProdVert200 ProdVert300 PetroPor PetroSwt PetroTOC PetroVClay PetroCarb PetroVQtz PetroPyr StagesPumped StimLength No..Screenouts TotalFluid.bbl AmtSlickwater AmtCrosslink AmtAcid AmtLinear CleanFluidTotal TotalProppant Amt100mesh AmtResin
Geol3DSpacing 1.00 0.06 -0.02 0.02 -0.04 0.02 0.00 0.03 -0.08 0.16 -0.10 0.12 0.02 0.01 0.00 -0.14 -0.07 -0.16 0.19 -0.13 0.02 -0.11 0.11 -0.10 0.26
GeolAPIGrav 0.06 1.00 -0.22 0.84 0.44 0.46 0.54 -0.03 0.28 0.08 -0.01 0.06 0.11 -0.41 -0.43 0.01 -0.40 -0.36 -0.16 -0.06 0.03 -0.32 -0.30 0.05 0.09
GeolX_Rel -0.02 -0.22 1.00 0.18 0.11 0.18 -0.08 0.09 -0.09 -0.01 -0.01 0.00 -0.02 0.02 0.22 -0.05 0.13 0.06 0.15 -0.16 0.07 0.19 0.09 -0.07 0.03
GeolY_Rel 0.02 0.84 0.18 1.00 0.51 0.53 0.54 0.06 0.22 0.09 -0.04 0.08 0.09 -0.38 -0.29 0.01 -0.33 -0.35 -0.01 -0.21 0.04 -0.18 -0.21 0.02 0.19
ProdVert200 -0.04 0.44 0.11 0.51 1.00 0.85 0.18 0.10 0.16 0.05 -0.03 0.05 0.11 -0.28 -0.14 0.01 -0.20 -0.16 -0.11 -0.04 0.04 -0.06 -0.17 -0.09 0.13
ProdVert300 0.02 0.46 0.18 0.53 0.85 1.00 0.27 0.00 0.06 -0.08 0.11 -0.07 0.08 -0.33 -0.15 -0.01 -0.21 -0.20 -0.07 -0.09 0.00 -0.04 -0.16 -0.12 0.15
PetroPor 0.00 0.54 -0.08 0.54 0.18 0.27 1.00 -0.62 0.02 -0.48 0.53 -0.47 -0.05 -0.27 -0.25 -0.01 -0.24 -0.26 0.02 -0.09 -0.04 -0.20 -0.06 0.05 0.25
PetroSwt 0.03 -0.03 0.09 0.06 0.10 0.00 -0.62 1.00 0.03 0.83 -0.86 0.81 0.01 0.04 0.19 -0.14 0.12 0.10 0.03 -0.03 0.09 0.19 0.04 -0.04 -0.06
PetroTOC -0.08 0.28 -0.09 0.22 0.16 0.06 0.02 0.03 1.00 0.16 -0.15 0.18 0.54 -0.01 -0.11 0.09 -0.09 -0.10 -0.04 0.03 0.23 -0.15 -0.09 0.08 -0.02
PetroVClay 0.16 0.08 -0.01 0.09 0.05 -0.08 -0.48 0.83 0.16 1.00 -0.95 0.95 0.17 0.04 0.05 -0.13 0.04 0.06 -0.08 0.01 0.18 0.06 -0.03 -0.07 -0.11
PetroCarb -0.10 -0.01 -0.01 -0.04 -0.03 0.11 0.53 -0.86 -0.15 -0.95 1.00 -0.95 -0.15 -0.06 -0.12 0.15 -0.07 -0.08 0.04 -0.01 -0.16 -0.08 -0.01 0.03 0.10
PetroVQtz 0.12 0.06 0.00 0.08 0.05 -0.07 -0.47 0.81 0.18 0.95 -0.95 1.00 0.20 0.02 0.05 -0.15 0.02 0.04 -0.07 0.01 0.17 0.04 -0.03 0.00 -0.09
PetroPyr 0.02 0.11 -0.02 0.09 0.11 0.08 -0.05 0.01 0.54 0.17 -0.15 0.20 1.00 0.02 -0.12 -0.03 -0.11 -0.10 -0.07 -0.03 0.09 -0.12 -0.14 -0.03 -0.03
StagesPumped 0.01 -0.41 0.02 -0.38 -0.28 -0.33 -0.27 0.04 -0.01 0.04 -0.06 0.02 0.02 1.00 0.63 0.06 0.65 0.51 0.36 0.11 0.27 0.30 0.52 0.06 -0.23
StimLength 0.00 -0.43 0.22 -0.29 -0.14 -0.15 -0.25 0.19 -0.11 0.05 -0.12 0.05 -0.12 0.63 1.00 -0.10 0.79 0.66 0.39 0.29 0.10 0.47 0.69 -0.04 -0.05
No..Screenouts -0.14 0.01 -0.05 0.01 0.01 -0.01 -0.01 -0.14 0.09 -0.13 0.15 -0.15 -0.03 0.06 -0.10 1.00 0.04 -0.02 0.09 0.10 0.20 -0.02 -0.01 -0.03 -0.09
TotalFluid.bbl -0.07 -0.40 0.13 -0.33 -0.20 -0.21 -0.24 0.12 -0.09 0.04 -0.07 0.02 -0.11 0.65 0.79 0.04 1.00 0.89 0.36 0.39 0.17 0.53 0.82 0.04 -0.26
AmtSlickwater -0.16 -0.36 0.06 -0.35 -0.16 -0.20 -0.26 0.10 -0.10 0.06 -0.08 0.04 -0.10 0.51 0.66 -0.02 0.89 1.00 -0.10 0.50 -0.02 0.55 0.55 0.06 -0.45
AmtCrosslink 0.19 -0.16 0.15 -0.01 -0.11 -0.07 0.02 0.03 -0.04 -0.08 0.04 -0.07 -0.07 0.36 0.39 0.09 0.36 -0.10 1.00 -0.17 0.24 0.05 0.66 -0.03 0.38
AmtAcid -0.13 -0.06 -0.16 -0.21 -0.04 -0.09 -0.09 -0.03 0.03 0.01 -0.01 0.01 -0.03 0.11 0.29 0.10 0.39 0.50 -0.17 1.00 0.00 0.10 0.24 -0.02 -0.22
AmtLinear 0.02 0.03 0.07 0.04 0.04 0.00 -0.04 0.09 0.23 0.18 -0.16 0.17 0.09 0.27 0.10 0.20 0.17 -0.02 0.24 0.00 1.00 -0.01 0.18 0.00 -0.07
CleanFluidTotal -0.11 -0.32 0.19 -0.18 -0.06 -0.04 -0.20 0.19 -0.15 0.06 -0.08 0.04 -0.12 0.30 0.47 -0.02 0.53 0.55 0.05 0.10 -0.01 1.00 0.31 0.01 -0.20
TotalProppant 0.11 -0.30 0.09 -0.21 -0.17 -0.16 -0.06 0.04 -0.09 -0.03 -0.01 -0.03 -0.14 0.52 0.69 -0.01 0.82 0.55 0.66 0.24 0.18 0.31 1.00 0.00 0.20
Amt100mesh -0.10 0.05 -0.07 0.02 -0.09 -0.12 0.05 -0.04 0.08 -0.07 0.03 0.00 -0.03 0.06 -0.04 -0.03 0.04 0.06 -0.03 -0.02 0.00 0.01 0.00 1.00 -0.08
AmtResin 0.26 0.09 0.03 0.19 0.13 0.15 0.25 -0.06 -0.02 -0.11 0.10 -0.09 -0.03 -0.23 -0.05 -0.09 -0.26 -0.45 0.38 -0.22 -0.07 -0.20 0.20 -0.08 1.00

In DGSA, parameters that show a trend in the low dimensional scatter plots (DGSA Scatter) are deemed influential on the analyzed functional responses. To quantify parameter importance we compute cluster specific CDFS (DGSA CDFS) and then we compute the L1 norm of each CDF from its prior CDF. The prior CDF’s are computed on the entire parameter data (both clusters). A parameter that shows a trend in the low dimensional scatter plot will also have cluster specific CDF’s that are significantly different from the prior CDF (Fenwick et al 2015). The parameters are ranked based on their average departure (L1 norm) from their prior CDF’s (DGSA Results).

DGSA also allows for quantification of parameter interactions. Parameter interactions were not computed in this exercise due to the fact that the input parameters were highly correlated. DGSA was designed for quantification of parameter importances in computer experiments where input parameters are uncorrelated by design. Correlation table shown in ‘Correlation Table’ tab clearly shows that many parameters are correlated, this would would have made DGSA parameter interaction analysis misleading, hence it was not performed in this exercise.

Based on the sensitivity analysis we can see that hydraulic fracturing parameters are the most influential on shale oil response. Furthermore, wells “Y location” also appears to be highly influential on oil rates. What is surprising is the low ranking of “Petro TOC”. This is because the “Petro TOC” data is most likely corrupted given that it takes one of only three values (see DGSA CDF tab). The same is most likely true for the “Petro Pyr” parameter.

Fitting a Universal Trace Kriging Model

The first step in fitting a forecasting model is parameter selection. In shale reservoir forecasting there are three guiding principles that can be followed. The first prinicple is to use the parameters that are influential on the forecasted response. The DGSA analysis performed previously clearly identified hydraulic fracturing parameters as the most influential and they will be carried forward into our forecasting models. The second principle considers only those parameters that are available at the time of the intended use of the fitted model. Clearly, API gravity and petrophysical parameters are not availabe to reservoir modelers when forecasting undrilled wells. Therefore, such parameters will not be carried forward into the forecasting analysis. The third principle is concerned with the type of the considered statistical model. In our case we aim to fit statistical models that take spatial correlations between the wells into consideration. For this reason, both X and Y well locations will be carried forward even though X location was not found to be influential on oil responses in our DGSA analysis.

The final set of selected forecasting parameters (predictors) is given below:

##  [1] "GeolX_Rel"      "GeolY_Rel"      "ProdVert200"    "ProdVert300"   
##  [5] "StagesPumped"   "TotalProppant"  "StimLength"     "TotalFluid.bbl"
##  [9] "AmtSlickwater"  "AmtCrosslink"   "AmtAcid"        "AmtLinear"     
## [13] "AmtResin"

Next we will fit a modified universal trace kriging forecasting model (UTrK, Menafoglio et al (2013)). The model is of the following form:

\[ \mathcal{X}_i(t) = m_i(t) + r_i(t)\]

Where:

  • The drift is modeled with a functional linear model \(m_i(t)=f(\boldsymbol{z}_i)=\sum_{l=0}^Lf_l(\boldsymbol{z})a_l(t)\). Here, vector \(\boldsymbol{z} \in R^n\) contains all predictors, spatial and non-spatial.
  • The residuals are assumed to be second order stationary and spatially correlated. In other words, \(cov(r_i(t), r_j(t))=\int_T r_i(t)r_j(t)dt=c(\|\boldsymbol{s}_i - \boldsymbol{s}_j\|)\). Here, vector \(\boldsymbol{s} \in R^2\) specifies well locations.

For theoretical details of this modeling paradigm please consult Menafoglio et al (2013) or Grujic (2017).

Before we proceed to fit a model we will rescale all of the selected well parameters. We choose to bring all parameters to [0,1] scale which is simply accomplished by subtracting the minimum and dividing by the scaled maximum. A special attention needs to be given to rescaling of spatial parameters. In order to maintain the proportion of the scaled and actual distances we will subtract the minimum of each coordinate but we will divide both coordinates with max(max(x), max(y)). The following code does the rescaling and plots the data:

COVARIATES <- WELL_PARAMETERS[,c(3,4,5,6,14,23,15,17,18,19,20,21,25)]
XY <- COVARIATES[,1:2]
XY[,1] <- XY[,1] - min(XY[,1])
XY[,2] <- XY[,2] - min(XY[,2])
maxXY <- max(XY)
XY <- XY / max(XY)
COVARIATES <- as.data.frame(apply(COVARIATES, 2, function(x) {
                                                x <- x-min(x)
                                                x <- x/max(x)
                                                }))
COVARIATES[,1:2] <- XY
Well Locations: Old and New Coordinates

Well Locations: Old and New Coordinates

The aforementioned forecasting model can be fitted with the fdagstat package. We start in the usual way by initializing the fstat data structure:

set.seed(1)
train <- sample(nrow(WELL_PARAMETERS), floor(0.7*nrow(WELL_PARAMETERS)), replace=FALSE)
g <- fstat(NULL, 'Oil', COVARIATES[train,], OilRates[,train])
g <- estimateDrift("~.", g, Intercept = TRUE)
g <- fvariogram("~GeolX_Rel + GeolY_Rel", g, Nlags = 100, LagMax = 1, ArgStep = 1, comments=FALSE)

There are two things that need to be noted here. The formula used in drift estimation is using all covariates ("~."), while the formula used in variogram estimation uses only the X and Y locations ("~GeolX_Rel + GeolY_Rel").

Variogram Analysis

Next we will fit a trace variogram and infer a GLS model for the drift term.

Raw Variogram (OLS)

Initial Trace Variogram

Initial Trace Variogram 

g <- fitVariograms(g, vgm(1, "Mat", 0.17, 250000, kappa=0.4),
                         fitRanges = FALSE, forceNugget = FALSE)
plotVariogram(g)

The fitted omni directional variogram has a range of 0.17 in rescaled spatial units. This value equals 15055.06 meters in the original scale.

Variogram 1 Iteration (GLS)

Initial Trace Variogram (GLS-1)

Initial Trace Variogram (GLS-1) 

# Rescaling:
g$model$omni$Oil$psill <- g$model$omni$Oil$psill/ sum(g$model$omni$Oil$psill)
g <- addCovariance(g)
g <- estimateDrift("~.", g, .type = "GLS", Intercept = TRUE)
g <- fvariogram("~GeolX_Rel + GeolY_Rel", g, Nlags=100, LagMax = 1, ArgStep = 1,
                      useResidual = TRUE, comments=FALSE)
g <- fitVariograms(g, vgm(1, "Mat", 0.17, 250000, kappa=0.4),
                         fitRanges = FALSE, forceNugget = FALSE)
plotVariogram(g)

Here, we fit a GLS model for the drift. Note that the fitted trace variogram had to be rescaled to produce a covariance matrix whose scale was similar to the scale of the COVARIATES matrix. This was necesarry to avoid numerical issues with the GLS computaions.

Variogram 2 Iterations (GLS)

Initial Trace Variogram (GLS-2)

Initial Trace Variogram (GLS-2) 

# Rescaling:
g$model$omni$Oil$psill <- g$model$omni$Oil$psill/ sum(g$model$omni$Oil$psill)
g <- addCovariance(g)
g <- estimateDrift("~.", g, .type = "GLS",  Intercept = TRUE)
g <- fvariogram("~GeolX_Rel + GeolY_Rel", g, Nlags=100, LagMax = 1, ArgStep = 1,
                      useResidual = TRUE, comments=FALSE)
g <- fitVariograms(g, vgm(1, "Mat", 0.17, 250000, kappa=0.4),
                         fitRanges = FALSE, forceNugget = FALSE)
plotVariogram(g)

Forecasting New Functions

All

True Functions vs. Forecasts

True Functions vs. Forecasts

Individual 1

Individual 2

Individual 3

Individual 4

Fitting a Universal Co-Kriging of FPC Scores Model

Similarly to the first vignette of the fdagstat package (link)we can formulate an alternative FPC-based model. The FPC modeling modeling paradigm is analogous to the first vignette and the aforementioned modified universal trace kriging.

We first start with functional principal component decomposition of the OilRates. Since in shale reservoir forecasting data interpretation is very important we will work with VARIMAX rotated fpc’s. In the figure given below we are plotting both the raw (top row) and VARIMAX rotated (bottom row) principal components. Notice that the first rotated fpc describes late time oil production while the second rotated fpc describes the variation in early oil production. For more details about VARIMAX rotated fpcs please consult Ramsay and Silverman (2005).

library(fda)
basis <- create.bspline.basis(c(1, 600), 30, 4)
expandedFunctions <- smooth.basis(1:600, as.matrix(OilRates[,train]), basis)
expandedFunctions.pca <- pca.fd(expandedFunctions$fd, 2)
expandedFunctions.vrmx <- varmx.pca.fd(expandedFunctions.pca)
MEAN <- eval.fd(1:600, expandedFunctions.pca$meanfd)
FPCS <- eval.fd(1:600, expandedFunctions.pca$harmonics)
FPCS.vrmx <- eval.fd(1:600, expandedFunctions.vrmx$harmonics)
Regular and VARIMAX rotated FPCS

Regular and VARIMAX rotated FPCS

As in the first vignette we will standardize the rotated functional principal component scores.

sd1 <- sd(expandedFunctions.vrmx$scores[,1])
mn1 <- mean( expandedFunctions.vrmx$scores[,1])
sd2 <- sd(expandedFunctions.vrmx$scores[,2])
mn2 <- mean( expandedFunctions.vrmx$scores[,2])

Oil_pc1 <- (expandedFunctions.vrmx$scores[,1] - mn1)/sd1
Oil_pc2 <- (expandedFunctions.vrmx$scores[,2] - mn2)/sd2

Variogram Modeling

Variogram (OLS)

g <- fstat(NULL, 'Oil-rFPC1', COVARIATES[train, ], as.data.frame(Oil_pc1), scalar=TRUE)
g <- fstat(g,    'Oil-rFPC2', COVARIATES[train, ], as.data.frame(Oil_pc2), scalar=TRUE)
g <- estimateDrift("~.", g, 'OLS')
g <- fvariogram("~GeolX_Rel + GeolY_Rel", g, Nlags = 100, LagMax = 1, useResidual =TRUE, comment = FALSE)
g <- fitVariograms(g, vgm(1, "Mat", 0.12, 0, kappa=0.4), fitRanges = FALSE, forceNugget = FALSE)
plotVariogram(g)

Variogram (GLS-1)

g <- addCovariance(g)
g <- estimateDrift("~.", g, .type = "GLS")
g <- fvariogram("~GeolX_Rel + GeolY_Rel", g, Nlags=100,
                     LagMax = 1, ArgStep = 1,
                     useResidual = TRUE, comment = FALSE)
g <- fitVariograms(g, vgm(1, "Mat", 0.12, 0, kappa=0.4), fitRanges = FALSE, forceNugget = FALSE)
plotVariogram(g)

Variogram (GLS-2)

g <- addCovariance(g)
g <- estimateDrift("~.", g, .type = "GLS")
g <- fvariogram("~GeolX_Rel + GeolY_Rel", g, Nlags=100,
                     LagMax = 1, ArgStep = 1,
                     useResidual = TRUE, comment = FALSE)
g <- fitVariograms(g, vgm(1, "Mat", 0.12, 0, kappa=0.4), fitRanges = FALSE, forceNugget = FALSE)
plotVariogram(g)

What is convenient in the case of universal co-kriging of fpc scores is that we can visualize the spatial component of the drift models. The following code evaluates and visualizes the spatial components of the drifts on a fine spatial grid

library(gridExtra)
library(ggplot2)

grid <- expand.grid(X=seq(0, 1, length.out = 30), Y=seq(0, 1, length.out = 30))
grid$rFpc1 <- as.matrix(grid) %*% as.matrix(g$drift$`Oil-rFPC1`$Betas[2:3])
grid$rFpc2 <- as.matrix(grid[,1:2]) %*% as.matrix(g$drift$`Oil-rFPC2`$Betas[2:3])
p1 <- ggplot(grid) + geom_point(aes(x=X, y=Y, colour=rFpc1), size=4) + scale_colour_gradientn(colours=rainbow(9))
p2 <- ggplot(grid) + geom_point(aes(x=X, y=Y, colour=rFpc2), size=4) + scale_colour_gradientn(colours=rainbow(9))

grid.arrange(p1, p2, ncol=2)

We can see that the first rotated fpc, that describes late time oil production, has a NW-SE spatial trend while the second fpc, that describes early production, has a NE-SW spatial trend. Geological interpretation of these trends is given in Grujic (2017).

Forecasting

All Data

forecast.rFPC1 <- predictFstat(g, .newCoordinates = COVARIATES[-train,], .what="Oil-rFPC1", .type="UcoK", algIndependent = TRUE)
forecast.rFPC2 <- predictFstat(g, .newCoordinates = COVARIATES[-train,], .what="Oil-rFPC2", .type="UcoK", algIndependent = TRUE)

# Bring back to the original scale:
forecast.rFPC1$Forecast <- (forecast.rFPC1$Forecast * sd1) + mn1
forecast.rFPC2$Forecast <- (forecast.rFPC2$Forecast * sd2) + mn2

# Construct functions:
UcoKforecast <- NULL
for(i in 1:length(forecast.rFPC1$Forecast)){
  UcoKforecast <- rbind(UcoKforecast, t(MEAN + FPCS.vrmx %*% c(forecast.rFPC1$Forecast[i], forecast.rFPC2$Forecast[i])))
}
matplot(OilRates[,-train], type="l", col="blue", xlab="Time (days)", ylab="Oil Rate (bbl/day)")
matplot(t(UcoKforecast), type = "l", col="red", add=TRUE)
legend('topright', c("True", "Forecasted"), lty=1, col=c("red", "blue"))

Individual 1

Individual 2

Individual 3

Individual 4