From 1afdd5c0c77422cd1380b23bb5087fda846cf5a7 Mon Sep 17 00:00:00 2001
From: w-bonelli <wesbonelli@gmail.com>
Date: Thu, 22 Dec 2022 15:19:19 -0500
Subject: [PATCH] feat(channel_utils): move open channel flow utils from
 modflow6

---
 autotest/test_channel_utils.py | 347 +++++++++++++++++++++
 flopy/utils/channel_utils.py   | 536 +++++++++++++++++++++++++++++++++
 2 files changed, 883 insertions(+)
 create mode 100644 autotest/test_channel_utils.py
 create mode 100644 flopy/utils/channel_utils.py

diff --git a/autotest/test_channel_utils.py b/autotest/test_channel_utils.py
new file mode 100644
index 0000000000..4db58dd3c3
--- /dev/null
+++ b/autotest/test_channel_utils.py
@@ -0,0 +1,347 @@
+import math
+
+import numpy as np
+import pytest
+
+from flopy.utils.channel_utils import (
+    get_depth,
+    get_discharge,
+    get_discharge_rect,
+    get_segment_wetted_area,
+    get_segment_wetted_perimeter,
+    get_segment_wetted_station,
+    get_wetted_area,
+    get_wetted_perimeter,
+)
+
+
+def test_get_segment_wetted_station_all_dry():
+    depth = 10
+    p0 = (0, 12)
+    p1 = (10, 11)
+    x0, x1 = get_segment_wetted_station(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+    # zero-length segment at x0
+    assert x0 == x1 == p0[0]
+
+
+def test_get_segment_wetted_station_partial():
+    depth = 10
+
+    # left bank (sloping downwards to the right)
+    p0 = (0, 12)
+    p1 = (10, 8)
+    x0, x1 = get_segment_wetted_station(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+    # left endpoint should be moved to the right
+    assert x0 != x1
+    assert x0 != p0[0]
+    assert x1 == p1[0]
+    assert x0 == (p1[0] - p0[0]) / 2
+
+    # right bank (sloping upwards to the right)
+    p0 = (0, 8)
+    p1 = (10, 12)
+    x0, x1 = get_segment_wetted_station(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+    # right endpoint should be moved to the left
+    assert x0 != x1
+    assert x0 == p0[0]
+    assert x1 != p1[0]
+    assert x1 == (p1[0] - p0[0]) / 2
+
+
+def test_get_segment_wetted_station_submerged():
+    depth = 13
+    p0 = (0, 12)
+    p1 = (10, 11)
+    x0, x1 = get_segment_wetted_station(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+    # entire segment should be returned
+    assert x0 == p0[0]
+    assert x1 == p1[0]
+
+
+def test_get_segment_wetted_perimeter_dry():
+    depth = 10
+    p0 = (0, 12)
+    p1 = (10, 11)
+    perim = get_segment_wetted_perimeter(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+    assert perim == 0
+
+
+point_lists = [
+    # single segments
+    [(0, -1), (1, -1)],
+    [(0, 0), (1, 1)],
+    [(0, 1), (2, -1)],
+    [(0, -1), (2, 1)],
+    [(0, 1), (1, 0)],
+    # channels with multiple segments
+    [(0, -1), (1, -1), (2, -1)],  # flat
+    [(0, 0), (1, -1), (2, 0)],  # triangular
+    [(0, -1), (1, -2), (2, -2), (3, -1)],  # trapezoidal
+    [(0, -1), (1, -2), (2, -4), (3, -4), (4, -1)],  # complex
+]
+
+
+@pytest.mark.parametrize(
+    "depth, p0, p1",
+    [
+        (0, point_lists[0][0], point_lists[0][1]),
+        (0, point_lists[1][0], point_lists[1][1]),
+        (0, point_lists[2][0], point_lists[2][1]),
+        (3, point_lists[3][0], point_lists[3][1]),
+    ],
+)
+def test_get_segment_wetted_perimeter(depth, p0, p1):
+    wp = get_segment_wetted_perimeter(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+
+    xlen = abs(p1[0] - p0[0])
+    hlen = abs(p1[1] - p0[1])
+    hmax = max([p0[1], p1[1]])
+    hmin = min([p0[1], p1[1]])
+    seg_len = math.sqrt(hlen**2 + xlen**2)
+
+    # if segment is fully wetted, wetted perimeter is just its length
+    if depth >= hmax:
+        # expect perimeter 0 if the water surface is level with a flat segment
+        if depth == hmin == hmax:
+            assert wp == 0
+        else:
+            assert wp == seg_len
+
+    # if segment is partially submerged, wetted perimeter should be
+    # less than the length of the segment but greater than zero
+    elif depth > hmin:
+        assert wp > 0
+        assert wp < seg_len
+
+    # if segment is completely dry, wetted perimeter should be zero
+    else:
+        assert wp == 0
+
+
+@pytest.mark.parametrize(
+    "depth, p0, p1",
+    [
+        (0, point_lists[0][0], point_lists[0][1]),
+        (0, point_lists[1][0], point_lists[1][1]),
+        (0, point_lists[2][0], point_lists[2][1]),
+        (3, point_lists[3][0], point_lists[3][1]),
+    ],
+)
+def test_get_segment_wetted_area(depth, p0, p1):
+    wa = get_segment_wetted_area(
+        x0=p0[0], x1=p1[0], h0=p0[1], h1=p1[1], depth=depth
+    )
+
+    xlen = abs(p1[0] - p0[0])
+    hlen = abs(p1[1] - p0[1])
+    hmax = max([p0[1], p1[1]])
+    hmin = min([p0[1], p1[1]])
+    seg_len = math.sqrt(hlen**2 + xlen**2)
+    tri_area = 0 if hlen == 0 else (0.5 * hlen * xlen)
+
+    # if segment is submerged wetted area is that of quadrilateral
+    # formed by the segment, the water surface, and vertical sides
+    if depth > hmax:
+        rect_area = xlen * (depth - hmax)
+        expected = rect_area + tri_area * (0.5 if p0[1] != p1[1] else 1)
+        assert wa == expected
+
+    # if segment is partially submerged, wetted area should be
+    # less than the area of the triangle formed by the segment
+    # and the water surface, with one vertical side
+    elif depth > hmin:
+        assert wa > 0
+        assert wa < tri_area
+
+    # if segment is completely dry, wetted area should be zero
+    else:
+        assert wa == 0
+
+
+def test_get_wetted_perimeter_rectangular():
+    depth = 0
+    points = np.array([[0, -0.2], [0, -1.4], [1.5, -1.4], [1.5, -0.2]])
+    perim = get_wetted_perimeter(points[:, 0], points[:, 1], depth)
+    assert perim == 1.5
+
+
+@pytest.mark.parametrize("points", [np.array(pts) for pts in point_lists])
+@pytest.mark.parametrize("depth", [0, 1, -1, -2])
+def test_get_wetted_perimeter(depth, points):
+    def total_perim(pts):
+        return sum(
+            [
+                math.sqrt(
+                    (pts[i][0] - pts[i - 1][0]) ** 2
+                    + (pts[i][1] - pts[i - 1][1]) ** 2
+                )
+                for i in range(1, len(pts))
+            ]
+        )
+
+    wp = get_wetted_perimeter(
+        x=points[:, 0], h=points[:, 1], depth=depth, verbose=True
+    )
+
+    hmax = max(points[:, 1])
+    hmin = min(points[:, 1])
+    total_perim = total_perim(points)
+
+    # if all segments are submerged, wetted perimeter is total perimeter
+    if depth >= hmax:
+        # expect perimeter 0 if the water surface is level with a flat channel
+        if all(p == depth for p in points[:, 1]):
+            assert wp == 0
+        else:
+            assert wp == total_perim
+
+    # if at least some segments are at least partially submerged...
+    elif depth > hmin:
+        assert wp > 0
+        assert wp < total_perim
+
+        def patch(x0, x1, h0, h1, depth):
+            # TODO: refactor get_segment_wetted_perimeter() to handle partial wetting
+            #  internally so separate get_segment_wetted_station() is no longer needed?
+
+            x0, x1 = get_segment_wetted_station(
+                x0=x0, x1=x1, h0=h0, h1=h1, depth=depth
+            )
+            return get_segment_wetted_perimeter(
+                x0=x0, x1=x1, h0=h0, h1=h1, depth=depth
+            )
+
+        assert np.isclose(
+            wp,
+            sum(
+                [
+                    patch(
+                        x0=points[i][0],
+                        x1=points[i + 1][0],
+                        h0=points[i][1],
+                        h1=points[i + 1][1],
+                        depth=depth,
+                    )
+                    for i in range(0, len(points) - 1)
+                ]
+            ),
+        )
+
+    # if all segments are dry, wetted perimeter is zero
+    else:
+        assert wp == 0
+
+
+def test_get_wetted_area_rectangular():
+    depth = 0
+    points = np.array([[0, -0.2], [0, -1.4], [1.5, -1.4], [1.5, -0.2]])
+    expected = 1.5 * 1.4
+    area = get_wetted_area(points[:, 0], points[:, 1], depth)
+    assert area == expected
+
+
+def test_get_discharge_rect():
+    # adapted from https://www.youtube.com/watch?v=R-Vhs_AH8mA
+
+    width = 0.5
+    depth = 0.25
+    n = 0.022
+    s = 0.005
+    k = 1.0
+    expected = 0.1
+    q = get_discharge_rect(width, depth, n, s, k)
+    assert np.isclose(expected, q, rtol=1e-2)
+
+
+def test_get_discharge_rectangular():
+    # adapted from https://www.youtube.com/watch?v=R-Vhs_AH8mA
+
+    # x and h as ints (width and height)
+    width = 0.5
+    depth = 0.25
+    n = 0.022
+    s = 0.005
+    k = 1.0
+    expected = 0.1
+    q = get_discharge(width, depth, depth, n, s, k)
+    assert np.isclose(expected, q, rtol=1e-2)
+
+    # x and h as arrays
+    x = np.array([0, 0, width, width])
+    h = np.array([depth, 0, 0, depth])
+    q = get_discharge(x, h, depth, n, s, k)
+    assert np.isclose(expected, q, rtol=1e-2)
+
+
+def test_get_discharge_trapezoidal():
+    # adapted from https://www.youtube.com/watch?v=ucLa9_DDWPA
+
+    x = np.array([0, 2, 4, 8, 10, 12])
+    h = np.array([2, 1, 0, 0, 1, 2])
+    n = 0.015
+    s = 0.002
+    k = 1.0
+    depth = 1.3
+    expected = 23.4
+    q = get_discharge(x, h, depth, n, s, k)
+    assert np.isclose(expected, q, rtol=0.1)
+
+
+@pytest.mark.xfail(reason="debug complex array input case")
+def test_get_discharge_compound_rectangular():
+    # adapted from https://www.youtube.com/watch?v=BJZ73WWEc3M
+
+    # manually sum over discharge for each rectangle
+    n = 0.016
+    s = 0.001
+    k = 1.0
+    expected = 3.3
+    q0 = get_discharge(3, 0.2, roughness=n, slope=s, conv=k)
+    q1 = get_discharge(1.5, 1.4, roughness=n, slope=s, conv=k)
+    q2 = get_discharge(3, 0.2, roughness=n, slope=s, conv=k)
+    sm = q0 + q1 + q2
+    assert np.isclose(expected, sm, rtol=1e-2)
+
+    # check single rectangular segment with array input for x and h
+    x = np.array([0, 0, 3, 3])
+    h = np.array([0, -0.2, -0.2, 0])
+    q = get_discharge(x, h, 0, n, s, k)
+    assert np.isclose(q, q0, rtol=1e-3)
+
+    # check multiple rectangles with array input
+    x = np.array([0, 0, 3, 3, 4.5, 4.5, 7.5, 7.5])
+    h = np.array([0, -0.2, -0.2, -1.4, -1.4, -0.2, -0.2, 0])
+    q = get_discharge(x, h, 0, n, s, k)
+    assert np.isclose(expected, q, rtol=1e-2)
+
+
+def test_get_depth():
+    # adapted from https://www.youtube.com/watch?v=t9ywTXEcScE
+
+    width = 5
+    x = np.array([0, 0, width, width])
+    h = np.array([2, 0, 0, 2])
+    q = 6.5
+    n = 0.02
+    s = 0.0005
+    k = 1.0
+    expected = 1.3
+    depth = get_depth(x, h, q, n, s, k)
+    assert np.isclose(expected, depth, rtol=1e-1)
+
+
+@pytest.mark.skip(reason="todo")
+def test_get_depths():
+    pass
diff --git a/flopy/utils/channel_utils.py b/flopy/utils/channel_utils.py
new file mode 100644
index 0000000000..9294f9bd01
--- /dev/null
+++ b/flopy/utils/channel_utils.py
@@ -0,0 +1,536 @@
+from typing import Optional, Union
+
+import numpy as np
+
+# power for Manning's hydraulic radius term
+_mpow = 2.0 / 3.0
+
+
+def get_segment_wetted_station(
+    x0,
+    x1,
+    h0,
+    h1,
+    depth,
+):
+    """
+    Computes the start and end stations of the wetted portion of a channel segment.
+    If the segment is fully submerged the start and end stations are the segment's
+    x-coordinates, unchanged. If the segment is partially submerged, the start/end
+    x-coordinates of the wetted subsegment are obtained using the segment's slope.
+    If the segment is wholly above the water surface, the end coordinate is set to
+    the value of the start coordinate (zero length in effect discards the segment).
+
+    Parameters
+    ----------
+    x0: x-coordinate of the segment's first point
+    x1: x-coordinate of the segment's second point
+    h0: elevation of the segment's first point
+    h1: elevation of the segment's second point
+    depth: elevation of the channel water surface
+
+    Returns
+    -------
+        The start and end stations (x-coordinates) of the wetted subsegment.
+    """
+
+    # calculate the minimum and maximum segment endpoint elevation
+    hmin = min(h0, h1)
+    hmax = max(h0, h1)
+
+    # if the water surface elevation is less than or equal to the
+    # minimum segment endpoint elevation, station length is zero.
+    if depth <= hmin:
+        x1 = x0
+
+    # if water surface is between hmin and hmax, station length is less than x1 - x0
+    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_segment_wetted_perimeter(
+    x0: float,
+    x1: float,
+    h0: float,
+    h1: float,
+    depth: float,
+):
+    """
+    Computes the wetted length of the channel cross-section segment given
+    the channel depth. See https://en.wikipedia.org/wiki/Wetted_perimeter.
+
+    Parameters
+    ----------
+    x0: x-coordinate of the segment's first point
+    x1: x-coordinate of the segment's second point
+    h0: elevation of the segment's first point
+    h1: elevation of the segment's second point
+    depth: elevation of the channel water surface
+
+    Returns
+    -------
+        The wetted length of the segment
+    """
+
+    # -- calculate the minimum and maximum elevation
+    hmin = min(h0, h1)
+    hmax = max(h0, h1)
+
+    if depth <= hmin:
+        return 0
+
+    # -- calculate the wetted perimeter for the segment
+    xlen = x1 - x0
+    if xlen < 0:
+        raise ValueError("x1 must be greater than x0")
+
+    if xlen == 0:
+        return 0
+
+    if xlen > 0.0:
+        if depth > hmax:
+            hlen = hmax - hmin
+        else:
+            hlen = depth - hmin
+    else:
+        if depth > hmin:
+            hlen = min(depth, hmax) - hmin
+        else:
+            hlen = 0.0
+
+    return np.sqrt(xlen**2.0 + hlen**2.0)
+
+
+def get_segment_wetted_area(
+    x0: float, x1: float, h0: float, h1: float, depth: float
+):
+    """
+    Computes the cross-sectional area above a segment of a channel with the
+    given geometry and filled to the given depth. Area above the segment is
+    defined as that of the quadrilateral formed by the channel segment, two
+    vertical sides, and the water surface.
+
+    Parameters
+    ----------
+    x0: x-coordinate of the segment's first point
+    x1: x-coordinate of the segment's second point
+    h0: elevation of the segment's first point
+    h1: elevation of the segment's second point
+    depth: elevation of the channel water surface
+
+    Returns
+    -------
+        The wetted area above the segment
+    """
+
+    # -- calculate the minimum and maximum elevation
+    hmin = min(h0, h1)
+    hmax = max(h0, h1)
+
+    if depth <= hmin:
+        return 0
+
+    # 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 hmax
+        area += (0.5 if hmax != hmin else 1) * (hmax - hmin)
+
+    return area
+
+
+def get_wetted_area(
+    x: np.ndarray,
+    h: np.ndarray,
+    depth: float,
+    verbose=False,
+):
+    """
+    Computes the cross-sectional wetted area of an open channel filled to a given depth.
+    The channel's geometry is provided as pairs of coordinates (x, h) defining boundary
+    segments. The segments are assumed to be connected in the order provided and arrays
+    for x and h must be the same length. See https://en.wikipedia.org/wiki/Wetted_area.
+
+    Parameters
+    ----------
+    x: x-coordinates of the channel's boundary segments
+    h: elevations of the channel's boundary segments
+    depth: elevation of the channel water surface
+    verbose: if True, print debugging information
+
+    Returns
+    -------
+        The wetted area of the channel
+    """
+
+    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_segment_wetted_station(x0, x1, h0, h1, depth)
+
+            # get wetted area
+            a = get_segment_wetted_area(x0, x1, h0, h1, depth)
+            area += a
+
+            if verbose:
+                print(
+                    f"{idx}->{idx + 1} ({x0},{x1}) - "
+                    f"perimeter={x1 - x0} - area={a}"
+                )
+
+    return area
+
+
+def get_wetted_perimeter(
+    x: np.ndarray,
+    h: np.ndarray,
+    depth: float,
+    verbose=False,
+):
+    """
+    Computes the cross-sectional wetted perimeter of an open channel filled to a given depth.
+    The channel's geometry is provided as pairs of coordinates (x, h) defining segments along
+    the channel's boundary. The segments are assumed connected in the order given, and arrays
+    x and h must share the same length. See https://en.wikipedia.org/wiki/Wetted_perimeter.
+
+    Parameters
+    ----------
+    x: x-coordinates of the channel's boundary segments
+    h: elevations of the channel's boundary segments
+    depth: elevation of the channel water surface
+    verbose: if True, print debugging information
+
+    Returns
+    -------
+        The wetted perimeter of the channel
+    """
+
+    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_segment_wetted_station(x0, x1, h0, h1, depth)
+
+            # get wetted perimeter
+            perimeter += get_segment_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 get_discharge_rect(
+    width: float,
+    depth: float,
+    roughness: float,
+    slope: float,
+    conv: float = 1.0,
+):
+    """
+    Calculates the discharge (volumetric flow rate) for a rectangular channel.
+
+    Parameters
+    ----------
+    width: the channel's width
+    depth: the channel's water depth
+    roughness: Manning's roughness coefficient
+    slope: the channel's slope
+    conv: the conversion factor
+
+    Returns
+    -------
+        The channel's discharge
+    """
+
+    if depth <= 0.0:
+        raise ValueError("Depth must be positive.")
+
+    if width <= 0.0:
+        raise ValueError("Width must be positive.")
+
+    if roughness <= 0:
+        raise ValueError(f"Roughness coefficient must be positive")
+
+    if slope < 0:
+        raise ValueError(f"Slope must be non-negative")
+
+    area = width * depth
+    perim = 2.0 * depth + width
+    radius = area / perim
+    # previously depth**_mpow was used below
+    return conv * area * radius**_mpow * slope**0.5 / roughness
+
+
+def get_discharge(
+    x: Union[int, float, np.ndarray],
+    h: Union[int, float, np.ndarray],
+    depth: Optional[float] = None,
+    roughness: float = 0.01,
+    slope: float = 0.001,
+    conv: float = 1.0,
+):
+    """
+    Calculates the discharge according to Manning's equation for a
+    channel with the provided cross-sectional geometry, slope, and
+    roughness. See https://en.wikipedia.org/wiki/Manning_formula.
+
+    Parameters
+    ----------
+    x: x-coordinates of the channel's boundary segments
+    h: elevations of the channel's boundary segments
+    depth: elevation of the channel water surface
+    roughness: Manning's roughness coefficient
+    slope: the channel's slope
+    conv: the conversion factor
+
+    Returns
+    -------
+        The channel's discharge
+    """
+
+    if roughness <= 0:
+        raise ValueError(f"Roughness coefficient must be positive")
+
+    if slope < 0:
+        raise ValueError(f"Slope must be non-negative")
+
+    if isinstance(x, int) or isinstance(x, float):
+        if not (isinstance(h, int) or isinstance(h, float)):
+            raise ValueError(
+                "Channel width and height must be integers or floats"
+            )
+
+        if x <= 0 or h <= 0:
+            raise ValueError(f"Channel width and height must be positive")
+
+        if depth is None:
+            depth = h
+
+        q = get_discharge_rect(x, depth, roughness, slope, conv)
+
+    elif isinstance(x, np.ndarray):
+        if not isinstance(h, np.ndarray):
+            raise ValueError("If x is an array, h must be an array")
+
+        if x.shape != h.shape:
+            raise ValueError(f"Arrays x and h must be the same shape")
+
+        if len(x) < 2:
+            raise ValueError(f"Arrays x and h must have at least 2 elements")
+
+        if isinstance(roughness, float):
+            roughness = (
+                np.ones(
+                    x.shape if isinstance(x, np.ndarray) else (1,), dtype=float
+                )
+                * roughness
+            )
+
+        q = 0.0
+        for i0 in range(x.shape[0] - 1):
+            in1 = i0 - 1
+            i1 = i0 + 1
+
+            perim = get_segment_wetted_perimeter(
+                x[i0], x[i1], h[i0], h[i1], depth
+            )
+            if perim == 0:
+                continue
+
+            # if left neighbor is vertical and descending, add its length
+            if i0 > 0:
+                left_len = x[in1] - x[i0]
+            else:
+                left_len = None
+            if left_len == 0 and h[in1] > h[i0]:
+                perim += min(depth, h[in1]) - h[i0]
+
+            # if right neighbor is vertical and ascending, add its length
+            if i0 < x.shape[0] - 2:
+                right_len = x[i1 + 1] - x[i1]
+                if right_len == 0 and h[i1 + 1] > h[i1]:
+                    perim += min(depth, h[i1 + 1]) - h[i1]
+
+            area = get_segment_wetted_area(x[i0], x[i1], h[i0], h[i1], depth)
+            radius = area / perim
+            q += conv * area * radius**_mpow * slope**0.5 / roughness[i0]
+
+    else:
+        raise ValueError(f"Invalid type for x: {type(x)}")
+
+    return q
+
+
+def get_depth(
+    x: np.ndarray,
+    h: np.ndarray,
+    q: float,
+    roughness: float = 0.01,
+    slope: float = 0.001,
+    conv: float = 1.0,
+    dd: float = 1e-4,
+    max_iter: int = 100,
+    verbose=False,
+):
+    """
+    Uses Manning's equation to approximate depth of flow for an open channel
+    given its cross-sectional geometry, a roughness coefficient, the rate of
+    flow, and the channel slope. An optional conversion factor may be given.
+    A single roughness coefficient or distinct coefficients for each segment
+    can also be provided. An iterative method is used to solve for the depth
+    and the number of iterations can be capped (by default, the max is 100).
+
+    Parameters
+    ----------
+    x: x-coordinates of the channel's boundary segments
+    h: elevations of the channel's boundary segments
+    q: the flow rate
+    roughness: Manning's roughness coefficient
+    slope: the channel's slope
+    conv: the conversion factor
+    dd: the increment used for approximating the derivative
+    verbose: if True, print debugging information
+
+    Returns
+    -------
+        The depth of flow
+    """
+
+    d0 = 0.0
+    q0 = get_discharge(
+        x,
+        h,
+        d0,
+        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 = get_discharge(
+            x,
+            h,
+            d0 + dd,
+            roughness=roughness,
+            slope=slope,
+            conv=conv,
+        )
+        dq = q1 - q0
+        if dq != 0.0:
+            derv = dd / (q1 - q0)
+        else:
+            derv = 0.0
+        d0 -= derv * r
+        q0 = get_discharge(
+            x,
+            h,
+            d0,
+            roughness=roughness,
+            slope=slope,
+            conv=conv,
+        )
+        r = q0 - q
+
+        iter += 1
+        if verbose:
+            print(f"iteration {iter:>2d} - residual={r}")
+        if iter > max_iter:
+            break
+    return d0
+
+
+def get_depths(
+    flows: Union[float, np.ndarray],
+    x: np.ndarray,
+    h: np.ndarray,
+    roughness: Union[float, np.ndarray] = 0.01,
+    slope: float = 0.001,
+    conv: float = 1.0,
+    dd: float = 1e-4,
+):
+    """
+    Uses Manning's equation to approximate depth of flow for an open channel
+    given its cross-sectional geometry, a roughness coefficient, one or more
+    flow rates, and the channel's slope. An optional conversion factor can be
+    provided, which is used to convert units of length from metric system to
+    another measurement system (e.g. for English units specify factor 1.49).
+    A single flow rate can be provided, or an array of flow rates. An array
+    is returned containing the depth of flow for each provided flow rate, in
+    the order given. A single roughness coefficient can be provided, or else
+    distinct values for each segment of the channel boundary can be provided.
+
+    Parameters
+    ----------
+    flows: the flow rate(s) to calculate depths for
+    x: x-coordinates of the channel's boundary segments
+    h: elevations of the channel's boundary segments
+    roughness: Manning's roughness coefficient
+    slope: the channel's slope
+    conv: the conversion factor for converting from metric to another unit
+    dd: the increment used for approximating the derivative
+
+    Returns
+    -------
+        An array of depths of flow (one element for each provided in flows)
+    """
+
+    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] = get_depth(
+            x,
+            h,
+            q,
+            roughness=roughness,
+            slope=slope,
+            conv=conv,
+            dd=dd,
+            verbose=False,
+        )
+
+    return depths