Hysteresis Methods

Individual hysteresis analysis modules implementing scientifically validated methods.

HARP Method

calculate_harp_metrics(df_obs, time_col='Etime', discharge_col='Q', concentration_col='C', intersection_method='auto')[source]

Calculate HARP metrics from discharge and concentration time series.

Parameters:

  • df_obs (pd.DataFrame) – Observed data with time, discharge, and concentration columns

  • time_col (str) – Column names for time, discharge, and concentration

  • discharge_col (str) – Column names for time, discharge, and concentration

  • concentration_col (str) – Column names for time, discharge, and concentration

  • intersection_method (str, default 'auto') – ‘auto’, ‘shapely’, ‘bentley-ottmann’, or ‘simple’

Returns:

  • metric_df (pd.DataFrame) – HARP metrics: - area: Total hysteresis area (magnitude) with sign of dominant portion

    Positive = diluting, Negative = enriching

    • residual: Change in scaled concentration (end - start)

      Positive = enriching, Negative = diluting

    • area_lower/area_upper: Areas of lower/upper portions (if intersection found)

    • peak_Q/peak_C: Scaled peak timings (0-1)

    • peaktime_Q/peaktime_C: Absolute peak timings (days)

    • radius_equiv: Equivalent circle radius

  • df_all (pd.DataFrame) – Processed data with scaled values and phase information

Calculate hysteresis metrics using the HARP (Hysteresis Analysis of Rising and falling Peaks) method developed by Roberts et al. (2023).

Method Overview:

HARP uses empirical classification based on:

  • Peak timing difference between Q and C

  • Loop area and shape

  • Residual (end-state deviation)

Returns:

Tuple of (metrics_dict, dataframe):

  • metrics_dict: Dictionary with area, residual, peak timings, classification

  • dataframe: Processed data with normalized values and classifications

Reference:

Roberts, M.E. et al. (2023). HARP: A new method for hysteresis analysis. Hydrological Processes.

harp_plot(df_processed, metrics)[source]

Visualize hysteresis loop df_processed :: df_all output from calculate_harp_metrics metrics :: metrics output from calculate_harp_metrics returns :: fig -> can be shown, saved, or further be processed

Create HARP visualization showing C-Q loop with peak timing markers.

Displays:

  • C-Q hysteresis loop

  • Peak Q and Peak C markers

  • Time progression coloring

  • Loop area shading

Zuecco Method

calculate_zuecco_metrics(df_obs, time_col='Etime', discharge_col='Q', concentration_col='C')[source]

Calculate Zuecco hysteresis index and classification.

Parameters:

  • df_obs (pd.DataFrame) – Observed data with time, discharge, and concentration columns

  • time_col (str) – Column names for time, discharge, and concentration

  • discharge_col (str) – Column names for time, discharge, and concentration

  • concentration_col (str) – Column names for time, discharge, and concentration

Returns:

  • metrics_df (pd.DataFrame) – Zuecco metrics: - h_index: Hysteresis index (sum of differential areas) - hyst_class: Hysteresis class (0-8) - min_diff_area: Minimum differential area - max_diff_area: Maximum differential area

  • df_all (pd.DataFrame) – Processed time series data with: - Q, C: Original discharge and concentration values - QS, CS: Normalized (0-1) discharge and concentration - Qphase, Cphase: Rising/falling phase indicators - switchpoints: Peak markers (gQ, gC, lQ, lC) Plus x_fixed interpolation stored in attrs: - attrs[‘x_fixed_interp’]: DataFrame with x_fixed, y_rise, y_fall, diff_area

Calculate hysteresis index using the Zuecco et al. (2016) integration method.

Method Overview:

Integration-based approach calculating differential areas between rising and falling limbs:

  • h-index: Sum of differential areas

  • 9-class classification system (0-8)

  • Quantifies magnitude and direction

Classification:

Class

Description

0

Linear / no hysteresis

1-4

Clockwise variants

5-8

Counter-clockwise variants

Returns:

Tuple of (metrics_dict, dataframe):

  • metrics_dict: h_index, hyst_class, min/max differential areas

  • dataframe: Processed data with cumulative areas and classifications

Reference:

Zuecco, G. et al. (2016). A versatile index to characterize hysteresis between hydrological variables at the runoff event timescale. Hydrological Processes, 30(9), 1449-1466.

zuecco_plot(df_all, metrics)[source]

Visualize Zuecco hysteresis analysis results in a two-panel plot: - Left: Hysteresis loop with rising/falling limbs - Right: Differential area plot

df_all :: Processed DataFrame from calculate_zuecco_metrics metrics :: Metrics DataFrame from calculate_zuecco_metrics returns :: fig -> can be shown, saved, or further be processed

Create Zuecco visualization showing C-Q loop with differential areas.

Displays:

  • C-Q hysteresis loop

  • Rising/falling limb differentiation

  • Differential area shading

  • h-index annotation

Lloyd/Lawler Methods

calculate_lawlerlloyd_metrics(df_obs, time_col='Etime', discharge_col='Q', concentration_col='C')[source]

Calculate Lawler-Lloyd hysteresis index.

Parameters:

  • df_obs (pd.DataFrame) – Observed data with time, discharge, and concentration columns

  • time_col (str) – Column names for time, discharge, and concentration

  • discharge_col (str) – Column names for time, discharge, and concentration

  • concentration_col (str) – Column names for time, discharge, and concentration

Returns:

  • metrics_df (pd.DataFrame) – Lawler-Lloyd metrics: - mean_HI: Mean hysteresis index across percentiles - median_HI: Median hysteresis index - HI_range: Range of HI values (max - min) - mean_abs_HI: Mean of absolute HI values

  • df_all (pd.DataFrame) – Processed time series data with: - Q, C: Original discharge and concentration values - QS, CS: Normalized (0-1) discharge and concentration - Qphase, Cphase: Rising/falling phase indicators - switchpoints: Peak markers (gQ, gC, lQ, lC) Plus HI percentile data stored in attrs: - attrs[‘HI_percentiles’]: DataFrame with Hi, RLi, FLi at each QS percentile

Calculate hysteresis indices using percentile-based methods from Lloyd et al. (2016) and Lawler et al. (2006).

Method Overview:

Two complementary approaches:

  1. HInew (Lloyd 2016) - RECOMMENDED

    Difference-based method:

    HI = (C_rise - C_fall) / C_mid

    • Symmetric range: [-1, 1]

    • Better for comparisons

    • More interpretable

  2. HIL (Lawler 2006) - ORIGINAL

    Ratio-based method:

    HI = (C_rise / C_fall) - 1  if C_rise > C_fall
HI = (-1 / (C_rise / C_fall)) + 1  otherwise

    • Asymmetric range

    • More sensitive at extremes

Both methods sample at 9 discharge percentiles (0.1, 0.2, …, 0.9).

Returns:

Tuple of (metrics_dict, dataframe):

  • metrics_dict: mean/median HInew, mean/median HIL, ranges

  • dataframe: Processed data with percentile samples and indices

References:

  • Lloyd, C.E.M. et al. (2016). Using hysteresis analysis of high-resolution water quality monitoring data. Hydrology and Earth System Sciences, 20, 2705-2719.

  • Lawler, D.M. et al. (2006). Turbidity dynamics during spring storm events. Science of the Total Environment, 360, 109-126.

lloyd_plot(df_all, metrics)[source]

Visualize Lawler-Lloyd hysteresis analysis in a two-panel plot: - Left: Hysteresis loop with rising/falling limbs - Right: HI index across discharge percentiles

df_all :: Processed DataFrame from calculate_lawlerlloyd_metrics metrics :: Metrics DataFrame from calculate_lawlerlloyd_metrics returns :: fig -> can be shown, saved, or further be processed

Create Lloyd/Lawler visualization showing percentile samples and HI values.

Displays:

  • C-Q hysteresis loop

  • Percentile sample points (rising/falling)

  • HInew and HIL indices

  • Hysteresis direction

Comparison and Selection

When to Use Each Method

HARP

  • Best for qualitative interpretation

  • Easy to understand (peak timing)

  • Good for identifying flushing vs. dilution

  • Use when process timing is important

Zuecco

  • Best for quantitative comparisons

  • Detects complex/mixed patterns

  • 9-class system captures subtle variations

  • Use for detailed pattern characterization

Lloyd/Lawler (HInew)

  • Best for standardized reporting

  • Symmetric range facilitates comparison

  • Standard in literature

  • Use for publication-quality quantification

Convergent Evidence

High Confidence - All three methods agree on direction and approximate magnitude

Moderate Confidence - Two methods agree, third differs slightly

Low Confidence - Methods disagree significantly

When methods disagree, this indicates:

  • Complex or mixed hysteresis pattern

  • Insufficient data points

  • Ambiguous event structure

  • Need for additional investigation

Usage Example

Compare All Three Methods:

import hygcs as gcs

# Calculate all methods
metrics = gcs.calculate_all_hysteresis_metrics(
    data,
    time_col='datetime',
    discharge_col='Q',
    concentration_col='C'
)

# Extract results
print("HARP:")
print(f"  Area: {metrics['harp_metrics']['area']:.4f}")
print(f"  Peak timing diff: {metrics['harp_metrics']['peaktime_C'] -
                               metrics['harp_metrics']['peaktime_Q']:.2f} days")

print("\nZuecco:")
print(f"  h-index: {metrics['zuecco_metrics']['h_index']:.4f}")
print(f"  Class: {metrics['zuecco_metrics']['hyst_class']}")

print("\nLloyd/Lawler:")
print(f"  HInew: {metrics['lloyd_metrics']['mean_HInew']:.4f}")
print(f"  HIL: {metrics['lloyd_metrics']['mean_HIL']:.4f}")
print(f"  Direction: {metrics['classifications']['lloyd_direction']}")

Visualize Each Method:

# HARP visualization
harp_m, harp_df = gcs.calculate_harp_metrics(data, ...)
fig_harp = gcs.harp_plot(harp_df, harp_m)
fig_harp.show()

# Zuecco visualization
zuecco_m, zuecco_df = gcs.calculate_zuecco_metrics(data, ...)
fig_zuecco = gcs.zuecco_plot(zuecco_df, zuecco_m)
fig_zuecco.show()

# Lloyd visualization
lloyd_m, lloyd_df = gcs.calculate_lawlerlloyd_metrics(data, ...)
fig_lloyd = gcs.lloyd_plot(lloyd_df, lloyd_m)
fig_lloyd.show()

Data Requirements

All three methods require:

  • Minimum 5 data points (10-15 recommended)

  • Complete Q and C time series

  • No missing values during event

  • Identifiable rising and falling limbs

For best results:

  • 20-30 points covering full event

  • Well-sampled peaks

  • Synchronized Q and C measurements

  • Quality-controlled data (outliers removed)

See Also