Methods

The page describes the methods used in pyBenD to analyze river meanders. It includes details on how pyBenD create Centerline objects, identifies meander bends and their characteristic points, track centerlines and bends through time, and calculates morphometric and kinematic parameters.

Centerline Creation

pyBenD provides methods to create Centerline objects from various input data formats (see Inputs section for more details). The Centerline class stores the channel centerline as a series of connected channel points (ClPoint). For subsequent analyzes, the centerline is optionally resampled and smoothed to reduce digitizing noise and to provide a more uniform point spacing prior to computing geometric attributes.

The main steps are:

Resampling and path smoothing

  1. Curvilinear abscissa computation

    The centerline cumulative distance (curvilinear abscissa) is computed from the input XY coordinates. This provides an estimate of total centerline length and is used to determine how many points are required after resampling.

  2. Resampling with a parametric spline

    If spacing > 0, the centerline is resampled with a parametric spline so that the output contains either:

    • a target spacing between consecutive points (default behavior), or

    • a fixed number of points when use_fix_nb_points=True. In this case, the number of points can be set with the spacing parameter (integer value).

    If spacing == 0, the centerline is not resampled and the original point positions are used.

  3. Smoothing of the resampled XY path

    When resampling is enabled (spacing > 0), the resampled X and Y coordinates are smoothed using a Savitzky-Golay filter (polynomial order 3, mode='nearest'). The smoothing window is derived from the requested smoothing distance:

    • window size (in points) is int(smooth_distance / spacing),

    • the window is forced to be odd and at least 5 points.

  4. Recompute derived geometry on the new points

    After resampling/smoothing, pyBenD recomputes (i) the curvilinear abscissa of the updated point set and (ii) the local normal vector components (NORMAL_X, NORMAL_Y) from neighboring points.

Property interpolation

If interpol_props=True, existing properties stored in the input DataFrame are transferred to the resampled points.

For each resampled point, pyBenD:

  • finds the two closest points on the original centerline (in the XY plane),

  • computes their distances to the resampled point,

  • interpolates each property using a distance-weighted average.

Some fields are intentionally excluded from interpolation because they are recomputed for the resampled geometry: curvilinear abscissa and XY coordinates.

Curvature calculation and curvature smoothing

If compute_curvature=True, curvature is computed along the (resampled and smoothed) centerline using a local 3-point method: the curvature at point i depends on the coordinates of points (i-1, i, i+1).

(See Curvature in the Glossary.)

To reduce local variability, a second curvature series is then computed by filtering the raw curvature. Two filters are available:

  • a moving uniform filter (window size in points), or

  • a Savitzky-Golay filter (polynomial order 2).

For the Savitzky-Golay curvature filter, the window size is forced to be odd and at least 5 points.

Meander Bend Identification

In pyBenD, meander bends are identified directly from a Centerline object once a filtered curvature series has been computed on the centerline points (Curvature_filtered property).

The workflow implemented in find_bends() is:

  1. Detect inflection points from curvature sign changes

    Inflection points are extracted from the filtered curvature array using find_inflection_points(). A point at index i is considered an inflection point when the curvature changes sign in its neighborhood (using a 3-point condition based on sums of consecutive curvatures):

    \[\mathrm{sign}(C_{i-1}+C_i) \neq \mathrm{sign}(C_i+C_{i+1})\]

    Very close inflection points are then filtered so that two consecutive inflection indices are separated by at least a small lag (in points). In the current implementation, lag=2 is used.

    The first and last centerline points are always added as bend delimiters if they were not detected as inflection points.

  2. Create bends as segments between consecutive inflection points

    Each bend is defined as the set of centerline points between two consecutive inflection points (upstream and downstream). Internally, this creates one Bend object per interval.

  3. Compute bend properties (side, validity, apex, and middle point)

    For each bend, pyBenD computes:

    • Bend side (BendSide): based on the sign of the sum of filtered curvature along the bend. Positive sum gives BendSide.UP, negative gives BendSide.DOWN.

    • Bend validity: based on a sinuosity criterion computed between the two inflection points.

      (See Sinuosity in the Glossary.)

      • bend arc-length L is measured as the curvilinear abscissa difference,

      • bend chord length D is the Euclidean distance between inflection-point

        coordinates.

      The bend sinuosity is then:

      \[S = \frac{L}{D}\]

      A bend is marked valid if S >= sinuo_thres and S < 10 (the upper bound is currently hard-coded).

    • Bend middle point (Bend.pt_center): the midpoint of the inflection-point chord (the bend chord) (XY midpoint), with Z set to the average Z of the two inflection points.

    • Bend apex index (Bend.index_apex): by default, the bend apex is determined from the filtered curvature distribution (see below).

When multiple CPU cores are available, bend creation and per-bend property computation can be parallelized using multiprocessing.

Bend apex identification

pyBenD provides multiple methods to locate the apex of each bend. All methods return an index along the full centerline (not a local bend index).

Curvature-distribution method (default)

The default apex definition used in find_bend_apex() is based on the median of the cumulative filtered curvature distribution within the bend.

Let \(\kappa_i\) be the filtered curvature magnitude within the bend. The cumulative weighting function is built from \(\kappa_i^n\) (n is an exponent parameter):

\[F(j) = \frac{\sum_{i \le j} |\kappa_i|^n}{\sum_i |\kappa_i|^n}\]

The apex is the first index where \(F(j) > 0.5\). Larger values of n put more emphasis on the highest-curvature part of the bend.

Probability-from-weights method

Alternatively, an apex probability can be computed with set_bend_apex_probability_user_weights() and used to detect the apex with find_bend_apex_from_weights().

The apex probability is a weighted combination of three normalized components:

  • filtered curvature magnitude,

  • bend amplitude (distance to the bend chord),

  • a length term that favors the center of the bend over its ends.

User weights (curvature_weight, amplitude_weight, length_weight) are renormalized so that their sum equals 1. Inflection points are forced to have zero probability.

To turn probability into an apex index, local maxima of the probability curve are detected (scipy.signal.find_peaks with a peak height threshold in [0.6, 1.0]). If multiple peaks are found, the apex is taken as the mean index of the peaks.

Azimuth change method

The method implemented in find_bend_apex_limaye() follows Limaye (2025) and defines the apex from half of the total channel path angle variation between bend endpoints if total change if less than 90°. If the total change is greater than 180°, the apex is located at 90° of change from the upstream inflection point.

Practically, the algorithm uses the local normal vectors along the bend (Normal_x and Normal_y properties) and finds the index where the angular deviation from the first normal reaches half of the overall deviation.

Morphometric Parameter Calculation

Morphometric parameters are computed once bends and their characteristic points (inflection points, apex, and chord midpoint) have been identified on a Centerline. In practice, the computations are implemented in Morphometry and exposed through compute_bends_morphometry().

Unless stated otherwise, all bend-scale metrics are computed on the subset of valid bends (see bend validity in the previous section). Internally, each bend stores indices to the upstream and downstream inflection points (Bend.index_inflex_up, Bend.index_inflex_down), the apex index (Bend.index_apex), and the chord midpoint (Bend.pt_center).

Morphometric parameters (bend scale)

Let o be the two inflection points delimiting a bend and x the bend apex (inflection points and bend apex). Let \(s\) denote the curvilinear abscissa along the centerline.

  • Bend arc length (MorphometricNames ARC_LENGTH):

    The bend arc length \(L_b\) is the curvilinear distance between the two inflection points:

    \[L_b = s_{\mathrm{down}} - s_{\mathrm{up}}\]

    This is computed from the Curvilinear_abscissa property stored on the centerline points within the bend.

  • Bend wavelength (MorphometricNames WAVELENGTH):

    The bend wavelength \(\lambda_b\) is the Euclidean distance between the upstream and downstream inflection points:

    \[\lambda_b = \|\mathbf{p}_{\mathrm{down}} - \mathbf{p}_{\mathrm{up}}\|\]

  • Bend sinuosity (MorphometricNames SINUOSITY):

    Computed as sinuosity \(S = L_b / \lambda_b\) when \(\lambda_b > 0\).

  • Bend amplitude (MorphometricNames AMPLITUDE):

    The bend amplitude \(A_b\) is the orthogonal (perpendicular) distance from the apex point to the bend chord (the line through the two inflection points):

    \[A_b = d_\perp\left(\mathbf{p}_{\mathrm{apex}},\ \mathbf{p}_{\mathrm{up}},\ \mathbf{p}_{\mathrm{down}}\right)\]

    In code this is computed with an orthogonal-distance routine applied to the apex and the two inflection-point coordinates.

  • Bend extension (MorphometricNames EXTENSION):

    The bend extension \(E_b\) is the distance from the apex to the midpoint of the chord (bend center):

    \[E_b = \|\mathbf{p}_{\mathrm{apex}} - \mathbf{p}_{\mathrm{center}}\|\]

    where \(\mathbf{p}_{\mathrm{center}}\) is the chord midpoint (XY midpoint; Z set to the mean Z of the two inflection points).

  • Radius of curvature (MorphometricNames RADIUS_CURVATURE):

    The bend radius of curvature is computed from the filtered curvature at the apex (see curvature filtering above):

    \[R = \frac{1}{|\kappa_{\mathrm{apex}}|}\]

    If \(|\kappa_{\mathrm{apex}}| = 0\), the radius is reported as NaN.

  • Asymmetry coefficient (MorphometricNames ASYMMETRY):

    Following Howard and Hemberger (1991), the asymmetry coefficient compares the curvilinear distances from the apex to the upstream and downstream inflection points:

    \[AC = \frac{L_{\mathrm{up}} - L_{\mathrm{down}}}{L_b}\]

    where \(L_{\mathrm{up}}\) is the arc length from the upstream inflection point to the apex, and \(L_{\mathrm{down}}\) is the arc length from the apex to the downstream inflection point.

  • Roundness (MorphometricNames ROUNDNESS):

    The roundness coefficient (Schwenk et al., 2015) is computed from the distribution of filtered curvature magnitude \(|\kappa|\) within the bend:

    \[\mathrm{Roundness} = \frac{\max(|\kappa|)}{\mathrm{mean}(|\kappa|)}\]

  • Skewness (MorphometricNames SKEWNESS):

    The skewness corresponds to Pearson’s skewness coefficient of the curvature spatial distribution. The curvilinear abscissa is first normalized to \(s' \in [0, 1]\) within the bend, and the skewness is evaluated as:

    \[\mathrm{Skewness} = \frac{3\,(\mu - s'_{\mathrm{apex}})}{\sigma}\]

    where \(\mu\) and \(\sigma\) are the mean and standard deviation of the curvature spatial distribution computed from \(|\kappa|\) and \(s'\), and \(s'_{\mathrm{apex}}\) is the normalized curvilinear abscissa at the apex index.

Morphometric parameters (meander scale, Leopold & Wolman)

For interior bends (i.e., bends that have both a preceding and a following bend), pyBenD also provides two meander-scale metrics consistent with Leopold and Wolman (1960). These use apex points of three consecutive bends:

  • Meander wavelength (WAVELENGTH_LEOPOLD): Euclidean distance between the apex of the previous bend and the apex of the next bend.

  • Meander amplitude (AMPLITUDE_LEOPOLD): orthogonal distance from the current bend apex to the straight line joining the previous and next bend apexes.

Centerline and Bend tracking

In pyBenD, kinematic metrics require linking geometric entities through time. This is done in two stages implemented in CenterlineCollection (a centerline collection):

  1. Centerline matching: associate channel points between consecutive

    centerlines (point-to-point tracking).

  2. Bend connection: associate bends between consecutive centerlines and

    group them into BendEvolution objects (bend-to-bend tracking; see bend evolution).

Centerline matching (channel-point tracking)

Channel-point tracking is computed with match_centerlines(). For each pair of consecutive centerlines (age prev_key and key), pyBenD builds a pairwise cost matrix and solves a Dynamic Time Warping (DTW) problem to obtain an alignment between point indices. Centerlines are tracked sequentially from the oldest to the youngest.

The methodology for channel point tracking is inspired from zsylvester/meanderpy.

Cost matrix components

For points \(i\) on the current centerline and \(j\) on the previous one, three non-negative distances can be combined:

  • Planform distance: Euclidean distance between \((x_i, y_i)\) and \((x_j, y_j)\). If the distance exceeds dmax, it is set to a very large penalty so that DTW avoids matching distant points. dmax is usually set to a meander wavelength.

  • Curvature distance: difference in filtered curvature magnitude \(\left|\,|\kappa_i| - |\kappa_j|\,\right|\).

  • Velocity-perturbation distance (optional): absolute difference of the VELOCITY_PERTURBATION property after Savitzky-Golay filtering (window_length=window, polyorder=2), when the property exists.

Each component is normalized by its maximum (when non-zero) so that the weights can be interpreted consistently. The final DTW cost is:

\[C = w_v\,C_v + w_\kappa\,C_\kappa + w_d\,C_d\]

with user weights vel_perturb_weight, curvature_weight and distance_weight.

DTW alignment and stored links

The DTW implementation is provided by the dtw-python package. DTW is applied to the cost matrix (dtw.dtw) with a configurable step pattern (pattern, default "asymmetric"). The resulting warping path is converted into a list of indices in the previous centerline (dtw.warp).

pyBenD then stores the matching in both directions:

  • On the current centerline: index_cl_pts_prev_centerline stores, for each current point index, the matched index in the previous centerline (or -1 when no acceptable match exists).

  • On the previous centerline: index_cl_pts_next_centerline stores, for each previous point index, the list of current indices that map to it (one-to-many links may occur due to DTW).

In addition, the per-point link lists on each channel point ClPoint are updated: cl_pt_index_prev and cl_pt_index_next.

Centerline tracking is generally time-consuming. When multiple CPU cores are available, matching of multiple consecutive pairs can be parallelized.

Bend connection (bend-to-bend tracking)

Once bends have been identified on each centerline (see previous sections), pyBenD can connect bends through time and build BendEvolution objects using connect_bends(). Three strategies are implemented (method argument):

1) Apex-distance connection

With method=BendConnectionMethod.APEX (BendConnectionMethod), bends are connected by minimizing the Euclidean distance between successive apex points (Bend.index_apex). The search is performed backward through time. For each valid bend at the current age:

  • candidate trajectories are restricted to those whose last bend belongs to the immediately previous age and has the same bend side (Bend.side),

  • the closest apex is selected if its distance is below dmax; otherwise a new trajectory is started.

2) Centroid-distance connection

With method=BendConnectionMethod.CENTROID (BendConnectionMethod), the same logic is applied but the distance is computed between bend centroids (Bend.pt_centroid) rather than apex points.

3) Connection from point matching

With method=BendConnectionMethod.MATCHING (BendConnectionMethod), bend association is derived from the channel-point tracking result and therefore requires that match_centerlines() has been run first.

For each bend of age key and its successive centerline next_key, pyBenD counts how many channel points inside the bend (inflection points excluded) connect to each bend in next_key via the DTW links (cl_pt_index_next). Connected bends are retained when the fraction of connected points exceeds a threshold (currently 0.2).

Reciprocal connections are stored on bends as lists of unique identifiers (bend_uid_next and bend_uid_prev). Finally, connected components of this directed/undirected connection graph are collected recursively to build the BendEvolution list.

Validity of BendEvolution objects

After bend connection, each BendEvolution contains a mapping from ages to bend indices (bend_indexes). A BendEvolution can be flagged as valid if it contains enough time steps; this is controlled by bend_evol_validity.

Section across channel belt

pyBenD can build 2D stratigraphic sections across a channel belt by intersecting a set of user-defined (or automatically generated) section lines with all centerlines stored in a CenterlineCollection.

The workflow is implemented in find_points_on_sections() and produces a list of Section objects, each storing the channel position through time in a section-aligned coordinate system.

Defining section lines

Sections are defined in plan view as straight line segments (shapely.LineString) stored in CenterlineCollection.section_lines. Two approaches are available:

  • Manual definition with set_section_lines(), by providing lists of start and end coordinates (pts_start and pts_end). This is the approach used in notebooks/bend_kinematics_analysis.ipynb where section endpoints are chosen so that each section intersects the centerlines of a given channel belt.

  • Automatic definition with create_section_lines(). The method uses the youngest (last) centerline and creates one section per valid interior bend (first and last bends are skipped). The section line depends on the selected CreateSectionMethod:

    • MIDDLE: line from bend apex to chord midpoint (Bend.pt_center).

    • CENTROID: line from bend apex to bend centroid (Bend.pt_centroid).

    • APEX: line from the bend apex to the midpoint between neighboring bend

      apexes (fallbacks use bend midpoints when apexes are undefined).

Finding channel location on a section

Given a section line, pyBenD searches for its intersections with every centerline polyline. For each age:

  1. A local search window is built around the section line to avoid testing all

    segments.

  2. Each candidate centerline segment (two consecutive

    ClPoint) is intersected with the section line.

  3. When an intersection exists, channel-point properties are interpolated at the

    intersection location using the relative distance along the segment.

Each intersection produces an ChannelCrossSection isoline (type CHANNEL) whose reference point is the interpolated ClPoint. An idealized cross-section geometry is then generated with complete_channel_shape() (parabolic shape whose width and depth come from the reference point properties).

Finally, if a section line intersects enough centerlines (controlled by the thres argument), a Section object is created and stored in CenterlineCollection.sections. During this process, the bend intersected by each channel position is also identified from the intersected point index and the bend is tagged with the section index (used later for kinematic summaries).

Section coordinate system

Within a Section, the channel position at each age is stored relative to the first intersected channel position according to the coordinates (d, dz), where:

  • \(d\) is the signed lateral position along the section, measured relative to the first channel occurrence on the section.

  • \(dz\) is the vertical offset relative to that first occurrence.

The sign of \(d\) is determined from the user-provided flow direction (flow_dir): the migration vector between the reference channel position and the current one is projected onto the direction perpendicular to flow_dir to assign the left/right side of migration. These \((d, dz)\) coordinates form Section.isolines_origin and are used to compute stacking-pattern classes (get_stacking_pattern_type()), apparent channel displacements (channel_apparent_displacements()), and section-averaged mobility metrics (section_averaged_channel_displacements()).

Stacking pattern classification

When a channel belt is observed on a 2D stratigraphic section, the successive intersections of the channel with the section can be interpreted as a trajectory of the channel (here, represented by the channel reference point of each intersected ChannelCrossSection).

pyBenD classifies this trajectory into stacking-pattern types using get_stacking_pattern_type(). The implemented classification follows the workflow used in:

Lemay, M., Grimaud, J. L., Cojan, I., Rivoirard, J., and Ors, F. (2024). Submarine channel stacking patterns controlled by the 3D kinematics of meander bends. Geological Society, London, Special Publications, SP540-2022-143. doi.org/10.1144/SP540-2022-143.

Step classification

The classification operates on the ordered list Section.isolines_origin of channel points for each time step containing pairs (d, dz). For each consecutive pair of time steps, pyBenD computes a signed lateral migration increment:

\[\Delta d_k = d_{k} - d_{k-1}\]

and discretizes it into a ternary sequence:

  • 0 if \(|\Delta d_k| < \mathtt{mig\_threshold}\) (treated as aggradation-only),

  • +1 if \(\Delta d_k \ge \mathtt{mig\_threshold}\),

  • -1 if \(\Delta d_k \le -\mathtt{mig\_threshold}\).

The parameter mig_threshold therefore controls when lateral motion is considered significant. In notebooks/bend_kinematics_analysis.ipynb it is set proportionally to channel width (mig_threshold = 0.0125 * width).

From steps to stacking-pattern types

Let \(f_0\), \(f_{+}\), and \(f_{-}\) be the fractions of steps classified as 0, +1 and -1 respectively. A first-order decision is made from these fractions using frac_threshold (default 0.95; also used in the notebook):

  • Aggradation Only (StackingPatternType.AGGRADATION) if \(f_0\) exceeds

    frac_threshold.

  • One Way (StackingPatternType.ONE_WAY) if either \(f_{+}\) or \(f_{-}\)

    exceeds frac_threshold.

When the trajectory contains both aggradation steps and a dominant migration direction such that \((f_0 + f_{+})\) or \((f_0 + f_{-})\) exceeds frac_threshold, pyBenD distinguishes between:

  • Aggradation then One Way (StackingPatternType.AGGRAD_ONE_WAY): a

    sustained early aggradation phase followed by a sustained migration phase.

  • One Way (StackingPatternType.ONE_WAY): migration dominates from early

    time steps.

In practice, the implementation groups consecutive steps of the dominant sign and of 0 into runs, and compares the run lengths to a small fraction of the record length (begin_threshold, default 0.1). This is the criterion used to decide whether the initial part of the trajectory is long enough to be treated as a distinct aggradation phase.

When neither direction dominates and aggradation is not overwhelming, the classification is based on the number of direction reversals after removing aggradation steps (0):

  • Two Way (StackingPatternType.TWO_WAYS) when the non-zero migration

    steps can be reduced to two main runs of opposite sign.

  • Multiple Ways (StackingPatternType.MULTIPLE_WAYS) when more than two

    alternating runs are required.

The resulting class is stored in Section.stacking_pattern_type and can be used downstream to compute “bend-scale” vs “whole-trajectory” section-averaged mobility metrics (see section_averaged_channel_displacements()).

Kinematic Parameter Calculation

In pyBenD, kinematic parameters quantify how the channel (or parts of it) move through time, once geometric entities have been linked between successive centerlines.

Two complementary families of kinematic measurements are commonly used:

  • Real (planform) kinematics: measured from successive mapped/simulated centerlines, after point-to-point tracking (DTW matching).

  • Apparent (section-based) kinematics: measured from channel positions recorded on stratigraphic sections, in a section-aligned coordinate system.

The notebooks notebooks/seine_river_migration_analysis.ipynb and notebooks/bend_kinematics_analysis.ipynb illustrate these computations.

Prerequisites

Real kinematics require:

  1. A CenterlineCollection built

    from multiple centerlines (ages/iterations).

  2. Channel-point tracking computed with

    match_centerlines().

Depending on the metric, you may also need bends to be detected on each centerline (find_bends=True) and/or cross-sections to be defined.

Channel-point migration (Seine river example)

The simplest kinematic metric is the migration distance of each channel point between two consecutive centerlines.

After DTW matching, each centerline of age \(t\) stores an index mapping index_cl_pts_prev_centerline such that for a point \(i\) on the current centerline, the matched point index on the previous centerline is \(j\).

Let \(\mathbf{p}_t(i)\) be the XY position of point \(i\) at age \(t\) and \(\mathbf{p}_{t-1}(j)\) the matched XY position at age \(t-1\). The planform migration vector and distance are:

\[\Delta \mathbf{p}(i) = \mathbf{p}_t(i) - \mathbf{p}_{t-1}(j)\]
\[d_{\mathrm{mig}}(i) = \|\Delta \mathbf{p}(i)\|\]

In notebooks/seine_river_migration_analysis.ipynb, this is implemented by looping over the matched indices and computing cpf.distance(pt1, pt2) for each valid match.

This point-wise migration can then be summarized (mean, median, histogram) or mapped along the centerline.

From trajectories to kinematic components

Many bend-scale metrics are computed by first building a trajectory: a time ordered list of positions of a tracked entity (apex, bend-average position, any point index, etc.).

Given a list of points l_pt (2D or 3D), compute_point_displacements() computes:

  • local_disp: incremental displacements between successive time steps, with columns (dX, dY, dZ, dMig).

  • whole_disp: displacement between first and last trajectory points, again as (DX, DY, DZ, DMig).

Internally, the XY displacement is projected into a user-defined local basis (a change of coordinates). This is what allows the notebooks to express displacements across or along a chosen direction.

Bend-average migration (bend_kinematics_analysis example)

The bend-averaged kinematics shown in notebooks/bend_kinematics_analysis.ipynb follow this logic:

  1. Select a bend on a reference centerline (often the youngest one).

  2. Build the list of point indices that belong to the bend

    (from upstream to downstream inflection points).

  3. Step through ages and use index_cl_pts_prev_centerline to follow each

    point to its match on the next centerline.

  4. For each consecutive pair of ages, compute the mean migration vector across

    all bend points (excluding unmatched points):

    \[\Delta \bar{\mathbf{p}} = \frac{1}{N}\sum_{k=1}^N \left(\mathbf{p}_t(k) - \mathbf{p}_{t-1}(k')\right)\]

    This produces one mean vector per time step; integrating those vectors over time yields a bend-average trajectory.

  5. Optionally smooth/resample that trajectory (the notebook uses a

    resampling step followed by a Savitzky-Golay filter).

  6. Compute incremental and total displacements with

    compute_point_displacements().

The notebook stores bend-scale metrics from the resulting whole_disp:

  • DxBend, DyBend: absolute values of the displacement components in the chosen local basis.

  • DmigBend: total planform migration magnitude (DMig).

Apex migration (bend_kinematics_analysis example)

The apex trajectory is computed similarly, but instead of averaging over all points of the bend, the notebook tracks only the apex index (Bend.index_apex) through time using the DTW index mapping. The resulting trajectory is then processed with the same displacement routine, yielding DxApex, DyApex, and DmigApex.

Section-based (apparent) kinematics and mobility numbers

When stratigraphic sections are available, pyBenD can compute apparent displacements along a section using channel_apparent_displacements().

Given successive channel positions on the section (stored as Section.isolines_origin), the method computes incremental lateral and vertical displacements:

\[\Delta x = x_{t} - x_{t-1},\qquad \Delta z = z_{t} - z_{t-1}\]

and the local (dimensionless) stratigraphic mobility number:

\[M_{\mathrm{sb}} = \frac{\Delta x}{\Delta z}\,\frac{H}{W}\]

where \(W\) is channel width and \(H\) is channel depth (see section_averaged_channel_displacements(), after Jobe et al., 2016).

In notebooks/bend_kinematics_analysis.ipynb, analogous mobility numbers are also reported for real planform kinematics by combining a planform migration magnitude (DmigApex or DmigBend) with a section-derived thickness estimate (DzApex / DzBend):

\[M = \frac{D_{\mathrm{mig}}}{D_z}\,\frac{H}{W}\]

This allows comparing planform-derived kinematics (“real”) to section-derived kinematics (“recorded”) on the same dimensionless scale.