diff --git a/flopy/utils/crosssection.py b/flopy/utils/crosssection.py
new file mode 100644
index 0000000000..725c1e250d
--- /dev/null
+++ b/flopy/utils/crosssection.py
@@ -0,0 +1,276 @@
+import numpy as np
+
+# power for Manning's hydraulic radius term
+_mpow = 2.0 / 3.0
+
+
+def calculate_rectchan_mannings_discharge(
+    conversion_factor, roughness, slope, width, depth
+):
+    """
+    Calculate Manning's discharge for a rectangular channel.
+
+    """
+    area = width * depth
+    return conversion_factor * area * depth**_mpow * slope**0.5 / roughness
+
+
+# n-point cross-section functions
+def get_wetted_station(
+    x0,
+    x1,
+    h0,
+    h1,
+    depth,
+):
+    """Get the wetted length in the x-direction"""
+    # -- calculate the minimum and maximum depth
+    hmin = min(h0, h1)
+    hmax = max(h0, h1)
+
+    # -- if depth is less than or equal to the minimum value the
+    #    station length (xlen) is zero
+    if depth <= hmin:
+        x1 = x0
+    # -- if depth is between hmin and hmax, station length is less
+    #    than h1 - h0
+    elif depth < hmax:
+        xlen = x1 - x0
+        dlen = h1 - h0
+        if abs(dlen) > 0.0:
+            slope = xlen / dlen
+        else:
+            slope = 0.0
+        if h0 > h1:
+            dx = (depth - h1) * slope
+            xt = x1 + dx
+            xt0 = xt
+            xt1 = x1
+        else:
+            dx = (depth - h0) * slope
+            xt = x0 + dx
+            xt0 = x0
+            xt1 = xt
+        x0 = xt0
+        x1 = xt1
+    return x0, x1
+
+
+def get_wetted_perimeter(
+    x0,
+    x1,
+    h0,
+    h1,
+    depth,
+):
+    # -- calculate the minimum and maximum depth
+    hmin = min(h0, h1)
+    hmax = max(h0, h1)
+
+    # -- calculate the wetted perimeter for the segment
+    xlen = x1 - x0
+    if xlen > 0.0:
+        if depth > hmax:
+            dlen = hmax - hmin
+        else:
+            dlen = depth - hmin
+    else:
+        if depth > hmin:
+            dlen = min(depth, hmax) - hmin
+        else:
+            dlen = 0.0
+    return np.sqrt(xlen**2.0 + dlen**2.0)
+
+
+def get_wetted_area(x0, x1, h0, h1, depth):
+    # -- calculate the minimum and maximum depth
+    hmin = min(h0, h1)
+    hmax = max(h0, h1)
+
+    # -- calculate the wetted area for the segment
+    xlen = x1 - x0
+    area = 0.0
+    if xlen > 0.0:
+        # -- add the area above hmax
+        if depth > hmax:
+            area = xlen * (depth - hmax)
+        # -- add the area below zmax
+        if hmax != hmin and depth > hmin:
+            area += 0.5 * (depth - hmin)
+    return area
+
+
+def wetted_area(
+    x,
+    h,
+    depth,
+    verbose=False,
+):
+    area = 0.0
+    if x.shape[0] == 1:
+        area = x[0] * depth
+    else:
+        for idx in range(0, x.shape[0] - 1):
+            x0, x1 = x[idx], x[idx + 1]
+            h0, h1 = h[idx], h[idx + 1]
+
+            # get station data
+            x0, x1 = get_wetted_station(x0, x1, h0, h1, depth)
+
+            # get wetted area
+            a = get_wetted_area(x0, x1, h0, h1, depth)
+            area += a
+
+            # write to screen
+            if verbose:
+                print(
+                    f"{idx}->{idx + 1} ({x0},{x1}) - "
+                    f"perimeter={x1 - x0} - area={a}"
+                )
+
+    return area
+
+
+def wetted_perimeter(
+    x,
+    h,
+    depth,
+    verbose=False,
+):
+    perimeter = 0.0
+    if x.shape[0] == 1:
+        perimeter = x[0]
+    else:
+        for idx in range(0, x.shape[0] - 1):
+            x0, x1 = x[idx], x[idx + 1]
+            h0, h1 = h[idx], h[idx + 1]
+
+            # get station data
+            x0, x1 = get_wetted_station(x0, x1, h0, h1, depth)
+
+            # get wetted perimeter
+            perimeter += get_wetted_perimeter(x0, x1, h0, h1, depth)
+
+            # write to screen
+            if verbose:
+                print(f"{idx}->{idx + 1} ({x0},{x1}) - perimeter={x1 - x0}")
+
+    return perimeter
+
+
+def manningsq(
+    x,
+    h,
+    depth,
+    roughness=0.01,
+    slope=0.001,
+    conv=1.0,
+):
+    if isinstance(roughness, float):
+        roughness = np.ones(x.shape, dtype=float) * roughness
+    if x.shape[0] > 1:
+        q = 0.0
+        for i0 in range(x.shape[0] - 1):
+            i1 = i0 + 1
+            perimeter = get_wetted_perimeter(x[i0], x[i1], h[i0], h[i1], depth)
+            area = get_wetted_area(x[i0], x[i1], h[i0], h[i1], depth)
+            if perimeter > 0.0:
+                radius = area / perimeter
+                q += (
+                        conv * area * radius ** _mpow * slope ** 0.5 / roughness[i0]
+                )
+    else:
+        perimeter = wetted_perimeter(x, h, depth)
+        area = wetted_area(x, h, depth)
+        radius = 0.0
+        if perimeter > 0.0:
+            radius = area / perimeter
+        q = conv * area * radius ** _mpow * slope ** 0.5 / roughness[0]
+    return q
+
+
+def get_depths(
+    flows,
+    x,
+    h,
+    roughness=0.01,
+    slope=0.001,
+    conv=1.0,
+    dd=1e-4,
+    verbose=False,
+):
+    if isinstance(flows, float):
+        flows = np.array([flows], dtype=float)
+    if isinstance(roughness, float):
+        roughness = np.ones(x.shape, dtype=float) * roughness
+    depths = np.zeros(flows.shape, dtype=float)
+    for idx, q in enumerate(flows):
+        depths[idx] = qtodepth(
+            x,
+            h,
+            q,
+            roughness=roughness,
+            slope=slope,
+            conv=conv,
+            dd=dd,
+            verbose=False,
+        )
+
+    return depths
+
+
+def qtodepth(
+    x,
+    h,
+    q,
+    roughness=0.01,
+    slope=0.001,
+    conv=1.0,
+    dd=1e-4,
+    verbose=False,
+):
+    h0 = 0.0
+    q0 = manningsq(
+        x,
+        h,
+        h0,
+        roughness=roughness,
+        slope=slope,
+        conv=conv,
+    )
+    r = q0 - q
+
+    iter = 0
+    if verbose:
+        print(f"iteration {iter:>2d} - residual={r}")
+    while abs(r) > 1e-12:
+        q1 = manningsq(
+            x,
+            h,
+            h0 + dd,
+            roughness=roughness,
+            slope=slope,
+            conv=conv,
+        )
+        dq = q1 - q0
+        if dq != 0.0:
+            derv = dd / (q1 - q0)
+        else:
+            derv = 0.0
+        h0 -= derv * r
+        q0 = manningsq(
+            x,
+            h,
+            h0,
+            roughness=roughness,
+            slope=slope,
+            conv=conv,
+        )
+        r = q0 - q
+
+        iter += 1
+        if verbose:
+            print(f"iteration {iter:>2d} - residual={r}")
+        if iter > 100:
+            break
+    return h0