LandscapeDNDC 1.37.0
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>

Inherits ldndc::MBE_LegacyModel.

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., More...
 
double balance_check (double *)
 checks each time step for correct water balance More...
 
void CalcRaining (unsigned int)
 apportion daily rainfall to subdaily rainfall
 
double rainfall_intensity ()
 
void CalcIrrigation (unsigned int)
 apportion daily irrigation to subdaily irrigation More...
 
lerr_t CalcIceContent ()
 Call of ice-content related processes.
 
lerr_t CalcPotEvapoTranspiration ()
 Potential evaporation. More...
 
void CalcInterception ()
 Calculates water quantity that is retained on leaf and stem bark surfaces. More...
 
double get_leaf_water ()
 Collects leaf water from all canopy layers. More...
 
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.
 

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.

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

◆ 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

◆ 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

◆ 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 .
Definition: ld_thornthwaite.cpp:41

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

Here is the call graph for this function:

◆ CalcSnowEvaporation()

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

◆ CalcSoilEvapoPercolation()

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

◆ CalcSurfaceFlux()

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

◆ CalcSurfaceWaterEvaporation()

void ldndc::WatercycleDNDC::CalcSurfaceWaterEvaporation ( )
private
Parameters
[in]None
[out]None
Returns
None
1998{
1999 if ( cbm::flt_greater_zero( wc_.surface_water) &&
2000 cbm::flt_greater_zero( hr_pot_evapotrans))
2001 {
2002 if ( hr_pot_evapotrans > wc_.surface_water)
2003 {
2004 wc_.accumulated_surfacewaterevaporation += wc_.surface_water;
2005 hr_pot_evapotrans -= wc_.surface_water;
2006 wc_.surface_water = 0.0;
2007 }
2008 else
2009 {
2010 wc_.surface_water -= hr_pot_evapotrans;
2011 wc_.accumulated_surfacewaterevaporation += hr_pot_evapotrans;
2012 hr_pot_evapotrans = 0.0;
2013 }
2014 }
2015}

◆ CalcTranspiration()

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

◆ CalcTranspirationCouvreur()

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

◆ 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}

◆ 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}

◆ 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 ldndc::climate::climate_info_t::rainfall_intensity.

◆ SoilWaterEvaporation()

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

◆ WatercycleDNDC_preferential_flow()

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

◆ 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

2377{
2378 double const sks( WaterFlowSaturated( _sl, _fact_impedance));
2379 double const fsl( get_fsl( _sl));
2380 double const fhr( get_time_step_duration());
2381
2383 double slope(( _sl < m_soillayers->soil_layers_in_litter_cnt()) ? m_sipar->SLOPE_FF() : m_sipar->SLOPE_MS());
2384
2385 double const travel_time( 1.0 - log10( wc_.sks_sl[_sl] * cbm::M_IN_CM * cbm::MIN_IN_DAY * fhr));
2386 double const travel_time_factor( cbm::flt_greater_zero( travel_time) ? (1.0 - exp(-1.0 / travel_time)) : 1.0);
2387 return cbm::bound_max( cbm::sqr(1.0 - (wc_.wc_fc_sl[_sl] / ((wc_.wc_sl[_sl] + wc_.ice_sl[_sl]) * slope)))
2388 * wc_.wc_sl[_sl] * sc_.h_sl[_sl] * fsl * travel_time_factor * _fact_impedance,
2389 sks);
2390}
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:2550

Member Data Documentation

◆ IMAX_W

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

number of iteration steps within a day