P Reduction — Equations and Routing

---

Overview

The P model is substantially more complex than the N model because: 1. P loss from fields occurs via three distinct physical pathways (erosion, macropore, matrix). 2. P is routed through a chain of lakes, not directly to the coast — with a transport matrix capturing retention in each lake. 3. Stream engineering measures (ochre, sand, re-meandering, raising, trees) reduce P at the watercourse level, not the field level. 4. P wetlands (PWET) are constrained to not exceed total catchment P loss and interact with stream measures.

---

P loss pathways at field level

For each field i and measure j, total P effectiveness is the sum of three components:

P_Effects_Total(i,j) = P_Effects_erosion(i,j) + P_Effects_macropore(i,j) + P_Effects_matrix(i,j)

See P loss pathways for the physical meaning of each pathway.

Matrix pathway

Slow subsurface transport through the soil matrix. Effect proportional to matrix_ha(i) (kg P/ha/yr background loss via matrix flow).

P_Effects_matrix(i,"NPB10") = PotV(i,"NPB10") × matrix_ha(i) × 0.05
P_Effects_matrix(i,"NPB20") = PotV(i,"NPB20") × matrix_ha(i) × 0.20
P_Effects_matrix(i,"PPC_NPB10") = 0.9 × PotV(i,"NPB10") × matrix_ha(i) × 0.05

(For constructed wetlands, MW uses its own matrix effect formula averaged over the retention area.)

Macropore pathway

Fast preferential flow through soil macropores — more important on heavy clay soils. Effect based on macropore_field(i) (kg P/yr total from macropore flow).

P_Effects_macropore(i,"NPB10") = macropore_field(i) × 0.05     [if NPB10 potential > 0]
P_Effects_macropore(i,"FO")    = macropore_field(i) × 0.35 × (PotV(i,"FO") / IniPotV(i))
P_Effects_macropore(i,"LRh")   = (1-lav(i)) × macropore_field(i) × 0.35

Erosion pathway

Surface erosion of P-rich topsoil into streams. Effect based on erosion_field(i) (kg P/yr total erosion loss).

P_Effects_erosion(i,"PPC")    = erosion_field(i) × 0.9
P_Effects_erosion(i,"OT")     = erosion_field(i) × 0.5
P_Effects_erosion(i,"BZ10")   = erosion_field(i) × 0.62   [if BZ10Pot(i) > 0.2]
P_Effects_erosion(i,"BZ20")   = erosion_field(i) × 0.75   [if BZ20Pot(i) > 0.2]
P_Effects_erosion(i,"LRl")    = 0   [set to zero — low-lying fields]
P_Effects_erosion(i,"Pwet")   = 15 × PotV(i,"Pwet")
P_Effects_erosion(i,"FO")     = erosion_field(i) × (PotV(i,"FO") / IniPotV(i))
P_Effects_erosion(i,"LRh")    = (1-lav(i)) × erosion_field(i)

Note on BZ threshold: A threshold of BZ10Pot(i) > 0.2 (20% of field within BZ potential) is required for P erosion effects. Reasoning: BZ potential is a fraction, so only a portion of the field needs to qualify for the measure to protect the whole field's erosion losses. ⚠️ The appropriate threshold is noted as a design decision still open in the code.

---

P reduction in upstream lake catchments (P_Red_bf_ret_eq_up)

P_Red_bf_ret_up(up_lakecatch) =e=
    Sum(ret$up_LakeRet(up_lakecatch,ret), MW_red_p(ret))                                          [MW constructed wetlands]
  + Sum(i$up_Lake_i(up_lakecatch,i), Sum(j$pho(j),  P_Effects_Total(i,j) * x(i,j)))              [P field measures]
  + Sum(i$up_Lake_i(up_lakecatch,i), Sum(j$nit(j),  P_Effects_Total(i,j) * x(i,j)))              [N measures with P co-benefit]
  + Sum(i$up_lake_i(up_lakecatch,i), Sum(j$comb(j), P_Effects_Total(i,j) * x(i,j)))              [combined measures]
  + Sum(i$up_lake_i(up_lakecatch,i), Sum(j$comb2(j),P_Effects_Total(i,j) * x(i,j)))              [combined with BZ]
  + WWT_P_Red(up_lakecatch)                                                                         [WWT upgrades]
  + Ochre_P_eff(up_lakecatch) + Sand_P_eff(up_lakecatch) + IBZ_P_eff(up_lakecatch)                [stream point measures]
  + Re_meandering_P_eff(up_lakecatch) + Raising_P_eff(up_lakecatch)                               [stream restoration]
  + Trees_P_eff(up_lakecatch) + OF_p_red(up_lakecatch)                                             [trees + overflow]

Important: N measures (nit set) contribute P co-benefits — BZ10, BZ20, FO, LRh, etc. all reduce P in addition to N. These P effects are also captured through P_Effects_Total.

---

Lake-chain routing (P_Red_aft_ret_eq)

P_Red_aft_ret(lakecatch) =e=
    Sum(up_lakecatch$lakechain(lakecatch,up_lakecatch),
        P_Red_bf_ret_up(up_lakecatch) × p_transp_matrix(lakecatch, up_lakecatch))

The p_transp_matrix(lakecatch, up_lakecatch) is a transport coefficient — the fraction of P reduced upstream that arrives at the downstream lake catchment after passing through any intermediate lakes. It is loaded from p_transp_matrix_NEW.csv.

A parallel matrix p_transp_matrix_k(k, up_lakecatch) routes P all the way to coastal catchments for coastal-level P aggregation.

---

P target constraint (P_target_eq)

P_target_eq(lakecatch)$(P_targets(lakecatch)>0)..
    P_Red_aft_ret(lakecatch) + exceed_p(lakecatch) =g= var_P × P_targets(lakecatch)

• Only active for lake catchments with a defined target (P_targets > 0).
• exceed_p(lakecatch) carries Penalty2 = 9.999×10¹³ DKK — an order of magnitude larger than the N penalty, reflecting the relative difficulty of meeting P targets.
• var_P = 0 deactivates all P targets.
• P targets loaded from P_targets_Dec2024.inc (443 records; 189 active; max 9,155.79 kg P/yr at lakecatch=456; total ~117,089 kg P/yr).

---

Stream measure P effects

Ochre traps

Ochre_P_eff(up_lakecatch) = Sum(w$up_lakew, w_okker(w) × P_Effects_Ochre(w) × v_w(w,"ochre"))
P_Effects_Ochre(w) = 140 kg P/yr   [if w_class = 1 or 2]

Sand traps

Sand_P_eff(up_lakecatch) = Sum(w$up_lakew, w_sand(w) × P_Effects_Sand(w) × v_w(w,"sand"))
P_Effects_Sand(w): 7 kg P/yr (class 1 or 2, geo zone 1)
                   26 kg P/yr (class 1 or 2, geo zone 2)
                   12 kg P/yr (class 1 or 2, geo zone 3)

Only watercourse classes 1 and 2 have sand trap potential and effects.

Re-meandering

Re_meandering_P_eff(up_lakecatch) = Sum(w$up_lakew, w_meandering(w) × w_P_red_re_meandering(w) × v_w(w,"re_meandering"))

Effect is watercourse-specific from w_P_red_re_meandering.inc.

Raising (bed raising)

Raising_P_eff(up_lakecatch) = Sum(w$up_lakew, w_raising(w) × w_P_red_raising(w) × v_w(w,"raising"))

Effect is watercourse-specific from w_P_red_raising.inc.

Trees on erosion stretches

Trees_P_eff(up_lakecatch) = Sum(ero_stretch$up_lake_ero, tree_eff(ero_stretch) × v_str(ero_stretch,"trees"))

Effect is stretch-specific from Tree_eff.inc.

IBZ (in-stream buffer zone)

IBZ_P_eff(up_lakecatch) = Sum(i$(up_lake_i AND PotV(i,"IBZ")>0), IBZ_eff(i) × x(i,"IBZ"))

Effect is field-specific from IBZ_eff.inc.

---

P wetland constraints

P wetland (PWET) is linked to specific stream reaches through a three-equation system:

Pwet1_eq(w)..          s(w) ≥ Sum(w_Pwet(w_m), v_w(w,w_m))
Pwet2_eq(w,i)$wi_pwet.. x(i,"Pwet") + s(w) ≤ 1
Pwet3_eq(w)..          Sum(w_Pwet(w_m), v_w(w,w_m)) ≥ s(w)

This ensures that a P wetland on field i and a stream measure (re-meandering or raising) on the adjacent watercourse w cannot both be selected. s(w) is a binary indicator variable.

Additionally:

P_Wet_limit_eq(up_lakecatch)..
    Sum(i$up_lake_i, P_Effects_Total(i,"Pwet") × x(i,"Pwet")) ≤ Total_P_loss(up_lakecatch)

Total P reduction from P wetlands cannot exceed the total P load of the catchment. This prevents the model from claiming more P reduction than physically possible.

---

Mini-wetland P effects

MW P effects are aggregated at retention area (ret) level:

MW_red_P(ret) = sum over MW sizes [area × (P_Effects_MW_Matrix(ret) + P_Effects_MW_macropore(ret)) × MW_binary]

Where: - P_Effects_MW_Matrix(ret) = 0.45 × weighted average matrix_ha for fields in ret with macropore loss - P_Effects_MW_macropore(ret) = 0.45 × weighted average macropore_field for fields in ret with macropore loss

(Note: Only fields with macropore loss > 0 are included in the average — per Hans Estrup, Nov 2024)

---

Open questions

1. WL P effect was set to 0 in Nov 2023. Is this permanent for VP3 or might it be revised?
2. Why is LRl P (erosion) = 0 while LRh has erosion effects? Is this because low-lying land is already wet?
3. The BZ20% threshold for P erosion effects — has this been validated or is it a pragmatic choice?
4. ✅ Total_P_loss(up_lakecatch) — confirmed: from Total_P_loss.inc, ~448 entries, baseline P load per upstream lake catchment (kg P/yr); source Total_P_loss.xlsx (confirmed 2026-04-05).

---

Related pages

Concepts: - P loss pathways — physical mechanisms for each pathway - Spatial Hierarchy — ID15 → upstream lake catchment → lake catchment → coast

Field measures with P effects: - PPC — 90% erosion reduction - OT — 50% erosion reduction - PWET — P wetland; sedimentation + catchment P cap constraint - NPB10, NPB20 — all three P pathways - BZ10, BZ20 — erosion P co-benefit - FO, LRh — erosion + macropore P co-benefits

Infrastructure measures: - IBZ — in-stream buffer zone - Stream measures — ochre, sand, re-meandering, raising, trees - MW — mini-wetlands; matrix + macropore P - WWT — wastewater P reduction