Coverage for wsimod\nodes\land.py: 10%
780 statements
« prev ^ index » next coverage.py v7.6.1, created at 2024-10-24 11:16 +0100
1# -*- coding: utf-8 -*-
2"""Created on Fri May 20 08:58:58 2022.
4@author: Barney
5"""
6import sys
7from bisect import bisect_left
8from math import exp, log, log10, sin
9from typing import Any, Dict
11from wsimod.core import constants
12from wsimod.nodes.nodes import Node
13from wsimod.nodes.nutrient_pool import NutrientPool
14from wsimod.nodes.tanks import DecayTank, ResidenceTank
17class Land(Node):
18 """"""
20 def __init__(
21 self,
22 name,
23 subsurface_residence_time=5,
24 percolation_residence_time=50,
25 surface_residence_time=1,
26 surfaces=[],
27 data_input_dict={},
28 ):
29 """An extensive node class that represents land processes (agriculture, soil,
30 subsurface flow, rural runoff, urban drainage, pollution deposition). The
31 expected use is that each distinctive type of land cover (different crop types,
32 gardens, forests, impervious urban drainage, etc.) each have a Surface object,
33 which is a subclass of Tank. The land node will iterate over its surfaces each
34 timestep, which will generally (except in the case of an impervious surface)
35 send water to three common Tanks: surface flow, subsurface flow and percolation.
36 These tanks will then send flows to rivers or groundwater.
38 (See wsimod/nodes/land.py/Surface and subclasses for currently available
39 surfaces)
41 Args:
42 name (str): node name. subsurface_residence_time (float, optional):
43 Residence time for
44 subsurface flow (see nodes.py/ResidenceTank). Defaults to 5.
45 percolation_residence_time (int, optional): Residence time for
46 percolation flow (see nodes.py/ResidenceTank). Defaults to 50.
47 surface_residence_time (int, optional): Residence time for surface flow
48 (see nodes.py/ResidenceTank). Defaults to 1.
49 surfaces (list, optional): list of dicts where each dict describes the
50 parameters of each surface in the Land node. Each dict also contains an
51 entry under 'type_' which describes which subclass of surface to use.
52 Defaults to [].
53 data_input_dict (dict, optional): Dictionary of data inputs relevant for
54 the node (generally, et0, precipitation and temperature). Keys are
55 tuples where first value is the name of the variable to read from the
56 dict and the second value is the time. Defaults to {}.
58 Functions intended to call in orchestration:
59 run apply_irrigation (if used)
61 Key assumptions:
62 - Percolation, surface runoff, and subsurface runoff, can be described with
63 a residence-time method.
64 - Flows to percolation, surface runoff, and subsurface runoff are
65 generated by different hydrological response units (subclasses of
66 `land.py/Surface`), but aggregated for a given land node.
67 - Flows to percolation are distributed to `storage.py/Groundwater`
68 nodes while surface/subsurface runoff to `nodes.py/Node` or
69 `storage.py/River` nodes.
70 - Input data associated with the land node (precipitation,
71 temperature, evapotranspiartion) are the same for every surface.
72 - Water received from `sewer.py/Sewer` objects is sent to the first
73 `land.py/ImperviousSurface` in the surfaces list.
75 Input data and parameter requirements:
76 - Precipitation and evapotranspiration are in the `data_input_dict`
77 at the model timestep. _Units_: metres/timestep
78 - Temperature in the `data_input_dict` at the model timestep.
79 _Units_: C
80 - Residence time of surface, subsurface and percolation flows.
81 _Units_: number of timesteps
82 """
83 # Assign parameters
84 self.subsurface_residence_time = subsurface_residence_time
85 self.percolation_residence_time = percolation_residence_time
86 self.surface_residence_time = surface_residence_time
87 self.data_input_dict = data_input_dict
89 super().__init__(name, data_input_dict=data_input_dict)
91 # This could be a deny but then you would have to know in advance whether a
92 # demand node has any gardening or not
93 self.push_check_handler[("Demand", "Garden")] = lambda x: self.empty_vqip()
94 self.push_set_handler[("Demand", "Garden")] = lambda x: self.empty_vqip()
96 # Create surfaces
97 self.irrigation_functions = [lambda: None]
99 surfaces_ = surfaces.copy()
100 surfaces = []
101 for surface in surfaces_:
102 # Assign parent (for data reading and to determine where to send flows to)
103 surface["parent"] = self
105 # Get surface type
106 type_ = surface["type_"]
107 del surface["type_"]
109 # Instantiate surface and store in list of surfaces
110 surfaces.append(getattr(sys.modules[__name__], type_)(**surface))
112 # Assign ds (mass balance checking)
113 self.mass_balance_ds.append(surfaces[-1].ds)
115 # Assign any irrigation functions
116 if isinstance(surfaces[-1], IrrigationSurface):
117 # TODO, this should probably be done in the Surface initialisation
118 self.irrigation_functions.append(surfaces[-1].irrigate)
120 # Assign garden surface functions
121 if isinstance(surfaces[-1], GardenSurface):
122 # TODO, this should probably be done in the Surface initialisation
123 self.push_check_handler[("Demand", "Garden")] = surfaces[
124 -1
125 ].calculate_irrigation_demand
126 self.push_set_handler[("Demand", "Garden")] = surfaces[
127 -1
128 ].receive_irrigation_demand
130 # Update handlers
131 self.push_set_handler["default"] = self.push_set_deny
132 self.push_check_handler["default"] = self.push_check_deny
133 self.push_set_handler["Sewer"] = self.push_set_sewer
135 # Create subsurface runoff, surface runoff and percolation tanks Can also do as
136 # timearea if this seems dodge (that is how it is done in IHACRES) TODO should
137 # these be decayresidencetanks?
138 self.subsurface_runoff = ResidenceTank(
139 residence_time=self.subsurface_residence_time,
140 capacity=constants.UNBOUNDED_CAPACITY,
141 )
142 self.percolation = ResidenceTank(
143 residence_time=self.percolation_residence_time,
144 capacity=constants.UNBOUNDED_CAPACITY,
145 )
146 self.surface_runoff = ResidenceTank(
147 residence_time=self.surface_residence_time,
148 capacity=constants.UNBOUNDED_CAPACITY,
149 )
151 # Store surfaces
152 self.surfaces = surfaces
154 # Mass balance checkign vqips and functions
155 self.running_inflow_mb = self.empty_vqip()
156 self.running_outflow_mb = self.empty_vqip()
158 self.mass_balance_in.append(lambda: self.running_inflow_mb)
159 self.mass_balance_out.append(lambda: self.running_outflow_mb)
160 self.mass_balance_ds.append(self.surface_runoff.ds)
161 self.mass_balance_ds.append(self.subsurface_runoff.ds)
162 self.mass_balance_ds.append(self.percolation.ds)
164 def apply_overrides(self, overrides=Dict[str, Any]):
165 """Apply overrides to the Land.
167 Enables a user to override any parameter of the residence_time and update
168 the residence_tank accordingly.
170 Args:
171 overrides (Dict[str, Any]): Dict describing which parameters should
172 be overridden (keys) and new values (values). Defaults to {}.
173 """
174 self.surface_residence_time = overrides.pop(
175 "surface_residence_time", self.surface_residence_time
176 )
177 self.subsurface_residence_time = overrides.pop(
178 "subsurface_residence_time", self.subsurface_residence_time
179 )
180 self.percolation_residence_time = overrides.pop(
181 "percolation_residence_time", self.percolation_residence_time
182 )
183 self.surface_runoff.residence_time = self.surface_residence_time
184 self.subsurface_runoff.residence_time = self.subsurface_residence_time
185 self.percolation.residence_time = self.percolation_residence_time
186 super().apply_overrides(overrides)
188 def apply_irrigation(self):
189 """Iterate over any irrigation functions (needs further testing..
191 maybe).
192 """
193 for f in self.irrigation_functions:
194 f()
196 def run(self):
197 """Call the run function in all surfaces, update surface/subsurface/ percolation
198 tanks, discharge to rivers/groundwater."""
199 # Run all surfaces
200 for surface in self.surfaces:
201 surface.run()
203 # Apply residence time to percolation
204 percolation = self.percolation.pull_outflow()
206 # Distribute percolation
207 reply = self.push_distributed(percolation, of_type=["Groundwater"])
209 if reply["volume"] > constants.FLOAT_ACCURACY:
210 # Update percolation 'tank'
211 _ = self.percolation.push_storage(reply, force=True)
213 # Apply residence time to subsurface/surface runoff
214 surface_runoff = self.surface_runoff.pull_outflow()
215 subsurface_runoff = self.subsurface_runoff.pull_outflow()
217 # Total runoff
218 total_runoff = self.sum_vqip(surface_runoff, subsurface_runoff)
219 if total_runoff["volume"] > 0:
220 # Send to rivers (or nodes, which are assumed to be junctions)
221 reply = self.push_distributed(total_runoff, of_type=["River", "Node"])
223 # Redistribute total_runoff not sent
224 if reply["volume"] > 0:
225 reply_surface = self.v_change_vqip(
226 reply,
227 reply["volume"] * surface_runoff["volume"] / total_runoff["volume"],
228 )
229 reply_subsurface = self.v_change_vqip(
230 reply,
231 reply["volume"]
232 * subsurface_runoff["volume"]
233 / total_runoff["volume"],
234 )
236 # Update surface/subsurface runoff 'tanks'
237 if reply_surface["volume"] > 0:
238 self.surface_runoff.push_storage(reply_surface, force=True)
239 if reply_subsurface["volume"] > 0:
240 self.subsurface_runoff.push_storage(reply_subsurface, force=True)
242 def push_set_sewer(self, vqip):
243 """Receive water from a sewer and send it to the first ImperviousSurface in
244 surfaces.
246 Args:
247 vqip (dict): A VQIP amount to be sent to the impervious surface
249 Returns:
250 vqip (dict): A VQIP amount of water that was not received
251 """
252 # TODO currently just push to the first impervious surface... not sure if people
253 # will be having multiple impervious surfaces. If people would be only having
254 # one then it would make sense to store as a parameter... this is probably fine
255 # for now
256 for surface in self.surfaces:
257 if isinstance(surface, ImperviousSurface):
258 vqip = self.surface.push_storage(vqip, force=True)
259 break
260 return vqip
262 def end_timestep(self):
263 """Update mass balance and end timestep of all tanks (and surfaces)."""
264 self.running_inflow_mb = self.empty_vqip()
265 self.running_outflow_mb = self.empty_vqip()
266 for tanks in self.surfaces + [
267 self.surface_runoff,
268 self.subsurface_runoff,
269 self.percolation,
270 ]:
271 tanks.end_timestep()
273 def get_surface(self, surface_):
274 """Return a surface from the list of surfaces by the 'surface' entry in the
275 surface. I.e., the name of the surface.
277 Args:
278 surface_ (str): Name of the surface
280 Returns:
281 surface (Surface): The first surface that matches the name
282 """
283 for surface in self.surfaces:
284 if surface.surface == surface_:
285 return surface
286 return None
288 def reinit(self):
289 """"""
290 self.end_timestep()
291 for surface in self.surfaces + [
292 self.surface_runoff,
293 self.subsurface_runoff,
294 self.percolation,
295 ]:
296 surface.reinit()
299class Surface(DecayTank):
300 """"""
302 def __init__(
303 self,
304 surface="",
305 area=1,
306 depth=1,
307 data_input_dict={},
308 pollutant_load={},
309 **kwargs,
310 ):
311 """A subclass of DecayTank. Each Surface is anticipated to represent a different
312 land cover type of a Land node. Besides functioning as a Tank, Surfaces have
313 three lists of functions (inflows, processes and outflows) where behaviour can
314 be added by appending new functions. We anticipate that customised surfaces
315 should be a subclass of Surface or its subclasses and add functions to these
316 lists. These lists are executed (inflows first, then processes, then outflows)
317 in the run function, which is called by the run function in Land. The lists must
318 return any model inflows or outflows as a VQIP for mass balance checking.
320 If a user wishes the DecayTank portion to be active, then can provide 'decays',
321 which are passed upwards (see wsimod/core/core.py/DecayObj for documentation)
323 Args:
324 surface (str, optional): String description of the surface type. Doesn't
325 serve a modelling purpose, just used for user reference. Defaults to ''.
326 area (float, optional): Area of surface. Defaults to 1. depth (float,
327 optional): Depth of tank (this has different physical
328 implications for different subclasses). Defaults to 1.
329 data_input_dict (dict, optional): Dictionary of data inputs relevant for
330 the surface (generally, deposition). Keys are tuples where first value
331 is the name of the variable to read from the dict and the second value
332 is the time. Note that this input should be specific to the surface, and
333 is not intended to be the same data input as for the land node. Also
334 note that with each surface having its own timeseries of data inputs,
335 this can take up a lot of memory, thus the default behavior is to have
336 this as monthly data where the time variable is a monthyear. Defaults to
337 {}.
338 pollutant_load (dict, optional): A dict of different pollutant amounts that
339 are deposited on the surface (units are mass per area per timestep).
340 Defaults to {}.
342 Key assumptions:
343 - Generic `Surface` that reads data and can apply simple forms of pollution
344 deposition.
345 - Formulated as a `Tank` object.
346 - Ammonia->Nitrite->Nitrate decay takes place if parameters describing this
347 process are provided in `decays` (see `core.py/DecayObj` for
348 transformation details).
350 Input data and parameter requirements:
351 - `data_input_dict` can contain a variety of pollutant deposition data.
352 `srp-dry` describes phosphate. `noy-dry` describes nitrogen as nitrates.
353 `nhx-dry` describes nitrogen as ammonia. `srp/noy/ nhx-wet` can also be
354 used to specify wet deposition. _Units_: kg/m2/timestep (data is read at
355 a monthly timestep)
356 """
357 # Assign parameters
358 self.depth = depth
359 self.data_input_dict = data_input_dict
360 self.surface = surface
361 self.pollutant_load = pollutant_load
362 # TODO this is a decaytank but growing surfaces don't have decay parameters...
363 # is it a problem.. we don't even take decays as an explicit argument and insert
364 # them in kwargs..
365 capacity = area * depth
366 # Parameters
367 super().__init__(capacity=capacity, area=area, **kwargs)
369 # Populate function lists TODO.. not sure why I have deposition but no
370 # precipitation here
371 if "nhx-dry" in set(x[0] for x in data_input_dict.keys()):
372 self.inflows = [self.atmospheric_deposition, self.precipitation_deposition]
373 else:
374 self.inflows = []
375 if len(self.pollutant_load) > 0:
376 self.inflows.append(self.simple_deposition)
377 self.processes = []
378 self.outflows = []
380 def apply_overrides(self, overrides=Dict[str, Any]):
381 """Override parameters.
383 Enables a user to override any of the following parameters:
384 area and depth (both will update the capacity), pollutant_load (the
385 entire dict does not need to be redefined, only changed values need to
386 be included).
388 Args:
389 overrides (Dict[str, Any]): Dict describing which parameters should
390 be overridden (keys) and new values (values). Defaults to {}.
391 """
392 self.surface = overrides.pop("surface", self.surface)
393 self.pollutant_load.update(overrides.pop("pollutant_load", {}))
395 self.area = overrides.pop("area", self.area)
396 self.depth = overrides.pop("depth", self.depth)
397 self.capacity = self.area * self.depth
399 if "capacity" in overrides.keys():
400 overrides.pop("capacity")
401 print(
402 "Warning: specifying capacity is depreciated in overrides for surface, \
403 please specify depth and area instead. capacity override value has been ignored"
404 )
406 # overrides data_input_dict
407 from wsimod.orchestration.model import read_csv
409 content = overrides.pop("data_input_dict", self.data_input_dict)
410 if isinstance(content, str):
411 self.data_input_dict = read_csv(content)
412 elif isinstance(content, dict):
413 self.data_input_dict = content
414 else:
415 raise ValueError(
416 f"{content.__class__} is not a recognised format for data_input_dict"
417 )
419 super().apply_overrides(overrides)
421 def run(self):
422 """Call run function (called from Land node)."""
423 if "nitrite" in constants.POLLUTANTS:
424 # Assume that if nitrite is modelled then nitrification is also modelled You
425 # will need ammonia->nitrite->nitrate decay to accurate simulate ammonia
426 # Thus, these decays are simulated here
428 # NOTE decay in a decaytank happens at start of timestep (confusingly) in
429 # the end_timestep function
430 self.storage["nitrate"] += self.total_decayed["nitrite"]
431 self.parent.running_inflow_mb["nitrate"] += self.total_decayed["nitrite"]
433 # Decayed ammonia becomes nitrite
434 self.storage["nitrite"] += self.total_decayed["ammonia"]
435 self.parent.running_inflow_mb["nitrite"] += self.total_decayed["ammonia"]
437 for f in self.inflows + self.processes + self.outflows:
438 # Iterate over function lists, updating mass balance
439 in_, out_ = f()
440 self.parent.running_inflow_mb = self.sum_vqip(
441 self.parent.running_inflow_mb, in_
442 )
443 self.parent.running_outflow_mb = self.sum_vqip(
444 self.parent.running_outflow_mb, out_
445 )
447 def get_data_input(self, var):
448 """Read data input from parent Land node (i.e., for precipitation/et0/temp).
450 Args:
451 var (str): Name of variable
453 Returns:
454 Data read
455 """
456 return self.parent.get_data_input(var)
458 def get_data_input_surface(self, var):
459 """Read data input from this surface's data_input_dict.
461 Args:
462 var (str): Name of variable
464 Returns:
465 Data read
466 """
467 return self.data_input_dict[(var, self.parent.monthyear)]
469 def dry_deposition_to_tank(self, vqip):
470 """Generic function for allocating dry pollution deposition to the surface.
471 Simply sends the pollution into the tank (some subclasses overwrite this
472 behaviour).
474 Args:
475 vqip (dict): A VQIP amount of dry deposition to send to tank
477 Returns:
478 vqip (dict): A VQIP amount of dry deposition that entered the tank (used
479 for mass balance checking)
480 """
481 # Default behaviour is to just enter the tank
482 _ = self.push_storage(vqip, force=True)
483 return vqip
485 def wet_deposition_to_tank(self, vqip):
486 """Generic function for allocating wet pollution deposition to the surface.
487 Simply sends the pollution into the tank (some subclasses overwrite this
488 behaviour).
490 Args:
491 vqip (dict): A VQIP amount of wet deposition to send to tank
493 Returns:
494 vqip (dict): A VQIP amount of wet deposition that entered the tank (used
495 for mass balance checking)
496 """
497 # Default behaviour is to just enter the tank
498 _ = self.push_storage(vqip, force=True)
499 return vqip
501 def simple_deposition(self):
502 """Inflow function to cause simple pollution deposition to occur, updating the
503 surface tank.
505 Returns:
506 (tuple): A tuple containing a VQIP amount for model inputs and outputs
507 for mass balance checking.
508 """
509 pollution = self.empty_vqip()
511 # Scale by area
512 for pol, item in self.pollutant_load.items():
513 if pol in constants.ADDITIVE_POLLUTANTS:
514 pollution[pol] = item * self.area
515 else:
516 pollution[pol] = item
517 pollution["volume"] = 0
519 # Update tank
520 _ = self.push_storage(pollution, force=True)
522 return (pollution, self.empty_vqip())
524 def atmospheric_deposition(self):
525 """Inflow function to cause dry atmospheric deposition to occur, updating the
526 surface tank.
528 Returns:
529 (tuple): A tuple containing a VQIP amount for model inputs and outputs
530 for mass balance checking.
531 """
532 # TODO double check units in preprocessing - is weight of N or weight of
533 # NHX/noy?
535 # Read data and scale
536 nhx = self.get_data_input_surface("nhx-dry") * self.area
537 noy = self.get_data_input_surface("noy-dry") * self.area
538 srp = self.get_data_input_surface("srp-dry") * self.area
540 # Assign pollutants
541 vqip = self.empty_vqip()
542 vqip["ammonia"] = nhx
543 vqip["nitrate"] = noy
544 vqip["phosphate"] = srp
546 # Update tank
547 in_ = self.dry_deposition_to_tank(vqip)
549 # Return mass balance
550 return (in_, self.empty_vqip())
552 def precipitation_deposition(self):
553 """Inflow function to cause wet precipitation deposition to occur, updating the
554 surface tank.
556 Returns:
557 (tuple): A tuple containing a VQIP amount for model inputs and outputs
558 for mass balance checking.
559 """
560 # TODO double check units - is weight of N or weight of NHX/noy?
562 # Read data and scale
563 nhx = self.get_data_input_surface("nhx-wet") * self.area
564 noy = self.get_data_input_surface("noy-wet") * self.area
565 srp = self.get_data_input_surface("srp-wet") * self.area
567 # Assign pollutants
568 vqip = self.empty_vqip()
569 vqip["ammonia"] = nhx
570 vqip["nitrate"] = noy
571 vqip["phosphate"] = srp
573 # Update tank
574 in_ = self.wet_deposition_to_tank(vqip)
576 # Return mass balance
577 return (in_, self.empty_vqip())
580class ImperviousSurface(Surface):
581 """"""
583 def __init__(self, pore_depth=0, et0_to_e=1, **kwargs):
584 """A surface to represent impervious surfaces that drain to storm sewers. Runoff
585 is generated by the surface tank overflowing, if a user wants all precipitation
586 to immediately go to runoff then they should reduce 'pore_depth', however
587 generally this is not what happens and a small (a few mm) depth should be
588 assigned to the tank. Also includes urban pollution deposition, though this will
589 only be mobilised if runoff occurs.
591 Note that the tank does not have a runoff coefficient because it doesn't make
592 sense from an integrated perspective. If a user wants to mimic runoff
593 coefficient-like behaviour, then they should reduce the ImperviousSurface tank
594 size, and increase other surfaces of the parent land node accordingly.
596 Args:
597 pore_depth (float, optional): The depth of the tank that must be exceeded
598 to generate runoff. Intended to represent the pores in ashpalt that
599 water accumulates in before flowing. Defaults to 0.
600 et0_to_e (float, optional): Multiplier applied to the parent's data
601 timeseries of et0 to determine how much evaporation takes place on the
602 ImperviousSurface. Defaults to 1.
603 """
604 # Assign parameters
605 self.et0_to_e = et0_to_e # Total evaporation
606 self.pore_depth = pore_depth
608 super().__init__(depth=pore_depth, **kwargs)
610 # Initialise state variables
611 self.evaporation = self.empty_vqip()
612 self.precipitation = self.empty_vqip()
614 # Populate function lists
615 self.inflows.append(self.precipitation_evaporation)
617 self.outflows.append(self.push_to_sewers)
619 def apply_overrides(self, overrides=Dict[str, Any]):
620 """Override parameters.
622 Enables a user to override any of the following parameters:
623 eto_to_e, pore_depth.
625 Args:
626 overrides (Dict[str, Any]): Dict describing which parameters should
627 be overridden (keys) and new values (values). Defaults to {}.
628 """
629 self.et0_to_e = overrides.pop("et0_to_e", self.et0_to_e)
630 if "depth" in overrides.keys():
631 overrides.pop("depth")
632 print(
633 "ERROR: specifying depth is depreciated in overrides for \
634 impervious surface, please specify pore_depth instead"
635 )
636 self.pore_depth = overrides.pop("pore_depth", self.pore_depth)
637 self.depth = self.pore_depth
638 self.capacity = self.area * self.depth
639 super().apply_overrides(overrides)
641 def precipitation_evaporation(self):
642 """Inflow function that is a simple rainfall-evaporation model, updating the.
644 surface tank. All precipitation that is not evaporated is forced into the tank
645 (even though some of that will later be pushed to sewers) - this enables runoff
646 to mix with the accumulated pollutants in the surface pores.
648 Returns:
649 (tuple): A tuple containing a VQIP amount for model inputs and outputs
650 for mass balance checking.
651 """
652 # Read data in length units
653 precipitation_depth = self.get_data_input("precipitation")
654 evaporation_depth = self.get_data_input("et0") * self.et0_to_e
656 if precipitation_depth < evaporation_depth:
657 # No effective precipitation
658 net_precipitation = 0
660 # Calculate how much should be evaporated from pores
661 evaporation_from_pores = evaporation_depth - precipitation_depth
663 # Scale
664 evaporation_from_pores *= self.area
666 # Pull from tank
667 evaporation_from_pores = self.evaporate(evaporation_from_pores)
669 # Scale to get actual evaporation
670 total_evaporation = evaporation_from_pores + precipitation_depth * self.area
671 else:
672 # Effective precipitation
673 net_precipitation = precipitation_depth - evaporation_depth
675 # Scale
676 net_precipitation *= self.area
677 net_precipitation = self.v_change_vqip(self.empty_vqip(), net_precipitation)
679 # Assign a temperature value TODO how hot is rain? No idea... just going to
680 # use surface air temperature
681 net_precipitation["temperature"] = self.get_data_input("temperature")
683 # Update tank
684 _ = self.push_storage(net_precipitation, force=True)
685 total_evaporation = evaporation_depth * self.area
687 # Converrt to VQIP
688 self.evaporation = self.v_change_vqip(self.empty_vqip(), total_evaporation)
689 self.precipitation = self.v_change_vqip(
690 self.empty_vqip(), precipitation_depth * self.area
691 )
693 return (self.precipitation, self.evaporation)
695 def push_to_sewers(self):
696 """Outflow function that distributes ponded water (i.e., surface runoff) to the
697 parent node's attached sewers.
699 Returns:
700 (tuple): A tuple containing a VQIP amount for model inputs and outputs
701 for mass balance checking.
702 """
703 # Get runoff
704 surface_runoff = self.pull_ponded()
706 # Distribute TODO in cwsd_partition this is done with timearea
707 reply = self.parent.push_distributed(surface_runoff, of_type=["Sewer"])
709 # Update tank (forcing, because if the water can't go to the sewer, where else
710 # can it go)
711 _ = self.push_storage(reply, force=True)
712 # TODO... possibly this could flow to attached river or land nodes.. or other
713 # surfaces? I expect this doesn't matter for large scale models.. but may need
714 # to be revisited for detailed sewer models
716 # Return empty mass balance because outflows are handled by parent
717 return (self.empty_vqip(), self.empty_vqip())
720class PerviousSurface(Surface):
721 """"""
723 def __init__(
724 self,
725 depth=0.75,
726 total_porosity=0.4,
727 field_capacity=0.3,
728 wilting_point=0.12,
729 infiltration_capacity=0.5,
730 surface_coefficient=0.05,
731 percolation_coefficient=0.75,
732 et0_coefficient=0.5,
733 ihacres_p=10,
734 **kwargs,
735 ):
736 """A generic pervious surface that represents hydrology with the IHACRES model.
738 Args:
739 depth (float, optional): Soil tank (i.e., root) depth. Defaults to 0.75.
740 total_porosity (float, optional): The total porosity IHACRES parameter
741 (i.e., defines the total porouse volume of the soil - the maximum volume
742 of soil pores can contain when saturated). Defaults to 0.4.
743 field_capacity (float, optional): The field capacity IHACRES parameter
744 (i.e., when water content in the soil tank is above this value - flows
745 of any kind can be generated). Defaults to 0.3.
746 wilting_point (float, optional): The wilting point IHACRES parameter (i.e.,
747 when water content content in the soil tank is above this value - plants
748 can uptake water and evaporation from the soil tank can occur). Defaults
749 to 0.12.
750 infiltration_capacity (float, optional): Depth of water per day that can
751 enter the soil tank. Non infiltrated water will pond and travel as
752 surface runoff from the parent Land node. Defaults to 0.5.
753 surface_coefficient (float, optional): If flow is generated, the proportion
754 of flow that goes to surface runoff. Defaults to 0.05.
755 percolation_coefficient (float, optional): If flow is generated, then the
756 proportion of water that does not go to surface runoff that goes to
757 percolation (i.e., groundwater) - the remainder goes to subsurface
758 runoff. Defaults to 0.75.
759 et0_coefficient (float, optional): Convert between the parent nodes data
760 timeseries et0 - and potential evaptranspiration per unit area for this
761 surface. Defaults to=0.5,
762 ihacres_p (float, optional): The IHACRES p parameter. Unless it is an
763 ephemeral stream this parameter probably can stay high. Defaults to 10.
765 Key assumptions:
766 - In IHACRES, the maximum infiltration per time step is controlled by an
767 infiltration capacity, beyond which the precipitation will flow directly
768 as surface runoff.
769 - Evapotranspiration and effective precipitation are calculated based on
770 soil moisture content.
771 - Effective precipitation is then divided into percolation, surface runoff,
772 and subsurface runoff by multiplying the corresponding coefficient.
773 - Percolation, surface runoff, and subsurface runoff are sent into the
774 corresponding residence tanks for rounting to downstream.
775 - The mass of pollutants in soil water tank proportionately leaves the
776 soil water tank into the routing residence tanks. Evapotranspiration can
777 only bring out water, with pollutants left in the soil tank.
779 Input data and parameter requirements:
780 - Field capacity and wilting point.
781 _Units_: -, both should in [0-1], with field capacity > wilting point
782 - Infiltration capacity.
783 _Units_: m/day
784 - Surface, percolation coefficient.
785 _Units_: -, both should in [0-1]
786 - et0 coefficient.
787 _Units_: -
788 - ihacres_p.
789 _Units_: -
790 """
791 # Assign parameters (converting field capacity and wilting point to depth)
792 self.field_capacity = field_capacity
793 self.field_capacity_m = field_capacity * depth
794 self.wilting_point = wilting_point
795 self.wilting_point_m = wilting_point * depth
796 self.infiltration_capacity = infiltration_capacity
797 self.surface_coefficient = surface_coefficient
798 self.percolation_coefficient = percolation_coefficient
799 self.et0_coefficient = et0_coefficient
800 self.ihacres_p = ihacres_p
801 self.total_porosity = total_porosity
803 # Parameters to determine how to calculate the temperature of outflowing water
804 # TODO what should these params be?
805 self.soil_temp_w_prev = 0.1 # previous timestep weighting
806 self.soil_temp_w_air = 0.6 # air temperature weighting
807 self.soil_temp_w_deep = 0.1 # deep soil temperature weighting
808 self.soil_temp_deep = 10 # deep soil temperature
810 # IHACRES is a deficit not a tank, so doesn't really have a capacity in this
811 # way... and if it did.. I don't know if it would be the root depth
812 super().__init__(depth=depth * total_porosity, **kwargs)
814 # Calculate subsurface coefficient
815 self.subsurface_coefficient = 1 - self.percolation_coefficient
817 # Initiliase state variables
818 self.infiltration_excess = self.empty_vqip()
819 self.subsurface_flow = self.empty_vqip()
820 self.percolation = self.empty_vqip()
821 self.tank_recharge = 0
822 self.evaporation = self.empty_vqip()
823 self.precipitation = self.empty_vqip()
825 # Populate function lists
826 self.inflows.append(self.ihacres) # work out runoff
828 # TODO interception if I hate myself enough?
829 self.processes.append(
830 self.calculate_soil_temperature
831 ) # Calculate soil temp + dependence factor
833 # self.processes.append(self.decay) #apply generic decay (currently handled by
834 # decaytank at end of timestep) TODO decaytank uses air temperature not soil
835 # temperature... probably need to just give it the decay function
837 self.outflows.append(self.route)
839 def apply_overrides(self, overrides=Dict[str, Any]):
840 """Override parameters.
842 Enables a user to override any of the following parameters:
843 field_capacity, wilting_point, total_porosity, infiltration_capacity,
844 surface_coefficient, percolation_coefficient, et0_coefficient, ihacres_p,
845 soil_temp_w_prev, soil_temp_w_air, soil_temp_w_deep, soil_temp_deep,
846 and the corresponding parameter values, including field_capacity_m,
847 wilting_point_m, depth, capacity, subsurface_coefficient.
849 Args:
850 overrides (Dict[str, Any]): Dict describing which parameters should
851 be overridden (keys) and new values (values). Defaults to {}.
852 """
854 self.depth /= self.total_porosity # restore the physical depth (root)
856 overwrite_params = [
857 "field_capacity",
858 "wilting_point",
859 "total_porosity",
860 "infiltration_capacity",
861 "surface_coefficient",
862 "percolation_coefficient",
863 "et0_coefficient",
864 "ihacres_p",
865 "soil_temp_w_prev",
866 "soil_temp_w_air",
867 "soil_temp_w_deep",
868 "soil_temp_deep",
869 ]
870 for param in overwrite_params:
871 setattr(self, param, overrides.pop(param, getattr(self, param)))
873 self.subsurface_coefficient = 1 - self.percolation_coefficient
875 super().apply_overrides(overrides)
876 # After the depth has been changed ...
877 self.field_capacity_m = self.field_capacity * self.depth
878 self.wilting_point_m = self.wilting_point * self.depth
879 self.depth *= self.total_porosity # update the simulation depth
880 self.capacity = self.depth * self.area
882 def get_cmd(self):
883 """Calculate moisture deficit (i.e., the tank excess converted to depth).
885 Returns:
886 (float): current moisture deficit
887 """
888 return self.get_excess()["volume"] / self.area
890 def get_smc(self):
891 """Calculate moisture content (i.e., the tank volume converted to depth).
893 Returns:
894 (float): soil moisture content
895 """
896 # Depth of soil moisture
897 return self.storage["volume"] / self.area
899 def get_climate(self):
900 """
902 Returns:
904 """
905 precipitation_depth = self.get_data_input("precipitation")
906 evaporation_depth = self.get_data_input("et0") * self.et0_coefficient
907 return precipitation_depth, evaporation_depth
909 def ihacres(self):
910 """Inflow function that runs the IHACRES model equations, updates tanks, and
911 store flows in state variables (which are later sent to the parent land node in
912 the route function).
914 Returns:
915 (tuple): A tuple containing a VQIP amount for model inputs and outputs
916 for mass balance checking.
917 """
918 # Read data (leave in depth units since that is what IHACRES equations are in)
919 precipitation_depth, evaporation_depth = self.get_climate()
920 temperature = self.get_data_input("temperature")
922 # Apply infiltration
923 infiltrated_precipitation = min(precipitation_depth, self.infiltration_capacity)
924 infiltration_excess = max(precipitation_depth - infiltrated_precipitation, 0)
926 # Get current moisture deficit
927 current_moisture_deficit_depth = self.get_cmd()
929 # IHACRES equations (we do (depth - wilting_point_m or field capacity) to
930 # convert from a deficit to storage tank)
931 evaporation = evaporation_depth * min(
932 1,
933 exp(
934 2
935 * (
936 1
937 - current_moisture_deficit_depth
938 / (self.depth - self.wilting_point_m)
939 )
940 ),
941 )
942 outflow = infiltrated_precipitation * (
943 1
944 - min(
945 1,
946 (current_moisture_deficit_depth / (self.depth - self.field_capacity_m))
947 ** self.ihacres_p,
948 )
949 )
951 # Can't evaporate more than available moisture (presumably the IHACRES equation
952 # prevents this ever being needed)
953 evaporation = min(evaporation, precipitation_depth + self.get_smc())
955 # Scale to volumes and apply proportions to work out percolation/surface
956 # runoff/subsurface runoff
957 surface = outflow * self.surface_coefficient * self.area
958 percolation = (
959 outflow
960 * (1 - self.surface_coefficient)
961 * self.percolation_coefficient
962 * self.area
963 )
964 subsurface_flow = (
965 outflow
966 * (1 - self.surface_coefficient)
967 * self.subsurface_coefficient
968 * self.area
969 )
970 tank_recharge = (infiltrated_precipitation - evaporation - outflow) * self.area
971 infiltration_excess *= self.area
972 infiltration_excess += surface
973 evaporation *= self.area
974 precipitation = precipitation_depth * self.area
976 # Mix in tank to calculate pollutant concentrations
977 total_water_passing_through_soil_tank = (
978 tank_recharge + subsurface_flow + percolation
979 )
981 if total_water_passing_through_soil_tank > 0:
982 # Net effective preipitation
983 total_water_passing_through_soil_tank = self.v_change_vqip(
984 self.empty_vqip(), total_water_passing_through_soil_tank
985 )
986 # Assign a temperature before sending into tank
987 total_water_passing_through_soil_tank["temperature"] = temperature
988 # Assign to tank
989 _ = self.push_storage(total_water_passing_through_soil_tank, force=True)
990 # Pull flows (which now have nonzero pollutant concentrations)
991 subsurface_flow = self.pull_storage({"volume": subsurface_flow})
992 percolation = self.pull_storage({"volume": percolation})
993 else:
994 # No net effective precipitation (instead evaporation occurs)
995 evap = self.evaporate(-total_water_passing_through_soil_tank)
996 subsurface_flow = self.empty_vqip()
997 percolation = self.empty_vqip()
999 if (
1000 abs(
1001 evap
1002 + infiltrated_precipitation * self.area
1003 - evaporation
1004 - infiltration_excess
1005 )
1006 > constants.FLOAT_ACCURACY
1007 ):
1008 print(
1009 "inaccurate evaporation calculation of {0}".format(
1010 abs(
1011 evap
1012 + infiltrated_precipitation * self.area
1013 - evaporation
1014 - infiltration_excess
1015 )
1016 )
1017 )
1019 # TODO saturation excess (think it should just be 'pull_ponded' presumably in
1020 # net effective precipitation? )
1022 # Convert to VQIPs
1023 infiltration_excess = self.v_change_vqip(self.empty_vqip(), infiltration_excess)
1024 infiltration_excess["temperature"] = temperature
1025 precipitation = self.v_change_vqip(self.empty_vqip(), precipitation)
1026 evaporation = self.v_change_vqip(self.empty_vqip(), evaporation)
1028 # Track flows (these are sent onwards in the route function)
1029 self.infiltration_excess = infiltration_excess
1030 self.subsurface_flow = subsurface_flow
1031 self.percolation = percolation
1032 self.tank_recharge = tank_recharge
1033 self.evaporation = evaporation
1034 self.precipitation = precipitation
1036 # Mass balance
1037 in_ = precipitation
1038 out_ = evaporation
1040 return (in_, out_)
1042 def route(self):
1043 """An outflow function that sends percolation, subsurface runoff and surface
1044 runoff to their respective tanks in the parent land node.
1046 Returns:
1047 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1048 for mass balance checking.
1049 """
1050 self.parent.surface_runoff.push_storage(self.infiltration_excess, force=True)
1051 self.parent.subsurface_runoff.push_storage(self.subsurface_flow, force=True)
1052 self.parent.percolation.push_storage(self.percolation, force=True)
1054 return (self.empty_vqip(), self.empty_vqip())
1056 def calculate_soil_temperature(self):
1057 """Process function that calculates soil temperature based on a weighted.
1059 average. This equation is from Lindstrom, Bishop & Lofvenius (2002),
1060 hydrological processes - but it is not clear what the parameters should be.
1062 Returns:
1063 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1064 for mass balance checking.
1065 """
1066 auto = self.storage["temperature"] * self.soil_temp_w_prev
1067 air = self.get_data_input("temperature") * self.soil_temp_w_air
1068 total_weight = (
1069 self.soil_temp_w_air + self.soil_temp_w_deep + self.soil_temp_w_prev
1070 )
1071 self.storage["temperature"] = (
1072 auto + air + self.soil_temp_deep * self.soil_temp_w_deep
1073 ) / total_weight
1075 return (self.empty_vqip(), self.empty_vqip())
1078class GrowingSurface(PerviousSurface):
1079 """"""
1081 def __init__(
1082 self,
1083 rooting_depth=1,
1084 ET_depletion_factor=1,
1085 crop_factor_stages=[1, 1],
1086 crop_factor_stage_dates=[0, 365],
1087 sowing_day=1,
1088 harvest_day=365,
1089 initial_soil_storage=None,
1090 **kwargs,
1091 ):
1092 """Extensive surface subclass that implements the CatchWat equations (Liu,
1093 Dobson & Mijic (2022) Science of the total environment), which in turn are
1094 primarily based on FAO document:
1095 https://www.fao.org/3/x0490e/x0490e0ehtm#soil%20water%20availability. This
1096 surface is a pervious surface that also has things that grow on it. This
1097 behaviour includes soil nutrient pools, crop planting/harvest calendars,
1098 erosion, crop behaviour.
1100 A key complexity of this surface is the nutrient pool (see wsimod/nodes/
1101 nutrient_pool.py), which is a class that tracks the amount of phosphorus and
1102 nitrogen in different states and performs transformations that occur in the
1103 phosphorus/nitrogen cycle. It is assumed that the phosphate/nitrate/nitrite/
1104 ammonia amounts in this Surface tank should track the dissolved inorganic pool
1105 in the nutrient pool. Meanwhile, the org-phosphorus/org-nitrogen amounts in this
1106 tank should track the dissolved organic pool in the nutrient pool. The total
1107 amount of pollutants that enter this tank may not be the same as the total
1108 amount that leave, because pollutants are transformed between inorganic/ organic
1109 and between wet/dry states - these transformations are accounted for in mass
1110 balance.
1112 For users to quickly enable/disable these nutrient processes, which are
1113 computationally intensive (in current case studies they account for about half
1114 of the total runtime), they are only active if 'nitrate' is one of the modelled
1115 pollutants. Note that the code will not check if nitrite/phosphate/
1116 org-phosphorus/org-nitrogen/ammonia are also included, but they should be if
1117 nitrate is included and otherwise the code will crash with a key error.
1119 Args:
1120 rooting_depth (float, optional): Depth of the soil tank (i.e., how deep do
1121 crop roots go). Defaults to 1.
1122 ET_depletion_factor (float, optional): Average fraction of soil that can be
1123 depleted from the root zone before moisture stress (reduction in ET)
1124 occurs. Defaults to 1.
1125 crop_factor_stages (list, optional): Crop factor is a multiplier on et0,
1126 more grown plants have higher transpiration and higher crop factors.
1127 This list shows changing crop factor at different times of year in
1128 relation to crop_factor_stage_dates. See wsimod/preprocessing/
1129 england_data_formatting.py/format_surfaces for further details on
1130 formulating these - since the interpolation used to find crop_factors in
1131 between the given values in the list is a bit involved. Defaults to
1132 [1,1].
1133 crop_factor_stage_dates (list, optional): Dates associated with
1134 crop_factor_stages. Defaults to [0, 365].
1135 sowing_day (int, optional): day of year that crops are sown. Defaults to 1.
1136 harvest_day (int, optional): day of year that crops are harvest. Defaults
1137 to 365.
1138 initial_soil_storage (dict or float, optional): Initial mass of solid
1139 pollutants
1140 in the soil nutrient pools (fast and adsorbed inorganic pools)
1142 Key assumptions:
1143 - In the soil water module, crop stages and crop coefficients control the
1144 evapotranspiration.
1145 - Fertiliser and manure application are the major source of soil nutrients,
1146 which are added
1147 into soil nutrient pools, including dissovled inorganic, dissolved
1148 organic, fast and humus for both nitrogen and phosphorus.
1149 - Nutrient transformation processes in soil are simulated, including fluxes
1150 between the soil
1151 nutrient pools, denitrification for nitrogen, adsorption/desorption for
1152 phosphorus. These processes are affected by temperature and soil
1153 moisture.
1154 - Crop uptake of nutrients are simulated based on crop stages, which is
1155 different for spring-sown
1156 and autumn-sown crops.
1157 - Soil erosion from the growing surface is simulated as one of the major
1158 sources of suspended solids
1159 in rivers, which is mainly affected by rainfall energy and crop/ground
1160 cover. Phosphorus will also be eroded along with the soil particles, in
1161 both adsorbed inorganic and humus form.
1163 Input data and parameter requirements:
1164 - `data_input_dict` can contain a variety of pollutant deposition data.
1165 `srp-fertiliser` describes phosphate. `noy-fertiliser` describes
1166 nitrogen as nitrates. `nhx-fertiliser` describes nitrogen as ammonia.
1167 `srp/noy/ nhx-manure` can also be used to specify manure application.
1168 _Units_: kg/m2/timestep (data is read at a monthly timestep)
1169 - Rooting depth.
1170 _Units_: m
1171 - Evapotranspiration depletion factor.
1172 _Units_: -
1173 - Sowing day, harvest day and crop calendars.
1174 _Units_: day number in Julian calendar
1175 - Crop factor.
1176 _Units_: -
1177 - Initial storage for solid pollutants.
1178 _Units_: kg
1180 """
1181 # Crop factors (set when creating object)
1182 self.ET_depletion_factor = (
1183 ET_depletion_factor # To do with water availability, p from FAOSTAT
1184 )
1185 self.rooting_depth = (
1186 rooting_depth # maximum depth that plants can absorb, Zr from FAOSTAT
1187 )
1188 depth = rooting_depth
1190 # Crop parameters
1191 self.crop_cover_max = 0.9 # [-] 0~1
1192 self.ground_cover_max = 0.3 # [-]
1193 # TODO... really I should just have this as an annual profile parameter and do
1194 # away with interpolation etc.
1195 self.crop_factor_stages = crop_factor_stages
1196 self.crop_factor_stage_dates = crop_factor_stage_dates
1197 self.sowing_day = sowing_day
1198 self.harvest_day = harvest_day
1200 # Soil moisture dependence parameters
1201 self.satact = 0.6 # [-] for calculating soil_moisture_dependence_factor
1202 self.thetaupp = 0.12 # [-] for calculating soil_moisture_dependence_factor
1203 self.thetalow = 0.08 # [-] for calculating soil_moisture_dependence_factor
1204 self.thetapow = 1 # [-] for calculating soil_moisture_dependence_factorself.
1205 # satact = 0.6 # [-] for calculating soil_moisture_dependence_factor
1207 # Crop uptake parameters
1208 self.uptake1 = (
1209 15 # [g/m2/y] shape factor for crop (Dissolved) Inorganic nitrogen uptake
1210 )
1211 self.uptake2 = (
1212 1 # [-] shape factor for crop (Dissolved) Inorganic nitrogen uptake
1213 )
1214 self.uptake3 = (
1215 0.02 # [1/day] shape factor for crop (Dissolved) Inorganic nitrogen uptake
1216 )
1217 self.uptake_PNratio = 1 / 7.2 # [-] P:N during crop uptake
1219 # Erosion parameters
1220 self.erodibility = 0.0025 # [g * d / (J * mm)]
1221 self.sreroexp = 1.2 # [-] surface runoff erosion exponent
1222 self.cohesion = 1 # [kPa]
1223 self.slope = 5 # [-] every 100
1224 self.srfilt = (
1225 0.7 # [-] ratio of eroded sediment left in surface runoff after filtration
1226 )
1227 self.macrofilt = 0.01 # [-] ratio of eroded sediment left in subsurface flow
1228 # after filtration
1230 # Denitrification parameters
1231 self.limpar = 0.7 # [-] above which denitrification begins
1232 self.exppar = 2.5 # [-] exponential parameter for
1233 # soil_moisture_dependence_factor_exp calculation
1234 self.hsatINs = 1 # [mg/l] for calculation of half-saturation concentration
1235 # dependence factor
1236 self.denpar = 0.015 # [-] denitrification rate coefficient
1238 # Adsorption parameters
1239 self.adosorption_nr_limit = 0.00001
1240 self.adsorption_nr_maxiter = 20
1241 self.kfr = 153.7 # [litter/kg] freundlich adsorption isoterm
1242 self.nfr = 1 / 2.6 # [-] freundlich exponential coefficient
1243 self.kadsdes = 0.03 # [1/day] adsorption/desorption coefficient
1245 # Other soil parameters
1246 self.bulk_density = 1300 # [kg/m3]
1247 super().__init__(depth=depth, **kwargs)
1249 (
1250 self.harvest_sow_calendar,
1251 self.ground_cover_stages,
1252 self.crop_cover_stages,
1253 self.autumn_sow,
1254 ) = self.infer_sow_harvest_calendar()
1256 # State variables
1257 self.days_after_sow = None
1258 self.crop_cover = 0
1259 self.ground_cover = 0
1260 self.crop_factor = 0
1261 self.et0_coefficient = 1
1263 (
1264 self.total_available_water,
1265 self.readily_available_water,
1266 ) = self.calculate_available_water()
1268 # Initiliase nutrient pools
1269 self.nutrient_pool = NutrientPool()
1271 self.inflows.insert(0, self.calc_crop_cover)
1272 if "nitrate" in constants.POLLUTANTS:
1273 # Populate function lists
1274 self.inflows.append(self.effective_precipitation_flushing)
1275 self.inflows.append(self.fertiliser)
1276 self.inflows.append(self.manure)
1277 # self.inflows.append(self.residue)
1279 self.processes.append(self.calc_temperature_dependence_factor)
1280 self.processes.append(self.calc_soil_moisture_dependence_factor)
1281 self.processes.append(self.soil_pool_transformation)
1282 self.processes.append(self.calc_crop_uptake)
1284 # TODO possibly move these into nutrient pool
1285 self.processes.append(self.erosion)
1286 self.processes.append(self.denitrification)
1287 self.processes.append(self.adsorption)
1289 # Reflect initial water concentration in dissolved nutrient pools
1290 self.nutrient_pool.dissolved_inorganic_pool.storage["P"] = self.storage[
1291 "phosphate"
1292 ]
1293 self.nutrient_pool.dissolved_inorganic_pool.storage["N"] = (
1294 self.storage["nitrate"]
1295 + self.storage["ammonia"]
1296 + self.storage["nitrite"]
1297 )
1298 self.nutrient_pool.dissolved_organic_pool.storage["P"] = self.storage[
1299 "org-phosphorus"
1300 ]
1301 self.nutrient_pool.dissolved_organic_pool.storage["N"] = self.storage[
1302 "org-nitrogen"
1303 ]
1304 if initial_soil_storage:
1305 self.initial_soil_storage = initial_soil_storage
1306 # Reflect initial nutrient stores in solid nutrient pools
1307 self.nutrient_pool.adsorbed_inorganic_pool.storage["P"] = (
1308 initial_soil_storage["phosphate"]
1309 )
1310 self.nutrient_pool.adsorbed_inorganic_pool.storage["N"] = (
1311 initial_soil_storage["ammonia"]
1312 + initial_soil_storage["nitrate"]
1313 + initial_soil_storage["nitrite"]
1314 )
1315 self.nutrient_pool.fast_pool.storage["N"] = initial_soil_storage[
1316 "org-nitrogen"
1317 ]
1318 self.nutrient_pool.fast_pool.storage["P"] = initial_soil_storage[
1319 "org-phosphorus"
1320 ]
1322 def infer_sow_harvest_calendar(self):
1323 """Infer basic sow/harvest calendar and indicate autumn-sown.
1324 Returns:
1325 (list): havest/sow calendar
1326 (list): ground cover stages
1327 (list): crop cover stages
1328 (boolean): indication for autumn-sown crops
1329 """
1330 # Infer basic sow/harvest calendar
1331 harvest_sow_calendar = [
1332 0,
1333 self.sowing_day,
1334 self.harvest_day,
1335 self.harvest_day + 1,
1336 365,
1337 ]
1338 ground_cover_stages = [0, 0, self.ground_cover_max, 0, 0]
1339 crop_cover_stages = [0, 0, self.crop_cover_max, 0, 0]
1341 # Use day number of 181 to indicate autumn-sown (from HYPE)
1342 if self.sowing_day > 181:
1343 autumn_sow = True
1344 else:
1345 autumn_sow = False
1347 return harvest_sow_calendar, ground_cover_stages, crop_cover_stages, autumn_sow
1349 def calculate_available_water(self):
1350 """Calculate total/readily available water based on capacity/wp.
1351 Returns:
1352 (float): total available water
1353 (float): readily available water
1354 """
1355 # Calculate parameters based on capacity/wp
1356 total_available_water = self.field_capacity_m - self.wilting_point_m
1357 if total_available_water < 0:
1358 print("warning: TAW < 0...")
1359 readily_available_water = total_available_water * self.ET_depletion_factor
1361 return total_available_water, readily_available_water
1363 def apply_overrides(self, overrides=Dict[str, Any]):
1364 """Override parameters.
1366 Enables a user to override any of the following parameters
1368 Args:
1369 overrides (Dict[str, Any]): Dict describing which parameters should
1370 be overridden (keys) and new values (values). Defaults to {}.
1371 """
1372 overwrite_params = [
1373 "ET_depletion_factor",
1374 "crop_cover_max",
1375 "ground_cover_max",
1376 "crop_factor_stages",
1377 "crop_factor_stage_dates",
1378 "sowing_day",
1379 "harvest_day",
1380 "satact",
1381 "thetaupp",
1382 "thetalow",
1383 "thetapow",
1384 "uptake1",
1385 "uptake2",
1386 "uptake3",
1387 "uptake_PNratio",
1388 "erodibility",
1389 "sreroexp",
1390 "cohesion",
1391 "slope",
1392 "srfilt",
1393 "macrofilt",
1394 "limpar",
1395 "exppar",
1396 "hsatINs",
1397 "denpar",
1398 "adosorption_nr_limit",
1399 "adsorption_nr_maxiter",
1400 "kfr",
1401 "nfr",
1402 "kadsdes",
1403 "bulk_density",
1404 ]
1405 for param in overwrite_params:
1406 setattr(self, param, overrides.pop(param, getattr(self, param)))
1408 if "depth" in overrides.keys():
1409 overrides.pop("depth")
1410 print(
1411 "ERROR: specifying depth is depreciated in overrides for \
1412 GrowingSurface, please specify rooting_depth instead"
1413 )
1414 self.rooting_depth = overrides.pop("rooting_depth", self.rooting_depth)
1415 overrides["depth"] = self.rooting_depth
1416 super().apply_overrides(overrides)
1418 (
1419 self.harvest_sow_calendar,
1420 self.ground_cover_stages,
1421 self.crop_cover_stages,
1422 self.autumn_sow,
1423 ) = self.infer_sow_harvest_calendar()
1424 (
1425 self.total_available_water,
1426 self.readily_available_water,
1427 ) = self.calculate_available_water()
1429 def pull_storage(self, vqip):
1430 """Pull water from the surface, updating the surface storage VQIP. Nutrient pool
1431 pollutants (nitrate/nitrite/ammonia/phosphate/org- phosphorus/ org-nitrogen) are
1432 removed in proportion to their amounts in the dissolved nutrient pools, if they
1433 are simulated. Other pollutants are removed in proportion to their amount in the
1434 surface tank.
1436 Args:
1437 vqip (dict): VQIP amount to be pulled, (only 'volume' key is needed)
1439 Returns:
1440 reply (dict): A VQIP amount successfully pulled from the tank
1441 """
1442 if self.storage["volume"] == 0:
1443 return self.empty_vqip()
1445 # Adjust based on available volume
1446 reply = min(vqip["volume"], self.storage["volume"])
1448 # Update reply to vqip (get concentration for non-nutrients)
1449 reply = self.v_change_vqip(self.storage, reply)
1451 if "nitrate" in constants.POLLUTANTS:
1452 # Update nutrient pool and get concentration for nutrients
1453 prop = reply["volume"] / self.storage["volume"]
1454 nutrients = self.nutrient_pool.extract_dissolved(prop)
1455 reply["nitrate"] = (
1456 nutrients["inorganic"]["N"]
1457 * self.storage["nitrate"]
1458 / (
1459 self.storage["nitrate"]
1460 + self.storage["ammonia"]
1461 + self.storage["nitrite"]
1462 )
1463 )
1464 reply["ammonia"] = (
1465 nutrients["inorganic"]["N"]
1466 * self.storage["ammonia"]
1467 / (
1468 self.storage["nitrate"]
1469 + self.storage["ammonia"]
1470 + self.storage["nitrite"]
1471 )
1472 )
1473 reply["nitrite"] = (
1474 nutrients["inorganic"]["N"]
1475 * self.storage["nitrite"]
1476 / (
1477 self.storage["nitrate"]
1478 + self.storage["ammonia"]
1479 + self.storage["nitrite"]
1480 )
1481 )
1482 reply["phosphate"] = nutrients["inorganic"]["P"]
1483 reply["org-phosphorus"] = nutrients["organic"]["P"]
1484 reply["org-nitrogen"] = nutrients["organic"]["N"]
1486 # Extract from storage
1487 self.storage = self.extract_vqip(self.storage, reply)
1489 return reply
1491 def quick_interp(self, x, xp, yp):
1492 """A simple version of np.interp to intepolate crop information on the fly.
1494 Args:
1495 x (int): Current time (i.e., day of year) xp (list): Predefined times (i.e.,
1496 list of days of year) yp (list): Predefined values associated with xp
1498 Returns:
1499 y (float): Interpolated value for current time
1500 """
1501 x_ind = bisect_left(xp, x)
1502 x_left = xp[x_ind - 1]
1503 x_right = xp[x_ind]
1504 dif = x - x_left
1505 y_left = yp[x_ind - 1]
1506 y_right = yp[x_ind]
1507 y = y_left + (y_right - y_left) * dif / (x_right - x_left)
1508 return y
1510 def calc_crop_cover(self):
1511 """Process function that calculates how much crop cover there is, assigns
1512 whether crops are sown/harvested, and calculates et0_coefficient based on growth
1513 stage of crops.
1515 Returns:
1516 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1517 for mass balance checking.
1518 """
1519 # Get current day of year
1520 doy = self.parent.t.dayofyear
1522 if self.parent.t.is_leap_year:
1523 # Hacky way to handle leap years
1524 if doy > 59:
1525 doy -= 1
1527 if self.days_after_sow is None:
1528 if self.parent.t.dayofyear == self.sowing_day:
1529 # sow
1530 self.days_after_sow = 0
1531 else:
1532 if self.parent.t.dayofyear == self.harvest_day:
1533 # harvest
1534 self.days_after_sow = None
1535 self.crop_factor = self.crop_factor_stages[0]
1536 self.crop_cover = 0
1537 self.ground_cover = 0
1538 else:
1539 # increment days since sow
1540 self.days_after_sow += 1
1542 # Calculate relevant parameters
1543 self.crop_factor = self.quick_interp(
1544 doy, self.crop_factor_stage_dates, self.crop_factor_stages
1545 )
1546 if self.days_after_sow:
1547 # Move outside of this if, if you want nonzero crop/ground cover outside of
1548 # season
1549 self.crop_cover = self.quick_interp(
1550 doy, self.harvest_sow_calendar, self.crop_cover_stages
1551 )
1552 self.ground_cover = self.quick_interp(
1553 doy, self.harvest_sow_calendar, self.ground_cover_stages
1554 )
1556 root_zone_depletion = max(self.field_capacity_m - self.get_smc(), 0)
1557 if root_zone_depletion < self.readily_available_water:
1558 crop_water_stress_coefficient = 1
1559 else:
1560 crop_water_stress_coefficient = max(
1561 0,
1562 (self.total_available_water - root_zone_depletion)
1563 / ((1 - self.ET_depletion_factor) * self.total_available_water),
1564 )
1566 self.et0_coefficient = crop_water_stress_coefficient * self.crop_factor
1568 return (self.empty_vqip(), self.empty_vqip())
1570 def adjust_vqip_to_liquid(self, vqip, deposition, in_):
1571 """Function to interoperate between surface tank and nutrient pool. Most
1572 depositions are given in terms of ammonia/nitrate/phosphate - they are then
1573 aggregated to total N or P to enter the nutrient pools. Depending on the source
1574 of deposition these may transform (e.g., some go to dissolved and some to
1575 solids) upon entering the nutrient pool. To reflect these transformations in the
1576 soil tank, the amounts entering the soil tank are adjusted proportionately.
1578 Args:
1579 vqip (dict): A VQIP amount of pollutants originally intended to enter the
1580 soil tank
1581 deposition (dict): A dict with nutrients (N and P) as keys, showing the
1582 total amount of nutrients entering the nutrient pool
1583 in_ (dict): A dict with nutrients as keys, showing the updated amount of
1584 nutrients that entered the nutrient pool as dissolved pollutants
1586 Returns:
1587 vqip (dict): A VQIP amount of pollutants that have been scaled to account
1588 for nutrient pool transformations
1589 """
1590 if "nitrate" in constants.POLLUTANTS:
1591 if deposition["N"] > 0:
1592 vqip["nitrate"] *= in_["N"] / deposition["N"]
1593 vqip["ammonia"] *= in_["N"] / deposition["N"]
1594 vqip["org-nitrogen"] *= in_["N"] / deposition["N"]
1595 if deposition["P"] > 0:
1596 vqip["phosphate"] *= in_["P"] / deposition["P"]
1597 vqip["org-phosphorus"] *= in_["P"] / deposition["P"]
1599 return vqip
1601 def effective_precipitation_flushing(self):
1602 """Remove the nutrients brought out by effective precipitation, which is surface
1603 runoff, subsurface runoff, and percolation, from the nutrients pool.
1605 Returns:
1606 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1607 for mass balance checking.
1608 """
1609 # inorganic
1610 out = self.nutrient_pool.get_empty_nutrient()
1611 out["N"] = (
1612 self.subsurface_flow["ammonia"]
1613 + self.subsurface_flow["nitrite"]
1614 + self.subsurface_flow["nitrate"]
1615 + self.percolation["ammonia"]
1616 + self.percolation["nitrite"]
1617 + self.percolation["nitrate"]
1618 + self.infiltration_excess["ammonia"]
1619 + self.infiltration_excess["nitrite"]
1620 + self.infiltration_excess["nitrate"]
1621 ) # TODO what happens if infiltration excess (the real part) has pollutants?
1622 out["P"] = (
1623 self.subsurface_flow["phosphate"]
1624 + self.percolation["phosphate"]
1625 + self.infiltration_excess["phosphate"]
1626 )
1627 self.nutrient_pool.dissolved_inorganic_pool.extract(out)
1629 # organic
1630 out = self.nutrient_pool.get_empty_nutrient()
1631 out["N"] = (
1632 self.subsurface_flow["org-nitrogen"]
1633 + self.percolation["org-nitrogen"]
1634 + self.infiltration_excess["org-nitrogen"]
1635 )
1636 out["P"] = (
1637 self.subsurface_flow["org-phosphorus"]
1638 + self.percolation["org-phosphorus"]
1639 + self.infiltration_excess["org-phosphorus"]
1640 )
1641 self.nutrient_pool.dissolved_organic_pool.extract(out)
1643 return (self.empty_vqip(), self.empty_vqip())
1645 def fertiliser(self):
1646 """Read, scale and allocate fertiliser, updating the tank.
1648 Returns:
1649 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1650 for mass balance checking.
1651 """
1652 # TODO tidy up fertiliser/manure/residue/deposition once preprocessing is sorted
1654 # Scale for surface
1655 nhx = self.get_data_input_surface("nhx-fertiliser") * self.area
1656 noy = self.get_data_input_surface("noy-fertiliser") * self.area
1657 srp = self.get_data_input_surface("srp-fertiliser") * self.area
1659 # Update as VQIP
1660 vqip = self.empty_vqip()
1661 vqip["ammonia"] = nhx
1662 vqip["nitrate"] = noy
1663 vqip["phosphate"] = srp
1665 # Enter nutrient pool
1666 deposition = self.nutrient_pool.get_empty_nutrient()
1667 deposition["N"] = vqip["nitrate"] + vqip["ammonia"]
1668 deposition["P"] = vqip["phosphate"]
1669 in_ = self.nutrient_pool.allocate_fertiliser(deposition)
1671 # Update tank
1672 vqip = self.adjust_vqip_to_liquid(vqip, deposition, in_)
1673 self.push_storage(vqip, force=True)
1675 return (vqip, self.empty_vqip())
1677 def manure(self):
1678 """Read, scale and allocate manure, updating the tank.
1680 Returns:
1681 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1682 for mass balance checking.
1683 """
1684 # Scale for surface
1685 nhx = self.get_data_input_surface("nhx-manure") * self.area
1686 noy = self.get_data_input_surface("noy-manure") * self.area
1687 srp = self.get_data_input_surface("srp-manure") * self.area
1689 # Formulate as VQIP
1690 vqip = self.empty_vqip()
1691 vqip["ammonia"] = nhx
1692 vqip["nitrate"] = noy
1693 vqip["phosphate"] = srp
1695 # Enter nutrient pool
1696 deposition = self.nutrient_pool.get_empty_nutrient()
1697 deposition["N"] = vqip["nitrate"] + vqip["ammonia"]
1698 deposition["P"] = vqip["phosphate"]
1699 in_ = self.nutrient_pool.allocate_manure(deposition)
1701 # Update tank
1702 vqip = self.adjust_vqip_to_liquid(vqip, deposition, in_)
1704 self.push_storage(vqip, force=True)
1706 return (vqip, self.empty_vqip())
1708 def residue(self):
1709 """Read, scale and allocate residue, updating the tank (NOT CURRENTLY USED
1710 BECAUSE NO DATA SOURCES FOR RESIDUE CAN BE IDENTIFIED).
1712 Returns:
1713 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1714 for mass balance checking.
1715 """
1716 nhx = self.get_data_input_surface("nhx-residue") * self.area
1717 noy = self.get_data_input_surface("noy-residue") * self.area
1718 srp = self.get_data_input_surface("srp-residue") * self.area
1720 vqip = self.empty_vqip()
1721 vqip["ammonia"] = nhx * self.nutrient_pool.fraction_residue_to_fast["N"]
1722 vqip["nitrate"] = noy * self.nutrient_pool.fraction_residue_to_fast["N"]
1723 vqip["org-nitrogen"] = (
1724 nhx + noy
1725 ) * self.nutrient_pool.fraction_residue_to_humus["N"]
1726 vqip["phosphate"] = srp * self.nutrient_pool.fraction_residue_to_fast["P"]
1727 vqip["org-phosphorus"] = srp * self.nutrient_pool.fraction_residue_to_humus["P"]
1729 deposition = self.nutrient_pool.get_empty_nutrient()
1730 deposition["N"] = vqip["nitrate"] + vqip["ammonia"] + vqip["org-nitrogen"]
1731 deposition["P"] = vqip["phosphate"] + vqip["org-phosphorus"]
1733 in_ = self.nutrient_pool.allocate_residue(deposition)
1734 vqip = self.adjust_vqip_to_liquid(vqip, deposition, in_)
1736 self.push_storage(vqip, force=True)
1738 return (vqip, self.empty_vqip())
1740 def soil_pool_transformation(self):
1741 """A process function that run transformation functions in the nutrient pool and
1742 updates the pollutant concentrations in the surface tank.
1744 Returns:
1745 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1746 for mass balance checking.
1747 """
1748 # Initialise mass balance tracking variables
1749 in_ = self.empty_vqip()
1750 out_ = self.empty_vqip()
1752 # Get proportion of nitrogen that is nitrate in the soil tank NOTE ignores
1753 # nitrite - couldn't find enough information on it
1754 nitrate_proportion = self.storage["nitrate"] / (
1755 self.storage["nitrate"] + self.storage["ammonia"]
1756 )
1758 # Run soil pool functions
1759 (
1760 increase_in_dissolved_inorganic,
1761 increase_in_dissolved_organic,
1762 ) = self.nutrient_pool.soil_pool_transformation()
1764 # Update tank and mass balance TODO .. there is definitely a neater way to write
1765 # this
1766 if increase_in_dissolved_inorganic["N"] > 0:
1767 # Increase in inorganic nitrogen, rescale back to nitrate and ammonia
1768 in_["nitrate"] = increase_in_dissolved_inorganic["N"] * nitrate_proportion
1769 in_["ammonia"] = increase_in_dissolved_inorganic["N"] * (
1770 1 - nitrate_proportion
1771 )
1772 else:
1773 # Decrease in inorganic nitrogen, rescale back to nitrate and ammonia
1774 out_["nitrate"] = -increase_in_dissolved_inorganic["N"] * nitrate_proportion
1775 out_["ammonia"] = -increase_in_dissolved_inorganic["N"] * (
1776 1 - nitrate_proportion
1777 )
1779 if increase_in_dissolved_organic["N"] > 0:
1780 # Increase in organic nitrogen
1781 in_["org-nitrogen"] = increase_in_dissolved_organic["N"]
1782 else:
1783 # Decrease in organic nitrogen
1784 out_["org-nitrogen"] = -increase_in_dissolved_organic["N"]
1786 if increase_in_dissolved_inorganic["P"] > 0:
1787 # Increase in inorganic phosphate
1788 in_["phosphate"] = increase_in_dissolved_inorganic["P"]
1789 else:
1790 # Decrease in inorganic phosphate
1791 out_["phosphate"] = -increase_in_dissolved_inorganic["P"]
1793 if increase_in_dissolved_organic["P"] > 0:
1794 # Increase in organic phosphorus
1795 in_["org-phosphorus"] = increase_in_dissolved_organic["P"]
1796 else:
1797 # Decrease in organic phosphorus
1798 out_["org-phosphorus"] = -increase_in_dissolved_organic["P"]
1800 # Update tank with inputs/outputs of pollutants
1801 _ = self.push_storage(in_, force=True)
1802 out2_ = self.pull_pollutants(out_)
1804 if not self.compare_vqip(out_, out2_):
1805 print("nutrient pool not tracking soil tank")
1807 return (in_, out_)
1809 def calc_temperature_dependence_factor(self):
1810 """Process function that calculates the temperature dependence factor for the
1811 nutrient pool (which impacts soil pool transformations).
1813 Returns:
1814 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1815 for mass balance checking.
1816 """
1817 # Parameters/equations from HYPE documentation
1818 if self.storage["temperature"] > 5:
1819 temperature_dependence_factor = 2 ** (
1820 (self.storage["temperature"] - 20) / 10
1821 )
1822 elif self.storage["temperature"] > 0:
1823 temperature_dependence_factor = self.storage["temperature"] / 5
1824 else:
1825 temperature_dependence_factor = 0
1826 self.nutrient_pool.temperature_dependence_factor = temperature_dependence_factor
1827 return (self.empty_vqip(), self.empty_vqip())
1829 def calc_soil_moisture_dependence_factor(self):
1830 """Process function that calculates the soil moisture dependence factor for the
1831 nutrient pool (which impacts soil pool transformations).
1833 Returns:
1834 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1835 for mass balance checking.
1836 """
1837 # Parameters/equations from HYPE documentation
1838 current_soil_moisture = self.get_smc()
1839 if current_soil_moisture >= self.field_capacity_m:
1840 self.nutrient_pool.soil_moisture_dependence_factor = self.satact
1841 elif current_soil_moisture <= self.wilting_point_m:
1842 self.nutrient_pool.soil_moisture_dependence_factor = 0
1843 else:
1844 fc_diff = self.field_capacity_m - current_soil_moisture
1845 fc_comp = (fc_diff / (self.thetaupp * self.rooting_depth)) ** self.thetapow
1846 fc_comp = (1 - self.satact) * fc_comp + self.satact
1847 wp_diff = current_soil_moisture - self.wilting_point_m
1848 wp_comp = (wp_diff / (self.thetalow * self.rooting_depth)) ** self.thetapow
1849 self.nutrient_pool.soil_moisture_dependence_factor = min(
1850 1, wp_comp, fc_comp
1851 )
1852 return (self.empty_vqip(), self.empty_vqip())
1854 def calc_crop_uptake(self):
1855 """Process function that calculates how much nutrient crops uptake and updates
1856 nutrient pool and surface tank.
1858 Returns:
1859 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1860 for mass balance checking.
1861 """
1862 # Parameters/equations from HYPE documentation
1864 # Initialise
1865 N_common_uptake = 0
1866 P_common_uptake = 0
1868 if self.days_after_sow:
1869 # If there are crops
1871 days_after_sow = self.days_after_sow
1873 if self.autumn_sow:
1874 temp_func = max(0, min(1, (self.storage["temperature"] - 5) / 20))
1875 days_after_sow -= 25 # Not sure why this is (but it's in HYPE)
1876 else:
1877 temp_func = 1
1879 # Calculate uptake
1880 uptake_par = (self.uptake1 - self.uptake2) * exp(
1881 -self.uptake3 * days_after_sow
1882 )
1883 if (uptake_par + self.uptake2) > 0:
1884 N_common_uptake = (
1885 self.uptake1
1886 * self.uptake2
1887 * self.uptake3
1888 * uptake_par
1889 / ((self.uptake2 + uptake_par) ** 2)
1890 )
1891 N_common_uptake *= temp_func * constants.G_M2_TO_KG_M2 * self.area # [kg]
1892 P_common_uptake = N_common_uptake * self.uptake_PNratio
1893 # calculate maximum available uptake
1894 N_maximum_available_uptake = (
1895 max(0, self.storage["volume"] - self.wilting_point_m * self.area)
1896 / self.storage["volume"]
1897 * self.nutrient_pool.dissolved_inorganic_pool.storage["N"]
1898 )
1899 P_maximum_available_uptake = (
1900 max(0, self.storage["volume"] - self.wilting_point_m * self.area)
1901 / self.storage["volume"]
1902 * self.nutrient_pool.dissolved_inorganic_pool.storage["P"]
1903 )
1905 uptake = {
1906 "P": min(P_common_uptake, P_maximum_available_uptake),
1907 "N": min(N_common_uptake, N_maximum_available_uptake),
1908 }
1909 crop_uptake = self.nutrient_pool.dissolved_inorganic_pool.extract(uptake)
1910 out_ = self.empty_vqip()
1912 # Assuming plants eat N and P as nitrate and phosphate
1913 out_["nitrate"] = crop_uptake["N"]
1914 out_["phosphate"] = crop_uptake["P"]
1916 out2_ = self.pull_pollutants(out_)
1917 if not self.compare_vqip(out_, out2_):
1918 print("nutrient pool not tracking soil tank")
1920 return (self.empty_vqip(), out_)
1921 else:
1922 return (self.empty_vqip(), self.empty_vqip())
1924 def erosion(self):
1925 """Outflow function that erodes adsorbed/humus phosphorus and sediment and sends
1926 onwards to percolation/surface runoff/subsurface runoff.
1928 Returns:
1929 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1930 for mass balance checking.
1931 """
1932 # Parameters/equations from HYPE documentation (which explains why my
1933 # documentation is a bit ambiguous - because theirs is too)
1935 # Convert precipitation to MM since all the equations assume that
1936 precipitation_depth = self.get_data_input("precipitation") * constants.M_TO_MM
1938 # Calculate how much rain is mobilising erosion
1939 if precipitation_depth > 5:
1940 rainfall_energy = 8.95 + 8.44 * log10(
1941 precipitation_depth
1942 * (
1943 0.257
1944 + sin(2 * constants.PI * ((self.parent.t.dayofyear - 70) / 365))
1945 * 0.09
1946 )
1947 * 2
1948 )
1949 rainfall_energy *= precipitation_depth
1950 mobilised_rain = rainfall_energy * (1 - self.crop_cover) * self.erodibility
1951 else:
1952 mobilised_rain = 0
1954 # Calculate if any infiltration is mobilising erosion
1955 if self.infiltration_excess["volume"] > 0:
1956 mobilised_flow = (
1957 self.infiltration_excess["volume"] / self.area * constants.M_TO_MM * 365
1958 ) ** self.sreroexp
1959 mobilised_flow *= (
1960 (1 - self.ground_cover)
1961 * (1 / (0.5 * self.cohesion))
1962 * sin(self.slope / 100)
1963 / 365
1964 )
1965 else:
1966 mobilised_flow = 0
1968 # Sum flows (not sure why surface runoff isn't included) TODO I'm pretty sure it
1969 # should be included here
1970 total_flows = (
1971 self.infiltration_excess["volume"]
1972 + self.subsurface_flow["volume"]
1973 + self.percolation["volume"]
1974 ) # m3/dt + self.tank_recharge['volume'] (guess not needed)
1976 # Convert to MM/M2
1977 erodingflow = total_flows / self.area * constants.M_TO_MM
1979 # Calculate eroded sediment
1980 transportfactor = min(1, (erodingflow / 4) ** 1.3)
1981 erodedsed = (
1982 1000 * (mobilised_flow + mobilised_rain) * transportfactor
1983 ) # [kg/km2]
1984 # TODO not sure what conversion this HYPE 1000 is referring to
1986 # soil erosion with adsorbed inorganic phosphorus and humus phosphorus (erodedP
1987 # as P in eroded sediments and effect of enrichment)
1988 if erodingflow > 4:
1989 enrichment = 1.5
1990 elif erodingflow > 0:
1991 enrichment = 4 - (4 - 1.5) * erodingflow / 4
1992 else:
1993 return (self.empty_vqip(), self.empty_vqip())
1995 # Get erodable phosphorus
1996 erodableP = (
1997 self.nutrient_pool.get_erodable_P() / self.area * constants.KG_M2_TO_KG_KM2
1998 )
1999 erodedP = (
2000 erodedsed
2001 * (
2002 erodableP
2003 / (
2004 self.rooting_depth
2005 * constants.M_TO_KM
2006 * self.bulk_density
2007 * constants.KG_M3_TO_KG_KM3
2008 )
2009 )
2010 * enrichment
2011 ) # [kg/km2]
2013 # Convert to kg
2014 erodedP *= self.area * constants.M2_TO_KM2 # [kg]
2015 erodedsed *= self.area * constants.M2_TO_KM2 # [kg]
2017 # Allocate to different flows
2018 surface_erodedP = (
2019 self.srfilt * self.infiltration_excess["volume"] / total_flows * erodedP
2020 ) # [kg]
2021 surface_erodedsed = (
2022 self.srfilt * self.infiltration_excess["volume"] / total_flows * erodedsed
2023 ) # [kg]
2025 subsurface_erodedP = (
2026 self.macrofilt * self.subsurface_flow["volume"] / total_flows * erodedP
2027 ) # [kg]
2028 subsurface_erodedsed = (
2029 self.macrofilt * self.subsurface_flow["volume"] / total_flows * erodedsed
2030 ) # [kg]
2032 percolation_erodedP = (
2033 self.macrofilt * self.percolation["volume"] / total_flows * erodedP
2034 ) # [kg]
2035 percolation_erodedsed = (
2036 self.macrofilt * self.percolation["volume"] / total_flows * erodedsed
2037 ) # [kg]
2039 # Track mass balance
2040 in_ = self.empty_vqip()
2042 # Total eroded phosphorus
2043 eff_erodedP = percolation_erodedP + surface_erodedP + subsurface_erodedP # [kg]
2044 if eff_erodedP > 0:
2045 # Update nutrient pool
2046 org_removed, inorg_removed = self.nutrient_pool.erode_P(eff_erodedP)
2047 total_removed = inorg_removed + org_removed
2049 if abs(total_removed - eff_erodedP) > constants.FLOAT_ACCURACY:
2050 print("weird nutrients")
2052 # scale flows to split between inorganic and organic eroded P
2053 self.infiltration_excess["org-phosphorus"] += (
2054 surface_erodedP * org_removed / eff_erodedP
2055 )
2056 self.subsurface_flow["org-phosphorus"] += (
2057 subsurface_erodedP * org_removed / eff_erodedP
2058 )
2059 self.percolation["org-phosphorus"] += (
2060 percolation_erodedP * org_removed / eff_erodedP
2061 )
2063 # TODO Leon reckons this is conceptually dodgy.. but i'm not sure where else
2064 # adsorbed inorganic phosphorus should go
2065 self.infiltration_excess["phosphate"] += (
2066 surface_erodedP * inorg_removed / eff_erodedP
2067 )
2068 self.subsurface_flow["phosphate"] += (
2069 subsurface_erodedP * inorg_removed / eff_erodedP
2070 )
2071 self.percolation["phosphate"] += (
2072 percolation_erodedP * inorg_removed / eff_erodedP
2073 )
2075 # Entering the model (no need to uptake surface tank because both adsorbed
2076 # inorganic pool and humus pool are solids and so no tracked in the soil
2077 # water tank)
2078 in_["phosphate"] = inorg_removed
2079 in_["org-phosphorus"] = org_removed
2080 else:
2081 pass
2083 # Track sediment as solids
2084 self.infiltration_excess["solids"] += surface_erodedsed
2085 self.subsurface_flow["solids"] += subsurface_erodedsed
2086 self.percolation["solids"] += percolation_erodedsed
2088 in_["solids"] = surface_erodedsed + subsurface_erodedsed + percolation_erodedsed
2090 return (in_, self.empty_vqip())
2092 def denitrification(self):
2093 """Outflow function that performs denitirication processes, updating nutrient
2094 pool and soil tank.
2096 Returns:
2097 (tuple): A tuple containing a VQIP amount for model inputs and outputs
2098 for mass balance checking.
2099 """
2100 # Parameters/equations from HYPE documentation TODO could more of this be moved
2101 # to NutrientPool Calculate soil moisture dependence of denitrification
2102 soil_moisture_content = self.get_smc()
2103 if soil_moisture_content > self.field_capacity_m:
2104 denitrifying_soil_moisture_dependence = 1
2105 elif soil_moisture_content / self.field_capacity_m > self.limpar:
2106 denitrifying_soil_moisture_dependence = (
2107 ((soil_moisture_content / self.field_capacity_m) - self.limpar)
2108 / (1 - self.limpar)
2109 ) ** self.exppar
2110 else:
2111 denitrifying_soil_moisture_dependence = 0
2112 return (self.empty_vqip(), self.empty_vqip())
2114 # Get dissolved inorg nitrogen as a concentration and calculate factor
2115 din_conc = (
2116 self.nutrient_pool.dissolved_inorganic_pool.storage["N"]
2117 / self.storage["volume"]
2118 ) # [kg/m3]
2119 din_conc *= constants.KG_M3_TO_MG_L
2120 half_saturation_concentration_dependence_factor = din_conc / (
2121 din_conc + self.hsatINs
2122 )
2124 # Calculate and extract dentrified nitrogen
2125 denitrified_N = (
2126 self.nutrient_pool.dissolved_inorganic_pool.storage["N"]
2127 * half_saturation_concentration_dependence_factor
2128 * denitrifying_soil_moisture_dependence
2129 * self.nutrient_pool.temperature_dependence_factor
2130 * self.denpar
2131 )
2132 denitrified_request = self.nutrient_pool.get_empty_nutrient()
2133 denitrified_request["N"] = denitrified_N
2134 denitrified_N = self.nutrient_pool.dissolved_inorganic_pool.extract(
2135 denitrified_request
2136 )
2138 # Leon reckons this should leave the model (though I think technically some
2139 # small amount goes to nitrite)
2140 out_ = self.empty_vqip()
2141 out_["nitrate"] = denitrified_N["N"]
2143 # Update tank
2144 out2_ = self.pull_pollutants(out_)
2145 if not self.compare_vqip(out_, out2_):
2146 print("nutrient pool not tracking soil tank")
2148 return (self.empty_vqip(), out_)
2150 def adsorption(self):
2151 """Outflow function that calculates phosphorus adsorption/desorptions and
2152 updates soil tank and nutrient pools.
2154 Returns:
2155 (tuple): A tuple containing a VQIP amount for model inputs and outputs
2156 for mass balance checking.
2157 """
2158 # Parameters/equations from HYPE documentation TODO could this be moved to the
2159 # nutrient pool?
2161 # Initialise mass balance checking
2162 in_ = self.empty_vqip()
2163 out_ = self.empty_vqip()
2165 # Get total phosphorus in pool available for adsorption/desorption
2166 limit = self.adosorption_nr_limit
2167 ad_de_P_pool = (
2168 self.nutrient_pool.adsorbed_inorganic_pool.storage["P"]
2169 + self.nutrient_pool.dissolved_inorganic_pool.storage["P"]
2170 ) # [kg]
2171 ad_de_P_pool /= self.area * constants.M2_TO_KM2 # [kg/km2]
2172 if ad_de_P_pool == 0:
2173 return (self.empty_vqip(), self.empty_vqip())
2175 # Calculate coefficient and concentration of adsorbed phosphorus
2176 soil_moisture_content = (
2177 self.get_smc() * constants.M_TO_MM
2178 ) # [mm] (not sure why HYPE has this in mm but whatever)
2179 conc_sol = (
2180 self.nutrient_pool.adsorbed_inorganic_pool.storage["P"]
2181 * constants.KG_TO_MG
2182 / (self.bulk_density * self.rooting_depth * self.area)
2183 ) # [mg P/kg soil]
2184 coeff = self.kfr * self.bulk_density * self.rooting_depth # [mm]
2186 # calculate equilibrium concentration
2187 if conc_sol <= 0:
2188 # Not sure how this would happen
2189 print("Warning: soil partP <=0. Freundlich will give error, take shortcut.")
2190 xn_1 = ad_de_P_pool / (soil_moisture_content + coeff) # [mg/l]
2191 ad_P_equi_conc = self.kfr * xn_1 # [mg/ kg]
2192 else:
2193 # Newton-Raphson method
2194 x0 = exp(
2195 (log(conc_sol) - log(self.kfr)) / self.nfr
2196 ) # initial guess of equilibrium liquid concentration
2197 fxn = x0 * soil_moisture_content + coeff * (x0**self.nfr) - ad_de_P_pool
2198 xn = x0
2199 xn_1 = xn
2200 j = 0
2201 while (
2202 abs(fxn) > limit and j < self.adsorption_nr_maxiter
2203 ): # iteration to calculate equilibrium concentations
2204 fxn = xn * soil_moisture_content + coeff * (xn**self.nfr) - ad_de_P_pool
2205 fprimxn = soil_moisture_content + self.nfr * coeff * (
2206 xn ** (self.nfr - 1)
2207 )
2208 dx = fxn / fprimxn
2209 if abs(dx) < (0.000001 * xn):
2210 # From HYPE... not sure what it means
2211 break
2212 xn_1 = xn - dx
2213 if xn_1 <= 0:
2214 xn_1 = 1e-10
2215 xn = xn_1
2216 j += 1
2217 ad_P_equi_conc = self.kfr * (xn_1**self.nfr)
2218 # print(ad_P_equi_conc, conc_sol)
2220 # Calculate new pool and concentration, depends on the equilibrium concentration
2221 if abs(ad_P_equi_conc - conc_sol) > 1e-6:
2222 request = self.nutrient_pool.get_empty_nutrient()
2224 # TODO not sure about this if statement, surely it would be triggered every
2225 # time
2226 adsdes = (ad_P_equi_conc - conc_sol) * (
2227 1 - exp(-self.kadsdes)
2228 ) # kinetic adsorption/desorption
2229 request["P"] = (
2230 adsdes
2231 * self.bulk_density
2232 * self.rooting_depth
2233 * (self.area * constants.M2_TO_KM2)
2234 ) # [kg]
2235 if request["P"] > 0:
2236 # Adsorption
2237 adsorbed = self.nutrient_pool.dissolved_inorganic_pool.extract(request)
2238 if (adsorbed["P"] - request["P"]) > constants.FLOAT_ACCURACY:
2239 print("Warning: freundlich flow adjusted, was larger than pool")
2240 self.nutrient_pool.adsorbed_inorganic_pool.receive(adsorbed)
2242 # Dissolved leaving the soil water tank and becoming solid
2243 out_["phosphate"] = adsorbed["P"]
2245 # Update tank
2246 out2_ = self.pull_pollutants(out_)
2247 if not self.compare_vqip(out_, out2_):
2248 print("nutrient pool not tracking soil tank")
2249 else:
2250 # Desorption
2251 request["P"] = -request["P"]
2252 desorbed = self.nutrient_pool.adsorbed_inorganic_pool.extract(request)
2253 if (desorbed["P"] - request["P"]) > constants.FLOAT_ACCURACY:
2254 print("Warning: freundlich flow adjusted, was larger than pool")
2255 self.nutrient_pool.dissolved_inorganic_pool.receive(desorbed)
2257 # Solid phosphorus becoming inorganic P in the soil water tank
2258 in_["phosphate"] = desorbed["P"]
2259 _ = self.push_storage(in_, force=True)
2261 return (in_, out_)
2263 def dry_deposition_to_tank(self, vqip):
2264 """Allocate dry deposition to surface tank, updating nutrient pool accordingly.
2266 Args:
2267 vqip (dict): A VQIP amount of dry deposition to send to tank
2269 Returns:
2270 vqip (dict): A VQIP amount of dry deposition that entered the tank (used
2271 for mass balance checking)
2272 """
2273 # Convert to nutrients
2274 deposition = self.nutrient_pool.get_empty_nutrient()
2275 deposition["N"] = vqip["nitrate"] + vqip["ammonia"]
2276 deposition["P"] = vqip["phosphate"]
2278 # Update nutrient pool
2279 in_ = self.nutrient_pool.allocate_dry_deposition(deposition)
2280 vqip = self.adjust_vqip_to_liquid(vqip, deposition, in_)
2282 # Update tank
2283 self.push_storage(vqip, force=True)
2284 return vqip
2286 def wet_deposition_to_tank(self, vqip):
2287 """Allocate wet deposition to surface tank, updating nutrient pool accordingly.
2289 Args:
2290 vqip (dict): A VQIP amount of dry deposition to send to tank
2292 Returns:
2293 vqip (dict): A VQIP amount of dry deposition that entered the tank (used
2294 for mass balance checking)
2295 """
2296 # Convert to nutrients
2297 deposition = self.nutrient_pool.get_empty_nutrient()
2298 deposition["N"] = vqip["nitrate"] + vqip["ammonia"]
2299 deposition["P"] = vqip["phosphate"]
2301 # Update nutrient pool
2302 in_ = self.nutrient_pool.allocate_wet_deposition(deposition)
2303 vqip = self.adjust_vqip_to_liquid(vqip, deposition, in_)
2305 # Update tank
2306 self.push_storage(vqip, force=True)
2307 return vqip
2310class IrrigationSurface(GrowingSurface):
2311 """"""
2313 def __init__(self, irrigation_coefficient=0.1, **kwargs):
2314 """A subclass of GrowingSurface that can calculate water demand for the crops
2315 that is not met by precipitation and use the parent node to acquire water. When
2316 the surface is created by the parent node, the irrigate function below is
2317 assigned.
2319 Args:
2320 irrigation_coefficient (float, optional): proportion area irrigated *
2321 proportion of demand met. Defaults to 0.1.
2322 """
2323 # Assign param
2324 self.irrigation_coefficient = irrigation_coefficient # proportion area
2325 # irrigated * proportion of demand met
2327 super().__init__(**kwargs)
2329 def apply_overrides(self, overrides=Dict[str, Any]):
2330 """Override parameters.
2332 Enables a user to override irrigation_coefficient
2334 Args:
2335 overrides (Dict[str, Any]): Dict describing which parameters should
2336 be overridden (keys) and new values (values). Defaults to {}.
2337 """
2338 self.irrigation_coefficient = overrides.pop(
2339 "irrigation_coefficient", self.irrigation_coefficient
2340 )
2341 super().apply_overrides(overrides)
2343 def irrigate(self):
2344 """Calculate water demand for crops and call parent node to acquire water,
2345 updating surface tank and nutrient pools."""
2346 if self.days_after_sow:
2347 # Irrigation is just difference between evaporation and precipitation amount
2348 irrigation_demand = (
2349 max(self.evaporation["volume"] - self.precipitation["volume"], 0)
2350 * self.irrigation_coefficient
2351 )
2352 if irrigation_demand > constants.FLOAT_ACCURACY:
2353 root_zone_depletion = self.get_cmd()
2354 if root_zone_depletion <= constants.FLOAT_ACCURACY:
2355 # TODO this isn't in FAO... but seems sensible
2356 irrigation_demand = 0
2358 # Pull water using parent node
2359 supplied = self.parent.pull_distributed(
2360 {"volume": irrigation_demand},
2361 of_type=["River", "Node", "Groundwater", "Reservoir"],
2362 )
2364 # update tank
2365 _ = self.push_storage(supplied, force=True)
2367 # update nutrient pools
2368 organic = {
2369 "N": supplied["org-nitrogen"],
2370 "P": supplied["org-phosphorus"],
2371 }
2372 inorganic = {
2373 "N": supplied["ammonia"] + supplied["nitrate"],
2374 "P": supplied["phosphate"],
2375 }
2376 self.nutrient_pool.allocate_organic_irrigation(organic)
2377 self.nutrient_pool.allocate_inorganic_irrigation(inorganic)
2380class GardenSurface(GrowingSurface):
2381 """"""
2383 # TODO - probably a simplier version of this is useful, building just on
2384 # pervioussurface
2385 def __init__(self, **kwargs):
2386 """A specific surface for gardens that treats the garden as a grass crop, but
2387 that can calculate/receive irrigation through functions that are assigned by the
2388 parent land node's handlers, which in turn are expected to be triggered by a
2389 query from an attached Demand node."""
2390 super().__init__(**kwargs)
2392 def calculate_irrigation_demand(self, ignore_vqip=None):
2393 """A check function (assigned by parent to push check from demand nodes) that
2394 calculations irrigation demand (i.e., difference between evaporation and
2395 preciptiation).
2397 Args:
2398 ignore_vqip (any, optional): Conventional push checks send an optional
2399 VQIP amount, however the intention with this check is to get the
2400 irrigation demand
2402 Returns:
2403 reply (dict): A VQIP amount of irrigation demand (note only 'volume' key
2404 is used)
2405 """
2406 # Calculate irrigation demand
2407 irrigation_demand = max(
2408 self.evaporation["volume"] - self.precipitation["volume"], 0
2409 )
2411 root_zone_depletion = self.get_cmd()
2412 if root_zone_depletion <= constants.FLOAT_ACCURACY:
2413 # TODO this isn't in FAO... but seems sensible
2414 irrigation_demand = 0
2416 # Reply as VQIP
2417 reply = self.empty_vqip()
2418 reply["volume"] = irrigation_demand
2419 return reply
2421 def receive_irrigation_demand(self, vqip):
2422 """A set function (assigned by parent to push set from demand nodes) that
2423 assigns irrigation water supply to the surface tank.
2425 Args:
2426 vqip (dict): A VQIP amount of irrigation to receive
2428 Returns:
2429 (dict): A VQIP amount of irrigation that was not received (should always
2430 be empty)
2431 """
2432 # update tank
2433 return self.push_storage(vqip, force=True)
2436class VariableAreaSurface(GrowingSurface):
2437 """"""
2439 def __init__(self, current_surface_area=0, **kwargs):
2440 """"""
2441 super().__init__(**kwargs)
2442 self.get_climate = self.get_climate_
2443 self.current_surface_area = current_surface_area
2445 def get_climate_(self):
2446 """
2448 Returns:
2450 """
2451 precipitation_depth = self.get_data_input("precipitation")
2452 evaporation_depth = self.get_data_input("et0") * self.et0_coefficient
2454 precipitation_depth *= self.current_surface_area / self.area
2455 evaporation_depth *= self.current_surface_area / self.area
2457 return precipitation_depth, evaporation_depth