Note

This notebook is located in the ./examples directory of the gwtransport repository.

Residence Time Distribution Analysis#

Learning Objectives#

Calculate residence time distributions from aquifer pore volume distributions
Understand the difference between infiltration to extraction and extraction to infiltration residence time perspectives
Analyze temporal variations in groundwater residence times
Apply residence time concepts to contaminant transport and water quality assessment

Overview#

This notebook demonstrates how to calculate residence time distributions from aquifer pore volume distributions and flow rates. Residence time quantifies how long water spends in the aquifer, which is crucial for understanding contaminant transport and water quality.

Two Perspectives#

Infiltration to Extraction: How long until infiltrating water is extracted?
Extraction to Infiltration: How long ago was extracted water infiltrated?

Key Assumptions#

Stationary pore volume distribution (steady streamlines)
Known aquifer heterogeneity (from Example 1)

Background Reading#

Residence Time — Time in aquifer between infiltration and extraction
Pore Volume Distribution — Aquifer heterogeneity modeling
Retardation Factor — Slower movement due to sorption

Theoretical Background#

Residence Time#

Residence time represents the travel time of water through the aquifer:

\[t_{residence} = \frac{V_{pore}}{Q} \cdot R_f\]

where \(V_{pore}\) is pore volume [m³], \(Q\) is flow rate [m³/day], and \(R_f\) is the retardation factor [-].

Residence times vary over time due to seasonal flow rate changes and extraction rate variations. See residence time and retardation factor for background.

[1]:

import matplotlib.pyplot as plt
import numpy as np

import gwtransport.residence_time as rt
from gwtransport import gamma as gamma_utils
from gwtransport.examples import generate_temperature_example_data
from gwtransport.utils import step_plot_coords

# Set random seed for reproducibility
np.random.seed(42)
plt.style.use("seaborn-v0_8-whitegrid")

print("Libraries imported successfully")

Libraries imported successfully

1. Data Setup and Aquifer Parameters#

We use aquifer parameters typically obtained from Example 1 (temperature-based characterization) and generate flow data for residence time analysis.

[2]:

df, tedges = generate_temperature_example_data(
    date_start="2020-01-01",
    date_end="2025-05-31",
    flow_mean=120.0,  # Base flow rate [m³/day]
    flow_amplitude=40.0,  # Seasonal flow variation [m³/day]
    flow_noise=5.0,  # Random daily fluctuations [m³/day]
    cin_method="soil_temperature",  # Use real soil temperature data
    aquifer_pore_volume_gamma_mean=8000.0,  # True mean pore volume [m³]
    aquifer_pore_volume_gamma_std=400.0,  # True standard deviation [m³]
)

print("Dataset Summary:")
print(f"Period: {df.index[0].date()} to {df.index[-1].date()}")
print(f"Mean flow: {df['flow'].mean():.1f} m³/day")
print(f"Flow range: {df['flow'].min():.1f} - {df['flow'].max():.1f} m³/day")
print(
    f"Aquifer parameters: {df.attrs['aquifer_pore_volume_gamma_mean']:.0f} ± {df.attrs['aquifer_pore_volume_gamma_std']:.0f} m³"
)

print("Dataset Summary:")
print(f"Period: {df.index[0].date()} to {df.index[-1].date()}")
print(f"Mean flow: {df['flow'].mean():.1f} m³/day")
print(f"Mean infiltration temperature: {df['cin'].mean():.1f} °C")
print(f"Mean extraction temperature: {df['cout'].mean():.1f} °C")
print(f"True mean pore volume: {df.attrs['aquifer_pore_volume_gamma_mean']:.1f} m³")
print(f"True std deviation: {df.attrs['aquifer_pore_volume_gamma_std']:.1f} m³")

Dataset Summary:
Period: 2020-01-01 to 2025-05-31
Mean flow: 111.8 m³/day
Flow range: 17.7 - 172.4 m³/day
Aquifer parameters: 8000 ± 400 m³
Dataset Summary:
Period: 2020-01-01 to 2025-05-31
Mean flow: 111.8 m³/day
Mean infiltration temperature: 11.8 °C
Mean extraction temperature: 11.8 °C
True mean pore volume: 8000.0 m³
True std deviation: 400.0 m³

[3]:

# Gamma pore-volume distribution -- consumed directly by the closed-form residence-time
# functions (no discretization). The residence-time percentile band evaluates `full` at the
# gamma quantile pore volumes, since residence time is monotone in pore volume.
pv_mean = df.attrs["aquifer_pore_volume_gamma_mean"]
pv_std = df.attrs["aquifer_pore_volume_gamma_std"]
quantiles = [1, 10, 90, 99]  # percentiles for the residence-time uncertainty band
pv_alpha, pv_beta = gamma_utils.mean_std_loc_to_alpha_beta(mean=pv_mean, std=pv_std)
pv_quantiles = gamma_utils.gamma_dist.ppf(np.array(quantiles) / 100, pv_alpha, scale=pv_beta)

print(f"Gamma pore volume: {pv_mean:.0f} +/- {pv_std:.0f} m3  (alpha={pv_alpha:.0f}, beta={pv_beta:.1f})")
print(f"Pore-volume {quantiles} percentiles: {np.round(pv_quantiles).astype(int)} m3")
print(f"Mean residence time (no retardation): {pv_mean / df['flow'].mean():.1f} days")
print(
    f"Mean residence time (with retardation): {pv_mean / df['flow'].mean() * df.attrs['retardation_factor']:.1f} days"
)

Gamma pore volume: 8000 +/- 400 m3  (alpha=400, beta=20.0)
Pore-volume [1, 10, 90, 99] percentiles: [7099 7492 8517 8960] m3
Mean residence time (no retardation): 71.5 days
Mean residence time (with retardation): 143.1 days

2. Infiltration to Extraction Residence Time Analysis#

Calculate how long infiltrating water takes to be extracted. We compute residence times for both water flow (conservative tracer) and thermal transport (with retardation).

[4]:

print("Computing infiltration to extraction residence times (closed-form gamma APVD)...")

# Mean residence time over the gamma pore-volume distribution (closed form; default warm-start).
rt_ite_water_mean = rt.gamma(
    flow=df.flow,
    tedges=tedges,
    cout_tedges=tedges,
    mean=pv_mean,
    std=pv_std,
    retardation_factor=1.0,  # Conservative tracer (water flow)
    direction="infiltration_to_extraction",
)
rt_ite_thermal_mean = rt.gamma(
    flow=df.flow,
    tedges=tedges,
    cout_tedges=tedges,
    mean=pv_mean,
    std=pv_std,
    retardation_factor=df.attrs["retardation_factor"],  # Heat transport (slower)
    direction="infiltration_to_extraction",
)

# Residence-time percentile band: tau is monotone in pore volume, so the q-th residence-time
# percentile is `full` at the q-th pore-volume quantile (water flow). Rows follow `quantiles`.
rt_ite_water_bands = rt.full(
    flow=df.flow,
    tedges=tedges,
    cout_tedges=tedges,
    aquifer_pore_volumes=pv_quantiles,
    retardation_factor=1.0,
    direction="infiltration_to_extraction",
)
print(
    f"Residence-time series: {rt_ite_water_mean.shape[0]} output bins; band = {rt_ite_water_bands.shape[0]} quantiles"
)

Computing infiltration to extraction residence times (closed-form gamma APVD)...
Residence-time series: 1978 output bins; band = 4 quantiles

[5]:

# Assemble the mean residence time and the percentile band into the dataframe
quantile_headers = [f"rt_infiltration_to_extraction_water_{q}%" for q in quantiles]

df["rt_infiltration_to_extraction_water_mean"] = rt_ite_water_mean
df["rt_infiltration_to_extraction_thermal_mean"] = rt_ite_thermal_mean
df[quantile_headers] = rt_ite_water_bands.T  # rows of `full` are the quantile residence times

print("Infiltration to extraction residence time statistics calculated")
print(f"Mean water residence time: {df['rt_infiltration_to_extraction_water_mean'].mean():.1f} days")
print(f"Mean thermal residence time: {df['rt_infiltration_to_extraction_thermal_mean'].mean():.1f} days")
print(
    f"Retardation effect: {df['rt_infiltration_to_extraction_thermal_mean'].mean() / df['rt_infiltration_to_extraction_water_mean'].mean():.1f}x"
)

Infiltration to extraction residence time statistics calculated
Mean water residence time: 75.4 days
Mean thermal residence time: 145.4 days
Retardation effect: 1.9x

Fraction Explained: Assessing Spin-Up#

During the initial period, some flow paths do not yet have sufficient flow history to compute a valid residence time. The :func:~gwtransport.residence_time.fraction_explained_gamma function quantifies, per output bin, the advective fraction of pore volumes that are out of spin-up (informed by the flow record). A value of 1.0 means all flow paths are fully informed. (This is an advective diagnostic; molecular diffusion and dispersion are not modelled here.)

[6]:

# Advective fraction of each output bin explained by the flow record (closed-form gamma APVD)
frac_explained_forward = rt.fraction_explained_gamma(
    flow=df.flow,
    tedges=tedges,
    cout_tedges=tedges,
    mean=pv_mean,
    std=pv_std,
    retardation_factor=1.0,
    direction="infiltration_to_extraction",
)

# First date the gamma distribution is (effectively) fully explained by the flow record
fully_explained_mask = frac_explained_forward >= 1.0
if fully_explained_mask.any():
    first_fully_explained = df.index[fully_explained_mask][0]
    print(f"Spin-up period ends: {first_fully_explained.date()}")
else:
    print("Warning: the gamma distribution is not fully explained within the data period")

print(f"Fraction explained at start: {frac_explained_forward[0]:.2f}")
print(f"Fraction explained at end: {frac_explained_forward[-1]:.2f}")

Spin-up period ends: 2020-01-01
Fraction explained at start: 1.00
Fraction explained at end: 0.00

3. Infiltration to Extraction Residence Time Visualization#

Visualize how residence times vary over time and compare water vs. thermal transport.

[7]:

fig, ax = plt.subplots(2, 1, figsize=(14, 10), sharex=True)

# Flow rate subplot - convert to step format
xstep_flow, ystep_flow = step_plot_coords(tedges, df.flow)
ax[0].plot(xstep_flow, ystep_flow, label="Discharge rate", color="steelblue", linewidth=1.2)
ax[0].set_ylabel("Flow Rate [m³/day]")
ax[0].set_title("Infiltration to Extraction Residence Time Analysis", fontsize=14, fontweight="bold")
ax[0].legend(loc="upper right")
ax[0].grid(True, alpha=0.3)

# Residence time subplot - convert all series to step format
xstep_rt, ystep_rt_water = step_plot_coords(tedges, df.rt_infiltration_to_extraction_water_mean)
_, ystep_rt_thermal = step_plot_coords(tedges, df.rt_infiltration_to_extraction_thermal_mean)

ax[1].plot(
    xstep_rt,
    ystep_rt_water,
    label="Mean (water, R=1.0)",
    color="blue",
    linewidth=1.5,
)
ax[1].plot(
    xstep_rt,
    ystep_rt_thermal,
    label=f"Mean (thermal, R={df.attrs['retardation_factor']:.1f})",
    color="red",
    linewidth=1.5,
)

# Add uncertainty bounds
for q in quantiles:
    alpha_val = 0.4 if q in {10, 90} else 0.3
    linestyle = "--" if q in {1, 99} else "-."
    _, ystep_quantile = step_plot_coords(tedges, df[f"rt_infiltration_to_extraction_water_{q}%"])
    ax[1].plot(
        xstep_rt,
        ystep_quantile,
        label=f"{q}% quantile (water)",
        color="blue",
        alpha=alpha_val,
        linewidth=0.8,
        linestyle=linestyle,
    )

ax[1].set_ylabel("Residence Time [days]")
ax[1].set_xlabel("Date")
ax[1].legend(loc="upper right", ncol=2, fontsize=9)
ax[1].grid(True, alpha=0.3)

# Add fraction explained on twin axis
ax1b = ax[1].twinx()
ax1b.plot(
    *step_plot_coords(tedges, frac_explained_forward * 100),
    color="gray",
    linewidth=1.5,
    linestyle=":",
    alpha=0.8,
    label="Fraction explained",
)
ax1b.set_ylabel("Fraction Explained [%]")
ax1b.set_ylim(0, 105)
ax1b.legend(loc="lower right", fontsize=9)

# Add explanatory text
ax[1].text(
    0.02,
    0.02,
    "Infiltration to extraction: How many days until infiltrating water is extracted?",
    ha="left",
    va="bottom",
    transform=ax[1].transAxes,
    fontsize=11,
    bbox={"boxstyle": "round,pad=0.3", "facecolor": "wheat", "alpha": 0.8},
)

plt.tight_layout()

plt.show()

../_images/examples_02_Residence_Time_Analysis_13_0.png

4. Extraction to Infiltration Residence Time Analysis#

Calculate when currently extracted water was originally infiltrated. This perspective is crucial for source identification and age dating.

[8]:

print("Computing extraction to infiltration residence times (closed-form gamma APVD)...")

rt_eti_water_mean = rt.gamma(
    flow=df.flow,
    tedges=tedges,
    cout_tedges=tedges,
    mean=pv_mean,
    std=pv_std,
    retardation_factor=1.0,
    direction="extraction_to_infiltration",
)
rt_eti_thermal_mean = rt.gamma(
    flow=df.flow,
    tedges=tedges,
    cout_tedges=tedges,
    mean=pv_mean,
    std=pv_std,
    retardation_factor=df.attrs["retardation_factor"],
    direction="extraction_to_infiltration",
)
rt_eti_water_bands = rt.full(
    flow=df.flow,
    tedges=tedges,
    cout_tedges=tedges,
    aquifer_pore_volumes=pv_quantiles,
    retardation_factor=1.0,
    direction="extraction_to_infiltration",
)
print(f"Residence-time series: {rt_eti_water_mean.shape[0]} output bins")

Computing extraction to infiltration residence times (closed-form gamma APVD)...
Residence-time series: 1978 output bins

[9]:

quantile_headers_extraction_to_infiltration = [f"rt_extraction_to_infiltration_water_{q}%" for q in quantiles]

df["rt_extraction_to_infiltration_water_mean"] = rt_eti_water_mean
df["rt_extraction_to_infiltration_thermal_mean"] = rt_eti_thermal_mean
df[quantile_headers_extraction_to_infiltration] = rt_eti_water_bands.T

# Advective fraction explained for the extraction to infiltration direction (closed-form gamma APVD)
frac_explained_backward = rt.fraction_explained_gamma(
    flow=df.flow,
    tedges=tedges,
    cout_tedges=tedges,
    mean=pv_mean,
    std=pv_std,
    retardation_factor=1.0,
    direction="extraction_to_infiltration",
)

print("Extraction to infiltration residence time statistics calculated")
print(
    f"Mean extraction to infiltration water residence time: {df['rt_extraction_to_infiltration_water_mean'].mean():.1f} days"
)
print(
    f"Mean extraction to infiltration thermal residence time: {df['rt_extraction_to_infiltration_thermal_mean'].mean():.1f} days"
)

fully_explained_backward_mask = frac_explained_backward >= 1.0
if fully_explained_backward_mask.any():
    print(f"Spin-up period ends: {df.index[fully_explained_backward_mask][0].date()}")

Extraction to infiltration residence time statistics calculated
Mean extraction to infiltration water residence time: 74.7 days
Mean extraction to infiltration thermal residence time: 143.2 days
Spin-up period ends: 2020-05-23

5. Extraction to Infiltration Residence Time Visualization#

Show how the age of extracted water varies over time, revealing the temporal dynamics of groundwater flow.

[10]:

fig, ax = plt.subplots(2, 1, figsize=(14, 10), sharex=True)

# Flow rate subplot - convert to step format
xstep_flow, ystep_flow = step_plot_coords(tedges, df.flow)
ax[0].plot(xstep_flow, ystep_flow, label="Discharge rate", color="steelblue", linewidth=1.2)
ax[0].set_ylabel("Flow Rate [m³/day]")
ax[0].set_title("Extraction to Infiltration Residence Time Analysis", fontsize=14, fontweight="bold")
ax[0].legend(loc="upper right")
ax[0].grid(True, alpha=0.3)

# Residence time subplot - convert all series to step format
xstep_rt, ystep_rt_water = step_plot_coords(tedges, df.rt_extraction_to_infiltration_water_mean)
_, ystep_rt_thermal = step_plot_coords(tedges, df.rt_extraction_to_infiltration_thermal_mean)

ax[1].plot(
    xstep_rt,
    ystep_rt_water,
    label="Mean (water, R=1.0)",
    color="green",
    linewidth=1.5,
)
ax[1].plot(
    xstep_rt,
    ystep_rt_thermal,
    label=f"Mean (thermal, R={df.attrs['retardation_factor']:.1f})",
    color="orange",
    linewidth=1.5,
)

# Add uncertainty bounds for water
for q in quantiles:
    alpha_val = 0.4 if q in {10, 90} else 0.3
    linestyle = "--" if q in {1, 99} else "-."
    _, ystep_quantile = step_plot_coords(tedges, df[f"rt_extraction_to_infiltration_water_{q}%"])
    ax[1].plot(
        xstep_rt,
        ystep_quantile,
        label=f"{q}% quantile (water)",
        color="green",
        alpha=alpha_val,
        linewidth=0.8,
        linestyle=linestyle,
    )

ax[1].set_ylabel("Residence Time [days]")
ax[1].set_xlabel("Date")
ax[1].legend(loc="upper right", ncol=2, fontsize=9)
ax[1].grid(True, alpha=0.3)

# Add fraction explained on twin axis
ax1b = ax[1].twinx()
ax1b.plot(
    *step_plot_coords(tedges, frac_explained_backward * 100),
    color="gray",
    linewidth=1.5,
    linestyle=":",
    alpha=0.8,
    label="Fraction explained",
)
ax1b.set_ylabel("Fraction Explained [%]")
ax1b.set_ylim(0, 105)
ax1b.legend(loc="lower right", fontsize=9)

# Add explanatory text
ax[1].text(
    0.02,
    0.02,
    "Extraction to infiltration: How many days ago was the extracted water infiltrated?",
    ha="left",
    va="bottom",
    transform=ax[1].transAxes,
    fontsize=11,
    bbox={"boxstyle": "round,pad=0.3", "facecolor": "lightblue", "alpha": 0.8},
)

plt.tight_layout()

plt.show()

../_images/examples_02_Residence_Time_Analysis_18_0.png

Results & Discussion#

Flow Rate Dependencies#

Residence times show strong inverse correlation with flow rates — higher flows yield shorter residence times. Both infiltration-to-extraction and extraction-to-infiltration perspectives yield comparable averages but differ in temporal patterns.

Retardation Effects#

Thermal retardation factor of 2.0 effectively doubles residence times, which must be considered in heat transport applications.

Practical Applications#

Infiltration to extraction: Use for contaminant arrival time predictions and spill response planning
Extraction to infiltration: Use for source identification and age dating
Seasonal variations: Consider flow variations when designing monitoring and management strategies

Key Takeaways#

Two Perspectives: Infiltration-to-extraction and extraction-to-infiltration residence times answer different but complementary questions
Flow Rate Control: Higher flows lead to shorter residence times and faster transport
Retardation Effects: Must account for tracer-specific retardation in transport analysis
Spin-up Assessment: Use :func:~gwtransport.residence_time.fraction_explained_gamma to identify when model output is fully informed by input signal

Residence Time Distribution Analysis

Contents

Residence Time Distribution Analysis#

Learning Objectives#

Overview#

Two Perspectives#

Key Assumptions#

Background Reading#

Theoretical Background#

Residence Time#

1. Data Setup and Aquifer Parameters#

2. Infiltration to Extraction Residence Time Analysis#

Fraction Explained: Assessing Spin-Up#

3. Infiltration to Extraction Residence Time Visualization#

4. Extraction to Infiltration Residence Time Analysis#

5. Extraction to Infiltration Residence Time Visualization#

Results & Discussion#

Flow Rate Dependencies#

Retardation Effects#

Practical Applications#

Key Takeaways#