LandscapeDNDC  1.36.0
ldndc::WatercycleDNDC Class Reference

Watercycle model WatercycleDNDC. More...

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 WaterFlowSaturated (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.

2490 {
2491  double balance = get_leaf_water() + wc_.surface_ice + wc_.surface_water;
2492  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2493  {
2494  balance += (wc_.wc_sl[sl] + wc_.ice_sl[sl]) * sc_.h_sl[sl];
2495  }
2496  balance += ( - m_wc.accumulated_irrigation_event_executed
2497  - wc_.accumulated_irrigation_automatic
2498  - wc_.accumulated_irrigation_ggcmi
2499  - wc_.accumulated_irrigation_reservoir_withdrawal
2500  - wc_.accumulated_precipitation
2501  - wc_.accumulated_groundwater_access
2502  + wc_.accumulated_percolation
2503  + wc_.accumulated_runoff
2504  + wc_.accumulated_wateruptake_sl.sum()
2505  + wc_.accumulated_interceptionevaporation
2506  + wc_.accumulated_soilevaporation
2507  + wc_.accumulated_surfacewaterevaporation);
2508 
2509  if ( _balance)
2510  {
2511  double const balance_delta( std::abs( *_balance - balance));
2512  if ( balance_delta > cbm::DIFFMAX)
2513  {
2514  KLOGWARN( "Water leakage: difference is ", balance - *_balance);
2515  }
2516  }
2517 
2518  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2519  {
2520  if ( cbm::flt_less( wc_.wc_sl[sl], 0.0))
2521  {
2522  KLOGWARN( "Negative water content of: ", wc_.wc_sl[sl], " in layer: ", sl, ". Water content set to zero.");
2523  wc_.wc_sl[sl] = 0.0;
2524  }
2525  if ( cbm::flt_less( wc_.ice_sl[sl], 0.0))
2526  {
2527  KLOGWARN( "Negative ice content of: ", wc_.ice_sl[sl], " in layer: ", sl, ". Ice content set to zero.");
2528  wc_.ice_sl[sl] = 0.0;
2529  }
2530  }
2531 
2532  return balance;
2533 }
double get_leaf_water()
Collects leaf water from all canopy layers.
Definition: watercycle-dndc.cpp:1161

◆ 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

643 {
644  //reset irrigation
645  m_wc.irrigation = 0.0;
646 
650  if ( cbm::flt_in_range_lu(0.0, m_cfg.automaticirrigation, 1.0))
651  {
652  for ( size_t l = 0; l < m_soillayers->soil_layer_cnt(); ++l)
653  {
654  // trigger irrigation events as soon as drought stress appears
655  double const wc_target( sc_.wc_wp_sl[l] + m_cfg.automaticirrigation * (sc_.wc_fc_sl[l] - sc_.wc_wp_sl[l]));
656  if ( cbm::flt_less( wc_.wc_sl[l], wc_target))
657  {
658  // irrigation assumed until target water content
659  double const irrigate_layer( (wc_target - wc_.wc_sl[l]) * sc_.h_sl[l]);
660  m_wc.irrigation += irrigate_layer;
661  wc_.accumulated_irrigation_automatic += irrigate_layer;
662  }
663  }
664  }
665 
666  if ( cbm::flt_in_range_lu( 0.0, m_eventflood.get_saturation_level(), 1.0) &&
667  cbm::is_valid( m_eventflood.get_irrigation_height()) &&
668  cbm::is_valid( m_eventflood.get_soil_depth()))
669  {
670  double wc_sum( 0.0);
671  double wc_fc_sum( 0.0);
672  double wc_wp_sum( 0.0);
673  for ( size_t l = 0; l < m_soillayers->soil_layer_cnt(); ++l)
674  {
675  wc_sum += wc_.wc_sl[l] * sc_.h_sl[l];
676  wc_fc_sum += sc_.wc_fc_sl[l] * sc_.h_sl[l];
677  wc_wp_sum += sc_.wc_wp_sl[l] * sc_.h_sl[l];
678 
679  if ( cbm::flt_greater_equal( sc_.depth_sl[l], m_eventflood.get_soil_depth()))
680  {
681  break;
682  }
683  }
684 
685  // trigger irrigation events as soon as drought stress appears
686  double const wc_target( wc_wp_sum + m_eventflood.get_saturation_level() * (wc_fc_sum - wc_wp_sum));
687  if ( cbm::flt_less( wc_sum, wc_target))
688  {
689  m_wc.irrigation += m_eventflood.get_irrigation_height();
690  wc_.accumulated_irrigation_automatic += m_eventflood.get_irrigation_height();
691  }
692  }
693 
701  if ( m_cfg.ggcmi_irrigation)
702  {
703  if( (lclock()->subday() == 1) && (_i == 0))
704  {
705  if( m_veg->size() > 0 ) //check crop is present
706  {
707  double water_supply( 0.0);
708  for ( size_t l = 0; l < m_soillayers->soil_layer_cnt(); ++l)
709  {
710  if ( cbm::flt_greater_zero( m_veg->mfrt_sl( l))) //only irrigate soil layers that contain roots
711  {
712  if( cbm::flt_less( wc_.wc_sl[l], sc_.wc_fc_sl[l]) &&
713  cbm::flt_equal_zero( wc_.ice_sl[l])) //only add water if no ice is present
714  {
715  water_supply += (sc_.wc_fc_sl[l] - wc_.wc_sl[l]) * sc_.h_sl[l]; //amount of water added (m)
716  wc_.wc_sl[l] = sc_.wc_fc_sl[l];
717  }
718  }
719  }
720  wc_.accumulated_irrigation_ggcmi += water_supply;
721  }
722  }
723  }
724 
728  if ( cbm::is_valid( m_eventflood.get_water_table()))
729  {
730  /* Dynamic flooding */
731  if ( cbm::is_valid( m_eventflood.get_irrigation_height()))
732  {
733  bool irrigate( false);
734 
735  /* Irrigate after surface water table drop */
736  if ( cbm::flt_greater_equal_zero( m_eventflood.get_water_table()))
737  {
738  if ( cbm::flt_greater( m_eventflood.get_water_table(), wc_.surface_water))
739  {
740  irrigate = true;
741  }
742  }
743  else
744  {
745  if ( !cbm::flt_greater_zero( wc_.surface_water))
746  {
747  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
748  {
749  if ( cbm::flt_greater_equal( sc_.depth_sl[sl], -m_eventflood.get_water_table()))
750  {
751  //check if water content below AWD_DEPTH is saturated
752  if ( cbm::flt_less( wc_.wc_sl[sl], 0.9 * sc_.poro_sl[sl]))
753  {
754  irrigate = true;
755  }
756  break;
757  }
758  }
759  }
760  }
761  // trigger AWD irrigation event
762  if ( irrigate)
763  {
764  double water_supply( cbm::bound_min( 0.0, m_eventflood.get_irrigation_height() - wc_.surface_water));
765  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
766  {
767  // set soil fully water saturated
768  if ( cbm::flt_greater( sc_.poro_sl[sl], wc_.wc_sl[sl]))
769  {
770  water_supply += ((sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl]);
771  }
772  }
773 
774  //remove rainfall (less irrigation water is supplied if it if raining)
775  water_supply = cbm::bound_min( 0.0, water_supply - m_mc.rainfall);
776 
777  if ( m_eventflood.have_unlimited_water())
778  {
779  wc_.accumulated_irrigation_automatic += water_supply;
780  }
781  else
782  {
783  if ( cbm::flt_greater( wc_.irrigation_reservoir, water_supply))
784  {
785  wc_.irrigation_reservoir -= water_supply;
786  }
787  else
788  {
789  water_supply = wc_.irrigation_reservoir;
790  wc_.irrigation_reservoir = 0.0;
791  }
792  wc_.accumulated_irrigation_reservoir_withdrawal += water_supply;
793  }
794  m_wc.irrigation += water_supply;
795  }
796  }
797  /* Static flooding */
798  else
799  {
800  if ( !cbm::flt_greater_zero( m_eventflood.get_water_table()))
801  {
802  double water_supply( 0.0);
803  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
804  {
805  if ( cbm::flt_greater( sc_.depth_sl[sl], -m_eventflood.get_water_table()) &&
806  cbm::flt_greater( sc_.poro_sl[sl], wc_.wc_sl[sl]))
807  {
808  water_supply += (sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl];
809  }
810  }
811 
812  //remove rainfall (less irrigation water is supplied if it if raining)
813  water_supply = cbm::bound_min( 0.0, water_supply - m_mc.rainfall);
814 
815  /* Surface water will be actively drained */
816  if ( cbm::flt_greater( wc_.surface_water, water_supply))
817  {
818  wc_.accumulated_runoff += wc_.surface_water - water_supply;
819  wc_.surface_water = water_supply;
820  }
821 
822  if ( m_eventflood.have_unlimited_water())
823  {
824  wc_.accumulated_irrigation_automatic += water_supply;
825  }
826  else
827  {
828  if ( cbm::flt_greater( wc_.irrigation_reservoir, water_supply))
829  {
830  wc_.irrigation_reservoir -= water_supply;
831  }
832  else
833  {
834  water_supply = wc_.irrigation_reservoir;
835  wc_.irrigation_reservoir = 0.0;
836  }
837  wc_.accumulated_irrigation_reservoir_withdrawal += water_supply;
838  }
839  m_wc.irrigation += water_supply;
840  }
844  else if ( cbm::flt_greater( m_eventflood.get_water_table(), wc_.surface_water))
845  {
846  //adjust surface and soil water in case of flooding event
847  //water required to fill up surface water
848  double water_supply( m_eventflood.get_water_table() - wc_.surface_water);
849 
850  //add water required to fill top x layers up to saturation
851  unsigned long no_layers(3);
852  if( no_layers > m_soillayers->soil_layer_cnt())
853  {
854  no_layers = m_soillayers->soil_layer_cnt();
855  }
856 
857  for ( size_t sl = 0; sl < no_layers; ++sl)
858  {
859  water_supply += ((sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl]);
860  }
861 
862  //remove rainfall (less irrigation water is supplied if it if raining)
863  water_supply = cbm::bound_min( 0.0, water_supply - m_mc.rainfall);
864 
865  if ( m_eventflood.have_unlimited_water())
866  {
867  wc_.accumulated_irrigation_automatic += water_supply;
868  }
869  else
870  {
871  if ( cbm::flt_greater( wc_.irrigation_reservoir, water_supply))
872  {
873  wc_.irrigation_reservoir -= water_supply;
874  }
875  else
876  {
877  water_supply = wc_.irrigation_reservoir;
878  wc_.irrigation_reservoir = 0.0;
879  }
880  wc_.accumulated_irrigation_reservoir_withdrawal += water_supply;
881  }
882  m_wc.irrigation += water_supply;
883  }
884  }
885  }
886 
890  double const irrigation_scheduled( cbm::bound_min( 0.0, wc_.accumulated_irrigation - m_wc.accumulated_irrigation_event_executed));
891 
892  if ( cbm::flt_greater( irrigation_scheduled, 24.0 * rainfall_intensity()))
893  {
894  m_wc.irrigation += irrigation_scheduled / 24.0;
895  m_wc.accumulated_irrigation_event_executed += irrigation_scheduled / 24.0;
896  }
897  else if ( cbm::flt_greater( irrigation_scheduled, (double)hours_per_time_step() * rainfall_intensity()))
898  {
899  m_wc.irrigation += (double)hours_per_time_step() * rainfall_intensity();
900  m_wc.accumulated_irrigation_event_executed += (double)hours_per_time_step() * rainfall_intensity();
901  }
902  else
903  {
904  m_wc.irrigation += irrigation_scheduled;
905  m_wc.accumulated_irrigation_event_executed += irrigation_scheduled;
906  }
907 }
double rainfall_intensity()
Definition: watercycle-dndc.cpp:115

◆ 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

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

1042 {
1043  day_potevapo = 0.0;
1044  if ( m_cfg.potentialevapotranspiration_method == "thornthwaite")
1045  {
1046  if ( lclock()->is_position( TMODE_PRE_YEARLY))
1047  {
1049  lclock()->days_in_year(),
1050  m_climate->annual_temperature_average(),
1051  m_climate->annual_temperature_amplitude());
1052  }
1053 
1054  double const daylength(ldndc::meteo::daylength(m_setup->latitude(), lclock()->yearday()));
1055  double const avg_dayl(12.0);
1056 
1057  day_potevapo = m_param->WCDNDC_INCREASE_POT_EVAPOTRANS()
1059  mc_.nd_airtemperature,
1060  daylength, avg_dayl, lclock()->days_in_month(),
1062  }
1063  else if ( (m_cfg.potentialevapotranspiration_method == "penman") ||
1064  (m_cfg.potentialevapotranspiration_method == "priestleytaylor"))
1065  {
1066  /* clock */
1067  cbm::sclock_t const & clk( this->lclock_ref());
1068 
1069  /* leaf area index */
1070  double lai( 0.0);
1071  for( PlantIterator vt = m_veg->begin(); vt != m_veg->end(); ++vt)
1072  {
1073  lai += (*vt)->lai();
1074  }
1075 
1076  /* vapour pressure [10^3:Pa] */
1077  double vps( 0.0);
1078  ldndc::meteo::vps( mc_.nd_airtemperature, &vps);
1079  double vp( vps - m_climate->vpd_day( clk));
1080 
1081  if( m_cfg.potentialevapotranspiration_method == "priestleytaylor")
1082  {
1083  day_potevapo = m_param->WCDNDC_INCREASE_POT_EVAPOTRANS()
1084  * ldndc::priestleytaylor( clk.yearday(),
1085  clk.days_in_year(),
1086  mc_.albedo,
1087  m_setup->latitude(),
1088  mc_.nd_airtemperature,
1089  mc_.nd_shortwaveradiation_in,
1090  vp,
1091  m_param->PT_ALPHA());
1092  }
1093  else if( m_cfg.potentialevapotranspiration_method == "penman")
1094  {
1095  ldndc::surface_type surface( cbm::flt_greater_zero( lai) ? short_grass : bare_soil);
1096 
1097  day_potevapo = m_param->WCDNDC_INCREASE_POT_EVAPOTRANS()
1098  * ldndc::penman( surface,
1099  clk.yearday(),
1100  clk.days_in_year(),
1101  mc_.albedo,
1102  m_setup->latitude(),
1103  mc_.nd_shortwaveradiation_in * cbm::SEC_IN_DAY,
1104  mc_.nd_airtemperature,
1105  mc_.nd_windspeed,
1106  vp);
1107  }
1108  }
1109  else
1110  {
1111  KLOGERROR("[BUG] huh :o how did you get here? ",
1112  "tell the maintainers refering to this error message");
1113  return LDNDC_ERR_FAIL;
1114  }
1115 
1116  if ( m_cfg.evapotranspiration_method == "uniform")
1117  {
1118  // hourly evaporation demand is calculated from daily demand equally throughout the day
1119  hr_pot_evapotrans = day_potevapo * nhr / cbm::HR_IN_DAY;
1120  }
1121  else if ( m_cfg.evapotranspiration_method == "previouspattern")
1122  {
1123 
1124  hr_pot_evapotrans = day_potevapo * pot_evapotranspiration_ts[lclock()->subday()-1] / nwc;
1125  }
1126  else //if ( m_cfg.evapotranspiration_method === "nice name") // could be for example 'daytime_only'
1127  {
1128  double const hour(lclock()->subday());
1129  double const day(lclock()->yearday());
1130  double const hsun(0.833 * cbm::PI / 180.0);
1131  double const lat(m_setup->latitude() * cbm::PI / 180.0);
1132  double const dec(-23.44 * (cbm::PI / 180.0) * std::cos(2.0 * cbm::PI * (day - 172.0) / lclock()->days_in_year()));
1133  double const help((std::sin(hsun) - std::sin(lat) * std::sin(dec)) / (std::cos(lat) * std::cos(dec)));
1134  double dayl(0.0);
1135  if (help < -0.999) dayl = 0.0;
1136  else if (help > 0.999) dayl = 24.0;
1137  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)));
1138  if (dayl < 1.0) dayl = 0.0;
1139  double const rise(12.0 - dayl * 0.5);
1140 
1141  double factor(0.0);
1142  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);
1143 
1144  double fday(0.0);
1145  if (hour > rise && hour < (24.0 - rise)) fday = factor;
1146 
1147  hr_pot_evapotrans = day_potevapo * fday;
1148  }
1149  // the evaporation limit for all timesteps is the daily demand
1150  wc_.accumulated_potentialevaporation += hr_pot_evapotrans;
1151 
1152  return LDNDC_ERR_OK;
1153 }
double LDNDC_API thornthwaite_heat_index(double, size_t, double, double)
Thornthwaite heat index.
Definition: ld_thornthwaite.cpp:100
double thornthwaite_heat_index
heat index for potential evaporation using approach after Thornthwaite
Definition: watercycle-dndc.h:137
double LDNDC_API thornthwaite(double, double, double, double, double)
Potential evapotranspiration after .
Definition: ld_thornthwaite.cpp:41
Here is the call graph for this function:

◆ CalcSnowEvaporation()

void ldndc::WatercycleDNDC::CalcSnowEvaporation ( )
private
Parameters
[in]None
[out]None
Returns
None
1870 {
1871  //Average limit of snowEvaporation is 0.00013m day-1 (Granger and Male 1978)
1872  double hr_potESnow( cbm::bound(0.0, hr_pot_evapotrans, 0.00013 * nhr/24.0));
1873 
1874  if ( cbm::flt_greater_zero( wc_.surface_ice) &&
1875  cbm::flt_greater_zero( hr_potESnow))
1876  {
1877  //If there is a snowpack and there is a potential to evaporate
1878  if ( wc_.surface_ice > hr_potESnow)
1879  {
1880  //Partial evaporation from the snow
1881  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - hr_potESnow);
1882  wc_.accumulated_surfacewaterevaporation += hr_potESnow;
1883  wc_.surface_ice -= hr_potESnow;
1884  }
1885  else
1886  {
1887  //Total evaporation of the snowpack
1888  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - wc_.surface_ice);;
1889  wc_.accumulated_surfacewaterevaporation += wc_.surface_ice;
1890  wc_.surface_ice = 0.0;
1891  }
1892  }
1893 }

◆ CalcSoilEvapoPercolation()

void ldndc::WatercycleDNDC::CalcSoilEvapoPercolation ( )
private
Parameters
[in]None
[out]None
Returns
None
2008 {
2009  /* Reset water fluxes */
2010  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2011  {
2012  m_wc.percolation_sl[sl] = 0.0;
2013  }
2014 
2015  /* Fill groundwater layer with water */
2016  for ( int sl = m_soillayers->soil_layer_cnt()-1; sl >= 0; --sl)
2017  {
2018  if ( cbm::flt_greater_equal( sc_.depth_sl[sl], ground_water_table))
2019  {
2020  /* ensure that ice volume is less than porosity */
2021  double const ice_vol( wc_.ice_sl[sl] / cbm::DICE);
2022  if ( cbm::flt_greater( sc_.poro_sl[sl], ice_vol))
2023  {
2024  //air filled porosity is filled with groundwater
2025  double const gw_access( cbm::bound_min( 0.0,
2026  (sc_.poro_sl[sl] - ice_vol - wc_.wc_sl[sl]) * sc_.h_sl[sl]));
2027  wc_.wc_sl[sl] += gw_access / sc_.h_sl[sl];
2028  wc_.accumulated_groundwater_access += gw_access;
2029  }
2030  }
2031  else
2032  {
2033  break;
2034  }
2035  }
2036 
2037  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2038  {
2039  if ( sl == 0)
2040  {
2041  wc_.accumulated_infiltration += wc_.surface_water;
2042  wc_.wc_sl[0] += wc_.surface_water / sc_.h_sl[0];
2043  wc_.surface_water = 0.0;
2044  }
2045  else
2046  {
2047  wc_.wc_sl[sl] += m_wc.percolation_sl[sl-1] / sc_.h_sl[sl];
2048  }
2049 
2050  if ( cbm::flt_greater( sc_.depth_sl[sl], ground_water_table))
2051  {
2052  //set constant flux for all ground water layers depending on last non-groundwater layer
2053  if ( sl > 0)
2054  {
2055  m_wc.percolation_sl[sl] = m_wc.percolation_sl[sl-1];
2056  }
2057 
2058  //in case of groundwater until soil surface first outflux is estimated by saturated water flow
2059  else
2060  {
2061  m_wc.percolation_sl[sl] = WaterFlowSaturated( sl, get_impedance_factor( sl));
2062  }
2063  }
2064  else if ( sl + 1 == m_soillayers->soil_layer_cnt())
2065  {
2066  m_wc.percolation_sl[sl] = 0.0;
2067  if ( cbm::flt_greater_zero( wc_.sks_sl[sl]))
2068  {
2069  if ( cbm::flt_greater( wc_.wc_sl[sl], sc_.poro_sl[sl]) &&
2070  cbm::flt_greater( wc_.wc_sl[sl], sc_.wc_fc_sl[sl])) //security, should be always true if first condition is met
2071  {
2072  //maximum allowed flow until field capacity
2073  double const max_flow( (wc_.wc_sl[sl] - sc_.wc_fc_sl[sl]) * sc_.h_sl[sl]);
2074  m_wc.percolation_sl[sl] = cbm::bound_max( WaterFlowSaturated( sl, get_impedance_factor( sl)),
2075  max_flow);
2076  }
2077 
2078  //water flux above field capacity
2079  else if ( cbm::flt_greater_equal( wc_.wc_sl[sl] + wc_.ice_sl[sl], sc_.wc_fc_sl[sl]))
2080  {
2081  m_wc.percolation_sl[sl] = WaterFlowAboveFieldCapacity( sl, get_impedance_factor( sl));
2082  }
2083 
2084  //water flux below field capacity and above wilting point
2085  else if ( cbm::flt_greater_zero( m_wc.percolation_sl[sl-1])
2086  && cbm::flt_greater( wc_.wc_sl[sl], sc_.wc_wp_sl[sl]))
2087  {
2088  m_wc.percolation_sl[sl] = m_param->FPERCOL() * m_wc.percolation_sl[sl-1];
2089  }
2090 
2091  //bound water flux by defined maximal percolation rate (e.g., during flooding events)
2092  m_wc.percolation_sl[sl] = cbm::bound_max( m_wc.percolation_sl[sl], max_percolation);
2093  }
2094  }
2095 
2096  //infiltrating water in current time step is exceeding available pore space of last time step
2097  //second condition: security, should be always true if first condition is met
2098  else if ( cbm::flt_greater( wc_.wc_sl[sl], sc_.poro_sl[sl]) &&
2099  cbm::flt_greater( wc_.wc_sl[sl], sc_.wc_fc_sl[sl]))
2100  {
2101  //maximum allowed flow until field capacity
2102  double const max_flow( (wc_.wc_sl[sl] - sc_.wc_fc_sl[sl]) * sc_.h_sl[sl]);
2103  m_wc.percolation_sl[sl] = cbm::bound_max( WaterFlowSaturated( sl, get_impedance_factor( sl)),
2104  max_flow);
2105  }
2106  else if ( cbm::flt_greater_equal( wc_.wc_sl[sl] + wc_.ice_sl[sl], sc_.wc_fc_sl[sl]))
2107  {
2108  m_wc.percolation_sl[sl] = WaterFlowAboveFieldCapacity( sl, get_impedance_factor( sl));
2109  }
2110  else if ( cbm::flt_greater( wc_.wc_sl[sl], sc_.wc_wp_sl[sl]))
2111  {
2112  m_wc.percolation_sl[sl] = WaterFlowCapillaryRise( sl, get_impedance_factor( sl));
2113  }
2114  else
2115  {
2116  m_wc.percolation_sl[sl] = 0.0;
2117  }
2118 
2119 
2120  // water content must be greater or equal wilting point
2121  if ( cbm::flt_less( (wc_.wc_sl[sl] + wc_.ice_sl[sl]) * sc_.h_sl[sl] - m_wc.percolation_sl[sl],
2122  sc_.wc_wp_sl[sl] * sc_.h_sl[sl]) &&
2123  cbm::flt_greater_zero( m_wc.percolation_sl[sl]))
2124  {
2125  m_wc.percolation_sl[sl] = cbm::bound( 0.0,
2126  (wc_.wc_sl[sl] + wc_.ice_sl[sl] - sc_.wc_wp_sl[sl]) * sc_.h_sl[sl],
2127  m_wc.percolation_sl[sl]);
2128  }
2129 
2130  // Update of the current soil layer water content
2131  wc_.wc_sl[sl] -= m_wc.percolation_sl[sl] / sc_.h_sl[sl];
2132 
2133  SoilWaterEvaporation( sl);
2134  }
2135 
2136  /* Correct water content and percolation if water content is above porosity
2137  * If the layer below cannot take up all of the water assigned for percolation,
2138  * so if the air filled pore space is not sufficient,
2139  * the water remains in the upper layer.
2140  */
2141  for ( size_t sl = m_soillayers->soil_layer_cnt()-1; sl > 0; sl--)
2142  {
2143  double const ice_vol( wc_.ice_sl[sl] / cbm::DICE);
2144  if ( cbm::flt_greater( wc_.wc_sl[sl] + ice_vol,
2145  sc_.poro_sl[sl]))
2146  {
2147  double delta_wc( wc_.wc_sl[sl] + ice_vol - sc_.poro_sl[sl]);
2148  if ( cbm::flt_greater_equal( ice_vol, sc_.poro_sl[sl]))
2149  {
2150  delta_wc = wc_.wc_sl[sl];
2151  wc_.wc_sl[sl] = 0.0;
2152  }
2153  else
2154  {
2155  wc_.wc_sl[sl] = sc_.poro_sl[sl] - ice_vol;
2156  }
2157  m_wc.percolation_sl[sl-1] -= delta_wc * sc_.h_sl[sl];
2158  wc_.wc_sl[sl-1] += delta_wc * sc_.h_sl[sl] / sc_.h_sl[sl-1];
2159  }
2160  }
2161 
2162  double const ice_vol( wc_.ice_sl[0] / cbm::DICE);
2163  if ( cbm::flt_greater( wc_.wc_sl[0] + ice_vol,
2164  sc_.poro_sl[0]))
2165  {
2166 
2167  double delta_wc( wc_.wc_sl[0] + ice_vol - sc_.poro_sl[0]);
2168  if ( cbm::flt_greater_equal( ice_vol, sc_.poro_sl[0]))
2169  {
2170  delta_wc = wc_.wc_sl[0];
2171  wc_.wc_sl[0] = 0.0;
2172  }
2173  else
2174  {
2175  wc_.wc_sl[0] = sc_.poro_sl[0] - ice_vol;
2176  }
2177  wc_.accumulated_infiltration -= delta_wc * sc_.h_sl[0];
2178  wc_.surface_water += delta_wc * sc_.h_sl[0];
2179  }
2180 
2181 
2182  CalcSurfaceFlux();
2183 
2184  wc_.accumulated_waterflux_sl += m_wc.percolation_sl;
2185  wc_.accumulated_percolation += m_wc.percolation_sl[ m_soillayers->soil_layer_cnt()-1];
2186 }
double WaterFlowSaturated(size_t, double const &)
Definition: watercycle-dndc.cpp:2278
double get_impedance_factor(size_t)
Reduces water fluxes out of a layer if ice lenses have formed.
Definition: watercycle-dndc.cpp:2354
void SoilWaterEvaporation(size_t)
Definition: watercycle-dndc.cpp:1935
double WaterFlowAboveFieldCapacity(size_t, double const &)
Definition: watercycle-dndc.cpp:2199
void CalcSurfaceFlux()
Definition: watercycle-dndc.cpp:2403
double WaterFlowCapillaryRise(size_t, double const &)
Similar to WaterFlowBelowFieldCapacity(...) but restricts water flow by water availability from soil ...
Definition: watercycle-dndc.cpp:2248

◆ CalcSurfaceFlux()

void ldndc::WatercycleDNDC::CalcSurfaceFlux ( )
private
Parameters
[in]None
[out]None
Returns
None
2404 {
2405  if ( cbm::flt_greater_zero( m_eventflood.get_bund_height()))
2406  {
2407  if ( cbm::flt_greater( wc_.surface_water, m_eventflood.get_bund_height()))
2408  {
2409  wc_.accumulated_runoff += wc_.surface_water - m_eventflood.get_bund_height();
2410  wc_.surface_water = m_eventflood.get_bund_height();
2411  }
2412  }
2413  else
2414  {
2415  /* runoff fraction per time step equals RO fraction per hour * length of a time step */
2416  double const runoff_fraction( m_param->FRUNOFF() * get_time_step_duration());
2417 
2418  // if frunoff_ts < 1
2419  if ( cbm::flt_less( runoff_fraction, 1.0))
2420  {
2421  double const runoff( runoff_fraction * wc_.surface_water);
2422  wc_.surface_water -= runoff;
2423  wc_.accumulated_runoff += runoff;
2424  }
2425  else
2426  {
2427  wc_.accumulated_runoff += wc_.surface_water;
2428  wc_.surface_water = 0.0;
2429  }
2430  }
2431 }
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2478

◆ CalcSurfaceWaterEvaporation()

void ldndc::WatercycleDNDC::CalcSurfaceWaterEvaporation ( )
private
Parameters
[in]None
[out]None
Returns
None
1839 {
1840  if ( cbm::flt_greater_zero( wc_.surface_water) &&
1841  cbm::flt_greater_zero( hr_pot_evapotrans))
1842  {
1843  if ( hr_pot_evapotrans > wc_.surface_water)
1844  {
1845  wc_.accumulated_surfacewaterevaporation += wc_.surface_water;
1846  hr_pot_evapotrans -= wc_.surface_water;
1847  wc_.surface_water = 0.0;
1848  }
1849  else
1850  {
1851  wc_.surface_water -= hr_pot_evapotrans;
1852  wc_.accumulated_surfacewaterevaporation += hr_pot_evapotrans;
1853  hr_pot_evapotrans = 0.0;
1854  }
1855  }
1856 }

◆ 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
1694 {
1695  // hourly potential transpiration
1696  double const hr_potTrans( hr_potTransinm * cbm::CM_IN_M);
1697 
1698  // hourly water uptake
1699  double hr_uptake( 0.0);
1700 
1701  // Calculate the root water uptake and hence the transpiration for every plant individually
1702  for( PlantIterator vt = m_veg->begin(); vt != m_veg->end(); ++vt)
1703  {
1704  // split up the potential transpiration equally to the various plants
1705  // frt mass of the current species
1706  double const mfrt_species( (*vt)->mFrt);
1707  // total frt mass
1708  double const mfrt_sum( m_veg->dw_frt());
1709  double const hr_potTransspecies = hr_potTrans * mfrt_species/mfrt_sum; // in cm
1710 
1711  // If some soil layers are frozen and contain only ice, no fluid water is present and the water potential diverges.
1712  // The uptake from these layers is prevented by considering only roots in layers with more water than ice.
1713  // In this case the root distribution needs to be rescaled (otherwise the non-contributing layers would effectively contribute a potential = 0).
1714  double fFrt_unfrozen( 1.0);
1715 
1716  /* TODO: include icetowatercrit as proper parameter */
1717  double icetowatercrit( 0.1);
1718  double effsoilwatpot( 0.0);
1719  for( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
1720  {
1721  // fine root conditions to prevent transpiration without uptake organs, wateruptake from frozen ground
1722  if( cbm::flt_greater_zero( (*vt)->fFrt_sl[sl] * mfrt_species))
1723  {
1724  if(cbm::flt_less( wc_.ice_sl[sl]/wc_.wc_sl[sl], icetowatercrit))
1725  {
1726  // overall water potential in this layer
1727  double const watpotl( Soillayerwaterhead( sl)); // negative
1728  effsoilwatpot += watpotl * (*vt)->fFrt_sl[sl]; // negative
1729  }
1730  else
1731  {
1732  fFrt_unfrozen -= (*vt)->fFrt_sl[sl];
1733  }
1734  }
1735  else
1736  {
1737  // roots are not wireless ;-)
1738  break;
1739  }
1740  }
1741  // normalize the root distribution in the potential, if some soil is frozen
1742  effsoilwatpot = effsoilwatpot / fFrt_unfrozen;
1743 
1744  // hydraulic conductivities of plant system
1745  double const Krs(vt->KRSINIT() * mfrt_species / vt->FRTMASSINIT() * nhr); // cm(water)/cm(from pressure)/timestep
1746  double const Kcomp(vt->KCOMPINIT() * mfrt_species / vt->FRTMASSINIT() * nhr); // cm(water)/cm(from pressure)/timestep
1747 
1748  // water potential in the leafs, depending on the potential transpiration
1749  double const leafwatpot( cbm::bound_min(vt->HLEAFCRIT(), effsoilwatpot - (hr_potTransspecies / Krs))); // cm
1750 
1751  if( cbm::flt_greater_zero( leafwatpot))
1752  {
1753  KLOGERROR( "leafwatpot is positive");
1754  return LDNDC_ERR_FAIL;
1755  }
1756  if( !cbm::flt_greater_zero( effsoilwatpot - leafwatpot))
1757  {
1758  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.");
1759  return LDNDC_ERR_FAIL;
1760  }
1761 
1762  double const effuptake( Krs * (effsoilwatpot - leafwatpot)); // cm/timestep
1763  double uptake_tot( 0.0);
1764  double uptakecomp_tot( 0.0);
1765 // double uptakecompabs_tot( 0.0);
1766 
1767  for( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
1768  {
1769  // fine root condition to prevent transpiration without uptake organs
1770  if( cbm::flt_greater_zero( (*vt)->fFrt_sl[sl] * mfrt_species))
1771  {
1772  if(cbm::flt_less( wc_.ice_sl[sl]/wc_.wc_sl[sl], icetowatercrit))
1773  {
1774  double const watpotl( Soillayerwaterhead( sl));
1775  double const compuptake( Kcomp * (watpotl - effsoilwatpot));
1776  double uptake_layer( (effuptake + compuptake) * (*vt)->fFrt_sl[sl]); // absolute value in cm
1777 
1778  wc_.accumulated_wateruptake_sl[sl] += uptake_layer * cbm::M_IN_CM;
1779  uptake_tot += uptake_layer; // absolute value in cm
1780  uptakecomp_tot += compuptake * (*vt)->fFrt_sl[sl]; // absolute value in cm
1781 // uptakecompabs_tot += fabs(compuptake * (*vt)->fFrt_sl[sl]); // absolute value in cm
1782 
1783  wc_.wc_sl[sl] -= uptake_layer * cbm::M_IN_CM /sc_.h_sl[sl]; // fraction of water in this layer, in m
1784 
1785  if( !cbm::flt_greater_zero( wc_.wc_sl[sl]))
1786  {
1787  KLOGERROR( "In soil layer ", sl, " the water content is negative: ", wc_.wc_sl[sl]);
1788  return LDNDC_ERR_FAIL;
1789  }
1790  }
1791  }
1792  else
1793  {
1794  // roots are not wireless ;-)
1795  break;
1796  }
1797  }
1798 
1799  hr_uptake += uptake_tot * cbm::M_IN_CM; // absolute value in m for all species
1800 
1801 
1802  if( cbm::flt_greater_zero( fabs(uptakecomp_tot))){
1803  KLOGWARN( "Compensatory water uptake ", uptakecomp_tot, " > 0.");
1804  }
1805  if( cbm::flt_greater( uptake_tot, hr_potTransspecies, 0.0000000001))
1806  {
1807  KLOGERROR( "Total water uptake ", uptake_tot, " larger than potential transpiration ", hr_potTransspecies);
1808  }
1809  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))
1810  {
1811  KLOGERROR( "T ", uptake_tot, "cm < Tpot ", hr_potTransspecies, "cm even though leafwatpot ", leafwatpot, " > HLEAFCRIT", vt->HLEAFCRIT());
1812  }
1813  }
1814 
1815  // total evaporation demand reduced to apply on other evaporating sources (ground, snow, surface water)
1816  hr_pot_evapotrans = cbm::bound_min(0.0, hr_pot_evapotrans - hr_uptake);
1817 
1818  return LDNDC_ERR_OK;
1819 }

◆ event_flood()

lerr_t ldndc::WatercycleDNDC::event_flood ( )
private

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

  • surface water table
  • bund height
613 {
614  lerr_t rc = m_eventflood.solve();
615  if ( rc != LDNDC_ERR_OK)
616  {
617  KLOGERROR("Irrigation event not successful!");
618  return rc;
619  }
620 
621  /* max_percolation is gradually decreased after end of flooding event to inital value */
622  double const maximum_percolation = m_eventflood.get_maximum_percolation();
623  if ( cbm::is_valid( maximum_percolation))
624  {
625  if ( cbm::flt_greater_zero( maximum_percolation))
626  {
627  max_percolation = maximum_percolation * nhr / 24.0;
628  }
629  }
630  else
631  {
632  /* time rate for gradual recovery of max_percolation */
633  double const time_rate( get_time_step_duration() * 0.1);
634  max_percolation -= (max_percolation - WaterFlowSaturated( m_soillayers->soil_layer_cnt()-1, 1.0)) * time_rate;
635  }
636 
637  return LDNDC_ERR_OK;
638 }
double WaterFlowSaturated(size_t, double const &)
Definition: watercycle-dndc.cpp:2278
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2478

◆ 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

References ldndc::site::site_info_t::pressure_at_field_capacity, and ldndc::rootdensityfortranspiration_sl().

1162 {
1163  if ( m_veg->size() > 0u)
1164  {
1165  return wc_.wc_fl.sum();
1166  }
1167  return 0.0;
1168 }
Here is the call graph for this function:

◆ rainfall_intensity()

double ldndc::WatercycleDNDC::rainfall_intensity ( )
private

fixed amount of maximum rainfall/irrigation triggered per time step

References ldndc::climate::climate_info_t::rainfall_intensity, and ldndc::thornthwaite_heat_index().

116 {
117  if ( m_iokcomm->get_input_class< input_class_climate_t >())
118  {
119  return m_iokcomm->get_input_class< climate::input_class_climate_t >()->station_info()->rainfall_intensity * cbm::M_IN_MM;
120  }
121  else
122  {
123  return ldndc::climate::climate_info_defaults.rainfall_intensity * cbm::M_IN_MM;
124  }
125 }
double rainfall_intensity
Definition: climatetypes.h:40
Here is the call graph for this function:

◆ SoilWaterEvaporation()

void ldndc::WatercycleDNDC::SoilWaterEvaporation ( size_t  _sl)
private
Parameters
[in]_slSoil layer
[out]None
Returns
None
1937 {
1938  //percolation and direct evaporation from the (upper) soil
1939  //if there is snow cover, evaporation from soil will be set to 0, and only sublimation is allowed.
1940  double const wc_min_evaporation( m_param->WCDNDC_EVALIM_FRAC_WCMIN() * sc_.wc_wp_sl[_sl]);
1941  if ( cbm::flt_less( sc_.depth_sl[_sl] - sc_.h_sl[_sl], m_param->EVALIM()) &&
1942  cbm::flt_greater( wc_.wc_sl[_sl], wc_min_evaporation) &&
1943  cbm::flt_greater_zero( hr_pot_evapotrans))
1944  {
1945  //limiting factor for soil water infiltration
1946  double const f_clay( 1.0 + m_param->SLOPE_CLAYF() * std::min(1.0 - sc_.fcorg_sl[_sl] / cbm::CCORG, sc_.clay_sl[_sl]));
1947  double const f_water( cbm::bound( 0.0,
1948  pow( (wc_.wc_sl[_sl] - wc_min_evaporation) / sc_.wc_fc_sl[_sl], f_clay),
1949  1.0));
1950  double const h_rest( std::min( sc_.h_sl[_sl], m_param->EVALIM() - ( sc_.depth_sl[_sl] - sc_.h_sl[_sl])));
1951  double const f_depth( h_rest / m_param->EVALIM() * std::max(0.0,
1952  1.0 - std::min((double)m_param->EVALIM(),
1953  sc_.depth_sl[_sl] - ( 0.5 * sc_.h_sl[_sl]))
1954  / m_param->EVALIM()));
1955  double const soilevaporation_max( cbm::bound_max( f_water * wc_.wc_sl[_sl] * sc_.h_sl[_sl], hr_pot_evapotrans));
1956 
1957  double const soilevaporation( cbm::bound( 0.0,
1958  hr_pot_evapotrans * f_depth * f_water,
1959  soilevaporation_max));
1960 
1961  //update of the current soil layer water content
1962  wc_.wc_sl[_sl] -= soilevaporation / sc_.h_sl[_sl];
1963  wc_.accumulated_soilevaporation += soilevaporation;
1964  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - soilevaporation);
1965  }
1966 }

◆ WatercycleDNDC_preferential_flow()

void ldndc::WatercycleDNDC::WatercycleDNDC_preferential_flow ( )
private
Parameters
[in]None
[out]None
Returns
None
2293 {
2294  //BY_PASSF is defined as the fraction of surface water that passes per hour (0.0 by default)
2295  if ( cbm::flt_greater_zero( m_param->BY_PASSF()) &&
2296  cbm::flt_greater_zero( wc_.surface_water))
2297  {
2298 
2299  double water_def_total( 0.0);
2300  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2301  {
2302  //water flux that flow through macro pores can maximally store
2303  water_def_total += cbm::bound_min( 0.0, (sc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]);
2304  }
2305 
2306  //duration of timestep (hr-1)
2307  double const fhr( get_time_step_duration());
2308 
2309 
2310  //fraction of surface water that is put into bypass flow (m timestep-1)
2311  double const bypass_fraction( cbm::bound_max( m_param->BY_PASSF() * fhr, 0.99));
2312  double const bypass_water( wc_.surface_water * bypass_fraction);
2313  wc_.surface_water -= bypass_water;
2314  wc_.accumulated_infiltration += bypass_water;
2315 
2316  if ( cbm::flt_greater_equal( water_def_total, bypass_water))
2317  {
2318 
2319  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2320  {
2321  double const water_def( cbm::bound_min( 0.0, (sc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]));
2322  double const water_add( water_def / water_def_total * bypass_water);
2323 
2324  wc_.wc_sl[sl] += water_add / sc_.h_sl[sl];
2325  }
2326  }
2327  else
2328  {
2329  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2330  {
2331  double const water_def( cbm::bound_min( 0.0, (sc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]));
2332  wc_.wc_sl[sl] += water_def / sc_.h_sl[sl];
2333  }
2334  /* add surplus to percolation out of last soil layer */
2335  wc_.accumulated_percolation += bypass_water - water_def_total;
2336  wc_.accumulated_waterflux_sl[m_soillayers->soil_layer_cnt()-1] += bypass_water - water_def_total;
2337  }
2338  }
2339 }
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2478

◆ 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

2202 {
2203  double const sks( WaterFlowSaturated( _sl, _fact_impedance));
2204  double const fsl( get_fsl( _sl));
2205  double const fhr( get_time_step_duration());
2206 
2208  double slope(( _sl < m_soillayers->soil_layers_in_litter_cnt()) ? m_param->SLOPE_FF() : m_param->SLOPE_MS());
2209 
2210  //max 0.005 means that no sks values up to 9.88 cm min-1 are allowed !!
2211  double const travelTime( std::max( 0.005, 1.0 - log10( wc_.sks_sl[_sl] * cbm::M_IN_CM * cbm::MIN_IN_DAY * fhr)));
2212  return( cbm::bound_max( cbm::sqr(1.0 - (sc_.wc_fc_sl[_sl] / ((wc_.wc_sl[_sl] + wc_.ice_sl[_sl]) * slope)))
2213  * 0.5 * wc_.wc_sl[_sl] * sc_.h_sl[_sl] * fsl * (1.0 - exp(-1.0 / travelTime)) * _fact_impedance,
2214  sks));
2215 }
double WaterFlowSaturated(size_t, double const &)
Definition: watercycle-dndc.cpp:2278
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2478
double get_fsl(size_t)
Returns factor accounting for soil layer depth. Originally, a constant soil layer depth of 0...
Definition: watercycle-dndc.cpp:2375

◆ WaterFlowSaturated()

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

Returns saturated water flow per time step accounting for impedance

2280 {
2281  double const sks( wc_.sks_sl[_sl] * cbm::M_IN_CM * cbm::MIN_IN_DAY * get_time_step_duration() * _fact_impedance);
2282  return sks;
2283 }
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2478

Member Data Documentation

◆ IMAX_W

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

number of iteration steps within a day