diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py
index bf5fa17ac..2ecd711c1 100644
--- a/python/ribasim/ribasim/model.py
+++ b/python/ribasim/ribasim/model.py
@@ -525,7 +525,7 @@ def plot(self, ax=None, indicate_subnetworks: bool = True) -> Any:
         self.plot_control_listen(ax)
         self.network.node.plot(ax=ax, zorder=3)
         if indicate_subnetworks:
-            self.network.node.plot_allocation_networks(ax=ax)
+            self.network.node.plot_allocation_networks(ax=ax, zorder=4)
 
         ax.legend(loc="lower left", bbox_to_anchor=(1, 0.5))
 
diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py
index 7de7eb347..652acc790 100644
--- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py
+++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py
@@ -7,8 +7,9 @@
 import ribasim_testmodels
 from ribasim_testmodels.allocation import (
     allocation_example_model,
-    # looped_subnetwork_model,
     fractional_flow_subnetwork_model,
+    # looped_subnetwork_model,
+    main_network_with_subnetworks_model,
     minimal_subnetwork_model,
     subnetwork_model,
     user_model,
@@ -85,6 +86,7 @@
     # Disable until this issue is resolved:
     # https://github.com/Deltares/Ribasim/issues/692
     # "looped_subnetwork_model",
+    "main_network_with_subnetworks_model",
 ]
 
 # provide a mapping from model name to its constructor, so we can iterate over all models
diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py
index c3416c4d4..4fa62a5e9 100644
--- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py
+++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py
@@ -123,7 +123,9 @@ def user_model():
 
 
 def subnetwork_model():
-    """Create a user testmodel representing a subnetwork."""
+    """Create a user testmodel representing a subnetwork.
+    This model is merged into main_network_with_subnetworks_model.
+    """
 
     # Setup the nodes:
     xy = np.array(
@@ -283,7 +285,9 @@ def subnetwork_model():
 
 
 def looped_subnetwork_model():
-    """Create a user testmodel representing a subnetwork containing a loop in the topology."""
+    """Create a user testmodel representing a subnetwork containing a loop in the topology.
+    This model is merged into main_network_with_subnetworks_model.
+    """
     # Setup the nodes:
     xy = np.array(
         [
@@ -670,7 +674,9 @@ def minimal_subnetwork_model():
 
 
 def fractional_flow_subnetwork_model():
-    """Create a small subnetwork that contains fractional flow nodes."""
+    """Create a small subnetwork that contains fractional flow nodes.
+    This model is merged into main_network_with_subnetworks_model.
+    """
 
     xy = np.array(
         [
@@ -1065,3 +1071,102 @@ def allocation_example_model():
     )
 
     return model
+
+
+def main_network_with_subnetworks_model():
+    """Generate a model with a main network with connected subnetworks which are other test models."""
+
+    # Main network
+    n_basins = 5
+    n_nodes = 2 * n_basins
+    xy = np.array([(3 * i, (-1) ** (i + 1)) for i in range(n_nodes)])
+    node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])
+
+    node_type = ["FlowBoundary", "Basin"] + (n_basins - 1) * [
+        "LinearResistance",
+        "Basin",
+    ]
+    node = ribasim.Node(
+        df=gpd.GeoDataFrame(
+            data={"type": node_type, "allocation_network_id": 1},
+            index=pd.Index(np.arange(len(xy)) + 1, name="fid"),
+            geometry=node_xy,
+            crs="EPSG:28992",
+        )
+    )
+
+    from_id = np.arange(1, n_nodes)
+    to_id = np.arange(2, n_nodes + 1)
+    allocation_network_id = len(from_id) * [None]
+    allocation_network_id[0] = 1
+    lines = node.geometry_from_connectivity(from_id, to_id)
+    edge = ribasim.Edge(
+        df=gpd.GeoDataFrame(
+            data={
+                "from_node_id": from_id,
+                "to_node_id": to_id,
+                "edge_type": len(from_id) * ["flow"],
+                "allocation_network_id": allocation_network_id,
+            },
+            geometry=lines,
+            crs="EPSG:28992",
+        )
+    )
+
+    # Setup the basins:
+    basin_node_ids = np.arange(2, n_nodes + 1, 2)
+    profile = pd.DataFrame(
+        data={
+            "node_id": np.repeat(basin_node_ids, [2] * len(basin_node_ids)),
+            "area": 1000.0,
+            "level": n_basins * [0.0, 1.0],
+        }
+    )
+
+    static = pd.DataFrame(
+        data={
+            "node_id": basin_node_ids,
+            "drainage": 0.0,
+            "potential_evaporation": 0.0,
+            "infiltration": 0.0,
+            "precipitation": 0.0,
+            "urban_runoff": 0.0,
+        }
+    )
+
+    state = pd.DataFrame(data={"node_id": basin_node_ids, "level": 1.0})
+
+    basin = ribasim.Basin(profile=profile, static=static, state=state)
+
+    flow_boundary = ribasim.FlowBoundary(
+        static=pd.DataFrame(
+            data={
+                "node_id": [1],
+                "flow_rate": 1.0,
+            }
+        )
+    )
+
+    linear_resistance = ribasim.LinearResistance(
+        static=pd.DataFrame(
+            data={"node_id": np.arange(3, n_nodes + 1, 2), "resistance": 1e-3}
+        )
+    )
+
+    model = ribasim.Model(
+        network=ribasim.Network(
+            node=node,
+            edge=edge,
+        ),
+        basin=basin,
+        flow_boundary=flow_boundary,
+        linear_resistance=linear_resistance,
+        starttime="2020-01-01 00:00:00",
+        endtime="2021-01-01 00:00:00",
+    )
+
+    model.merge_model(subnetwork_model(), offset_spacial=(0.0, 3.0))
+    model.merge_model(fractional_flow_subnetwork_model(), offset_spacial=(14.0, 3.0))
+    model.merge_model(looped_subnetwork_model(), offset_spacial=(26.0, 3.0))
+
+    return model