## 3 Solutions
### Solution to the Complete Model (1)
#### 3.2.1 Front Evolution in Time  
Let us apply the direct double integration technique to (1) (with a transformed diffusion term):  
\[\int\limits_{0}^{\delta}\int\limits_{x}^{\delta}\frac{\partial u}{\partial t} dxdx=\frac{D_{0}}{m+1}\int\limits_{0}^{\delta}\int\limits_{x}^{\delta}\frac{ \partial^{2}u^{m+1}}{\partial x^{2}}dxdx+b\int\limits_{0}^{\delta}\int\limits _{x}^{\delta}u^{p}dxdx-c\int\limits_{0}^{\delta}\int\limits_{x}^{\delta}u^{q }dxdx \tag{15}\]  
Figure 1: Barenblatt profile (function) vs. the dimensionless spatial coordinate \(X=x/\delta_{u}\) with a shape controlled by the exponent n. (**a**) Concave profiles for \(n\geq 1\); (**b**) Convex profiles for \(n\leq 1\). Note: For ease of comparison with plots in panel (**a**), the exponents in panel (**b**) are specifically displayed as \(1/n\). Now, applying the Leibniz rule to the integral on the left-hand side, we have  
\[\frac{d}{dt}\int\limits_{0}^{\delta}\int\limits_{x}^{\delta}u(x,t)dxdx=\frac{D_{0 }}{m+1}u^{m+1}(0,t)+b\int\limits_{0}^{\delta}\int\limits_{x}^{\delta}u^{p}dxdxx- \int\limits_{0}^{\delta}\int\limits_{x}^{\delta}u^{q}dxdx \tag{16}\]  
At this step, we did not touch the absorption and source terms since their integration, strongly depends on the choice of the assumed profile. Then, using Equation (15) and replacing \(u(x,t)\) by \(u(\frac{x}{\delta},t)\), the integration yields  
\[\frac{1}{(n+1)(n+2)}\frac{d\delta^{2}}{dt}=\frac{D_{0}}{m+1}u_{s}^{m+1}\int \limits_{0}^{\delta}\int\limits_{x}^{\delta}\frac{\partial^{2}u^{m+1}}{ \partial x^{2}}+b\frac{\delta^{2}}{(np+1)(np+2)}-c\frac{\delta^{2}}{(nq+1)(nq+ 2)} \tag{17}\]  
A rearrangement of Equation (17) allows us to see the equation about \(\delta(t)\) as  
\[\frac{d\delta^{2}}{dt}=D_{0}\frac{(n+1)(n+2)}{m+1}+\delta^{2}b\bigg{(}\frac{1 }{(np+1)(np+2)}-\frac{c}{b}\frac{1}{(nq+1)(nq+2)}\bigg{)}(n+1)(n+2) \tag{18}\]  
In a compact form, we have  
\[\frac{d\delta^{2}}{dt}=A+B\delta^{2} \tag{19}\]  
where  
\[ \begin{split}& A=D_{0}\frac{(n+1)(n+2)}{m+1}=D_{0}N(n,m),\\ & B=b\bigg{(}\frac{1}{(np+1)(np+2)}-\frac{c}{b}\frac{1}{(nq+1)( nq+2)}\bigg{)}(n+1)(n+2)=bM(n,p,q)\end{split}  \tag{20}\]  
Hereafter, for the sake of simplicity of the expression, we will use only the notations \(N\) and \(M\). We have a separable equation with a solution (taking into account the physically based condition \(\delta(t=0)=0\):  
\[\delta=\sqrt{\frac{A}{B}}\sqrt{e^{Bt}-1} \tag{21}\]  
Now, taking into account that \(A\equiv D_{0}\) and \(B\equiv b\), we have that \(\delta\equiv\sqrt{D_{0}/b}\) that makes the result Equation (21) physically correct. Moreover, the scaled depth is \(\delta/\sqrt{D_{0}/b}=f(Fo)\) because of the product \(Bt\equiv bt=Fo\). In a more detailed form, Equation (21) can be presented as  
\[\delta(t)=\sqrt{\frac{D_{0}}{b}}\sqrt{\frac{N}{M}}\sqrt{e^{MFo}-1} \tag{22}\]  
For a short time, i.e., for \(t_{0+}\), we may approximate \(e^{MFo}\approx 1+MFo\), and in such a case \(\delta(t)\approx\sqrt{D_{0}t}\sqrt{N}\) corresponds to a Gaussian diffusion. For a long time when \(e^{MFo}\gg 1\), we have \(\delta(t)\approx\sqrt{\frac{D_{0}}{b}}\sqrt{\frac{N}{M}}\sqrt{e^{MFo}}\)  
**Remark 3** (How the degeneracy should be accounted for).