# Simulating a partial differential equation, uniform vegetation and bare soil system (2 species)

I want to simulate a system of two plants ,woody and herbaceous species; with starting condition that half of the grid (square grid) is filled with woody specie, the herbaceous specie exist on the line between the woody specie and the bare soil.

The system model is the following:

$$\frac{\partial B_1}{\partial t}= \Lambda_1 WB_1(1-B_1/K_1)(1+E_1B_1)^2-M_1B_1+D_{B_1}\nabla^2B_1 ,$$

$$\frac{\partial B_2}{\partial t}= \Lambda_2 WB_2(1-B_2/K_2)(1+E_2B_2)^2-M_2B_2+D_{B_2}\nabla^2B_2 ,$$

$$\frac{\partial W}{\partial t} = I H - \frac{N W}{1+RB_1/K_1} -\Gamma_1WB_1(1+E_1B_1)^2-\Gamma_2WB_2(1+E_2B_2)^2+D_W\nabla^2W,$$

$$\frac{\partial H}{\partial t}=P-IH +D_H\nabla^2(H^2)$$

where

$$I=A\frac{Y_1B_1+Y_2B_2+qf}{Y_1B_1+Y_2B_2+q}.$$

subject to Newman Boundary Conditions and Initial Conditions i mentioned above, I wrote the following code, but with no success.

L = 1;(*length of square*)
T = 8;(*Time integration*)
pts = 100;
f = 0.2; (*The infiltration contrast between vegetated and bare-soil \
areas*)
A = 40.0; (*The maximal infiltration rate, obtained in densely \
vegetated areas*)
n = 30.0; (*the soil-water evaporation rate*)
R = 30.0; (*the shading effect*)
E1 = 2.0; (*root to shoot ratio of B1*)
E2 = 0.0;(*root to shoot ratio of B2*)
M1 = 0.75; (*the mortality rate of B1*)
M2 = 7.5; (*the mortality rate of B2*)
K1 = 5.0; (*the maximum standing biomass of B1*)
K2 = 0.5; (*the maximum standing biomass of B2*)
\[CapitalLambda]1 = 0.05; (*growth rate of B1*)
\[CapitalLambda]2 = 1.0; (*growth rate of B2*)
\[CapitalGamma]1 = 10.0; (*the water uptake rate of B1*)
\[CapitalGamma]2 = 15.0; (*the water uptake rate of B2*)
q = 5.0; (*reference biomass value beyond which he effect of \
vegetation on infiltration becomes significant*)
Y1 = 1.0; (*relative contribution of the woody specie (B1)*)
Y2 = 100; (*relative contribution of the herbaceous specie (B2)*)
DB1 = 0.1; (*the biomass expansion rate of B1*)
DB2 = 0.05; (*the biomass expansion rate of B2*)
DW = 1.0; (*the soil-water diffusion*)
DH = 1.0; (*the rate of surface-water spread*)
i = A (Y1 B1[t, x, y] + Y2 B2[t, x, y] + q f)/(
Y1 B1[t, x, y] + Y2 B2[t, x, y] + q );
P = 2;
(*system of nonlinear PDE*)
pde = {D[B1[t, x, y],
t] == \[CapitalLambda]1 W[t, x, y] B1[t, x,
y] (1 - B1[t, x, y]/K1) (1 + E1 B1[t, x, y])^2 -
M1 B1[t, x, y] +
DB1 (D[B1[t, x, y], x, x] + D[B1[t, x, y], y, y]),
D[B2[t, x, y],
t] == \[CapitalLambda]2 W[t, x, y] B2[t, x,
y] (1 - B2[t, x, y]/K2) (1 + E2 B2[t, x, y])^2 -
M2 B2[t, x, y] +
DB2 (D[B2[t, x, y], x, x] + D[B2[t, x, y], y, y]),
D[W[t, x, y], t] ==
i H - n W[t, x, y]/(
1 + R B1[t, x, y]/K1 ) - \[CapitalGamma]1 W[t, x, y] B1[t, x,
y] (1 + E1 B1[t, x, y])^2 - \[CapitalGamma]2 W[t, x, y] B2[t,
x, y] (1 + E2 B2[t, x, y])^2 +
DW (D[W[t, x, y], x, x] + D[W[t, x, y], y, y]),
D[H[t, x, y], t] ==
P - i H[t, x, y] +
DH (D[H[t, x, y], x, x]^2 + D[H[t, x, y], y, y]^2)};
(*Newman boundary condition*)
bc = {(D[B1[t, x, y], x] /. x -> -L) ==
0, (D[B1[t, x, y], x] /. x -> L) ==
0, (D[B1[t, x, y], y] /. y -> -L) ==
0, (D[B1[t, x, y], y] /. y -> L) ==
0, (D[B2[t, x, y], x] /. x -> -L) ==
0, (D[B2[t, x, y], x] /. x -> L) ==
0, (D[B2[t, x, y], y] /. y -> -L) ==
0, (D[B2[t, x, y], y] /. y -> L) ==
0, (D[W[t, x, y], x] /. x -> -L) ==
0, (D[W[t, x, y], x] /. x -> L) ==
0, (D[W[t, x, y], y] /. y -> -L) ==
0, (D[W[t, x, y], y] /. y -> L) == 0,
(D[H[t, x, y], x] /. x -> -L) == 0, (D[H[t, x, y], x] /. x -> L) ==
0, (D[H[t, x, y], y] /. y -> -L) ==
0, (D[H[t, x, y], y] /. y -> L) == 0};

ic = {B1[0, x, y] ==
Interpolation[
Flatten[Table[{x, y, 1}, {x, -L, L, 2/pts}, {y, -L, 0, 2/pts}],
1]][x, y],
B1[0, x, y] ==
Interpolation[
Flatten[Table[{x, y, 0}, {x, -L, L, 2/pts}, {y, 0, L, 2/pts}],
1]][x, y],
B2[0, x, y] ==
Interpolation[
Flatten[Table[{x, y, 1}, {x, -L, L, 2/pts}, {y, L, L, 2/pts}],
1]][x, y],
W[0, x, y] ==
Interpolation[
Flatten[Table[{x, y, RandomReal[]}, {x, -L, L, 2/pts}, {y, -L,
L, 2/pts}], 1]][x, y],
H[0, x, y] ==
Interpolation[
Flatten[Table[{x, y, RandomReal[]}, {x, -L, L, 2/pts}, {y, -L,
L, 2/pts}], 1]][x, y]};

eqns = Flatten@{pde, bc, ic};

sol = NDSolve[
eqns, {B1[t, x, y], B2[t, x, y], W[t, x, y], H[t, x, y]}, {t, 0,
T}, {x, -L, L}, {y, -L, L},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> pts, "MaxPoints" -> pts},
Method -> {"IDA", "ImplicitSolver" -> {"GMRES"}}}];

GraphicsGrid[Table[time = t0 + 3*t1;
DensityPlot[B1[time, x, y] /. sol, {x, -L, L}, {y, -L, L},
ColorFunction -> "SunsetColors",
PlotLabel -> "t=" <> ToString[time], Ticks -> False], {t1, 0,
2}, {t0, 0, 2}], ImageSize -> 600]


Can someone help me? Thanks in advance.

• A few things. First, your equation for dW/dt has H instead of H[t, x, y]. Second, the initial conditions for B1 are repeated. The error messages you got could help find those problems, don't ignore them! Less crucial, it would be easier to define the initial conditions as B1[0, x, y] == 1, B2[0, x, y] == 1. Also, your implementation of the nonlinear diffusion of H does not seem to match your equations. After fixing all those things, it runs but takes a long time because the equations may be stiff due to parameter values. BTW, do you have a reference for this model? Feb 18, 2020 at 17:43
• first of all, thank you! secondly, but if choose B1[0,x,y]==1, B2[0,x,y]==1 if wouldnt take into account that just half of space have them (B1), and that B2 is only located on the line between the dryland and the vegetation. This model is used in this paper. Feb 18, 2020 at 18:28
• Whoops, sorry I misinterpreted your intention for the initial conditions. Try B1[0, x, y] == If[y < 0, 1, 0], B2[0, x, y] == If[Abs[y] < 0.1, 1, 0] (I think you need a finite thickness of the B2 band). Feb 18, 2020 at 18:51