diff --git a/doc/body.tex b/doc/body.tex
index fe1e0bc5..bd217ab9 100644
--- a/doc/body.tex
+++ b/doc/body.tex
@@ -300,4 +300,8 @@
 \insection
 \input{sections/ex-gwe-prt.tex}
 
+\clearpage
+\insection
+\input{sections/ex-gwt-adv-schemes.tex}
+
 %MODFLOW API GWT examples
diff --git a/doc/sections/ex-gwt-adv-schemes.tex b/doc/sections/ex-gwt-adv-schemes.tex
new file mode 100644
index 00000000..cf3f1cf4
--- /dev/null
+++ b/doc/sections/ex-gwt-adv-schemes.tex
@@ -0,0 +1,114 @@
+\section{Advection schemes in MODFLOW6}
+
+This example demonstrates the performance of different numerical advection schemes when solving the groundwater transport equation under pure advection conditions. The example provides a comparison of four advection schemes (upstream, central difference, TVD, and UTVD) across three different grid types (structured, triangular, and Voronoi) using three inflow concentration profiles with varying numerical challenges.
+
+\subsection{Example description}
+
+The problem solves the pure advection equation:
+
+\begin{equation}
+\frac{\partial \left(S_w \theta C\right)}{\partial t} = -\vec{\nabla} \cdot \left(\vec{q} \cdot C\right)
+\label{eq:pureadvection}
+\end{equation}
+
+\noindent where $C$ is the volumetric concentration ($M/L^3$), $S_w$ is water saturation (dimensionless), $\theta$ is porosity (dimensionless), and $\mathbf{q}$ is the specific discharge field ($L/T$). The problem is configured with no dispersion or diffusion terms, making it ideal for testing numerical scheme performance since an analytical solution exists.
+
+The model domain consists of a 100 cm $\times$ 100 cm square. A uniform inflow (flux) of magnitude 0.5 cm/s is applied to the inflow boundaries (left and bottom edges) resulting in a specific discharge field of $\vec{q} = (q_x, q_y) = (0.354, 0.354)$ cm/s, which is oriented at 45 degrees. Prescribed concentrations are applied on the inflow boundaries,while CHD boundary conditions are used on the outflow boundaries which extract any outflowing concentration.
+
+The spatial discretization uses three different grid types to test scheme performance across different grid geometries:
+\begin{itemize}
+\item \textbf{Structured grid}: Grid using 50 $\times$ 50 rectangular cells (2 cm $\times$ 2 cm)
+\item \textbf{Triangular grid}: Grid using triangular cells with maximum area constraint
+\item \textbf{Voronoi grid}: Grid using Voronoi cells derived from triangular cells
+\end{itemize}
+
+\noindent The simulation time spans 300 seconds using adaptive time stepping with an initial time step of 5 seconds and a Courant number of 0.7. The Courant number is a dimensionless quantity representing the fraction of a cell that fluid travels in a single time step. A value less or equal to 1 is required for stability. A relatively high Courant number is chosen here specifically to challenge the central difference scheme and demonstrate its potential for oscillatory behavior on discontinuous functions. Four advection schemes are compared:
+\begin{itemize}
+\item \textbf{Upstream}: First-order accurate, stable but diffusive scheme
+\item \textbf{Central difference}: Second-order accurate but prone to oscillations on discontinuities
+\item \textbf{TVD}: Total Variation Diminishing scheme that handles sharp fronts well on structured grids
+\item \textbf{UTVD}: Unstructured TVD scheme extending TVD capabilities to unstructured grids
+\end{itemize}
+
+\noindent Three inflow concentration profiles are used to probe different numerical challenges:
+\begin{itemize}
+\item \textbf{Sin² wave}: Smooth function testing second-order accuracy
+\item \textbf{Step wave}: Sharp transition testing boundedness properties
+\item \textbf{Square wave}: Sharp rectangular discontinuity testing whether schemes can capture the original peak value of the analytical solution
+\end{itemize}
+
+\noindent These inflow concentration profiles are applied as concentration profiles on the inflow boundaries (left and bottom edges), where the concentration varies with position along the edge. In this pure advection problem these profiles should be transported along the streamlines without any change in shape. By comparing the simulation results to this expected behavior, we can evaluate how well each numerical scheme handles different challenges, such as the smooth gradients of the sin² wave versus the sharp discontinuities of the step and square waves. Any observed numerical diffusion or oscillations errors highlight the limitations of a given scheme.
+
+The model parameters are summarized in table~\ref{tab:ex-gwt-adv-schemes-01}.
+
+\input{../tables/ex-gwt-adv-schemes-01.tex}
+
+\subsection{Example Results}
+
+The example generates 39 total simulations: 3 flow simulations (one per grid type) and 36 transport simulations models (3 grids $\times$ 4 schemes $\times$ 3 inflow concentration profiles). Results are presented in three complementary views to highlight different aspects of scheme performance.
+
+\subsubsection{Flow Field}
+
+Figure~\ref{fig:ex-gwt-adv-schemes-flow} shows the computed head fields for all three grid types. The results verify that the flow field produces the expected uniform gradient with a 45° angle across all grid geometries, providing the necessary foundation for the pure advection transport simulations. Head contours are parallel straight lines as expected for uniform flow, confirming the correct implementation of the boundary conditions and hydraulic properties.
+
+\begin{StandardFigure}{
+    Computed head fields for the three grid types used in the advection scheme comparison. \textit{Left}: Structured rectangular grid. \textit{Center}: Triangular grid. \textit{Right}: Voronoi polygon grid. All grids produce the expected uniform head gradient at 45° angle.
+}{fig:ex-gwt-adv-schemes-flow}{../figures/ex-gwt-adv-schemes-flow.png}
+\end{StandardFigure}
+
+\subsubsection{Concentration Field Results}
+
+Figures~\ref{fig:ex-gwt-adv-schemes-sin²}, \ref{fig:ex-gwt-adv-schemes-step}, and \ref{fig:ex-gwt-adv-schemes-square} show the final concentration distributions at the end of the simulation, which are near the steady-state solution, for the three inflow concentration profiles across all scheme and grid combinations. These 2D concentration maps reveal the characteristic behavior of each numerical scheme, though the specific wave profiles are more clearly visualized in the cross-section analysis presented later. A mask has been applied to filter out all values below 0 and above 1. This has been done to make it easier to compare the results between different schemes.
+
+\paragraph{Sin² wave results} (Figure~\ref{fig:ex-gwt-adv-schemes-sin²}) reveal important differences in scheme performance for this smooth inflow concentration profile. While all schemes perform well on the structured grid, the upstream scheme is noticeable more diffusive. The central difference scheme becomes unstable on the triangular and Voronoi grids. On the structured grid, the central difference scheme is the least diffusive but exhibits undershoots resulting in negative concentrations due to the relatively high Courant number (0.7); lowering the Courant number eliminates this behavior. The UTVD scheme demonstrates robust performance across all grid types, maintaining both accuracy and stability.
+\begin{StandardFigure}{
+    Final concentration distributions for the sin² wave inflow concentration profile. Rows represent the four advection schemes (upstream, central difference, TVD, UTVD) and columns represent the three grid types (structured, triangular, Voronoi).
+}{fig:ex-gwt-adv-schemes-sin²}{../figures/ex-gwt-adv-schemes-sin²-wave-conc.png}
+\end{StandardFigure}
+
+\paragraph{Step wave results} (Figure~\ref{fig:ex-gwt-adv-schemes-step}) show similar trends to the sin² wave results. On the structured grid, the central difference scheme outperforms UTVD in preserving the overall shape of the wave. However, while it only exhibited undershoot for the smooth sin² wave, the sharp discontinuity of the step wave causes both undershoot and overshoot. For the step wave, the central difference scheme becomes unstable on the triangular and Voronoi grids. The TVD scheme outperforms upstream on the structured grid but performs similarly on the triangular and Voronoi grids.The UTVD scheme emerges as the least dissipative scheme that maintains stability across all grid types.
+\begin{StandardFigure}{
+    Final concentration distributions for the step wave inflow concentration profile. The central difference scheme begins to show oscillatory behavior near the sharp transition, while TVD and UTVD schemes maintain stable solutions without spurious oscillations.
+}{fig:ex-gwt-adv-schemes-step}{../figures/ex-gwt-adv-schemes-step-wave-conc.png}
+\end{StandardFigure}
+
+\paragraph{Square wave results} (Figure~\ref{fig:ex-gwt-adv-schemes-square}) exhibit similar behavior to the step wave results. The TVD scheme outperforms upstream on the structured grid but performs similarly on the triangular and Voronoi grids. The central difference scheme shows oscillations on the structured grid and blows up on the other grids. The UTVD scheme works well on all grids and is the only stable scheme that actually comes close to the maximum value of the analytical solution, demonstrating its capability for handling sharp discontinuities across all grid types.
+\begin{StandardFigure}{
+    Final concentration distributions for the square wave inflow concentration profile. This test with sharp discontinuities demonstrates the oscillatory behavior of the central difference scheme and the ability of the TVD/UTVD schemes to maintain stable solutions without spurious oscillations.
+}{fig:ex-gwt-adv-schemes-square}{../figures/ex-gwt-adv-schemes-square-wave-conc.png}
+\end{StandardFigure}
+
+\subsubsection{Cross-Section Analysis}
+
+Figure~\ref{fig:ex-gwt-adv-schemes-conc-cross-section} presents concentration profiles along the centerline (y = 50 cm) for all inflow concentration profiles and grid types, comparing numerical solutions against the analytical solution. These 1D cross-sections provide quantitative assessment of scheme accuracy and enable detailed comparison of numerical behavior:
+
+\textbf{Accuracy assessment}: The cross-sections reveal that the UTVD scheme follows the analytical solution best for discontinuous functions, while maintaining computational stability. The TVD scheme only slightly outperforms the upstream scheme in terms of accuracy. The central difference scheme works well for the smooth function on a structured grid but its Courant condition is restricting, limiting its practical applicability.
+
+\textbf{Oscillation identification}: Clear oscillations in the central difference scheme are visible as concentration values exceeding the physical bounds [0,1] or exhibiting spurious undershoot/overshoot behavior.
+
+\textbf{Grid sensitivity}: The comparison across grid types demonstrates that UTVD maintains consistent performance regardless of grid geometry, while TVD may show some degradation on unstructured grids.
+
+\textbf{Diffusion quantification}: The upstream scheme shows characteristic smearing of sharp features, providing a baseline for assessing the improvement offered by higher-order schemes.
+
+\begin{StandardFigure}{
+    Concentration cross-sections at y = 50 cm comparing numerical solutions (dashed lines with markers) against analytical solutions (solid lines) for all inflow concentration profiles and grid types. \textit{Rows}: Inflow concentration profiles (sin² wave, step wave, square wave). \textit{Columns}: Grid types (structured, triangular, Voronoi). The restricted y-axis range [-0.1, 1.1] helps visualize oscillations and overshoots.
+}{fig:ex-gwt-adv-schemes-conc-cross-section}{../figures/ex-gwt-adv-schemes-conc-cross-section.png}
+\end{StandardFigure}
+
+\subsubsection{Key Findings}
+
+The results from these simulations provide several key findings that can be used as guidance for advection scheme selection in groundwater transport modeling. While this comparison is based on solid testing, it is also limited in scope and the findings should be interpreted accordingly:
+
+\begin{enumerate}
+\item \textbf{Grid-dependent performance}: The central difference scheme performed reasonably on the structured grid but became unstable on the triangular and Voronoi grids, particularly for the discontinuous input concentration profiles. This demonstrates the critical importance of considering grid type when selecting numerical schemes.
+
+\item \textbf{Courant number sensitivity}: The central difference scheme's stability was highly sensitive to the Courant number. At Courant = 0.7, it exhibited undershoots leading to negative concentrations on smooth functions and oscillartory instability for the discontinuous inflow concentration profiles. Reducing the Courant number mitigated these issues for the smooth inflow concentration profile but at the cost of computational efficiency.
+
+\item \textbf{UTVD scheme as overall best choice}: The UTVD scheme consistently outperformed all others across grid types and inflow concentration profiles .It was the least dissipative scheme that maintained stability, and is the only stable scheme that approached the maximum values of analytical solutions for sharp discontinuities.
+
+\item \textbf{TVD limitations on unstructured grids}: While TVD performed well on the structured grid and outperformed the upstream scheme for the discontinuous input concentration profiles, its advantage diminished on the unstructured grids, where it performed similarly to the upstream scheme.
+
+\item \textbf{Upstream scheme reliability}: The upstream scheme provided a stable baseline across all grid types and inflow concentration profiles, though with significant diffusion. It can serve as a reliable fallback option when stability is paramount.
+
+\item \textbf{Practical recommendations}: For unstructured grid applications involving transport with potential sharp fronts, the tests conducted here suggest that UTVD is the best choice. For structured grids, TVD offers a good balance of accuracy and stability. The central difference scheme should be used with caution and requires careful Courant number selection, particularly on unstructured grids. If computational speed is prioritized and users can accept more diffusive solutions, the upstream scheme provides a fast and stable option. Similarly, if users can tolerate potential instabilities in exchange for speed and higher accuracy on structured grids, the central difference scheme may be considered with appropriate Courant number restrictions.
+\end{enumerate}
diff --git a/environment.yml b/environment.yml
index 240ae418..72c8518a 100644
--- a/environment.yml
+++ b/environment.yml
@@ -14,6 +14,7 @@ dependencies:
   - flaky
   - fprettify
   - fortran-language-server
+  - geopandas 
   - gitpython
   - Jinja2>=3.1.5,<4
   - jupytext
diff --git a/scripts/ex-gwt-adv-schemes.py b/scripts/ex-gwt-adv-schemes.py
new file mode 100644
index 00000000..2eb41477
--- /dev/null
+++ b/scripts/ex-gwt-adv-schemes.py
@@ -0,0 +1,930 @@
+# %% [markdown]
+# ## Advection schemes in MODFLOW 6
+#
+# This example demonstrates the performance of different numerical advection schemes when solving the groundwater transport equation under pure advection conditions. We solve the pure advection equation:
+#
+# $$\frac{\partial \left(S_w \theta C\right)}{\partial t} = -\nabla \left(\mathbf{q} \cdot C\right)$$
+#
+# where C is concentration [g/cm³] and **q** is the specific discharge field = (qx, qy) = (0.354, 0.354) cm/s at 45°. The problem is configured with no dispersion or diffusion terms, making it a perfect test case for numerical scheme performance since an analytical solution exists.
+#
+# **Problem Setup:**
+# - Domain: 100cm x 100cm square with uniform flow at 45° angle
+# - Boundary conditions: Prescribed concentrations on inflow boundaries
+# - Time: 300 seconds with adaptive time stepping (initial dt = 5s)
+# - Physics: Pure advection without mixing processes (analytical solution available)
+#
+# **Four Advection Schemes Tested:**
+# - **Upstream**: 1st-order accurate, stable but diffusive
+# - **Central Difference (CD)**: 2nd-order accurate but can oscillate on discontinuities
+# - **TVD** (Total Variation Diminishing): Handles sharp fronts well, works reliably only on structured grids
+# - **UTVD** (Unstructured TVD): TVD extended with unstructured grid support, maintains TVD-quality performance on all grid types
+#
+# **Three Test Functions (probing different numerical challenges):**
+# - **sin² wave**: Smooth function testing 2nd-order accuracy
+# - **Square wave**: Sharp discontinuity testing stability
+# - **Step wave**: Sharp transition testing boundedness
+#
+# **Three Grid Types:**
+# - **Structured**: Regular rectangular cells
+# - **Triangle**: Triangular mesh elements
+# - **Voronoi**: Voronoi polygon cells
+#
+# **Expected Results:**
+# - CD scheme should oscillate/fail on discontinuous functions
+# - TVD should work well on structured grids but may have issues on unstructured grids
+# - UTVD should handle discontinuities without oscillation across all grid types
+# - Different grid geometries may show different accuracy characteristics for the same numerical scheme
+
+# %% [markdown]
+# ### Initial setup
+#
+# Import dependencies, define the example name and workspace, and read settings from environment variables.
+
+# %%
+import collections.abc
+import itertools
+import math
+from pathlib import Path
+
+import flopy
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+from flopy.discretization.vertexgrid import VertexGrid
+from flopy.plot.styles import styles
+from flopy.utils import GridIntersect
+from flopy.utils.triangle import Triangle
+from flopy.utils.voronoi import VoronoiGrid
+from modflow_devtools.misc import get_env, timed
+from scipy.interpolate import LinearNDInterpolator
+from shapely.geometry import LineString
+
+try:
+    import git
+except ImportError:
+    git = None
+
+# Example name and workspace paths. If this example is running
+# in the git repository, use the folder structure described in
+# the README. Otherwise just use the current working directory.
+sim_name = "ex-gwt-adv-schemes"
+try:
+    if git is not None:
+        root = Path(git.Repo(".", search_parent_directories=True).working_dir)
+    else:
+        root = None
+except Exception:
+    root = None
+workspace = root / "examples" if root else Path.cwd()
+figs_path = root / "figures" if root else Path.cwd()
+
+# Settings from environment variables
+write = get_env("WRITE", True)
+run = get_env("RUN", True)
+plot = get_env("PLOT", True)
+plot_show = get_env("PLOT_SHOW", True)
+plot_save = get_env("PLOT_SAVE", True)
+
+# %% [markdown]
+# ## Define parameters
+#
+# Define model units, parameters and other settings.
+
+# %%
+# Model units
+length_units = "centimeters"  # Length units
+time_units = "seconds"  # Time units
+
+# Model parameters
+nper = 1  # Number of periods
+nlay = 1  # Number of layers for the structured grid
+ncol = 50  # Number of columns for the structured grid
+nrow = 50  # Number of rows for the structured grid
+Length = 100.0  # Domain length (cm)
+Width = 100.0  # Domain width (cm)
+
+delr = 2  # Column width (cm) for the structured grid
+delc = 2  # Row width (cm) for the structured grid
+top = 1.0  # Top elevation of the model (cm)
+botm = 0  # Bottom elevation of the model (cm)
+
+specific_discharge = 0.5  # Specific discharge (cm/s)
+hydraulic_conductivity = 0.01  # Hydraulic conductivity (cm/s)
+angledeg = 45  # Flow direction (°)
+profile_width = 20.0  # Width of the inflow concentration profiles (cm)
+
+total_time = 300.0  # Total simulation time (s)
+dt = 5  # Initial time step (s)
+percel = 0.7  # Courant number (-)
+
+# Solver parameters
+nouter = 100  # Maximum outer iterations
+ninner = 300  # Maximum inner iterations
+hclose = 1e-6  # Head closure criterion
+cclose = 1e-6  # Concentration closure criterion
+rclose = 1e-6  # Residual closure criterion
+
+# Grid and scheme definitions
+grids = ["structured", "triangle", "voronoi"]  # 3 grid types (36 total simulations)
+schemes = ["upstream", "central", "tvd", "utvd"]  # 4 advection schemes
+wave_functions = ["sin²-wave", "step-wave", "square-wave"]  # 3 test functions
+
+# Compute discharge components
+angle = math.radians(angledeg)
+qx = specific_discharge * math.cos(angle)  # x-component of specific discharge (cm/s)
+qy = specific_discharge * math.cos(angle)  # y-component of specific discharge (cm/s)
+
+AXES_FRACTION = "axes fraction"
+OFFSET_POINTS = "offset points"
+
+# %% [markdown]
+# # Analytical Solution
+#
+# Exact solution of the concentration field for pure advection.
+#
+# The analytical solution works by:
+# 1. Rotating coordinates to align with flow direction
+# 2. Applying the 1D inlet signal in the rotated coordinate system
+# 3. This works because advection simply translates the inlet pattern
+#
+# For uniform flow at angle θ, any point (x,y) maps to:
+# - Rotated coordinate: y' = sin(-θ)·x + cos(-θ)·y
+# - Concentration: C(x,y,t) = inlet_signal(y' - v·t)
+# Note: At t=300s, the pattern has fully advected through the domain
+
+
+# %%
+def exact_solution_concentration(x, y, analytical):
+    """Calculate exact concentration at any point in the domain.
+
+    The analytical solution works by:
+    1. Rotating coordinates to align with flow direction
+    2. Applying the 1D inlet signal in the rotated coordinate system
+    3. This works because advection simply translates the inlet pattern
+
+    Args:
+        x, y: Spatial coordinates (cm)
+        analytical: Signal type ('sin²-wave', 'step-wave', 'square-wave')
+
+    Returns:
+        Concentration values [dimensionless, 0-1]
+    """
+    # Rotate to 1d solution space
+    reverse_angle = -angle
+    rotated_y = math.sin(reverse_angle) * x + math.cos(reverse_angle) * y
+
+    # Compute concentration
+    return inlet_signal(rotated_y, analytical)
+
+
+def inlet_signal(y, signal_name):
+    """Generate inlet signal based on the signal type.
+
+    Args:
+        y: y-coordinate values
+        signal_name: Type of signal ('step-wave', 'square-wave', 'sin²-wave')
+
+    Returns:
+        Concentration values based on the signal type
+    """
+    clipped_y = np.clip(y, -profile_width / 2, profile_width / 2)
+    match signal_name:
+        case "step-wave":
+            return np.where(y < 0, np.ones(y.shape), np.zeros(y.shape))
+        case "square-wave":
+            return np.where(
+                np.abs(y) < profile_width / 2, np.ones(y.shape), np.zeros(y.shape)
+            )
+        case "sin²-wave":
+            return np.power(np.cos(math.pi * clipped_y / profile_width), 2)
+        case _:
+            raise ValueError("Unknown signal name")
+
+
+# %% [markdown]
+# # Grid Helper methods
+
+
+# %%
+def grid_triangulator(itri, delr, delc):
+    nrow, ncol = itri.shape
+    if np.isscalar(delr):
+        delr = delr * np.ones(ncol)
+    if np.isscalar(delc):
+        delc = delc * np.ones(nrow)
+    regular_grid = flopy.discretization.StructuredGrid(delc, delr)
+    vertdict = {}
+    icell = 0
+    for i in range(nrow):
+        for j in range(ncol):
+            vs = regular_grid.get_cell_vertices(i, j)
+            if itri[i, j] == 0:
+                vertdict[icell] = [vs[0], vs[1], vs[2], vs[3], vs[0]]
+                icell += 1
+            elif itri[i, j] == 1:
+                vertdict[icell] = [vs[0], vs[1], vs[3], vs[0]]
+                icell += 1
+                vertdict[icell] = [vs[3], vs[1], vs[2], vs[3]]
+                icell += 1
+            elif itri[i, j] == 2:
+                vertdict[icell] = [vs[0], vs[2], vs[3], vs[0]]
+                icell += 1
+                vertdict[icell] = [vs[0], vs[1], vs[2], vs[0]]
+                icell += 1
+            else:
+                raise ValueError(f"Unknown itri value: {itri[i, j]}")
+    verts, iverts = flopy.utils.cvfdutil.to_cvfd(vertdict)
+    return verts, iverts
+
+
+def cvfd_to_cell2d(verts, iverts):
+    vertices = []
+    for i in range(verts.shape[0]):
+        x = verts[i, 0]
+        y = verts[i, 1]
+        vertices.append([i, x, y])
+    cell2d = []
+    for icell2d, vs in enumerate(iverts):
+        points = [tuple(verts[ip]) for ip in vs]
+        xc, yc = flopy.utils.cvfdutil.centroid_of_polygon(points)
+        cell2d.append([icell2d, xc, yc, len(vs), *vs])
+    return vertices, cell2d
+
+
+def grid_intersector(vgrid):
+    """Create a grid intersector for the given vertex grid."""
+    return flopy.utils.GridIntersect(vgrid)
+
+
+def get_boundary(gi, line):
+    line = LineString(line)
+    cells = gi.intersect(line)["cellids"]
+    cells = np.array(list(cells))
+
+    return cells
+
+
+def create_bc(boundary_ids, value, auxvalue=None):
+    if auxvalue is None:
+        if isinstance(value, collections.abc.Sequence | np.ndarray):
+            return [
+                [(0, cell_id), value[idx]] for idx, cell_id in enumerate(boundary_ids)
+            ]
+        else:
+            return [[(0, cell_id), value] for idx, cell_id in enumerate(boundary_ids)]
+    else:
+        if isinstance(auxvalue, collections.abc.Sequence | np.ndarray):
+            return [
+                [(0, cell_id), value, auxvalue[idx]]
+                for idx, cell_id in enumerate(boundary_ids)
+            ]
+        else:
+            return [
+                [(0, cell_id), value, auxvalue]
+                for idx, cell_id in enumerate(boundary_ids)
+            ]
+
+
+def merge_bc_sum(boundaries):
+    df = pd.DataFrame(boundaries)
+    summed = df.groupby([0], as_index=False)[1].sum()
+
+    return summed.values.tolist()
+
+
+def merge_bc_mean(boundaries):
+    df = pd.DataFrame(boundaries)
+    mean = df.groupby([0], as_index=False)[1].mean()
+
+    return mean.values.tolist()
+
+
+def create_grid(grid_type):
+    """Create grid based on the specified type.
+
+    Args:
+        grid_type: Type of grid ('structured', 'triangle', 'voronoi')
+
+    Returns:
+        Tuple of (ncpl, nvert, vertices, cell2d)
+    """
+    if grid_type == "structured":
+        itri = np.zeros((nrow, ncol), dtype=int)
+        verts, iverts = grid_triangulator(itri, delr, delc)
+        vertices, cell2d = cvfd_to_cell2d(verts, iverts)
+
+        ncpl = len(cell2d)
+        nvert = len(verts)
+
+        return ncpl, nvert, vertices, cell2d
+    elif grid_type == "triangle":
+        active_domain = [(0, 0), (Length, 0), (Length, Width), (0, Width)]
+        tri = Triangle(angle=30, maximum_area=delc * delr * 1.5, model_ws=workspace)
+        tri.add_polygon(active_domain)
+        tri.build()
+
+        cell2d = tri.get_cell2d()
+        vertices = tri.get_vertices()
+        ncpl = tri.ncpl
+        nvert = tri.nvert
+
+        return ncpl, nvert, vertices, cell2d
+    elif grid_type == "voronoi":
+        active_domain = [(0, 0), (Length, 0), (Length, Width), (0, Width)]
+        tri = Triangle(
+            angle=30, maximum_area=delc * delr / 1.5 * 1.2, model_ws=workspace
+        )
+        tri.add_polygon(active_domain)
+        tri.build()
+
+        vor = VoronoiGrid(tri)
+        disv_gridprops = vor.get_gridprops_vertexgrid()
+
+        cell2d = disv_gridprops["cell2d"]
+        vertices = disv_gridprops["vertices"]
+        ncpl = disv_gridprops["ncpl"]
+        nvert = len(vertices)
+
+        return ncpl, nvert, vertices, cell2d
+
+    else:
+        raise ValueError(f"grid of type '{grid_type}' is not supported.")
+
+
+def axis_aligned_segment_length(polygon, axis="y", value=0):
+    """Calculate the total length of segments aligned with the specified axis at the given value."""
+    total_length = 0.0
+    coords = list(polygon.exterior.coords)
+
+    for i in range(len(coords) - 1):
+        p1, p2 = coords[i], coords[i + 1]
+        is_aligned = (axis == "y" and p1[1] == value and p2[1] == value) or (
+            axis == "x" and p1[0] == value and p2[0] == value
+        )
+        if is_aligned:
+            segment = LineString([p1, p2])
+            total_length += segment.length
+
+    return total_length
+
+
+# %% [markdown]
+# # Model setup
+#
+# Define functions to build models, write input files, and run the simulation.
+
+
+# %%
+def build_mf6gwf(grid_type):
+    gwfname = f"flow_{grid_type}"
+    sim_ws = workspace / sim_name / Path(gwfname)
+    sim = flopy.mf6.MFSimulation(sim_name=sim_name, sim_ws=sim_ws, exe_name="mf6")
+
+    tdis_ds = ((total_time, 1, 1.0),)
+    flopy.mf6.ModflowTdis(sim, nper=nper, perioddata=tdis_ds, time_units=time_units)
+
+    flopy.mf6.ModflowIms(
+        sim,
+        print_option="ALL",
+        linear_acceleration="bicgstab",
+        outer_maximum=nouter,
+        outer_dvclose=hclose,
+        inner_maximum=ninner,
+        inner_dvclose=hclose,
+        rcloserecord=rclose,
+    )
+
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=gwfname, save_flows=True)
+
+    flopy.mf6.ModflowGwfnpf(
+        gwf,
+        save_specific_discharge=True,
+        save_saturation=True,
+        export_array_ascii=True,
+        xt3doptions=True,
+        icelltype=0,
+        k=hydraulic_conductivity,
+    )
+
+    flopy.mf6.ModflowGwfic(gwf, strt=0.0, export_array_ascii=True)
+
+    ncpl, nvert, vertices, cell2d = create_grid(grid_type)
+    flopy.mf6.ModflowGwfdisv(
+        gwf,
+        nlay=nlay,
+        ncpl=ncpl,
+        nvert=nvert,
+        top=top,
+        botm=botm,
+        vertices=vertices,
+        cell2d=cell2d,
+        filename=f"{gwfname}.disv",
+    )
+
+    head_filerecord = f"{gwfname}.hds"
+    budget_filerecord_csv = f"{gwfname}.bud.csv"
+    budget_filerecord = f"{gwfname}.bud"
+    flopy.mf6.ModflowGwfoc(
+        gwf,
+        head_filerecord=head_filerecord,
+        budgetcsv_filerecord=budget_filerecord_csv,
+        budget_filerecord=budget_filerecord,
+        saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+    )
+
+    # Boundary lines
+    bottom_edge = [(0, 0), (Length, 0)]
+    right_edge = [(Length, 0), (Length, Width)]
+    top_edge = [(Length, Width), (0, Width)]
+    left_edge = [(0, Width), (0, 0)]
+
+    # Identify boundary ids
+    vgrid = VertexGrid(vertices=vertices, cell2d=cell2d, ncpl=ncpl, nlay=1)
+    gi = grid_intersector(vgrid)
+
+    bottom_boundary_ids = get_boundary(gi, bottom_edge)
+    right_boundary_ids = get_boundary(gi, right_edge)
+    top_boundary_ids = get_boundary(gi, top_edge)
+    left_boundary_ids = get_boundary(gi, left_edge)
+
+    # Set top right element to a chd 0. This makes the solution unique
+    geometry = vgrid.geo_dataframe.geometry
+
+    topright_id, _, _ = np.intersect1d(
+        right_boundary_ids, top_boundary_ids, return_indices=True
+    )
+
+    if not topright_id.any():
+        topright_id = [
+            geometry.loc[top_boundary_ids].get_coordinates()["x"].idxmax().tolist()
+        ]
+
+    flopy.mf6.ModflowGwfchd(gwf, stress_period_data=create_bc(topright_id, 0))
+
+    # Set boundary conditions
+    cell_area = geometry.area.values
+    inflow_area_left = (
+        geometry.loc[left_boundary_ids]
+        .apply(lambda poly: axis_aligned_segment_length(poly, axis="x", value=0))
+        .values
+    )
+    inflow_area_bot = (
+        geometry.loc[bottom_boundary_ids]
+        .apply(lambda poly: axis_aligned_segment_length(poly, axis="y", value=0))
+        .values
+    )
+    outflow_area_right = (
+        geometry.loc[right_boundary_ids]
+        .apply(lambda poly: axis_aligned_segment_length(poly, axis="x", value=Length))
+        .values
+    )
+    outflow_area_top = (
+        geometry.loc[top_boundary_ids]
+        .apply(lambda poly: axis_aligned_segment_length(poly, axis="y", value=Width))
+        .values
+    )
+
+    flopy.mf6.ModflowGwfrch(
+        gwf,
+        stress_period_data=merge_bc_sum(
+            create_bc(
+                left_boundary_ids, qx * inflow_area_left / cell_area[left_boundary_ids]
+            )
+            + create_bc(
+                bottom_boundary_ids,
+                qy * inflow_area_bot / cell_area[bottom_boundary_ids],
+            )
+            + create_bc(
+                right_boundary_ids,
+                -qx * outflow_area_right / cell_area[right_boundary_ids],
+            )
+            + create_bc(
+                top_boundary_ids, -qy * outflow_area_top / cell_area[top_boundary_ids]
+            )
+        ),
+    )
+
+    return sim
+
+
+def build_mf6gwt(grid_type, scheme, wave_func):
+    pathname = f"trans_{grid_type}_{wave_func}_{scheme}"
+    gwtname = "trans"
+    gwfname = f"flow_{grid_type}"
+    sim_ws = workspace / sim_name / Path(pathname)
+    sim = flopy.mf6.MFSimulation(sim_name=sim_name, sim_ws=sim_ws, exe_name="mf6")
+
+    nsteps = int(total_time / dt)
+    tdis_ds = ((total_time, nsteps, 1.0),)
+
+    tdis = flopy.mf6.ModflowTdis(
+        sim, nper=nper, perioddata=tdis_ds, time_units=time_units
+    )
+
+    dtmin = 1e-5
+    dtmax = 10
+    dtadj = 2
+    dtfailadj = dtadj
+    ats_filerecord = gwtname + ".ats"
+    atsperiod = [(0, dt, dtmin, dtmax, dtadj, dtfailadj)]
+    tdis.ats.initialize(
+        maxats=len(atsperiod),
+        perioddata=atsperiod,
+        filename=ats_filerecord,
+    )
+
+    flopy.mf6.ModflowIms(
+        sim,
+        linear_acceleration="bicgstab",
+        print_option="SUMMARY",
+        outer_maximum=nouter,
+        outer_dvclose=cclose,
+        inner_maximum=ninner,
+        inner_dvclose=cclose,
+        rcloserecord=rclose,
+    )
+
+    gwt = flopy.mf6.ModflowGwt(sim, modelname=gwtname, save_flows=True)
+
+    flopy.mf6.ModflowGwtmst(gwt, porosity=1.0)
+
+    flopy.mf6.ModflowGwtssm(gwt)
+
+    flopy.mf6.ModflowGwtadv(gwt, scheme=scheme, ats_percel=percel)
+
+    packagedata = [
+        ("GWFHEAD", f"../{gwfname}/{gwfname}.hds", None),
+        ("GWFBUDGET", f"../{gwfname}/{gwfname}.bud", None),
+    ]
+    flopy.mf6.ModflowGwtfmi(gwt, packagedata=packagedata, save_flows=True)
+
+    flopy.mf6.ModflowGwtoc(
+        gwt,
+        budget_filerecord=f"{gwtname}.cbc",
+        concentration_filerecord=f"{gwtname}.ucn",
+        saverecord=[("CONCENTRATION", "LAST"), ("BUDGET", "LAST")],
+        printrecord=[("CONCENTRATION", "LAST"), ("BUDGET", "LAST")],
+    )
+
+    ncpl, nvert, vertices, cell2d = create_grid(grid_type)
+    flopy.mf6.ModflowGwtdisv(
+        gwt,
+        nlay=nlay,
+        ncpl=ncpl,
+        nvert=nvert,
+        top=top,
+        botm=botm,
+        vertices=vertices,
+        cell2d=cell2d,
+        filename=f"{gwtname}.disv",
+    )
+
+    # Boundary lines
+    bottom_edge = [(0, 0), (Length, 0)]
+    left_edge = [(0, Width), (0, 0)]
+
+    # Identify boundary ids
+    vgrid = VertexGrid(vertices=vertices, cell2d=cell2d, ncpl=ncpl, nlay=1)
+    gi = grid_intersector(vgrid)
+
+    bottom_boundary_ids = get_boundary(gi, bottom_edge)
+    left_boundary_ids = get_boundary(gi, left_edge)
+
+    # Compute exact solution for the inlet boundaries
+    cell2d_df = pd.DataFrame(cell2d)
+    xc, yc = cell2d_df.loc[bottom_boundary_ids, 1:2].T.to_numpy()
+    conc_bottom = exact_solution_concentration(xc, yc, wave_func)
+    xc, yc = cell2d_df.loc[left_boundary_ids, 1:2].T.to_numpy()
+    conc_left = exact_solution_concentration(xc, yc, wave_func)
+
+    # Set inlet concentrations
+    flopy.mf6.ModflowGwtcnc(
+        gwt,
+        stress_period_data=merge_bc_mean(
+            create_bc(bottom_boundary_ids, conc_bottom)
+            + create_bc(left_boundary_ids, conc_left)
+        ),
+    )
+
+    # Set initial condition equal to zero
+    flopy.mf6.ModflowGwtic(gwt, strt=0)
+
+    return sim
+
+
+# %%
+def write_models(sim, silent=True):
+    sim.write_simulation(silent=silent)
+
+
+@timed
+def run_models(sim, silent=True):
+    success, buff = sim.run_simulation(silent=silent)
+    assert success, buff
+
+
+# %% [markdown]
+# # Plotting results
+#
+# Define functions to plot model results.
+
+
+# %%
+def plot_results(gwf_sims, gwt_sims):
+    plot_flows(gwf_sims)
+    plot_concentrations(gwt_sims)
+    plot_concentration_cross_sections(gwt_sims)
+
+
+def plot_flows(gwf_sims):
+    with styles.USGSPlot():
+        figure_size = (7, 7 / len(gwf_sims))
+        fig, axs = plt.subplots(
+            1, len(gwf_sims), dpi=300, figsize=figure_size, tight_layout=True
+        )
+        fig.suptitle("Head - flow angle 45")
+
+        for idx, (grid, sim) in enumerate(gwf_sims.items()):
+            plot_flow(sim, axs[idx])
+
+        pad = 5  # in points
+        for ax, col in zip(axs, gwf_sims.keys()):
+            ax.annotate(
+                col,
+                xy=(0.5, 1),
+                xytext=(0, pad),
+                xycoords=AXES_FRACTION,
+                textcoords=OFFSET_POINTS,
+                size="large",
+                ha="center",
+                va="baseline",
+            )
+
+    if plot_show:
+        plt.show()
+    if plot_save:
+        fname = f"{sim_name}-flow.png"
+        fpth = figs_path / fname
+        fig.savefig(fpth)
+
+
+def plot_flow(sim, ax):
+    gwf = sim.get_model()
+    head = gwf.output.head().get_data()
+    bud = gwf.output.budget()
+    spdis = bud.get_data(text="DATA-SPDIS")[0]
+    vmin = head.min()
+    vmax = head.max()
+
+    pmv = flopy.plot.PlotMapView(gwf, ax=ax)
+    pmv.plot_grid(colors="k", alpha=0.1)
+    pc = pmv.plot_array(head, vmin=vmin, vmax=vmax, alpha=0.5)
+    pmv.contour_array(head, levels=np.linspace(vmin, vmax, 10))
+    plt.colorbar(pc)
+
+    ax.set_xlabel("x [cm]")
+    ax.set_ylabel("y [cm]")
+    ax.set_aspect("equal")
+
+
+def plot_concentrations(gwt_sims):
+    for wave_func in wave_functions:
+        with styles.USGSPlot():
+            fig, axs = plt.subplots(
+                nrows=len(schemes),
+                ncols=len(grids),
+                dpi=300,
+                figsize=(7, 7 * len(schemes) / (len(grids) + 1)),
+                tight_layout=True,
+            )
+            fig.suptitle(f"Concentration - {wave_func}")
+
+            for idx_scheme, scheme in enumerate(schemes):
+                for idx_grid, grid in enumerate(grids):
+                    sim = gwt_sims[(grid, scheme, wave_func)]
+                    plot_concentration(sim, axs[idx_scheme, idx_grid])
+
+            pad = 5  # in points
+            for ax, col in zip(axs[0], grids):
+                ax.annotate(
+                    col,
+                    xy=(0.5, 1),
+                    xytext=(0, pad),
+                    xycoords=AXES_FRACTION,
+                    textcoords=OFFSET_POINTS,
+                    size="large",
+                    ha="center",
+                    va="baseline",
+                )
+
+            for ax, row in zip(axs[:, 0], schemes):
+                ax.annotate(
+                    row,
+                    xy=(0, 0.5),
+                    xytext=(-ax.yaxis.labelpad - pad, 0),
+                    xycoords=ax.yaxis.label,
+                    textcoords=OFFSET_POINTS,
+                    size="large",
+                    ha="right",
+                    va="center",
+                )
+
+            fig.subplots_adjust(left=0.25, top=0.95)
+
+            if plot_show:
+                plt.show()
+            if plot_save:
+                fname = f"{sim_name}-{wave_func}-conc.png"
+                fpth = figs_path / fname
+                fig.savefig(fpth)
+
+
+def plot_concentration(sim, ax):
+    gwt = sim.get_model()
+    ucnobj_mf6 = gwt.output.concentration()
+    conc = ucnobj_mf6.get_data(totim=ucnobj_mf6.get_times()[-1]).flatten()
+
+    tol = 1e-2
+    vmin = -0.0 - tol
+    vmax = 1.0 + tol
+
+    masked_conc = np.ma.masked_outside(conc, vmin, vmax)
+    pmv = flopy.plot.PlotMapView(gwt, ax=ax)
+    pc = pmv.plot_array(masked_conc, vmin=vmin, vmax=vmax, alpha=1)
+    plt.colorbar(pc)
+
+    ax.set_xlabel("x [cm]")
+    ax.set_ylabel("y [cm]")
+    ax.set_aspect("equal")
+
+
+def plot_concentration_cross_sections(gwt_sims):
+    with styles.USGSPlot():
+        fig, axs = plt.subplots(
+            nrows=len(wave_functions),
+            ncols=len(grids),
+            dpi=300,
+            figsize=(7, 7 * len(wave_functions) / (len(grids) + 1)),
+            tight_layout=True,
+        )
+        fig.suptitle("Concentration cross-section")
+
+        for wave_idx, wave_func in enumerate(wave_functions):
+            for grid_idx, grid in enumerate(grids):
+                ax = axs[wave_idx][grid_idx]
+                plot_concentration_analytical(wave_func, ax)
+                for scheme in schemes:
+                    sim = gwt_sims[(grid, scheme, wave_func)]
+                    plot_concentration_cross_section(sim, scheme, ax)
+
+                ax.legend(fontsize="small")
+                ax.set_xlabel("x [cm]")
+                ax.set_ylabel("C [g/cm³]")
+                ax.set_ylim(-0.1, 1.1)
+
+        pad = 5  # in points
+        for ax, col in zip(axs[0], grids):
+            ax.annotate(
+                col,
+                xy=(0.5, 1),
+                xytext=(0, pad),
+                xycoords=AXES_FRACTION,
+                textcoords=OFFSET_POINTS,
+                size="large",
+                ha="center",
+                va="baseline",
+            )
+
+        for ax, row in zip(axs[:, 0], wave_functions):
+            ax.annotate(
+                row,
+                xy=(0, 0.5),
+                xytext=(-ax.yaxis.labelpad - pad, 0),
+                xycoords=ax.yaxis.label,
+                textcoords=OFFSET_POINTS,
+                size="large",
+                ha="right",
+                va="center",
+            )
+
+        fig.subplots_adjust(left=0.25, top=0.95)
+
+    if plot_show:
+        plt.show()
+    if plot_save:
+        fname = f"{sim_name}-conc-cross-section.png"
+        fpth = figs_path / fname
+        fig.savefig(fpth)
+
+
+def plot_concentration_analytical(analytical_func, ax):
+    x = np.linspace(0, Length, 100)
+    y = Width / 2
+
+    conc = exact_solution_concentration(x, y, analytical_func)
+
+    ax.plot(
+        x, conc, linestyle="-", mfc="none", markersize="3", linewidth=1, label="exact"
+    )
+
+
+def plot_concentration_cross_section(sim, scheme, ax):
+    gwt = sim.get_model()
+    ucnobj_mf6 = gwt.output.concentration()
+    conc = ucnobj_mf6.get_data(totim=ucnobj_mf6.get_times()[-1]).flatten()
+
+    grid = gwt.modelgrid
+    xc = grid.xcellcenters
+    yc = grid.ycellcenters
+
+    ix = GridIntersect(grid)
+    x_line = LineString([(0, Width / 2), (Length, Width / 2)])
+    xi = ix.intersect(x_line)
+    x_along_line = sorted(xc[xi.cellids.tolist()])
+
+    interp = LinearNDInterpolator(list(zip(xc, yc)), conc)
+    conc_interp = interp([(i, Width / 2) for i in x_along_line])
+
+    ax.plot(
+        x_along_line,
+        conc_interp,
+        linestyle="--",
+        marker="o",
+        mfc="none",
+        markersize="3",
+        linewidth="0.5",
+        markeredgewidth="0.5",
+        label=scheme,
+    )
+
+
+# %% [markdown]
+# # Running the example
+#
+# Define and invoke a function to run the example scenario, then plot results.
+
+
+# %%
+def scenario(silent=True):
+    print("=" * 60)
+    print("MODFLOW 6 Advection Schemes Comparison")
+    print("=" * 60)
+    print(
+        f"Total simulations: {len(grids)} grids x {len(schemes)} schemes x {len(wave_functions)} functions = {len(grids) * len(schemes) * len(wave_functions)} transport models"
+    )
+    print(f"Plus {len(grids)} flow models")
+    print()
+
+    gwf_sims = {}
+
+    # Build and run the gwt models on different grids
+    for grid in grids:
+        gwf_sims[grid] = build_mf6gwf(grid)
+
+    # Build the gwt models
+    gwt_sims = {}
+    combinations = list(itertools.product(*[grids, schemes, wave_functions]))
+    for combination in combinations:
+        gwt_sims[combination] = build_mf6gwt(
+            combination[0], combination[1], combination[2]
+        )
+
+    if write:
+        print("\nWriting flow models:")
+        for grid, sim in gwf_sims.items():
+            print(f"- {grid} grid")
+            write_models(sim, silent=silent)
+
+        print("\nWriting transport models:")
+        for key, sim in gwt_sims.items():
+            grid, scheme, wave = key
+            print(f"- {grid} grid \t {scheme} scheme \t {wave} function")
+            write_models(sim, silent=silent)
+
+    if run:
+        print("\nRunning flow models:")
+        for grid, sim in gwf_sims.items():
+            print(f"- {grid} grid")
+            run_models(sim, silent=silent)
+
+        print("\nRunning transport models:")
+        for key, sim in gwt_sims.items():
+            grid, scheme, wave = key
+            print(f"- {grid} grid \t {scheme} scheme \t {wave} function")
+            run_models(sim, silent=silent)
+    if plot:
+        print("\n" + "=" * 40)
+        print("PLOTTING RESULTS")
+        print("=" * 40)
+        print("Generating 3 figure sets:")
+        print("1. Flow fields (3 subplots)")
+        print(
+            f"2. Concentration maps ({len(schemes)} rows x {len(grids)} cols per wave func = {len(schemes) * len(grids) * len(wave_functions)} total subplots)"
+        )
+        print(
+            f"3. Cross-section comparisons ({len(wave_functions)} rows x {len(grids)} cols = {len(wave_functions) * len(grids)} subplots)"
+        )
+        plot_results(gwf_sims, gwt_sims)
+
+
+scenario()