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.

2477 {
2478  double balance = get_leaf_water() + wc_.surface_ice + wc_.surface_water;
2479  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2480  {
2481  balance += (wc_.wc_sl[sl] + wc_.ice_sl[sl]) * sc_.h_sl[sl];
2482  }
2483  balance += ( - m_wc.accumulated_irrigation_event_executed
2484  - wc_.accumulated_irrigation_automatic
2485  - wc_.accumulated_irrigation_ggcmi
2486  - wc_.accumulated_irrigation_reservoir_withdrawal
2487  - wc_.accumulated_precipitation
2488  - wc_.accumulated_groundwater_access
2489  + wc_.accumulated_percolation
2490  + wc_.accumulated_runoff
2491  + wc_.accumulated_wateruptake_sl.sum()
2492  + wc_.accumulated_interceptionevaporation
2493  + wc_.accumulated_soilevaporation
2494  + wc_.accumulated_surfacewaterevaporation);
2495 
2496  if ( _balance)
2497  {
2498  double const balance_delta( std::abs( *_balance - balance));
2499  if ( balance_delta > cbm::DIFFMAX)
2500  {
2501  KLOGWARN( "Water leakage: difference is ", balance - *_balance);
2502  }
2503  }
2504 
2505  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2506  {
2507  if ( cbm::flt_less( wc_.wc_sl[sl], 0.0))
2508  {
2509  KLOGWARN( "Negative water content of: ", wc_.wc_sl[sl], " in layer: ", sl, ". Water content set to zero.");
2510  wc_.wc_sl[sl] = 0.0;
2511  }
2512  if ( cbm::flt_less( wc_.ice_sl[sl], 0.0))
2513  {
2514  KLOGWARN( "Negative ice content of: ", wc_.ice_sl[sl], " in layer: ", sl, ". Ice content set to zero.");
2515  wc_.ice_sl[sl] = 0.0;
2516  }
2517  }
2518 
2519  return balance;
2520 }
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
1857 {
1858  //Average limit of snowEvaporation is 0.00013m day-1 (Granger and Male 1978)
1859  double hr_potESnow( cbm::bound(0.0, hr_pot_evapotrans, 0.00013 * nhr/24.0));
1860 
1861  if ( cbm::flt_greater_zero( wc_.surface_ice) &&
1862  cbm::flt_greater_zero( hr_potESnow))
1863  {
1864  //If there is a snowpack and there is a potential to evaporate
1865  if ( wc_.surface_ice > hr_potESnow)
1866  {
1867  //Partial evaporation from the snow
1868  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - hr_potESnow);
1869  wc_.accumulated_surfacewaterevaporation += hr_potESnow;
1870  wc_.surface_ice -= hr_potESnow;
1871  }
1872  else
1873  {
1874  //Total evaporation of the snowpack
1875  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - wc_.surface_ice);;
1876  wc_.accumulated_surfacewaterevaporation += wc_.surface_ice;
1877  wc_.surface_ice = 0.0;
1878  }
1879  }
1880 }

◆ CalcSoilEvapoPercolation()

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

◆ CalcSurfaceFlux()

void ldndc::WatercycleDNDC::CalcSurfaceFlux ( )
private
Parameters
[in]None
[out]None
Returns
None
2391 {
2392  if ( cbm::flt_greater_zero( m_eventflood.get_bund_height()))
2393  {
2394  if ( cbm::flt_greater( wc_.surface_water, m_eventflood.get_bund_height()))
2395  {
2396  wc_.accumulated_runoff += wc_.surface_water - m_eventflood.get_bund_height();
2397  wc_.surface_water = m_eventflood.get_bund_height();
2398  }
2399  }
2400  else
2401  {
2402  /* runoff fraction per time step equals RO fraction per hour * length of a time step */
2403  double const runoff_fraction( m_param->FRUNOFF() * get_time_step_duration());
2404 
2405  // if frunoff_ts < 1
2406  if ( cbm::flt_less( runoff_fraction, 1.0))
2407  {
2408  double const runoff( runoff_fraction * wc_.surface_water);
2409  wc_.surface_water -= runoff;
2410  wc_.accumulated_runoff += runoff;
2411  }
2412  else
2413  {
2414  wc_.accumulated_runoff += wc_.surface_water;
2415  wc_.surface_water = 0.0;
2416  }
2417  }
2418 }
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2465

◆ CalcSurfaceWaterEvaporation()

void ldndc::WatercycleDNDC::CalcSurfaceWaterEvaporation ( )
private
Parameters
[in]None
[out]None
Returns
None
1826 {
1827  if ( cbm::flt_greater_zero( wc_.surface_water) &&
1828  cbm::flt_greater_zero( hr_pot_evapotrans))
1829  {
1830  if ( hr_pot_evapotrans > wc_.surface_water)
1831  {
1832  wc_.accumulated_surfacewaterevaporation += wc_.surface_water;
1833  hr_pot_evapotrans -= wc_.surface_water;
1834  wc_.surface_water = 0.0;
1835  }
1836  else
1837  {
1838  wc_.surface_water -= hr_pot_evapotrans;
1839  wc_.accumulated_surfacewaterevaporation += hr_pot_evapotrans;
1840  hr_pot_evapotrans = 0.0;
1841  }
1842  }
1843 }

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

◆ 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:2265
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2465

◆ 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
1924 {
1925  //percolation and direct evaporation from the (upper) soil
1926  //if there is snow cover, evaporation from soil will be set to 0, and only sublimation is allowed.
1927  double const wc_min_evaporation( m_param->WCDNDC_EVALIM_FRAC_WCMIN() * sc_.wc_wp_sl[_sl]);
1928  if ( cbm::flt_less( sc_.depth_sl[_sl] - sc_.h_sl[_sl], m_param->EVALIM()) &&
1929  cbm::flt_greater( wc_.wc_sl[_sl], wc_min_evaporation) &&
1930  cbm::flt_greater_zero( hr_pot_evapotrans))
1931  {
1932  //limiting factor for soil water infiltration
1933  double const f_clay( 1.0 + m_param->SLOPE_CLAYF() * std::min(1.0 - sc_.fcorg_sl[_sl] / cbm::CCORG, sc_.clay_sl[_sl]));
1934  double const f_water( cbm::bound( 0.0,
1935  pow( (wc_.wc_sl[_sl] - wc_min_evaporation) / sc_.wc_fc_sl[_sl], f_clay),
1936  1.0));
1937  double const h_rest( std::min( sc_.h_sl[_sl], m_param->EVALIM() - ( sc_.depth_sl[_sl] - sc_.h_sl[_sl])));
1938  double const f_depth( h_rest / m_param->EVALIM() * std::max(0.0,
1939  1.0 - std::min((double)m_param->EVALIM(),
1940  sc_.depth_sl[_sl] - ( 0.5 * sc_.h_sl[_sl]))
1941  / m_param->EVALIM()));
1942  double const soilevaporation_max( cbm::bound_max( f_water * wc_.wc_sl[_sl] * sc_.h_sl[_sl], hr_pot_evapotrans));
1943 
1944  double const soilevaporation( cbm::bound( 0.0,
1945  hr_pot_evapotrans * f_depth * f_water,
1946  soilevaporation_max));
1947 
1948  //update of the current soil layer water content
1949  wc_.wc_sl[_sl] -= soilevaporation / sc_.h_sl[_sl];
1950  wc_.accumulated_soilevaporation += soilevaporation;
1951  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - soilevaporation);
1952  }
1953 }

◆ WatercycleDNDC_preferential_flow()

void ldndc::WatercycleDNDC::WatercycleDNDC_preferential_flow ( )
private
Parameters
[in]None
[out]None
Returns
None
2280 {
2281  //BY_PASSF is defined as the fraction of surface water that passes per hour (0.0 by default)
2282  if ( cbm::flt_greater_zero( m_param->BY_PASSF()) &&
2283  cbm::flt_greater_zero( wc_.surface_water))
2284  {
2285 
2286  double water_def_total( 0.0);
2287  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2288  {
2289  //water flux that flow through macro pores can maximally store
2290  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]);
2291  }
2292 
2293  //duration of timestep (hr-1)
2294  double const fhr( get_time_step_duration());
2295 
2296 
2297  //fraction of surface water that is put into bypass flow (m timestep-1)
2298  double const bypass_fraction( cbm::bound_max( m_param->BY_PASSF() * fhr, 0.99));
2299  double const bypass_water( wc_.surface_water * bypass_fraction);
2300  wc_.surface_water -= bypass_water;
2301  wc_.accumulated_infiltration += bypass_water;
2302 
2303  if ( cbm::flt_greater_equal( water_def_total, bypass_water))
2304  {
2305 
2306  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2307  {
2308  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]));
2309  double const water_add( water_def / water_def_total * bypass_water);
2310 
2311  wc_.wc_sl[sl] += water_add / sc_.h_sl[sl];
2312  }
2313  }
2314  else
2315  {
2316  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2317  {
2318  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]));
2319  wc_.wc_sl[sl] += water_def / sc_.h_sl[sl];
2320  }
2321  /* add surplus to percolation out of last soil layer */
2322  wc_.accumulated_percolation += bypass_water - water_def_total;
2323  wc_.accumulated_waterflux_sl[m_soillayers->soil_layer_cnt()-1] += bypass_water - water_def_total;
2324  }
2325  }
2326 }
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2465

◆ 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

2189 {
2190  double const sks( WaterFlowSaturated( _sl, _fact_impedance));
2191  double const fsl( get_fsl( _sl));
2192  double const fhr( get_time_step_duration());
2193 
2195  double slope(( _sl < m_soillayers->soil_layers_in_litter_cnt()) ? m_param->SLOPE_FF() : m_param->SLOPE_MS());
2196 
2197  //max 0.005 means that no sks values up to 9.88 cm min-1 are allowed !!
2198  double const travelTime( std::max( 0.005, 1.0 - log10( wc_.sks_sl[_sl] * cbm::M_IN_CM * cbm::MIN_IN_DAY * fhr)));
2199  return( cbm::bound_max( cbm::sqr(1.0 - (sc_.wc_fc_sl[_sl] / ((wc_.wc_sl[_sl] + wc_.ice_sl[_sl]) * slope)))
2200  * 0.5 * wc_.wc_sl[_sl] * sc_.h_sl[_sl] * fsl * (1.0 - exp(-1.0 / travelTime)) * _fact_impedance,
2201  sks));
2202 }
double WaterFlowSaturated(size_t, double const &)
Definition: watercycle-dndc.cpp:2265
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2465
double get_fsl(size_t)
Returns factor accounting for soil layer depth. Originally, a constant soil layer depth of 0...
Definition: watercycle-dndc.cpp:2362

◆ WaterFlowSaturated()

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

Returns saturated water flow per time step accounting for impedance

2267 {
2268  double const sks( wc_.sks_sl[_sl] * cbm::M_IN_CM * cbm::MIN_IN_DAY * get_time_step_duration() * _fact_impedance);
2269  return sks;
2270 }
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2465

Member Data Documentation

◆ IMAX_W

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

number of iteration steps within a day