1"""modvege_lib.py
2
3Main ModVege functions.
4"""
5
6import numpy as np
7
8np.seterr("raise")
9
10
[docs]
11def leaf_area_index(
12 ts_vals: dict[str, float], params: dict[str, float]
13) -> float:
14 """Calculate the leaf area index (LAI)
15
16 Equation (12) in Jouven et al. (2006a)
17
18 Parameters
19 ----------
20 params : dict[str, float]
21 A dictionary containing these model parameters:
22 - ``pct_lam``: Percentage of laminae in the green vegetative (GV)
23 biomass compartment; default is 0.68 (%LAM) [dimensionless]
24 - ``sla``: Specific leaf area; default is 0.033 (SLA) [m² g⁻¹]
25 ts_vals : dict[str, float]
26 A dictionary with intermediate time series values for:
27 - ``bm_gv``: Standing biomass of the green vegetative (GV)
28 compartment (BM_GV) [kg DM ha⁻¹]
29 - ``bm_gr``: Standing biomass of the green reproductive (GR)
30 compartment (BM_GR) [kg DM ha⁻¹]
31
32 Returns
33 -------
34 float
35 Leaf area index (LAI) [dimensionless]
36 """
37
38 # use the sum of both green compartments
39 return (
40 params["sla"]
41 * (ts_vals["bm_gv"] + ts_vals["bm_gr"])
42 / 10.0
43 * params["pct_lam"]
44 )
45 # return (
46 # params["sla"] *
47 # ts_vals["bm_gv"] / 10.0 *
48 # params["pct_lam"]
49 # )
50
51
[docs]
52def actual_evapotranspiration(pet: float, ts_vals: dict[str, float]) -> float:
53 """Calculate the actual evapotranspiration (AET)
54
55 AET is equivalent to potential evapotranspiration (PET) when the cover
56 intercepts approximately 0.95 of the incident photosynthetically active
57 radiation (PAR), i.e., when the leaf area index (LAI) > 3,
58 based on Johnson and Parsons (1985).
59 AET is proportional to LAI when the proportion of intercepted radiation
60 is lower than 0.95, i.e. LAI < 3.
61
62 See Equation (14) in Jouven et al. (2006a)
63
64 Parameters
65 ----------
66 pet : float
67 Potential evapotranspiration (PET) [mm]
68 ts_vals : dict[str, float]
69 A dictionary with intermediate time series values for:
70 - ``lai``: Leaf area index (LAI) [dimensionless]
71
72 Returns
73 -------
74 float
75 Actual evapotranspiration (AET) [mm]
76 """
77
78 return min(pet, pet * ts_vals["lai"] / 3.0)
79
80
[docs]
81def potential_growth(
82 par_i: float, ts_vals: dict[str, float], params: dict[str, float]
83) -> float:
84 """Calculate potential growth (PGRO)
85
86 See Equation (12) in Jouven et al. (2006a)
87
88 Based on Schapendonk et al. (1998).
89
90 The model extinction coefficient is set to a constant value of 0.6
91 according to Schapendonk et al. (1998) and Bonesmo and Bélanger (2002).
92
93 The maximum radiation use efficiency is 3 g DM MJ⁻¹ based on
94 Schapendonk et al. (1998).
95
96 Parameters
97 ----------
98 par_i : float
99 Incident photosynthetically active radiation (PAR_i) [MJ m⁻²]
100 params : dict[str, float]
101 A dictionary containing model parameters:
102 - ``rue_max``: Maximum radiation use efficiency (RUE_max); default
103 is 3.0 [g DM MJ⁻¹]
104 ts_vals : dict[str, float]
105 A dictionary with intermediate time series values for:
106 - ``lai``: Leaf area index (LAI) [dimensionless]
107
108 Returns
109 -------
110 float
111 Potential growth (PGRO) [kg DM ha⁻¹]
112 """
113
114 return (
115 par_i
116 * params["rue_max"]
117 * (1.0 - np.exp(-0.6 * ts_vals["lai"]))
118 * 10.0
119 )
120
121
[docs]
122def par_function(par_i: float) -> float:
123 """Incident photosynthetically active radiation
124
125 Incident photosynthetically active radiation (PAR_i) function
126 (f(PAR_i)) needed to calculate the environmental limitation of growth
127 (ENV).
128
129 The definition has been derived from Schapendonk et al. (1998).
130 This function accounts for the decrease in radiation use efficiency (RUE)
131 at light intensities higher than 5 MJ m⁻².
132
133 See Figure 2(a), Equation (13), and the section on "Growth functions" in
134 Jouven et al. (2006a).
135
136 Parameters
137 ----------
138 par_i : float
139 Incident photosynthetically active radiation (PAR_i) [MJ m⁻²]
140
141 Returns
142 -------
143 float
144 PAR_i function (f(PAR_i)) [dimensionless]
145 """
146
147 if par_i < 5.0:
148 val = 1.0
149 else:
150 # linear gradient
151 gradient = 1.0 / (5.0 - (25.0 + 30.0) / 2.0)
152 intercept = 1.0 - gradient * 5.0
153 val = max(gradient * par_i + intercept, 0.0)
154 return val
155
156
[docs]
157def sum_of_temperatures(
158 params: dict[str, float],
159 ts_vals: dict[str, float],
160 t_ts: list[float],
161 day: int,
162) -> float:
163 """Sum of temperature threshold
164
165 Return the sum of temperatures for each day of the year above the minimum
166 temperature for growth (ST)
167
168 Parameters
169 ----------
170 t_ts : list[float]
171 Temperature (T) field of the input time series data (temperature
172 should be in °C)
173 day : int
174 Day number
175 params : dict[str, float]
176 A dictionary containing model parameters:
177 - ``t_0``: Minimum temperature for growth (T₀); default is 4 [°C]
178 ts_vals : dict[str, float]
179 A dictionary with intermediate time series values for:
180 - ``st``: Sum of temperatures value for the previous data row (ST)
181 [°C d]
182
183 Returns
184 -------
185 float
186 Sum of temperatures above T₀ corresponding to each day of the year
187 (ST) [°C d]
188
189 Notes
190 -----
191 - Degree days are measures of how cold or warm a location is
192 - A degree day compares the mean (the average of the high and low)
193 outdoor temperatures recorded for a location to a
194 standard temperature
195 - Also known as heat units or thermal units
196 - All species of plants have a cutoff temperature below which no
197 development occurs (developmental threshold)
198 - Degree days are accumulated whenever the temperature exceeds the
199 predetermined developmental threshold
200 - Calculate degree days by subtracting the developmental threshold from
201 the average daily temperature
202 - If the average degree day value for a given day is less than zero, just
203 record zero, not a negative number
204
205 References
206 ----------
207 - https://hort.extension.wisc.edu/articles/degree-day-calculation/
208 - https://www.eia.gov/energyexplained/units-and-calculators/degree-days.php
209 """
210
211 for i in range(day):
212 if t_ts[i] > params["t_0"]:
213 val = ts_vals["st"] + t_ts[i] - params["t_0"]
214 else:
215 val = ts_vals["st"]
216 return val
217
218
[docs]
219def ten_day_moving_avg_temperature(
220 day: int, t_ts: list[float], t_init: list[float] = None
221) -> float:
222 """Calculate the 10-d moving average temperature.
223
224 See sec. "Growth functions", par. above Equation (13) in Jouven et al.
225 (2006a).
226
227 Parameters
228 ----------
229 t_ts : list[float]
230 Temperature (T) field of the input time series data (temperature
231 should be in °C)
232 day : int
233 Day number
234
235 Returns
236 -------
237 float
238 10-d moving average temperature [°C]
239 """
240
241 if day < 10:
242 if t_init is None:
243 # Note: USING THE TEMP, NOT 10-d MOVING AVG!
244 val = t_ts[day - 1]
245 else:
246 t_list = list(t_init) + list(t_ts)
247 val = np.mean([t_list[(day + 8) - j] for j in range(9, -1, -1)])
248 else:
249 val = np.mean([t_ts[(day - 1) - j] for j in range(9, -1, -1)])
250 return val
251
252
[docs]
253def temperature_function(
254 ts_vals: dict[str, float], params: dict[str, float]
255) -> float:
256 """Temperature function, f(T)
257
258 See Figure 2(b) of Jouven et al. (2006a) and the accompanying text for
259 more info; f(T) has been derived based on Schapendonk et al. (1998)
260
261 Assume no growth takes place after a maximum temperature
262
263 Parameters
264 ----------
265 ts_vals : dict[str, float]
266 A dictionary with intermediate time series values for:
267 - ``t_m10``: 10-d moving average temperature [°C]
268 params : dict[str, float]
269 A dictionary containing these model parameters:
270 - ``t_0``: Minimum temperature for growth (T₀); default is 4 [°C]
271 - ``t_1``: Minimum temperature for optimal growth; default is 10
272 [°C]
273 - ``t_2``: Maximum temperature for optimal growth; default is 20
274 [°C]
275 - ``t_max``: Maximum temperature for growth; default is 40 [°C]
276
277 Returns
278 -------
279 float
280 Temperature function (f(T)) [dimensionless]
281 """
282
283 if (
284 ts_vals["t_m10"] <= params["t_0"]
285 or ts_vals["t_m10"] >= params["t_max"]
286 ):
287 val = 0.0
288 elif params["t_0"] < ts_vals["t_m10"] < params["t_1"]:
289 # linear relationship
290 gradient = 1.0 / (params["t_1"] - params["t_0"])
291 intercept = 1.0 - gradient * params["t_1"]
292 val = gradient * ts_vals["t_m10"] + intercept
293 elif params["t_1"] <= ts_vals["t_m10"] <= params["t_2"]:
294 val = 1.0
295 elif params["t_2"] < ts_vals["t_m10"] < params["t_max"]:
296 # linear relationship
297 gradient = 1.0 / (params["t_2"] - params["t_max"])
298 intercept = 1.0 - gradient * params["t_2"]
299 val = gradient * ts_vals["t_m10"] + intercept
300 return val
301
302
[docs]
303def seasonal_effect(
304 # ts_vals: dict[str, float],
305 params: dict[str, float]
306) -> float:
307 """Seasonal effect
308
309 Calculate seasonal effect (SEA) on growth, driven by the sum of
310 temperatures.
311
312 Note: A constant value (the average of minSEA and maxSEA) is used for
313 the SEA if the sum of temperatures at the beginning of the reproductive
314 period is lower than the sum of temperatures at the onset of reproductive
315 growth.
316
317 SEA > 1 indicates above-ground stimulation by mobilisation of reserves;
318 SEA < 1 indicates growth limitation by storage of reserves
319
320 SEA = minSEA when ST < 200°C d, then increases and reaches maxSEA when
321 (ST₁ - 200) < ST < (ST₁ - 100) (ST = ST₁ at the beginning of the
322 reproductive period); during summer, SEA decreases, returning to minSEA at
323 ST₂ (ST = ST₂ at the end of the reproductive period)
324
325 "T-sum 200 date" had its origins in the Netherlands and reflects the
326 amount of heat absorbed and hence the amount of energy available to
327 promote grass growth, and is used by farmers as an indication of
328 conditions suitable for nitrogen application to grass swards
329 (Collins and Cummins, 1998).
330
331 See Figure 3 of Jouven et al. (2006a) and the accompanying paragraphs for
332 more info
333
334 minSEA and maxSEA are functional traits arranged symmetrically around 1:
335 (minSEA + maxSEA) / 2 = 1
336
337 Parameters
338 ----------
339 params : dict[str, float]
340 A dictionary containing these model parameters:
341 - ``max_sea``: Maximum seasonal effect (maxSEA); default is 1.2
342 [dimensionless]
343 - ``min_sea``: Minimum seasonal effect (minSEA); default is 0.8
344 [dimensionless]
345 - ``st_1``: Sum of temperatures at the beginning of the
346 reproductive period (ST₁) [°C d]
347 - ``st_2``: Sum of temperatures at the end of the reproductive
348 period (ST₂) [°C d]
349 ts_vals : dict[str, float]
350 A dictionary with intermediate time series values for:
351 - ``st``: Sum of temperatures (ST) [°C d]
352
353 Returns
354 -------
355 float
356 Seasonal effect [dimensionless]
357 """
358
359 # if params["st_1"] <= 200.0:
360 # # use a constant value
361 # val = np.mean([params["max_sea"], params["min_sea"]])
362 # else:
363 # if ts_vals["st"] <= 200.0 or ts_vals["st"] >= params["st_2"]:
364 # val = params["min_sea"]
365 # elif (
366 # params["st_1"] - 200.0 <=
367 # ts_vals["st"] <=
368 # params["st_1"] - 100.0
369 # ):
370 # val = params["max_sea"]
371 # elif 200.0 < ts_vals["st"] < (params["st_1"] - 200.0):
372 # # assume SEA increases linearly from minSEA at the onset of
373 # # growth to maxSEA
374 # gradient = (
375 # (params["max_sea"] - params["min_sea"]) /
376 # ((params["st_1"] - 200.0) - 200.0)
377 # )
378 # intercept = params["min_sea"] - gradient * 200.0
379 # val = max(
380 # gradient * ts_vals["st"] + intercept, params["min_sea"]
381 # )
382 # elif params["st_1"] - 100.0 < ts_vals["st"] < params["st_2"]:
383 # # SEA decreases linearly from maxSEA to minSEA at ST_2
384 # gradient = (
385 # (params["max_sea"] - params["min_sea"]) /
386 # ((params["st_1"] - 100.0) - params["st_2"])
387 # )
388 # intercept = params["min_sea"] - gradient * params["st_2"]
389 # val = max(
390 # gradient * ts_vals["st"] + intercept, params["min_sea"]
391 # )
392 # return val
393 return np.mean([params["max_sea"], params["min_sea"]])
394
395
[docs]
396def water_reserves(
397 ts_vals: dict[str, float], params: dict[str, float], precipitation: float
398) -> float:
399 """Calculate the water reserves (WR).
400
401 WR vary between zero and the soil water-holding capacity (WHC).
402 Precipitation (PP) fill the WHC, increasing WR, while actual
403 evapotranspiration (AET) empties it.
404
405 See Equation (14) in Jouven et al. (2006a).
406
407 Parameters
408 ----------
409 precipitation : float
410 Precipitation (PP) [mm]
411 ts_vals : dict[str, float]
412 A dictionary with intermediate time series values for:
413 - ``wr``: Water reserves (WR) [mm]
414 - ``aet``: Actual evapotranspiration (AET) [mm]
415 params : dict[str, float]
416 A dictionary containing model parameters:
417 - ``whc``: Soil water-holding capacity (WHC) [mm]
418
419 Returns
420 -------
421 float
422 Water reserves (WR) [mm]
423 """
424
425 return min(
426 max(0.0, ts_vals["wr"] + precipitation - ts_vals["aet"]), params["whc"]
427 )
428
429
[docs]
430def water_stress(ts_vals: dict[str, float], params: dict[str, float]) -> float:
431 """Calculate the water stress (W).
432
433 See Equation (14) in Jouven et al. (2006a)
434
435 Parameters
436 ----------
437 ts_vals : dict[str, float]
438 A dictionary with intermediate time series values for:
439 - ``wr``: Water reserves (WR) [mm]
440 params : dict[str, float]
441 A dictionary containing model parameters:
442 - ``whc``: Soil water-holding capacity (WHC) [mm]
443
444 Returns
445 -------
446 float
447 Water stress (W) [dimensionless]
448 """
449
450 try:
451 val = min(ts_vals["wr"] / params["whc"], 1.0)
452 except ZeroDivisionError:
453 val = 0.0
454 return val
455
456
[docs]
457def water_stress_function(ts_vals: dict[str, float], pet: float) -> float:
458 """Water stress function (f(W)).
459
460 See Figure 2(c) and Equation (14) of Jouven et al. (2006a).
461
462 Based on McCall and Bishop-Hurley (2003).
463
464 Parameters
465 ----------
466 ts_vals : dict[str, float]
467 A dictionary with intermediate time series values for:
468 - ``w``: Water stress (W) [dimensionless]
469 pet : float
470 Potential evapotranspiration (PET) [mm]
471
472 Returns
473 -------
474 float
475 Water stress function (f(W)) [dimensionless]
476 """
477
478 if pet < 3.8:
479 # linear gradients
480 if ts_vals["w"] < 0.2:
481 gradient = 0.8 / 0.2
482 val = gradient * ts_vals["w"]
483 elif ts_vals["w"] < 0.4:
484 gradient = (0.95 - 0.8) / (0.4 - 0.2)
485 intercept = 0.8 - gradient * 0.2
486 val = gradient * ts_vals["w"] + intercept
487 elif ts_vals["w"] < 0.6:
488 gradient = (1.0 - 0.95) / (0.6 - 0.4)
489 intercept = 1.0 - gradient * 0.6
490 val = gradient * ts_vals["w"] + intercept
491 else:
492 val = 1.0
493 elif pet <= 6.5:
494 if ts_vals["w"] < 0.2:
495 gradient = 0.4 / 0.2
496 val = gradient * ts_vals["w"]
497 elif ts_vals["w"] < 0.4:
498 gradient = (0.7 - 0.4) / (0.4 - 0.2)
499 intercept = 0.4 - gradient * 0.2
500 val = gradient * ts_vals["w"] + intercept
501 elif ts_vals["w"] < 0.6:
502 gradient = (0.9 - 0.7) / (0.6 - 0.4)
503 intercept = 0.9 - gradient * 0.6
504 val = gradient * ts_vals["w"] + intercept
505 elif ts_vals["w"] < 0.8:
506 gradient = (1.0 - 0.9) / (0.8 - 0.6)
507 intercept = 1.0 - gradient * 0.8
508 val = 0.5 * ts_vals["w"] + 0.6
509 else:
510 val = 1.0
511 else:
512 val = ts_vals["w"]
513 return val
514
515
[docs]
516def reproductive_function(
517 params: dict[str, float], ts_vals: dict[str, float]
518) -> float:
519 """Reproductive function (REP).
520
521 REP is zero when there is a cut due to grazing or harvesting.
522 REP is also zero before and after the period of reproductive growth.
523
524 See Equation (15) in Jouven et al. (2006a)
525
526 Parameters
527 ----------
528 ts_vals : dict[str, float]
529 A dictionary with intermediate time series values for:
530 - ``st``: Sum of temperatures [°C d]
531 - ``h_bm``: The total harvested biomass amount [kg DM ha⁻¹]
532 - ``i_bm``: The total ingested biomass amount [kg DM ha⁻¹]
533 params : dict[str, float]
534 A dictionary containing model parameters:
535 - ``ni``: Nitrogen nutritional index (NI) [dimensionless]
536 - ``st_1``: Sum of temperatures at the beginning of the
537 reproductive period [°C d]
538
539 Returns
540 -------
541 float
542 Reproductive function [dimensionless]
543 """
544
545 if (
546 ts_vals["st"] < params["st_1"]
547 or ts_vals["i_bm"] > 0.0
548 or ts_vals["h_bm"] > 0.0
549 ):
550 val = 0.0
551 else:
552 val = 0.25 + ((1.0 - 0.25) * (params["ni"] - 0.35)) / (1.0 - 0.35)
553 return val
554
555
[docs]
556def environmental_limitation(
557 ts_vals: dict[str, float], params: dict[str, float], par_i: float
558) -> float:
559 """Environmental limitation of growth (ENV).
560
561 See Equation (13) of Jouven et al. (2006a).
562
563 Parameters
564 ----------
565 ts_vals : dict[str, float]
566 A dictionary with intermediate time series values for:
567 - ``f_t``: temperature function (f(T)) [dimensionless]
568 - ``f_w``: Water stress function (f(W)) [dimensionless]
569 params : dict[str, float]
570 A dictionary containing model parameters:
571 - ``ni``: Nutritional index of pixel (NNI) [dimensionless]
572 par_i : float
573 Incident photosynthetically active radiation (PAR_i) [MJ m⁻²]
574
575 Returns
576 -------
577 float
578 Environmental limitation of growth (ENV) [dimensionless]
579 """
580
581 return (
582 ts_vals["f_t"]
583 * params["ni"]
584 * par_function(par_i=par_i)
585 * ts_vals["f_w"]
586 )
587
588
[docs]
589def total_growth(ts_vals: dict[str, float]) -> float:
590 """Calculate the total biomass growth (GRO)
591
592 See Equation (11) in Jouven et al. (2006a)
593
594 Parameters
595 ----------
596 ts_vals : dict[str, float]
597 A dictionary with intermediate time series values for:
598 - ``pgro``: Potential growth (PGRO) [kg DM ha⁻¹]
599 - ``env``: Environmental limitation of growth (ENV) [dimensionless]
600 - ``sea``: Seasonal effect (SEA) [dimensionless]
601
602 Returns
603 -------
604 float
605 Total biomass growth (GRO) [kg DM ha⁻¹]
606 """
607
608 return ts_vals["pgro"] * ts_vals["env"] * ts_vals["sea"]
609
610
[docs]
611def abscission(
612 ts_vals: dict[str, float], params: dict[str, float], temperature: float
613) -> dict[str, float]:
614 """Abscission biomass.
615
616 Compute abscission biomass for the dead vegetative (DV) and dead
617 reproductive (DR) compartments.
618 See Equation (18) and Figure 4(c) and (d) in in Jouven et al. (2006a).
619
620 Note that abscission only occurs when T > 0.
621
622 Parameters
623 ----------
624 params : dict[str, float]
625 A dictionary containing model parameters:
626 - ``lls``: Leaf lifespan (LLS) [500 °C d]
627 - ``kl_dv``: Basic abscission rate for the dead vegetative;
628 compartment default is 0.001 (Kl_DV) [dimensionless]
629 - ``kl_dr``: Basic abscission rate for the dead reproductive
630 compartment; default is 0.0005 (Kl_DR) [dimensionless]
631 - ``st_1``: Sum of temperatures at the beginning of the
632 reproductive period; (ST₁) [°C d]
633 - ``st_2``: Sum of temperatures at the end of the reproductive
634 period; (ST₂) [°C d]
635 temperature : float
636 Mean daily temperature (T) [°C]
637 ts_vals : dict[str, float]
638 A dictionary with intermediate time series values for:
639 - ``bm_dv``: DV biomass (BM_DV) [kg DM ha⁻¹]
640 - ``age_dv``: Age of the DV compartment (AGE_DV) [°C d]
641 - ``bm_dr``: DR biomass (BM_DR) [kg DM ha⁻¹]
642 - ``age_dr``: Age of the DR compartment (AGE_DR) [°C d]
643
644 Returns
645 -------
646 dict[str, float]
647 An updated `ts_vals` dictionary with:
648 - ``abs_dv``: Abscission of the DV biomass (ABS_DV) [kg DM ha⁻¹]
649 - ``abs_dr``: Abscission of the DR biomass (ABS_DR) [kg DM ha⁻¹]
650 """
651
652 # abscission of the DV biomass
653 if ts_vals["age_dv"] / params["lls"] < 1.0 / 3.0:
654 ts_vals["f_age_dv"] = 1.0
655 elif ts_vals["age_dv"] / params["lls"] < 2.0 / 3.0:
656 ts_vals["f_age_dv"] = 2.0
657 else:
658 ts_vals["f_age_dv"] = 3.0
659 if temperature > 0.0:
660 ts_vals["abs_dv"] = (
661 params["kl_dv"]
662 * ts_vals["bm_dv"]
663 * temperature
664 * ts_vals["f_age_dv"]
665 )
666 else:
667 ts_vals["abs_dv"] = 0.0
668
669 # abscission of the DR biomass
670 if ts_vals["age_dr"] / (params["st_2"] - params["st_1"]) < 1.0 / 3.0:
671 ts_vals["f_age_dr"] = 1.0
672 elif ts_vals["age_dr"] / (params["st_2"] - params["st_1"]) < 2.0 / 3.0:
673 ts_vals["f_age_dr"] = 2.0
674 else:
675 ts_vals["f_age_dr"] = 3.0
676 if temperature > 0.0:
677 ts_vals["abs_dr"] = (
678 params["kl_dr"]
679 * ts_vals["bm_dr"]
680 * temperature
681 * ts_vals["f_age_dr"]
682 )
683 else:
684 ts_vals["abs_dr"] = 0.0
685
686
[docs]
687def senescence(
688 ts_vals: dict[str, float], params: dict[str, float], temperature: float
689) -> dict[str, float]:
690 """Senescing biomass for the GV and GR compartments.
691
692 See Equations (16) and (17) and Figure 4(a) and (b) in Jouven et al.
693 (2006a).
694
695 No senescence occurs when T is between zero and T₀.
696 When T drops below zero, senescence is driven by freezing effects and is
697 proportional to :math:`(|T|)`.
698
699 Parameters
700 ----------
701 params : dict[str, float]
702 A dictionary containing model parameters:
703 - ``k_gv``: Basic senescence rate for the green vegetative
704 compartment; default is 0.002 (K_GV) [dimensionless]
705 - ``k_gr``: Basic senescence rate for the green reproductive
706 compartment; default is 0.001 (K_GR) [dimensionless]
707 - ``t_0``: Minimum temperature for growth; default is 4 (T₀) [°C]
708 - ``lls``: Leaf lifespan; default is 500 (LLS) [°C d]
709 - ``st_1`` : Sum of temperatures at the beginning of the
710 reproductive period; (ST₁) [°C d]
711 - ``st_2`` : Sum of temperatures at the end of the reproductive
712 period; (ST₂) [°C d]
713 ts_vals : dict[str, float]
714 A dictionary with intermediate time series values for:
715 - ``bm_gv``: GV biomass (BM_GV) [kg DM ha⁻¹]
716 - ``age_gv``: Age of the GV compartment (AGE_GV) [°C d]
717 - ``bm_gr``: Biomass available for GR (BM_GR) [kg DM ha⁻¹]
718 - ``age_gr``: Age of the GR compartment (AGE_GR) [°C d]
719 temperature : float
720 Mean daily temperature (T) [°C]
721
722 Returns
723 -------
724 dict[str, float]
725 An updated `ts_vals` dictionary with:
726 - Senescing GV biomass (SEN_GV) [kg DM ha⁻¹]
727 - Senescing GR biomass (SEN_GR) [kg DM ha⁻¹]
728 """
729
730 # senescing GV biomass
731 if ts_vals["age_gv"] / params["lls"] < 1.0 / 3.0:
732 ts_vals["f_age_gv"] = 1.0
733 elif ts_vals["age_gv"] / params["lls"] < 1.0:
734 # linear gradient
735 gradient = (3.0 - 1.0) / (1.0 - 1.0 / 3.0)
736 intercept = 3.0 - gradient * 1.0
737 ts_vals["f_age_gv"] = (
738 gradient * ts_vals["age_gv"] / params["lls"] + intercept
739 )
740 else:
741 ts_vals["f_age_gv"] = 3.0
742 if temperature > params["t_0"]:
743 ts_vals["sen_gv"] = (
744 params["k_gv"]
745 * ts_vals["bm_gv"]
746 * temperature
747 * ts_vals["f_age_gv"]
748 )
749 elif temperature < 0.0:
750 ts_vals["sen_gv"] = (
751 params["k_gv"] * ts_vals["bm_gv"] * abs(temperature)
752 )
753 else:
754 ts_vals["sen_gv"] = 0.0
755
756 # senescing GR biomass
757 if ts_vals["age_gr"] / (params["st_2"] - params["st_1"]) < 1.0 / 3.0:
758 ts_vals["f_age_gr"] = 1.0
759 elif ts_vals["age_gr"] / (params["st_2"] - params["st_1"]) < 1.0:
760 # linear gradient
761 gradient = (3.0 - 1.0) / (1.0 - 1.0 / 3.0)
762 intercept = 3.0 - gradient * 1.0
763 ts_vals["f_age_gr"] = (
764 gradient * ts_vals["age_gr"] / (params["st_2"] - params["st_1"])
765 + intercept
766 )
767 else:
768 ts_vals["f_age_gr"] = 3.0
769 if temperature > params["t_0"]:
770 ts_vals["sen_gr"] = (
771 params["k_gr"]
772 * ts_vals["bm_gr"]
773 * temperature
774 * ts_vals["f_age_gr"]
775 )
776 elif temperature < 0.0:
777 ts_vals["sen_gr"] = (
778 params["k_gr"] * ts_vals["bm_gr"] * abs(temperature)
779 )
780 else:
781 ts_vals["sen_gr"] = 0.0
782
783
[docs]
784def biomass_growth(ts_vals: dict[str, float]) -> tuple[float, float]:
785 """Calculate the growth of the GV and GR biomass compartments.
786
787 See Equations (1) and (2) in Jouven et al. (2006a)
788
789 Parameters
790 ----------
791 ts_vals : dict[str, float]
792 A dictionary with intermediate time series values for:
793 - ``gro``: Total growth (GRO) [kg DM ha⁻¹]
794 - ``rep``: Reproductive function (REP) [dimensionless]
795
796 Returns
797 -------
798 tuple[float, float]
799 A tuple of:
800 - Growth of the GV compartment (GRO_GV) [kg DM ha⁻¹]
801 - Growth of the GR compartment (GRO_GR) [kg DM ha⁻¹]
802 """
803
804 return (
805 ts_vals["gro"] * (1.0 - ts_vals["rep"]),
806 ts_vals["gro"] * ts_vals["rep"],
807 )
808
809
[docs]
810def standing_biomass(
811 ts_vals: dict[str, float], params: dict[str, float]
812) -> dict[str, float]:
813 """Update the standing biomass for each compartment.
814
815 See Equations (1), (2), (3), and (4) in Jouven et al. (2006a).
816
817 Parameters
818 ----------
819 params : dict[str, float]
820 A dictionary containing model parameters:
821 - ``sigma_gv``: Rate of biomass loss with respiration for GV
822 [dimensionless]
823 - ``sigma_gr``: Rate of biomass loss with respiration for GR
824 [dimensionless]
825 ts_vals : dict[str, float]
826 A dictionary with intermediate time series values for:
827 - ``bm_gv``: GV biomass (BM_GV) [kg DM ha⁻¹]
828 - ``sen_gv``: Senescence of GV compartment (SEN_GV) [kg DM ha⁻¹]
829 - ``bm_gr``: GR biomass (BM_GR) [kg DM ha⁻¹]
830 - ``sen_gr``: Senescence of GR compartment (SEN_GR) [kg DM ha⁻¹]
831 - ``bm_dv``: DV biomass (BM_DV) [kg DM ha⁻¹]
832 - ``abs_dv``: Abscission of the DV compartment (ABS_DV)
833 [kg DM ha⁻¹]
834 - ``bm_dr``: DR biomass (BM_DR) [kg DM ha⁻¹]
835 - ``abs_dr``: Abscission of the DR compartment (ABS_DR)
836 [kg DM ha⁻¹]
837 - ``gro``: Total growth (GRO) [kg DM ha⁻¹]
838 - ``rep``: Reproductive function (REP) [dimensionless]
839
840 Returns
841 -------
842 dict[str, float]
843 An updated `ts_vals` dictionary with:
844 - ``bm_gv``: GV biomass (BM_GV) [kg DM ha⁻¹]
845 - ``bm_gr``: GR biomass (BM_GR) [kg DM ha⁻¹]
846 - ``bm_dv``: DV biomass (BM_DV) [kg DM ha⁻¹]
847 - ``bm_dr``: DR biomass (BM_DR) [kg DM ha⁻¹]
848 """
849
850 # GV compartment
851 ts_vals["bm_gv"] = (
852 ts_vals["bm_gv"]
853 + biomass_growth(ts_vals=ts_vals)[0]
854 - ts_vals["sen_gv"]
855 )
856
857 # GR compartment
858 ts_vals["bm_gr"] = (
859 ts_vals["bm_gr"]
860 + biomass_growth(ts_vals=ts_vals)[1]
861 - ts_vals["sen_gr"]
862 )
863
864 # DV compartment
865 ts_vals["bm_dv"] = (
866 ts_vals["bm_dv"]
867 + (1.0 - params["sigma_gv"]) * ts_vals["sen_gv"]
868 - ts_vals["abs_dv"]
869 )
870
871 # DR compartment
872 ts_vals["bm_dr"] = (
873 ts_vals["bm_dr"]
874 + (1.0 - params["sigma_gr"]) * ts_vals["sen_gr"]
875 - ts_vals["abs_dr"]
876 )
877
878
[docs]
879def biomass_age(
880 temperature: float, params: dict[str, float], ts_vals: dict[str, float]
881) -> dict[str, float]:
882 """Update the age of each biomass compartment.
883
884 See Equations (5), (6), (7), and (8) in Jouven et al. (2006a).
885
886 The age of the residual biomass is increased daily by the mean daily
887 temperature, if this temperature is positive.
888
889 Parameters
890 ----------
891 params : dict[str, float]
892 A dictionary containing model parameters:
893 - ``sigma_gv``: Rate of biomass loss with respiration for GV
894 [dimensionless]
895 - ``sigma_gr``: Rate of biomass loss with respiration for GR
896 [dimensionless]
897 ts_vals : dict[str, float]
898 A dictionary with intermediate time series values for:
899 - ``bm_gv``: GV biomass (BM_GV) [kg DM ha⁻¹]
900 - ``age_gv``: Age of the GV compartment (AGE_GV) [°C d]
901 - ``sen_gv``: Senescence of GV compartment (SEN_GV) [kg DM ha⁻¹]
902 - ``bm_gr``: GR biomass (BM_GR) [kg DM ha⁻¹]
903 - ``age_gr``: Age of the GR compartment (AGE_GR) [°C d]
904 - ``sen_gr``: Senescence of GR compartment (SEN_GR) [kg DM ha⁻¹]
905 - ``bm_dv``: DV biomass (BM_DV) [kg DM ha⁻¹]
906 - ``age_dv``: Age of the DV compartment (AGE_DV) [°C d]
907 - ``abs_dv``: Abscission of the DV compartment (ABS_DV)
908 [kg DM ha⁻¹]
909 - ``bm_dr``: DR biomass (BM_DR) [kg DM ha⁻¹]
910 - ``age_dr``: Age of the DR compartment (AGE_DR) [°C d]
911 - ``abs_dr``: Abscission of the DR compartment (ABS_DR)
912 [kg DM ha⁻¹]
913 - ``gro``: Total growth (GRO) [kg DM ha⁻¹]
914 - ``rep``: Reproductive function (REP) [dimensionless]
915 temperature : float
916 Mean daily temperature (T) [°C]
917
918 Returns
919 -------
920 dict[str, float]
921 An updated `ts_vals` dictionary with:
922 - ``age_gv``: Age of the GV compartment (AGE_GV) [°C d]
923 - ``age_gr``: Age of the GR compartment (AGE_GR) [°C d]
924 - ``age_dv``: Age of the DV compartment (AGE_DV) [°C d]
925 - ``age_dr``: Age of the DR compartment (AGE_DR) [°C d]
926 """
927
928 # GV compartment
929 if temperature > 0.0 and ts_vals["bm_gv"] > 0.0:
930 ts_vals["age_gv"] = (
931 (ts_vals["bm_gv"] - ts_vals["sen_gv"])
932 / (
933 ts_vals["bm_gv"]
934 - ts_vals["sen_gv"]
935 + biomass_growth(ts_vals=ts_vals)[0]
936 )
937 * (ts_vals["age_gv"] + temperature)
938 )
939
940 # GR compartment
941 if temperature > 0.0 and ts_vals["bm_gr"] > 0.0:
942 ts_vals["age_gr"] = (
943 (ts_vals["bm_gr"] - ts_vals["sen_gr"])
944 / (
945 ts_vals["bm_gr"]
946 - ts_vals["sen_gr"]
947 + biomass_growth(ts_vals=ts_vals)[1]
948 )
949 * (ts_vals["age_gr"] + temperature)
950 )
951
952 # DV compartment
953 if temperature > 0.0 and ts_vals["bm_dv"] > 0.0:
954 ts_vals["age_dv"] = (
955 (ts_vals["bm_dv"] - ts_vals["abs_dv"])
956 / (
957 ts_vals["bm_dv"]
958 - ts_vals["abs_dv"]
959 + (1.0 - params["sigma_gv"]) * ts_vals["sen_gv"]
960 )
961 * (ts_vals["age_dv"] + temperature)
962 )
963
964 # DR compartment
965 if temperature > 0.0 and ts_vals["bm_dr"] > 0.0:
966 ts_vals["age_dr"] = (
967 (ts_vals["bm_dr"] - ts_vals["abs_dr"])
968 / (
969 ts_vals["bm_dr"]
970 - ts_vals["abs_dr"]
971 + (1.0 - params["sigma_gr"]) * ts_vals["sen_gr"]
972 )
973 * (ts_vals["age_dr"] + temperature)
974 )