Explimpl eq (#9)

* add equations

* add explicit and implicit in progress

* finished rill implicit description

reference_manual/config/usepackage.tex

@@ -142,6 +142,7 @@
 
 
 
+\usepackage{breqn}
 
 %%%% upraveni sekce
 % \titlespacing*{\section}

reference_manual/text_en/hillhyd.tex

@@ -110,7 +110,7 @@
 
         \begin{equation} 
         q = n^{-1} s^Y h^b,
-        \label{eq:powerlaw}
+        \label{eq:powerlaw2}
         \end{equation}
         where n is roughness,  Y  empirical parameter and s is the surface
         slope ($\mathrm{L.L^{-1}}$).
@@ -264,7 +264,7 @@
             %
             %P stands for the size of the raster cell. 
 
-        \subsubsection{Rill hydraulics (JJ)}
+        \subsubsection{Rill hydraulics (JJ)}\label{sec:rill}
 
             The rill is simplified as a small channel at the soil surface with
             a rectangular cross section. The rectangle has width $b_{rill}$ and
@@ -449,6 +449,7 @@
         \subsubsection{Kinematic}
         \subsubsection{Diffuse}
 
+    \FloatBarrier
     \subsection{Implicit /Explicit approach}
 
         \subsubsection{Explicit}
@@ -464,11 +465,155 @@
             \ref{eq:bilance} which incorporates the sum \ref{eq:d8} is calculated with the
             explicit time discretisation scheme for cell i as:
             \begin{equation} 
-            h_{i,t} =h_{i,t} + \mathrm{d}t (es_{i,t-1} + \sum_j^m q^{out}_{j,t-1}-
+            h_{i,t} =h_{i,t-1} + \mathrm{d}t (es_{i,t-1} + \sum_j^m q^{out}_{j,t-1}-
             inf_{i,t-1} - q^{out}_{i,t-1}),
-            \label{eq:bilance}
+            \label{eq:bilanceexpl}
             \end{equation}
+
+
+           After inserting the equations~(\ref{eq:infiltration})
+           and~(\ref{eq:powerlaw2}) into equation~(\ref{eq:bilanceexpl}) the
+           final explicit form balance equation can be written as:
+            \begin{dmath}
+              h_{i,t} = h_{i,t-1} +
+              es_{i,t-1}  + exf_{i,t-1} +  \sum_{j}^{inflows}
+              1/n_{sheet,j}I_j^{y_j} h_{j,t-1}^{b_j} -
+              1/n_{sheet,i}I_i^{y_i}h_{i,t-1}^{b_i} - ( 
+              \frac{1}{2}S_it_1^{-1/2}+K_{s,i}).
+            \end{dmath}
+
+
         \subsubsection{Implicit}
+        Alternatively the equation (\ref{eq:bilance}) solved using implicit
+        scheme. The right hand side can be expressed in for time $t$ as:
+
+        $$
+          \frac{h_{i,t} - h_{i,t-1} }{\mathrm{d} t} =  \left(es_{i,t} +
+          \sum_j^{inflows} q^{in}_{j,t} - inf_{i,t} - q^{out}_{i,t}\right).
+        $$
+
+        After inserting the power-law equation (\ref{eq:powerlaw}) and several rearrangements the formula above can be written 
+        as follows: 
+
+
+
+        $$
+          h_{i,t} + \mathrm{d} t a_ih^{b_{i}-1}_{i,t} h_{i,t} - \mathrm{d} t \sum_j^{inflows}
+          a_jh^{b_{j}-1}_{j,t} h_{j,t} = h_{i,t-1} + \mathrm{d} t es_{i,t} - \mathrm{d}
+          t inf_{i,t}.
+        $$
+        Here the known part of the balance equation (rainfall and infiltration) are at the right hand side. To 
+        extract the unknown $h$ form the power law a following rearrangement was used:
+
+
+
+        \begin{equation}
+          (1+\mathrm{d} t a_ih^{b_{i}-1}_{i,t})h_{i,t} - \mathrm{d} t \sum_j^{inflows}
+          a_jh^{b_{j}-1}_{j,t} h_{j,t} = h_{i,t-1} + \mathrm{d} t es_{i,t} -
+          \mathrm{d} t inf_{i,t}
+          \label{eq:impl}
+        \end{equation}
+
+        From the equation (\ref{eq:impl}) the system of linear equations can be constructed. 
+
+        As example, lets show the equation (\ref{eq:impl}) for
+        $
+        i=5\ \mathrm{and}\ inflows\in\{7,8,9\}
+        $:
+
+
+        \begin{dmath}
+          (1+\Delta T a_5h^{b_{5}-1}_{5,t})h_{5,t} - \Delta T a_7h^{b_{7}-1}_{7,t} h_{7,t} + \Delta T a_8h^{b_{8}-1}_{8,t} h_{8,t} + \Delta T a_9 h^{b_{9}-1}_{9,t} h_{9,t}  = h_{5,t-1} + \Delta T es_{5,t} - \Delta T inf_{5,t}.
+        \end{dmath}
+
+        In matrix form the above equation can be written as:
+
+
+
+        \begin{dmath}
+        \begin{bmatrix}
+        \hdots& & & & & &  \\
+        \hdots  & 1+\Delta T a_5h^{b_{5}-1}_{5,t} & \hdots &  \Delta T a_7h^{b_{7}-1}_{7,t} &  \Delta T a_8h^{b_{8}-1}_{8,t} & \Delta T a_9 h^{b_{9}-1}_{9,t} & \hdots \\
+        \hdots&  & & & & &  \\
+        \hdots& && &  &  & \\
+        \hdots& && &  &   & \\
+        \hdots& && &  &   & \\
+
+        \end{bmatrix} 
+        % 
+        \begin{bmatrix}
+         \vdots \\
+         h_{5,t} \\
+         \vdots \\
+         h_{7,t} \\
+         h_{8,t} \\
+         h_{9,t} \\
+         \vdots \\
+        %  
+        \end{bmatrix}
+        =
+        \begin{bmatrix}
+         \vdots \\
+        h_{5,t-1} + \Delta T es_{5,t} - \Delta T inf_{5,t} \\
+         \vdots \\
+         \vdots \\
+         \vdots \\
+         \vdots \\
+         \vdots \\
+        %  
+        \end{bmatrix}
+        \label{eq:sys}
+        \end{dmath}
+
+        The system (\ref{eq:sys}) is solved using \texttt{scipy} package method \texttt{root}
+        using \texttt{krigov} methof that shown the best convergence.
+
+
+        \paragraph{Implicit solution with rill flow} To calculated rill flow the water level in rill and needs to be calculated as 
+        \begin{equation}
+        h_{rill} = \text{max}(h-h_{crit},0),
+            \label{eq:hrill2}
+        \end{equation}
+        where 
+        \begin{equation}
+        h_{sheet} = \text{min}(h,h_{crit}).
+            \label{eq:hsheet2}
+        \end{equation}
+
+        Mannings formula is used to calculated the steady state flow in
+        the rill assuming rill is a channel (description in section \ref{sec:rill}.  
+        $$
+        q_{rill} = 1/n R(h_{rill})^{2/3} i^{1/2}
+        $$
+
+
+        The balance equation becomes
+        $$
+          \frac{h_{i,t} - h_{i,t-1} }{\Delta T} = 
+          es_{i,t} + \sum_j^{inflows} a_jh^{b_{j}}_{sh,j,t}  + \sum_k^{inflows} 1/n_k R_k(h_{rl,k,t})^{2/3} i_k^{1/2} - inf_{i,t} - a_ih^{b_{i}}_{sh,i,t} - 1/n R(h_{rl,i,t})^{2/3} i^{1/2},
+        $$
+        where $h_{sheet}$ and $h_{rill}$ are defined as (\ref{eq:hsheet2}) and  (\ref{eq:hrill2}).
+
+        This formula is rearranged and, for practical purposes, the linear equation is the system is calculated with condition \\
+        for  $h<=h_{crit}$ as
+        \begin{equation}
+            \left(\frac{1}{\Delta T}+a_ih^{b_{i}-1}_{i,t}\right)h_{i,t} -  \sum_j^{inflows} a_jh^{b_{j}-1}_{j,t} h_{j,t} = \frac{h_{i,t-1}}{\Delta} +  es_{i,t} - inf_{i,t}
+        \end{equation}
+        % 
+        % 
+        and for  $h>h_{crit}$ as
+        \begin{multline}
+          \left(\frac{1}{\Delta T}
+          + 1/n_i R_i(h_{i,t}-h_{crit,i})^{2/3} i_i^{1/2} \frac{1}{h_{i,t}}\right)h_{i,t}
+            - \sum_k^{inflows} \left( 1/n_k R_k(h_{k,t}-h_{crit,k})^{2/3} i_k^{1/2}  \frac{1}{h_{k,t}}\right)h_{k,t}
+          =  \\
+          \frac{h_{i,t-1} }{\Delta T}
+          + es_{i,t} 
+          + \sum_j^{inflows} a_j h^{b_{j}}_{crit,j}  
+          - inf_{i,t} 
+          - a_i h^{b_{i}}_{crit,i} 
+        \end{multline}
+    The later system of equations is constructed in similar fashion as the case of (\ref{eq:sys}) and solved with the same solver.
 
     \subsection{Impute data requirements and description}