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).

11.1. Division of pools

The soil P pool is divided into:

An inert P pool with no change over time,
A stable P pool, S, changing slowly according to a rate limited reaction with dissolved P in the soil solution, and
A labile P pool, L, changing rapidly according to an equilibrium reaction with dissolved P in soil solution.

11.2. Initialization

The sum of the labile and stable pool is assumed to be approximated by oxalate extractable P (P_ox), further denoted as the reactive P pool. The initial size of P_ox must be provided by user. The maximum size of P_ox 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). The labile pool, L, is at the start assumed to be 1/3 of P_ox, and the stable pool, S, thus being equal to 2/3 of P_ox.

Equation 11.1

where:

	are the initial sizes (mmol kg^–1) of the labile and stable P pools, respectively.
	is the initial size of the oxalate extractable P (mmol kg^–1).

11.3. Soil P concentrations

The concentration of inorganic P in soil solution, [P_i], 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). See Equation 11.11.

is a Langmuir adsorption constant (m³ 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, i.e., mg L^–1). The value of 1000 is used to convert result from kg m^–3 to g m^–3.

The calculated values of L and [P_i] are constrained by the following boundaries:

L_min ≤ L ≤ L_max
0.001 ≤ [P_i] ≤ 90

See Minimum and maximum pool sizes for further explanation.

The total (both inorganic and organic) P concentration, P_t, 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, i.e., mg L^–1).
	is the inorganic P concentration (g m^–3, i.e., mg L^–1). See Equation 11.2.
	is the base for natural logarithms, approximately 2.71828.

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

Equation 11.4

where:

is the P transferred between the soil solution and the stable pool (mmol kg^–1 day^–1).

The value of P_trans can be positive or negative.

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

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).

L_max and S_max are the maximum sizes of the labile and stable P pool, respectively (mmol kg^–1). See Equation 11.11 and Equation 11.12.

is the inorganic P concentration (g m^–3, i.e., mg L^–1). See Equation 11.2.

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

is the size of the stable pool (mmol kg^–1), as calculated in Equation 11.9.

11.4. P surplus

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

Equation 11.5

where:

	is the difference between P input and crop P uptake (g P ha^–1 day^–1).
	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).

11.5. Runoff and leaching

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

Equation 11.6

where:

	is the P loss via surface runoff (g P ha^–1 day^–1).
	is the P leached below root zone (g P ha^–1 day^–1).
	is the combined P loss via runoff and leaching (g P ha^–1 day^–1).
	is the total P concentration (g m^–3, i.e., mg L^–1). See Equation 11.3.
	is the daily interflow (subsurface runoff) (m³ ha^–1 day^–1). See 5. Water Fluxes.
	is the daily leaching effluent to shallow groundwater (m³ ha^–1 day^–1). See 5. Water Fluxes.

11.6. P accumulation

P accumulation (P_acc) is calculated as:

Equation 11.7

where:

	is the P accumulation in soil (g P ha^–1 day^–1).
	the P surplus (g P ha^–1 day^–1). See Equation 11.5.
	is the loss of P by leaching and runoff (g P ha^–1 day^–1). See Equation 11.6.

11.7. 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 sizes of the labile P pools (mmol kg^–1) at time t and time t–1 (the previous time step), respectively.

is the P accumulation (g ha^–1 day^–1). See Equation 11.7.

To convert the unit of 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).

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

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

The change in the stable P pool is calculated as:

Equation 11.9

where:

	are the sizes of the stable P pools (mmol kg^–1) at time t and time t–1 (the previous time step), 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, P_ox, can be determined as:

Equation 11.10

where:

	are the sizes of the reactive P pool (mmol kg^–1) at time t and time t–1 (the previous time step), respectively.
	are the sizes of the labile and stable P pools (mmol kg^–1) at time t, respectively. See Equation 11.8 and Equation 11.9 .
	is the P accumulation (g ha^–1 day^–1). See Equation 11.7.
	is the length of the time step (1 day).

Minimum and maximum pool sizes

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

The minimum and maximum boundaries L_min and L_max may be determined as:

Equation 11.11

where:

	is the theoretical maximum size of the adsorbed labile P pool (mmol kg^–1).
	is the oxalate extractable Al and Fe contents (mmol kg^–1).
	are the applicable minimum and maximum sizes of the labile P pools (mmol kg^–1), respectively.
	is a Langmuir adsorption constant (m³ kg^–1). See Equation 11.2.

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

Equation 11.12

where:

	is the maximum size of the stable P pool (mmol kg^–1).
	is the oxalate extractable Al and Fe contents (mmol kg^–1).
	are the maximum size of the labile P pool (mmol kg^–1). See Equation 11.11.

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

Equation 11.13

where:

	are the sizes of the labile and stable P pools (mmol kg^–1) at time t, respectively. See Equation 11.8 and Equation 11.9.
	are the applicable minimum and maximum sizes of the labile P pools (mmol kg^–1), respectively. See Equation 11.11.
	is the maximum size of the stable P pool (mmol kg^–1). See Equation 11.12.

Buffer pool for excess P

Following Equation 11.13, in the case of L_t > L_max or S_t > S_max, L_t and S_t are clipped to their respective maximums. The excess P (L_t − L_max, or S_t − S_max) 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 adsorbed 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 adsorbed stable P pool.
B = ...        # Buffer precipitation pool.
dB = ...       # Change in the B pool.
S = ...        # Adsorbed 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	A negative P_trans indicates transfer from adsorbed stable pool S to dissolved P. The buffer pool B is mined before the S pool.
2	If the buffer pool size B is smaller than what’s required for P_trans, B is depleted (i.e., ΔB = –B), and P_trans is offset by B (P_trans = P_trans + B, where P_trans < 0 and B > 0).
3	If the size of B is larger than the requirement by P_trans, P_trans is set to 0 as all transfer to dissolved P will come from the buffer and none from S, and B is reduced by the amount of P_trans (ΔB = P_trans, where P_trans < 0).
4	A positive P_trans indicates transfer from dissolved P to adsorbed stable pool S. B acts as a buffer for any excess P.
5	If the size of S after receiving the transfer (S + P_trans) is less than S_max, there is no excess P going into the buffer B (i.e., ΔB = 0).
6	If the size of S after receiving the transfer (S + P_trans) is greater than S_max, S will be constrained to S_max. The excess P (S + P_trans − S_max) will go into B. There is no need to adjust the sizes of P_trans or S, since S + P_trans will automatically become S_max in later steps.
7	Finally, update the sizes of B and S.

P_acc = ...    # P accumulation.
P_trans = ...  # Transfer flow between dissolved P and adsorbed stable P pool.
B = ...        # Buffer precipitation pool.
dB = ...       # Change in the B pool.
L = ...        # Adsorbed 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	P_trans > P_acc > 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	P_acc > P_trans and/or P_trans < 0, and the size of L will increase to a provisional new size L'.
6	When L' < L_min, L' will be padded to L_min, and the difference ( \|L' − L_min\| ) 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' − L_min\| − B ) will appear to be coming out of nowhere. However, since this value will be very small, it may be safely neglected.
8	B is larger than the difference, B is reduced by the amount of the difference.
9	When L' > L_max, B acts as a buffer pool for excess P.
10	When L_min < L' < L_max, B remains unchanged.
11	Update the sizes of B and L.