update doc

docs/tech_guide/dep.rst

@@ -42,7 +42,7 @@ Once the volumes of these depressions have been obtained, their subsequent filli
 
   Illustration of the two cases that may arise depending on the volume of sediment entering an internally drained depression (panel **a**). The red line shows the limit of the depression at the minimal spillover elevation. **b)** The volume of sediment (:math:`\mathrm{V_s^{in}}`) is lower than the depression volume :math:`\mathrm{V_{pit}}`. In this case all sediments are deposited and no additional calculation is required. **c)** If :math:`\mathrm{V_s^{in}\ge V_{pit}}`, the depression is filled up to depression filling elevation (priority-flood + :math:`\mathrm{\epsilon}`), the flow calculation needs to be recalculated and the excess sediment flux (:math:`\mathrm{Q_s^{ex}}`) is transported to downstream nodes.
 
-In cases where the incoming sediment volume is lower than the depression volume (Fig. 4b), all sediments are deposited and the elevation at node $i$ in the depression is increased by a thickness :math:`\mathrm{\delta_i}` such that:
+In cases where the incoming sediment volume is lower than the depression volume (Fig. 4b), all sediments are deposited and the elevation at node :math:`i` in the depression is increased by a thickness :math:`\mathrm{\delta_i}` such that:
 
 .. math::
 

docs/tech_guide/ero.rst

@@ -198,24 +198,24 @@ The approach considers the local balance between erosion and deposition and is b
 
 .. math::
 
-   \mathrm{Q_s^{out}} = (1.0 - G \Omega / Q) \mathrm{Q_s^{in} + E \Omega}
+	\mathrm{\frac{\eta_i^{t+\Delta t}-\eta_i^t}{\Delta t}} =  \mathrm{-\kappa P^d_i \sqrt{Q_i} \frac{\eta_i^{t+\Delta t} - \eta_{rcv}^{t+\Delta t}}{\lambda_{i,rcv}}} + \mathrm{G' Q_{s_i} / \Omega_i}
 
-This equation is solved similarly to the one detailed above using the same matrix system where the weighted in sparse matrix are multiplied by :math:`G \Omega / Q` based on local cell areas and water fluxes.
+where :math:`\mathrm{\lambda_{i,rcv}}` is the length of the edges connecting the considered vertex to its receiver and :math:`\mathrm{\Omega_i}` is the area (voronoi) of the node :math:`i`. 
 
-.. note::
-
-  As our previous technique, this implicit method is unconditionally stable ans allows to study the dynamics of fluvial systems including the transition from detachment-limited to transport-limited behavior.
+:math:`\mathrm{Q_{s_i}}` is the upstream incoming sediment flux in m3/yr and :math:`\mathrm{G'}` is equal to :math:`\mathrm{G \Omega_i / \bar{P}A}`.
 
-In turn, the deposition flux (:math:`Q_{d}`) at any point in the mesh is given by:
+The upstream incoming sediment flux is obtained from the total sediment flux :math:`\mathrm{Q_{t_i}}` where:
 
 .. math::
 
-   \mathrm{Q_d} = \Upsilon \mathrm{Q_s^{\star}}
+	\mathrm{Q_{t_i}^{t+\Delta t} - \sum_{ups} w_{i,j} Q_{t_u}^{t+\Delta t}}= \mathrm{(\eta_i^{t} - \eta_i^{t+\Delta t}) \frac{\Delta t}{\Omega_i}}
 
-where :math:`\Upsilon` is expressed as: 
+which gives:
 
 .. math::
 
-   \Upsilon = \frac{G \Omega / Q}{1.0 - G \Omega / Q} 
+	\mathrm{Q_{s_i}} = \mathrm{Q_{t_i}} - \mathrm{(\eta_i^{t} - \eta_i^{t+\Delta t}) \frac{\Delta t}{\Omega_i}}
+
+This system of coupled equations is solved implicitly using PETSc by assembling the matrix and vectors using the nested submatrix and subvectors and by using the ``fieldsplit`` preconditioner combining two separate preconditioners for the collections of variables. 
 
-and local :math:`\mathrm{Q_s^{\star}}` equals :math:`\mathrm{Q_s^{out}} - E \Omega`.
+The ``TFQMR`` (transpose-free QMR (quasi minimal residual)) KSP solver is used to solve the coupled system with sub KSPs set to ``preonly`` and preconditioner set to ``hypre``. (See PETSC documentation for more details about the solver and preconditoner options and settings).  

docs/tech_guide/hill.rst

@@ -9,7 +9,7 @@ Hillslope: soil creep
 
 Hillslope processes are known to strongly influence catchment morphology and drainage density and several formulations of hillslope transport laws have been proposed. Most of these formulations are based on a mass conservation equation and assume that a layer of soil available for transport is always present (*i.e.* precluding  case of bare exposed bedrock) and that dissolution and mass transport in solution can be neglected.
 
-Under such assumptions and via the Exner's law, the mass conservation equation widely applied in landscape modelling is of the form:
+Under such assumptions, the mass conservation equation widely applied in landscape modelling is of the form:
 
 .. math::
 
@@ -48,20 +48,16 @@ In goSPL, these parameters remain fixed  during a model run and therefore :math:
 Marine deposition
 --------------------
 
-In the marine realm, a nonlinear diffusion model is used for sediment-transport by rivers. When the dual lithology is activated, :mod:`gospl`  accounts for distinct transport coefficients for the two different grain sizes (`Rivenaes, 1992 <https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1365-2117.1992.tb00136.x>`_).
-
-Sediment transport is modelled through nonlinear diffusion equation. The rate of elevation change in the marine environment is governed by:
+In the marine realm, sediment transport is modelled through nonlinear diffusion equation. The rate of elevation change in the marine environment is governed by:
 
 
 .. math::
 
-  \mathrm{\frac{\partial \eta}{\partial t}} = \mathrm{\nabla \cdot \left( K_M(\eta) \nabla \eta \right)} + Q_{sr}
+  \mathrm{\frac{\partial \eta}{\partial t}} = \mathrm{\nabla \cdot \left( K_m(\eta) \nabla \eta \right)} + Q_{sr}
   
 
-where :math:`\mathrm{K_M}` is the marine sediment transport coefficient (m2/yr), and :math:`\mathrm{Q_{sr}}` is the sediment flux coming at the river mouth. As the model progresses over time so does the shoreline position due to both offshore sedimentation and prescribed eustatic sea-level variations.
-
-As already mentioned, the distinct transport efficiency of different grain sizes is included in our algorithm for marine sediment transport and deposition by distinct transport coefficients, :math:`\mathrm{K_{M1}}` and :math:`\mathrm{K_{M2}}` for coarse and fine sediments, respectively. goSPL considers that these transport coefficients are uniform in space and in time. When using equation above, :math:`\mathrm{Q_{sr}}` is divided in two components :math:`\mathrm{Q_{sr1}}` and :math:`\mathrm{Q_{sr2}}` which are the fully uncompacted flux of coarse and fine coming from the continental domain of the model, available for transport and deposition.
+where :math:`\mathrm{K_m}` is the marine sediment transport coefficient (m2/yr), and :math:`\mathrm{Q_{sr}}` is the sediment flux coming at the river mouth. As the model progresses over time so does the shoreline position due to both offshore sedimentation and prescribed eustatic sea-level variations.
 
-During a single time step, marine sediment transport is performed until all available sediment transported by rivers to the ocean have been diffused and the accumulations on these specific nodes remains below water depth or below a prescribed slope computed based on local water depth and distance to the nearest coastline.
+During a single time step, marine sediment transport is performed until all available sediment transported by rivers to the ocean have been diffused and the accumulations on these specific nodes remains below water depth or below a prescribed slope (i.e., clinoform slope) computed based on local water depth and distance to the nearest coastline.
 
-Like for inland deposition, the coarser sediments are deposited first, followed by the finer ones. It allows for finer sediments to be deposited further and reproduce the standard behaviour observed in stratigraphic architectures. The implicit finite volume approach already presented above is implemented and the solution for elevation changes are obtained in a similar fashion as for the solution of the other matrix systems using `PETSc <https://www.mcs.anl.gov/petsc/>`_ *Richardson solver* and *block Jacobi* preconditioning.
+The implicit finite volume approach is implemented using `PETSc <https://www.mcs.anl.gov/petsc/>`_ SNES functionality. The nonlinear system at each time step (``arkimex``) is solved iteratively with time stepping and the SNES solution is based on a Nonlinear Generalized Minimum Residual method (``NGMRES``) and the linear solver uses a ``richardson`` KSP with a ``bjacobi`` preconditioner.  (See PETSC documentation for more details about the solver and preconditoner options and settings). 

docs/tech_guide/strat.rst

@@ -46,7 +46,7 @@ where :math:`\mathrm{\phi_0}` is the surface porosity of sediments, and :math:`\
 
 As shown in the figure above, porosity decreases for burial depth of the order of a few thousand meters (:math:`\mathrm{z_0}`), accordingly associated compaction increases substantially (`Sclater and Christie, 1980 <https://agupubs.onlinelibrary.wiley.com/doi/10.1029/JB085iB07p03711>`_). We can also see that the porosity and rate of compaction between fine and coarse sediments are significantly different. As a result, goSPL uses the proportion of coarse versus fine at all depths within the underlying stratigraphy column to properly estimate compaction and induced elevation changes.
 
-For a given stratigraphic layer :math:`\mathrm{i}`, the associated porosity fis obtained at the centre of the layer for a specific depth :math:`\mathrm{\bar{z}_i}` by:
+For a given stratigraphic layer :math:`\mathrm{i}`, the associated porosity is obtained at the centre of the layer for a specific depth :math:`\mathrm{\bar{z}_i}` by:
 
 .. math::
 

docs/tech_guide/tecto.rst

@@ -60,7 +60,7 @@ Flexural isostasy
 
 
 .. figure:: ../images/flex.png
-  :scale: 40 %
+  :scale: 50 %
   :align: center
 
   Flexural isostasy can be produced in response to a range of geological loads (from `Wickert, 2016 <https://gmd.copernicus.org/articles/9/997/2016/gmd-9-997-2016.pdf>`_).

gospl/flow/flowplex.py

@@ -618,13 +618,37 @@ def _eroMats(self, hOldArray):
         return eMat, PA
 
     def _coupledEDSystem(self, eMat):
-        """
+        r"""
         Setup matrix for the coupled linear system in which the SPL model takes into account sediment deposition.
 
         .. note::
 
             The approach follows `Yuan et al, 2019 <https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018JF004867>`_, where the deposition flux depends on a deposition coefficient :math:`G` and is proportional to the ratio between cell area :math:`A` and flow accumulation :math:`FA`.
 
+        The approach considers the local balance between erosion and deposition and is based on sediment flux resulting from net upstream erosion.
+
+        .. math::
+
+            \mathrm{\frac{\eta_i^{t+\Delta t}-\eta_i^t}{\Delta t}} =  \mathrm{-\kappa P^d_i \sqrt{Q_i} \frac{\eta_i^{t+\Delta t} - \eta_{rcv}^{t+\Delta t}}{\lambda_{i,rcv}}} + \mathrm{G' Q_{s_i} / \Omega_i}
+
+        where :math:`\mathrm{\lambda_{i,rcv}}` is the length of the edges connecting the considered vertex to its receiver and :math:`\mathrm{\Omega_i}` is the area (voronoi) of the node :math:`i`.
+
+        :math:`\mathrm{Q_{s_i}}` is the upstream incoming sediment flux in m3/yr and :math:`\mathrm{G'}` is equal to :math:`\mathrm{G \Omega_i / \bar{P}A}`.
+
+        The upstream incoming sediment flux is obtained from the total sediment flux :math:`\mathrm{Q_{t_i}}` where:
+
+        .. math::
+
+            \mathrm{Q_{t_i}^{t+\Delta t} - \sum_{ups} w_{i,j} Q_{t_u}^{t+\Delta t}}= \mathrm{(\eta_i^{t} - \eta_i^{t+\Delta t}) \frac{\Delta t}{\Omega_i}}
+
+        which gives:
+
+        .. math::
+
+            \mathrm{Q_{s_i}} = \mathrm{Q_{t_i}} - \mathrm{(\eta_i^{t} - \eta_i^{t+\Delta t}) \frac{\Delta t}{\Omega_i}}
+
+        This system of coupled equations is solved implicitly using PETSc by assembling the matrix and vectors using the nested submatrix and subvectors and by using the ``fieldsplit`` preconditioner combining two separate preconditioners for the collections of variables.
+
         :arg eMat: erosion matrix (from the simple SPL model)
         """
 
@@ -724,13 +748,7 @@ def _getEroDepRate(self):
 
         The erosion rate is solved by an implicit time integration method, the matrix system is based on the receiver distributions and is assembled from local Compressed Sparse Row (**CSR**) matrices into a global PETSc matrix. The PETSc *scalable linear equations solvers* (**KSP**) is used with both an iterative method and a preconditioner and erosion rate solution is obtained using PETSc Richardson solver (`richardson`) with block Jacobian preconditioning (`bjacobi`).
 
-        Once the erosion rate solution has been obtained, local sediment flux depends on upstream fluxes, local eroded flux and local deposition flux. We assume that local deposition depends on user-defined forced deposition value (:math:`fDep`), the local water flux and the cell area. Here again the sediment flux is determined implicitly and corresponding deposition flux are calculated subsequently once the local total flux is known.
-
-        .. math::
-
-          Qs_{in} = Qs_e + (1 - GA/FA) \sum_{j \in upstream} w_{j,i} Qs_{in_j}
-
-        where the incoming sediment flux adds the local erosion flux obtained from the stream power law equation and the upstream sediment flux minus the deposition flux. Following `Yuan et al, 2019 <https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018JF004867>`_, the deposition flux depends on a deposition coefficient :math:`G` and is proportional to the ratio between cell area :math:`A` and flow accumulation :math:`FA`.
+        An alternative method to the detachment-limited approach proposed above consists in accounting for the role played by sediment in modulating erosion and deposition rates. It follows the model of `Yuan et al, 2019 <https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018JF004867>`_, whereby the deposition flux depends on a deposition coefficient :math:`G` and is proportional to the ratio between cell area :math:`\mathrm{\Omega}` and water discharge :math:`\mathrm{Q}=\bar{P}A`.
         """
 
         t0 = process_time()

gospl/sed/seaplex.py

@@ -203,7 +203,9 @@ def _matOcean(self):
 
     def _evalFunction(self, ts, t, x, xdot, f):
         """
-        The nonlinear system at each time step is solved iteratively using PETSc time stepping and SNES method and is based on a Newton/GMRES method. Here we define the function for the nonlinear solve.
+        The nonlinear system at each time step is solved iteratively using PETSc time stepping and SNES solution and is based on a Nonlinear Generalized Minimum Residual method (``NGMRES``) .
+        
+        Here we define the function for the nonlinear solve.
         """
 
         self.dm.globalToLocal(x, self.hl)
@@ -218,7 +220,9 @@ def _evalFunction(self, ts, t, x, xdot, f):
 
     def _evalJacobian(self, ts, t, x, xdot, a, J, P):
         """
-        The nonlinear system at each time step is solved iteratively using PETSc time stepping and SNES method and is based on a Newton/GMRES method. Here we define the Jacobian for the nonlinear solve.
+        The nonlinear system at each time step is solved iteratively using PETSc time stepping and SNES solution and is based on a Nonlinear Generalized Minimum Residual method (``NGMRES``) .
+         
+        Here we define the Jacobian for the nonlinear solve.
         """
 
         self.dm.globalToLocal(x, self.hl)
@@ -260,7 +264,6 @@ def _diffuseOcean(self, dh):
 
             PETSc SNES and time stepping TS approaches are used to solve the nonlinear equation above over the considered time step.
 
-        The nonlinear diffusion equation is ran for the coarser sediment first and for the finest ones afterwards. This mimicks the standard behaviour observed in stratigraphic architectures where the fine fraction are generally transported over longer distances.
 
         :arg dh: numpy array of incoming marine depositional thicknesses