LandscapeDNDC 1.37.0
Loading...
Searching...
No Matches
Public Member Functions | Static Public Attributes | Private Member Functions | Private Attributes | List of all members
ldndc::WatercycleDNDC Class Reference

Watercycle model WatercycleDNDC. More...

#include <models/watercycle/dndc/watercycle-dndc.h>

Public Member Functions

lerr_t solve ()
 Kicks off computation for one time step.
 

Static Public Attributes

static const unsigned int IMAX_W = 24
 

Private Member Functions

lerr_t event_irrigate ()
 considers water input due to irrigation and/or liquid manure evetns
 
lerr_t event_flood ()
 sets hydrologic conditions during flooding events, e.g.,
 
double balance_check (double *)
 checks each time step for correct water balance
 
void CalcRaining (unsigned int)
 apportion daily rainfall to subdaily rainfall
 
double rainfall_intensity ()
 
void CalcIrrigation (unsigned int)
 apportion daily irrigation to subdaily irrigation
 
lerr_t CalcIceContent ()
 Call of ice-content related processes.
 
lerr_t CalcPotEvapoTranspiration ()
 Potential evaporation.
 
void CalcInterception ()
 Calculates water quantity that is retained on leaf and stem bark surfaces.
 
double get_leaf_water ()
 Collects leaf water from all canopy layers.
 
double get_impedance_factor (size_t)
 Reduces water fluxes out of a layer if ice lenses have formed.
 
lerr_t CalcTranspiration ()
 
lerr_t CalcTranspirationCouvreur (double const &hr_potTransinm)
 
void CalcSnowEvaporation ()
 
void CalcSurfaceWaterEvaporation ()
 
void CalcSurfaceFlux ()
 
void WatercycleDNDC_preferential_flow ()
 
void CalcSoilEvapoPercolation ()
 
void SoilWaterEvaporation (size_t)
 
double WaterFlowCapillaryRise (size_t, double const &)
 Similar to WaterFlowBelowFieldCapacity(...) but restricts water flow by water availability from soil layer below.
 
double WaterFlowAboveFieldCapacity (size_t, double const &)
 
double get_fsl (size_t)
 Returns factor accounting for soil layer depth. Originally, a constant soil layer depth of 0.02 m has been assumed. Due to flexible layer depth, an additional factor has been introduced that accounts for a higher resistance of larger layers.
 
double get_time_step_duration ()
 Time fraction with regard to daily rates.
 
double get_clay_limitation (size_t _sl, double _water_avail, double _water_ref)
 Restrains water availability due to clay content.
 

Private Attributes

double thornthwaite_heat_index
 heat index for potential evaporation using approach after Thornthwaite
 

Detailed Description

Watercycle model WatercycleDNDC.

Member Function Documentation

◆ balance_check()

double ldndc::WatercycleDNDC::balance_check ( double * _balance)
private

checks each time step for correct water balance

Balance check of water.

2660{
2661 double balance = get_leaf_water() + wc_.surface_ice + wc_.surface_water;
2662 for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2663 {
2664 balance += (wc_.wc_sl[sl] + wc_.ice_sl[sl]) * sc_.h_sl[sl];
2665 }
2666 balance += ( - m_wc.accumulated_irrigation_event_executed
2667 - wc_.accumulated_irrigation_automatic
2668 - wc_.accumulated_irrigation_ggcmi
2669 - wc_.accumulated_irrigation_reservoir_withdrawal
2670 - wc_.accumulated_precipitation
2671 - wc_.accumulated_groundwater_access
2672 + wc_.accumulated_percolation
2673 + wc_.accumulated_runoff
2674 + wc_.accumulated_wateruptake_sl.sum()
2675 + wc_.accumulated_interceptionevaporation
2676 + wc_.accumulated_soilevaporation
2677 + wc_.accumulated_surfacewaterevaporation);
2678
2679 if ( _balance)
2680 {
2681 double const balance_delta( std::abs( *_balance - balance));
2682 if ( balance_delta > cbm::DIFFMAX)
2683 {
2684 KLOGWARN( "Water leakage: difference is ", balance - *_balance);
2685 }
2686 }
2687
2688 for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2689 {
2690 if ( cbm::flt_less( wc_.wc_sl[sl], 0.0))
2691 {
2692 KLOGWARN( "Negative water content of: ", wc_.wc_sl[sl], " in layer: ", sl, ". Water content set to zero.");
2693 wc_.wc_sl[sl] = 0.0;
2694 }
2695 if ( cbm::flt_less( wc_.ice_sl[sl], 0.0))
2696 {
2697 KLOGWARN( "Negative ice content of: ", wc_.ice_sl[sl], " in layer: ", sl, ". Ice content set to zero.");
2698 wc_.ice_sl[sl] = 0.0;
2699 }
2700 }
2701
2702 return balance;
2703}
ldndc::WatercycleDNDC::get_leaf_water
double get_leaf_water()
Collects leaf water from all canopy layers.
Definition watercycle-dndc.cpp:1110

References balance_check(), and get_leaf_water().

Referenced by balance_check(), and solve().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcInterception()

void ldndc::WatercycleDNDC::CalcInterception ( )
private

Calculates water quantity that is retained on leaf and stem bark surfaces.

  • Rainfall and interception are assumed to happen at the end of the timestep and evaporation is realised at once at the start of the timestep.
  • The model assumes a homogeneous horizontal distribution of LAI although a non-linear dependency to crown cover (relative foliage clumping) has been reported (e.g. Teklehaimanot et al. 1991).
  • The default water-interception for foliage almost equal to the largest species specific value

References get_impedance_factor(), and get_leaf_water().

Referenced by solve().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcIrrigation()

void ldndc::WatercycleDNDC::CalcIrrigation ( unsigned int _i)
private

apportion daily irrigation to subdaily irrigation

Automatic irrigation to prevent water stress

GGCMI style irrigation

  • fill up the water content of each soil layer to field capacity at the start of each day
  • water supplied just injected into soil layer
  • no interception, no run off, no surface-water
  • only irrigate if a crop is present

Irrigation due to flooding

irrigation triggered by flooding management

Irrigation due to explicit irrigation event

805{
806 //reset irrigation
807 m_wc.irrigation = 0.0;
808
812 if (cbm::flt_in_range_lu(0.0, m_cfg.automaticirrigation, 1.0))
813 {
814 if (m_veg->size() > 0) //check crop is present
815 {
816 if (cbm::flt_greater(mc_.nd_airtemperature,10.0)) // only add water if temperature > 10C)
817 {
818 for (size_t l = 0; l < m_soillayers->soil_layer_cnt(); ++l)
819 {
820 if (cbm::flt_greater_zero(m_veg->mfrt_sl(l))) //only count water need for soil layers that contain roots
821 {
822 // trigger irrigation events as soon as drought stress appears
823 double const wc_target(wc_.wc_wp_sl[l] + m_cfg.automaticirrigation * (wc_.wc_fc_sl[l] - wc_.wc_wp_sl[l]));
824 if (cbm::flt_less(wc_.wc_sl[l], wc_target) &&
825 cbm::flt_equal_zero(wc_.ice_sl[l])) //only add water if no ice is present
826 {
827 // irrigation assumed until target water content
828 double const irrigate_layer((wc_target - wc_.wc_sl[l]) * sc_.h_sl[l]);
829 m_wc.irrigation += irrigate_layer;
830 wc_.accumulated_irrigation_automatic += irrigate_layer;
831 }
832 }
833 }
834
835 }
836 else
837 {
838 }
839
840 }
841 }
842
843 if ( cbm::flt_in_range_lu( 0.0, m_eventflood.get_saturation_level(), 1.0) &&
844 cbm::is_valid( m_eventflood.get_irrigation_height()) &&
845 cbm::is_valid( m_eventflood.get_soil_depth()))
846 {
847 double wc_sum( 0.0);
848 double wc_fc_sum( 0.0);
849 double wc_wp_sum( 0.0);
850 for ( size_t l = 0; l < m_soillayers->soil_layer_cnt(); ++l)
851 {
852 wc_sum += wc_.wc_sl[l] * sc_.h_sl[l];
853 wc_fc_sum += wc_.wc_fc_sl[l] * sc_.h_sl[l];
854 wc_wp_sum += wc_.wc_wp_sl[l] * sc_.h_sl[l];
855
856 if ( cbm::flt_greater_equal( sc_.depth_sl[l], m_eventflood.get_soil_depth()))
857 {
858 break;
859 }
860 }
861
862 // trigger irrigation events as soon as drought stress appears
863 double const wc_target( wc_wp_sum + m_eventflood.get_saturation_level() * (wc_fc_sum - wc_wp_sum));
864 if ( cbm::flt_less( wc_sum, wc_target))
865 {
866 m_wc.irrigation += m_eventflood.get_irrigation_height();
867 wc_.accumulated_irrigation_automatic += m_eventflood.get_irrigation_height();
868 }
869 }
870
878 if ( m_cfg.ggcmi_irrigation)
879 {
880 if( (lclock()->subday() == 1) && (_i == 0))
881 {
882 if( m_veg->size() > 0 ) //check crop is present
883 {
884 double water_supply( 0.0);
885 for ( size_t l = 0; l < m_soillayers->soil_layer_cnt(); ++l)
886 {
887 if ( cbm::flt_greater_zero( m_veg->mfrt_sl( l))) //only irrigate soil layers that contain roots
888 {
889 if( cbm::flt_less( wc_.wc_sl[l], wc_.wc_fc_sl[l]) &&
890 cbm::flt_equal_zero( wc_.ice_sl[l])) //only add water if no ice is present
891 {
892 water_supply += (wc_.wc_fc_sl[l] - wc_.wc_sl[l]) * sc_.h_sl[l]; //amount of water added (m)
893 wc_.wc_sl[l] = wc_.wc_fc_sl[l];
894 }
895 }
896 }
897 wc_.accumulated_irrigation_ggcmi += water_supply;
898 }
899 }
900 }
901
905 if ( cbm::is_valid( m_eventflood.get_water_table()))
906 {
907 /* Dynamic flooding */
908 if ( cbm::is_valid( m_eventflood.get_irrigation_height()))
909 {
910 bool irrigate( false);
911
912 /* Irrigate after surface water table drop */
913 if ( cbm::flt_greater_equal_zero( m_eventflood.get_water_table()))
914 {
915 if ( cbm::flt_greater( m_eventflood.get_water_table(), wc_.surface_water))
916 {
917 irrigate = true;
918 }
919 }
920 else
921 {
922 if ( !cbm::flt_greater_zero( wc_.surface_water))
923 {
924 double water_tot_sum( 0.0);
925 double water_fc_sum( 0.0);
926 for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
927 {
928 water_tot_sum += wc_.wc_sl[sl] * sc_.h_sl[sl];
929 water_fc_sum += wc_.wc_fc_sl[sl] * sc_.h_sl[sl];
930 if ( cbm::flt_greater_equal( sc_.depth_sl[sl], -m_eventflood.get_water_table()))
931 {
932 break;
933 }
934 }
935 // Note: The tipping bucket approach may block percolation into deeper soil layers
936 // if the upper layers have higher water content. This can lead to water reduction
937 // from evapotranspiration in deeper layers without replenishment from above.
938 // Workaround: AWD is triggered when the integrated topsoil water drops below
939 // 90% of the integrated topsoil field capacity.
940 if ( cbm::flt_less( water_tot_sum, 0.9 * water_fc_sum))
941 {
942 irrigate = true;
943 }
944 }
945 }
946 // trigger AWD irrigation event
947 if ( irrigate)
948 {
949 double water_supply( cbm::bound_min( 0.0, m_eventflood.get_irrigation_height() - wc_.surface_water));
950 for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
951 {
952 // set soil fully water saturated
953 if ( cbm::flt_greater( sc_.poro_sl[sl], wc_.wc_sl[sl]))
954 {
955 water_supply += ((sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl]);
956 }
957 }
958
959 //remove rainfall (less irrigation water is supplied if it if raining)
960 water_supply = cbm::bound_min( 0.0, water_supply - m_mc.precipitation);
961
962 if ( m_eventflood.have_unlimited_water())
963 {
964 wc_.accumulated_irrigation_automatic += water_supply;
965 }
966 else
967 {
968 if ( cbm::flt_greater( wc_.irrigation_reservoir, water_supply))
969 {
970 wc_.irrigation_reservoir -= water_supply;
971 }
972 else
973 {
974 water_supply = wc_.irrigation_reservoir;
975 wc_.irrigation_reservoir = 0.0;
976 }
977 wc_.accumulated_irrigation_reservoir_withdrawal += water_supply;
978 }
979 m_wc.irrigation += water_supply;
980 }
981 }
982 /* Static flooding */
983 else
984 {
985 if ( !cbm::flt_greater_zero( m_eventflood.get_water_table()))
986 {
987 double water_supply( 0.0);
988 for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
989 {
990 if ( cbm::flt_greater( sc_.depth_sl[sl], -m_eventflood.get_water_table()) &&
991 cbm::flt_greater( sc_.poro_sl[sl], wc_.wc_sl[sl]))
992 {
993 double const scale( cbm::flt_less( sc_.depth_sl[sl] - sc_.h_sl[sl], -m_eventflood.get_water_table()) ?
994 (-m_eventflood.get_water_table() - (sc_.depth_sl[sl] - sc_.h_sl[sl])) / sc_.h_sl[sl] :
995 1.0);
996
997 water_supply += scale * (sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl];
998 }
999 }
1000
1001 //remove rainfall (less irrigation water is supplied if it is raining)
1002 water_supply = cbm::bound( 0.0, water_supply - m_mc.precipitation, wc_.sks_sl[0]);
1003
1004 /* Surface water will be actively drained */
1005 if ( cbm::flt_greater_equal( wc_.surface_water, water_supply))
1006 {
1007 wc_.accumulated_runoff += wc_.surface_water - water_supply;
1008 wc_.surface_water = water_supply;
1009 water_supply = 0.0;
1010 }
1011 else if ( cbm::flt_greater_zero( wc_.surface_water))
1012 {
1013 water_supply -= wc_.surface_water;
1014 }
1015
1016 if ( m_eventflood.have_unlimited_water())
1017 {
1018 wc_.accumulated_irrigation_automatic += water_supply;
1019 }
1020 else
1021 {
1022 if ( cbm::flt_greater( wc_.irrigation_reservoir, water_supply))
1023 {
1024 wc_.irrigation_reservoir -= water_supply;
1025 }
1026 else
1027 {
1028 water_supply = wc_.irrigation_reservoir;
1029 wc_.irrigation_reservoir = 0.0;
1030 }
1031 wc_.accumulated_irrigation_reservoir_withdrawal += water_supply;
1032 }
1033
1034 m_wc.irrigation += water_supply;
1035 }
1039 else if ( cbm::flt_greater( m_eventflood.get_water_table(), wc_.surface_water))
1040 {
1041 //adjust surface and soil water in case of flooding event
1042 //water required to fill up surface water
1043 double water_supply( m_eventflood.get_water_table() - wc_.surface_water);
1044
1045 //add water required to fill top x layers up to saturation
1046 unsigned long no_layers(3);
1047 if( no_layers > m_soillayers->soil_layer_cnt())
1048 {
1049 no_layers = m_soillayers->soil_layer_cnt();
1050 }
1051
1052 for ( size_t sl = 0; sl < no_layers; ++sl)
1053 {
1054 water_supply += ((sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl]);
1055 }
1056
1057 //remove rainfall (less irrigation water is supplied if it if raining)
1058 water_supply = cbm::bound_min( 0.0, water_supply - m_mc.precipitation);
1059
1060 if ( m_eventflood.have_unlimited_water())
1061 {
1062 wc_.accumulated_irrigation_automatic += water_supply;
1063 }
1064 else
1065 {
1066 if ( cbm::flt_greater( wc_.irrigation_reservoir, water_supply))
1067 {
1068 wc_.irrigation_reservoir -= water_supply;
1069 }
1070 else
1071 {
1072 water_supply = wc_.irrigation_reservoir;
1073 wc_.irrigation_reservoir = 0.0;
1074 }
1075 wc_.accumulated_irrigation_reservoir_withdrawal += water_supply;
1076 }
1077 m_wc.irrigation += water_supply;
1078 }
1079 }
1080 }
1081
1085 double const irrigation_scheduled( cbm::bound_min( 0.0, wc_.accumulated_irrigation - m_wc.accumulated_irrigation_event_executed));
1086
1087 if ( cbm::flt_greater( irrigation_scheduled, 24.0 * rainfall_intensity()))
1088 {
1089 m_wc.irrigation += irrigation_scheduled / 24.0;
1090 m_wc.accumulated_irrigation_event_executed += irrigation_scheduled / 24.0;
1091 }
1092 else if ( cbm::flt_greater( irrigation_scheduled, (double)hours_per_time_step() * rainfall_intensity()))
1093 {
1094 m_wc.irrigation += (double)hours_per_time_step() * rainfall_intensity();
1095 m_wc.accumulated_irrigation_event_executed += (double)hours_per_time_step() * rainfall_intensity();
1096 }
1097 else
1098 {
1099 m_wc.irrigation += irrigation_scheduled;
1100 m_wc.accumulated_irrigation_event_executed += irrigation_scheduled;
1101 }
1102}
ldndc::WatercycleDNDC::rainfall_intensity
double rainfall_intensity()
Definition watercycle-dndc.cpp:103

References CalcIrrigation(), and rainfall_intensity().

Referenced by CalcIrrigation(), and solve().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcPotEvapoTranspiration()

lerr_t ldndc::WatercycleDNDC::CalcPotEvapoTranspiration ( )
private

Potential evaporation.


there are two methods of calculating potential evaporation (choose via macro POTENTIAL_EVAPORATION above)

  • [0] Procedure according to Thornthwaite 1948 with modifications by Camargo et al. 1999 and Peireira and Pruitt 2004 (implemented from Pereira and Pruitt 2004) to better account for dry environments and daylength.
    NOTE:
    • monthly temperatures are assumed to be equal to daily average temperature of the middle of the month, calculated from continous equations equal to those used in the weather generator.
    • the model uses always daily average temperature as input even if it is run in sub-daily mode
    • the original PnET-N-DNDC code which contains modification that could not be found in the literature has been outcommented.
  • [1] uses the Priestley Taylor equation (Priestly and Taylor 1972) to calculate daily and hourly potential evapotranspiration (implemented from internet) vps calculation according to Allen et al. 1998, cit. in Cai et al. 2007
537{
538 // Daily evaportation demand is derived from the Thornthwaite, the Penman, or the Priestley-Taylor method (standard is 'thornthwaite).
539 day_potevapo = 0.0;
540 if ( m_cfg.potentialevapotranspiration_method == "thornthwaite")
541 {
542 if ( lclock()->is_position( TMODE_PRE_YEARLY))
543 {
544 thornthwaite_heat_index = ldndc::thornthwaite_heat_index( m_setup->latitude(),
545 lclock()->days_in_year(),
546 m_climate->annual_temperature_average(),
547 m_climate->annual_temperature_amplitude());
548 }
549
550 double const daylength(ldndc::meteo::daylength(m_setup->latitude(), lclock()->yearday()));
551 double const avg_dayl(12.0);
552
553 day_potevapo = m_sipar->WCDNDC_INCREASE_POT_EVAPOTRANS()
554 * ldndc::thornthwaite(
555 mc_.nd_airtemperature,
556 daylength, avg_dayl, lclock()->days_in_month(),
557 thornthwaite_heat_index);
558 }
559 else if ( (m_cfg.potentialevapotranspiration_method == "penman") ||
560 (m_cfg.potentialevapotranspiration_method == "priestleytaylor"))
561 {
562 cbm::sclock_t const & clk( lclock_ref());
563
564 /* vapour pressure [10^3:Pa] */
565 double vps( 0.0);
566 ldndc::meteo::vps( mc_.nd_airtemperature, &vps);
567 double vp( vps - m_climate->vpd_day( clk));
568
569 if( m_cfg.potentialevapotranspiration_method == "priestleytaylor")
570 {
571 day_potevapo = m_sipar->WCDNDC_INCREASE_POT_EVAPOTRANS()
572 * ldndc::priestleytaylor( clk.yearday(),
573 clk.days_in_year(),
574 mc_.albedo,
575 m_setup->latitude(),
576 mc_.nd_airtemperature,
577 mc_.nd_shortwaveradiation_in,
578 vp,
579 m_sipar->PT_ALPHA());
580 }
581 else if( m_cfg.potentialevapotranspiration_method == "penman")
582 {
583 /* leaf area index */
584 double const lai( m_veg->lai());
585 ldndc::surface_type surface( cbm::flt_greater_zero( lai) ? short_grass : bare_soil);
586
587 day_potevapo = m_sipar->WCDNDC_INCREASE_POT_EVAPOTRANS()
588 * ldndc::penman( surface,
589 clk.yearday(),
590 clk.days_in_year(),
591 mc_.albedo,
592 m_setup->latitude(),
593 mc_.nd_shortwaveradiation_in * cbm::SEC_IN_DAY,
594 mc_.nd_airtemperature,
595 mc_.nd_windspeed,
596 vp);
597 }
598 }
599 else
600 {
601 KLOGERROR("[BUG] huh : it seems that this potentialevapotranspiration method is not registered",
602 "tell the maintainers refering to this error message");
603 return LDNDC_ERR_FAIL;
604 }
605
606 // The daily evaporation demand is distributed throughout the day (standard is 'previouspattern)
607 // TODO: design a routine that distributes evaporation demand according to temperature or humidity
608 if ( m_cfg.evapotranspiration_method == "uniform")
609 {
610 // hourly evaporation demand is calculated from daily demand equally throughout the day
611 hr_pot_evapotrans = day_potevapo * nhr / cbm::HR_IN_DAY;
612 }
613 else if ( m_cfg.evapotranspiration_method == "previouspattern")
614 {
615 // hourly evaporation demand is derived from the weighted expression of previous values
616 hr_pot_evapotrans = day_potevapo * pot_evapotranspiration_ts[lclock()->subday()-1] / nwc;
617 }
618 else if (m_cfg.evapotranspiration_method == "daytimeonly")
619 {
620 // hourly evaporation demand is distributed throughout the sunlight hours only (similar as radiation)
621 double const hour(lclock()->subday());
622 double const day(lclock()->yearday());
623 double const hsun(0.833 * cbm::PI / 180.0);
624 double const lat(m_setup->latitude() * cbm::PI / 180.0);
625 double const dec(-23.44 * (cbm::PI / 180.0) * std::cos(2.0 * cbm::PI * (day - 172.0) / lclock()->days_in_year()));
626 double const help((std::sin(hsun) - std::sin(lat) * std::sin(dec)) / (std::cos(lat) * std::cos(dec)));
627 double dayl(0.0);
628 if (help < -0.999) dayl = 0.0;
629 else if (help > 0.999) dayl = 24.0;
630 else dayl = 24.0 - 24.0 / cbm::PI * std::acos((std::sin(hsun) - std::sin(lat) * std::sin(dec)) / (std::cos(lat) * std::cos(dec)));
631 if (dayl < 1.0) dayl = 0.0;
632 double const rise(12.0 - dayl * 0.5);
633
634 double factor(0.0);
635 if (dayl > 0.0) factor = (-std::cos(2.0 * cbm::PI * (hour - rise) / (cbm::HR_IN_DAY - 2.0 * rise)) + 1.0) / (cbm::HR_IN_DAY - 2.0 * rise);
636
637 double fday(0.0);
638 if (hour > rise && hour < (24.0 - rise)) fday = factor;
639
640 hr_pot_evapotrans = day_potevapo * fday;
641 }
642 else
643 {
644 KLOGERROR("[BUG] huh : it seems that this evapotranspiration method option is not registered ",
645 "tell the maintainers refering to this error message");
646 return LDNDC_ERR_FAIL;
647 }
648
649 // the evaporation limit for all timesteps is the daily demand
650 wc_.accumulated_potentialevaporation += hr_pot_evapotrans;
651
652 return LDNDC_ERR_OK;
653}
ldndc::WatercycleDNDC::thornthwaite_heat_index
double thornthwaite_heat_index
heat index for potential evaporation using approach after Thornthwaite
Definition watercycle-dndc.h:134
ldndc::thornthwaite_heat_index
double LDNDC_API thornthwaite_heat_index(double, size_t, double, double)
Thornthwaite heat index.
Definition ld_thornthwaite.cpp:100
ldndc::thornthwaite
double LDNDC_API thornthwaite(double, double, double, double, double)
Potential evapotranspiration after thornthwaite:1948a.
Definition ld_thornthwaite.cpp:41

References CalcPotEvapoTranspiration(), ldndc::thornthwaite(), ldndc::thornthwaite_heat_index(), and thornthwaite_heat_index.

Referenced by CalcPotEvapoTranspiration(), and solve().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcSnowEvaporation()

void ldndc::WatercycleDNDC::CalcSnowEvaporation ( )
private
Parameters
[in]None
[out]None
Returns
None
2023{
2024 //Average limit of snowEvaporation is 0.00013m day-1 (Granger and Male 1978)
2025 double hr_potESnow( cbm::bound(0.0, hr_pot_evapotrans, 0.00013 * nhr/24.0));
2026
2027 if ( cbm::flt_greater_zero( wc_.surface_ice) &&
2028 cbm::flt_greater_zero( hr_potESnow))
2029 {
2030 //If there is a snowpack and there is a potential to evaporate
2031 if ( wc_.surface_ice > hr_potESnow)
2032 {
2033 //Partial evaporation from the snow
2034 hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - hr_potESnow);
2035 wc_.accumulated_surfacewaterevaporation += hr_potESnow;
2036 wc_.surface_ice -= hr_potESnow;
2037 }
2038 else
2039 {
2040 //Total evaporation of the snowpack
2041 hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - wc_.surface_ice);;
2042 wc_.accumulated_surfacewaterevaporation += wc_.surface_ice;
2043 wc_.surface_ice = 0.0;
2044 }
2045 }
2046}

References CalcSnowEvaporation().

Referenced by CalcSnowEvaporation(), CalcTranspiration(), and solve().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcSoilEvapoPercolation()

void ldndc::WatercycleDNDC::CalcSoilEvapoPercolation ( )
private
Parameters
[in]None
[out]None
Returns
None
2160{
2161 // reset all water fluxes
2162 for (size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2163 {
2164 m_wc.percolation_sl[sl] = 0.0;
2165 }
2166
2167 // fill groundwater layer with water
2168 for (int sl = m_soillayers->soil_layer_cnt() - 1; sl >= 0; --sl)
2169 {
2170 if (cbm::flt_greater_equal(sc_.depth_sl[sl], ground_water_table))
2171 {
2172 // ice volume is excepted from groundwater
2173 double const ice_vol(wc_.ice_sl[sl] / cbm::DICE);
2174 if (cbm::flt_greater(sc_.poro_sl[sl], ice_vol))
2175 {
2176 // new water in the layer is generated from groundwater
2177 double const gw_access(cbm::bound_min(0.0, (sc_.poro_sl[sl] - ice_vol - wc_.wc_sl[sl]) * sc_.h_sl[sl] ));
2178 wc_.wc_sl[sl] += gw_access / sc_.h_sl[sl];
2179 wc_.accumulated_groundwater_access += gw_access;
2180 }
2181 }
2182 else
2183 {
2184 break;
2185 }
2186 }
2187
2188 for (size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2189 {
2190 // the first soil (litter) layer is filled with surface water up to its capacity
2191 if (sl == 0)
2192 {
2193 double delta_wc(0.0);
2194 double const wc_old(wc_.wc_sl[sl]);
2195 double const wc_cap( (sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl]); // uptake capacity in first soil layer [mm]
2196
2197 if (cbm::flt_greater(wc_.surface_water, wc_cap))
2198 {
2199 wc_.wc_sl[sl] = sc_.poro_sl[sl];
2200 delta_wc = cbm::bound_max(wc_.surface_water, (cbm::bound_min(0.0, wc_.wc_sl[sl] - wc_old) * sc_.h_sl[sl]));
2201 wc_.surface_water -= delta_wc;
2202 }
2203 else
2204 {
2205 wc_.wc_sl[sl] = (wc_.wc_sl[sl] * sc_.h_sl[sl] + wc_.surface_water) / sc_.h_sl[sl];
2206 delta_wc = cbm::bound_max(wc_.surface_water, (cbm::bound_min(0.0, wc_.wc_sl[sl] - wc_old) * sc_.h_sl[sl]));
2207 wc_.surface_water = 0.0;
2208 }
2209
2210 wc_.accumulated_infiltration += delta_wc;
2211 }
2212 // other soil layers get updated by percolation rates from the layer above
2213 else
2214 {
2215 wc_.wc_sl[sl] += m_wc.percolation_sl[sl - 1] / sc_.h_sl[sl];
2216 }
2217
2218 // percolation for all layers below groundwater
2219 if (cbm::flt_greater(sc_.depth_sl[sl], ground_water_table))
2220 {
2221
2222 // in case groundwater reaches the soil surface, the first outflux is estimated by saturated water flow
2223 if (sl == 0)
2224 {
2225 m_wc.percolation_sl[sl] = WaterFlowSaturated(sl, get_impedance_factor(sl));
2226 }
2227 // other soil layers are assumed to percolate by the same rate as the last non-groundwater layer
2228 else
2229 {
2230 m_wc.percolation_sl[sl] = m_wc.percolation_sl[sl - 1];
2231 }
2232
2233 }
2234 // in all other cases percolation
2235 else if (sl + 1 == m_soillayers->soil_layer_cnt())
2236 {
2237 m_wc.percolation_sl[sl] = 0.0;
2238 if (cbm::flt_greater_zero(wc_.sks_sl[sl]))
2239 {
2240 if (cbm::flt_greater(wc_.wc_sl[sl], sc_.poro_sl[sl]) &&
2241 cbm::flt_greater(wc_.wc_sl[sl], wc_.wc_fc_sl[sl])) //security, should be always true if first condition is met
2242 {
2243 //maximum allowed flow until field capacity
2244 double const max_flow((wc_.wc_sl[sl] - wc_.wc_fc_sl[sl]) * sc_.h_sl[sl]);
2245 m_wc.percolation_sl[sl] = cbm::bound_max(WaterFlowSaturated(sl, get_impedance_factor(sl)), max_flow);
2246 }
2247
2248 //water flux above field capacity
2249 else if (cbm::flt_greater_equal(wc_.wc_sl[sl] + wc_.ice_sl[sl], wc_.wc_fc_sl[sl]))
2250 {
2251 m_wc.percolation_sl[sl] = WaterFlowAboveFieldCapacity(sl, get_impedance_factor(sl));
2252 }
2253
2254 //water flux below field capacity and above wilting point
2255 else if (cbm::flt_greater_zero(m_wc.percolation_sl[sl - 1]) &&
2256 cbm::flt_greater(wc_.wc_sl[sl], wc_.wc_wp_sl[sl]))
2257 {
2258 m_wc.percolation_sl[sl] = m_sipar->FPERCOL() * m_wc.percolation_sl[sl - 1];
2259 }
2260
2261 //bound water flux by defined maximal percolation rate (e.g., during flooding events)
2262 m_wc.percolation_sl[sl] = cbm::bound_max(m_wc.percolation_sl[sl], max_percolation);
2263 }
2264 }
2265
2266 //infiltrating water in current time step is exceeding available pore space of last time step
2267 //second condition: security, should be always true if first condition is met
2268 else if (cbm::flt_greater(wc_.wc_sl[sl], sc_.poro_sl[sl]) &&
2269 cbm::flt_greater(wc_.wc_sl[sl], wc_.wc_fc_sl[sl]))
2270 {
2271 //maximum allowed flow until field capacity
2272 double const max_flow((wc_.wc_sl[sl] - wc_.wc_fc_sl[sl]) * sc_.h_sl[sl]);
2273 m_wc.percolation_sl[sl] = cbm::bound_max(WaterFlowSaturated(sl, get_impedance_factor(sl)),
2274 max_flow);
2275 }
2276 else if (cbm::flt_greater_equal(wc_.wc_sl[sl] + wc_.ice_sl[sl], wc_.wc_fc_sl[sl]))
2277 {
2278 m_wc.percolation_sl[sl] = WaterFlowAboveFieldCapacity(sl, get_impedance_factor(sl));
2279 }
2280 else if (cbm::flt_greater(wc_.wc_sl[sl], wc_.wc_wp_sl[sl]))
2281 {
2282 m_wc.percolation_sl[sl] = WaterFlowCapillaryRise(sl, get_impedance_factor(sl));
2283 }
2284 else
2285 {
2286 m_wc.percolation_sl[sl] = 0.0;
2287 }
2288
2289
2290 // water content must be greater or equal wilting point
2291 if (cbm::flt_less((wc_.wc_sl[sl] + wc_.ice_sl[sl]) * sc_.h_sl[sl] - m_wc.percolation_sl[sl], wc_.wc_wp_sl[sl] * sc_.h_sl[sl]) &&
2292 cbm::flt_greater_zero(m_wc.percolation_sl[sl]))
2293 {
2294 m_wc.percolation_sl[sl] = cbm::bound(0.0,
2295 (wc_.wc_sl[sl] + wc_.ice_sl[sl] - wc_.wc_wp_sl[sl]) * sc_.h_sl[sl],
2296 m_wc.percolation_sl[sl]);
2297 }
2298
2299 // Update of the current soil layer water content
2300 wc_.wc_sl[sl] -= m_wc.percolation_sl[sl] / sc_.h_sl[sl];
2301
2302 SoilWaterEvaporation(sl);
2303 }
2304
2305 /* Correct water content and percolation if water content is above porosity
2306 * If the layer below cannot take up all of the water assigned for percolation,
2307 * so if the air filled pore space is not sufficient,
2308 * the water remains in the upper layer.
2309 */
2310 for (size_t sl = m_soillayers->soil_layer_cnt() - 1; sl > 0; sl--)
2311 {
2312 double const ice_vol(wc_.ice_sl[sl] / cbm::DICE);
2313 if (cbm::flt_greater(wc_.wc_sl[sl] + ice_vol, sc_.poro_sl[sl]))
2314 {
2315 double delta_wc(wc_.wc_sl[sl] + ice_vol - sc_.poro_sl[sl]);
2316 if (cbm::flt_greater_equal(ice_vol, sc_.poro_sl[sl]))
2317 {
2318 delta_wc = wc_.wc_sl[sl];
2319 wc_.wc_sl[sl] = 0.0;
2320 }
2321 else
2322 {
2323 wc_.wc_sl[sl] = sc_.poro_sl[sl] - ice_vol;
2324 }
2325 m_wc.percolation_sl[sl - 1] -= delta_wc * sc_.h_sl[sl];
2326 wc_.wc_sl[sl - 1] += delta_wc * sc_.h_sl[sl] / sc_.h_sl[sl - 1];
2327 }
2328 }
2329
2330 double const ice_vol(wc_.ice_sl[0] / cbm::DICE);
2331 if (cbm::flt_greater(wc_.wc_sl[0] + ice_vol, sc_.poro_sl[0]))
2332 {
2333
2334 double delta_wc(wc_.wc_sl[0] + ice_vol - sc_.poro_sl[0]);
2335 if (cbm::flt_greater_equal(ice_vol, sc_.poro_sl[0]))
2336 {
2337 delta_wc = wc_.wc_sl[0];
2338 wc_.wc_sl[0] = 0.0;
2339 }
2340 else
2341 {
2342 wc_.wc_sl[0] = sc_.poro_sl[0] - ice_vol;
2343 }
2344 wc_.accumulated_infiltration -= delta_wc * sc_.h_sl[0];
2345 wc_.surface_water += delta_wc * sc_.h_sl[0];
2346 }
2347
2348 if (cbm::flt_greater_zero(wc_.surface_water))
2349 {
2350 CalcSurfaceFlux();
2351 }
2352
2353 wc_.accumulated_waterflux_sl += m_wc.percolation_sl;
2354 wc_.accumulated_percolation += m_wc.percolation_sl[ m_soillayers->soil_layer_cnt()-1];
2355}
ldndc::WatercycleDNDC::WaterFlowAboveFieldCapacity
double WaterFlowAboveFieldCapacity(size_t, double const &)
Definition watercycle-dndc.cpp:2368
ldndc::WatercycleDNDC::WaterFlowCapillaryRise
double WaterFlowCapillaryRise(size_t, double const &)
Similar to WaterFlowBelowFieldCapacity(...) but restricts water flow by water availability from soil ...
Definition watercycle-dndc.cpp:2417
ldndc::WatercycleDNDC::CalcSurfaceFlux
void CalcSurfaceFlux()
Definition watercycle-dndc.cpp:2572
ldndc::WatercycleDNDC::SoilWaterEvaporation
void SoilWaterEvaporation(size_t)
Definition watercycle-dndc.cpp:2088
ldndc::WatercycleDNDC::get_impedance_factor
double get_impedance_factor(size_t)
Reduces water fluxes out of a layer if ice lenses have formed.
Definition watercycle-dndc.cpp:2523

References CalcSoilEvapoPercolation(), CalcSurfaceFlux(), get_impedance_factor(), SoilWaterEvaporation(), WaterFlowAboveFieldCapacity(), and WaterFlowCapillaryRise().

Referenced by CalcSoilEvapoPercolation(), CalcTranspiration(), and solve().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcSurfaceFlux()

void ldndc::WatercycleDNDC::CalcSurfaceFlux ( )
private
Parameters
[in]None
[out]None
Returns
None
2573{
2574 if ( cbm::flt_greater_zero( m_eventflood.get_bund_height()))
2575 {
2576 if ( cbm::flt_greater( wc_.surface_water, m_eventflood.get_bund_height()))
2577 {
2578 wc_.accumulated_runoff += wc_.surface_water - m_eventflood.get_bund_height();
2579 wc_.surface_water = m_eventflood.get_bund_height();
2580 }
2581 }
2582 else
2583 {
2584 /* runoff fraction per time step equals RO fraction per hour * length of a time step */
2585 double const runoff_fraction(m_sipar->FRUNOFF() * get_time_step_duration());
2586
2587 // - not all surface water is removed within one timestep
2588 if ( cbm::flt_less( runoff_fraction, 1.0))
2589 {
2590 double const runoff( runoff_fraction * wc_.surface_water);
2591 wc_.surface_water -= runoff;
2592 wc_.accumulated_runoff += runoff;
2593 }
2594 // - all surface water is immediatly removed
2595 else
2596 {
2597 wc_.accumulated_runoff += wc_.surface_water;
2598 wc_.surface_water = 0.0;
2599 }
2600 }
2601}
ldndc::WatercycleDNDC::get_time_step_duration
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition watercycle-dndc.cpp:2648

References CalcSurfaceFlux(), and get_time_step_duration().

Referenced by CalcSoilEvapoPercolation(), CalcSurfaceFlux(), and CalcTranspiration().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcSurfaceWaterEvaporation()

void ldndc::WatercycleDNDC::CalcSurfaceWaterEvaporation ( )
private
Parameters
[in]None
[out]None
Returns
None
1992{
1993 if ( cbm::flt_greater_zero( wc_.surface_water) &&
1994 cbm::flt_greater_zero( hr_pot_evapotrans))
1995 {
1996 if ( hr_pot_evapotrans > wc_.surface_water)
1997 {
1998 wc_.accumulated_surfacewaterevaporation += wc_.surface_water;
1999 hr_pot_evapotrans -= wc_.surface_water;
2000 wc_.surface_water = 0.0;
2001 }
2002 else
2003 {
2004 wc_.surface_water -= hr_pot_evapotrans;
2005 wc_.accumulated_surfacewaterevaporation += hr_pot_evapotrans;
2006 hr_pot_evapotrans = 0.0;
2007 }
2008 }
2009}

References CalcSurfaceWaterEvaporation().

Referenced by CalcSurfaceWaterEvaporation(), CalcTranspiration(), and solve().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcTranspiration()

lerr_t ldndc::WatercycleDNDC::CalcTranspiration ( )
private
Parameters
[in]None
[out]None
Returns
None

References CalcSnowEvaporation(), CalcSoilEvapoPercolation(), CalcSurfaceFlux(), CalcSurfaceWaterEvaporation(), CalcTranspirationCouvreur(), get_clay_limitation(), get_fsl(), get_time_step_duration(), SoilWaterEvaporation(), WatercycleDNDC_preferential_flow(), WaterFlowAboveFieldCapacity(), and WaterFlowCapillaryRise().

Referenced by solve().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalcTranspirationCouvreur()

lerr_t ldndc::WatercycleDNDC::CalcTranspirationCouvreur ( double const & hr_potTransinm)
private
Parameters
[in]None
[out]None
Returns
None
1847{
1848 // hourly potential transpiration
1849 double const hr_potTrans( hr_potTransinm * cbm::CM_IN_M);
1850
1851 // hourly water uptake
1852 double hr_uptake( 0.0);
1853
1854 // Calculate the root water uptake and hence the transpiration for every plant individually
1855 for( PlantIterator vt = m_veg->begin(); vt != m_veg->end(); ++vt)
1856 {
1857 // split up the potential transpiration equally to the various plants
1858 // frt mass of the current species
1859 double const mfrt_species( (*vt)->mFrt);
1860 // total frt mass
1861 double const mfrt_sum( m_veg->dw_frt());
1862 double const hr_potTransspecies = hr_potTrans * mfrt_species/mfrt_sum; // in cm
1863
1864 // If some soil layers are frozen and contain only ice, no fluid water is present and the water potential diverges.
1865 // The uptake from these layers is prevented by considering only roots in layers with more water than ice.
1866 // In this case the root distribution needs to be rescaled (otherwise the non-contributing layers would effectively contribute a potential = 0).
1867 double fFrt_unfrozen( 1.0);
1868
1869 /* TODO: include icetowatercrit as proper parameter */
1870 double icetowatercrit( 0.1);
1871 double effsoilwatpot( 0.0);
1872 for( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
1873 {
1874 // fine root conditions to prevent transpiration without uptake organs, wateruptake from frozen ground
1875 if( cbm::flt_greater_zero( (*vt)->fFrt_sl[sl] * mfrt_species))
1876 {
1877 if(cbm::flt_less( wc_.ice_sl[sl]/wc_.wc_sl[sl], icetowatercrit))
1878 {
1879 // overall water potential in this layer
1880 double const watpotl( Soillayerwaterhead( sl)); // negative
1881 effsoilwatpot += watpotl * (*vt)->fFrt_sl[sl]; // negative
1882 }
1883 else
1884 {
1885 fFrt_unfrozen -= (*vt)->fFrt_sl[sl];
1886 }
1887 }
1888 else
1889 {
1890 // roots are not wireless ;-)
1891 break;
1892 }
1893 }
1894 // normalize the root distribution in the potential, if some soil is frozen
1895 effsoilwatpot = effsoilwatpot / fFrt_unfrozen;
1896
1897 // hydraulic conductivities of plant system
1898 double const Krs(vt->KRSINIT() * mfrt_species / vt->FRTMASSINIT() * nhr); // cm(water)/cm(from pressure)/timestep
1899 double const Kcomp(vt->KCOMPINIT() * mfrt_species / vt->FRTMASSINIT() * nhr); // cm(water)/cm(from pressure)/timestep
1900
1901 // water potential in the leafs, depending on the potential transpiration
1902 double const leafwatpot( cbm::bound_min(vt->HLEAFCRIT(), effsoilwatpot - (hr_potTransspecies / Krs))); // cm
1903
1904 if( cbm::flt_greater_zero( leafwatpot))
1905 {
1906 KLOGERROR( "leafwatpot is positive");
1907 return LDNDC_ERR_FAIL;
1908 }
1909 if( !cbm::flt_greater_zero( effsoilwatpot - leafwatpot))
1910 {
1911 KLOGERROR( "effsoilwatpot - leafwatpot is smaller 0: ", effsoilwatpot - leafwatpot, ", effsoilwatpot = ", effsoilwatpot, ", leafwatpot = ", leafwatpot, "\n-> Hourly transpiration ", hr_uptake, " is negative!\nPlant dies :-( Try with a smaller KRSINIT or different KCOMPINIT or different EXP_ROOT_DISTRIBUTION.\nOr wc_wp is to close to the residual water content.. increase it.");
1912 return LDNDC_ERR_FAIL;
1913 }
1914
1915 double const effuptake( Krs * (effsoilwatpot - leafwatpot)); // cm/timestep
1916 double uptake_tot( 0.0);
1917 double uptakecomp_tot( 0.0);
1918// double uptakecompabs_tot( 0.0);
1919
1920 for( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
1921 {
1922 // fine root condition to prevent transpiration without uptake organs
1923 if( cbm::flt_greater_zero( (*vt)->fFrt_sl[sl] * mfrt_species))
1924 {
1925 if(cbm::flt_less( wc_.ice_sl[sl]/wc_.wc_sl[sl], icetowatercrit))
1926 {
1927 double const watpotl( Soillayerwaterhead( sl));
1928 double const compuptake( Kcomp * (watpotl - effsoilwatpot));
1929 double uptake_layer( (effuptake + compuptake) * (*vt)->fFrt_sl[sl]); // absolute value in cm
1930
1931 wc_.accumulated_wateruptake_sl[sl] += uptake_layer * cbm::M_IN_CM;
1932 uptake_tot += uptake_layer; // absolute value in cm
1933 uptakecomp_tot += compuptake * (*vt)->fFrt_sl[sl]; // absolute value in cm
1934// uptakecompabs_tot += fabs(compuptake * (*vt)->fFrt_sl[sl]); // absolute value in cm
1935
1936 wc_.wc_sl[sl] -= uptake_layer * cbm::M_IN_CM /sc_.h_sl[sl]; // fraction of water in this layer, in m
1937
1938 if( !cbm::flt_greater_zero( wc_.wc_sl[sl]))
1939 {
1940 KLOGERROR( "In soil layer ", sl, " the water content is negative: ", wc_.wc_sl[sl]);
1941 return LDNDC_ERR_FAIL;
1942 }
1943 }
1944 }
1945 else
1946 {
1947 // roots are not wireless ;-)
1948 break;
1949 }
1950 }
1951
1952 hr_uptake += uptake_tot * cbm::M_IN_CM; // absolute value in m for all species
1953
1954
1955 if( cbm::flt_greater_zero( fabs(uptakecomp_tot))){
1956 KLOGWARN( "Compensatory water uptake ", uptakecomp_tot, " > 0.");
1957 }
1958 if( cbm::flt_greater( uptake_tot, hr_potTransspecies, 0.0000000001))
1959 {
1960 KLOGERROR( "Total water uptake ", uptake_tot, " larger than potential transpiration ", hr_potTransspecies);
1961 }
1962 if(cbm::flt_equal_eps(fFrt_unfrozen, 1.0, 0.0000000001) && !cbm::flt_equal_eps( uptake_tot, hr_potTransspecies, 0.0000000001) && !cbm::flt_equal_eps( leafwatpot, vt->HLEAFCRIT(), 0.0000000001))
1963 {
1964 KLOGERROR( "T ", uptake_tot, "cm < Tpot ", hr_potTransspecies, "cm even though leafwatpot ", leafwatpot, " > HLEAFCRIT", vt->HLEAFCRIT());
1965 }
1966 }
1967
1968 // total evaporation demand reduced to apply on other evaporating sources (ground, snow, surface water)
1969 hr_pot_evapotrans = cbm::bound_min(0.0, hr_pot_evapotrans - hr_uptake);
1970
1971 return LDNDC_ERR_OK;
1972}

References CalcTranspirationCouvreur().

Referenced by CalcTranspiration(), and CalcTranspirationCouvreur().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ event_flood()

lerr_t ldndc::WatercycleDNDC::event_flood ( )
private

sets hydrologic conditions during flooding events, e.g.,

  • surface water table
  • bund height
775{
776 lerr_t rc = m_eventflood.solve();
777 if ( rc != LDNDC_ERR_OK)
778 {
779 KLOGERROR("Irrigation event not successful!");
780 return rc;
781 }
782
783 /* max_percolation is gradually decreased after end of flooding event to inital value */
784 double const maximum_percolation = m_eventflood.get_maximum_percolation();
785 if ( cbm::is_valid( maximum_percolation))
786 {
787 if ( cbm::flt_greater_zero( maximum_percolation))
788 {
789 max_percolation = maximum_percolation * nhr / 24.0;
790 }
791 }
792 else
793 {
794 /* time rate for gradual recovery of max_percolation */
795 double const time_rate( get_time_step_duration() * 0.1);
796 max_percolation -= (max_percolation - WaterFlowSaturated( m_soillayers->soil_layer_cnt()-1, 1.0)) * time_rate;
797 }
798
799 return LDNDC_ERR_OK;
800}

References event_flood(), and get_time_step_duration().

Referenced by event_flood(), and solve().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_leaf_water()

double ldndc::WatercycleDNDC::get_leaf_water ( )
private

Collects leaf water from all canopy layers.

Parameters
[in]None
[out]None
Returns
Leaf water
1111{
1112 if ( m_veg->size() > 0u)
1113 {
1114 return wc_.wc_fl.sum();
1115 }
1116 return 0.0;
1117}

References get_leaf_water().

Referenced by balance_check(), CalcInterception(), and get_leaf_water().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ rainfall_intensity()

double ldndc::WatercycleDNDC::rainfall_intensity ( )
private

fixed amount of maximum rainfall/irrigation triggered per time step

104{
105 if ( m_iokcomm->get_input_class< input_class_climate_t >())
106 {
107 return m_iokcomm->get_input_class< climate::input_class_climate_t >()->station_info()->rainfall_intensity * cbm::M_IN_MM;
108 }
109 else
110 {
111 return ldndc::climate::climate_info_defaults.rainfall_intensity * cbm::M_IN_MM;
112 }
113}
ldndc::climate::climate_info_t::rainfall_intensity
double rainfall_intensity
Definition climatetypes.h:40

References rainfall_intensity().

Referenced by CalcIrrigation(), CalcRaining(), and rainfall_intensity().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ SoilWaterEvaporation()

void ldndc::WatercycleDNDC::SoilWaterEvaporation ( size_t _sl)
private
Parameters
[in]_slSoil layer
[out]None
Returns
None
2090{
2091 //percolation and direct evaporation from the (upper) soil
2092 //if there is snow cover, evaporation from soil will be set to 0, and only sublimation is allowed.
2093 double const wc_min_evaporation(m_sipar->WCDNDC_EVALIM_FRAC_WCMIN() * wc_.wc_wp_sl[_sl]);
2094 if ( cbm::flt_less( sc_.depth_sl[_sl] - sc_.h_sl[_sl], m_sipar->EVALIM()) &&
2095 cbm::flt_greater( wc_.wc_sl[_sl], wc_min_evaporation) &&
2096 cbm::flt_greater_zero( hr_pot_evapotrans))
2097 {
2098 //limiting factor for soil water infiltration
2099 double const f_clay( 1.0 + m_sipar->SLOPE_CLAYF() * std::min(1.0 - sc_.fcorg_sl[_sl] / cbm::CCORG, sc_.clay_sl[_sl]));
2100 double const f_water( cbm::bound( 0.0,
2101 pow( (wc_.wc_sl[_sl] - wc_min_evaporation) / wc_.wc_fc_sl[_sl], f_clay),
2102 1.0));
2103 double const h_rest( std::min( sc_.h_sl[_sl], m_sipar->EVALIM() - ( sc_.depth_sl[_sl] - sc_.h_sl[_sl])));
2104 double const f_depth( h_rest / m_sipar->EVALIM() * std::max(0.0,
2105 1.0 - std::min((double)m_sipar->EVALIM(),
2106 sc_.depth_sl[_sl] - ( 0.5 * sc_.h_sl[_sl]))
2107 / m_sipar->EVALIM()));
2108 double const soilevaporation_max( cbm::bound_max( f_water * wc_.wc_sl[_sl] * sc_.h_sl[_sl], hr_pot_evapotrans));
2109
2110 double const soilevaporation( cbm::bound( 0.0,
2111 hr_pot_evapotrans * f_depth * f_water,
2112 soilevaporation_max));
2113
2114 //update of the current soil layer water content
2115 wc_.wc_sl[_sl] -= soilevaporation / sc_.h_sl[_sl];
2116 wc_.accumulated_soilevaporation += soilevaporation;
2117 hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - soilevaporation);
2118 }
2119}

References SoilWaterEvaporation().

Referenced by CalcSoilEvapoPercolation(), CalcTranspiration(), and SoilWaterEvaporation().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ WatercycleDNDC_preferential_flow()

void ldndc::WatercycleDNDC::WatercycleDNDC_preferential_flow ( )
private
Parameters
[in]None
[out]None
Returns
None
2462{
2463 //BY_PASSF is defined as the fraction of surface water that passes per hour (0.0 by default)
2464 if ( cbm::flt_greater_zero(m_sipar->BY_PASSF()) &&
2465 cbm::flt_greater_zero( wc_.surface_water))
2466 {
2467
2468 double water_def_total( 0.0);
2469 for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2470 {
2471 //water flux that flow through macro pores can maximally store
2472 water_def_total += cbm::bound_min( 0.0, (wc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]);
2473 }
2474
2475 //duration of timestep (hr-1)
2476 double const fhr( get_time_step_duration());
2477
2478
2479 //fraction of surface water that is put into bypass flow (m timestep-1)
2480 double const bypass_fraction( cbm::bound_max(m_sipar->BY_PASSF() * fhr, 0.99));
2481 double const bypass_water( wc_.surface_water * bypass_fraction);
2482 wc_.surface_water -= bypass_water;
2483 wc_.accumulated_infiltration += bypass_water;
2484
2485 if ( cbm::flt_greater_equal( water_def_total, bypass_water))
2486 {
2487
2488 for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2489 {
2490 double const water_def( cbm::bound_min( 0.0, (wc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]));
2491 double const water_add( water_def / water_def_total * bypass_water);
2492
2493 wc_.wc_sl[sl] += water_add / sc_.h_sl[sl];
2494 }
2495 }
2496 else
2497 {
2498 for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2499 {
2500 double const water_def( cbm::bound_min( 0.0, (wc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]));
2501 wc_.wc_sl[sl] += water_def / sc_.h_sl[sl];
2502 }
2503 /* add surplus to percolation out of last soil layer */
2504 wc_.accumulated_percolation += bypass_water - water_def_total;
2505 wc_.accumulated_waterflux_sl[m_soillayers->soil_layer_cnt()-1] += bypass_water - water_def_total;
2506 }
2507 }
2508}

References get_time_step_duration(), and WatercycleDNDC_preferential_flow().

Referenced by CalcTranspiration(), solve(), and WatercycleDNDC_preferential_flow().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ WaterFlowAboveFieldCapacity()

double ldndc::WatercycleDNDC::WaterFlowAboveFieldCapacity ( size_t _sl,
double const & _fact_impedance )
private

factor accounting for different water travelling characteristics in litter and mineral soil

2371{
2372 double const sks( WaterFlowSaturated( _sl, _fact_impedance));
2373 double const fsl( get_fsl( _sl));
2374 double const fhr( get_time_step_duration());
2375
2377 double slope(( _sl < m_soillayers->soil_layers_in_litter_cnt()) ? m_sipar->SLOPE_FF() : m_sipar->SLOPE_MS());
2378
2379 double const travel_time( 1.0 - log10( wc_.sks_sl[_sl] * cbm::M_IN_CM * cbm::MIN_IN_DAY * fhr));
2380 double const travel_time_factor( cbm::flt_greater_zero( travel_time) ? (1.0 - exp(-1.0 / travel_time)) : 1.0);
2381 return cbm::bound_max( cbm::sqr(1.0 - (wc_.wc_fc_sl[_sl] / ((wc_.wc_sl[_sl] + wc_.ice_sl[_sl]) * slope)))
2382 * wc_.wc_sl[_sl] * sc_.h_sl[_sl] * fsl * travel_time_factor * _fact_impedance,
2383 sks);
2384}
ldndc::WatercycleDNDC::get_fsl
double get_fsl(size_t)
Returns factor accounting for soil layer depth. Originally, a constant soil layer depth of 0....
Definition watercycle-dndc.cpp:2544

References get_fsl(), get_time_step_duration(), and WaterFlowAboveFieldCapacity().

Referenced by CalcSoilEvapoPercolation(), CalcTranspiration(), and WaterFlowAboveFieldCapacity().

Here is the call graph for this function:
Here is the caller graph for this function:

Member Data Documentation

◆ IMAX_W

const int unsigned ldndc::WatercycleDNDC::IMAX_W = 24
static

number of iteration steps within a day