Merge branch 'ESCOMP:main' into v3.0.1-26-g85bbe71f_tracer

Author: nmizukami

docs/source/tech_note/images/overall_comp_workflow.png

@@ -1,7 +1,7 @@
 Lake water balance computation
 ================================
 
-mizuRoute include one method-Doll-to compute natural lake discharge, and two methods-Hanasaki and HYPE- to compute discharge for managed lake.
+MizuRoute include one method-Doll-to compute natural lake discharge, and two methods-Hanasaki and HYPE- to compute discharge for managed lake.
 
 Doll model
 --------------------------

docs/source/tech_note/images/uh_convolution.png

@@ -4,16 +4,66 @@
 Overall workflow
 ======================
 
-#. Compute HRU (Hydrologic Response Unit, or simply catchment) mean runoff [m/s]
+Overall computation workflow is shown in :numref:`Figure_overall_comp_workflow`. Starting with runoff depth from netCDF or coupler (e.g. CTSM coupling),
 
-#. Perform overland routing with HRU mean runoff as an input to compute lateral runoff (:math:`R_{lat}`) [m/s]
+#. Remap runoff depth [m/s] to river network HRU (Hydrologic Response Unit or simply catchment), :math:`R_{lat}`, if runoff is given at hydrologic model HRU
 
-#. Convert :math:`R_{lat}` from depth unit to volume (lateral runoff times HRU area) to get lateral discharge (:math:`q_{lat}`) [m3/s]
+#. Convert :math:`R_{lat}` from depth unit to volume (:math:`R_{lat}` times HRU area) to get lateral runoff volume (:math:`q_{lat}`) [m3/s]
 
-#. Route lateral discharge at each river reach through river network.
+#. Perform hillslope routing to delay lateral runoff volume, if travel time of runoff is not counted outside mizuRoute.
 
-The overland routing is optional, and currently simple unit hydrograph based on gamma distribution is used to delay instantaneous runoff.
+#. Route inflow from upstream and add delayed lateral discharge to routed inflow at each river reach outlet.
 
+The hillslope routing method currently uses a simple unit hydrograph based on gamma distribution (only one method available) to delay instantaneous runoff.
+See :ref:`Hillslope routing scheme <Hillslope_routing_scheme>`).
+
+The river reach routing needs to be performed in the order of upstream-to-downstream to complete the routing in the entire network.
+This routing order is internally computed in mizuRoute based on network topology information, ``downSegId`` (immediate downstream reach ID)
+
+.. _Figure_overall_comp_workflow:
+
+.. figure:: images/overall_comp_workflow.png
+ :width: 700
+ :height: 400
+
+ Overall routing procedures from runoff input into model to streamflow computation.
+
+.. _Hillslope_routing_scheme:
+
+Hillslope routing scheme
+========================
+
+MizuRoute uses a similar concept to a unit-hydrograph commonly used in engineering hydrology to account for travel time of instantaneous runoff to the river reach (delay and attenuate runoff).
+A unit-hydrograph is defined by a hydrograph (time series of discharge) that is derived from a unit depth of excess rainfall on a drainage area within a specific time period.
+Here, this concept is applied to directly runoff volume over a HRU, instead of rainfall excess over a drainage area.
+Therefore, a unit-hydrograph represents a time series of lateral discharge into a river reach from a HRU derived from a unit volume of runoff.
+Here, a probability density function (PDF) is used as a unit-hydrograph, so that cumulative sum of PDF is 1.0.
+This means that runoff volume at a current time step is just distributed in the future, with the sum of future distributed runoff is equal to the current runoff volume (i.e., volume conserved)
+
+To get actual delayed lateral flow series to the river reach, unit-hydrograph convolution is performed as below:
+
+.. _Figure_uh_convolution:
+
+.. figure:: images/uh_convolution.png
+ :width: 700
+ :height: 500
+
+ Illustration of discrete Unit hydrograph convolution.
+
+In mizuRoute, gamma distribution is used for PDF-based unit-hydrograph and written as:
+
+.. math::
+   :label: gamma_distribution
+
+   f(t; a, \theta) = \frac{1}{\Gamma(a)\theta^{a}}t^{a - 1} e^{-\frac{t}{\theta}},
+   \quad t > 0
+
+where *t* is time [sec], *a* is a shape parameter [–] (a > 0), and :math:`\theta` is a timescale parameter [sec].
+Both the shape and timescale parameters affect the peak time of (mode of the distribution: :math:`(a - 1)\theta` and flashiness (variance of the distribution: :math:`a\theta^2` of the unit-hydrograph (UH).
+These UH should depend on the physical HRU characteristics. Continuous gamm distribution is descritized per time step before the convolution performed.
+These shape and scale parameters are currently specified as a spatially constant parameter (see :ref:`Spatially-constant parameter namelist <namelist_file>`), though they could be provided as spatially distributed parameters from river data netCDF (potential future implementation)
+
+Also, please see section 3.1 in :ref:`Mizukami et al. 2016 <Mizukami2016>` for hillslope routing theory.
 
 .. _River_routing_schemes:
 
@@ -96,14 +146,37 @@ The numerical implementations of Euler kinematic wave and Diffusive wave are ess
 Impulse response function
 --------------------------
 
+The impulse response function (IRF) is derived from 1-D diffusive wave equation (Eq. :eq:`0.5`). Eq. :eq:`0.5` can be solved using convolution integrals:
+
+.. math::
+   :label: 0.8
+
+   Q(x,t) = \int_{0}^t Q_{in}(t-s)h(x,s)ds
+
+where :math:`Q_{in}` is a pulse or here input runoff volume in the river reach (e.g., inflow from HRU, or/and upstream reach), and :math:`h(x,s)` is the green function or *impulse response function*, and *x* is the distance from the inflow point in the reach to the reach location where discharge desires to be computed (i.e., HRU outlet).
 
+.. math::
+   :label: 0.9
+
+   h(x,t) = \frac{x}{2t \sqrt{\pi t D}} \exp(-\frac{(Ct-x)^2}{4Dt})
+
+where *C* (Eq. :eq:`0.6`) is a wave celerity [m/s] and *D* (Eq. :eq:`0.7`) is a diffusivity [m\ :sup:`2`\/s].
+
+This is the same concept as river routing described in Lohmann et al., (1996). Also, please see section 3.2.2 in :ref:`Mizukami et al. 2016 <Mizukami2016>`.
+
+Lohmann et al., (1996) and original mizuRoute (version 1) use the *source-to-sink* approach where unique IRFs are developed for each HRU outlet (sink) to its all upstream HRUs (sources), the IRFs are applied with lateral flows from the corresponding HRUs,
+then routed flow from all the source HRUs are summed up to the river discharge at the target HRU outlet. This process is repeated for all the HRUs and each time step.
+
+One modification made from Lohmann et al., 1996 and mizuRoute version 1 was to develop reach-specific IRF based on reach length, *C* and *D*, then apply it with inflow from the immediate upstream reach(es) to move water from upstream end to the outlet.
+The computed river discharge at each river reach is used as inflow of the next downstream reach. Therefore, the IRF computation is performed from upstream HRUs to downstream in the order.
+This approach is consistent with the other routing methods used in mizuRoute, and also requires less memory usage of IRF ordinates than source-to-sink approach, which requires storing IRFs of all the upstream HRUs for all the HRUs.
 
 .. _Lagrangian_kinematic_wave:
 
 Lagrangian kinmatic wave
 --------------------------
 
-
+Also, please see section 3.2.1 in Mizukami et al. 2016
 
 .. _Euler_kinematic_wave:
 

docs/source/tech_note/lake_routing.rst

@@ -97,7 +97,7 @@ Here assume that all the libaries are installed using homebrew 🍺, and gnu com
 
      Execute::
 
-         mk -p libraries/parallelio
+         mkdir -p libraries/parallelio
          ./bin/git-fleximod -g .gitmodules update
 
      See mizuRoute/README_EXTERNALS.rst for more details.