diff --git a/docs/source/refs.bib b/docs/source/refs.bib index a33ec3607..58a9aa600 100644 --- a/docs/source/refs.bib +++ b/docs/source/refs.bib @@ -1,4 +1,89 @@ +@article{hodnett_marked_2002, + title = {Marked differences between van Genuchten soil water-retention parameters for temperate and tropical soils: a new water-retention pedo-transfer functions developed for tropical soils}, + volume = {108}, + issn = {0016-7061}, + url = {https://www.sciencedirect.com/science/article/pii/S0016706102001052}, + doi = {https://doi.org/10.1016/S0016-7061(02)00105-2}, + pages = {155--180}, + number = {3}, + journal = {Geoderma}, + author = {Hodnett, M. G. and Tomasella, J.}, + year = {2002}, + keywords = {Andosols, Ferralsols, Pedo-transfer function, Tropical soils, van Genuchten parameters, Water-retention curves}, +} + +@article{su_aerodynamic_2021, + title = {Aerodynamic resistance and Bowen ratio explain the biophysical effects of forest cover on understory air and soil temperatures at the global scale}, + volume = {308-309}, + issn = {0168-1923}, + url = {https://www.sciencedirect.com/science/article/pii/S0168192321003014}, + doi = {https://doi.org/10.1016/j.agrformet.2021.108615}, + pages = {108615}, + journal = {Agricultural and Forest Meteorology}, + author = {Su, Yongxian and Zhang, Chaoqun and Chen, Xiuzhi and Liu, Liyang and Ciais, Philippe and Peng, Jian and Wu, Shengbiao and Wu, Jianping and Shang, Jiali and Wang, Yingping and Yuan, Wenping and Yang, Yuanzhi and Wu, Zhifeng and Lafortezza, Raffaele}, + year = {2021} +} + +@article{saldarriaga_solar_1991, + title={Solar energy conversion efficiencies during succession of a tropical rain forest in Amazonia}, + volume={7}, + DOI={10.1017/S0266467400005393}, + number={2}, + journal={Journal of Tropical Ecology}, + author={Saldarriaga, Juan Guillermo and Luxmoore, Robert John}, + year={1991}, + pages={233-242} + } + +@article{tao_simplified_2021, + author = {Tao, Gaoliang and Wu, Zhijia and Li, Wentao and Li, Yi and Dong, Heming}, + title = {Simplified Relation Model of Soil Saturation Permeability Coefficient and Air-Entry Value and Its Application}, + journal = {Fractal and Fractional}, + volume = {5}, + year = {2021}, + number = {4}, + article-number = {180}, + url = {https://www.mdpi.com/2504-3110/5/4/180}, + issn = {2504-3110}, + doi = {10.3390/fractalfract5040180} +} + + +@misc{brink_modelling_2009, + month = {October}, + year = {2009}, + title = {Modelling the discharge of the Cidanau River in West Java with the HBV model}, + author = {F. van den {Brink}}, + url = {http://essay.utwente.nl/68862/} +} + +@article{reichardt_hydraulic_1993, + title = {Hydraulic variability in space and time in a dark red latosol of the tropics}, + volume = {60}, + issn = {0016-7061}, + url = {https://www.sciencedirect.com/science/article/pii/001670619390024F}, + doi = {https://doi.org/10.1016/0016-7061(93)90024-F}, + pages = {159--168}, + number = {1}, + journal = {Geoderma}, + author = {Reichardt, Klaus and Bacchi, Osny O. S. and Villagra, Maria de las Mercedes and Turatti, Ariovaldo L. and Pedrosa, Zildene O.}, + year = {1993} +} + +@article{gupta_global_2022, + title = {Global Soil Hydraulic Properties dataset based on legacy site observations and robust parameterization}, + volume = {9}, + issn = {2052-4463}, + url = {https://doi.org/10.1038/s41597-022-01481-5}, + doi = {10.1038/s41597-022-01481-5}, + pages = {444}, + number = {1}, + journal = {Scientific Data}, + author = {Gupta, Surya and Papritz, Andreas and Lehmann, Peter and Hengl, Tomislav and Bonetti, Sara and Or, Dani}, + year = {2022} +} + @article{geary_guide_2020, author = {Geary, William and Bode, Michael and Doherty, Tim and Fulton, Elizabeth and Nimmo, Dale and Tulloch, Ayesha and Tulloch, Vivitskaia and Ritchie, Euan}, year = {2020}, diff --git a/tests/models/hydrology/test_above_ground.py b/tests/models/hydrology/test_above_ground.py index 0d16bc656..5be20d7b4 100644 --- a/tests/models/hydrology/test_above_ground.py +++ b/tests/models/hydrology/test_above_ground.py @@ -126,7 +126,7 @@ def test_calculate_soil_evaporation(dens_air, latvap): pyrealm_const=PyrealmConst, ) - exp_evap = np.array([2.466861, 0.612504, 0.110356]) + exp_evap = np.array([2.18791, 0.521941, 0.090352]) np.testing.assert_allclose(result["soil_evaporation"], exp_evap, rtol=0.01) exp_ra = np.array([5.0, 10.0, 50.0]) np.testing.assert_allclose( diff --git a/tests/models/hydrology/test_below_ground.py b/tests/models/hydrology/test_below_ground.py index feb0f72de..dd4efe9ab 100644 --- a/tests/models/hydrology/test_below_ground.py +++ b/tests/models/hydrology/test_below_ground.py @@ -92,7 +92,7 @@ def test_calculate_matric_potential(): ) constants = HydroConsts() - expected_potentials = np.repeat(-17.320508, 3) + expected_potentials = np.repeat(-68.197326, 3) actual_potentials = calculate_matric_potential( effective_saturation=np.repeat(0.5, 3), air_entry_potential_inverse=constants.air_entry_potential_inverse, @@ -123,9 +123,11 @@ def test_update_groundwater_storage(dummy_climate_data): reservoir_const_lower_groundwater=HydroConsts.reservoir_const_lower_groundwater, ) - exp_groundwater = np.array([[453, 457, 459, 459], [450.0, 450.0, 450.0, 450]]) - exp_upper_flow = np.array([22.65, 22.85, 22.95, 22.95]) - exp_lower_flow = np.array([22.5, 22.5, 22.5, 22.5]) - np.testing.assert_allclose(result["groundwater_storage"], exp_groundwater) - np.testing.assert_allclose(result["subsurface_flow"], exp_upper_flow) - np.testing.assert_allclose(result["baseflow"], exp_lower_flow) + exp_groundwat = np.array( + [[451.3, 455.3, 457.3, 457.3], [451.7, 451.7, 451.7, 451.7]] + ) + exp_upper_flow = np.array([22.565, 22.765, 22.865, 22.865]) + exp_lower_flow = np.array([22.585, 22.585, 22.585, 22.585]) + np.testing.assert_allclose(result["groundwater_storage"], exp_groundwat, rtol=1e-05) + np.testing.assert_allclose(result["subsurface_flow"], exp_upper_flow, rtol=1e-05) + np.testing.assert_allclose(result["baseflow"], exp_lower_flow, rtol=1e-5) diff --git a/tests/models/hydrology/test_hydrology_model.py b/tests/models/hydrology/test_hydrology_model.py index 21c48a6dc..5c961c6e2 100644 --- a/tests/models/hydrology/test_hydrology_model.py +++ b/tests/models/hydrology/test_hydrology_model.py @@ -303,12 +303,12 @@ def test_setup( # Test 2d variables expected_2d = { "soil_moisture": [ - [246.698563, 246.698536, 246.698563, 246.698563], - [218.996509, 218.996509, 218.996509, 218.996509], + [247.25352, 247.25352, 247.25352, 247.25352], + [218.999041, 218.999041, 218.999041, 218.999041], ], "matric_potential": [ - [-76.910901, -76.910929, -76.910901, -76.910901], - [-105.707405, -105.707405, -105.707405, -105.707405], + [-66.550207, -66.550207, -66.550207, -66.550207], + [-217.596714, -217.596714, -217.596714, -217.596714], ], } @@ -325,11 +325,11 @@ def test_setup( # Test one dimensional variables expected_1d = { - "vertical_flow": [0.012935, 0.012935, 0.012935, 0.012935], + "vertical_flow": [6.916e-05, 6.916e-05, 6.916e-05, 6.916e-05], "total_river_discharge": [0, 0, 69035, 22660], - "surface_runoff": [106.654852, 106.654852, 106.654852, 106.654852], - "surface_runoff_accumulated": [0, 0, 9150, 2730], - "soil_evaporation": [7.172842, 7.172842, 7.172842, 7.172842], + "surface_runoff": [122.432417, 122.432417, 122.432417, 122.432417], + "surface_runoff_accumulated": [0, 0, 10530, 3300], + "soil_evaporation": [5.870856, 5.870856, 5.870856, 5.870856], } for var_name, expected_vals in expected_1d.items(): diff --git a/virtual_ecosystem/models/hydrology/constants.py b/virtual_ecosystem/models/hydrology/constants.py index c11a6719f..b90926aee 100644 --- a/virtual_ecosystem/models/hydrology/constants.py +++ b/virtual_ecosystem/models/hydrology/constants.py @@ -3,9 +3,8 @@ :mod:`~virtual_ecosystem.models.hydrology.hydrology_model`. These parameters are constants in that they should not be changed during a particular simulation. -TODO Soil parameters vary strongly with soil type and will require literature search and -sensitivity analysis to produce meaningful results. The current default values are just -examples within reasonable bounds. +Note that soil parameters vary strongly with soil type and can change over time. The +current default values are average best estimates within reasonable bounds. """ # noqa: D205 from dataclasses import dataclass @@ -17,73 +16,84 @@ class HydroConsts(ConstantsDataclass): """Dataclass to store all constants for the `hydrology` model.""" - soil_moisture_residual: float = 0.1 + soil_moisture_residual: float = 0.175 """Residual soil moisture, unitless. The residual soil moisture refers to the water that remains in the soil after prolonged drainage due to the force of gravity. It is the water that is tightly held by soil particles and is not easily available for plant roots to extract. The value is soil specific, the format here is volumentic relative water content (unitless, - between 0 and 1). + between 0 and 1). Average value of different soil textures across tropical regions + :cite:p:`hodnett_marked_2002`. """ - soil_moisture_saturation: float = 0.6 - """Soil moisture saturation, [%]. + soil_moisture_saturation: float = 0.51 + """Soil moisture saturation, unitless. Maximum amount of water a soil can hold when all its pores are completely filled with water — that is, the soil is fully saturated and contains no air in the pore - spaces. + spaces. Average value of different soil textures across tropical regions + :cite:p:`hodnett_marked_2002` + . """ - saturated_hydraulic_conductivity: float = 0.001 + saturated_hydraulic_conductivity: float = 3.5e-5 """Saturated hydraulic conductivity, [m s-1]. The saturated hydraulic conductivity is the measure of a soil's ability to transmit water through its pores. More specifically, is defined as the volumetric flow rate of water passing through a unit cross-sectional area of soil under a unit hydraulic - gradient (pressure difference). + gradient (pressure difference). Value for average tropical rainforest taken from + :cite:t:`gupta_global_2022`. """ - hydraulic_gradient: float = 0.01 - """The hydraulic gradient, [m m-1]. + hydraulic_gradient: float = 1.31 + """The hydraulic gradient, [m]. The hydraulic gradient is a measure of the change in hydraulic head (pressure) per unit of distance in a particular direction within a fluid or porous medium, such as soil or an aquifer. It represents the driving force behind the - movement of water and indicates the direction in which water will flow. + movement of water and indicates the direction in which water will flow. Value for + subtropical regions, taken from :cite:t:`reichardt_hydraulic_1993`; depends on the + soil type, permeability, slope, and water table depth. """ - van_genuchten_nonlinearily_parameter: float = 2.0 + van_genuchten_nonlinearily_parameter: float = 1.598 """Nonlinearity parameter n (dimensionless) in Mualem-van Genuchten model. This parameter is a fitting shape parameters of soil water retention curve, see - :cite:p:`van_genuchten_closed-form_1980`.""" + :cite:p:`van_genuchten_closed-form_1980`. Average value of different soil textures + across tropical regions + is taken from :cite:t:`hodnett_marked_2002`. + """ stream_flow_capacity: float = 5000.0 - """Stream flow capacity, [mm per timestep]. + """Stream flow capacity, [mm per day]. This parameter represents the maximum capacity of an average stream in the model. - At the moment, this is set as an arbitrary value, but could be used in the future to - flag flood events.""" + At the moment, this is set as an arbitrary value; we are working on getting a best + estimate.""" intercept_parameters: tuple[float, float, float] = (0.935, 0.498, 0.00575) - """Interception parameters. + """Interception parameters, unitless. Parameters in equation that estimates maximum canopy interception capacity after :cite:t:`von_hoyningen-huene_interzeption_1981`.""" veg_density_param: float = 0.046 - """Parameter to estimate vegetation density for maximum canopy interception. + """Parameter to estimate vegetation density for maximum interception, unitless. This parameter is used to estimate the water holding capacity of a canopy after - :cite:t:`von_hoyningen-huene_interzeption_1981`.""" + :cite:t:`von_hoyningen-huene_interzeption_1981`. The value is taken from + :cite:t:`van_der_knijff_lisflood_2010`.""" groundwater_capacity: float = 500 """Ground water storage capacity, [mm]. This parameter indicates how much water can be stored in the ground water reservoir which affects the vertical flow of water and the horizontal sub-surface flow. This - parameter is currently set to an arbitrary value and might.""" + parameter is currently set to an arbitrary value; we are working on getting a best + estimate.""" bypass_flow_coefficient: float = 1.0 """Empirical bypass flow coefficient, unitless. @@ -93,38 +103,30 @@ class HydroConsts(ConstantsDataclass): 0 means all surface water goes directly to groundwater, a value of 1 gives a linear relation between soil moisture and bypass flow.""" - air_entry_water_potential: float = -3.815 + air_entry_water_potential: float = -3.648 """Water potential at which soil pores begin to aerate, [kPa]. The constant is used to estimate soil water potential from soil moisture. As this estimation is a stopgap this constant probably shouldn't become a core constant. The - value is the average across all soil types found in - :cite:t:`cosby_statistical_1984`. In future, this could be calculated based on soil - texture. - """ - - campbell_pore_size_distribution: float = -7.22 - """Curvature of the water retention curve as indicator of pore size distribution. - - This constant is used to convert soil moisture to matric potential following - :cite:t:`campbell_simple_1974`. The value is the average across all soil types found - in :cite:t:`cosby_statistical_1984`; see documentation for - :attr:`air_entry_water_potential` for further details. + value is the average across soil types found in :cite:t:`tao_simplified_2021`. """ - extinction_coefficient_global_radiation: float = 0.7 + extinction_coefficient_global_radiation: float = 0.74 """Extinction coefficient for global radiation, [unitless]. This constant is used to reduce potential evaporation for bare soil to maximum shaded evaporation in :func:~virtual_ecosystem.models.hydrology.above_ground.calculate_soil_evaporation`. Typical values are 0.4 to 0.7 for monocotyledons and 0.65 to 1.1 for broad leaved - dicotyledons :cite:t:`monteith_light_1969`. The extinction coefficient can be + dicotyledons :cite:t:`monteith_light_1969`. The value for tropical forest is taken + from :cite:t:`saldarriaga_solar_1991`. The extinction coefficient can be estimated from measurements of PAR above and below a canopy with a known LAI. """ - max_percolation_rate_uzlz: float = 1 - """Maximum percolation rate between upper and lower groundwater zone, [mm d-1]""" + max_percolation_rate_uzlz: float = 2.7 + """Maximum percolation rate between upper and lower groundwater zone, [mm d-1]. + + Values for tropical rainforest are taken from :cite:t:`brink_modelling_2009`.""" groundwater_loss: float = 1 """Constant ground water loss, [mm]. @@ -134,16 +136,40 @@ class HydroConsts(ConstantsDataclass): """ reservoir_const_upper_groundwater: float = 20 - """Reservoir constant for the upper groundwater layer, [days]""" + """Reservoir constant for the upper groundwater layer, [days]. + + This parameter defines the residence time (in days) of water in the upper + groundwater zone before contributing to streamflow, with typical values for tropical + catchments ranging from 5 to 30 days depending on soil permeability and slope. + This parameter is currently set to an arbitrary value; we are working on getting a + better estimate.""" reservoir_const_lower_groundwater: float = 20 - """Reservoir constant for the lower groundwater layer, [days]""" + """Reservoir constant for the lower groundwater layer, [days]. + + This reservoir constant, measured in days, determines the residence time of water in + the lower groundwater zone. It influences how quickly water exits the lower zone as + baseflow. Typical values range from 10 to 5000 depending on catchment + characteristics. This parameter is currently set to an arbitrary value; we are + working on getting a better estimate.""" initial_aerodynamic_resistance_surface: float = 12.5 - """Initial aerodynamic resistance at the soil surface, [s m-1].""" + """Initial aerodynamic resistance at the soil surface, [s m-1]. + + This parameter is an initial estimate of the resistance to the transfer of momentum, + heat, and water vapour between the soil surface and the atmosphere. The value is + based on Australian evergreen forest, taken from :cite:t:`su_aerodynamic_2021`; + note that this assumes a dense canopy. + """ - initial_aerodynamic_resistance_canopy: float = 12.5 - """Initial aerodynamic resistance of the canopy, [s m-1].""" + initial_aerodynamic_resistance_canopy: float = 12.1 + """Initial aerodynamic resistance of the canopy, [s m-1]. + + This parameter is an initial estimate of the resistance to the transfer of momentum, + heat, and water vapour between the leaf surface and the atmosphere. The value is + based on Australian evergreen forest, taken from :cite:t:`su_aerodynamic_2021`; + note that this assumes a dense canopy. + """ drag_coefficient_evaporation: float = 0.2 """Drag coefficient for evaporation, dimensionless. @@ -152,22 +178,33 @@ class HydroConsts(ConstantsDataclass): the atmosphere.""" intercept_residence_time: float = 86400.0 - """Intecept residence time. + """Intecept residence time, [s]. The amount of time that water sits on the leaves before it evaporates or falls to - the ground.""" + the ground. We currently assume that at the end of each day, all water has either + evaporated or fallen to the ground.""" initial_stomatal_conductance: float = 1000.0 - """Initial stomatal conductance, [mmol m-2 s-1]""" + """Initial stomatal conductance, [mmol m-2 s-1]. + + Initial estimate of the rate at which water vapor and carbon dioxide pass + through the stomata of plant leaves, reflecting how open the stomata are and + regulating both transpiration and gas exchange. + """ pore_connectivity_parameter: float = 0.5 - """Pore connectivity parameter. + """Pore connectivity parameter, dimensionless. Dimensionless parameter used in van Genuchten-Mualem model to calculate unsaturated hydraulic conductivity.""" - air_entry_potential_inverse: float = 0.1 - """Inverse of air entry potential (parameter alpha in van Genuchten), [m-1].""" + air_entry_potential_inverse: float = 0.042 + """Inverse of air entry potential (parameter alpha in van Genuchten), [m-1]. + + The inverse of air entry potential describes how quickly soil water retention + decreases with increasing soil suction, with higher values indicating coarser soils + that drain more readily. Average value of different soil textures across tropical + regions :cite:p:`hodnett_marked_2002`.""" m_to_kpa: float = 9.804 """Factor to convert matric potential from m to kPa."""