diff --git a/changelog-entries/684.md b/changelog-entries/684.md
new file mode 100644
index 000000000..dc4f8adf8
--- /dev/null
+++ b/changelog-entries/684.md
@@ -0,0 +1 @@
+- Added PCE-base surrogate for micro-dumux in two-scale heat conduction [#684](https://github.com/precice/tutorials/pull/684)
diff --git a/tools/cleaning-tools.sh b/tools/cleaning-tools.sh
index eaa4758fd..b8003eee7 100755
--- a/tools/cleaning-tools.sh
+++ b/tools/cleaning-tools.sh
@@ -199,6 +199,7 @@ clean_dumux() {
         echo "- Cleaning up DuMuX case in $(pwd)"
         rm -fv ./*.vtu
         rm -fv ./*.pvd
+	rm -fv ./*.hdf5
         clean_precice_logs .
         clean_case_logs .
    )
diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/clean.sh b/two-scale-heat-conduction/micro-dumux-surrogate/clean.sh
new file mode 100755
index 000000000..3e9fd4e51
--- /dev/null
+++ b/two-scale-heat-conduction/micro-dumux-surrogate/clean.sh
@@ -0,0 +1,6 @@
+#!/usr/bin/env sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+clean_dumux .
diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/micro-manager-model-switching-config.json b/two-scale-heat-conduction/micro-dumux-surrogate/micro-manager-model-switching-config.json
new file mode 100644
index 000000000..860327c02
--- /dev/null
+++ b/two-scale-heat-conduction/micro-dumux-surrogate/micro-manager-model-switching-config.json
@@ -0,0 +1,22 @@
+{
+   "micro_file_names": ["micro_sim_sur", "micro_sim"],
+   "coupling_params": {
+      "precice_config_file_name": "../precice-config.xml",
+      "macro_mesh_name": "Macro-Mesh",
+      "write_data_names": ["K00", "K11", "Porosity"],
+      "read_data_names": ["Concentration"]
+   },
+   "simulation_params": {
+      "micro_dt": 0.01,
+      "macro_domain_bounds": [0.0, 1.0, 0.0, 0.5],
+      "decomposition": [2, 1],
+      "adaptivity": false,
+      "model_adaptivity": true,
+      "model_adaptivity_settings": {
+         "switching_function": "switch-model"
+      }
+    },
+    "diagnostics": {
+      "data_from_micro_sims": ["grain_size"]
+    }
+}
diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/micro-manager-snapshot-config.json b/two-scale-heat-conduction/micro-dumux-surrogate/micro-manager-snapshot-config.json
new file mode 100644
index 000000000..3a527daad
--- /dev/null
+++ b/two-scale-heat-conduction/micro-dumux-surrogate/micro-manager-snapshot-config.json
@@ -0,0 +1,15 @@
+{
+    "micro_file_names": ["micro_sim"],
+    "coupling_params": {
+            "parameter_file_name": "macro-concentration-samples.hdf5",
+            "write_data_names": ["K00", "K11", "Porosity"],
+            "read_data_names": ["Concentration"]
+    },
+    "simulation_params": {
+      "micro_dt": 0.01
+    },
+    "snapshot_params": {
+      "output_file_name": "micro-dumux-snapshots",
+      "initialize_once": false
+    }
+}
diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/micro_sim.pc.in b/two-scale-heat-conduction/micro-dumux-surrogate/micro_sim.pc.in
new file mode 100644
index 000000000..566aba5e4
--- /dev/null
+++ b/two-scale-heat-conduction/micro-dumux-surrogate/micro_sim.pc.in
@@ -0,0 +1,15 @@
+prefix=@prefix@
+exec_prefix=@exec_prefix@
+libdir=@libdir@
+includedir=@includedir@
+CXX=@CXX@
+CC=@CC@
+DEPENDENCIES=@REQUIRES@
+
+Name: @PACKAGE_NAME@
+Version: @VERSION@
+Description: micro_sim module
+URL: http://dune-project.org/
+Requires: dumux-phasefield dumux-precice
+Libs: -L${libdir}
+Cflags: -I${includedir}
diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/micro_sim_sur.py b/two-scale-heat-conduction/micro-dumux-surrogate/micro_sim_sur.py
new file mode 100644
index 000000000..0a6a6d19e
--- /dev/null
+++ b/two-scale-heat-conduction/micro-dumux-surrogate/micro_sim_sur.py
@@ -0,0 +1,86 @@
+"""
+Micro simulation Surrogate, requrie previous computation of surrogate model
+"""
+import joblib
+import numpy as np
+import math
+
+
+class MicroSimulation:
+
+    def __init__(self, sim_id):
+        """
+        Get the micro-scale model from the BayesValidRox surrogate model.
+
+        Parameters
+        ----------
+        sim_id : int
+            The simulation ID for the micro-scale simulation.
+        """
+        self._sim_id = sim_id
+        self._state = None
+
+        self._model = None
+        with open('micro-dumux-surrogate.pkl', 'rb') as input:
+            self._model = joblib.load(input)
+        if self._model is None:
+            raise RuntimeError("Failed to load model.")
+
+    def initialize(self):
+        output_data = dict()
+        output_data["K00"] = 0.4912490635619572
+        output_data["K11"] = 0.4912490635989945
+        output_data["Porosity"] = 0.4933482661391027
+
+        if self._sim_id == 0:
+            output_data["K00"] = 0.4912490640081466
+            output_data["K11"] = 0.4912490640081367
+
+        return output_data
+
+    def get_state(self):
+        """
+        Get the current state of the micro-scale simulation.
+
+        Returns
+        -------
+        state : dict
+            The current state of the micro-scale simulation.
+        """
+        return self._state
+
+    def set_state(self, state):
+        """
+        Set the current state of the micro-scale simulation.
+
+        Parameters
+        ----------
+        state : dict
+            The state to set for the micro-scale simulation.
+        """
+        self._state = state
+
+    def solve(self, macro_data, dt):
+        """
+        Solve the micro-scale simulation using the surrogate model.
+
+        Parameters
+        ----------
+        macro_data : dict
+            The macro-scale data required for the micro-scale simulation.
+        dt : float
+            The time step for the micro-scale simulation.
+
+        Returns
+        -------
+        output_data : dict
+            The output data from the micro-scale simulation.
+        """
+        model_eval, _ = self._model.eval_metamodel(np.array([macro_data["Concentration"]])[:, np.newaxis])
+        output_data = dict()
+        output_data["K00"] = model_eval["K00"][0][0]
+        output_data["K11"] = model_eval["K11"][0][0]
+        output_data["Porosity"] = model_eval["Porosity"][0][0]
+        output_data["grain_size"] = math.sqrt((1 - model_eval["Porosity"][0][0]) / math.pi)
+
+        return output_data
diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/model.py b/two-scale-heat-conduction/micro-dumux-surrogate/model.py
new file mode 100644
index 000000000..5284f7485
--- /dev/null
+++ b/two-scale-heat-conduction/micro-dumux-surrogate/model.py
@@ -0,0 +1,5 @@
+# Dummy model function for the micro-dumux surrogate.
+# We do not wrap the original DuMuX model because we will directly provide
+# snapshots (computed by the Micro Manager) to BayesValidRox.
+def model(samples):
+    return None
diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/params.input b/two-scale-heat-conduction/micro-dumux-surrogate/params.input
new file mode 100644
index 000000000..964a6896d
--- /dev/null
+++ b/two-scale-heat-conduction/micro-dumux-surrogate/params.input
@@ -0,0 +1,25 @@
+[Assembly]
+Multithreading = false
+
+[TimeLoop]
+TEnd = 0.25                 # end time of the simulation
+DtInitial = 0.01            # initial time step size
+MaxTimeStepSize = 0.01      # maximal time step size
+
+[Grid]
+LowerLeft = 0.0 0.0         # lower left (front) corner of the domain (keep this fixed at 0 0!)
+UpperRight = 1.0 1.0        # upper right (back) corner of the domain
+Cells = 80 80               # grid resolution in each coordinate direction
+Periodic = 1 1              # Periodic Boundary conditions in both dimensions
+
+[Problem]
+xi = 0.08                   # phasefield parameter (lambda, set to around 4/Ncells) 
+omega = 0.01                # phasefield diffusivity/surface tension parameter (gamma)
+kt = 1.0                    # constant deciding speed of expansion/contraction
+eqconc = 0.5                # equilibrium concentration 
+ks = 1.0                    # conductivity of sand material
+kg = 0.0                    # conductivity of void material
+Name = cell_phase           # base name for VTK output files
+Radius = 0.4                # initial radius of the grain
+PhasefieldICScaling = 4.0   # factor in initial phasefield function
+MaxPorosity = 0.9686        # porosity cap
diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/requirements.txt b/two-scale-heat-conduction/micro-dumux-surrogate/requirements.txt
new file mode 100644
index 000000000..8970d34da
--- /dev/null
+++ b/two-scale-heat-conduction/micro-dumux-surrogate/requirements.txt
@@ -0,0 +1,4 @@
+numpy
+bayesvalidrox
+pyprecice
+micro-manager-precice
diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/run-surrogate-workflow.sh b/two-scale-heat-conduction/micro-dumux-surrogate/run-surrogate-workflow.sh
new file mode 100755
index 000000000..ddf636773
--- /dev/null
+++ b/two-scale-heat-conduction/micro-dumux-surrogate/run-surrogate-workflow.sh
@@ -0,0 +1,9 @@
+#!/usr/bin/env bash
+set -e -u
+
+python3 -m venv .venv
+. .venv/bin/activate
+
+pip install -r requirements.txt
+
+python surrogate_workflow.py
diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/run.sh b/two-scale-heat-conduction/micro-dumux-surrogate/run.sh
new file mode 100755
index 000000000..7736ee581
--- /dev/null
+++ b/two-scale-heat-conduction/micro-dumux-surrogate/run.sh
@@ -0,0 +1,29 @@
+#!/usr/bin/env bash
+set -e -u
+
+. ../../tools/log.sh
+exec > >(tee --append "$LOGFILE") 2>&1
+
+usage() { echo "Usage: cmd [-s] [-p n]" 1>&2; exit 1; }
+
+# Check if no input argument was provided
+if [ -z "$*" ] ; then
+  echo "No input argument provided. Micro Manager is launched in serial"
+  micro-manager-precice micro-manager-model-switching-config.json
+fi
+
+while getopts ":sp" opt; do
+  case ${opt} in
+  s)
+    micro-manager-precice micro-manager-model-switching-config.json
+    ;;
+  p)
+    mpiexec -n "$2" micro-manager-precice micro-manager-model-switching-config.json
+    ;;
+  *)
+    usage
+    ;;
+  esac
+done
+
+close_log
diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/surrogate_workflow.py b/two-scale-heat-conduction/micro-dumux-surrogate/surrogate_workflow.py
new file mode 100644
index 000000000..69581a747
--- /dev/null
+++ b/two-scale-heat-conduction/micro-dumux-surrogate/surrogate_workflow.py
@@ -0,0 +1,194 @@
+import os
+import subprocess
+from bayesvalidrox import PyLinkForwardModel, Input, PCE, ExpDesigns, Engine
+import h5py
+import joblib
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def create_snapshots() -> None:
+    """
+    Create snapshots of the DuMuX model in micro-dumux/ using the Micro Manager.
+    """
+    # What effect do uniformly distributed concentration samples have on the quality of the surrogate model?
+    concentration_samples = np.linspace(0.0, 0.9, 500)
+
+    # The name of the HDF5 file should be the same as the parameter_file_name in micro-manager-snapshot-config.json
+    with h5py.File("macro-concentration-samples.hdf5", "w") as f:
+        f.create_dataset("Concentration", data=concentration_samples)
+
+    # Run the Micro Manager to create snapshot database
+    subprocess.run(
+        ["mpiexec", "-n", "2", "micro-manager-precice", "--snapshot", "micro-manager-snapshot-config.json"],
+        check=True,
+    )
+
+
+def read_snapshots() -> tuple:
+    """
+    Read the snapshots created by the Micro Manager and return the inputs and outputs.
+    The inputs are the concentration samples, and the outputs are the porosity and conductivity matrix values.
+
+    Returns
+    -------
+    tuple
+        A tuple containing the inputs (concentration samples) and outputs (porosity and conductivity values).
+    """
+    with h5py.File("micro-dumux-snapshots.hdf5", "r") as f:
+        concentration_data = f["Concentration"][:]
+        porosity_data = f["Porosity"][:]
+        k_00_data = f["K00"][:]
+        k_11_data = f["K11"][:]
+
+    concentration = np.swapaxes(np.array([concentration_data]), 0, 1)
+    porosity = np.swapaxes(np.array([porosity_data]), 0, 1)
+    k_00 = np.swapaxes(np.array([k_00_data]), 0, 1)
+    k_11 = np.swapaxes(np.array([k_11_data]), 0, 1)
+
+    outputs = {"Porosity": porosity, "K00": k_00, "K11": k_11, "x_values": np.array([0])}
+    return concentration, outputs
+
+
+def split_samples(X, y, n_valid):
+    """
+    Split the samples and evaluations into training and validation/test data.
+    The split is performed randomly.
+
+    Parameters
+    ----------
+    X : np.ndarray
+        Samples, shape (#samples, #parameters)
+    y : dict
+        Corresponding model evaluations. Expected to match BVR output format.
+    n_valid : int
+        Number of samples to keep for validation.
+
+    Returns
+    -------
+    X_train, y_train,
+    X_valid, y_valid
+    """
+    n_samples = X.shape[0]
+    if n_valid >= n_samples:
+        raise AttributeError('The set number of validation points is invalid.')
+
+    # Random split
+    n_train = n_samples - n_valid
+    choice = np.random.choice(
+        range(n_samples), size=(n_train,), replace=False
+    )
+    ind = np.zeros(n_samples, dtype=bool)
+    ind[choice] = True
+
+    # Split samples
+    X_train = X[ind]
+    X_valid = X[~ind]
+
+    # Split outputs
+    y_train = {}
+    y_valid = {}
+    for key in y:
+        if key != "x_values":
+            y_train[key] = y[key][ind]
+            y_valid[key] = y[key][~ind]
+
+    return X_train, y_train, X_valid, y_valid
+
+
+def create_surrogate() -> tuple:
+    """
+    Create a surrogate model from the Micro Manager snapshots.
+
+    Returns
+    -------
+    tuple
+        A tuple containing the validation samples and corresponding model evaluations.
+    """
+    # We create a fake model from model.py because we directly provide the
+    # input and outputs from the Micro Manager snapshots.
+    model = PyLinkForwardModel()
+    model.py_file = "model"
+    model.name = "micro-dumux-surrogate"
+    model.link_type = "function"
+    model.output_names = ["Porosity", "K00", "K11"]
+
+    x, y = read_snapshots()
+
+    # Split the samples into training and validation sets
+    n_valid = 200
+    x_train, y_train, x_valid, y_valid = split_samples(x, y, n_valid)
+
+    inputs = Input()
+    inputs.add_marginals(name="Concentration", dist_type="unif", parameters=[0, 0.5])
+
+    exp_design = ExpDesigns(inputs)
+    exp_design.x = x_train
+    exp_design.y = y_train
+
+    # Create the surrogate model
+    meta_model = PCE(inputs)
+    meta_model.meta_model_type = "aPCE"
+    meta_model.pce_reg_method = "FastARD"
+    meta_model.pce_deg = 5
+
+    # Train the surrogate model
+    engine = Engine(meta_model, model, exp_design)
+    engine.train_normal()
+
+    with open(f'{model.name}.pkl', 'wb') as output:
+        joblib.dump(engine, output, 2)
+
+    return x_valid, y_valid
+
+
+def validate_surrogate(x_valid, y_valid, model_name="micro-dumux-surrogate.pkl") -> None:
+    """
+    Validate the surrogate model using the validation samples and outputs.
+
+    Parameters
+    ----------
+    x_valid : np.ndarray
+        Validation samples.
+    y_valid : dict
+        Corresponding model evaluations for validation samples.
+    model_name : str
+        Name of the surrogate model file.
+    """
+    with open(model_name, 'rb') as input:
+        engine = joblib.load(input)
+
+    y_metamod, _ = engine.eval_metamodel(x_valid)
+
+    # engine.plot_adapt(y_valid, y_metamod, y_metamod_std, x_valid)
+
+    # Compare predictions with true values
+    plt.figure()
+    plt.scatter(y_valid["Porosity"], y_metamod["Porosity"])
+    plt.xlabel("True Values")
+    plt.ylabel("Predictions")
+    plt.title(f"Validation: porosity")
+    plt.plot([0.4, 1.1], [0.4, 1.1], "k--")
+    plt.xlim(0.4, 1.1)
+    plt.ylim(0.4, 1.1)
+    plt.savefig("surrogate-validation-porosity.pdf", bbox_inches="tight")
+
+    plt.figure()
+    plt.scatter(x_valid, y_valid["Porosity"])
+    plt.xlabel("Validation Concentration")
+    plt.ylabel("Validation Porosity")
+    plt.title("IO: Concentration to Porosity")
+    plt.savefig("surrogate-validation-concentration-porosity.pdf", bbox_inches="tight")
+
+
+def main():
+    """
+    Create snaphots of the DuMuX model, create a surrogate model from the snapshots, and validate the surrogate model.
+    """
+    create_snapshots()
+    x_valid, y_valid = create_surrogate()
+    validate_surrogate(x_valid, y_valid)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/switch-model.py b/two-scale-heat-conduction/micro-dumux-surrogate/switch-model.py
new file mode 100644
index 000000000..eb2d99b0d
--- /dev/null
+++ b/two-scale-heat-conduction/micro-dumux-surrogate/switch-model.py
@@ -0,0 +1,48 @@
+import numpy as np
+from typing import Dict
+
+
+def switching_function(model_res: int, location: np.ndarray, t: float, input: Dict, prev_output: Dict):
+    """
+    Determine which model to use based on the resolution, location, time, input, and previous output.
+
+    Parameters
+    ----------
+    model_res : int
+        The resolution level of the model (0 for full-order model, >0 for reduced-order model).
+    location : np.ndarray
+        The macro-scale coordinates corresponding to the micro-scale simulation.
+    t : float
+        The current time step.
+    input : Dict
+        The input data for the model, including concentration.
+    prev_output : dict
+        The output from the previous time step.
+
+    Returns
+    -------
+    model_res_offset : int
+        The offset to change the model used from the provided hierarchy of models.
+    """
+    model_res_offset = 0
+
+    # Only use the full-order model for the first time step
+    if t == 0.0:
+        if model_res > 0:
+            model_res_offset = -1
+        else:
+            model_res_offset = 0
+
+    # After the first time step, the model is selected dynamically based on
+    # the concentration value and the model resolution.
+    else:
+        concentration = input["Concentration"]
+        is_valid_range = 0.2 < concentration < 0.3
+        is_fom = model_res == 0
+
+        if is_fom and is_valid_range:
+            model_res_offset = 1
+        if not is_fom and not is_valid_range:
+            model_res_offset = -1
+
+    return model_res_offset
diff --git a/two-scale-heat-conduction/precice-config.xml b/two-scale-heat-conduction/precice-config.xml
index 028d0bf31..8e85953c2 100644
--- a/two-scale-heat-conduction/precice-config.xml
+++ b/two-scale-heat-conduction/precice-config.xml
@@ -16,6 +16,7 @@
   <data:scalar name="Grain-Size" />
   <data:scalar name="Active-State" />
   <data:scalar name="Active-Steps" />
+  <data:scalar name="model_resolution" />
 
   <mesh name="Macro-Mesh" dimensions="2">
     <use-data name="K00" />
@@ -25,6 +26,7 @@
     <use-data name="Grain-Size" />
     <use-data name="Active-State" />
     <use-data name="Active-Steps" />
+    <use-data name="model_resolution" />
   </mesh>
 
   <participant name="Macro">
@@ -44,6 +46,7 @@
     <write-data name="Grain-Size" mesh="Macro-Mesh" />
     <write-data name="Active-State" mesh="Macro-Mesh" />
     <write-data name="Active-Steps" mesh="Macro-Mesh" />
+    <write-data name="model_resolution" mesh="Macro-Mesh" />
     <export:vtu directory="precice-exports" />
   </participant>