update input

docs/user_guide/inputfile.rst

@@ -39,7 +39,7 @@ Initial mesh definition and simulation declaration
         a. the initial spherical or 2D surface mesh ``npdata`` (**.npz** file). This file contains the following keys: ``v`` the mesh vertices coordinates, ``z`` the vertice elevations, ``c`` the mesh cells.
         b. the flow direction method to be used ``flowdir`` that takes an integer value between 1 (for SFD) and 6 (for MFD)
         c. the exponent (``flowexp``) used in the flow direction approach. Default value is set to 0.42.
-        d. boundary conditions (``bc``) when not running a global model. Each integer corresponds to an edge defined in the following order: south, east, north, and west. The integer is set to either 0 for open or to 1 for closed boundaries.
+        d. boundary conditions (``bc``) when not running a global model. Each integer corresponds to an edge defined in the following order: south, east, north, and west. The integer is set to either 0 for open or to 1 for fixed boundaries.
 
         In addition the following optional parameters could be set:
 
@@ -54,7 +54,7 @@ Initial mesh definition and simulation declaration
 
 .. note::
 
-  It is also possible to have only one **.npz** file containing the required keys. For example you could set an input file containing the following keys: ``v`` the mesh vertices coordinates, ``z`` the vertice elevations, ``c`` the mesh cells, ``ed`` erosion deposition and call it for both the ``npdata`` and ``nperodep`` parameters.
+  It is also possible to have only one **.npz** file containing the required keys. For example you could set an input file containing the following keys: ``v`` the mesh vertices coordinates, ``z`` the vertice elevations, ``c`` the mesh cells, and ``ed`` the erosion deposition node values and call it for both the ``npdata`` and ``nperodep`` parameters.
 
 
 Setting model temporal evolution

docs/user_guide/optfile1.rst

@@ -8,6 +8,10 @@ Climate related parameters
 Sea-level (eustatic) forcing
 -----------------------------
 
+.. note::
+
+    By default, the sea-level position in goSPL is set to 0 m. If you wish to set it to another position you can use the ``position`` key that changes the sea-level to a new value relative to sea-level. Another option consists in defining your own sea-level curve (``curve``) or using a published one (e.g. Haq curve for example). 
+
 .. grid:: 1
     :padding: 3
 
@@ -24,9 +28,12 @@ Sea-level (eustatic) forcing
 
         The sea-level declaration is defined with 2 optional parameters:
 
-        a. the relative sea-level ``position`` in meters (optional),
-        b. a sea-level ``curve`` *e.g.* a file containing 2 columns (time and sea-level position).
+        a. the relative sea-level ``position`` in meters (optional), its default value is set to 0.0
+        b. a sea-level ``curve`` *e.g.* a file containing 2 columns (time and sea-level position). Not required in case no sea-level fluctuations needs to be specified. 
+
+.. important::
 
+    The sea-level curve is defined as a 2 columns **csv** file containing in the first column the time in years (it doesn't need to be regularly temporally spaced) and in the second the sea-level position for the given time. When goSPL interprets this file, it will interpolate linearly between the defined times to find the position of the sea-level for every time step.
 
 Climatic (rainfall) forcing conditions
 ----------------------------------------
@@ -42,15 +49,31 @@ Climatic (rainfall) forcing conditions
 
             climate:
                 - start: -20000000.
-                map: ['input8/rain20Ma','r']
+                  map: ['input8/rain20Ma','r']
                 - start: -15000000.
-                uniform: 1.
+                  uniform: 1.
+
+
+        The climatic forcing is defined in a similar fashion as the tectonic one. It is based on a sequence of events with each event starting at a given time (``start`` in years) and corresponding to a given precipitation condition. This could either be an uniform rainfall over the entire mesh (``uniform``) or a precipitation mesh ``map``. The rainfall values have to be in metres per year and the precipitation is updated at every time step (defined by ``dt``).
+
+.. important::
+
+    When defining a precipitation grid, one needs to use the **npz** format and needs to specify the key corresponding to the precipitation variable in the file. In the above example this key is ``'r'``. The precipitation grid needs to define values for all vertices in the mesh.
+
+
+Orographic precipitation definition
+------------------------------------
+
+goSPL implements the Linear Theory of Orographic Precipitation following `Smith & Barstad (2004) <https://journals.ametsoc.org/view/journals/atsc/61/12/1520-0469_2004_061_1377_altoop_2.0.co_2.xml>`_.
 
+.. note::
+
+    The model includes airflow dynamics, condensed water advection, and downslope evaporation. It consists of two vertically-integrated steady-state advection equations describing: 
 
-        The climatic forcing is defined in a similar fashion as the tectonic one with again a sequence of events by a starting time (``start``) and either an uniform rainfall over the entire mesh (``uniform``) or with a precipitation mesh ``map``. The rainfall values have to be in metres per year.
+    1. the cloud water density and 
+    2. the hydrometeor density. 
 
-Orographic rain definition
----------------------------
+    Solving these equations using Fourier transform techniques, derives a single formula relating terrain and precipitation.
 
 .. grid:: 1
     :padding: 3
@@ -62,6 +85,7 @@ Orographic rain definition
         .. code:: python
 
             orography:
+                regdx: 200.
                 latitude: 40.0  
                 wind_speed: 10.0 
                 wind_dir: 0 
@@ -78,16 +102,17 @@ Orographic rain definition
             
         This part of the input file define the parameters for the orographic rain:
 
-        a. ``latitude``: average latitude used to compute the Coriolis factors [degrees btw -90 and 90]; default 0
-        b. ``wind_speed``: wind speed in m/s; default 10
-        c. ``wind_dir``: wind direction [0: north, 270: west]; default 0
-        d. ``nm``: moist stability frequency [1/s]; default 0.01
-        e. ``env_lapse_rate``: environmental lapse rate [degrees Celsius/km]; default -4.0
-        f. ``moist_lapse_rate``: moist adiabatic lapse rate [degrees Celsius/km]; default -7.0
-        g. ``ref_density``: reference saturation water vapor density [kg/m^3]; default 7.4e-3
-        h. ``hw``:  water vapor scale height [m]; default 3400
-        i. ``conv_time``: cloud water to hydrometeor conversion time [s]; default 1000
-        j. ``fall_time``: hydrometeor fallout time [s]; default 1000
-        k. ``oro_precip_base``: non-orographic, uniform precipitation rate [mm/h]; default 7.
-        l. ``oro_precip_min``: minimum precipitation [mm/h] when precipitation rate <= 0; default 0.01
-        m. ``rainfall_frequency``: number of storm of 1 hour duration per day; default 1
+        a. ``regdx``: the resolution of the regular grid used to perform the orographic rain calculation.
+        b. ``latitude``: average latitude used to compute the Coriolis factors [degrees btw -90 and 90]; default 0
+        c. ``wind_speed``: wind speed in m/s; default 10
+        d. ``wind_dir``: wind direction [0: north, 270: west]; default 0
+        e. ``nm``: moist stability frequency [1/s]; default 0.01
+        f. ``env_lapse_rate``: environmental lapse rate [degrees Celsius/km]; default -4.0
+        g. ``moist_lapse_rate``: moist adiabatic lapse rate [degrees Celsius/km]; default -7.0
+        h. ``ref_density``: reference saturation water vapor density [kg/m^3]; default 7.4e-3
+        i. ``hw``:  water vapor scale height [m]; default 3400
+        j. ``conv_time``: cloud water to hydrometeor conversion time [s]; default 1000
+        k. ``fall_time``: hydrometeor fallout time [s]; default 1000
+        l. ``oro_precip_base``: non-orographic, uniform precipitation rate [mm/h]; default 7.
+        m. ``oro_precip_min``: minimum precipitation [mm/h] when precipitation rate <= 0; default 0.01
+        n. ``rainfall_frequency``: number of storm of 1 hour duration per day; default 1

docs/user_guide/optfile2.rst

@@ -19,32 +19,41 @@ Tectonic forcing parameters
 
             tectonic:
                 - start: -20000000.
-                end: -19000000.
-                mapH: 'input8/disp20Ma'
+                  end: -19000000.
+                  mapH: 'input8/disp20Ma'
                 - start: -19000000.
-                end: -18000000.
-                mapH: 'input8/disp19Ma'
+                  end: -18000000.
+                  mapH: 'input8/disp19Ma'
                 - start: -18000000.
-                end: -17000000.
-                mapH: 'input8/disp18Ma'
+                  end: -17000000.
+                  mapH: 'input8/disp18Ma'
                 - start: -17000000.
-                end: -16000000.
-                mapH: 'input8/disp17Ma'
-                mapV: 'input8/dispv17Ma'
+                  end: -16000000.
+                  mapH: 'input8/disp17Ma'
+                  mapV: 'input8/dispv17Ma'
                 - start: -16000000.
-                end: -15000000.
-                mapV: 'input8/dispv16Ma'
+                  end: -15000000.
+                  mapV: 'input8/dispv16Ma'
 
-        Follows the tectonic forcing conditions with a sequence of events defined by a starting time (``start``) and either a vertical only forcing (*e.g.* uplift and/or subsidence defined with ``mapV``) or a fully 3D displacement mesh ``mapH``. These displacements are set in metres per year.
+        It defines the tectonic forcing conditions from a sequence of events defined by a starting and ending time (``start`` and ``end``) and either a vertical only forcing (*e.g.* uplift and/or subsidence defined with ``mapV``) or a fully 3D displacement mesh ``mapH``. **These displacement rates are set in metres per year**.
 
 
 .. important::
 
-  As mentioned above and for the next key parameter as well, these forcing files are defined as numpy zip array (**.npz**).
+  Here again, these forcing files are defined as numpy zip array (**.npz**). These files use specific keys to identify the tectonic forcing. For vertical only condition, the key ``z`` specified the displacements. For the horizontal condition, the ``xyz`` key is used where `xyz` is a 3D array containing the displacement rates along the x, y and z axis (in m/yr). 
 
-Forcing paleo-topography definition
+.. note::
+
+  There is no requirement to impose continuous tectonics forcing and you might chose to have periods without displacement by making discontinuous events using the ``start`` and ``end`` keys. 
+
+Flexural isostasy definition
 -----------------------------------
 
+This function computes the flexural isostasy equilibrium based on topographic change. It is a simple routine that accounts for flexural isostatic rebound associated with erosional loading/unloading.
+
+.. warning::
+    This function assumes a value of 1011 Pa for Young's modulus, 0.25 for Poisson's ratio and 9.81 m/s2 for g, the gravitational acceleration.
+
 .. grid:: 1
     :padding: 3
 
@@ -54,15 +63,15 @@ Forcing paleo-topography definition
 
         .. code:: python
 
-            forcepaleo:
-                dir: 'output-backward'
-                steps: [5,10,5]
+            flexure: 
+                regdx: 200.
+                thick: 20.e3
+                rhoc: 2800.0
+                rhoa: 3150.0
 
-        For simulations that require to be forced with paleo-topography maps obtained from backward models, you will also have to set this key composed of 2 parameters:
-
-        a. ``dir`` the directory containing the outputs of the backward model,
-        b. ``steps`` the steps from the model outputs that will be used to force the forward model topography.
-
-.. important::
+        Used to consider flexural isostasy in 2D simulation (*i.e.* not global scale).  paleo-topography maps obtained from backward models, you will also have to set this key composed of 2 parameters:
 
-  The ``steps`` often correspond to the time where you have a paleotopography dataset that you want to match for example from a Scotese paleotopography map.
+        a. ``regdx``: the resolution of the regular grid used to perform the flexural isostasy calculation,
+        b. ``thick`` effective elastic plate thickness in m,
+        c. ``rhoc`` crust density in kg/m3,
+        d. ``rhoa`` asthenospheric density in kg/m3.

docs/user_guide/surfproc.rst

@@ -28,14 +28,17 @@ Stream Power Law parameters
 
         a. ``K`` representing the erodibility coefficient which is scale-dependent and its value depend on lithology and mean precipitation rate, channel width, flood frequency, channel hydraulics. It is used in the SPL law: :math:`E = K (\bar{P}A)^m S^n`
 
-        .. warning::
-            It is worth noting that the coefficient *n* is fixed and take the value *1*.
+        The following parameters are **optional**:
 
-        b. Studies have shown that the physical strength of bedrock which varies with the degree of chemical weathering, increases systematically with local rainfall rate. Following `Murphy et al. (2016) <https://doi.org/10.1038/nature17449>`_, the stream power equation is adapted to explicitly incorporate the effect of local mean annual precipitation rate, P, on erodibility: :math:`E = (K_i P^d) (\bar{P}A)^m S^n`. ``d`` (:math:`d` in the equation) is a positive exponent that has been estimated from field-based relationships to 0.42. Its default value is set to 0.
+        b. Studies have shown that the physical strength of bedrock which varies with the degree of chemical weathering, increases systematically with local rainfall rate. Following `Murphy et al. (2016) <https://doi.org/10.1038/nature17449>`_, the stream power equation could be adapted to explicitly incorporate the effect of local mean annual precipitation rate, P, on erodibility: :math:`E = (K_i P^d) (\bar{P}A)^m S^n`. ``d`` (:math:`d` in the equation) is a positive exponent that has been estimated from field-based relationships to 0.42. Its default value is set to 0.0
         c. ``m`` is the coefficient from the SPL law: :math:`E = K (\bar{P}A)^m S^n` and takes the default value of 0.5.
-        d. ``fDa`` dimensionless deposition coefficient for continental domain
-        e. ``fDm`` dimensionless deposition coefficient for marine domain
-        f. ``mthd`` chosen approach to account for sediment deposition (should be either 1 or 2)
+
+        .. note::
+            It is worth noting that the coefficient *n* in the SPL is fixed and take the value *1*.
+
+        d. ``fDa`` dimensionless deposition coefficient for continental domain when accounting for sedimentation rate in the SPL following the model of `Yuan et al, 2019 <https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018JF004867>`_. The default value is 0.0 (purely detachment-limited model).
+        e. ``fDm`` dimensionless deposition coefficient for marine domain. The default value is 40.0.
+        f. ``mthd`` chosen approach to account for sediment deposition (should be either 1 or 2). The default value is set to 1 (other choice is 2). While the first method uses the approach from `Yuan et al, 2019 <https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018JF004867>`_, the second provides a faster calculation but might not be conservative.
 
 Hillslope and marine deposition parameters
 -------------------------------------------
@@ -60,9 +63,8 @@ Hillslope and marine deposition parameters
         a. ``hillslopeKa`` is the diffusion coefficient for the aerial domain,
         b. ``hillslopeKm`` is the diffusion coefficient for the marine domain,
         c. ``smthDep`` is the transport coefficient of freshly deposited sediments entering the ocean from rivers,
-        d. ``clinSlp`` is the maximum slope of clinoforms (needs to be positive), this slope is then used to estimate the top of the marine deposition based on distance to shore.        
-
-
+        d. ``clinSlp`` is the maximum slope of clinoforms (needs to be positive), this slope is then used to estimate the top of the marine deposition based on distance to shore. 
+        d. ``diffNb`` is the number of steps used to distribute the sediment fluxes in the marine domain. Default value is set to 1.        
 
 Compaction & porosity variables definition
 ------------------------------------------
@@ -82,6 +84,6 @@ Compaction & porosity variables definition
 
         We assume  a depth-porosity relationship for the sediment compaction based on the following parameters:
 
-        a. porosity at the surface ``phis``,
-        b. e-folding depth ``z0s`` (in metres)
+        a. porosity at the surface ``phis``, default value is set to 0.49,       
+        b. e-folding depth ``z0s`` (in metres), default value is set to 3700.       
 

gospl/tools/inputparser.py

@@ -372,7 +372,7 @@ def _readSPL(self):
             try:
                 self.fDepm = splDict["fDm"]
             except KeyError:
-                self.fDepm = 0.0
+                self.fDepm = 40.0
             try:
                 self.spl_m = splDict["m"]
             except KeyError:
@@ -385,7 +385,7 @@ def _readSPL(self):
             self.K = 1.0e-12
             self.coeffd = 0.0
             self.fDepa = 0.0
-            self.fDepm = 30.0
+            self.fDepm = 40.0
             self.spl_m = 0.5
             self.dmthd = 1