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.

2585 {
2586  double balance = get_leaf_water() + wc_.surface_ice + wc_.surface_water;
2587  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2588  {
2589  balance += (wc_.wc_sl[sl] + wc_.ice_sl[sl]) * sc_.h_sl[sl];
2590  }
2591  balance += ( - m_wc.accumulated_irrigation_event_executed
2592  - wc_.accumulated_irrigation_automatic
2593  - wc_.accumulated_irrigation_ggcmi
2594  - wc_.accumulated_irrigation_reservoir_withdrawal
2595  - wc_.accumulated_precipitation
2596  - wc_.accumulated_groundwater_access
2597  + wc_.accumulated_percolation
2598  + wc_.accumulated_runoff
2599  + wc_.accumulated_wateruptake_sl.sum()
2600  + wc_.accumulated_interceptionevaporation
2601  + wc_.accumulated_soilevaporation
2602  + wc_.accumulated_surfacewaterevaporation);
2603 
2604  if ( _balance)
2605  {
2606  double const balance_delta( std::abs( *_balance - balance));
2607  if ( balance_delta > cbm::DIFFMAX)
2608  {
2609  KLOGWARN( "Water leakage: difference is ", balance - *_balance);
2610  }
2611  }
2612 
2613  for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
2614  {
2615  if ( cbm::flt_less( wc_.wc_sl[sl], 0.0))
2616  {
2617  KLOGWARN( "Negative water content of: ", wc_.wc_sl[sl], " in layer: ", sl, ". Water content set to zero.");
2618  wc_.wc_sl[sl] = 0.0;
2619  }
2620  if ( cbm::flt_less( wc_.ice_sl[sl], 0.0))
2621  {
2622  KLOGWARN( "Negative ice content of: ", wc_.ice_sl[sl], " in layer: ", sl, ". Ice content set to zero.");
2623  wc_.ice_sl[sl] = 0.0;
2624  }
2625  }
2626 
2627  return balance;
2628 }
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
1948 {
1949  //Average limit of snowEvaporation is 0.00013m day-1 (Granger and Male 1978)
1950  double hr_potESnow( cbm::bound(0.0, hr_pot_evapotrans, 0.00013 * nhr/24.0));
1951 
1952  if ( cbm::flt_greater_zero( wc_.surface_ice) &&
1953  cbm::flt_greater_zero( hr_potESnow))
1954  {
1955  //If there is a snowpack and there is a potential to evaporate
1956  if ( wc_.surface_ice > hr_potESnow)
1957  {
1958  //Partial evaporation from the snow
1959  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - hr_potESnow);
1960  wc_.accumulated_surfacewaterevaporation += hr_potESnow;
1961  wc_.surface_ice -= hr_potESnow;
1962  }
1963  else
1964  {
1965  //Total evaporation of the snowpack
1966  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - wc_.surface_ice);;
1967  wc_.accumulated_surfacewaterevaporation += wc_.surface_ice;
1968  wc_.surface_ice = 0.0;
1969  }
1970  }
1971 }

◆ CalcSoilEvapoPercolation()

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

◆ CalcSurfaceFlux()

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

◆ CalcSurfaceWaterEvaporation()

void ldndc::WatercycleDNDC::CalcSurfaceWaterEvaporation ( )
private
Parameters
[in]None
[out]None
Returns
None
1917 {
1918  if ( cbm::flt_greater_zero( wc_.surface_water) &&
1919  cbm::flt_greater_zero( hr_pot_evapotrans))
1920  {
1921  if ( hr_pot_evapotrans > wc_.surface_water)
1922  {
1923  wc_.accumulated_surfacewaterevaporation += wc_.surface_water;
1924  hr_pot_evapotrans -= wc_.surface_water;
1925  wc_.surface_water = 0.0;
1926  }
1927  else
1928  {
1929  wc_.surface_water -= hr_pot_evapotrans;
1930  wc_.accumulated_surfacewaterevaporation += hr_pot_evapotrans;
1931  hr_pot_evapotrans = 0.0;
1932  }
1933  }
1934 }

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

◆ 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:2573

◆ 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
2015 {
2016  //percolation and direct evaporation from the (upper) soil
2017  //if there is snow cover, evaporation from soil will be set to 0, and only sublimation is allowed.
2018  double const wc_min_evaporation(m_sipar->WCDNDC_EVALIM_FRAC_WCMIN() * wc_.wc_wp_sl[_sl]);
2019  if ( cbm::flt_less( sc_.depth_sl[_sl] - sc_.h_sl[_sl], m_sipar->EVALIM()) &&
2020  cbm::flt_greater( wc_.wc_sl[_sl], wc_min_evaporation) &&
2021  cbm::flt_greater_zero( hr_pot_evapotrans))
2022  {
2023  //limiting factor for soil water infiltration
2024  double const f_clay( 1.0 + m_sipar->SLOPE_CLAYF() * std::min(1.0 - sc_.fcorg_sl[_sl] / cbm::CCORG, sc_.clay_sl[_sl]));
2025  double const f_water( cbm::bound( 0.0,
2026  pow( (wc_.wc_sl[_sl] - wc_min_evaporation) / wc_.wc_fc_sl[_sl], f_clay),
2027  1.0));
2028  double const h_rest( std::min( sc_.h_sl[_sl], m_sipar->EVALIM() - ( sc_.depth_sl[_sl] - sc_.h_sl[_sl])));
2029  double const f_depth( h_rest / m_sipar->EVALIM() * std::max(0.0,
2030  1.0 - std::min((double)m_sipar->EVALIM(),
2031  sc_.depth_sl[_sl] - ( 0.5 * sc_.h_sl[_sl]))
2032  / m_sipar->EVALIM()));
2033  double const soilevaporation_max( cbm::bound_max( f_water * wc_.wc_sl[_sl] * sc_.h_sl[_sl], hr_pot_evapotrans));
2034 
2035  double const soilevaporation( cbm::bound( 0.0,
2036  hr_pot_evapotrans * f_depth * f_water,
2037  soilevaporation_max));
2038 
2039  //update of the current soil layer water content
2040  wc_.wc_sl[_sl] -= soilevaporation / sc_.h_sl[_sl];
2041  wc_.accumulated_soilevaporation += soilevaporation;
2042  hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - soilevaporation);
2043  }
2044 }

◆ WatercycleDNDC_preferential_flow()

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

◆ 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

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

Member Data Documentation

◆ IMAX_W

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

number of iteration steps within a day