diff --git a/python/metrics/src/hydrotools/metrics/signature.py b/python/metrics/src/hydrotools/metrics/signature.py
index f0de6233..855e511c 100644
--- a/python/metrics/src/hydrotools/metrics/signature.py
+++ b/python/metrics/src/hydrotools/metrics/signature.py
@@ -3,6 +3,7 @@
 from dataclasses import dataclass
 import numpy as np
 import numpy.typing as npt
+from arch.bootstrap import IIDBootstrap
 
 @dataclass
 class FlowDurationCurve:
@@ -18,14 +19,26 @@ class FlowDurationCurve:
         Array of sorted values that correspond to probabilities. Typically
         the y-axis (vertical axis) in plots. Often depicted in log-scale for
         streamflow.
+    quantiles: numpy.ndarray, optional
+        Optional probability or sequence of probabilities that indicate what
+        the additional traces in `quantile_values` represent.
+    quantile_values: numpy.ndarray, optional
+        Optional array containing additional traces that represent alternative
+        values to plot on the y-axis. Typically used to represent the
+        boundaries of a confidence interval. Shape will be (len(quantiles), 
+        len(values)).
     """
     probabilities: np.ndarray
     values: np.ndarray
+    quantiles: np.ndarray | None = None
+    quantile_values: np.ndarray | None = None
 
 def compute_fdc(
         series: npt.ArrayLike,
         alpha: float = 0.0,
-        beta: float = 0.0
+        beta: float = 0.0,
+        quantiles: npt.ArrayLike | None = None,
+        boostrap_samples: int = 1000
         ) -> FlowDurationCurve:
     """
     Computes an empirical flow duration curve derived from series. The flow
@@ -44,6 +57,25 @@ def compute_fdc(
     beta: float, optional, default 0.0
         Plotting position parameter used to adjust plotting positions. Valid
         range is [0,1].
+    quantiles: array-like, optional
+        Probability or sequence of probabilities of quantiles to compute.
+        Values must be between 0 and 1 inclusive. Uses bootstrap resampling to
+        estimate additional quantiles of the ordered series. Passed directly
+        to numpy.quantile. Often used to construct confidence intervals. For
+        example, setting `quantiles=[0.025, 0.975]` will result in two
+        additional traces of values that represent the 95% confidence
+        interval.
+    bootstrap_samples: int, optional, default 1000
+        Used with `quantiles`. Indicates the number of bootstrap samples
+        to generate for the posterior distribution of values.
+    
+    Returns
+    -------
+    FlowDurationCurve dataclass. Probabilities are accessible through
+    FlowDurationCurve.probabilities. Sorted values are available through
+    FlowDurationCurve.values. If the `quantiles` option is used, additional
+    traces will be available through FlowDurationCurve.quantile_values, which
+    directly correspond to the values in FlowDurationCurve.quantiles.
     """
     # Validate
     length = len(series)
@@ -53,7 +85,16 @@ def compute_fdc(
 
     # Build FDC
     p = (np.arange(1, length+1) - alpha) / (length+1-alpha-beta)
-    return FlowDurationCurve(
+    fdc = FlowDurationCurve(
         p,
         np.sort(series)[::-1]
     )
+    if quantiles is None:
+        return fdc
+    
+    # Generate additional quantiles
+    bs = IIDBootstrap(series, seed=2025)
+    posterior = bs.apply(lambda a: np.sort(a)[::-1], boostrap_samples)
+    fdc.quantiles = quantiles
+    fdc.quantile_values = np.quantile(posterior, quantiles, axis=0)
+    return fdc
diff --git a/python/metrics/tests/test_signature.py b/python/metrics/tests/test_signature.py
index 2fa2865d..387a860b 100644
--- a/python/metrics/tests/test_signature.py
+++ b/python/metrics/tests/test_signature.py
@@ -8,6 +8,7 @@ def test_compute_fdc():
     fdc = signature.compute_fdc(y)
     assert np.all(fdc.probabilities <= 1.0)
     assert np.all(fdc.values == 1)
+    assert len(fdc.values) == 100
 
     y = np.linspace(-1.0, 100.0, 102)
     fdc = signature.compute_fdc(y)
@@ -15,3 +16,26 @@ def test_compute_fdc():
     assert fdc.values[-1] < fdc.values[0]
     assert np.min(fdc.values) == fdc.values[-1]
     assert np.max(fdc.values) == fdc.values[0]
+    assert fdc.quantiles == None
+    assert fdc.quantile_values == None
+
+    fdc = signature.compute_fdc(y, alpha=0.5, beta=0.5)
+    assert np.all(fdc.probabilities <= 1.0)
+    assert fdc.values[-1] < fdc.values[0]
+    assert np.min(fdc.values) == fdc.values[-1]
+    assert np.max(fdc.values) == fdc.values[0]
+
+    fdc = signature.compute_fdc(y, quantiles=[0.025, 0.975])
+    assert fdc.quantiles[0] == 0.025
+    assert fdc.quantiles[1] == 0.975
+    assert fdc.quantile_values.shape[0] == 2
+    assert fdc.quantile_values.shape[1] == 102
+    assert fdc.values[-1] < fdc.values[0]
+    assert np.min(fdc.values) == fdc.values[-1]
+    assert np.max(fdc.values) == fdc.values[0]
+    assert fdc.quantile_values[0][-1] < fdc.quantile_values[0][0]
+    assert np.min(fdc.quantile_values[0]) == fdc.quantile_values[0][-1]
+    assert np.max(fdc.quantile_values[0]) == fdc.quantile_values[0][0]
+    assert fdc.quantile_values[1][-1] < fdc.quantile_values[1][0]
+    assert np.min(fdc.quantile_values[1]) == fdc.quantile_values[1][-1]
+    assert np.max(fdc.quantile_values[1]) == fdc.quantile_values[1][0]