LandscapeDNDC  1.36.0
ldndc::EcHy Class Reference

Watercycle model EcosystemHydrology - EcHy. More...

Inherits ldndc::MBE_LegacyModel.

Public Member Functions

lerr_t solve ()
 

Private Member Functions

lerr_t EcHyIrrigation ()
 Irrigation.
 
lerr_t EcHyFlood ()
 sets hydrologic conditions during flooding events, e.g., More...
 
lerr_t EcHyPercolation (size_t, size_t)
 Calculates water percolation within the soil profile.
 
lerr_t EcHyBypassFlow (double, double &)
 Calculates water percolation within the soil profile.
 
lerr_t EcHyEvapotranspiration ()
 Calculates evapotranspiration within the soil profile.
 
lerr_t EcHySubsl2 (double const &, double const &, int const &, double const &, double &)
 ... More...
 
lerr_t EcHyBalanceCheck (double &)
 ... More...
 
void EcHyreset ()
 ...
 
lerr_t EcHyGroundwater ()
 ... More...
 
lerr_t EcHyIntegration ()
 ... More...
 
lerr_t EcHyStepInit ()
 ...
 
lerr_t EcHyStepExit ()
 ...
 
double EcHyGetInterceptionCapacity ()
 ...
 
lerr_t EcHyCalculateLeafWaterDistribution (lvector_t< double > &)
 ...
 
lerr_t EcHyPotentialEvapotranspiration ()
 Calculates potential evapotranspiration. Specific concept can be given as model option.
 
lerr_t EcHySnowIce ()
 Calls SnowDNDC for the calculation of snowpack and soil ice formation.
 
double EcHyGetMinimumWater (size_t)
 ...
 
double EcHyGetWiltingPoint (size_t)
 ...
 
double EcHyGetAvailableWaterTranspiration (size_t)
 ...
 
double EcHyGetAvailableWaterEvaporation (size_t)
 ...
 
double EcHyGetRootLimitation (double)
 ...
 
double EcHyGetWaterLimitationTranspiration (double, size_t)
 ...
 
double EcHySoilWaterChange (size_t)
 

Private Attributes

double kst_bottom
 
double gw_depth_static
 
double accumulated_potentialtranspiration_old
 
double accumulated_irrigation_old
 
lvector_t< double > trwl_sl
 
lvector_t< double > wlfl_sl
 
lvector_t< double > cr_fill_sl
 
lvector_t< double > gw_fill_sl
 
lvector_t< double > bypass_fill_sl
 
lvector_t< double > wc_sat_sl
 
lvector_t< double > wc_res_sl
 
lvector_t< double > evsws_sl
 
lvector_t< double > kst_sl
 
lvector_t< double > wl_sl
 
lvector_t< double > wlfc_sl
 
lvector_t< double > wlwp_sl
 
lvector_t< double > wlst_sl
 
double ev_leaf
 
double gw_fill_surface
 
double cr_fill_groundwater
 
double runoff
 
double snowfall
 
double throughfall
 
double interception_water
 
double thornthwaite_heat_index
 
double daily_potential_evapotranspiration
 
double daily_potential_leaf_evaporation
 
double daily_potential_soil_evaporation
 
double daily_potential_transpiration
 
WaterCycleSnowDNDC::IceContentStateIn m_icecontent_in
 
cbm::string_t evapotranspiration_method
 
cbm::string_t runoff_method
 

Detailed Description

Watercycle model EcosystemHydrology - EcHy.

Member Function Documentation

◆ EcHyBalanceCheck()

lerr_t ldndc::EcHy::EcHyBalanceCheck ( double &  _balance)
private

...

Checks balance between all incoming and outgoing water fluxes.

References cr_fill_groundwater, ev_leaf, evsws_sl, gw_fill_sl, gw_fill_surface, interception_water, runoff, snowfall, throughfall, trwl_sl, wl_sl, and wlfl_sl.

Referenced by solve().

455 {
456  double balance( interception_water + surface_water + surface_ice);
457  for ( size_t sl = 0; sl < soillayers_in->soil_layer_cnt(); ++sl)
458  {
459  balance += wl_sl[sl] + wc_.ice_sl[sl] * sc_.h_sl[sl];
460  }
461 
462  if ( _balance > 0.0)
463  {
464  balance += wlfl_sl[soillayers_in->soil_layer_cnt()]
465  + ev_leaf + evsws_sl.sum() + trwl_sl.sum() + runoff
467  - irrigation - throughfall - snowfall;
468 
469  double const balance_delta( std::abs( _balance - balance));
470  if ( cbm::flt_greater( balance_delta, 1.0e-4))
471  {
472  KLOGWARN( "Water leakage in: ", name(),
473  " Difference: ", balance - _balance);
474  return LDNDC_ERR_FAIL;
475  }
476  }
477  else
478  {
479  _balance = balance;
480  }
481 
482  return LDNDC_ERR_OK;
483 }
double gw_fill_surface
Definition: echy.h:128
double ev_leaf
Definition: echy.h:125
lvector_t< double > evsws_sl
Definition: echy.h:107
double cr_fill_groundwater
Definition: echy.h:131
double interception_water
Definition: echy.h:143
lvector_t< double > gw_fill_sl
Definition: echy.h:95
lvector_t< double > wl_sl
Definition: echy.h:113
double snowfall
Definition: echy.h:137
double runoff
Definition: echy.h:134
lvector_t< double > wlfl_sl
Definition: echy.h:89
double throughfall
Definition: echy.h:140
lvector_t< double > trwl_sl
Definition: echy.h:81
Here is the caller graph for this function:

◆ EcHyFlood()

lerr_t ldndc::EcHy::EcHyFlood ( )
private

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

Flooding.

  • surface water table
  • bund height

References kst_bottom, and kst_sl.

Referenced by solve().

18 {
19  lerr_t rc = m_eventflood.solve();
20  if ( rc != LDNDC_ERR_OK){ return rc; }
21 
22  bund_height = m_eventflood.get_bund_height();
23 
24  /* todo */
25  //irrigation_amount = m_eventflood.get_irrigation_amount();
26 
27  /* todo */
28  //double const irrigation_height = m_eventflood.get_irrigation_height();
29  double const max_percolation = m_eventflood.get_maximum_percolation();
30  double const water_table_flooding = m_eventflood.get_water_table();
31 
32 
33  irrigation_switch = NONE;
34  minimum_watertable_height = 0.0;
35  if ( cbm::is_valid( water_table_flooding))
36  {
37  if ( cbm::flt_greater_zero( water_table_flooding))
38  {
39  irrigation_switch = CONSTANT_POSITIVE_WATER_TABLE;
40  minimum_watertable_height = water_table_flooding;
41  bund_height = water_table_flooding;
42  }
43  else
44  {
45  irrigation_switch = CONSTANT_NEGATIVE_WATER_TABLE;
46  minimum_watertable_height = water_table_flooding;
47  bund_height = 0.0;
48  }
49  }
50 
51 
52  size_t sl_max( soillayers_in->soil_layer_cnt()-1);
53  if ( cbm::is_valid( max_percolation))
54  {
55  kst_sl[sl_max] = max_percolation / lclock()->time_resolution();
56  }
57  else
58  {
59  double const time_rate( 0.1 / lclock()->time_resolution());
60  kst_sl[sl_max] = cbm::bound_min( 0.0,
61  kst_sl[sl_max] - (kst_sl[sl_max] - kst_bottom) * time_rate);
62  }
63 
64  return LDNDC_ERR_OK;
65 }
lvector_t< double > kst_sl
Definition: echy.h:110
double kst_bottom
Definition: echy.h:67
Here is the caller graph for this function:

◆ EcHyGroundwater()

lerr_t ldndc::EcHy::EcHyGroundwater ( )
private

...

Negative groundwater table represents water on the soil surface

Negative groundwater table represents water on the soil surface

References cr_fill_groundwater, cr_fill_sl, EcHySoilWaterChange(), gw_depth_static, gw_fill_sl, gw_fill_surface, kst_sl, wc_res_sl, wc_sat_sl, wl_sl, and wlst_sl.

Referenced by solve().

21 {
22  cbm::invalidate( kst_bottom_gw_vertical);
23  if ( cbm::flt_greater( sc_.depth_sl[soillayers_in->soil_layer_cnt()-1], gw_depth_static))
24  {
29  if ( cbm::flt_greater( -gw_depth_static, surface_water))
30  {
31  gw_fill_surface = -gw_depth_static - surface_water;
32  }
33 
34  //groundwater access
35  {
36  size_t sl( 0);
37  double const layer_midpoint_depth( 0.5 * sc_.h_sl[sl]);
38  if ( cbm::flt_greater_equal( layer_midpoint_depth, gw_depth_static))
39  {
40  gw_fill_sl[sl] = cbm::bound_min( 0.0,
41  wlst_sl[sl] - (wl_sl[sl] + EcHySoilWaterChange( sl)));
42  wc_.accumulated_groundwater_access += gw_fill_sl[sl];
43  }
44  else
45  {
46  double const wl_new( wl_sl[sl] + EcHySoilWaterChange( sl));
47  double const wc_new( cbm::bound( wc_res_sl[sl], wl_new / sc_.h_sl[sl], wc_sat_sl[sl]));
48  kst_bottom_gw_vertical = cbm::bound_min( 0.0,
49  ldndc::hydraulic_conductivity(
50  wc_new,
51  wc_.vgm_sl[sl],
52  wc_.mvg_tau_sl[sl],
53  wc_sat_sl[sl],
54  wc_res_sl[sl],
55  kst_sl[sl]));
56  gw_fill_sl[sl] = 0.0;
57  }
58  }
59 
60  for (size_t sl = 1; sl < soillayers_in->soil_layer_cnt()-1; sl++)
61  {
62  double const layer_midpoint_depth( sc_.depth_sl[sl] - 0.5 * sc_.h_sl[sl]);
63  if ( cbm::flt_greater_equal( layer_midpoint_depth, gw_depth_static))
64  {
65  gw_fill_sl[sl] = cbm::bound_min( 0.0,
66  wlst_sl[sl] - (wl_sl[sl] + EcHySoilWaterChange( sl)));
67  wc_.accumulated_groundwater_access += gw_fill_sl[sl];
68  }
69  else
70  {
71  double const wl_new( wl_sl[sl] + EcHySoilWaterChange( sl));
72  double const wc_new( cbm::bound( wc_res_sl[sl], wl_new / sc_.h_sl[sl], wc_sat_sl[sl]));
73  kst_bottom_gw_vertical = cbm::bound_min( 0.0,
74  ldndc::hydraulic_conductivity(
75  wc_new,
76  wc_.vgm_sl[sl],
77  wc_.mvg_tau_sl[sl],
78  wc_sat_sl[sl],
79  wc_res_sl[sl],
80  kst_sl[sl]));
81  gw_fill_sl[sl] = 0.0;
82  }
83  }
84 
85  {
86  size_t sl( soillayers_in->soil_layer_cnt()-1);
87  double const layer_midpoint_depth( sc_.depth_sl[sl] - 0.5 * sc_.h_sl[sl]);
88  if ( cbm::flt_greater_equal( layer_midpoint_depth, gw_depth_static))
89  {
90  gw_fill_sl[sl] = cbm::bound_min( 0.0,
91  wlst_sl[sl] - (wl_sl[sl] + EcHySoilWaterChange( sl)));
92  wc_.accumulated_groundwater_access += gw_fill_sl[sl];
93  }
94  else
95  {
96  double const wl_new( wl_sl[sl] + EcHySoilWaterChange( sl));
97  double const wc_new( cbm::bound( wc_res_sl[sl], wl_new / sc_.h_sl[sl], wc_sat_sl[sl]));
98  kst_bottom_gw_vertical = cbm::bound_min( 0.0,
99  ldndc::hydraulic_conductivity(
100  wc_new,
101  wc_.vgm_sl[sl],
102  wc_.mvg_tau_sl[sl],
103  wc_sat_sl[sl],
104  wc_res_sl[sl],
105  kst_sl[sl]));
106  gw_fill_sl[sl] = 0.0;
107  }
108  }
109  }
110 
111  //capillary rise
112  if ( m_param->WCDNDC_HAVE_CAPILLARY_ACTION())
113  {
114  size_t const STEPS( 10);
115  for (size_t steps = 0; steps < STEPS; steps++)
116  {
117  for (size_t sl = 0; sl < soillayers_in->soil_layer_cnt()-1; sl++)
118  {
119  double const layer_midpoint_depth( sc_.depth_sl[sl] - 0.5 * sc_.h_sl[sl]);
120  if ( cbm::flt_greater_equal( layer_midpoint_depth, gw_depth_static))
121  {
122  break;
123  }
124 
125  /* Water contents with until now calculated water balance */
126  double const wl_new_0( wl_sl[sl] + EcHySoilWaterChange( sl));
127  double const wl_new_1( wl_sl[sl+1] + EcHySoilWaterChange( sl+1));
128 
129  double const wc_vg_0( cbm::bound( wc_res_sl[sl], wl_new_0 / sc_.h_sl[sl], wc_sat_sl[sl]));
130  double const wc_vg_1( cbm::bound( wc_res_sl[sl+1], wl_new_1 / sc_.h_sl[sl+1], wc_sat_sl[sl+1]));
131 
132  /* Unsaturated hydraulic conductivity with until now calculated water balance */
133  double const kust_0( cbm::bound_min( 0.0,
134  ldndc::hydraulic_conductivity(
135  wc_vg_0, wc_.vgm_sl[sl], wc_.mvg_tau_sl[sl],
136  wc_sat_sl[sl], wc_res_sl[sl],
137  1.0/double(STEPS)*kst_sl[sl])));
138  double const kust_1( cbm::bound_min( 0.0,
139  ldndc::hydraulic_conductivity(
140  wc_vg_1, wc_.vgm_sl[sl+1], wc_.mvg_tau_sl[sl],
141  wc_sat_sl[sl+1], wc_res_sl[sl+1],
142  1.0/double(STEPS)*kst_sl[sl+1])));
143 
144  double const capillary_pressure_0( ldndc::capillary_pressure(
145  wl_new_0/sc_.h_sl[sl], wc_.vga_sl[sl],
146  wc_.vgn_sl[sl], wc_.vgm_sl[sl],
147  wc_sat_sl[sl], wc_res_sl[sl]));
148  double const capillary_pressure_1( ldndc::capillary_pressure(
149  wl_new_1/sc_.h_sl[sl+1], wc_.vga_sl[sl+1],
150  wc_.vgn_sl[sl+1], wc_.vgm_sl[sl+1],
151  wc_sat_sl[sl], wc_res_sl[sl+1]));
152 
153  double const kust( cbm::flt_greater( capillary_pressure_0, capillary_pressure_1) ? kust_1 : kust_0 );
154  double const flow( (capillary_pressure_1 - capillary_pressure_0) / (0.5 * (sc_.h_sl[sl] + sc_.h_sl[sl+1])) * kust);
155 
156  /* downward flow*/
157  if ( cbm::flt_greater_zero( flow))
158  {
159  double const left( wl_sl[sl] + EcHySoilWaterChange( sl) - (wc_res_sl[sl] * sc_.h_sl[sl]));
160  double const space( wlst_sl[sl+1] - (wl_sl[sl+1] + EcHySoilWaterChange( sl+1)));
161 
162  if ( cbm::flt_greater_zero( left) &&
163  cbm::flt_greater_zero( space))
164  {
165  double const bound_flow( 1.0 / static_cast<double>(STEPS) * cbm::bound_max( flow, std::min( left, space)));
166  cr_fill_sl[sl] -= bound_flow;
167  cr_fill_sl[sl+1] += bound_flow;
168  }
169  }
170  /* upward flow */
171  else if ( cbm::flt_less( flow, 0.0))
172  {
173  double const left( wl_sl[sl+1] + EcHySoilWaterChange( sl+1) - (wc_res_sl[sl+1] * sc_.h_sl[sl+1]));
174  double const space( wlst_sl[sl] - (wl_sl[sl] + EcHySoilWaterChange( sl)));
175 
176  if ( cbm::flt_greater_zero( left) &&
177  cbm::flt_greater_zero( space))
178  {
179  double const bound_flow( cbm::bound_max( -flow, 1.0 / static_cast<double>(STEPS) * std::min( left, space)));
180  cr_fill_sl[sl] += bound_flow;
181  cr_fill_sl[sl+1] -= bound_flow;
182  }
183  }
184  }
185 
186  size_t sl( soillayers_in->soil_layer_cnt()-1);
187  double const layer_midpoint_depth( sc_.depth_sl[sl] - 0.5 * sc_.h_sl[sl]);
188  if ( cbm::flt_greater( gw_depth_static, layer_midpoint_depth))
189  {
190  /* Water contents with until now calculated water balance */
191  double const wl_new( wl_sl[sl] + EcHySoilWaterChange( sl));
192 
193  double const wc_vg( cbm::bound( wc_res_sl[sl], wl_new / sc_.h_sl[sl], wc_sat_sl[sl]));
194 
195  /* Unsaturated hydraulic conductivity with until now calculated water balance */
196  double const kust( cbm::bound_min( 0.0,
197  ldndc::hydraulic_conductivity(
198  wc_vg, wc_.vgm_sl[sl], wc_.mvg_tau_sl[sl],
199  wc_sat_sl[sl], wc_res_sl[sl],
200  1.0/double(STEPS)*kst_sl[sl])));
201 
202  double const capillary_pressure( ldndc::capillary_pressure(
203  wl_new/sc_.h_sl[sl], wc_.vga_sl[sl],
204  wc_.vgn_sl[sl], wc_.vgm_sl[sl],
205  wc_sat_sl[sl], wc_res_sl[sl]));
206 
207  /* upward flow */
208  double const flow( -capillary_pressure / (gw_depth_static - layer_midpoint_depth) * kust);
209  if ( cbm::flt_less( flow, 0.0))
210  {
211  double const space( wlst_sl[sl] - (wl_sl[sl] + EcHySoilWaterChange( sl)));
212  if ( cbm::flt_greater_zero( space))
213  {
214  double const bound_flow( cbm::bound_max( -flow, 1.0 / static_cast<double>(STEPS) * space));
215  cr_fill_sl[sl] += bound_flow;
216  cr_fill_groundwater += bound_flow;
217  }
218  }
219  }
220  }
221  }
222 
223  return LDNDC_ERR_OK;
224 }
double gw_fill_surface
Definition: echy.h:128
double gw_depth_static
Definition: echy.h:72
lvector_t< double > kst_sl
Definition: echy.h:110
double cr_fill_groundwater
Definition: echy.h:131
lvector_t< double > gw_fill_sl
Definition: echy.h:95
lvector_t< double > wl_sl
Definition: echy.h:113
lvector_t< double > wc_res_sl
Definition: echy.h:104
double EcHySoilWaterChange(size_t)
Definition: echy.cpp:417
lvector_t< double > wc_sat_sl
Definition: echy.h:101
lvector_t< double > wlst_sl
Definition: echy.h:122
lvector_t< double > cr_fill_sl
Definition: echy.h:92
Here is the call graph for this function:
Here is the caller graph for this function:

◆ EcHyIntegration()

lerr_t ldndc::EcHy::EcHyIntegration ( )
private

...

Integrates all water state variables.

References bypass_fill_sl, EcHySoilWaterChange(), gw_fill_surface, runoff, throughfall, wc_res_sl, wc_sat_sl, wl_sl, and wlfl_sl.

Referenced by solve().

431 {
432  surface_water = cbm::bound_min( 0.0,
433  surface_water +
434  (throughfall + gw_fill_surface + irrigation - wlfl_sl[0] - runoff - bypass_fill_sl.sum()));
435 
436  for (size_t sl = 0; sl < soillayers_in->soil_layer_cnt(); sl++)
437  {
438  wl_sl[sl] = cbm::bound_min( 0.0, wl_sl[sl] + EcHySoilWaterChange( sl));
439  wc_.mskpa_sl[sl] = ldndc::hydrology::capillary_pressure( wl_sl[sl]/sc_.h_sl[sl], wc_.vga_sl[sl],
440  wc_.vgn_sl[sl], wc_.vgm_sl[sl], wc_sat_sl[sl], wc_res_sl[sl]);
441  }
442 
443  return LDNDC_ERR_OK;
444 }
double gw_fill_surface
Definition: echy.h:128
lvector_t< double > wl_sl
Definition: echy.h:113
double runoff
Definition: echy.h:134
lvector_t< double > wlfl_sl
Definition: echy.h:89
lvector_t< double > wc_res_sl
Definition: echy.h:104
double throughfall
Definition: echy.h:140
double EcHySoilWaterChange(size_t)
Definition: echy.cpp:417
lvector_t< double > wc_sat_sl
Definition: echy.h:101
lvector_t< double > bypass_fill_sl
Definition: echy.h:98
Here is the call graph for this function:
Here is the caller graph for this function:

◆ EcHySoilWaterChange()

double ldndc::EcHy::EcHySoilWaterChange ( size_t  _sl)
private

Integrates all water state variables.

References cr_fill_sl, evsws_sl, gw_fill_sl, trwl_sl, and wlfl_sl.

Referenced by EcHyGetAvailableWaterEvaporation(), EcHyGetAvailableWaterTranspiration(), EcHyGroundwater(), EcHyIntegration(), and EcHyIrrigation().

419 {
420  return wlfl_sl[_sl] - wlfl_sl[_sl+1] + cr_fill_sl[_sl] - trwl_sl[_sl] - evsws_sl[_sl] + gw_fill_sl[_sl];
421 }
lvector_t< double > evsws_sl
Definition: echy.h:107
lvector_t< double > gw_fill_sl
Definition: echy.h:95
lvector_t< double > wlfl_sl
Definition: echy.h:89
lvector_t< double > trwl_sl
Definition: echy.h:81
lvector_t< double > cr_fill_sl
Definition: echy.h:92
Here is the caller graph for this function:

◆ EcHySubsl2()

lerr_t ldndc::EcHy::EcHySubsl2 ( double const &  _cp,
double const &  _delta_z,
int const &  _sl,
double const &  _saturation,
double &  _flow 
)
private

...

This routine calculates the rate of capillary flow or percolation between groundwater table and root zone. The stationary flow is found by integration of dZL = K.d(MH)/(K + FLW), where Z= height above groundwater, MH= matric head, K= conductivity and FLW= chosen flow. In an iteration loop the correct flow is found. The integration goes at most over four intervals: [0,45], [45,170], [170,330] and [330,MH-rootzone] (last one on logarithmic scale).

Author
C. Rappoldt M. Wopereis
Date
January 1986, revised June 1990

Chapter 15 in documentation WOFOST Version 4.1 (1988)

References kst_sl.

255 {
256  //calculation of pf from matrix head
257  double pf( std::log10( _cp));
258 
259  //in case of small matric head (high water contents)
260  if ( cbm::flt_greater( 1.0, pf))
261  {
262  _flow = 0.0;
263  return LDNDC_ERR_OK;
264  }
265 
266  double elog10( 2.302585);
267  double logst4( 2.518514);
268 
269  double start[4] = {0.0, 45.0, 170.0, 330.0};
270 
271  //number and width of integration intervals
272  int iint( 0);
273  double del[4] = { 0.0, 0.0, 0.0, 0.0};
274  for (int i1 = 0; i1 < 4; i1++)
275  {
276  if (i1 < 3)
277  {
278  del[i1] = std::min( start[i1+1], _cp) - start[i1];
279  }
280  else
281  {
282  del[i1] = pf - logst4;
283  }
284  if (del[i1] <= 0.0)
285  {
286  break;
287  }
288  iint += 1;
289  }
290 
291  //preparation of three-point gaussian integration
292  double hulp[12] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
293  double conduc[12] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
294  for (int i1 = 0; i1 < iint; i1++)
295  {
296  for (int i2 = 0; i2 < 3; i2++)
297  {
298  double pgau[3] = {0.1127016654, 0.5, 0.8872983346};
299  double wgau[3] = {0.2777778, 0.4444444, 0.2777778};
300 
301  int i3( 3 * i1 + i2);
302 
303  //the three points in the full-width intervals are standard
304  double pfstan[9] = {0.705143, 1.352183, 1.601282, 1.771497, 2.031409, 2.192880, 2.274233, 2.397940, 2.494110};
305  double pfgau( pfstan[i3]);
306 
307  //the three points in the last interval are calculated
308  if ( i1 == (iint-1))
309  {
310  if ( iint < 4)
311  {
312  pfgau = std::log10( start[iint-1] + pgau[i2] * del[iint-1]);
313  }
314  else
315  {
316  pfgau = logst4 + pgau[i2] * del[iint-1];
317  }
318  }
319 
320  double const wcl_loc( ldndc::hydrology::water_content( std::exp(elog10 * pfgau), wc_.vga_sl[_sl], wc_.vgn_sl[_sl], wc_.vgm_sl[_sl], _saturation, 0.1));
321  conduc[i3] = ldndc::hydraulic_conductivity( wcl_loc, wc_.vgm_sl[_sl], wc_.mvg_tau_sl[_sl],
322  _saturation, sc_.wc_wp_sl[_sl], kst_sl[_sl] * cbm::CM_IN_M);
323  hulp[i3] = del[i1] * wgau[i2] * conduc[i3];
324  if (i3 > 8)
325  {
326  hulp[i3] = hulp[i3] * elog10 * std::exp( elog10 * pfgau);
327  }
328  }
329  }
330 
331 
332  //setting upper and lower limit
333  double const wcl_loc( ldndc::hydrology::water_content( _cp, wc_.vga_sl[_sl], wc_.vgn_sl[_sl], wc_.vgm_sl[_sl], _saturation, 0.1));
334  double const kms( ldndc::hydraulic_conductivity( wcl_loc, wc_.vgm_sl[_sl], wc_.mvg_tau_sl[_sl],
335  _saturation, sc_.wc_wp_sl[_sl], kst_sl[_sl] * cbm::CM_IN_M));
336  double fu( (_cp <= _delta_z) ? 0.0 : 1.27);
337  double fl( (_cp >= _delta_z) ? 0.0 : -kms);
338 
339 
340  if ( !cbm::flt_equal( _cp, _delta_z))
341  {
342  //iteration loop
343  int imax( 3 * iint);
344  for ( int i1 = 0; i1 < 14; i1++)
345  {
346  double flw( (fu + fl) / 2.0);
347  double df( (fu - fl) / 2.0);
348  if ( (df < 0.01) &&
349  ((df / std::abs(flw)) < 0.1))
350  {
351  break;
352  }
353 
354  double z( 0.0);
355  for ( int i2 = 0; i2 < imax; i2++)
356  {
357  z = z + hulp[i2] / (conduc[i2] + flw);
358  }
359  if (z >= _delta_z)
360  {
361  fl = flw;
362  }
363  if (z <= _delta_z)
364  {
365  fu = flw;
366  }
367  }
368  }
369 
370  //flow output in mm
371  _flow = 10.0 * (fu + fl) / 2.0;
372 
373  return LDNDC_ERR_OK;
374 }
lvector_t< double > kst_sl
Definition: echy.h:110

◆ solve()

lerr_t ldndc::EcHy::solve ( )

Kicks off computation for one time step.

References EcHyBalanceCheck(), EcHyBypassFlow(), EcHyEvapotranspiration(), EcHyFlood(), EcHyGroundwater(), EcHyIntegration(), EcHyIrrigation(), EcHyPercolation(), EcHyreset(), EcHySnowIce(), EcHyStepExit(), EcHyStepInit(), kst_bottom, kst_sl, runoff, runoff_method, throughfall, and wlfl_sl.

46 {
47  /* reset of all fluxes */
48  EcHyreset();
49 
50  /* flooding events */
51  lerr_t rc = EcHyFlood();
52  if ( rc != LDNDC_ERR_OK)
53  {
54  KLOGERROR("Handling flooding event failed in: ", name());
55  return rc;
56  }
57 
58  /* tasks:
59  * - update internal state (e.g., throughfall, ...)
60  * - calculate evapotranspiration on demand
61  */
62  rc = EcHyStepInit();
63  if ( rc != LDNDC_ERR_OK)
64  { return rc; }
65 
66  /* begin water balance */
67  double balance( -1.0);
68  EcHyBalanceCheck( balance);
69 
70  /* snow and ice calculated */
71  EcHySnowIce();
72 
73  /*
74  * apply irrigation triggered by
75  * - irrigation events
76  * - flooding events
77  */
79 
80  double const potential_infiltration( surface_water + irrigation + throughfall);
81  double bypass_flow( m_param->CRACK_FRACTION() * potential_infiltration);
82 
83  /* bypass flow */
84  EcHyBypassFlow( m_param->CRACK_DEPTH(), bypass_flow);
85  double const infiltration( potential_infiltration - bypass_flow);
86 
87  /* downward water flow */
88  wlfl_sl[0] = cbm::bound_max( infiltration, kst_sl[0]);
89  for ( size_t sl = 0; sl < soillayers_in->soil_layer_cnt(); ++sl)
90  {
91  EcHyPercolation( 1, sl);
92  }
93 
94  if ( cbm::is_valid( kst_bottom))
95  {
96  wlfl_sl[soillayers_in->soil_layer_cnt()] = cbm::bound_max( wlfl_sl[soillayers_in->soil_layer_cnt()],
97  kst_bottom);
98  }
99  if ( cbm::is_valid( kst_bottom_gw_lateral))
100  {
101  wlfl_sl[soillayers_in->soil_layer_cnt()] = cbm::bound_max( wlfl_sl[soillayers_in->soil_layer_cnt()],
102  kst_bottom_gw_lateral);
103  }
104  if ( cbm::is_valid( kst_bottom_gw_vertical))
105  {
106  wlfl_sl[soillayers_in->soil_layer_cnt()] = cbm::bound_max( wlfl_sl[soillayers_in->soil_layer_cnt()],
107  kst_bottom_gw_vertical);
108  }
109 
111 
112  /* upward water flow */
113  for ( int sl = soillayers_in->soil_layer_cnt()-1; sl >= 0; --sl)
114  {
115  EcHyPercolation( 2, sl);
116  }
117 
118  /* runoff */
119  if ( runoff_method == "curvenumber")
120  {
121  double const surface_water_in_mm( surface_water * cbm::MM_IN_M);
122  double const curvenumber( 70.0);
123 
124  //retention parameter [mm]
125  double const retention_factor( 25.4 * (1000.0 / curvenumber - 10));
126 
127  //runoff [mm]
128  double const q_surf( cbm::sqr( surface_water_in_mm) /
129  ((surface_water_in_mm + retention_factor) * lclock()->time_resolution()));
130 
131  runoff = cbm::bound( 0.0,
132  q_surf * cbm::M_IN_MM,
133  0.99 * (surface_water + throughfall + irrigation - wlfl_sl[0] - bypass_flow));
134  }
135  else
136  {
137  double const f_runoff( cbm::bound_max( m_param->FRUNOFF() / lclock()->time_resolution(), 1.0));
138  runoff = cbm::bound_min( 0.0, f_runoff * (surface_water + throughfall + irrigation - wlfl_sl[0] - bypass_flow - bund_height));
139  }
140 
141  /* groundwater interaction:
142  * - set soil layers water-saturated
143  * - capillary rise
144  */
145  EcHyGroundwater();
146 
147  /* integration of state variables */
148  EcHyIntegration();
149 
150  /* update external state */
151  EcHyStepExit();
152 
153  /* perform water balance */
154  rc = EcHyBalanceCheck( balance);
155  if ( rc){ return rc; }
156 
157  rc = EcHy_check_for_negative_value("exit");
158  if ( rc){ return rc; }
159 
160  return LDNDC_ERR_OK;
161 }
lvector_t< double > kst_sl
Definition: echy.h:110
lerr_t EcHySnowIce()
Calls SnowDNDC for the calculation of snowpack and soil ice formation.
lerr_t EcHyFlood()
sets hydrologic conditions during flooding events, e.g.,
Definition: echy-management.cpp:17
void EcHyreset()
...
Definition: echy.cpp:170
lerr_t EcHyBypassFlow(double, double &)
Calculates water percolation within the soil profile.
lerr_t EcHyStepInit()
...
Definition: echy.cpp:202
lerr_t EcHyIntegration()
...
Definition: echy.cpp:430
lerr_t EcHyEvapotranspiration()
Calculates evapotranspiration within the soil profile.
cbm::string_t runoff_method
Definition: echy.h:194
double kst_bottom
Definition: echy.h:67
double runoff
Definition: echy.h:134
lvector_t< double > wlfl_sl
Definition: echy.h:89
double throughfall
Definition: echy.h:140
lerr_t EcHyGroundwater()
...
Definition: echy-groundwater.cpp:20
lerr_t EcHyStepExit()
...
Definition: echy.cpp:305
lerr_t EcHyPercolation(size_t, size_t)
Calculates water percolation within the soil profile.
lerr_t EcHyBalanceCheck(double &)
...
Definition: echy.cpp:453
lerr_t EcHyIrrigation()
Irrigation.
Definition: echy-management.cpp:74
Here is the call graph for this function:

Member Data Documentation

◆ accumulated_irrigation_old

double ldndc::EcHy::accumulated_irrigation_old
private

Accumulated irrigation of last time step [m]

Referenced by EcHyIrrigation().

◆ accumulated_potentialtranspiration_old

double ldndc::EcHy::accumulated_potentialtranspiration_old
private

Accumulated potential transpiration of last time step [m]

Referenced by EcHyPotentialEvapotranspiration().

◆ bypass_fill_sl

lvector_t< double > ldndc::EcHy::bypass_fill_sl
private

Soillayer water input due to bypass flow

Referenced by EcHyIntegration(), and EcHyreset().

◆ cr_fill_groundwater

double ldndc::EcHy::cr_fill_groundwater
private

Water addition to last soil layer from groundwater water by cappillary rise (only used for water balance calculation)

Referenced by EcHyBalanceCheck(), EcHyGroundwater(), and EcHyreset().

◆ cr_fill_sl

lvector_t< double > ldndc::EcHy::cr_fill_sl
private

Soillayer water input due to capillary rise

Referenced by EcHyGroundwater(), EcHyreset(), and EcHySoilWaterChange().

◆ daily_potential_evapotranspiration

double ldndc::EcHy::daily_potential_evapotranspiration
private

◆ daily_potential_leaf_evaporation

double ldndc::EcHy::daily_potential_leaf_evaporation
private

◆ daily_potential_soil_evaporation

double ldndc::EcHy::daily_potential_soil_evaporation
private

◆ daily_potential_transpiration

double ldndc::EcHy::daily_potential_transpiration
private

◆ ev_leaf

double ldndc::EcHy::ev_leaf
private

Evaporation from plant surface

Referenced by EcHyBalanceCheck(), EcHyPotentialEvapotranspiration(), and EcHyStepExit().

◆ evapotranspiration_method

cbm::string_t ldndc::EcHy::evapotranspiration_method
private

◆ evsws_sl

lvector_t< double > ldndc::EcHy::evsws_sl
private

◆ gw_depth_static

double ldndc::EcHy::gw_depth_static
private

Depth of groundwater table [m]

Referenced by EcHyGroundwater(), and EcHyStepInit().

◆ gw_fill_sl

lvector_t< double > ldndc::EcHy::gw_fill_sl
private

Soillayer water input due to groundwater flow

Referenced by EcHyBalanceCheck(), EcHyGroundwater(), EcHyreset(), and EcHySoilWaterChange().

◆ gw_fill_surface

double ldndc::EcHy::gw_fill_surface
private

Water addition to surface water from groundwater water boundary condition

Referenced by EcHyBalanceCheck(), EcHyGroundwater(), EcHyIntegration(), and EcHyreset().

◆ interception_water

double ldndc::EcHy::interception_water
private

◆ kst_bottom

double ldndc::EcHy::kst_bottom
private

Saturated hydraulic conductivity below last soil layer [cm:min-1]

Referenced by EcHyFlood(), EcHyStepInit(), and solve().

◆ kst_sl

lvector_t< double > ldndc::EcHy::kst_sl
private

Saturated hydraulic conductivity

Referenced by EcHyFlood(), EcHyGroundwater(), EcHySubsl2(), and solve().

◆ m_icecontent_in

WaterCycleSnowDNDC::IceContentStateIn ldndc::EcHy::m_icecontent_in
private

...

◆ runoff

double ldndc::EcHy::runoff
private

Lateral surface runoff

Referenced by EcHyBalanceCheck(), EcHyIntegration(), EcHyreset(), EcHyStepExit(), and solve().

◆ runoff_method

cbm::string_t ldndc::EcHy::runoff_method
private

...

Referenced by solve().

◆ snowfall

double ldndc::EcHy::snowfall
private

Snowfall

Referenced by EcHyBalanceCheck(), and EcHyreset().

◆ thornthwaite_heat_index

double ldndc::EcHy::thornthwaite_heat_index
private

◆ throughfall

double ldndc::EcHy::throughfall
private

Throughfall to soil

Referenced by EcHyBalanceCheck(), EcHyIntegration(), EcHyStepExit(), EcHyStepInit(), and solve().

◆ trwl_sl

lvector_t< double > ldndc::EcHy::trwl_sl
private

Water withdrawl by transpiration [m d−1]

Referenced by EcHyBalanceCheck(), EcHyPotentialEvapotranspiration(), EcHyreset(), EcHySoilWaterChange(), and EcHyStepExit().

◆ wc_res_sl

lvector_t< double > ldndc::EcHy::wc_res_sl
private

Minimum water content

Referenced by EcHyGetMinimumWater(), EcHyGroundwater(), and EcHyIntegration().

◆ wc_sat_sl

lvector_t< double > ldndc::EcHy::wc_sat_sl
private

Maximum water content

Referenced by EcHyGroundwater(), and EcHyIntegration().

◆ wl_sl

lvector_t< double > ldndc::EcHy::wl_sl
private

◆ wlfc_sl

lvector_t< double > ldndc::EcHy::wlfc_sl
private

Amount of water in soil layer at field capacity

Referenced by EcHyGetAvailableWaterTranspiration(), EcHyGetWaterLimitationTranspiration(), and EcHyIrrigation().

◆ wlfl_sl

lvector_t< double > ldndc::EcHy::wlfl_sl
private

Water flux at boundaries of soil layer [m d-1] the flux of water between soil layer is tracked by wlfl (mm d−1). in total, there are nl+1 flow rates, where nl is the number of soil layers, wlfl[1] is the flow rate between the ponded water layer and the soil surface, wlfl[2] is the flow rate between soil layer 1 and soil layer 2, etc.

Referenced by EcHyBalanceCheck(), EcHyIntegration(), EcHyreset(), EcHySoilWaterChange(), EcHyStepExit(), and solve().

◆ wlst_sl

lvector_t< double > ldndc::EcHy::wlst_sl
private

Amount of water in soil layer at saturation

Referenced by EcHyGroundwater(), and EcHyIrrigation().

◆ wlwp_sl

lvector_t< double > ldndc::EcHy::wlwp_sl
private

Amount of water in soil layer at wilting point

Referenced by EcHyGetWiltingPoint().