Forest or Farmland?
Clearing forest for crops is tempting. Farmland feeds more people right away, and a growing population needs food. But forests also act as a sponge for rainfall: they recharge groundwater and keep rivers flowing. Remove too much of them and the water that irrigates those fields starts to disappear decades later.
In this article we turn that dilemma into a concrete optimization problem and solve it with Python. We search for the land-use policy that maximizes a combination of population and environmental resource value 50 years from now.
1. The Problem Setting
Consider a region with a total area of $L = 1000\ \mathrm{km^2}$. At the start:
- Forest: $F_0 = 700\ \mathrm{km^2}$
- Farmland: $A_0 = 300\ \mathrm{km^2}$
- Population: $P_0 = 30$ thousand people
- Water stock (aquifer and reservoirs): $W_0 = 800$ million $\mathrm{m^3}$, with capacity $W_{\max} = 1000$ million $\mathrm{m^3}$
Every year we may convert some forest into farmland. We restrict the decision to a simple two-parameter policy:
- $u$: the fraction of the remaining forest converted each year
- $T_c$: the number of years during which conversion continues, after which it stops
The conversion in year $t$ is
$$
c_t =
\begin{cases}
u,F_t & (t < T_c)\
0 & (t \ge T_c)
\end{cases}
$$
and the land-use dynamics are
$$
F_{t+1} = F_t - c_t, \qquad A_{t+1} = A_t + c_t .
$$
2. The Mathematical Model
Food production. Farmland yields $y_0$ thousand people’s worth of food per $\mathrm{km^2}$ at full water supply, scaled by a water factor:
$$
Y_t = y_0, A_t \sqrt{\frac{W_t}{W_{\max}}}
$$
The food ratio $r_t$ is production divided by population:
$$
r_t = \frac{Y_t}{P_t}
$$
Population dynamics. Population grows when $r_t > 1$ and shrinks when $r_t < 1$:
$$
P_{t+1} = P_t\left[1 + g_{\max}\tanh\bigl(2,(r_t - 1)\bigr)\right]
$$
Water dynamics. Forest increases recharge, while farmland irrigation and domestic use withdraw water. A fraction $\lambda$ is lost naturally each year:
$$
W_{t+1} = \mathrm{clip}\Bigl(W_t + a_0 + a_1 F_t - w_a A_t - w_p P_t - \lambda W_t,\ 0,\ W_{\max}\Bigr)
$$
This is the key mechanism. Converting forest cuts recharge $a_1 F_t$ and raises irrigation demand $w_a A_t$ at the same time, so the water stock declines slowly and the damage shows up late.
Objective function. We want to maximize the value at the horizon $T = 50$:
$$
J = w,\frac{P_T}{P_0} + (1-w)\left[\kappa,\frac{F_T}{F_0} + (1-\kappa),\frac{W_T}{W_{\max}}\right]
$$
The first term is the normalized population. The bracket is the normalized environmental value, a mix of forest area and water stock. We use $w = 0.5$ and $\kappa = 0.5$. The optimization problem is
$$
\max_{u,,T_c}\ J(u, T_c)
$$
Parameter values used in the code:
| Symbol | Value | Meaning |
|---|---|---|
| $y_0$ | 0.11 | People fed per $\mathrm{km^2}$ of farmland (thousand) |
| $a_0$ | 25 | Base recharge (million $\mathrm{m^3}$/yr) |
| $a_1$ | 0.30 | Extra recharge per $\mathrm{km^2}$ of forest |
| $\lambda$ | 0.04 | Natural water loss rate |
| $w_a$ | 0.25 | Irrigation withdrawal per $\mathrm{km^2}$ of farmland |
| $w_p$ | 0.60 | Domestic withdrawal per thousand people |
| $g_{\max}$ | 0.03 | Maximum population growth rate |
3. Solution Strategy
The decision space has only two dimensions, so we can evaluate it exhaustively with a grid search:
- $u$: 81 values from 0% to 10% per year
- $T_c$: 51 values from 0 to 50 years
That gives $81 \times 51 = 4131$ policies. Running a 50-year simulation for each one in a Python loop would be slow. Instead, we write the simulation so that every state variable is a NumPy array of shape $(81, 51)$. A single 50-step loop then advances all 4131 policies at the same time.
We also run a sensitivity analysis on the population weight $w$. Because the final states of all policies are stored, we can re-evaluate $J$ for any $w$ without simulating again.
4. Source Code
1 | import numpy as np |
5. Code Walkthrough
Section 1: Parameters. All physical constants and weights are defined at the top, so every experiment in this article can be reproduced by changing a single line. The area is in $\mathrm{km^2}$, population in thousands, and water in millions of $\mathrm{m^3}$.
Section 2: The simulate function. This is the core of the program.
np.broadcast_arrays(u, tc)makes the policy parameters share one shape. We can pass either scalars (a single policy) or a 2D grid (thousands of policies) and use the same code.- Each state variable (
F,A,P,W) is an array of that shape, so onefor t in range(T)loop simulates every policy at once. This is the speed-up: it replaces roughly four thousand 50-step Python loops with one. np.where(t < tc, u * F, 0.0)implements the rule “convert a fraction $u$ of the forest while $t < T_c$, then stop”, elementwise for every policy.- The food ratio uses
np.sqrt(np.clip(W, 0.0, WMAX) / WMAX). The clip keeps the square root safe, andnp.maximum(P, 1e-9)avoids division by zero. - The water update is the equation from Section 2. Forest adds recharge, farmland and people subtract withdrawals, and the whole stock is clipped to $[0, W_{\max}]$.
- Recharge and withdrawal are computed from the land state before this year’s conversion is applied, so the effect of deforestation reaches the water stock with a delay.
- With
keep_path=True, the function also returns the full time series. We use this later to draw trajectories.
Section 3: Grid search. np.meshgrid(..., indexing="ij") builds the policy grid, and a single call to simulate(U, TC) returns the final state of every policy. np.argmax followed by np.unravel_index finds the best cell. The program then prints the optimal policy and compares it against three reference policies: no conversion, a moderate one (3% per year for 20 years), and an aggressive one (8% per year for 50 years).
Section 4: Sensitivity analysis. The objective is cheap to re-evaluate because the final states Ff, Wf, Pf are already stored. For 41 values of the population weight $w$ from 0 to 1, we find the best policy and record the remaining forest share and the population it achieves.
Section 5: Visualization. Nine panels are combined into one figure with GridSpec:
- (a) A 3D surface of the objective $J$ over the policy space, with the optimum marked.
- (b) The same function as a contour map.
- (c) A 3D trajectory of each policy in (forest, water, population) space.
- (d)–(f) Time series of population, forest area, and water stock.
- (g) A 3D bar chart of final population over a coarse policy grid.
- (h)–(i) The trade-off between population and forest as the weight $w$ changes.
6. Execution Results
=== Optimal land-use policy === Annual conversion rate u* : 6.50 % of remaining forest per year Conversion period tc* : 4 years Forest after 50 years : 535.0 km^2 (53.5 % of land) Farmland after 50 years : 465.0 km^2 Population after 50 years : 49.8 thousand Water stock after 50 years : 997.5 million m^3 Objective J* : 1.2700 === Policy comparison (year 50) === Policy Forest Water Pop J No conversion 700.0 1000.0 32.8 1.0475 Optimal 535.0 997.5 49.8 1.2700 Moderate (3%/yr, 20 yr) 380.7 0.0 23.3 0.5238 Aggressive (8%/yr, 50 yr) 10.8 0.0 14.5 0.2460
7. Interpreting the Results
The optimal policy is an early, short burst of conversion. The grid search selects a conversion rate of about 6.5% per year for only 4 years. That converts roughly 165 $\mathrm{km^2}$ of forest, so about 53% of the land stays forested at year 50. Farmland grows from 300 to about 465 $\mathrm{km^2}$, and the population reaches about 49.8 thousand, against about 32.8 thousand if nothing is converted. Meanwhile the water stock stays essentially full, at about 997 million $\mathrm{m^3}$.
Panel (a) and (b): a narrow ridge, not a single peak. The objective surface shows a bright ridge running along the curve where the total amount of forest removed, roughly $u \times T_c$, is moderate. Many combinations of $u$ and $T_c$ that clear about the same area score about the same. The decision that matters is how much forest to remove in total, and the exact schedule matters much less. The surface then falls steeply toward the far corner, where large amounts of forest are cleared.
Panel (d)–(f): the delayed collapse. The aggressive and moderate policies show the trap clearly. Their populations climb faster than the optimal one for the first 15 to 25 years, because more farmland means more food. Their forest keeps shrinking, so recharge falls and irrigation demand keeps rising. Eventually the water stock crashes to zero, in panel (f) after roughly 15 years for the aggressive policy and 30 years for the moderate one. The water factor in food production then collapses, and the population falls, ending below even the no-conversion case. In the end, the short-term gain is far outweighed by the long-term loss.
The comparison at year 50 makes this concrete:
| Policy | Forest ($\mathrm{km^2}$) | Water (Mm$^3$) | Population (thousand) | $J$ |
|---|---|---|---|---|
| No conversion | 700.0 | 1000.0 | 32.8 | 1.0475 |
| Optimal | 535.0 | 997.5 | 49.8 | 1.2700 |
| Moderate (3%/yr, 20 yr) | 380.7 | 0.0 | 23.3 | 0.5238 |
| Aggressive (8%/yr, 50 yr) | 10.8 | 0.0 | 14.5 | 0.2460 |
Panel (c): state-space trajectories. In the 3D view, the optimal and no-conversion paths stay near the high-water plane. The moderate and aggressive paths first move toward lower forest while population rises, then turn and plunge to zero water. The picture shows that the destructive policies pass through a high-population state that cannot be sustained.
Panel (g): a cliff, not a slope. The final-population bars show that the population is not a smooth function of the policy. Cautious policies give moderate populations, a narrow band of mild conversion gives the tallest bars, and heavy conversion causes a sudden drop. This is typical of systems with a tipping point, where the water stock hits its lower bound.
Panel (h) and (i): how values change the answer. As the population weight $w$ approaches 0, the optimum keeps about 70% of the land as forest, which is simply the do-nothing policy. At around $w \approx 0.15$ to $0.2$ the optimum switches to a limited conversion that leaves a little over half the land forested, and it stays there for all larger $w$. Even a society that values only population should not clear everything, because a policy that destroys the water supply also destroys the population. This is why the environmental and population objectives are less opposed than they first appear, at least up to the tipping point.
8. Summary
- We modeled forest-to-farmland conversion with land, water, food, and population dynamics, and searched for the policy that maximizes $J$ over 50 years.
- Vectorizing the simulation lets a single 50-step loop evaluate 4131 policies at once.
- The best policy is a brief, early conversion of about 165 $\mathrm{km^2}$ of forest, which raises the final population by about 50% while keeping the water stock intact.
- Policies that clear more forest look better for 15 to 25 years and then collapse as the water supply fails.
- The timing of damage matters as much as its size. Because deforestation affects water with a delay, the true cost of a decision only becomes visible long after it is made.
The same framework extends naturally to richer models, for example adding soil erosion, a climate-driven rainfall trend, or a continuous-time control law solved with dynamic programming.
