11. Soil P Dynamics

Apart from input and uptake, the P concentration in the topsoil and subsoil is affected by P accumulation. Modelling of P accumulation or release is included by using a Langmuir equilibrium, supplemented with rate limited diffusion, based on the approach used in the INITIATOR model (van der Salm et al., 2016#; #De Vries et al., 2023).

Division of pools

The soil P pool is divided into:

  1. An inert P pool with no change over time,

  2. A stable P pool, S, changing slowly according to a rate limited reaction with dissolved P in the soil solution, and

  3. A labile P pool, L, changing rapidly according to an equilibrium reaction with dissolved P in soil solution.

Initialization

The sum of the labile and stable pool is assumed to be approximated by oxalate extractable P (Pox), further denoted as the reactive P pool. The initial size of Pox must be provided by user. The labile pool, L, is at the start assumed to be 1/3 of Pox, and the stable pool, S, thus being equal to 2/3 of Pox.

The maximum size of Pox is assumed to be 1/2 of (Al + Fe)ox, with (Al + Fe)ox being the oxalate extractable Al and Fe contents (mmol kg–1).

Equation 11.1

Soil P concentrations

The concentration of inorganic P in soil solution, [Pi], is calculated from the size of the labile P pool (L) according to a Langmuir equation:

Equation 11.2

where:

is the size of the adsorbed labile P pool (mmol kg–1).

is the theoretical maximum adsorption (mmol kg–1), which is 1/6 of (Al+Fe)ox.

is a Langmuir adsorption constant (m3 kg–1).

Soil Texture

Peat

Sand

Loam

Clay

Soil Texture Class

1

2

3

500

1000

1500

2000

is the inorganic P concentration (g m–3, or mg L–1). The value of 1000 is used to convert result from kg m–3 to g m–3 (mg L–1).

  • The size of L is bounded by Lmin and Lmax.

  • [Pi] is bounded to [0.001, 90].

If any calculated value falls outside the boundaries, the value must be clipped. See Minimum and maximum pool sizes for further explanation.

The total (both inorganic and organic) P concentration, Pt, is then derived using an exponential relationship between total P and inorganic P based on measurements in soil solution, drainage water and surface water, as reported by Chardon et al. (2007) (in Dutch):

Equation 11.3

where:

is the total P concentration (g m–3, or mg L–1).

is the inorganic P concentration (g m–3, or mg L–1).

The transfer between P in the soil solution and the stable pool, Ptrans, is described by a rate-limited Freundlich equation:

Equation 11.4

where:

is the rate constant (day–1) for the transfer from soil solution to the stable pool, with a default value of 0.0014.

is the rate constant (day–1) for the transfer from the stable pool to soil solution, with a default value of 2 × 10–6 for sand, and 44 × 10–6 for clay and peat soils.

is the Freundlich constant of the stable pool (mmol kg–1 (mg L–1)–n).

is the inorganic P concentration (g m–3, or mg L–1).

is the exponent of stable pool, with a default value of 0.26.

is the size of the stable pool (mmol kg–1).

The value of Ptrans can be positive or negative.

When Ptrans > 0, dissolved P transfers to S; and when Ptrans < 0, S transfers to dissolved P.

P surplus

P surplus is simply calculated as the difference between P input and crop P uptake:

Equation 11.5

where:

is the annual P deposition (kg P ha–1 yr–1).

is the annual P input from mineral fertilisers (kg P ha–1 yr–1).

is the annual P input from animal manure (kg P ha–1 yr–1).

is the annual P removed from soil by crop uptake (kg P ha–1 yr–1).

Runoff and leaching

P runoff and leaching are calculated based on the total P concentration and the waterflux:

Equation 11.6

where:

is the total P concentration (g m–3, or mg L–1).

is the daily interflow (subsurface runoff) (m3 ha–1 day–1).

is the daily leaching effluent to shallow groundwater (m3 ha–1 day–1).
For calculations of Qint and Qeff, refer to 5. Water Fluxes.

P accumulation

P accumulation (Pacc) is calculated as:

Equation 11.7

where:

the P surplus (g ha–1 day–1).

is the loss of P by leaching and runoff (g ha–1 day–1).

Pool dynamics

The changes in the labile and stable pools are calculated on a daily basis.

The change in the labile P pool is calculated by a mass balance:

Equation 11.8

where:

are the size of the labile P pool (mmol kg–1) at time t + 1 and time t, respectively.

is the P accumulation (g ha–1 day–1).

is the transfer flux from the soil solution to the stable pool of adsorbed P (mmol kg–1 day–1).

is the length of the time step (1 day).

To convert the unit of and from g ha–1 day–1 to mmol kg–1 day–1:

where:

is the thickness of the soil layer (cm).

the bulk density of the soil (g cm–3).

is the molar mass of P (31 g mol–1).

The change in the stable P pool is calculated as:

Equation 11.9

where:

are the size of the stable P pool (mmol kg–1) at time t + 1 and time t, respectively.

is the transfer flux from the soil solution to the stable pool of adsorbed P (mmol kg–1 day–1).

is the length of the time step (1 day).

The change in the reactive P pool, Pox, can be determined as:

Equation 11.10

Minimum and maximum pool sizes

From Equation 11.2 it can be deduced that when L approaches L'max, [Pi] can become infinitely high. To maintain [Pi] in a reasonable range, we set for [Pi] a minimum value of 0.001 mg L–1, and a maximum value of 90 mg L–1.

The minimum and maximum boundaries Lmin and Lmax may be determined as:

Equation 11.11

The upper boundary of S pool is equal to Smax, which is defined as:

Equation 11.12

The calculated size of L and S in each time step is checked so that they do not exceed the boundaries:

Equation 11.13

Buffer pool for excess P

Following Equation 11.13, in the case of Lt > Lmax or St > Smax, Lt and St are clipped to their respective maximums. The excess P (Lt − Lmax or St − Smax) must go to a provisional buffer pool B, which acts as a precipitation pool to temporarily store these excess P.

The buffer pool also have precedence over L and S in P mining. When B > 0, any transfer of adsporped P to soil solution will first come from B, and any remainder will then come from L and S.

The following pseudocode in Python demonstrates how the buffer pool is used. The dynamics of S is calculated before the dynamics of L.

P_trans = ...  # Transfer flow between dissolved P and adsorped stable P pool.
B = ...        # Buffer precipitation pool.
dB = ...       # Change in the B pool.
S = ...        # Adsorped stable P pool.
S_max = ...    # Maximum size of the S pool.

if P_trans < 0:  (1)
    if abs(P_trans) > B:  (2)
        dB = -B
        P_trans = P_trans + B
    else:  (3)
        dB = P_trans
        P_trans = 0

else:  (4)
    if S + P_trans < S_max:  (5)
        dB = 0
    else:  (6)
        dB = S + P_trans - S_max

(7)
B = B + dB
S = S + P_trans
S = min(S, S_max)
1 Transfer from adsorped stable pool to dissolved P. B is used before S.
2 B is smaller than Ptrans, B is depleted, and Ptrans is offset by B.
3 B is larger than Ptrans, Ptrans is set to 0 since no P will be coming from S, and B is reduced by Ptrans.
4 Transfer from dissolved P to adsorped stable pool. B acts as a buffer pool for excess P.
5 No excess P. B remains unchanged.
6 Excess P goes into B. No need to adjust the sizes of Ptrans or S, since S + Ptrans will be clipped to Smax.
7 Update the sizes of B and S. 
P_acc = ...    # P accumulation.
P_trans = ...  # Transfer flow between dissolved P and adsorped stable P pool.
B = ...        # Buffer precipitation pool.
dB = ...       # Change in the B pool.
L = ...        # Adsorped labile P pool.
L_min = ...    # Minimum size of the L pool.
L_max = ...    # Maximum size of the L pool.

dL = P_acc - P_trans  (1)

if dL < 0:  (2)
    if abs(dL) > B:  (3)
        dB = -B
        dL = dL + B
    else:  (4)
        dB = dL
        dL = 0

else:  (5)
    _L = L + dL

    if _L < L_min:  (6)
        if abs(_L - L_min) > B:  (7)
            dB = -B
        else:  (8)
            dB = _L - L_min
    elif _L > L_max:  (9)
        dB = _L - L_max
    else:  (10)
        dB = 0

(11)
B = B + dB
L = L + dL
L = max(L_min, min(L, L_max))
1 Change in size of the L pool (ΔL) is determined according to Equation 11.8.
2 Ptrans > Pacc > 0, B is used before L.
3 B is smaller than ΔL, B is depleted, and ΔL is offset by B.
4 B is larger than ΔL, ΔL is set to 0 since no P will be coming from L, and B is reduced by ΔL.
5 Pacc > Ptrans and/or Ptrans < 0, and the size of L will increase to a provisional new size L'.
6 When L' < Lmin, L' will be padded to Lmin, and the difference ( |L' − Lmin| ) will come from B.
7 B is smaller than the difference, B is depleted. Note that in this case a small amount of P ( |L' − Lmin| − B ) will appear to be coming out of nowhere. However, since this value will be very small, it may be safely negelected.
8 B is larger than the difference, B is reduced by the amount of the difference.
9 When L' > Lmax, B acts as a buffer pool for excess P.
10 When Lmin < L' < Lmax, B remains unchanged.
11 Update the sizes of B and L.