Kriging using well logs to constrain permeability and porosity parame…

…ter distributions in flow tubes (#155)

* Add changelog template and CHANGELOG.md

* Add changelog item

* Updated contributing md

* Initial commit to add well_log information

* Incorporate kriging results for permeability in ert setup

* Changes due to new black version

* Added CHANGELOG entry

* Implemented on porosity as well

* Made mypy happy

* Minor updates in docstring and coding style

* Updates in code style due to update of Black

* Updated test_one_dimensional black

* Fixed too many eqlnum/satum bug

* Black

* Add kriging parameters to config

* Include kriging option in config

* Make sure empty dict is None

* Fix typo in simulationTime()

* Fix error in usage of variance instead of (two) standard deviations

* black

* Typo in doc

* Remove defaults in description

* Set semeio to 0.5.6 to test

* Ooops, wrong PR

* Change well log reading to use ecl2df

* Test CI

* black

* Add List type import

* Revert CI

CHANGELOG.md

@@ -11,6 +11,7 @@ This project adheres to [Semantic Versioning](https://semver.org/).
 - [#197](https://github.com/equinor/flownet/pull/197) Added opening/closing of connections based on straddling/plugging/perforations through time. This PR also adds a new perforation strategy `multiple_based_on_workovers` which models the well connections with as few connections as possible taking into account opening/closing of groups of connections. 
 - [#188](https://github.com/equinor/flownet/pull/188) The possibility to extract regions from an input simulation model extended to also include SATNUM regions. For relative permeability, the scheme keyword can be set to 'regions_from_sim' in the configuration file.
 - [#189](https://github.com/equinor/flownet/pull/189) User can now provide both a _base_ configuration file, and an optional extra configuration file which will be used to override the base settings.
+- [#155](https://github.com/equinor/flownet/pull/155) Adds reading of simulation 'well logs' to condition the priors of permeability and porosity based using kriging
 
 ### Changed
 - [#199](https://github.com/equinor/flownet/pull/199) Removed deprecated parameters in pyscal ('krowend', 'krogend') from config file. Added 'kroend' to config file. 

setup.py

@@ -14,6 +14,7 @@
     "opm>=2020.10",
     "pandas~=1.0",
     "pyarrow>=0.14",
+    "pykrige>=1.5",
     "pyscal~=0.6.1",
     "pyvista>=0.23",
     "pyyaml>=5.2",

src/flownet/ahm/_run_ahm.py

@@ -1,5 +1,5 @@
 import argparse
-from typing import Union, List, Optional
+from typing import Union, List, Optional, Dict
 import json
 import pathlib
 from operator import add
@@ -14,6 +14,7 @@
 from ..network_model import NetworkModel
 from ..network_model import create_connections
 from ._assisted_history_matching import AssistedHistoryMatching
+from ..utils import kriging
 
 from ..parameters import (
     PorvPoroTrans,
@@ -215,6 +216,99 @@ def _get_distribution(
     return df
 
 
+def _constrain_using_well_logs(
+    porv_poro_trans_dist_values: pd.DataFrame,
+    data: np.ndarray,
+    network: NetworkModel,
+    measurement_type: str,
+    config: ConfigSuite.snapshot,
+) -> pd.DataFrame:
+    """
+    Function to constrain permeability and porosity distributions of flow tubes by using 3D kriging of
+    porosity and permeability values from well logs.
+
+    Args:
+        porv_poro_trans_dist_values: pre-constraining dataframe
+        data: well log data
+        network: FlowNet network model
+        measurement_type: 'prorosity' or 'permeability' (always log10)
+        config: FlowNet configparser snapshot
+
+    Returns:
+        Well-log constrained "porv_poro_trans_dist_values" DataFrame
+
+    """
+    n = config.flownet.constraining.kriging.n
+    n_lags = config.flownet.constraining.kriging.n_lags
+    anisotropy_scaling_z = config.flownet.constraining.kriging.anisotropy_scaling_z
+    variogram_model = config.flownet.constraining.kriging.variogram_model
+
+    if measurement_type == "permeability":
+        data[:, 3] = np.log10(data[:, 3])
+
+    variogram_parameters: Optional[Dict] = None
+    if measurement_type == "permeability":
+        variogram_parameters = dict(
+            config.flownet.constraining.kriging.permeability_variogram_parameters
+        )
+    elif measurement_type == "porosity":
+        variogram_parameters = dict(
+            config.flownet.constraining.kriging.porosity_variogram_parameters
+        )
+
+    if not variogram_parameters:
+        variogram_parameters = None
+
+    k3d3_interpolator, ss3d_interpolator = kriging.execute(
+        data,
+        n=n,
+        n_lags=n_lags,
+        variogram_model=variogram_model,
+        variogram_parameters=variogram_parameters,
+        anisotropy_scaling_z=anisotropy_scaling_z,
+    )
+
+    parameter_min_kriging = k3d3_interpolator(
+        network.connection_midpoints
+    ) - 2 * np.sqrt(ss3d_interpolator(network.connection_midpoints))
+
+    parameter_max_kriging = k3d3_interpolator(
+        network.connection_midpoints
+    ) + 2 * np.sqrt(ss3d_interpolator(network.connection_midpoints))
+
+    if measurement_type == "permeability":
+        parameter_min_kriging = np.power(10, parameter_min_kriging)
+        parameter_max_kriging = np.power(10, parameter_max_kriging)
+
+    parameter_min = np.maximum(
+        np.minimum(
+            parameter_min_kriging,
+            porv_poro_trans_dist_values[f"maximum_{measurement_type}"].values,
+        ),
+        porv_poro_trans_dist_values[f"minimum_{measurement_type}"].values,
+    )
+    parameter_max = np.minimum(
+        np.maximum(
+            parameter_max_kriging,
+            porv_poro_trans_dist_values[f"minimum_{measurement_type}"].values,
+        ),
+        porv_poro_trans_dist_values[f"maximum_{measurement_type}"].values,
+    )
+
+    # Set NaN's to the full original range as set in the config
+    parameter_min[np.isnan(parameter_min)] = porv_poro_trans_dist_values[
+        f"minimum_{measurement_type}"
+    ].values[0]
+    parameter_max[np.isnan(parameter_max)] = porv_poro_trans_dist_values[
+        f"maximum_{measurement_type}"
+    ].values[0]
+
+    porv_poro_trans_dist_values[f"minimum_{measurement_type}"] = parameter_min
+    porv_poro_trans_dist_values[f"maximum_{measurement_type}"] = parameter_max
+
+    return porv_poro_trans_dist_values
+
+
 def update_distribution(
     parameters: List[Parameter], ahm_case: pathlib.Path
 ) -> List[Parameter]:
@@ -322,6 +416,13 @@ def run_flownet_history_matching(
     df_production_data: pd.DataFrame = field_data.production
     df_well_connections: pd.DataFrame = field_data.well_connections
 
+    # Load log data if required
+    df_well_logs: Optional[pd.DataFrame] = (
+        field_data.well_logs
+        if config.flownet.data_source.simulation.well_logs
+        else None
+    )
+
     # Load fault data if required
     df_fault_planes: Optional[pd.DataFrame] = (
         field_data.faults if config.model_parameters.fault_mult else None
@@ -365,6 +466,23 @@ def run_flownet_history_matching(
         network.grid.model.unique(),
     )
 
+    if df_well_logs is not None:
+        # Use well logs to constrain priors.
+
+        perm_data = df_well_logs[["X", "Y", "Z", "PERM"]].values
+        poro_data = df_well_logs[["X", "Y", "Z", "PORO"]].values
+
+        porv_poro_trans_dist_values = _constrain_using_well_logs(
+            porv_poro_trans_dist_values,
+            perm_data,
+            network,
+            "permeability",
+            config=config,
+        )
+        porv_poro_trans_dist_values = _constrain_using_well_logs(
+            porv_poro_trans_dist_values, poro_data, network, "porosity", config=config
+        )
+
     #########################################
     # Relative Permeability                 #
     #########################################

src/flownet/config_parser/_config_parser.py

@@ -56,7 +56,7 @@ def _to_lower(input_data: Union[List[str], str]) -> Union[List[str], str]:
     @configsuite.transformation_msg("Convert input string to absolute path")
     def _to_abs_path(path: Optional[str]) -> str:
         """
-        Helper function for the configsuite. Take in a path as a string and
+        Helper function for the configsuite. Takes in a path as a string and
         attempts to convert it to an absolute path.
 
         Args:
@@ -93,6 +93,10 @@ def _to_abs_path(path: Optional[str]) -> str:
                                         MK.Transformation: _to_abs_path,
                                         MK.Description: "Simulation input case to be used as data source for FlowNet",
                                     },
+                                    "well_logs": {
+                                        MK.Type: types.Bool,
+                                        MK.AllowNone: True,
+                                    },
                                     "vectors": {
                                         MK.Type: types.NamedDict,
                                         MK.Description: "Which vectors to use as observation data sources",
@@ -209,6 +213,62 @@ def _to_abs_path(path: Optional[str]) -> str:
                             "concave_hull": {MK.Type: types.Bool, MK.AllowNone: True},
                         },
                     },
+                    "constraining": {
+                        MK.Type: types.NamedDict,
+                        MK.Content: {
+                            "kriging": {
+                                MK.Type: types.NamedDict,
+                                MK.Content: {
+                                    "enabled": {
+                                        MK.Type: types.Bool,
+                                        MK.Default: False,
+                                        MK.Description: "Switch to enable or disable kriging on well log data.",
+                                    },
+                                    "n": {
+                                        MK.Type: types.Integer,
+                                        MK.Default: 20,
+                                        MK.Description: "Number of kriged values in each direct. E.g, n = 10 -> "
+                                        "10x10x10 = 1000 values.",
+                                    },
+                                    "n_lags": {
+                                        MK.Type: types.Integer,
+                                        MK.Default: 6,
+                                        MK.Description: "Number of averaging bins for the semivariogram.",
+                                    },
+                                    "anisotropy_scaling_z": {
+                                        MK.Type: types.Number,
+                                        MK.Default: 10,
+                                        MK.Description: "Scalar stretching value to take into account anisotropy. ",
+                                    },
+                                    "variogram_model": {
+                                        MK.Type: types.String,
+                                        MK.Default: "spherical",
+                                        MK.Description: "Specifies which variogram model to use. See PyKridge "
+                                        "documentation for valid options.",
+                                    },
+                                    "permeability_variogram_parameters": {
+                                        MK.Type: types.Dict,
+                                        MK.Description: "Parameters that define the specified variogram model. "
+                                        "Permeability model sill and nugget are in log scale. See "
+                                        "PyKridge documentation for valid options.",
+                                        MK.Content: {
+                                            MK.Key: {MK.Type: types.String},
+                                            MK.Value: {MK.Type: types.Number},
+                                        },
+                                    },
+                                    "porosity_variogram_parameters": {
+                                        MK.Type: types.Dict,
+                                        MK.Description: "Parameters that define the specified variogram model. See "
+                                        "PyKridge documentation for valid options.",
+                                        MK.Content: {
+                                            MK.Key: {MK.Type: types.String},
+                                            MK.Value: {MK.Type: types.Number},
+                                        },
+                                    },
+                                },
+                            },
+                        },
+                    },
                     "pvt": {
                         MK.Type: types.NamedDict,
                         MK.Content: {
@@ -1187,4 +1247,13 @@ def parse_config(
             "'min', 'max' and 'log_unif'. Currently one or more parameters are missing."
         )
 
+    if (
+        config.flownet.constraining.kriging.enabled
+        and not config.flownet.data_source.simulation.well_logs
+    ):
+        raise ValueError(
+            "Ambiguous configuration input: well log data needs to be loaded (from the simulation model) in order "
+            "to allow for enabling of kriging."
+        )
+
     return config

src/flownet/data/from_flow.py

@@ -1,6 +1,6 @@
 import warnings
 from pathlib import Path
-from typing import Union
+from typing import Union, List
 
 import numpy as np
 import pandas as pd
@@ -34,6 +34,7 @@ def __init__(
 
         self._input_case: Path = Path(input_case)
         self._eclsum = EclSum(str(self._input_case))
+        self._init = EclFile(str(self._input_case.with_suffix(".INIT")))
         self._grid = EclGrid(str(self._input_case.with_suffix(".EGRID")))
         self._restart = EclFile(str(self._input_case.with_suffix(".UNRST")))
         self._init = EclInitFile(self._grid, str(self._input_case.with_suffix(".INIT")))
@@ -88,6 +89,44 @@ def _well_connections(self) -> pd.DataFrame:
 
         return perforation_strategy_method(df).sort_values(["DATE"])
 
+    def _well_logs(self) -> pd.DataFrame:
+        """
+        Function to extract well log information from a Flow simulation.
+
+        Returns:
+            columns: WELL_NAME, X, Y, Z, PERM (mD), PORO (-)
+
+        """
+        coords: List = []
+
+        for well_name in self._wells["WELL"].unique():
+            unique_connections = self._wells[
+                self._wells["WELL"] == well_name
+            ].drop_duplicates(subset=["I", "J", "K1", "K2"])
+            for _, connection in unique_connections.iterrows():
+                ijk = (connection["I"] - 1, connection["J"] - 1, connection["K1"] - 1)
+                xyz = self._grid.get_xyz(ijk=ijk)
+
+                perm_kw = self._init.iget_named_kw("PERMX", 0)
+                poro_kw = self._init.iget_named_kw("PORO", 0)
+
+                coords.append(
+                    [
+                        well_name,
+                        *xyz,
+                        perm_kw[
+                            self._grid.cell(i=ijk[0], j=ijk[1], k=ijk[2]).active_index
+                        ],
+                        poro_kw[
+                            self._grid.cell(i=ijk[0], j=ijk[1], k=ijk[2]).active_index
+                        ],
+                    ]
+                )
+
+        return pd.DataFrame(
+            coords, columns=["WELL_NAME", "X", "Y", "Z", "PERM", "PORO"]
+        )
+
     def _production_data(self) -> pd.DataFrame:
         """
         Function to read production data for all producers and injectors from an
@@ -295,6 +334,11 @@ def well_connections(self) -> pd.DataFrame:
         """dataframe with all well connection coordinates"""
         return self._well_connections()
 
+    @property
+    def well_logs(self) -> pd.DataFrame:
+        """dataframe with all well log"""
+        return self._well_logs()
+
     @property
     def grid(self) -> EclGrid:
         """the simulation grid with properties"""

src/flownet/data/from_source.py

@@ -29,3 +29,10 @@ def faults(self) -> Optional[pd.DataFrame]:
         raise NotImplementedError(
             "The faults property is required to be implemented in a FromSource class."
         )
+
+    @property
+    @abstractmethod
+    def well_logs(self) -> Optional[pd.DataFrame]:
+        raise NotImplementedError(
+            "The well_log property is required to be implemented in a FromSource class."
+        )

src/flownet/network_model/_network_model.py

@@ -97,6 +97,28 @@ def get_connection_midpoints(self, i: Optional[int] = None) -> np.ndarray:
 
         return (coordinates_start + coordinates_end) / 2
 
+    @property
+    def connection_midpoints(self) -> np.ndarray:
+        """
+        Returns a numpy array with the midpoint of each connection in the network
+
+        Returns:
+            (Nx3) np.ndarray with connection midpoint coordinates.
+
+        """
+        coordinates_start = self._df_entity_connections[
+            ["xstart", "ystart", "zstart"]
+        ].values
+        coordinates_end = self._df_entity_connections[
+            [
+                "xend",
+                "yend",
+                "zend",
+            ]
+        ].values
+
+        return (coordinates_start + coordinates_end) / 2
+
     @property
     def aquifers_xyz(self) -> List[Coordinate]:
         """