LandscapeDNDC  1.37.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 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.

2577 {
2578  double balance = get_leaf_water() + wc_.surface_ice + wc_.surface_water;
2579  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2580  {
2581  balance += (wc_.wc_sl[sl] + wc_.ice_sl[sl]) * sc_.h_sl[sl];
2582  }
2583  balance += ( - m_wc.accumulated_irrigation_event_executed
2584  - wc_.accumulated_irrigation_automatic
2585  - wc_.accumulated_irrigation_ggcmi
2586  - wc_.accumulated_irrigation_reservoir_withdrawal
2587  - wc_.accumulated_precipitation
2588  - wc_.accumulated_groundwater_access
2589  + wc_.accumulated_percolation
2590  + wc_.accumulated_runoff
2591  + wc_.accumulated_wateruptake_sl.sum()
2592  + wc_.accumulated_interceptionevaporation
2593  + wc_.accumulated_soilevaporation
2594  + wc_.accumulated_surfacewaterevaporation);
2595 
2596  if ( _balance)
2597  {
2598  double const balance_delta( std::abs( *_balance - balance));
2599  if ( balance_delta > cbm::DIFFMAX)
2600  {
2601  KLOGWARN( "Water leakage: difference is ", balance - *_balance);
2602  }
2603  }
2604 
2605  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2606  {
2607  if ( cbm::flt_less( wc_.wc_sl[sl], 0.0))
2608  {
2609  KLOGWARN( "Negative water content of: ", wc_.wc_sl[sl], " in layer: ", sl, ". Water content set to zero.");
2610  wc_.wc_sl[sl] = 0.0;
2611  }
2612  if ( cbm::flt_less( wc_.ice_sl[sl], 0.0))
2613  {
2614  KLOGWARN( "Negative ice content of: ", wc_.ice_sl[sl], " in layer: ", sl, ". Ice content set to zero.");
2615  wc_.ice_sl[sl] = 0.0;
2616  }
2617  }
2618 
2619  return balance;
2620 }
double get_leaf_water()
Collects leaf water from all canopy layers.
Definition: watercycle-dndc.cpp:1117

◆ 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

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

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

550 {
551  // Daily evaportation demand is derived from the Thornthwaite, the Penman, or the Priestley tayler method (standard is 'thornthwaite).
552  day_potevapo = 0.0;
553  if ( m_cfg.potentialevapotranspiration_method == "thornthwaite")
554  {
555  if ( lclock()->is_position( TMODE_PRE_YEARLY))
556  {
558  lclock()->days_in_year(),
559  m_climate->annual_temperature_average(),
560  m_climate->annual_temperature_amplitude());
561  }
562 
563  double const daylength(ldndc::meteo::daylength(m_setup->latitude(), lclock()->yearday()));
564  double const avg_dayl(12.0);
565 
566  day_potevapo = m_sipar->WCDNDC_INCREASE_POT_EVAPOTRANS()
568  mc_.nd_airtemperature,
569  daylength, avg_dayl, lclock()->days_in_month(),
571  }
572  else if ( (m_cfg.potentialevapotranspiration_method == "penman") ||
573  (m_cfg.potentialevapotranspiration_method == "priestleytaylor"))
574  {
575  cbm::sclock_t const & clk( lclock_ref());
576 
577  /* vapour pressure [10^3:Pa] */
578  double vps( 0.0);
579  ldndc::meteo::vps( mc_.nd_airtemperature, &vps);
580  double vp( vps - m_climate->vpd_day( clk));
581 
582  if( m_cfg.potentialevapotranspiration_method == "priestleytaylor")
583  {
584  day_potevapo = m_sipar->WCDNDC_INCREASE_POT_EVAPOTRANS()
585  * ldndc::priestleytaylor( clk.yearday(),
586  clk.days_in_year(),
587  mc_.albedo,
588  m_setup->latitude(),
589  mc_.nd_airtemperature,
590  mc_.nd_shortwaveradiation_in,
591  vp,
592  m_sipar->PT_ALPHA());
593  }
594  else if( m_cfg.potentialevapotranspiration_method == "penman")
595  {
596  /* leaf area index */
597  double const lai( m_veg->lai());
598  ldndc::surface_type surface( cbm::flt_greater_zero( lai) ? short_grass : bare_soil);
599 
600  day_potevapo = m_sipar->WCDNDC_INCREASE_POT_EVAPOTRANS()
601  * ldndc::penman( surface,
602  clk.yearday(),
603  clk.days_in_year(),
604  mc_.albedo,
605  m_setup->latitude(),
606  mc_.nd_shortwaveradiation_in * cbm::SEC_IN_DAY,
607  mc_.nd_airtemperature,
608  mc_.nd_windspeed,
609  vp);
610  }
611  }
612  else
613  {
614  KLOGERROR("[BUG] huh :o how did you get here? ",
615  "tell the maintainers refering to this error message");
616  return LDNDC_ERR_FAIL;
617  }
618 
619  // The daily evaporation demand is distributed throughout the day (standard is 'previouspattern)
620  // TODO: design a routine that distributes evaporation demand according to temperature or humidity
621  if ( m_cfg.evapotranspiration_method == "uniform")
622  {
623  // hourly evaporation demand is calculated from daily demand equally throughout the day
624  hr_pot_evapotrans = day_potevapo * nhr / cbm::HR_IN_DAY;
625  }
626  else if ( m_cfg.evapotranspiration_method == "previouspattern")
627  {
628  // hourly evaporation demand is derived from the weighted expression of previous values
629  hr_pot_evapotrans = day_potevapo * pot_evapotranspiration_ts[lclock()->subday()-1] / nwc;
630  }
631  else
632  {
633  // hourly evaporation demand is distributed throughout the sunlight hours only (similar as radiation)
634  double const hour(lclock()->subday());
635  double const day(lclock()->yearday());
636  double const hsun(0.833 * cbm::PI / 180.0);
637  double const lat(m_setup->latitude() * cbm::PI / 180.0);
638  double const dec(-23.44 * (cbm::PI / 180.0) * std::cos(2.0 * cbm::PI * (day - 172.0) / lclock()->days_in_year()));
639  double const help((std::sin(hsun) - std::sin(lat) * std::sin(dec)) / (std::cos(lat) * std::cos(dec)));
640  double dayl(0.0);
641  if (help < -0.999) dayl = 0.0;
642  else if (help > 0.999) dayl = 24.0;
643  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)));
644  if (dayl < 1.0) dayl = 0.0;
645  double const rise(12.0 - dayl * 0.5);
646 
647  double factor(0.0);
648  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);
649 
650  double fday(0.0);
651  if (hour > rise && hour < (24.0 - rise)) fday = factor;
652 
653  hr_pot_evapotrans = day_potevapo * fday;
654  }
655 
656  // the evaporation limit for all timesteps is the daily demand
657  wc_.accumulated_potentialevaporation += hr_pot_evapotrans;
658 
659  return LDNDC_ERR_OK;
660 }
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:134
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
1940 {
1941  //Average limit of snowEvaporation is 0.00013m day-1 (Granger and Male 1978)
1942  double hr_potESnow( cbm::bound(0.0, hr_pot_evapotrans, 0.00013 * nhr/24.0));
1943 
1944  if ( cbm::flt_greater_zero( wc_.surface_ice) &&
1945  cbm::flt_greater_zero( hr_potESnow))
1946  {
1947  //If there is a snowpack and there is a potential to evaporate
1948  if ( wc_.surface_ice > hr_potESnow)
1949  {
1950  //Partial evaporation from the snow
1951  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - hr_potESnow);
1952  wc_.accumulated_surfacewaterevaporation += hr_potESnow;
1953  wc_.surface_ice -= hr_potESnow;
1954  }
1955  else
1956  {
1957  //Total evaporation of the snowpack
1958  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - wc_.surface_ice);;
1959  wc_.accumulated_surfacewaterevaporation += wc_.surface_ice;
1960  wc_.surface_ice = 0.0;
1961  }
1962  }
1963 }

◆ CalcSoilEvapoPercolation()

void ldndc::WatercycleDNDC::CalcSoilEvapoPercolation ( )
private
Parameters
[in]None
[out]None
Returns
None
2077 {
2078  // reset all water fluxes
2079  for (size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2080  {
2081  m_wc.percolation_sl[sl] = 0.0;
2082  }
2083 
2084  // fill groundwater layer with water
2085  for (int sl = m_soillayers->soil_layer_cnt() - 1; sl >= 0; --sl)
2086  {
2087  if (cbm::flt_greater_equal(sc_.depth_sl[sl], ground_water_table))
2088  {
2089  // ice volume is excepted from groundwater
2090  double const ice_vol(wc_.ice_sl[sl] / cbm::DICE);
2091  if (cbm::flt_greater(sc_.poro_sl[sl], ice_vol))
2092  {
2093  // new water in the layer is generated from groundwater
2094  double const gw_access(cbm::bound_min(0.0, (sc_.poro_sl[sl] - ice_vol - wc_.wc_sl[sl]) * sc_.h_sl[sl] ));
2095  wc_.wc_sl[sl] += gw_access / sc_.h_sl[sl];
2096  wc_.accumulated_groundwater_access += gw_access;
2097  }
2098  }
2099  else
2100  {
2101  break;
2102  }
2103  }
2104 
2105  for (size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2106  {
2107  // the first soil (litter) layer is filled with surface water up to its capacity
2108  if (sl == 0)
2109  {
2110  double delta_wc(0.0);
2111  double const wc_old(wc_.wc_sl[sl]);
2112  double const wc_cap( (sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl]); // uptake capacity in first soil layer [mm]
2113 
2114  if (cbm::flt_greater(wc_.surface_water, wc_cap))
2115  {
2116  wc_.wc_sl[sl] = sc_.poro_sl[sl];
2117  delta_wc = cbm::bound_max(wc_.surface_water, (cbm::bound_min(0.0, wc_.wc_sl[sl] - wc_old) * sc_.h_sl[sl]));
2118  wc_.surface_water -= delta_wc;
2119  }
2120  else
2121  {
2122  wc_.wc_sl[sl] = (wc_.wc_sl[sl] * sc_.h_sl[sl] + wc_.surface_water) / sc_.h_sl[sl];
2123  delta_wc = cbm::bound_max(wc_.surface_water, (cbm::bound_min(0.0, wc_.wc_sl[sl] - wc_old) * sc_.h_sl[sl]));
2124  wc_.surface_water = 0.0;
2125  }
2126 
2127  wc_.accumulated_infiltration += delta_wc;
2128  }
2129  // other soil layers get updated by percolation rates from the layer above
2130  else
2131  {
2132  wc_.wc_sl[sl] += m_wc.percolation_sl[sl - 1] / sc_.h_sl[sl];
2133  }
2134 
2135  // percolation for all layers below groundwater
2136  if (cbm::flt_greater(sc_.depth_sl[sl], ground_water_table))
2137  {
2138 
2139  // in case groundwater reaches the soil surface, the first outflux is estimated by saturated water flow
2140  if (sl == 0)
2141  {
2142  m_wc.percolation_sl[sl] = WaterFlowSaturated(sl, get_impedance_factor(sl));
2143  }
2144  // other soil layers are assumed to percolate by the same rate as the last non-groundwater layer
2145  else
2146  {
2147  m_wc.percolation_sl[sl] = m_wc.percolation_sl[sl - 1];
2148  }
2149 
2150  }
2151  // in all other cases percolation
2152  else if (sl + 1 == m_soillayers->soil_layer_cnt())
2153  {
2154  m_wc.percolation_sl[sl] = 0.0;
2155  if (cbm::flt_greater_zero(wc_.sks_sl[sl]))
2156  {
2157  if (cbm::flt_greater(wc_.wc_sl[sl], sc_.poro_sl[sl]) &&
2158  cbm::flt_greater(wc_.wc_sl[sl], wc_.wc_fc_sl[sl])) //security, should be always true if first condition is met
2159  {
2160  //maximum allowed flow until field capacity
2161  double const max_flow((wc_.wc_sl[sl] - wc_.wc_fc_sl[sl]) * sc_.h_sl[sl]);
2162  m_wc.percolation_sl[sl] = cbm::bound_max(WaterFlowSaturated(sl, get_impedance_factor(sl)), max_flow);
2163  }
2164 
2165  //water flux above field capacity
2166  else if (cbm::flt_greater_equal(wc_.wc_sl[sl] + wc_.ice_sl[sl], wc_.wc_fc_sl[sl]))
2167  {
2168  m_wc.percolation_sl[sl] = WaterFlowAboveFieldCapacity(sl, get_impedance_factor(sl));
2169  }
2170 
2171  //water flux below field capacity and above wilting point
2172  else if (cbm::flt_greater_zero(m_wc.percolation_sl[sl - 1]) &&
2173  cbm::flt_greater(wc_.wc_sl[sl], wc_.wc_wp_sl[sl]))
2174  {
2175  m_wc.percolation_sl[sl] = m_sipar->FPERCOL() * m_wc.percolation_sl[sl - 1];
2176  }
2177 
2178  //bound water flux by defined maximal percolation rate (e.g., during flooding events)
2179  m_wc.percolation_sl[sl] = cbm::bound_max(m_wc.percolation_sl[sl], max_percolation);
2180  }
2181  }
2182 
2183  //infiltrating water in current time step is exceeding available pore space of last time step
2184  //second condition: security, should be always true if first condition is met
2185  else if (cbm::flt_greater(wc_.wc_sl[sl], sc_.poro_sl[sl]) &&
2186  cbm::flt_greater(wc_.wc_sl[sl], wc_.wc_fc_sl[sl]))
2187  {
2188  //maximum allowed flow until field capacity
2189  double const max_flow((wc_.wc_sl[sl] - wc_.wc_fc_sl[sl]) * sc_.h_sl[sl]);
2190  m_wc.percolation_sl[sl] = cbm::bound_max(WaterFlowSaturated(sl, get_impedance_factor(sl)),
2191  max_flow);
2192  }
2193  else if (cbm::flt_greater_equal(wc_.wc_sl[sl] + wc_.ice_sl[sl], wc_.wc_fc_sl[sl]))
2194  {
2195  m_wc.percolation_sl[sl] = WaterFlowAboveFieldCapacity(sl, get_impedance_factor(sl));
2196  }
2197  else if (cbm::flt_greater(wc_.wc_sl[sl], wc_.wc_wp_sl[sl]))
2198  {
2199  m_wc.percolation_sl[sl] = WaterFlowCapillaryRise(sl, get_impedance_factor(sl));
2200  }
2201  else
2202  {
2203  m_wc.percolation_sl[sl] = 0.0;
2204  }
2205 
2206 
2207  // water content must be greater or equal wilting point
2208  if (cbm::flt_less((wc_.wc_sl[sl] + wc_.ice_sl[sl]) * sc_.h_sl[sl] - m_wc.percolation_sl[sl], wc_.wc_wp_sl[sl] * sc_.h_sl[sl]) &&
2209  cbm::flt_greater_zero(m_wc.percolation_sl[sl]))
2210  {
2211  m_wc.percolation_sl[sl] = cbm::bound(0.0,
2212  (wc_.wc_sl[sl] + wc_.ice_sl[sl] - wc_.wc_wp_sl[sl]) * sc_.h_sl[sl],
2213  m_wc.percolation_sl[sl]);
2214  }
2215 
2216  // Update of the current soil layer water content
2217  wc_.wc_sl[sl] -= m_wc.percolation_sl[sl] / sc_.h_sl[sl];
2218 
2220  }
2221 
2222  /* Correct water content and percolation if water content is above porosity
2223  * If the layer below cannot take up all of the water assigned for percolation,
2224  * so if the air filled pore space is not sufficient,
2225  * the water remains in the upper layer.
2226  */
2227  for (size_t sl = m_soillayers->soil_layer_cnt() - 1; sl > 0; sl--)
2228  {
2229  double const ice_vol(wc_.ice_sl[sl] / cbm::DICE);
2230  if (cbm::flt_greater(wc_.wc_sl[sl] + ice_vol, sc_.poro_sl[sl]))
2231  {
2232  double delta_wc(wc_.wc_sl[sl] + ice_vol - sc_.poro_sl[sl]);
2233  if (cbm::flt_greater_equal(ice_vol, sc_.poro_sl[sl]))
2234  {
2235  delta_wc = wc_.wc_sl[sl];
2236  wc_.wc_sl[sl] = 0.0;
2237  }
2238  else
2239  {
2240  wc_.wc_sl[sl] = sc_.poro_sl[sl] - ice_vol;
2241  }
2242  m_wc.percolation_sl[sl - 1] -= delta_wc * sc_.h_sl[sl];
2243  wc_.wc_sl[sl - 1] += delta_wc * sc_.h_sl[sl] / sc_.h_sl[sl - 1];
2244  }
2245  }
2246 
2247  double const ice_vol(wc_.ice_sl[0] / cbm::DICE);
2248  if (cbm::flt_greater(wc_.wc_sl[0] + ice_vol, sc_.poro_sl[0]))
2249  {
2250 
2251  double delta_wc(wc_.wc_sl[0] + ice_vol - sc_.poro_sl[0]);
2252  if (cbm::flt_greater_equal(ice_vol, sc_.poro_sl[0]))
2253  {
2254  delta_wc = wc_.wc_sl[0];
2255  wc_.wc_sl[0] = 0.0;
2256  }
2257  else
2258  {
2259  wc_.wc_sl[0] = sc_.poro_sl[0] - ice_vol;
2260  }
2261  wc_.accumulated_infiltration -= delta_wc * sc_.h_sl[0];
2262  wc_.surface_water += delta_wc * sc_.h_sl[0];
2263  }
2264 
2265  if (cbm::flt_greater_zero(wc_.surface_water))
2266  {
2267  CalcSurfaceFlux();
2268  }
2269 
2270  wc_.accumulated_waterflux_sl += m_wc.percolation_sl;
2271  wc_.accumulated_percolation += m_wc.percolation_sl[ m_soillayers->soil_layer_cnt()-1];
2272 }
double get_impedance_factor(size_t)
Reduces water fluxes out of a layer if ice lenses have formed.
Definition: watercycle-dndc.cpp:2440
void SoilWaterEvaporation(size_t)
Definition: watercycle-dndc.cpp:2005
double WaterFlowAboveFieldCapacity(size_t, double const &)
Definition: watercycle-dndc.cpp:2285
void CalcSurfaceFlux()
Definition: watercycle-dndc.cpp:2489
double WaterFlowCapillaryRise(size_t, double const &)
Similar to WaterFlowBelowFieldCapacity(...) but restricts water flow by water availability from soil ...
Definition: watercycle-dndc.cpp:2334

◆ CalcSurfaceFlux()

void ldndc::WatercycleDNDC::CalcSurfaceFlux ( )
private
Parameters
[in]None
[out]None
Returns
None
2490 {
2491  if ( cbm::flt_greater_zero( m_eventflood.get_bund_height()))
2492  {
2493  if ( cbm::flt_greater( wc_.surface_water, m_eventflood.get_bund_height()))
2494  {
2495  wc_.accumulated_runoff += wc_.surface_water - m_eventflood.get_bund_height();
2496  wc_.surface_water = m_eventflood.get_bund_height();
2497  }
2498  }
2499  else
2500  {
2501  /* runoff fraction per time step equals RO fraction per hour * length of a time step */
2502  double const runoff_fraction(m_sipar->FRUNOFF() * get_time_step_duration());
2503 
2504  // - not all surface water is removed within one timestep
2505  if ( cbm::flt_less( runoff_fraction, 1.0))
2506  {
2507  double const runoff( runoff_fraction * wc_.surface_water);
2508  wc_.surface_water -= runoff;
2509  wc_.accumulated_runoff += runoff;
2510  }
2511  // - all surface water is immediatly removed
2512  else
2513  {
2514  wc_.accumulated_runoff += wc_.surface_water;
2515  wc_.surface_water = 0.0;
2516  }
2517  }
2518 }
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2565

◆ CalcSurfaceWaterEvaporation()

void ldndc::WatercycleDNDC::CalcSurfaceWaterEvaporation ( )
private
Parameters
[in]None
[out]None
Returns
None
1909 {
1910  if ( cbm::flt_greater_zero( wc_.surface_water) &&
1911  cbm::flt_greater_zero( hr_pot_evapotrans))
1912  {
1913  if ( hr_pot_evapotrans > wc_.surface_water)
1914  {
1915  wc_.accumulated_surfacewaterevaporation += wc_.surface_water;
1916  hr_pot_evapotrans -= wc_.surface_water;
1917  wc_.surface_water = 0.0;
1918  }
1919  else
1920  {
1921  wc_.surface_water -= hr_pot_evapotrans;
1922  wc_.accumulated_surfacewaterevaporation += hr_pot_evapotrans;
1923  hr_pot_evapotrans = 0.0;
1924  }
1925  }
1926 }

◆ 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
1764 {
1765  // hourly potential transpiration
1766  double const hr_potTrans( hr_potTransinm * cbm::CM_IN_M);
1767 
1768  // hourly water uptake
1769  double hr_uptake( 0.0);
1770 
1771  // Calculate the root water uptake and hence the transpiration for every plant individually
1772  for( PlantIterator vt = m_veg->begin(); vt != m_veg->end(); ++vt)
1773  {
1774  // split up the potential transpiration equally to the various plants
1775  // frt mass of the current species
1776  double const mfrt_species( (*vt)->mFrt);
1777  // total frt mass
1778  double const mfrt_sum( m_veg->dw_frt());
1779  double const hr_potTransspecies = hr_potTrans * mfrt_species/mfrt_sum; // in cm
1780 
1781  // If some soil layers are frozen and contain only ice, no fluid water is present and the water potential diverges.
1782  // The uptake from these layers is prevented by considering only roots in layers with more water than ice.
1783  // In this case the root distribution needs to be rescaled (otherwise the non-contributing layers would effectively contribute a potential = 0).
1784  double fFrt_unfrozen( 1.0);
1785 
1786  /* TODO: include icetowatercrit as proper parameter */
1787  double icetowatercrit( 0.1);
1788  double effsoilwatpot( 0.0);
1789  for( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
1790  {
1791  // fine root conditions to prevent transpiration without uptake organs, wateruptake from frozen ground
1792  if( cbm::flt_greater_zero( (*vt)->fFrt_sl[sl] * mfrt_species))
1793  {
1794  if(cbm::flt_less( wc_.ice_sl[sl]/wc_.wc_sl[sl], icetowatercrit))
1795  {
1796  // overall water potential in this layer
1797  double const watpotl( Soillayerwaterhead( sl)); // negative
1798  effsoilwatpot += watpotl * (*vt)->fFrt_sl[sl]; // negative
1799  }
1800  else
1801  {
1802  fFrt_unfrozen -= (*vt)->fFrt_sl[sl];
1803  }
1804  }
1805  else
1806  {
1807  // roots are not wireless ;-)
1808  break;
1809  }
1810  }
1811  // normalize the root distribution in the potential, if some soil is frozen
1812  effsoilwatpot = effsoilwatpot / fFrt_unfrozen;
1813 
1814  // hydraulic conductivities of plant system
1815  double const Krs(vt->KRSINIT() * mfrt_species / vt->FRTMASSINIT() * nhr); // cm(water)/cm(from pressure)/timestep
1816  double const Kcomp(vt->KCOMPINIT() * mfrt_species / vt->FRTMASSINIT() * nhr); // cm(water)/cm(from pressure)/timestep
1817 
1818  // water potential in the leafs, depending on the potential transpiration
1819  double const leafwatpot( cbm::bound_min(vt->HLEAFCRIT(), effsoilwatpot - (hr_potTransspecies / Krs))); // cm
1820 
1821  if( cbm::flt_greater_zero( leafwatpot))
1822  {
1823  KLOGERROR( "leafwatpot is positive");
1824  return LDNDC_ERR_FAIL;
1825  }
1826  if( !cbm::flt_greater_zero( effsoilwatpot - leafwatpot))
1827  {
1828  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.");
1829  return LDNDC_ERR_FAIL;
1830  }
1831 
1832  double const effuptake( Krs * (effsoilwatpot - leafwatpot)); // cm/timestep
1833  double uptake_tot( 0.0);
1834  double uptakecomp_tot( 0.0);
1835 // double uptakecompabs_tot( 0.0);
1836 
1837  for( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
1838  {
1839  // fine root condition to prevent transpiration without uptake organs
1840  if( cbm::flt_greater_zero( (*vt)->fFrt_sl[sl] * mfrt_species))
1841  {
1842  if(cbm::flt_less( wc_.ice_sl[sl]/wc_.wc_sl[sl], icetowatercrit))
1843  {
1844  double const watpotl( Soillayerwaterhead( sl));
1845  double const compuptake( Kcomp * (watpotl - effsoilwatpot));
1846  double uptake_layer( (effuptake + compuptake) * (*vt)->fFrt_sl[sl]); // absolute value in cm
1847 
1848  wc_.accumulated_wateruptake_sl[sl] += uptake_layer * cbm::M_IN_CM;
1849  uptake_tot += uptake_layer; // absolute value in cm
1850  uptakecomp_tot += compuptake * (*vt)->fFrt_sl[sl]; // absolute value in cm
1851 // uptakecompabs_tot += fabs(compuptake * (*vt)->fFrt_sl[sl]); // absolute value in cm
1852 
1853  wc_.wc_sl[sl] -= uptake_layer * cbm::M_IN_CM /sc_.h_sl[sl]; // fraction of water in this layer, in m
1854 
1855  if( !cbm::flt_greater_zero( wc_.wc_sl[sl]))
1856  {
1857  KLOGERROR( "In soil layer ", sl, " the water content is negative: ", wc_.wc_sl[sl]);
1858  return LDNDC_ERR_FAIL;
1859  }
1860  }
1861  }
1862  else
1863  {
1864  // roots are not wireless ;-)
1865  break;
1866  }
1867  }
1868 
1869  hr_uptake += uptake_tot * cbm::M_IN_CM; // absolute value in m for all species
1870 
1871 
1872  if( cbm::flt_greater_zero( fabs(uptakecomp_tot))){
1873  KLOGWARN( "Compensatory water uptake ", uptakecomp_tot, " > 0.");
1874  }
1875  if( cbm::flt_greater( uptake_tot, hr_potTransspecies, 0.0000000001))
1876  {
1877  KLOGERROR( "Total water uptake ", uptake_tot, " larger than potential transpiration ", hr_potTransspecies);
1878  }
1879  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))
1880  {
1881  KLOGERROR( "T ", uptake_tot, "cm < Tpot ", hr_potTransspecies, "cm even though leafwatpot ", leafwatpot, " > HLEAFCRIT", vt->HLEAFCRIT());
1882  }
1883  }
1884 
1885  // total evaporation demand reduced to apply on other evaporating sources (ground, snow, surface water)
1886  hr_pot_evapotrans = cbm::bound_min(0.0, hr_pot_evapotrans - hr_uptake);
1887 
1888  return LDNDC_ERR_OK;
1889 }

◆ event_flood()

lerr_t ldndc::WatercycleDNDC::event_flood ( )
private

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

  • surface water table
  • bund height
782 {
783  lerr_t rc = m_eventflood.solve();
784  if ( rc != LDNDC_ERR_OK)
785  {
786  KLOGERROR("Irrigation event not successful!");
787  return rc;
788  }
789 
790  /* max_percolation is gradually decreased after end of flooding event to inital value */
791  double const maximum_percolation = m_eventflood.get_maximum_percolation();
792  if ( cbm::is_valid( maximum_percolation))
793  {
794  if ( cbm::flt_greater_zero( maximum_percolation))
795  {
796  max_percolation = maximum_percolation * nhr / 24.0;
797  }
798  }
799  else
800  {
801  /* time rate for gradual recovery of max_percolation */
802  double const time_rate( get_time_step_duration() * 0.1);
803  max_percolation -= (max_percolation - WaterFlowSaturated( m_soillayers->soil_layer_cnt()-1, 1.0)) * time_rate;
804  }
805 
806  return LDNDC_ERR_OK;
807 }
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2565

◆ 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
1118 {
1119  if ( m_veg->size() > 0u)
1120  {
1121  return wc_.wc_fl.sum();
1122  }
1123  return 0.0;
1124 }

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

117 {
118  if ( m_iokcomm->get_input_class< input_class_climate_t >())
119  {
120  return m_iokcomm->get_input_class< climate::input_class_climate_t >()->station_info()->rainfall_intensity * cbm::M_IN_MM;
121  }
122  else
123  {
124  return ldndc::climate::climate_info_defaults.rainfall_intensity * cbm::M_IN_MM;
125  }
126 }
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
2007 {
2008  //percolation and direct evaporation from the (upper) soil
2009  //if there is snow cover, evaporation from soil will be set to 0, and only sublimation is allowed.
2010  double const wc_min_evaporation(m_sipar->WCDNDC_EVALIM_FRAC_WCMIN() * wc_.wc_wp_sl[_sl]);
2011  if ( cbm::flt_less( sc_.depth_sl[_sl] - sc_.h_sl[_sl], m_sipar->EVALIM()) &&
2012  cbm::flt_greater( wc_.wc_sl[_sl], wc_min_evaporation) &&
2013  cbm::flt_greater_zero( hr_pot_evapotrans))
2014  {
2015  //limiting factor for soil water infiltration
2016  double const f_clay( 1.0 + m_sipar->SLOPE_CLAYF() * std::min(1.0 - sc_.fcorg_sl[_sl] / cbm::CCORG, sc_.clay_sl[_sl]));
2017  double const f_water( cbm::bound( 0.0,
2018  pow( (wc_.wc_sl[_sl] - wc_min_evaporation) / wc_.wc_fc_sl[_sl], f_clay),
2019  1.0));
2020  double const h_rest( std::min( sc_.h_sl[_sl], m_sipar->EVALIM() - ( sc_.depth_sl[_sl] - sc_.h_sl[_sl])));
2021  double const f_depth( h_rest / m_sipar->EVALIM() * std::max(0.0,
2022  1.0 - std::min((double)m_sipar->EVALIM(),
2023  sc_.depth_sl[_sl] - ( 0.5 * sc_.h_sl[_sl]))
2024  / m_sipar->EVALIM()));
2025  double const soilevaporation_max( cbm::bound_max( f_water * wc_.wc_sl[_sl] * sc_.h_sl[_sl], hr_pot_evapotrans));
2026 
2027  double const soilevaporation( cbm::bound( 0.0,
2028  hr_pot_evapotrans * f_depth * f_water,
2029  soilevaporation_max));
2030 
2031  //update of the current soil layer water content
2032  wc_.wc_sl[_sl] -= soilevaporation / sc_.h_sl[_sl];
2033  wc_.accumulated_soilevaporation += soilevaporation;
2034  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - soilevaporation);
2035  }
2036 }

◆ WatercycleDNDC_preferential_flow()

void ldndc::WatercycleDNDC::WatercycleDNDC_preferential_flow ( )
private
Parameters
[in]None
[out]None
Returns
None
2379 {
2380  //BY_PASSF is defined as the fraction of surface water that passes per hour (0.0 by default)
2381  if ( cbm::flt_greater_zero(m_sipar->BY_PASSF()) &&
2382  cbm::flt_greater_zero( wc_.surface_water))
2383  {
2384 
2385  double water_def_total( 0.0);
2386  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2387  {
2388  //water flux that flow through macro pores can maximally store
2389  water_def_total += cbm::bound_min( 0.0, (wc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]);
2390  }
2391 
2392  //duration of timestep (hr-1)
2393  double const fhr( get_time_step_duration());
2394 
2395 
2396  //fraction of surface water that is put into bypass flow (m timestep-1)
2397  double const bypass_fraction( cbm::bound_max(m_sipar->BY_PASSF() * fhr, 0.99));
2398  double const bypass_water( wc_.surface_water * bypass_fraction);
2399  wc_.surface_water -= bypass_water;
2400  wc_.accumulated_infiltration += bypass_water;
2401 
2402  if ( cbm::flt_greater_equal( water_def_total, bypass_water))
2403  {
2404 
2405  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2406  {
2407  double const water_def( cbm::bound_min( 0.0, (wc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]));
2408  double const water_add( water_def / water_def_total * bypass_water);
2409 
2410  wc_.wc_sl[sl] += water_add / sc_.h_sl[sl];
2411  }
2412  }
2413  else
2414  {
2415  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2416  {
2417  double const water_def( cbm::bound_min( 0.0, (wc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]));
2418  wc_.wc_sl[sl] += water_def / sc_.h_sl[sl];
2419  }
2420  /* add surplus to percolation out of last soil layer */
2421  wc_.accumulated_percolation += bypass_water - water_def_total;
2422  wc_.accumulated_waterflux_sl[m_soillayers->soil_layer_cnt()-1] += bypass_water - water_def_total;
2423  }
2424  }
2425 }
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2565

◆ 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

2288 {
2289  double const sks( WaterFlowSaturated( _sl, _fact_impedance));
2290  double const fsl( get_fsl( _sl));
2291  double const fhr( get_time_step_duration());
2292 
2294  double slope(( _sl < m_soillayers->soil_layers_in_litter_cnt()) ? m_sipar->SLOPE_FF() : m_sipar->SLOPE_MS());
2295 
2296  double const travel_time( 1.0 - log10( wc_.sks_sl[_sl] * cbm::M_IN_CM * cbm::MIN_IN_DAY * fhr));
2297  double const travel_time_factor( cbm::flt_greater_zero( travel_time) ? (1.0 - exp(-1.0 / travel_time)) : 1.0);
2298  return cbm::bound_max( cbm::sqr(1.0 - (wc_.wc_fc_sl[_sl] / ((wc_.wc_sl[_sl] + wc_.ice_sl[_sl]) * slope)))
2299  * wc_.wc_sl[_sl] * sc_.h_sl[_sl] * fsl * travel_time_factor * _fact_impedance,
2300  sks);
2301 }
double get_time_step_duration()
Time fraction with regard to daily rates.
Definition: watercycle-dndc.cpp:2565
double get_fsl(size_t)
Returns factor accounting for soil layer depth. Originally, a constant soil layer depth of 0...
Definition: watercycle-dndc.cpp:2461

Member Data Documentation

◆ IMAX_W

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

number of iteration steps within a day