Numerical methods

Cross-section construction

Conveyance-storage separation

In 1D hydrodynamic models, flow through a cross-section is resolved assuming a cross-sectionally average velocity. This assumed that the entire cross-section is available to for conveyance. However in reality some parts of the cross-section do not contribute to flow. For example, sections of a river behind a levee where water is stagnant contribute to storage (volume), but not flow.

SOBEK enables distinction between 'flow area' and 'storage area'. fm2prof implements methods to resolve from 2D model output which cells add to the 'flow volume' within a control volume and which to the storage volume.

fm2prof implements two methods. The configuration parameter ConveyanceDetectionMethod is used to determine which method is used.

max_method A cell is considered flowing if the velocity magnitude is more than the average of the three higher flow velocities per outputmap multiplied by the relative velocity threshold OR if the flow velocity meets the absolute threshold absolute velocity threshold

mean_method Not recommended. Legacy method.

Parameters:

Name	Type	Description	Default
waterdepth	DataFrame	dataframe of a control volume with waterdepths per cel per output map	required
velocity	DataFrame	dataframe of a control volume with velocity magnitude per cel per output map	required

Returns:

Name	Type	Description
flow_mask	DataFrame	dataframe of a control volume with the flow condition per cel per output map. True means flowing, False storage.

Source code in fm2prof\CrossSection.py

def _distinguish_conveyance_from_storage(self, waterdepth: pd.DataFrame, velocity: pd.DataFrame) -> pd.DataFrame:
    """
    In 1D hydrodynamic models, flow through a cross-section is resolved assuming a 
    cross-sectionally average velocity. This assumed that the entire cross-section
    is available to for conveyance. However in reality some parts of the cross-section
    do not contribute to flow. For example, sections of a river behind a levee where
    water is stagnant contribute to storage (volume), but not flow. 

    SOBEK enables distinction between 'flow area' and 'storage area'. `fm2prof` implements
    methods to resolve from 2D model output which cells add to the 'flow volume' within a
    [control volume](glossary.md#control-volume) and which to the storage volume. 

    `fm2prof` implements two methods. The configuration parameter [`ConveyanceDetectionMethod`](configuration.md#exec-1--conveyancedetectionmethod) is used
    to determine which method is used.

    **`max_method`**
    A cell is considered flowing if the velocity magnitude is more than the average
    of the three higher flow velocities per outputmap multiplied by the 
    [`relative velocity threshold`](configuration.md#exec-1--relativevelocitythreshold) OR
    if the flow velocity meets the absolute threshold [`absolute velocity threshold`](configuration.md#exec-1--absolutevelocitythreshold)

    **`mean_method`**
    Not recommended. Legacy method.

    Parameters:
        waterdepth: dataframe of a control volume with waterdepths per cel per output map
        velocity:  dataframe of a control volume with velocity magnitude per cel per output map

    Returns:
        flow_mask: dataframe of a control volume with the flow condition per cel per output map. `True` means flowing, `False` storage. 
    """
    @staticmethod
    def max_velocity_method(waterdepth: pd.DataFrame, velocity: pd.DataFrame) -> pd.DataFrame:
        """
        This method was added in version 2.3 because the mean_velocity_method
        led to unreasonably high conveyance if the river was connected to 
        an inland harbour. 
        """
        # This condition may be redundant
        waterdepth_condition = waterdepth > 0

        # Determine maximum as the average of the top 3 flow velocities
        maxv = velocity.max()
        for i in velocity:
            maxv[i] = velocity[i].sort_values().iloc[-3:].mean()

        # Relative to max condition
        relative_velocity_condition = velocity > maxv*self.get_parameter(self.__cs_parameter_relative_threshold)

        # Absolute flow condition
        absolute_velocity_condition =  velocity > self.get_parameter(self.__cs_parameter_velocity_threshold)

        # Flow mask determines which cells are conveyance (TRUE)
        flow_mask = waterdepth_condition & (relative_velocity_condition | absolute_velocity_condition)

        return flow_mask

    @staticmethod
    def mean_velocity_method(waterdepth, velocity):
        """
        This was the default method < 2.3. This method leads to unreasonably 
        high conveyance if the river was connected to an inland harbour. 
        """
        # apply rolling average over the velocities
        # to smooth out extreme values
        velocity = velocity.rolling(
            window=10, min_periods=1, center=True).mean()

        flow_mask = (
            (waterdepth > 0)
            & (velocity > self.get_parameter(self.__cs_parameter_velocity_threshold))
            & (
                velocity
                > self.get_parameter(self.__cs_parameter_relative_threshold)
                * np.mean(velocity)
            )
        )

        return flow_mask

    match self.get_inifile().get_parameter(self.__cs_parameter_conveyance_detection_method):
        case 0:
            return mean_velocity_method(waterdepth, velocity)    
        case 1:
            return max_velocity_method(waterdepth, velocity)
        case _:
            self.set_logger_message('Invalid conveyance method. Defaulting to [1]', 'warning')
            return max_velocity_method(waterdepth, velocity)

Simplification

The cross-section geometry generated by fm2prof contains one point per output timestep in the 2D map file. This resolution is often too high given the complexity of the cross-sections, and results in very large input files for the 1D model. Therefore fm2prof includes a simplification algorithm that reduces the number of points while preservering the shape of the geometry. This algorithm reduces as many points until the number specified in MaximumPointsInProfile is reached.

We use the Visvalingam-Whyatt method of poly-line vertex reduction¹. The total width is leading for the simplification of the geometry meaning that the choice for which points to remove to simplify the geometry is based on the total width. Subsequently, the corresponding point are removed from the flow width.

---

1. Visvalingam, M and Whyatt J D (1993) "Line Generalisation by Repeated Elimination of Points", Cartographic J., 30 (1), 46 - 51 URL: http://web.archive.org/web/20100428020453/http://www2.dcs.hull.ac.uk/CISRG/publications/DPs/DP10/DP10.html Implemented vertex reduction methods: ↩

Parameters:

Name	Type	Description	Default
count_after	int	number of points in cross-section after application of this function	20

Source code in fm2prof\CrossSection.py

def reduce_points(self, count_after: int = 20) -> None:
    """
    The cross-section geometry generated by `fm2prof` contains one point per output
    timestep in the 2D map file. This resolution is often too high given the
    complexity of the cross-sections, and results in very large input files for the
    1D model. Therefore `fm2prof` includes a simplification algorithm that reduces
    the number of points while preservering the shape of the geometry. This algorithm
    reduces as many points until the number specified in
    `MaximumPointsInProfile` is reached.

    We use the Visvalingam-Whyatt method of poly-line vertex reduction[^1].
    The [total width](glossary.md#total-width) is leading for the simplification of the geometry meaning
    that the choice for which points to remove to simplify the geometry is based on
    the total width. Subsequently, the corresponding point are removed from the [flow width](glossary.md#flow-width).

    [^1]:
        Visvalingam, M and Whyatt J D (1993) "Line Generalisation by Repeated Elimination of Points", Cartographic J., 30 (1), 46 - 51 URL: http://web.archive.org/web/20100428020453/http://www2.dcs.hull.ac.uk/CISRG/publications/DPs/DP10/DP10.html
            Implemented vertex reduction methods:


    Parameters:
        count_after: number of points in cross-section after application of this function

    """

    n_before_reduction = self.get_number_of_vertices()

    points = np.array(
        [
            [self._css_z[i], self._css_total_width[i]]
            for i in range(n_before_reduction)
        ]
    )

    # The number of points is equal to n, it cannot be further reduced
    reduced_index = np.array([True] * n_before_reduction)

    if n_before_reduction > count_after:
        try:
            simplifier = PS.VWSimplifier(points)
            reduced_index = simplifier.from_number_index(count_after)
        except Exception as e:
            self.set_logger_message(
                "Exception thrown while using polysimplify: "
                + "{}".format(str(e)),
                "error",
            )

    # Write to attributes
    self.z = self._css_z[reduced_index]
    self.total_width = self._css_total_width[reduced_index]
    self.flow_width = self._css_flow_width[reduced_index]

    self.set_logger_message(
        "Cross-section reduced "
        + "from {} ".format(n_before_reduction)
        + "to {} points".format(len(self.total_width))
    )

    self._css_is_reduced = True

Lake identification

This algorithms determines whether a 2D cell should be marked as Lake.

Cells are marked as lake if the following conditions are both met: - the waterdepth on timestep LakeTimeSteps is positive - the waterdepth on timestep LakeTimeSteps is at least 1 cm higher than the waterlevel on timestep 0.

Next, the following steps are taken

• It is determined at what timestep the waterlevel in the lake starts rising. From that point on the lake counts as regular geometry and counts toward the total volume. A cell is considered active if its waterlevel has risen by 1 mm.
• A correction matrix is built that contains the 'lake water level' for each lake cell. This matrix is subtracted from the waterdepth to compute volumes.

Parameters:

Name	Type	Description	Default
waterdepth	DataFrame	a DataFrame containing all waterdepth output in the control volume	required

Returns:

Name	Type	Description
lake_mask	ndarray	mask of all cells that are a 'lake'
wet_not_lake_mask	ndarray	mask of all cells that are wet, but not a lake
lake_depth_correction	ndarray	the depth of a lake at the start of the 2D computation

Source code in fm2prof\CrossSection.py

def _identify_lakes(self, waterdepth:pd.DataFrame) -> np.ndarray:
    """
    This algorithms determines whether a 2D cell should
    be marked as [Lake](glossary.md#Lakes).

    Cells are marked as lake if the following conditions are both met:
    - the waterdepth on timestep [LakeTimeSteps](configuration.md#exec-1--laketimesteps) is positive
    - the waterdepth on timestep [LakeTimeSteps](configuration.md#exec-1--laketimesteps) is at least 1 cm higher than the waterlevel on timestep 0.

    Next, the following steps are taken

    - It is determined at what timestep the waterlevel in the lake starts rising. From that point on the lake counts as regular geometry and counts toward the total volume. A cell is considered active if its waterlevel has risen by 1 mm.
    - A correction matrix is built that contains the 'lake water level' for each lake cell. This matrix is subtracted from the waterdepth to compute volumes.


    Parameters:
        waterdepth: a DataFrame containing all waterdepth output in the [control volume](glossary.md#control-volume)

    Returns:
        lake_mask: mask of all cells that are a 'lake'
        wet_not_lake_mask: mask of all cells that are wet, but not a lake
        lake_depth_correction: the depth of a lake at the start of the 2D computation

    """
    # preallocate arrays
    plassen_depth_correction = np.zeros(waterdepth.shape, dtype=float)

    # check for non-rising waterlevels
    waterdepth_diff = np.diff(waterdepth, n=1, axis=-1)

    # find all wet cells
    wet_mask = waterdepth > 0

    # find all lakes
    lake_mask = (
        waterdepth.T.iloc[self.get_parameter(self.__cs_parameter_plassen_timesteps)]
        > 0
    ) & (
        np.abs(
            waterdepth.T.iloc[
                self.get_parameter(self.__cs_parameter_plassen_timesteps)
            ]
            - waterdepth.T.iloc[0]
        )
        <= 0.01
    )

    self.plassen_mask = lake_mask

    # Plassen_mask_time is to determine at whata timestep the lake starts rising again.
    plassen_mask_time = np.zeros((len(waterdepth.T), len(lake_mask)), dtype=bool)

    # At t=0, all lakes are inactive
    plassen_mask_time[0, :] = lake_mask

    # walk through dataframe in time, for each timestep check
    # when to unmask a plassen cell
    for i, diff in enumerate(waterdepth_diff.T):
        final_mask = reduce(
            np.logical_and, [(diff <= 0.001), (plassen_mask_time[i] == True)]
        )
        plassen_mask_time[i + 1, :] = final_mask

    plassen_mask_time = pd.DataFrame(plassen_mask_time).T

    # The depth of a lake is the waterdepth at t=0
    for i, depths in enumerate(waterdepth):
        plassen_depth_correction[lake_mask, i] = -waterdepth.T.iloc[0][lake_mask]

    # correct wet cells for plassen
    wet_not_plas_mask = reduce(
        np.logical_and, [(wet_mask == True), np.asarray(plassen_mask_time == False)]
    )

    return lake_mask, wet_not_plas_mask, plassen_depth_correction