From 4da61b738138890e023abdef4230faeb030905ae Mon Sep 17 00:00:00 2001
From: "Troy Frever, NOAA Affiliate" <troy.frever@noaa.gov>
Date: Mon, 27 Jan 2025 13:57:49 -0800
Subject: [PATCH] scale down the flow speed at blind channels and impoundments
 (#17)

---
 src/hydro.cpp | 35 +++++++++++++++++++++++++++++++----
 src/hydro.h   |  2 ++
 src/map.cpp   |  9 +++++++++
 src/map.h     |  4 ++++
 4 files changed, 46 insertions(+), 4 deletions(-)

diff --git a/src/hydro.cpp b/src/hydro.cpp
index 3a3471e..2600a4a 100644
--- a/src/hydro.cpp
+++ b/src/hydro.cpp
@@ -104,7 +104,7 @@ float HydroModel::getFlowSpeedAlong(Edge &edge) {
     // Calculate how colinear the flow vector and the given edge are
     float d = sqrt(dx * dx + dy * dy);
     float scalar_proj = u*(dx/d) + v*(dy/d);
-    return scalar_proj;
+    return scaledFlowSpeed(scalar_proj, *(edge.source));
 }
 
 // Get the current horizontal (E/W) flow velocity in m/s at the given node
@@ -118,13 +118,40 @@ float HydroModel::getCurrentV(DistribHydroNode &hydroNode) {
 }
 
 // Get the total flow velocity in m/s at the given node
+float HydroModel::getFlowSpeedAtHydroNode(DistribHydroNode &hydroNode) {
+    float currU = this->getCurrentU(hydroNode);
+    float currV = this->getCurrentV(hydroNode);
+    return sqrt(currU*currU + currV*currV);
+}
+
+float HydroModel::scaledFlowSpeed(const float velocity, const MapNode &node) {
+    if (!isBlindChannel(node.type) && !isImpoundment(node.type)) {
+        return velocity;
+    }
+    constexpr double M_TO_CM_CONV = 100.0;
+    const double hydroVelocity = this->getFlowSpeedAtHydroNode(this->hydroNodes[node.nearestHydroNodeID]);
+    const double hydroWidth = 0.00332119 * pow(hydroVelocity * M_TO_CM_CONV, 125.0 / 48.0);
+    const double blindChannelWidth = sqrt(node.area);
+    double scalar = blindChannelWidth / hydroWidth;
+
+    if (scalar > 1.0) {
+        scalar = 1.0;
+    }
+    if (isImpoundment(node.type)) {
+        constexpr double IMPOUNDMENT_MIN_FLOW_ADDL_SCALAR = 0.1;
+        scalar = IMPOUNDMENT_MIN_FLOW_ADDL_SCALAR * scalar;
+    }
+
+    const double scaledFlowSpeed = scalar * velocity;
+    return static_cast<float>(scaledFlowSpeed);
+}
+
 float HydroModel::getFlowSpeedAt(MapNode &node) {
     if (this->useSimData) {
         return isDistributary(node.type) ? this->simDistFlow / (this->getDepth(node) * sqrt(node.area)) : 0.0f;
     }
-    float currU = this->getCurrentU(this->hydroNodes[node.nearestHydroNodeID]);
-    float currV = this->getCurrentV(this->hydroNodes[node.nearestHydroNodeID]);
-    return sqrt(currU*currU + currV*currV);
+    const float velocity = this->getFlowSpeedAtHydroNode(this->hydroNodes[node.nearestHydroNodeID]);
+    return scaledFlowSpeed(velocity, node);
 }
 
 float limitWaterTemp(float waterTemp, HabitatType nodeType) {
diff --git a/src/hydro.h b/src/hydro.h
index 63b71f8..69a4a5c 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -34,6 +34,8 @@ class HydroModel {
     float getFlowSpeedAlong(Edge &edge);
     // Return the flow speed in m/s at a given location
     float getFlowSpeedAt(MapNode &node);
+    float getFlowSpeedAtHydroNode(DistribHydroNode &hydroNode);
+    float scaledFlowSpeed(float velocity, const MapNode &node);
     // Return the temperature in degrees C at a given location
     float getTemp(MapNode &node);
     // Return the water depth in meters at a given location
diff --git a/src/map.cpp b/src/map.cpp
index 85263e4..c9bebff 100644
--- a/src/map.cpp
+++ b/src/map.cpp
@@ -32,6 +32,15 @@ bool isNearshore(HabitatType t) {
     }
 }
 
+bool isBlindChannel(HabitatType t) {
+    return t == HabitatType::BlindChannel;
+}
+
+bool isImpoundment(HabitatType t) {
+    return t == HabitatType::Impoundment;
+}
+
+
 float habTypeMortalityConst(HabitatType t) {
     switch(t){
     case HabitatType::Distributary:
diff --git a/src/map.h b/src/map.h
index 6dbf7ce..3b4c3eb 100644
--- a/src/map.h
+++ b/src/map.h
@@ -45,6 +45,10 @@ bool isDistributaryOrHarbor(HabitatType t);
 // Returns true if t == Distributary or t == isNearshore
 bool isNearshore(HabitatType t);
 
+bool isBlindChannel(HabitatType t);
+
+bool isImpoundment(HabitatType t);
+
 // Returns the mortality constant for a given habitat type
 float habTypeMortalityConst(HabitatType t);