Coverage for wsimod/nodes/land.py: 9%
713 statements
« prev ^ index » next coverage.py v7.3.2, created at 2024-01-11 16:39 +0000
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
10from wsimod.core import constants
11from wsimod.nodes.nodes import DecayTank, Node, ResidenceTank
12from wsimod.nodes.nutrient_pool import NutrientPool
15class Land(Node):
16 """"""
18 def __init__(
19 self,
20 name,
21 subsurface_residence_time=5,
22 percolation_residence_time=50,
23 surface_residence_time=1,
24 surfaces=[],
25 data_input_dict={},
26 ):
27 """An extensive node class that represents land processes (agriculture, soil,
28 subsurface flow, rural runoff, urban drainage, pollution deposition). The
29 expected use is that each distinctive type of land cover (different crop types,
30 gardens, forests, impervious urban drainage, etc.) each have a Surface object,
31 which is a subclass of Tank. The land node will iterate over its surfaces each
32 timestep, which will generally (except in the case of an impervious surface)
33 send water to three common Tanks: surface flow, subsurface flow and percolation.
34 These tanks will then send flows to rivers or groundwater.
36 (See wsimod/nodes/land.py/Surface and subclasses for currently available
37 surfaces)
39 Args:
40 name (str): node name. subsurface_residence_time (float, optional):
41 Residence time for
42 subsurface flow (see nodes.py/ResidenceTank). Defaults to 5.
43 percolation_residence_time (int, optional): Residence time for
44 percolation flow (see nodes.py/ResidenceTank). Defaults to 50.
45 surface_residence_time (int, optional): Residence time for surface flow
46 (see nodes.py/ResidenceTank). Defaults to 1.
47 surfaces (list, optional): list of dicts where each dict describes the
48 parameters of each surface in the Land node. Each dict also contains an
49 entry under 'type_' which describes which subclass of surface to use.
50 Defaults to [].
51 data_input_dict (dict, optional): Dictionary of data inputs relevant for
52 the node (generally, et0, precipitation and temperature). Keys are
53 tuples where first value is the name of the variable to read from the
54 dict and the second value is the time. Defaults to {}.
56 Functions intended to call in orchestration:
57 run apply_irrigation (if used)
59 Key assumptions:
60 - Percolation, surface runoff, and subsurface runoff, can be described with
61 a residence-time method.
62 - Flows to percolation, surface runoff, and subsurface runoff are
63 generated by different hydrological response units (subclasses of
64 `land.py/Surface`), but aggregated for a given land node.
65 - Flows to percolation are distributed to `storage.py/Groundwater`
66 nodes while surface/subsurface runoff to `nodes.py/Node` or
67 `storage.py/River` nodes.
68 - Input data associated with the land node (precipitation,
69 temperature, evapotranspiartion) are the same for every surface.
70 - Water received from `sewer.py/Sewer` objects is sent to the first
71 `land.py/ImperviousSurface` in the surfaces list.
73 Input data and parameter requirements:
74 - Precipitation and evapotranspiration are in the `data_input_dict`
75 at the model timestep. _Units_: metres/timestep
76 - Temperature in the `data_input_dict` at the model timestep.
77 _Units_: C
78 - Residence time of surface, subsurface and percolation flows.
79 _Units_: number of timesteps
80 """
81 # Assign parameters
82 self.subsurface_residence_time = subsurface_residence_time
83 self.percolation_residence_time = percolation_residence_time
84 self.surface_residence_time = surface_residence_time
85 self.data_input_dict = data_input_dict
87 super().__init__(name, data_input_dict=data_input_dict)
89 # This could be a deny but then you would have to know in advance whether a
90 # demand node has any gardening or not
91 self.push_check_handler[("Demand", "Garden")] = lambda x: self.empty_vqip()
92 self.push_set_handler[("Demand", "Garden")] = lambda x: self.empty_vqip()
94 # Create surfaces
95 self.irrigation_functions = [lambda: None]
97 surfaces_ = surfaces.copy()
98 surfaces = []
99 for surface in surfaces_:
100 # Assign parent (for data reading and to determine where to send flows to)
101 surface["parent"] = self
103 # Get surface type
104 type_ = surface["type_"]
105 del surface["type_"]
107 # Instantiate surface and store in list of surfaces
108 surfaces.append(getattr(sys.modules[__name__], type_)(**surface))
110 # Assign ds (mass balance checking)
111 self.mass_balance_ds.append(surfaces[-1].ds)
113 # Assign any irrigation functions
114 if isinstance(surfaces[-1], IrrigationSurface):
115 # TODO, this should probably be done in the Surface initialisation
116 self.irrigation_functions.append(surfaces[-1].irrigate)
118 # Assign garden surface functions
119 if isinstance(surfaces[-1], GardenSurface):
120 # TODO, this should probably be done in the Surface initialisation
121 self.push_check_handler[("Demand", "Garden")] = surfaces[
122 -1
123 ].calculate_irrigation_demand
124 self.push_set_handler[("Demand", "Garden")] = surfaces[
125 -1
126 ].receive_irrigation_demand
128 # Update handlers
129 self.push_set_handler["default"] = self.push_set_deny
130 self.push_check_handler["default"] = self.push_check_deny
131 self.push_set_handler["Sewer"] = self.push_set_sewer
133 # Create subsurface runoff, surface runoff and percolation tanks Can also do as
134 # timearea if this seems dodge (that is how it is done in IHACRES) TODO should
135 # these be decayresidencetanks?
136 self.subsurface_runoff = ResidenceTank(
137 residence_time=self.subsurface_residence_time,
138 capacity=constants.UNBOUNDED_CAPACITY,
139 )
140 self.percolation = ResidenceTank(
141 residence_time=self.percolation_residence_time,
142 capacity=constants.UNBOUNDED_CAPACITY,
143 )
144 self.surface_runoff = ResidenceTank(
145 residence_time=self.surface_residence_time,
146 capacity=constants.UNBOUNDED_CAPACITY,
147 )
149 # Store surfaces
150 self.surfaces = surfaces
152 # Mass balance checkign vqips and functions
153 self.running_inflow_mb = self.empty_vqip()
154 self.running_outflow_mb = self.empty_vqip()
156 self.mass_balance_in.append(lambda: self.running_inflow_mb)
157 self.mass_balance_out.append(lambda: self.running_outflow_mb)
158 self.mass_balance_ds.append(self.surface_runoff.ds)
159 self.mass_balance_ds.append(self.subsurface_runoff.ds)
160 self.mass_balance_ds.append(self.percolation.ds)
162 def apply_irrigation(self):
163 """Iterate over any irrigation functions (needs further testing..
165 maybe).
166 """
167 for f in self.irrigation_functions:
168 f()
170 def run(self):
171 """Call the run function in all surfaces, update surface/subsurface/ percolation
172 tanks, discharge to rivers/groundwater."""
173 # Run all surfaces
174 for surface in self.surfaces:
175 surface.run()
177 # Apply residence time to percolation
178 percolation = self.percolation.pull_outflow()
180 # Distribute percolation
181 reply = self.push_distributed(percolation, of_type=["Groundwater"])
183 if reply["volume"] > constants.FLOAT_ACCURACY:
184 # Update percolation 'tank'
185 _ = self.percolation.push_storage(reply, force=True)
187 # Apply residence time to subsurface/surface runoff
188 surface_runoff = self.surface_runoff.pull_outflow()
189 subsurface_runoff = self.subsurface_runoff.pull_outflow()
191 # Total runoff
192 total_runoff = self.sum_vqip(surface_runoff, subsurface_runoff)
193 if total_runoff["volume"] > 0:
194 # Send to rivers (or nodes, which are assumed to be junctions)
195 reply = self.push_distributed(total_runoff, of_type=["River", "Node"])
197 # Redistribute total_runoff not sent
198 if reply["volume"] > 0:
199 reply_surface = self.v_change_vqip(
200 reply,
201 reply["volume"] * surface_runoff["volume"] / total_runoff["volume"],
202 )
203 reply_subsurface = self.v_change_vqip(
204 reply,
205 reply["volume"]
206 * subsurface_runoff["volume"]
207 / total_runoff["volume"],
208 )
210 # Update surface/subsurface runoff 'tanks'
211 if reply_surface["volume"] > 0:
212 self.surface_runoff.push_storage(reply_surface, force=True)
213 if reply_subsurface["volume"] > 0:
214 self.subsurface_runoff.push_storage(reply_subsurface, force=True)
216 def push_set_sewer(self, vqip):
217 """Receive water from a sewer and send it to the first ImperviousSurface in
218 surfaces.
220 Args:
221 vqip (dict): A VQIP amount to be sent to the impervious surface
223 Returns:
224 vqip (dict): A VQIP amount of water that was not received
225 """
226 # TODO currently just push to the first impervious surface... not sure if people
227 # will be having multiple impervious surfaces. If people would be only having
228 # one then it would make sense to store as a parameter... this is probably fine
229 # for now
230 for surface in self.surfaces:
231 if isinstance(surface, ImperviousSurface):
232 vqip = self.surface.push_storage(vqip, force=True)
233 break
234 return vqip
236 def end_timestep(self):
237 """Update mass balance and end timestep of all tanks (and surfaces)."""
238 self.running_inflow_mb = self.empty_vqip()
239 self.running_outflow_mb = self.empty_vqip()
240 for tanks in self.surfaces + [
241 self.surface_runoff,
242 self.subsurface_runoff,
243 self.percolation,
244 ]:
245 tanks.end_timestep()
247 def get_surface(self, surface_):
248 """Return a surface from the list of surfaces by the 'surface' entry in the
249 surface. I.e., the name of the surface.
251 Args:
252 surface_ (str): Name of the surface
254 Returns:
255 surface (Surface): The first surface that matches the name
256 """
257 for surface in self.surfaces:
258 if surface.surface == surface_:
259 return surface
260 return None
262 def reinit(self):
263 """"""
264 self.end_timestep()
265 for surface in self.surfaces + [
266 self.surface_runoff,
267 self.subsurface_runoff,
268 self.percolation,
269 ]:
270 surface.reinit()
273class Surface(DecayTank):
274 """"""
276 def __init__(
277 self,
278 surface="",
279 area=1,
280 depth=1,
281 data_input_dict={},
282 pollutant_load={},
283 **kwargs,
284 ):
285 """A subclass of DecayTank. Each Surface is anticipated to represent a different
286 land cover type of a Land node. Besides functioning as a Tank, Surfaces have
287 three lists of functions (inflows, processes and outflows) where behaviour can
288 be added by appending new functions. We anticipate that customised surfaces
289 should be a subclass of Surface or its subclasses and add functions to these
290 lists. These lists are executed (inflows first, then processes, then outflows)
291 in the run function, which is called by the run function in Land. The lists must
292 return any model inflows or outflows as a VQIP for mass balance checking.
294 If a user wishes the DecayTank portion to be active, then can provide 'decays',
295 which are passed upwards (see wsimod/core/core.py/DecayObj for documentation)
297 Args:
298 surface (str, optional): String description of the surface type. Doesn't
299 serve a modelling purpose, just used for user reference. Defaults to ''.
300 area (float, optional): Area of surface. Defaults to 1. depth (float,
301 optional): Depth of tank (this has different physical
302 implications for different subclasses). Defaults to 1.
303 data_input_dict (dict, optional): Dictionary of data inputs relevant for
304 the surface (generally, deposition). Keys are tuples where first value
305 is the name of the variable to read from the dict and the second value
306 is the time. Note that this input should be specific to the surface, and
307 is not intended to be the same data input as for the land node. Also
308 note that with each surface having its own timeseries of data inputs,
309 this can take up a lot of memory, thus the default behavior is to have
310 this as monthly data where the time variable is a monthyear. Defaults to
311 {}.
312 pollutant_load (dict, optional): A dict of different pollutant amounts that
313 are deposited on the surface (units are mass per area per timestep).
314 Defaults to {}.
316 Key assumptions:
317 - Generic `Surface` that reads data and can apply simple forms of pollution
318 deposition.
319 - Formulated as a `Tank` object.
320 - Ammonia->Nitrite->Nitrate decay takes place if parameters describing this
321 process are provided in `decays` (see `core.py/DecayObj` for
322 transformation details).
324 Input data and parameter requirements:
325 - `data_input_dict` can contain a variety of pollutant deposition data.
326 `srp-dry` describes phosphate. `noy-dry` describes nitrogen as nitrates.
327 `nhx-dry` describes nitrogen as ammonia. `srp/noy/ nhx-wet` can also be
328 used to specify wet deposition. _Units_: kg/m2/timestep (data is read at
329 a monthly timestep)
330 """
331 # Assign parameters
332 self.depth = depth
333 self.data_input_dict = data_input_dict
334 self.surface = surface
335 self.pollutant_load = pollutant_load
336 # TODO this is a decaytank but growing surfaces don't have decay parameters...
337 # is it a problem.. we don't even take decays as an explicit argument and insert
338 # them in kwargs..
339 capacity = area * depth
340 # Parameters
341 super().__init__(capacity=capacity, area=area, **kwargs)
343 # Populate function lists TODO.. not sure why I have deposition but no
344 # precipitation here
345 if "nhx-dry" in set(x[0] for x in data_input_dict.keys()):
346 self.inflows = [self.atmospheric_deposition, self.precipitation_deposition]
347 else:
348 self.inflows = []
349 if len(self.pollutant_load) > 0:
350 self.inflows.append(self.simple_deposition)
351 self.processes = []
352 self.outflows = []
354 def run(self):
355 """Call run function (called from Land node)."""
356 if "nitrite" in constants.POLLUTANTS:
357 # Assume that if nitrite is modelled then nitrification is also modelled You
358 # will need ammonia->nitrite->nitrate decay to accurate simulate ammonia
359 # Thus, these decays are simulated here
361 # NOTE decay in a decaytank happens at start of timestep (confusingly) in
362 # the end_timestep function
363 self.storage["nitrate"] += self.total_decayed["nitrite"]
364 self.parent.running_inflow_mb["nitrate"] += self.total_decayed["nitrite"]
366 # Decayed ammonia becomes nitrite
367 self.storage["nitrite"] += self.total_decayed["ammonia"]
368 self.parent.running_inflow_mb["nitrite"] += self.total_decayed["ammonia"]
370 for f in self.inflows + self.processes + self.outflows:
371 # Iterate over function lists, updating mass balance
372 in_, out_ = f()
373 self.parent.running_inflow_mb = self.sum_vqip(
374 self.parent.running_inflow_mb, in_
375 )
376 self.parent.running_outflow_mb = self.sum_vqip(
377 self.parent.running_outflow_mb, out_
378 )
380 def get_data_input(self, var):
381 """Read data input from parent Land node (i.e., for precipitation/et0/temp).
383 Args:
384 var (str): Name of variable
386 Returns:
387 Data read
388 """
389 return self.parent.get_data_input(var)
391 def get_data_input_surface(self, var):
392 """Read data input from this surface's data_input_dict.
394 Args:
395 var (str): Name of variable
397 Returns:
398 Data read
399 """
400 return self.data_input_dict[(var, self.parent.monthyear)]
402 def dry_deposition_to_tank(self, vqip):
403 """Generic function for allocating dry pollution deposition to the surface.
404 Simply sends the pollution into the tank (some subclasses overwrite this
405 behaviour).
407 Args:
408 vqip (dict): A VQIP amount of dry deposition to send to tank
410 Returns:
411 vqip (dict): A VQIP amount of dry deposition that entered the tank (used
412 for mass balance checking)
413 """
414 # Default behaviour is to just enter the tank
415 _ = self.push_storage(vqip, force=True)
416 return vqip
418 def wet_deposition_to_tank(self, vqip):
419 """Generic function for allocating wet pollution deposition to the surface.
420 Simply sends the pollution into the tank (some subclasses overwrite this
421 behaviour).
423 Args:
424 vqip (dict): A VQIP amount of wet deposition to send to tank
426 Returns:
427 vqip (dict): A VQIP amount of wet deposition that entered the tank (used
428 for mass balance checking)
429 """
430 # Default behaviour is to just enter the tank
431 _ = self.push_storage(vqip, force=True)
432 return vqip
434 def simple_deposition(self):
435 """Inflow function to cause simple pollution deposition to occur, updating the
436 surface tank.
438 Returns:
439 (tuple): A tuple containing a VQIP amount for model inputs and outputs
440 for mass balance checking.
441 """
442 pollution = self.empty_vqip()
444 # Scale by area
445 for pol, item in self.pollutant_load.items():
446 if pol in constants.ADDITIVE_POLLUTANTS:
447 pollution[pol] = item * self.area
448 else:
449 pollution[pol] = item
450 pollution["volume"] = 0
452 # Update tank
453 _ = self.push_storage(pollution, force=True)
455 return (pollution, self.empty_vqip())
457 def atmospheric_deposition(self):
458 """Inflow function to cause dry atmospheric deposition to occur, updating the
459 surface tank.
461 Returns:
462 (tuple): A tuple containing a VQIP amount for model inputs and outputs
463 for mass balance checking.
464 """
465 # TODO double check units in preprocessing - is weight of N or weight of
466 # NHX/noy?
468 # Read data and scale
469 nhx = self.get_data_input_surface("nhx-dry") * self.area
470 noy = self.get_data_input_surface("noy-dry") * self.area
471 srp = self.get_data_input_surface("srp-dry") * self.area
473 # Assign pollutants
474 vqip = self.empty_vqip()
475 vqip["ammonia"] = nhx
476 vqip["nitrate"] = noy
477 vqip["phosphate"] = srp
479 # Update tank
480 in_ = self.dry_deposition_to_tank(vqip)
482 # Return mass balance
483 return (in_, self.empty_vqip())
485 def precipitation_deposition(self):
486 """Inflow function to cause wet precipitation deposition to occur, updating the
487 surface tank.
489 Returns:
490 (tuple): A tuple containing a VQIP amount for model inputs and outputs
491 for mass balance checking.
492 """
493 # TODO double check units - is weight of N or weight of NHX/noy?
495 # Read data and scale
496 nhx = self.get_data_input_surface("nhx-wet") * self.area
497 noy = self.get_data_input_surface("noy-wet") * self.area
498 srp = self.get_data_input_surface("srp-wet") * self.area
500 # Assign pollutants
501 vqip = self.empty_vqip()
502 vqip["ammonia"] = nhx
503 vqip["nitrate"] = noy
504 vqip["phosphate"] = srp
506 # Update tank
507 in_ = self.wet_deposition_to_tank(vqip)
509 # Return mass balance
510 return (in_, self.empty_vqip())
513class ImperviousSurface(Surface):
514 """"""
516 def __init__(self, pore_depth=0, et0_to_e=1, **kwargs):
517 """A surface to represent impervious surfaces that drain to storm sewers. Runoff
518 is generated by the surface tank overflowing, if a user wants all precipitation
519 to immediately go to runoff then they should reduce 'pore_depth', however
520 generally this is not what happens and a small (a few mm) depth should be
521 assigned to the tank. Also includes urban pollution deposition, though this will
522 only be mobilised if runoff occurs.
524 Note that the tank does not have a runoff coefficient because it doesn't make
525 sense from an integrated perspective. If a user wants to mimic runoff
526 coefficient-like behaviour, then they should reduce the ImperviousSurface tank
527 size, and increase other surfaces of the parent land node accordingly.
529 Args:
530 pore_depth (float, optional): The depth of the tank that must be exceeded
531 to generate runoff. Intended to represent the pores in ashpalt that
532 water accumulates in before flowing. Defaults to 0.
533 et0_to_e (float, optional): Multiplier applied to the parent's data
534 timeseries of et0 to determine how much evaporation takes place on the
535 ImperviousSurface. Defaults to 1.
536 """
537 # Assign parameters
538 self.et0_to_e = et0_to_e # Total evaporation
539 self.pore_depth = pore_depth
541 super().__init__(depth=pore_depth, **kwargs)
543 # Initialise state variables
544 self.evaporation = self.empty_vqip()
545 self.precipitation = self.empty_vqip()
547 # Populate function lists
548 self.inflows.append(self.precipitation_evaporation)
550 self.outflows.append(self.push_to_sewers)
552 def precipitation_evaporation(self):
553 """Inflow function that is a simple rainfall-evaporation model, updating the.
555 surface tank. All precipitation that is not evaporated is forced into the tank
556 (even though some of that will later be pushed to sewers) - this enables runoff
557 to mix with the accumulated pollutants in the surface pores.
559 Returns:
560 (tuple): A tuple containing a VQIP amount for model inputs and outputs
561 for mass balance checking.
562 """
563 # Read data in length units
564 precipitation_depth = self.get_data_input("precipitation")
565 evaporation_depth = self.get_data_input("et0") * self.et0_to_e
567 if precipitation_depth < evaporation_depth:
568 # No effective precipitation
569 net_precipitation = 0
571 # Calculate how much should be evaporated from pores
572 evaporation_from_pores = evaporation_depth - precipitation_depth
574 # Scale
575 evaporation_from_pores *= self.area
577 # Pull from tank
578 evaporation_from_pores = self.evaporate(evaporation_from_pores)
580 # Scale to get actual evaporation
581 total_evaporation = evaporation_from_pores + precipitation_depth * self.area
582 else:
583 # Effective precipitation
584 net_precipitation = precipitation_depth - evaporation_depth
586 # Scale
587 net_precipitation *= self.area
588 net_precipitation = self.v_change_vqip(self.empty_vqip(), net_precipitation)
590 # Assign a temperature value TODO how hot is rain? No idea... just going to
591 # use surface air temperature
592 net_precipitation["temperature"] = self.get_data_input("temperature")
594 # Update tank
595 _ = self.push_storage(net_precipitation, force=True)
596 total_evaporation = evaporation_depth * self.area
598 # Converrt to VQIP
599 self.evaporation = self.v_change_vqip(self.empty_vqip(), total_evaporation)
600 self.precipitation = self.v_change_vqip(
601 self.empty_vqip(), precipitation_depth * self.area
602 )
604 return (self.precipitation, self.evaporation)
606 def push_to_sewers(self):
607 """Outflow function that distributes ponded water (i.e., surface runoff) to the
608 parent node's attached sewers.
610 Returns:
611 (tuple): A tuple containing a VQIP amount for model inputs and outputs
612 for mass balance checking.
613 """
614 # Get runoff
615 surface_runoff = self.pull_ponded()
617 # Distribute TODO in cwsd_partition this is done with timearea
618 reply = self.parent.push_distributed(surface_runoff, of_type=["Sewer"])
620 # Update tank (forcing, because if the water can't go to the sewer, where else
621 # can it go)
622 _ = self.push_storage(reply, force=True)
623 # TODO... possibly this could flow to attached river or land nodes.. or other
624 # surfaces? I expect this doesn't matter for large scale models.. but may need
625 # to be revisited for detailed sewer models
627 # Return empty mass balance because outflows are handled by parent
628 return (self.empty_vqip(), self.empty_vqip())
631class PerviousSurface(Surface):
632 """"""
634 def __init__(
635 self,
636 depth=0.75,
637 total_porosity=0.4,
638 field_capacity=0.3,
639 wilting_point=0.12,
640 infiltration_capacity=0.5,
641 surface_coefficient=0.05,
642 percolation_coefficient=0.75,
643 et0_coefficient=0.5,
644 ihacres_p=10,
645 **kwargs,
646 ):
647 """A generic pervious surface that represents hydrology with the IHACRES model.
649 Args:
650 depth (float, optional): Soil tank (i.e., root) depth. Defaults to 0.75.
651 total_porosity (float, optional): The total porosity IHACRES parameter
652 (i.e., defines the total porouse volume of the soil - the maximum volume
653 of soil pores can contain when saturated). Defaults to 0.4.
654 field_capacity (float, optional): The field capacity IHACRES parameter
655 (i.e., when water content in the soil tank is above this value - flows
656 of any kind can be generated). Defaults to 0.3.
657 wilting_point (float, optional): The wilting point IHACRES parameter (i.e.,
658 when water content content in the soil tank is above this value - plants
659 can uptake water and evaporation from the soil tank can occur). Defaults
660 to 0.12.
661 infiltration_capacity (float, optional): Depth of water per day that can
662 enter the soil tank. Non infiltrated water will pond and travel as
663 surface runoff from the parent Land node. Defaults to 0.5.
664 surface_coefficient (float, optional): If flow is generated, the proportion
665 of flow that goes to surface runoff. Defaults to 0.05.
666 percolation_coefficient (float, optional): If flow is generated, then the
667 proportion of water that does not go to surface runoff that goes to
668 percolation (i.e., groundwater) - the remainder goes to subsurface
669 runoff. Defaults to 0.75.
670 et0_coefficient (float, optional): Convert between the parent nodes data
671 timeseries et0 - and potential evaptranspiration per unit area for this
672 surface. Defaults to=0.5,
673 ihacres_p (float, optional): The IHACRES p parameter. Unless it is an
674 ephemeral stream this parameter probably can stay high. Defaults to 10.
676 Key assumptions:
677 - In IHACRES, the maximum infiltration per time step is controlled by an
678 infiltration capacity, beyond which the precipitation will flow directly
679 as surface runoff.
680 - Evapotranspiration and effective precipitation are calculated based on
681 soil moisture content.
682 - Effective precipitation is then divided into percolation, surface runoff,
683 and subsurface runoff by multiplying the corresponding coefficient.
684 - Percolation, surface runoff, and subsurface runoff are sent into the
685 corresponding residence tanks for rounting to downstream.
686 - The mass of pollutants in soil water tank proportionately leaves the
687 soil water tank into the routing residence tanks. Evapotranspiration can
688 only bring out water, with pollutants left in the soil tank.
690 Input data and parameter requirements:
691 - Field capacity and wilting point.
692 _Units_: -, both should in [0-1], with field capacity > wilting point
693 - Infiltration capacity.
694 _Units_: m/day
695 - Surface, percolation coefficient.
696 _Units_: -, both should in [0-1]
697 - et0 coefficient.
698 _Units_: -
699 - ihacres_p.
700 _Units_: -
701 """
702 # Assign parameters (converting field capacity and wilting point to depth)
703 self.field_capacity = field_capacity
704 self.field_capacity_m = field_capacity * depth
705 self.wilting_point = wilting_point
706 self.wilting_point_m = wilting_point * depth
707 self.infiltration_capacity = infiltration_capacity
708 self.surface_coefficient = surface_coefficient
709 self.percolation_coefficient = percolation_coefficient
710 self.et0_coefficient = et0_coefficient
711 self.ihacres_p = ihacres_p
712 self.total_porosity = total_porosity
714 # Parameters to determine how to calculate the temperature of outflowing water
715 # TODO what should these params be?
716 self.soil_temp_w_prev = 0.1 # previous timestep weighting
717 self.soil_temp_w_air = 0.6 # air temperature weighting
718 self.soil_temp_w_deep = 0.1 # deep soil temperature weighting
719 self.soil_temp_deep = 10 # deep soil temperature
721 # IHACRES is a deficit not a tank, so doesn't really have a capacity in this
722 # way... and if it did.. I don't know if it would be the root depth
723 super().__init__(depth=depth * total_porosity, **kwargs)
725 # Calculate subsurface coefficient
726 self.subsurface_coefficient = 1 - self.percolation_coefficient
728 # Initiliase state variables
729 self.infiltration_excess = self.empty_vqip()
730 self.subsurface_flow = self.empty_vqip()
731 self.percolation = self.empty_vqip()
732 self.tank_recharge = 0
733 self.evaporation = self.empty_vqip()
734 self.precipitation = self.empty_vqip()
736 # Populate function lists
737 self.inflows.append(self.ihacres) # work out runoff
739 # TODO interception if I hate myself enough?
740 self.processes.append(
741 self.calculate_soil_temperature
742 ) # Calculate soil temp + dependence factor
744 # self.processes.append(self.decay) #apply generic decay (currently handled by
745 # decaytank at end of timestep) TODO decaytank uses air temperature not soil
746 # temperature... probably need to just give it the decay function
748 self.outflows.append(self.route)
750 def get_cmd(self):
751 """Calculate moisture deficit (i.e., the tank excess converted to depth).
753 Returns:
754 (float): current moisture deficit
755 """
756 return self.get_excess()["volume"] / self.area
758 def get_smc(self):
759 """Calculate moisture content (i.e., the tank volume converted to depth).
761 Returns:
762 (float): soil moisture content
763 """
764 # Depth of soil moisture
765 return self.storage["volume"] / self.area
767 def get_climate(self):
768 """
770 Returns:
772 """
773 precipitation_depth = self.get_data_input("precipitation")
774 evaporation_depth = self.get_data_input("et0") * self.et0_coefficient
775 return precipitation_depth, evaporation_depth
777 def ihacres(self):
778 """Inflow function that runs the IHACRES model equations, updates tanks, and
779 store flows in state variables (which are later sent to the parent land node in
780 the route function).
782 Returns:
783 (tuple): A tuple containing a VQIP amount for model inputs and outputs
784 for mass balance checking.
785 """
786 # Read data (leave in depth units since that is what IHACRES equations are in)
787 precipitation_depth, evaporation_depth = self.get_climate()
788 temperature = self.get_data_input("temperature")
790 # Apply infiltration
791 infiltrated_precipitation = min(precipitation_depth, self.infiltration_capacity)
792 infiltration_excess = max(precipitation_depth - infiltrated_precipitation, 0)
794 # Get current moisture deficit
795 current_moisture_deficit_depth = self.get_cmd()
797 # IHACRES equations (we do (depth - wilting_point_m or field capacity) to
798 # convert from a deficit to storage tank)
799 evaporation = evaporation_depth * min(
800 1,
801 exp(
802 2
803 * (
804 1
805 - current_moisture_deficit_depth
806 / (self.depth - self.wilting_point_m)
807 )
808 ),
809 )
810 outflow = infiltrated_precipitation * (
811 1
812 - min(
813 1,
814 (current_moisture_deficit_depth / (self.depth - self.field_capacity_m))
815 ** self.ihacres_p,
816 )
817 )
819 # Can't evaporate more than available moisture (presumably the IHACRES equation
820 # prevents this ever being needed)
821 evaporation = min(evaporation, precipitation_depth + self.get_smc())
823 # Scale to volumes and apply proportions to work out percolation/surface
824 # runoff/subsurface runoff
825 surface = outflow * self.surface_coefficient * self.area
826 percolation = (
827 outflow
828 * (1 - self.surface_coefficient)
829 * self.percolation_coefficient
830 * self.area
831 )
832 subsurface_flow = (
833 outflow
834 * (1 - self.surface_coefficient)
835 * self.subsurface_coefficient
836 * self.area
837 )
838 tank_recharge = (infiltrated_precipitation - evaporation - outflow) * self.area
839 infiltration_excess *= self.area
840 infiltration_excess += surface
841 evaporation *= self.area
842 precipitation = precipitation_depth * self.area
844 # Mix in tank to calculate pollutant concentrations
845 total_water_passing_through_soil_tank = (
846 tank_recharge + subsurface_flow + percolation
847 )
849 if total_water_passing_through_soil_tank > 0:
850 # Net effective preipitation
851 total_water_passing_through_soil_tank = self.v_change_vqip(
852 self.empty_vqip(), total_water_passing_through_soil_tank
853 )
854 # Assign a temperature before sending into tank
855 total_water_passing_through_soil_tank["temperature"] = temperature
856 # Assign to tank
857 _ = self.push_storage(total_water_passing_through_soil_tank, force=True)
858 # Pull flows (which now have nonzero pollutant concentrations)
859 subsurface_flow = self.pull_storage({"volume": subsurface_flow})
860 percolation = self.pull_storage({"volume": percolation})
861 else:
862 # No net effective precipitation (instead evaporation occurs)
863 evap = self.evaporate(-total_water_passing_through_soil_tank)
864 subsurface_flow = self.empty_vqip()
865 percolation = self.empty_vqip()
867 if (
868 abs(
869 evap
870 + infiltrated_precipitation * self.area
871 - evaporation
872 - infiltration_excess
873 )
874 > constants.FLOAT_ACCURACY
875 ):
876 print(
877 "inaccurate evaporation calculation of {0}".format(
878 abs(
879 evap
880 + infiltrated_precipitation * self.area
881 - evaporation
882 - infiltration_excess
883 )
884 )
885 )
887 # TODO saturation excess (think it should just be 'pull_ponded' presumably in
888 # net effective precipitation? )
890 # Convert to VQIPs
891 infiltration_excess = self.v_change_vqip(self.empty_vqip(), infiltration_excess)
892 infiltration_excess["temperature"] = temperature
893 precipitation = self.v_change_vqip(self.empty_vqip(), precipitation)
894 evaporation = self.v_change_vqip(self.empty_vqip(), evaporation)
896 # Track flows (these are sent onwards in the route function)
897 self.infiltration_excess = infiltration_excess
898 self.subsurface_flow = subsurface_flow
899 self.percolation = percolation
900 self.tank_recharge = tank_recharge
901 self.evaporation = evaporation
902 self.precipitation = precipitation
904 # Mass balance
905 in_ = precipitation
906 out_ = evaporation
908 return (in_, out_)
910 def route(self):
911 """An outflow function that sends percolation, subsurface runoff and surface
912 runoff to their respective tanks in the parent land node.
914 Returns:
915 (tuple): A tuple containing a VQIP amount for model inputs and outputs
916 for mass balance checking.
917 """
918 self.parent.surface_runoff.push_storage(self.infiltration_excess, force=True)
919 self.parent.subsurface_runoff.push_storage(self.subsurface_flow, force=True)
920 self.parent.percolation.push_storage(self.percolation, force=True)
922 return (self.empty_vqip(), self.empty_vqip())
924 def calculate_soil_temperature(self):
925 """Process function that calculates soil temperature based on a weighted.
927 average. This equation is from Lindstrom, Bishop & Lofvenius (2002),
928 hydrological processes - but it is not clear what the parameters should be.
930 Returns:
931 (tuple): A tuple containing a VQIP amount for model inputs and outputs
932 for mass balance checking.
933 """
934 auto = self.storage["temperature"] * self.soil_temp_w_prev
935 air = self.get_data_input("temperature") * self.soil_temp_w_air
936 total_weight = (
937 self.soil_temp_w_air + self.soil_temp_w_deep + self.soil_temp_w_prev
938 )
939 self.storage["temperature"] = (
940 auto + air + self.soil_temp_deep * self.soil_temp_w_deep
941 ) / total_weight
943 return (self.empty_vqip(), self.empty_vqip())
946class GrowingSurface(PerviousSurface):
947 """"""
949 def __init__(
950 self,
951 rooting_depth=1,
952 ET_depletion_factor=1,
953 crop_factor_stages=[1, 1],
954 crop_factor_stage_dates=[0, 365],
955 sowing_day=1,
956 harvest_day=365,
957 initial_soil_storage=None,
958 **kwargs,
959 ):
960 """Extensive surface subclass that implements the CatchWat equations (Liu,
961 Dobson & Mijic (2022) Science of the total environment), which in turn are
962 primarily based on FAO document:
963 https://www.fao.org/3/x0490e/x0490e0ehtm#soil%20water%20availability. This
964 surface is a pervious surface that also has things that grow on it. This
965 behaviour includes soil nutrient pools, crop planting/harvest calendars,
966 erosion, crop behaviour.
968 A key complexity of this surface is the nutrient pool (see wsimod/nodes/
969 nutrient_pool.py), which is a class that tracks the amount of phosphorus and
970 nitrogen in different states and performs transformations that occur in the
971 phosphorus/nitrogen cycle. It is assumed that the phosphate/nitrate/nitrite/
972 ammonia amounts in this Surface tank should track the dissolved inorganic pool
973 in the nutrient pool. Meanwhile, the org-phosphorus/org-nitrogen amounts in this
974 tank should track the dissolved organic pool in the nutrient pool. The total
975 amount of pollutants that enter this tank may not be the same as the total
976 amount that leave, because pollutants are transformed between inorganic/ organic
977 and between wet/dry states - these transformations are accounted for in mass
978 balance.
980 For users to quickly enable/disable these nutrient processes, which are
981 computationally intensive (in current case studies they account for about half
982 of the total runtime), they are only active if 'nitrate' is one of the modelled
983 pollutants. Note that the code will not check if nitrite/phosphate/
984 org-phosphorus/org-nitrogen/ammonia are also included, but they should be if
985 nitrate is included and otherwise the code will crash with a key error.
987 Args:
988 rooting_depth (float, optional): Depth of the soil tank (i.e., how deep do
989 crop roots go). Defaults to 1.
990 ET_depletion_factor (float, optional): Average fraction of soil that can be
991 depleted from the root zone before moisture stress (reduction in ET)
992 occurs. Defaults to 1.
993 crop_factor_stages (list, optional): Crop factor is a multiplier on et0,
994 more grown plants have higher transpiration and higher crop factors.
995 This list shows changing crop factor at different times of year in
996 relation to crop_factor_stage_dates. See wsimod/preprocessing/
997 england_data_formatting.py/format_surfaces for further details on
998 formulating these - since the interpolation used to find crop_factors in
999 between the given values in the list is a bit involved. Defaults to
1000 [1,1].
1001 crop_factor_stage_dates (list, optional): Dates associated with
1002 crop_factor_stages. Defaults to [0, 365].
1003 sowing_day (int, optional): day of year that crops are sown. Defaults to 1.
1004 harvest_day (int, optional): day of year that crops are harvest. Defaults
1005 to 365.
1006 initial_soil_storage (dict or float, optional): Initial mass of solid
1007 pollutants
1008 in the soil nutrient pools (fast and adsorbed inorganic pools)
1010 Key assumptions:
1011 - In the soil water module, crop stages and crop coefficients control the
1012 evapotranspiration.
1013 - Fertiliser and manure application are the major source of soil nutrients,
1014 which are added
1015 into soil nutrient pools, including dissovled inorganic, dissolved
1016 organic, fast and humus for both nitrogen and phosphorus.
1017 - Nutrient transformation processes in soil are simulated, including fluxes
1018 between the soil
1019 nutrient pools, denitrification for nitrogen, adsorption/desorption for
1020 phosphorus. These processes are affected by temperature and soil
1021 moisture.
1022 - Crop uptake of nutrients are simulated based on crop stages, which is
1023 different for spring-sown
1024 and autumn-sown crops.
1025 - Soil erosion from the growing surface is simulated as one of the major
1026 sources of suspended solids
1027 in rivers, which is mainly affected by rainfall energy and crop/ground
1028 cover. Phosphorus will also be eroded along with the soil particles, in
1029 both adsorbed inorganic and humus form.
1031 Input data and parameter requirements:
1032 - `data_input_dict` can contain a variety of pollutant deposition data.
1033 `srp-fertiliser` describes phosphate. `noy-fertiliser` describes
1034 nitrogen as nitrates. `nhx-fertiliser` describes nitrogen as ammonia.
1035 `srp/noy/ nhx-manure` can also be used to specify manure application.
1036 _Units_: kg/m2/timestep (data is read at a monthly timestep)
1037 - Rooting depth.
1038 _Units_: m
1039 - Evapotranspiration depletion factor.
1040 _Units_: -
1041 - Sowing day, harvest day and crop calendars.
1042 _Units_: day number in Julian calendar
1043 - Crop factor.
1044 _Units_: -
1045 - Initial storage for solid pollutants.
1046 _Units_: kg
1048 """
1049 # Crop factors (set when creating object)
1050 self.ET_depletion_factor = (
1051 ET_depletion_factor # To do with water availability, p from FAOSTAT
1052 )
1053 self.rooting_depth = (
1054 rooting_depth # maximum depth that plants can absorb, Zr from FAOSTAT
1055 )
1056 depth = rooting_depth
1058 # Crop parameters
1059 self.crop_cover_max = 0.9 # [-] 0~1
1060 self.ground_cover_max = 0.3 # [-]
1061 # TODO... really I should just have this as an annual profile parameter and do
1062 # away with interpolation etc.
1063 self.crop_factor_stages = crop_factor_stages
1064 self.crop_factor_stage_dates = crop_factor_stage_dates
1065 self.sowing_day = sowing_day
1066 self.harvest_day = harvest_day
1068 # Soil moisture dependence parameters
1069 self.satact = 0.6 # [-] for calculating soil_moisture_dependence_factor
1070 self.thetaupp = 0.12 # [-] for calculating soil_moisture_dependence_factor
1071 self.thetalow = 0.08 # [-] for calculating soil_moisture_dependence_factor
1072 self.thetapow = 1 # [-] for calculating soil_moisture_dependence_factorself.
1073 # satact = 0.6 # [-] for calculating soil_moisture_dependence_factor
1075 # Crop uptake parameters
1076 self.uptake1 = (
1077 15 # [g/m2/y] shape factor for crop (Dissolved) Inorganic nitrogen uptake
1078 )
1079 self.uptake2 = (
1080 1 # [-] shape factor for crop (Dissolved) Inorganic nitrogen uptake
1081 )
1082 self.uptake3 = (
1083 0.02 # [1/day] shape factor for crop (Dissolved) Inorganic nitrogen uptake
1084 )
1085 self.uptake_PNratio = 1 / 7.2 # [-] P:N during crop uptake
1087 # Erosion parameters
1088 self.erodibility = 0.0025 # [g * d / (J * mm)]
1089 self.sreroexp = 1.2 # [-] surface runoff erosion exponent
1090 self.cohesion = 1 # [kPa]
1091 self.slope = 5 # [-] every 100
1092 self.srfilt = (
1093 0.7 # [-] ratio of eroded sediment left in surface runoff after filtration
1094 )
1095 self.macrofilt = 0.01 # [-] ratio of eroded sediment left in subsurface flow
1096 # after filtration
1098 # Denitrification parameters
1099 self.limpar = 0.7 # [-] above which denitrification begins
1100 self.exppar = 2.5 # [-] exponential parameter for
1101 # soil_moisture_dependence_factor_exp calculation
1102 self.hsatINs = 1 # [mg/l] for calculation of half-saturation concentration
1103 # dependence factor
1104 self.denpar = 0.015 # [-] denitrification rate coefficient
1106 # Adsorption parameters
1107 self.adosorption_nr_limit = 0.00001
1108 self.adsorption_nr_maxiter = 20
1109 self.kfr = 153.7 # [litter/kg] freundlich adsorption isoterm
1110 self.nfr = 1 / 2.6 # [-] freundlich exponential coefficient
1111 self.kadsdes = 0.03 # [1/day] adsorption/desorption coefficient
1113 # Other soil parameters
1114 self.bulk_density = 1300 # [kg/m3]
1115 super().__init__(depth=depth, **kwargs)
1117 # Infer basic sow/harvest calendar
1118 self.harvest_sow_calendar = [
1119 0,
1120 self.sowing_day,
1121 self.harvest_day,
1122 self.harvest_day + 1,
1123 365,
1124 ]
1125 self.ground_cover_stages = [0, 0, self.ground_cover_max, 0, 0]
1126 self.crop_cover_stages = [0, 0, self.crop_cover_max, 0, 0]
1128 # Use day number of 181 to indicate autumn-sown (from HYPE)
1129 if self.sowing_day > 181:
1130 self.autumn_sow = True
1131 else:
1132 self.autumn_sow = False
1134 # State variables
1135 self.days_after_sow = None
1136 self.crop_cover = 0
1137 self.ground_cover = 0
1138 self.crop_factor = 0
1139 self.et0_coefficient = 1
1141 # Calculate parameters based on capacity/wp
1142 self.total_available_water = self.field_capacity_m - self.wilting_point_m
1143 if self.total_available_water < 0:
1144 print("warning: TAW < 0...")
1145 self.readily_available_water = (
1146 self.total_available_water * self.ET_depletion_factor
1147 )
1149 # Initiliase nutrient pools
1150 self.nutrient_pool = NutrientPool()
1152 self.inflows.insert(0, self.calc_crop_cover)
1153 if "nitrate" in constants.POLLUTANTS:
1154 # Populate function lists
1155 self.inflows.append(self.effective_precipitation_flushing)
1156 self.inflows.append(self.fertiliser)
1157 self.inflows.append(self.manure)
1158 # self.inflows.append(self.residue)
1160 self.processes.append(self.calc_temperature_dependence_factor)
1161 self.processes.append(self.calc_soil_moisture_dependence_factor)
1162 self.processes.append(self.soil_pool_transformation)
1163 self.processes.append(self.calc_crop_uptake)
1165 # TODO possibly move these into nutrient pool
1166 self.processes.append(self.erosion)
1167 self.processes.append(self.denitrification)
1168 self.processes.append(self.adsorption)
1170 # Reflect initial water concentration in dissolved nutrient pools
1171 self.nutrient_pool.dissolved_inorganic_pool.storage["P"] = self.storage[
1172 "phosphate"
1173 ]
1174 self.nutrient_pool.dissolved_inorganic_pool.storage["N"] = (
1175 self.storage["nitrate"]
1176 + self.storage["ammonia"]
1177 + self.storage["nitrite"]
1178 )
1179 self.nutrient_pool.dissolved_organic_pool.storage["P"] = self.storage[
1180 "org-phosphorus"
1181 ]
1182 self.nutrient_pool.dissolved_organic_pool.storage["N"] = self.storage[
1183 "org-nitrogen"
1184 ]
1185 if initial_soil_storage:
1186 self.initial_soil_storage = initial_soil_storage
1187 # Reflect initial nutrient stores in solid nutrient pools
1188 self.nutrient_pool.adsorbed_inorganic_pool.storage[
1189 "P"
1190 ] = initial_soil_storage["phosphate"]
1191 self.nutrient_pool.adsorbed_inorganic_pool.storage["N"] = (
1192 initial_soil_storage["ammonia"]
1193 + initial_soil_storage["nitrate"]
1194 + initial_soil_storage["nitrite"]
1195 )
1196 self.nutrient_pool.fast_pool.storage["N"] = initial_soil_storage[
1197 "org-nitrogen"
1198 ]
1199 self.nutrient_pool.fast_pool.storage["P"] = initial_soil_storage[
1200 "org-phosphorus"
1201 ]
1203 def pull_storage(self, vqip):
1204 """Pull water from the surface, updating the surface storage VQIP. Nutrient pool
1205 pollutants (nitrate/nitrite/ammonia/phosphate/org- phosphorus/ org-nitrogen) are
1206 removed in proportion to their amounts in the dissolved nutrient pools, if they
1207 are simulated. Other pollutants are removed in proportion to their amount in the
1208 surface tank.
1210 Args:
1211 vqip (dict): VQIP amount to be pulled, (only 'volume' key is needed)
1213 Returns:
1214 reply (dict): A VQIP amount successfully pulled from the tank
1215 """
1216 if self.storage["volume"] == 0:
1217 return self.empty_vqip()
1219 # Adjust based on available volume
1220 reply = min(vqip["volume"], self.storage["volume"])
1222 # Update reply to vqip (get concentration for non-nutrients)
1223 reply = self.v_change_vqip(self.storage, reply)
1225 if "nitrate" in constants.POLLUTANTS:
1226 # Update nutrient pool and get concentration for nutrients
1227 prop = reply["volume"] / self.storage["volume"]
1228 nutrients = self.nutrient_pool.extract_dissolved(prop)
1229 reply["nitrate"] = (
1230 nutrients["inorganic"]["N"]
1231 * self.storage["nitrate"]
1232 / (
1233 self.storage["nitrate"]
1234 + self.storage["ammonia"]
1235 + self.storage["nitrite"]
1236 )
1237 )
1238 reply["ammonia"] = (
1239 nutrients["inorganic"]["N"]
1240 * self.storage["ammonia"]
1241 / (
1242 self.storage["nitrate"]
1243 + self.storage["ammonia"]
1244 + self.storage["nitrite"]
1245 )
1246 )
1247 reply["nitrite"] = (
1248 nutrients["inorganic"]["N"]
1249 * self.storage["nitrite"]
1250 / (
1251 self.storage["nitrate"]
1252 + self.storage["ammonia"]
1253 + self.storage["nitrite"]
1254 )
1255 )
1256 reply["phosphate"] = nutrients["inorganic"]["P"]
1257 reply["org-phosphorus"] = nutrients["organic"]["P"]
1258 reply["org-nitrogen"] = nutrients["organic"]["N"]
1260 # Extract from storage
1261 self.storage = self.extract_vqip(self.storage, reply)
1263 return reply
1265 def quick_interp(self, x, xp, yp):
1266 """A simple version of np.interp to intepolate crop information on the fly.
1268 Args:
1269 x (int): Current time (i.e., day of year) xp (list): Predefined times (i.e.,
1270 list of days of year) yp (list): Predefined values associated with xp
1272 Returns:
1273 y (float): Interpolated value for current time
1274 """
1275 x_ind = bisect_left(xp, x)
1276 x_left = xp[x_ind - 1]
1277 x_right = xp[x_ind]
1278 dif = x - x_left
1279 y_left = yp[x_ind - 1]
1280 y_right = yp[x_ind]
1281 y = y_left + (y_right - y_left) * dif / (x_right - x_left)
1282 return y
1284 def calc_crop_cover(self):
1285 """Process function that calculates how much crop cover there is, assigns
1286 whether crops are sown/harvested, and calculates et0_coefficient based on growth
1287 stage of crops.
1289 Returns:
1290 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1291 for mass balance checking.
1292 """
1293 # Get current day of year
1294 doy = self.parent.t.dayofyear
1296 if self.parent.t.is_leap_year:
1297 # Hacky way to handle leap years
1298 if doy > 59:
1299 doy -= 1
1301 if self.days_after_sow is None:
1302 if self.parent.t.dayofyear == self.sowing_day:
1303 # sow
1304 self.days_after_sow = 0
1305 else:
1306 if self.parent.t.dayofyear == self.harvest_day:
1307 # harvest
1308 self.days_after_sow = None
1309 self.crop_factor = self.crop_factor_stages[0]
1310 self.crop_cover = 0
1311 self.ground_cover = 0
1312 else:
1313 # increment days since sow
1314 self.days_after_sow += 1
1316 # Calculate relevant parameters
1317 self.crop_factor = self.quick_interp(
1318 doy, self.crop_factor_stage_dates, self.crop_factor_stages
1319 )
1320 if self.days_after_sow:
1321 # Move outside of this if, if you want nonzero crop/ground cover outside of
1322 # season
1323 self.crop_cover = self.quick_interp(
1324 doy, self.harvest_sow_calendar, self.crop_cover_stages
1325 )
1326 self.ground_cover = self.quick_interp(
1327 doy, self.harvest_sow_calendar, self.ground_cover_stages
1328 )
1330 root_zone_depletion = max(self.field_capacity_m - self.get_smc(), 0)
1331 if root_zone_depletion < self.readily_available_water:
1332 crop_water_stress_coefficient = 1
1333 else:
1334 crop_water_stress_coefficient = max(
1335 0,
1336 (self.total_available_water - root_zone_depletion)
1337 / ((1 - self.ET_depletion_factor) * self.total_available_water),
1338 )
1340 self.et0_coefficient = crop_water_stress_coefficient * self.crop_factor
1342 return (self.empty_vqip(), self.empty_vqip())
1344 def adjust_vqip_to_liquid(self, vqip, deposition, in_):
1345 """Function to interoperate between surface tank and nutrient pool. Most
1346 depositions are given in terms of ammonia/nitrate/phosphate - they are then
1347 aggregated to total N or P to enter the nutrient pools. Depending on the source
1348 of deposition these may transform (e.g., some go to dissolved and some to
1349 solids) upon entering the nutrient pool. To reflect these transformations in the
1350 soil tank, the amounts entering the soil tank are adjusted proportionately.
1352 Args:
1353 vqip (dict): A VQIP amount of pollutants originally intended to enter the
1354 soil tank
1355 deposition (dict): A dict with nutrients (N and P) as keys, showing the
1356 total amount of nutrients entering the nutrient pool
1357 in_ (dict): A dict with nutrients as keys, showing the updated amount of
1358 nutrients that entered the nutrient pool as dissolved pollutants
1360 Returns:
1361 vqip (dict): A VQIP amount of pollutants that have been scaled to account
1362 for nutrient pool transformations
1363 """
1364 if "nitrate" in constants.POLLUTANTS:
1365 if deposition["N"] > 0:
1366 vqip["nitrate"] *= in_["N"] / deposition["N"]
1367 vqip["ammonia"] *= in_["N"] / deposition["N"]
1368 vqip["org-nitrogen"] *= in_["N"] / deposition["N"]
1369 if deposition["P"] > 0:
1370 vqip["phosphate"] *= in_["P"] / deposition["P"]
1371 vqip["org-phosphorus"] *= in_["P"] / deposition["P"]
1373 return vqip
1375 def effective_precipitation_flushing(self):
1376 """Remove the nutrients brought out by effective precipitation, which is surface
1377 runoff, subsurface runoff, and percolation, from the nutrients pool.
1379 Returns:
1380 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1381 for mass balance checking.
1382 """
1383 # inorganic
1384 out = self.nutrient_pool.get_empty_nutrient()
1385 out["N"] = (
1386 self.subsurface_flow["ammonia"]
1387 + self.subsurface_flow["nitrite"]
1388 + self.subsurface_flow["nitrate"]
1389 + self.percolation["ammonia"]
1390 + self.percolation["nitrite"]
1391 + self.percolation["nitrate"]
1392 + self.infiltration_excess["ammonia"]
1393 + self.infiltration_excess["nitrite"]
1394 + self.infiltration_excess["nitrate"]
1395 ) # TODO what happens if infiltration excess (the real part) has pollutants?
1396 out["P"] = (
1397 self.subsurface_flow["phosphate"]
1398 + self.percolation["phosphate"]
1399 + self.infiltration_excess["phosphate"]
1400 )
1401 self.nutrient_pool.dissolved_inorganic_pool.extract(out)
1403 # organic
1404 out = self.nutrient_pool.get_empty_nutrient()
1405 out["N"] = (
1406 self.subsurface_flow["org-nitrogen"]
1407 + self.percolation["org-nitrogen"]
1408 + self.infiltration_excess["org-nitrogen"]
1409 )
1410 out["P"] = (
1411 self.subsurface_flow["org-phosphorus"]
1412 + self.percolation["org-phosphorus"]
1413 + self.infiltration_excess["org-phosphorus"]
1414 )
1415 self.nutrient_pool.dissolved_organic_pool.extract(out)
1417 return (self.empty_vqip(), self.empty_vqip())
1419 def fertiliser(self):
1420 """Read, scale and allocate fertiliser, updating the tank.
1422 Returns:
1423 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1424 for mass balance checking.
1425 """
1426 # TODO tidy up fertiliser/manure/residue/deposition once preprocessing is sorted
1428 # Scale for surface
1429 nhx = self.get_data_input_surface("nhx-fertiliser") * self.area
1430 noy = self.get_data_input_surface("noy-fertiliser") * self.area
1431 srp = self.get_data_input_surface("srp-fertiliser") * self.area
1433 # Update as VQIP
1434 vqip = self.empty_vqip()
1435 vqip["ammonia"] = nhx
1436 vqip["nitrate"] = noy
1437 vqip["phosphate"] = srp
1439 # Enter nutrient pool
1440 deposition = self.nutrient_pool.get_empty_nutrient()
1441 deposition["N"] = vqip["nitrate"] + vqip["ammonia"]
1442 deposition["P"] = vqip["phosphate"]
1443 in_ = self.nutrient_pool.allocate_fertiliser(deposition)
1445 # Update tank
1446 vqip = self.adjust_vqip_to_liquid(vqip, deposition, in_)
1447 self.push_storage(vqip, force=True)
1449 return (vqip, self.empty_vqip())
1451 def manure(self):
1452 """Read, scale and allocate manure, updating the tank.
1454 Returns:
1455 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1456 for mass balance checking.
1457 """
1458 # Scale for surface
1459 nhx = self.get_data_input_surface("nhx-manure") * self.area
1460 noy = self.get_data_input_surface("noy-manure") * self.area
1461 srp = self.get_data_input_surface("srp-manure") * self.area
1463 # Formulate as VQIP
1464 vqip = self.empty_vqip()
1465 vqip["ammonia"] = nhx
1466 vqip["nitrate"] = noy
1467 vqip["phosphate"] = srp
1469 # Enter nutrient pool
1470 deposition = self.nutrient_pool.get_empty_nutrient()
1471 deposition["N"] = vqip["nitrate"] + vqip["ammonia"]
1472 deposition["P"] = vqip["phosphate"]
1473 in_ = self.nutrient_pool.allocate_manure(deposition)
1475 # Update tank
1476 vqip = self.adjust_vqip_to_liquid(vqip, deposition, in_)
1478 self.push_storage(vqip, force=True)
1480 return (vqip, self.empty_vqip())
1482 def residue(self):
1483 """Read, scale and allocate residue, updating the tank (NOT CURRENTLY USED
1484 BECAUSE NO DATA SOURCES FOR RESIDUE CAN BE IDENTIFIED).
1486 Returns:
1487 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1488 for mass balance checking.
1489 """
1490 nhx = self.get_data_input_surface("nhx-residue") * self.area
1491 noy = self.get_data_input_surface("noy-residue") * self.area
1492 srp = self.get_data_input_surface("srp-residue") * self.area
1494 vqip = self.empty_vqip()
1495 vqip["ammonia"] = nhx * self.nutrient_pool.fraction_residue_to_fast["N"]
1496 vqip["nitrate"] = noy * self.nutrient_pool.fraction_residue_to_fast["N"]
1497 vqip["org-nitrogen"] = (
1498 nhx + noy
1499 ) * self.nutrient_pool.fraction_residue_to_humus["N"]
1500 vqip["phosphate"] = srp * self.nutrient_pool.fraction_residue_to_fast["P"]
1501 vqip["org-phosphorus"] = srp * self.nutrient_pool.fraction_residue_to_humus["P"]
1503 deposition = self.nutrient_pool.get_empty_nutrient()
1504 deposition["N"] = vqip["nitrate"] + vqip["ammonia"] + vqip["org-nitrogen"]
1505 deposition["P"] = vqip["phosphate"] + vqip["org-phosphorus"]
1507 in_ = self.nutrient_pool.allocate_residue(deposition)
1508 vqip = self.adjust_vqip_to_liquid(vqip, deposition, in_)
1510 self.push_storage(vqip, force=True)
1512 return (vqip, self.empty_vqip())
1514 def soil_pool_transformation(self):
1515 """A process function that run transformation functions in the nutrient pool and
1516 updates the pollutant concentrations in the surface tank.
1518 Returns:
1519 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1520 for mass balance checking.
1521 """
1522 # Initialise mass balance tracking variables
1523 in_ = self.empty_vqip()
1524 out_ = self.empty_vqip()
1526 # Get proportion of nitrogen that is nitrate in the soil tank NOTE ignores
1527 # nitrite - couldn't find enough information on it
1528 nitrate_proportion = self.storage["nitrate"] / (
1529 self.storage["nitrate"] + self.storage["ammonia"]
1530 )
1532 # Run soil pool functions
1533 (
1534 increase_in_dissolved_inorganic,
1535 increase_in_dissolved_organic,
1536 ) = self.nutrient_pool.soil_pool_transformation()
1538 # Update tank and mass balance TODO .. there is definitely a neater way to write
1539 # this
1540 if increase_in_dissolved_inorganic["N"] > 0:
1541 # Increase in inorganic nitrogen, rescale back to nitrate and ammonia
1542 in_["nitrate"] = increase_in_dissolved_inorganic["N"] * nitrate_proportion
1543 in_["ammonia"] = increase_in_dissolved_inorganic["N"] * (
1544 1 - nitrate_proportion
1545 )
1546 else:
1547 # Decrease in inorganic nitrogen, rescale back to nitrate and ammonia
1548 out_["nitrate"] = -increase_in_dissolved_inorganic["N"] * nitrate_proportion
1549 out_["ammonia"] = -increase_in_dissolved_inorganic["N"] * (
1550 1 - nitrate_proportion
1551 )
1553 if increase_in_dissolved_organic["N"] > 0:
1554 # Increase in organic nitrogen
1555 in_["org-nitrogen"] = increase_in_dissolved_organic["N"]
1556 else:
1557 # Decrease in organic nitrogen
1558 out_["org-nitrogen"] = -increase_in_dissolved_organic["N"]
1560 if increase_in_dissolved_inorganic["P"] > 0:
1561 # Increase in inorganic phosphate
1562 in_["phosphate"] = increase_in_dissolved_inorganic["P"]
1563 else:
1564 # Decrease in inorganic phosphate
1565 out_["phosphate"] = -increase_in_dissolved_inorganic["P"]
1567 if increase_in_dissolved_organic["P"] > 0:
1568 # Increase in organic phosphorus
1569 in_["org-phosphorus"] = increase_in_dissolved_organic["P"]
1570 else:
1571 # Decrease in organic phosphorus
1572 out_["org-phosphorus"] = -increase_in_dissolved_organic["P"]
1574 # Update tank with inputs/outputs of pollutants
1575 _ = self.push_storage(in_, force=True)
1576 out2_ = self.pull_pollutants(out_)
1578 if not self.compare_vqip(out_, out2_):
1579 print("nutrient pool not tracking soil tank")
1581 return (in_, out_)
1583 def calc_temperature_dependence_factor(self):
1584 """Process function that calculates the temperature dependence factor for the
1585 nutrient pool (which impacts soil pool transformations).
1587 Returns:
1588 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1589 for mass balance checking.
1590 """
1591 # Parameters/equations from HYPE documentation
1592 if self.storage["temperature"] > 5:
1593 temperature_dependence_factor = 2 ** (
1594 (self.storage["temperature"] - 20) / 10
1595 )
1596 elif self.storage["temperature"] > 0:
1597 temperature_dependence_factor = self.storage["temperature"] / 5
1598 else:
1599 temperature_dependence_factor = 0
1600 self.nutrient_pool.temperature_dependence_factor = temperature_dependence_factor
1601 return (self.empty_vqip(), self.empty_vqip())
1603 def calc_soil_moisture_dependence_factor(self):
1604 """Process function that calculates the soil moisture dependence factor for the
1605 nutrient pool (which impacts soil pool transformations).
1607 Returns:
1608 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1609 for mass balance checking.
1610 """
1611 # Parameters/equations from HYPE documentation
1612 current_soil_moisture = self.get_smc()
1613 if current_soil_moisture >= self.field_capacity_m:
1614 self.nutrient_pool.soil_moisture_dependence_factor = self.satact
1615 elif current_soil_moisture <= self.wilting_point_m:
1616 self.nutrient_pool.soil_moisture_dependence_factor = 0
1617 else:
1618 fc_diff = self.field_capacity_m - current_soil_moisture
1619 fc_comp = (fc_diff / (self.thetaupp * self.rooting_depth)) ** self.thetapow
1620 fc_comp = (1 - self.satact) * fc_comp + self.satact
1621 wp_diff = current_soil_moisture - self.wilting_point_m
1622 wp_comp = (wp_diff / (self.thetalow * self.rooting_depth)) ** self.thetapow
1623 self.nutrient_pool.soil_moisture_dependence_factor = min(
1624 1, wp_comp, fc_comp
1625 )
1626 return (self.empty_vqip(), self.empty_vqip())
1628 def calc_crop_uptake(self):
1629 """Process function that calculates how much nutrient crops uptake and updates
1630 nutrient pool and surface tank.
1632 Returns:
1633 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1634 for mass balance checking.
1635 """
1636 # Parameters/equations from HYPE documentation
1638 # Initialise
1639 N_common_uptake = 0
1640 P_common_uptake = 0
1642 if self.days_after_sow:
1643 # If there are crops
1645 days_after_sow = self.days_after_sow
1647 if self.autumn_sow:
1648 temp_func = max(0, min(1, (self.storage["temperature"] - 5) / 20))
1649 days_after_sow -= 25 # Not sure why this is (but it's in HYPE)
1650 else:
1651 temp_func = 1
1653 # Calculate uptake
1654 uptake_par = (self.uptake1 - self.uptake2) * exp(
1655 -self.uptake3 * days_after_sow
1656 )
1657 if (uptake_par + self.uptake2) > 0:
1658 N_common_uptake = (
1659 self.uptake1
1660 * self.uptake2
1661 * self.uptake3
1662 * uptake_par
1663 / ((self.uptake2 + uptake_par) ** 2)
1664 )
1665 N_common_uptake *= temp_func * constants.G_M2_TO_KG_M2 * self.area # [kg]
1666 P_common_uptake = N_common_uptake * self.uptake_PNratio
1667 # calculate maximum available uptake
1668 N_maximum_available_uptake = (
1669 max(0, self.storage["volume"] - self.wilting_point_m * self.area)
1670 / self.storage["volume"]
1671 * self.nutrient_pool.dissolved_inorganic_pool.storage["N"]
1672 )
1673 P_maximum_available_uptake = (
1674 max(0, self.storage["volume"] - self.wilting_point_m * self.area)
1675 / self.storage["volume"]
1676 * self.nutrient_pool.dissolved_inorganic_pool.storage["P"]
1677 )
1679 uptake = {
1680 "P": min(P_common_uptake, P_maximum_available_uptake),
1681 "N": min(N_common_uptake, N_maximum_available_uptake),
1682 }
1683 crop_uptake = self.nutrient_pool.dissolved_inorganic_pool.extract(uptake)
1684 out_ = self.empty_vqip()
1686 # Assuming plants eat N and P as nitrate and phosphate
1687 out_["nitrate"] = crop_uptake["N"]
1688 out_["phosphate"] = crop_uptake["P"]
1690 out2_ = self.pull_pollutants(out_)
1691 if not self.compare_vqip(out_, out2_):
1692 print("nutrient pool not tracking soil tank")
1694 return (self.empty_vqip(), out_)
1695 else:
1696 return (self.empty_vqip(), self.empty_vqip())
1698 def erosion(self):
1699 """Outflow function that erodes adsorbed/humus phosphorus and sediment and sends
1700 onwards to percolation/surface runoff/subsurface runoff.
1702 Returns:
1703 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1704 for mass balance checking.
1705 """
1706 # Parameters/equations from HYPE documentation (which explains why my
1707 # documentation is a bit ambiguous - because theirs is too)
1709 # Convert precipitation to MM since all the equations assume that
1710 precipitation_depth = self.get_data_input("precipitation") * constants.M_TO_MM
1712 # Calculate how much rain is mobilising erosion
1713 if precipitation_depth > 5:
1714 rainfall_energy = 8.95 + 8.44 * log10(
1715 precipitation_depth
1716 * (
1717 0.257
1718 + sin(2 * constants.PI * ((self.parent.t.dayofyear - 70) / 365))
1719 * 0.09
1720 )
1721 * 2
1722 )
1723 rainfall_energy *= precipitation_depth
1724 mobilised_rain = rainfall_energy * (1 - self.crop_cover) * self.erodibility
1725 else:
1726 mobilised_rain = 0
1728 # Calculate if any infiltration is mobilising erosion
1729 if self.infiltration_excess["volume"] > 0:
1730 mobilised_flow = (
1731 self.infiltration_excess["volume"] / self.area * constants.M_TO_MM * 365
1732 ) ** self.sreroexp
1733 mobilised_flow *= (
1734 (1 - self.ground_cover)
1735 * (1 / (0.5 * self.cohesion))
1736 * sin(self.slope / 100)
1737 / 365
1738 )
1739 else:
1740 mobilised_flow = 0
1742 # Sum flows (not sure why surface runoff isn't included) TODO I'm pretty sure it
1743 # should be included here
1744 total_flows = (
1745 self.infiltration_excess["volume"]
1746 + self.subsurface_flow["volume"]
1747 + self.percolation["volume"]
1748 ) # m3/dt + self.tank_recharge['volume'] (guess not needed)
1750 # Convert to MM/M2
1751 erodingflow = total_flows / self.area * constants.M_TO_MM
1753 # Calculate eroded sediment
1754 transportfactor = min(1, (erodingflow / 4) ** 1.3)
1755 erodedsed = (
1756 1000 * (mobilised_flow + mobilised_rain) * transportfactor
1757 ) # [kg/km2]
1758 # TODO not sure what conversion this HYPE 1000 is referring to
1760 # soil erosion with adsorbed inorganic phosphorus and humus phosphorus (erodedP
1761 # as P in eroded sediments and effect of enrichment)
1762 if erodingflow > 4:
1763 enrichment = 1.5
1764 elif erodingflow > 0:
1765 enrichment = 4 - (4 - 1.5) * erodingflow / 4
1766 else:
1767 return (self.empty_vqip(), self.empty_vqip())
1769 # Get erodable phosphorus
1770 erodableP = (
1771 self.nutrient_pool.get_erodable_P() / self.area * constants.KG_M2_TO_KG_KM2
1772 )
1773 erodedP = (
1774 erodedsed
1775 * (
1776 erodableP
1777 / (
1778 self.rooting_depth
1779 * constants.M_TO_KM
1780 * self.bulk_density
1781 * constants.KG_M3_TO_KG_KM3
1782 )
1783 )
1784 * enrichment
1785 ) # [kg/km2]
1787 # Convert to kg
1788 erodedP *= self.area * constants.M2_TO_KM2 # [kg]
1789 erodedsed *= self.area * constants.M2_TO_KM2 # [kg]
1791 # Allocate to different flows
1792 surface_erodedP = (
1793 self.srfilt * self.infiltration_excess["volume"] / total_flows * erodedP
1794 ) # [kg]
1795 surface_erodedsed = (
1796 self.srfilt * self.infiltration_excess["volume"] / total_flows * erodedsed
1797 ) # [kg]
1799 subsurface_erodedP = (
1800 self.macrofilt * self.subsurface_flow["volume"] / total_flows * erodedP
1801 ) # [kg]
1802 subsurface_erodedsed = (
1803 self.macrofilt * self.subsurface_flow["volume"] / total_flows * erodedsed
1804 ) # [kg]
1806 percolation_erodedP = (
1807 self.macrofilt * self.percolation["volume"] / total_flows * erodedP
1808 ) # [kg]
1809 percolation_erodedsed = (
1810 self.macrofilt * self.percolation["volume"] / total_flows * erodedsed
1811 ) # [kg]
1813 # Track mass balance
1814 in_ = self.empty_vqip()
1816 # Total eroded phosphorus
1817 eff_erodedP = percolation_erodedP + surface_erodedP + subsurface_erodedP # [kg]
1818 if eff_erodedP > 0:
1819 # Update nutrient pool
1820 org_removed, inorg_removed = self.nutrient_pool.erode_P(eff_erodedP)
1821 total_removed = inorg_removed + org_removed
1823 if abs(total_removed - eff_erodedP) > constants.FLOAT_ACCURACY:
1824 print("weird nutrients")
1826 # scale flows to split between inorganic and organic eroded P
1827 self.infiltration_excess["org-phosphorus"] += (
1828 surface_erodedP * org_removed / eff_erodedP
1829 )
1830 self.subsurface_flow["org-phosphorus"] += (
1831 subsurface_erodedP * org_removed / eff_erodedP
1832 )
1833 self.percolation["org-phosphorus"] += (
1834 percolation_erodedP * org_removed / eff_erodedP
1835 )
1837 # TODO Leon reckons this is conceptually dodgy.. but i'm not sure where else
1838 # adsorbed inorganic phosphorus should go
1839 self.infiltration_excess["phosphate"] += (
1840 surface_erodedP * inorg_removed / eff_erodedP
1841 )
1842 self.subsurface_flow["phosphate"] += (
1843 subsurface_erodedP * inorg_removed / eff_erodedP
1844 )
1845 self.percolation["phosphate"] += (
1846 percolation_erodedP * inorg_removed / eff_erodedP
1847 )
1849 # Entering the model (no need to uptake surface tank because both adsorbed
1850 # inorganic pool and humus pool are solids and so no tracked in the soil
1851 # water tank)
1852 in_["phosphate"] = inorg_removed
1853 in_["org-phosphorus"] = org_removed
1854 else:
1855 pass
1857 # Track sediment as solids
1858 self.infiltration_excess["solids"] += surface_erodedsed
1859 self.subsurface_flow["solids"] += subsurface_erodedsed
1860 self.percolation["solids"] += percolation_erodedsed
1862 in_["solids"] = surface_erodedsed + subsurface_erodedsed + percolation_erodedsed
1864 return (in_, self.empty_vqip())
1866 def denitrification(self):
1867 """Outflow function that performs denitirication processes, updating nutrient
1868 pool and soil tank.
1870 Returns:
1871 (tuple): A tuple containing a VQIP amount for model inputs and outputs
1872 for mass balance checking.
1873 """
1874 # Parameters/equations from HYPE documentation TODO could more of this be moved
1875 # to NutrientPool Calculate soil moisture dependence of denitrification
1876 soil_moisture_content = self.get_smc()
1877 if soil_moisture_content > self.field_capacity_m:
1878 denitrifying_soil_moisture_dependence = 1
1879 elif soil_moisture_content / self.field_capacity_m > self.limpar:
1880 denitrifying_soil_moisture_dependence = (
1881 ((soil_moisture_content / self.field_capacity_m) - self.limpar)
1882 / (1 - self.limpar)
1883 ) ** self.exppar
1884 else:
1885 denitrifying_soil_moisture_dependence = 0
1886 return (self.empty_vqip(), self.empty_vqip())
1888 # Get dissolved inorg nitrogen as a concentration and calculate factor
1889 din_conc = (
1890 self.nutrient_pool.dissolved_inorganic_pool.storage["N"]
1891 / self.storage["volume"]
1892 ) # [kg/m3]
1893 din_conc *= constants.KG_M3_TO_MG_L
1894 half_saturation_concentration_dependence_factor = din_conc / (
1895 din_conc + self.hsatINs
1896 )
1898 # Calculate and extract dentrified nitrogen
1899 denitrified_N = (
1900 self.nutrient_pool.dissolved_inorganic_pool.storage["N"]
1901 * half_saturation_concentration_dependence_factor
1902 * denitrifying_soil_moisture_dependence
1903 * self.nutrient_pool.temperature_dependence_factor
1904 * self.denpar
1905 )
1906 denitrified_request = self.nutrient_pool.get_empty_nutrient()
1907 denitrified_request["N"] = denitrified_N
1908 denitrified_N = self.nutrient_pool.dissolved_inorganic_pool.extract(
1909 denitrified_request
1910 )
1912 # Leon reckons this should leave the model (though I think technically some
1913 # small amount goes to nitrite)
1914 out_ = self.empty_vqip()
1915 out_["nitrate"] = denitrified_N["N"]
1917 # Update tank
1918 out2_ = self.pull_pollutants(out_)
1919 if not self.compare_vqip(out_, out2_):
1920 print("nutrient pool not tracking soil tank")
1922 return (self.empty_vqip(), out_)
1924 def adsorption(self):
1925 """Outflow function that calculates phosphorus adsorption/desorptions and
1926 updates soil tank and nutrient pools.
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 TODO could this be moved to the
1933 # nutrient pool?
1935 # Initialise mass balance checking
1936 in_ = self.empty_vqip()
1937 out_ = self.empty_vqip()
1939 # Get total phosphorus in pool available for adsorption/desorption
1940 limit = self.adosorption_nr_limit
1941 ad_de_P_pool = (
1942 self.nutrient_pool.adsorbed_inorganic_pool.storage["P"]
1943 + self.nutrient_pool.dissolved_inorganic_pool.storage["P"]
1944 ) # [kg]
1945 ad_de_P_pool /= self.area * constants.M2_TO_KM2 # [kg/km2]
1946 if ad_de_P_pool == 0:
1947 return (self.empty_vqip(), self.empty_vqip())
1949 # Calculate coefficient and concentration of adsorbed phosphorus
1950 soil_moisture_content = (
1951 self.get_smc() * constants.M_TO_MM
1952 ) # [mm] (not sure why HYPE has this in mm but whatever)
1953 conc_sol = (
1954 self.nutrient_pool.adsorbed_inorganic_pool.storage["P"]
1955 * constants.KG_TO_MG
1956 / (self.bulk_density * self.rooting_depth * self.area)
1957 ) # [mg P/kg soil]
1958 coeff = self.kfr * self.bulk_density * self.rooting_depth # [mm]
1960 # calculate equilibrium concentration
1961 if conc_sol <= 0:
1962 # Not sure how this would happen
1963 print("Warning: soil partP <=0. Freundlich will give error, take shortcut.")
1964 xn_1 = ad_de_P_pool / (soil_moisture_content + coeff) # [mg/l]
1965 ad_P_equi_conc = self.kfr * xn_1 # [mg/ kg]
1966 else:
1967 # Newton-Raphson method
1968 x0 = exp(
1969 (log(conc_sol) - log(self.kfr)) / self.nfr
1970 ) # initial guess of equilibrium liquid concentration
1971 fxn = x0 * soil_moisture_content + coeff * (x0**self.nfr) - ad_de_P_pool
1972 xn = x0
1973 xn_1 = xn
1974 j = 0
1975 while (
1976 abs(fxn) > limit and j < self.adsorption_nr_maxiter
1977 ): # iteration to calculate equilibrium concentations
1978 fxn = (
1979 xn * soil_moisture_content + coeff * (xn**self.nfr) - ad_de_P_pool
1980 )
1981 fprimxn = soil_moisture_content + self.nfr * coeff * (
1982 xn ** (self.nfr - 1)
1983 )
1984 dx = fxn / fprimxn
1985 if abs(dx) < (0.000001 * xn):
1986 # From HYPE... not sure what it means
1987 break
1988 xn_1 = xn - dx
1989 if xn_1 <= 0:
1990 xn_1 = 1e-10
1991 xn = xn_1
1992 j += 1
1993 ad_P_equi_conc = self.kfr * (xn_1**self.nfr)
1994 # print(ad_P_equi_conc, conc_sol)
1996 # Calculate new pool and concentration, depends on the equilibrium concentration
1997 if abs(ad_P_equi_conc - conc_sol) > 1e-6:
1998 request = self.nutrient_pool.get_empty_nutrient()
2000 # TODO not sure about this if statement, surely it would be triggered every
2001 # time
2002 adsdes = (ad_P_equi_conc - conc_sol) * (
2003 1 - exp(-self.kadsdes)
2004 ) # kinetic adsorption/desorption
2005 request["P"] = (
2006 adsdes
2007 * self.bulk_density
2008 * self.rooting_depth
2009 * (self.area * constants.M2_TO_KM2)
2010 ) # [kg]
2011 if request["P"] > 0:
2012 # Adsorption
2013 adsorbed = self.nutrient_pool.dissolved_inorganic_pool.extract(request)
2014 if (adsorbed["P"] - request["P"]) > constants.FLOAT_ACCURACY:
2015 print("Warning: freundlich flow adjusted, was larger than pool")
2016 self.nutrient_pool.adsorbed_inorganic_pool.receive(adsorbed)
2018 # Dissolved leaving the soil water tank and becoming solid
2019 out_["phosphate"] = adsorbed["P"]
2021 # Update tank
2022 out2_ = self.pull_pollutants(out_)
2023 if not self.compare_vqip(out_, out2_):
2024 print("nutrient pool not tracking soil tank")
2025 else:
2026 # Desorption
2027 request["P"] = -request["P"]
2028 desorbed = self.nutrient_pool.adsorbed_inorganic_pool.extract(request)
2029 if (desorbed["P"] - request["P"]) > constants.FLOAT_ACCURACY:
2030 print("Warning: freundlich flow adjusted, was larger than pool")
2031 self.nutrient_pool.dissolved_inorganic_pool.receive(desorbed)
2033 # Solid phosphorus becoming inorganic P in the soil water tank
2034 in_["phosphate"] = desorbed["P"]
2035 _ = self.push_storage(in_, force=True)
2037 return (in_, out_)
2039 def dry_deposition_to_tank(self, vqip):
2040 """Allocate dry deposition to surface tank, updating nutrient pool accordingly.
2042 Args:
2043 vqip (dict): A VQIP amount of dry deposition to send to tank
2045 Returns:
2046 vqip (dict): A VQIP amount of dry deposition that entered the tank (used
2047 for mass balance checking)
2048 """
2049 # Convert to nutrients
2050 deposition = self.nutrient_pool.get_empty_nutrient()
2051 deposition["N"] = vqip["nitrate"] + vqip["ammonia"]
2052 deposition["P"] = vqip["phosphate"]
2054 # Update nutrient pool
2055 in_ = self.nutrient_pool.allocate_dry_deposition(deposition)
2056 vqip = self.adjust_vqip_to_liquid(vqip, deposition, in_)
2058 # Update tank
2059 self.push_storage(vqip, force=True)
2060 return vqip
2062 def wet_deposition_to_tank(self, vqip):
2063 """Allocate wet deposition to surface tank, updating nutrient pool accordingly.
2065 Args:
2066 vqip (dict): A VQIP amount of dry deposition to send to tank
2068 Returns:
2069 vqip (dict): A VQIP amount of dry deposition that entered the tank (used
2070 for mass balance checking)
2071 """
2072 # Convert to nutrients
2073 deposition = self.nutrient_pool.get_empty_nutrient()
2074 deposition["N"] = vqip["nitrate"] + vqip["ammonia"]
2075 deposition["P"] = vqip["phosphate"]
2077 # Update nutrient pool
2078 in_ = self.nutrient_pool.allocate_wet_deposition(deposition)
2079 vqip = self.adjust_vqip_to_liquid(vqip, deposition, in_)
2081 # Update tank
2082 self.push_storage(vqip, force=True)
2083 return vqip
2086class IrrigationSurface(GrowingSurface):
2087 """"""
2089 def __init__(self, irrigation_coefficient=0.1, **kwargs):
2090 """A subclass of GrowingSurface that can calculate water demand for the crops
2091 that is not met by precipitation and use the parent node to acquire water. When
2092 the surface is created by the parent node, the irrigate function below is
2093 assigned.
2095 Args:
2096 irrigation_coefficient (float, optional): proportion area irrigated *
2097 proportion of demand met. Defaults to 0.1.
2098 """
2099 # Assign param
2100 self.irrigation_coefficient = irrigation_coefficient # proportion area
2101 # irrigated * proportion of demand met
2103 super().__init__(**kwargs)
2105 def irrigate(self):
2106 """Calculate water demand for crops and call parent node to acquire water,
2107 updating surface tank and nutrient pools."""
2108 if self.days_after_sow:
2109 # Irrigation is just difference between evaporation and precipitation amount
2110 irrigation_demand = (
2111 max(self.evaporation["volume"] - self.precipitation["volume"], 0)
2112 * self.irrigation_coefficient
2113 )
2114 if irrigation_demand > constants.FLOAT_ACCURACY:
2115 root_zone_depletion = self.get_cmd()
2116 if root_zone_depletion <= constants.FLOAT_ACCURACY:
2117 # TODO this isn't in FAO... but seems sensible
2118 irrigation_demand = 0
2120 # Pull water using parent node
2121 supplied = self.parent.pull_distributed(
2122 {"volume": irrigation_demand},
2123 of_type=["River", "Node", "Groundwater", "Reservoir"],
2124 )
2126 # update tank
2127 _ = self.push_storage(supplied, force=True)
2129 # update nutrient pools
2130 organic = {
2131 "N": supplied["org-nitrogen"],
2132 "P": supplied["org-phosphorus"],
2133 }
2134 inorganic = {
2135 "N": supplied["ammonia"] + supplied["nitrate"],
2136 "P": supplied["phosphate"],
2137 }
2138 self.nutrient_pool.allocate_organic_irrigation(organic)
2139 self.nutrient_pool.allocate_inorganic_irrigation(inorganic)
2142class GardenSurface(GrowingSurface):
2143 """"""
2145 # TODO - probably a simplier version of this is useful, building just on
2146 # pervioussurface
2147 def __init__(self, **kwargs):
2148 """A specific surface for gardens that treats the garden as a grass crop, but
2149 that can calculate/receive irrigation through functions that are assigned by the
2150 parent land node's handlers, which in turn are expected to be triggered by a
2151 query from an attached Demand node."""
2152 super().__init__(**kwargs)
2154 def calculate_irrigation_demand(self, ignore_vqip=None):
2155 """A check function (assigned by parent to push check from demand nodes) that
2156 calculations irrigation demand (i.e., difference between evaporation and
2157 preciptiation).
2159 Args:
2160 ignore_vqip (any, optional): Conventional push checks send an optional
2161 VQIP amount, however the intention with this check is to get the
2162 irrigation demand
2164 Returns:
2165 reply (dict): A VQIP amount of irrigation demand (note only 'volume' key
2166 is used)
2167 """
2168 # Calculate irrigation demand
2169 irrigation_demand = max(
2170 self.evaporation["volume"] - self.precipitation["volume"], 0
2171 )
2173 root_zone_depletion = self.get_cmd()
2174 if root_zone_depletion <= constants.FLOAT_ACCURACY:
2175 # TODO this isn't in FAO... but seems sensible
2176 irrigation_demand = 0
2178 # Reply as VQIP
2179 reply = self.empty_vqip()
2180 reply["volume"] = irrigation_demand
2181 return reply
2183 def receive_irrigation_demand(self, vqip):
2184 """A set function (assigned by parent to push set from demand nodes) that
2185 assigns irrigation water supply to the surface tank.
2187 Args:
2188 vqip (dict): A VQIP amount of irrigation to receive
2190 Returns:
2191 (dict): A VQIP amount of irrigation that was not received (should always
2192 be empty)
2193 """
2194 # update tank
2195 return self.push_storage(vqip, force=True)
2198class VariableAreaSurface(GrowingSurface):
2199 """"""
2201 def __init__(self, current_surface_area=0, **kwargs):
2202 """"""
2203 super().__init__(**kwargs)
2204 self.get_climate = self.get_climate_
2205 self.current_surface_area = current_surface_area
2207 def get_climate_(self):
2208 """
2210 Returns:
2212 """
2213 precipitation_depth = self.get_data_input("precipitation")
2214 evaporation_depth = self.get_data_input("et0") * self.et0_coefficient
2216 precipitation_depth *= self.current_surface_area / self.area
2217 evaporation_depth *= self.current_surface_area / self.area
2219 return precipitation_depth, evaporation_depth