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.
 
double get_clay_limitation (size_t, double, double)
 
double rootdensityfortranspiration_sl (double, double, double, bool)
 
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.

2450 {
2451  double balance = get_leaf_water() + wc_.surface_ice + wc_.surface_water + plant_waterdeficit_compensation;
2452  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2453  {
2454  balance += (wc_.wc_sl[sl] + wc_.ice_sl[sl]) * sc_.h_sl[sl];
2455  }
2456  balance += ( - m_wc.accumulated_irrigation_event_executed
2457  - wc_.accumulated_irrigation_automatic
2458  - wc_.accumulated_irrigation_ggcmi
2459  - wc_.accumulated_irrigation_reservoir_withdrawal
2460  - wc_.accumulated_precipitation
2461  - wc_.accumulated_groundwater_access
2462  + wc_.accumulated_percolation
2463  + wc_.accumulated_runoff
2464  + wc_.accumulated_wateruptake_sl.sum()
2465  + wc_.accumulated_interceptionevaporation
2466  + wc_.accumulated_soilevaporation
2467  + wc_.accumulated_surfacewaterevaporation);
2468 
2469  if ( _balance)
2470  {
2471  double const balance_delta( std::abs( *_balance - balance));
2472  if ( balance_delta > cbm::DIFFMAX)
2473  {
2474  KLOGWARN( "Water leakage: difference is ", balance - *_balance);
2475  }
2476  }
2477 
2478  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2479  {
2480  if ( cbm::flt_less( wc_.wc_sl[sl], 0.0))
2481  {
2482  KLOGWARN( "Negative water content of: ", wc_.wc_sl[sl], " in layer: ", sl, ". Water content set to zero.");
2483  wc_.wc_sl[sl] = 0.0;
2484  }
2485  if ( cbm::flt_less( wc_.ice_sl[sl], 0.0))
2486  {
2487  KLOGWARN( "Negative ice content of: ", wc_.ice_sl[sl], " in layer: ", sl, ". Ice content set to zero.");
2488  wc_.ice_sl[sl] = 0.0;
2489  }
2490  }
2491 
2492  return balance;
2493 }
double get_leaf_water()
Collects leaf water from all canopy layers.
Definition: watercycle-dndc.cpp:1163

◆ 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

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

◆ 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().

1044 {
1045  day_potevapo = 0.0;
1046  if ( m_cfg.potentialevapotranspiration_method == "thornthwaite")
1047  {
1048  if ( lclock()->is_position( TMODE_PRE_YEARLY))
1049  {
1051  lclock()->days_in_year(),
1052  m_climate->annual_temperature_average(),
1053  m_climate->annual_temperature_amplitude());
1054  }
1055 
1056  double const daylength(ldndc::meteo::daylength(m_setup->latitude(), lclock()->yearday()));
1057  double const avg_dayl(12.0);
1058 
1059  day_potevapo = m_param->WCDNDC_INCREASE_POT_EVAPOTRANS()
1061  mc_.nd_airtemperature,
1062  daylength, avg_dayl, lclock()->days_in_month(),
1064  }
1065  else if ( (m_cfg.potentialevapotranspiration_method == "penman") ||
1066  (m_cfg.potentialevapotranspiration_method == "priestleytaylor"))
1067  {
1068  /* clock */
1069  cbm::sclock_t const & clk( this->lclock_ref());
1070 
1071  /* leaf area index */
1072  double lai( 0.0);
1073  for( PlantIterator vt = m_veg->begin(); vt != m_veg->end(); ++vt)
1074  {
1075  lai += (*vt)->lai();
1076  }
1077 
1078  /* vapour pressure [10^3:Pa] */
1079  double vps( 0.0);
1080  ldndc::meteo::vps( mc_.nd_airtemperature, &vps);
1081  double vp( vps - m_climate->vpd_day( clk));
1082 
1083  if( m_cfg.potentialevapotranspiration_method == "priestleytaylor")
1084  {
1085  day_potevapo = m_param->WCDNDC_INCREASE_POT_EVAPOTRANS()
1086  * ldndc::priestleytaylor( clk.yearday(),
1087  clk.days_in_year(),
1088  mc_.albedo,
1089  m_setup->latitude(),
1090  mc_.nd_airtemperature,
1091  mc_.nd_shortwaveradiation_in,
1092  vp,
1093  m_param->PT_ALPHA());
1094  }
1095  else if( m_cfg.potentialevapotranspiration_method == "penman")
1096  {
1097  ldndc::surface_type surface( cbm::flt_greater_zero( lai) ? short_grass : bare_soil);
1098 
1099  day_potevapo = m_param->WCDNDC_INCREASE_POT_EVAPOTRANS()
1100  * ldndc::penman( surface,
1101  clk.yearday(),
1102  clk.days_in_year(),
1103  mc_.albedo,
1104  m_setup->latitude(),
1105  mc_.nd_shortwaveradiation_in * cbm::SEC_IN_DAY,
1106  mc_.nd_airtemperature,
1107  mc_.nd_windspeed,
1108  vp);
1109  }
1110  }
1111  else
1112  {
1113  KLOGERROR("[BUG] huh :o how did you get here? ",
1114  "tell the maintainers refering to this error message");
1115  return LDNDC_ERR_FAIL;
1116  }
1117 
1118  if ( m_cfg.evapotranspiration_method == "uniform")
1119  {
1120  // hourly evaporation demand is calculated from daily demand equally throughout the day
1121  hr_pot_evapotrans = day_potevapo * nhr / cbm::HR_IN_DAY;
1122  }
1123  else if ( m_cfg.evapotranspiration_method == "previouspattern")
1124  {
1125 
1126  hr_pot_evapotrans = day_potevapo * pot_evapotranspiration_ts[lclock()->subday()-1] / nwc;
1127  }
1128  else //if ( m_cfg.evapotranspiration_method === "nice name") // could be for example 'daytime_only'
1129  {
1130  double const hour(lclock()->subday());
1131  double const day(lclock()->yearday());
1132  double const hsun(0.833 * cbm::PI / 180.0);
1133  double const lat(m_setup->latitude() * cbm::PI / 180.0);
1134  double const dec(-23.44 * (cbm::PI / 180.0) * std::cos(2.0 * cbm::PI * (day - 172.0) / lclock()->days_in_year()));
1135  double const help((std::sin(hsun) - std::sin(lat) * std::sin(dec)) / (std::cos(lat) * std::cos(dec)));
1136  double dayl(0.0);
1137  if (help < -0.999) dayl = 0.0;
1138  else if (help > 0.999) dayl = 24.0;
1139  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)));
1140  if (dayl < 1.0) dayl = 0.0;
1141  double const rise(12.0 - dayl * 0.5);
1142 
1143  double factor(0.0);
1144  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);
1145 
1146  double fday(0.0);
1147  if (hour > rise && hour < (24.0 - rise)) fday = factor;
1148 
1149  hr_pot_evapotrans = day_potevapo * fday;
1150  }
1151  // the evaporation limit for all timesteps is the daily demand
1152  wc_.accumulated_potentialevaporation += hr_pot_evapotrans;
1153 
1154  return LDNDC_ERR_OK;
1155 }
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:138
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
1830 {
1831  //Average limit of snowEvaporation is 0.00013m day-1 (Granger and Male 1978)
1832  double hr_potESnow( cbm::bound(0.0, hr_pot_evapotrans, 0.00013 * nhr/24.0));
1833 
1834  if ( cbm::flt_greater_zero( wc_.surface_ice) &&
1835  cbm::flt_greater_zero( hr_potESnow))
1836  {
1837  //If there is a snowpack and there is a potential to evaporate
1838  if ( wc_.surface_ice > hr_potESnow)
1839  {
1840  //Partial evaporation from the snow
1841  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - hr_potESnow);
1842  wc_.accumulated_surfacewaterevaporation += hr_potESnow;
1843  wc_.surface_ice -= hr_potESnow;
1844  }
1845  else
1846  {
1847  //Total evaporation of the snowpack
1848  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - wc_.surface_ice);;
1849  wc_.accumulated_surfacewaterevaporation += wc_.surface_ice;
1850  wc_.surface_ice = 0.0;
1851  }
1852  }
1853 }

◆ CalcSoilEvapoPercolation()

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

◆ CalcSurfaceFlux()

void ldndc::WatercycleDNDC::CalcSurfaceFlux ( )
private
Parameters
[in]None
[out]None
Returns
None
2364 {
2365  if ( cbm::flt_greater_zero( m_eventflood.get_bund_height()))
2366  {
2367  if ( cbm::flt_greater( wc_.surface_water, m_eventflood.get_bund_height()))
2368  {
2369  wc_.accumulated_runoff += wc_.surface_water - m_eventflood.get_bund_height();
2370  wc_.surface_water = m_eventflood.get_bund_height();
2371  }
2372  }
2373  else
2374  {
2375  /* runoff fraction per time step equals RO fraction per hour * length of a time step */
2376  double const runoff_fraction( m_param->FRUNOFF() * get_time_step_duration());
2377 
2378  // if frunoff_ts < 1
2379  if ( cbm::flt_less( runoff_fraction, 1.0))
2380  {
2381  double const runoff( runoff_fraction * wc_.surface_water);
2382  wc_.surface_water -= runoff;
2383  wc_.accumulated_runoff += runoff;
2384  }
2385  else
2386  {
2387  wc_.accumulated_runoff += wc_.surface_water;
2388  wc_.surface_water = 0.0;
2389  }
2390  }
2391 }
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2438

◆ CalcSurfaceWaterEvaporation()

void ldndc::WatercycleDNDC::CalcSurfaceWaterEvaporation ( )
private
Parameters
[in]None
[out]None
Returns
None
1799 {
1800  if ( cbm::flt_greater_zero( wc_.surface_water) &&
1801  cbm::flt_greater_zero( hr_pot_evapotrans))
1802  {
1803  if ( hr_pot_evapotrans > wc_.surface_water)
1804  {
1805  wc_.accumulated_surfacewaterevaporation += wc_.surface_water;
1806  hr_pot_evapotrans -= wc_.surface_water;
1807  wc_.surface_water = 0.0;
1808  }
1809  else
1810  {
1811  wc_.surface_water -= hr_pot_evapotrans;
1812  wc_.accumulated_surfacewaterevaporation += hr_pot_evapotrans;
1813  hr_pot_evapotrans = 0.0;
1814  }
1815  }
1816 }

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

◆ event_flood()

lerr_t ldndc::WatercycleDNDC::event_flood ( )
private

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

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

◆ get_clay_limitation()

double ldndc::WatercycleDNDC::get_clay_limitation ( size_t  _sl,
double  _water_avail,
double  _water_ref 
)
private
Parameters
[in]_slSoil layer
[in]_water_availAvailable water content
[in]_water_refReference water content
[out]None
Returns
Clay limitation
1315 {
1316  // layer specific water availability
1317  // organic fraction [%]
1318  double const orgf( sc_.fcorg_sl[_sl] / cbm::CCORG);
1319  // clay factor "1" [-]
1320  double const clay_f1( pow(10.0, (-sc_.clay_sl[_sl] - m_param->RCLAY())));
1321  // clay factor "2" [-]
1322  double const clay_f2( std::max(0.0, 1.0 + std::min(1.0 - orgf, sc_.clay_sl[_sl])));
1323 
1324  // fraction of max uptake that can be used according to soil conditions [%]
1325  return cbm::bound(0.0, clay_f1 + pow( _water_avail / _water_ref, clay_f2), 1.0);
1326 }

◆ 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
1164 {
1165  if ( m_veg->size() > 0u)
1166  {
1167  return wc_.wc_fl.sum();
1168  }
1169  return 0.0;
1170 }

◆ 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().

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

◆ rootdensityfortranspiration_sl()

double ldndc::WatercycleDNDC::rootdensityfortranspiration_sl ( double  m_sl,
double  l_sl,
double  h_sl,
bool  root_length_h2o_up 
)
private
Parameters
[in]m_slroot mass in soil layer
[in]l_slroot length in soil layer
[in]h_slheight of soil layer
[in]root_length_h2o_up,determinesif root mass or length is used
Returns
root density in this soil layer, either by mass or length
1335 {
1336  double rootdensity_sl = -99.99;
1337  if( root_length_h2o_up == true)
1338  {
1339  rootdensity_sl = l_sl / h_sl;
1340  }
1341  else
1342  {
1343  rootdensity_sl = m_sl / h_sl;
1344  }
1345 
1346  return rootdensity_sl;
1347 }

◆ SoilWaterEvaporation()

void ldndc::WatercycleDNDC::SoilWaterEvaporation ( size_t  _sl)
private
Parameters
[in]_slSoil layer
[out]None
Returns
None
1897 {
1898  //percolation and direct evaporation from the (upper) soil
1899  //if there is snow cover, evaporation from soil will be set to 0, and only sublimation is allowed.
1900  double const wc_min_evaporation( m_param->WCDNDC_EVALIM_FRAC_WCMIN() * sc_.wc_wp_sl[_sl]);
1901  if ( cbm::flt_less( sc_.depth_sl[_sl] - sc_.h_sl[_sl], m_param->EVALIM()) &&
1902  cbm::flt_greater( wc_.wc_sl[_sl], wc_min_evaporation) &&
1903  cbm::flt_greater_zero( hr_pot_evapotrans))
1904  {
1905  //limiting factor for soil water infiltration
1906  double const f_clay( 1.0 + m_param->SLOPE_CLAYF() * std::min(1.0 - sc_.fcorg_sl[_sl] / cbm::CCORG, sc_.clay_sl[_sl]));
1907  double const f_water( cbm::bound( 0.0,
1908  pow( (wc_.wc_sl[_sl] - wc_min_evaporation) / sc_.wc_fc_sl[_sl], f_clay),
1909  1.0));
1910  double const h_rest( std::min( sc_.h_sl[_sl], m_param->EVALIM() - ( sc_.depth_sl[_sl] - sc_.h_sl[_sl])));
1911  double const f_depth( h_rest / m_param->EVALIM() * std::max(0.0,
1912  1.0 - std::min((double)m_param->EVALIM(),
1913  sc_.depth_sl[_sl] - ( 0.5 * sc_.h_sl[_sl]))
1914  / m_param->EVALIM()));
1915  double const soilevaporation_max( cbm::bound_max( f_water * wc_.wc_sl[_sl] * sc_.h_sl[_sl], hr_pot_evapotrans));
1916 
1917  double const soilevaporation( cbm::bound( 0.0,
1918  hr_pot_evapotrans * f_depth * f_water,
1919  soilevaporation_max));
1920 
1921  //update of the current soil layer water content
1922  wc_.wc_sl[_sl] -= soilevaporation / sc_.h_sl[_sl];
1923  wc_.accumulated_soilevaporation += soilevaporation;
1924  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - soilevaporation);
1925  }
1926 }

◆ WatercycleDNDC_preferential_flow()

void ldndc::WatercycleDNDC::WatercycleDNDC_preferential_flow ( )
private
Parameters
[in]None
[out]None
Returns
None
2253 {
2254  //BY_PASSF is defined as the fraction of surface water that passes per hour (0.0 by default)
2255  if ( cbm::flt_greater_zero( m_param->BY_PASSF()) &&
2256  cbm::flt_greater_zero( wc_.surface_water))
2257  {
2258 
2259  double water_def_total( 0.0);
2260  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2261  {
2262  //water flux that flow through macro pores can maximally store
2263  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]);
2264  }
2265 
2266  //duration of timestep (hr-1)
2267  double const fhr( get_time_step_duration());
2268 
2269 
2270  //fraction of surface water that is put into bypass flow (m timestep-1)
2271  double const bypass_fraction( cbm::bound_max( m_param->BY_PASSF() * fhr, 0.99));
2272  double const bypass_water( wc_.surface_water * bypass_fraction);
2273  wc_.surface_water -= bypass_water;
2274  wc_.accumulated_infiltration += bypass_water;
2275 
2276  if ( cbm::flt_greater_equal( water_def_total, bypass_water))
2277  {
2278 
2279  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2280  {
2281  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]));
2282  double const water_add( water_def / water_def_total * bypass_water);
2283 
2284  wc_.wc_sl[sl] += water_add / sc_.h_sl[sl];
2285  }
2286  }
2287  else
2288  {
2289  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2290  {
2291  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]));
2292  wc_.wc_sl[sl] += water_def / sc_.h_sl[sl];
2293  }
2294  /* add surplus to percolation out of last soil layer */
2295  wc_.accumulated_percolation += bypass_water - water_def_total;
2296  wc_.accumulated_waterflux_sl[m_soillayers->soil_layer_cnt()-1] += bypass_water - water_def_total;
2297  }
2298  }
2299 }
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2438

◆ 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

2162 {
2163  double const sks( WaterFlowSaturated( _sl, _fact_impedance));
2164  double const fsl( get_fsl( _sl));
2165  double const fhr( get_time_step_duration());
2166 
2168  double slope(( _sl < m_soillayers->soil_layers_in_litter_cnt()) ? m_param->SLOPE_FF() : m_param->SLOPE_MS());
2169 
2170  //max 0.005 means that no sks values up to 9.88 cm min-1 are allowed !!
2171  double const travelTime( std::max( 0.005, 1.0 - log10( sc_.sks_sl[_sl] * cbm::M_IN_CM * cbm::MIN_IN_DAY * fhr)));
2172  return( cbm::bound_max( cbm::sqr(1.0 - (sc_.wc_fc_sl[_sl] / ((wc_.wc_sl[_sl] + wc_.ice_sl[_sl]) * slope)))
2173  * 0.5 * wc_.wc_sl[_sl] * sc_.h_sl[_sl] * fsl * (1.0 - exp(-1.0 / travelTime)) * _fact_impedance,
2174  sks));
2175 }
double WaterFlowSaturated(size_t, double const &)
Definition: watercycle-dndc.cpp:2238
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2438
double get_fsl(size_t)
Returns factor accounting for soil layer depth. Originally, a constant soil layer depth of 0...
Definition: watercycle-dndc.cpp:2335

◆ WaterFlowSaturated()

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

Returns saturated water flow per time step accounting for impedance

2240 {
2241  double const sks( sc_.sks_sl[_sl] * cbm::M_IN_CM * cbm::MIN_IN_DAY * get_time_step_duration() * _fact_impedance);
2242  return sks;
2243 }
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2438

Member Data Documentation

◆ IMAX_W

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

number of iteration steps within a day