ldndc::WatercycleDNDC Class Reference

Watercycle model WatercycleDNDC. More...

Inherits ldndc::MBE_LegacyModel.

Public Member Functions
lerr_t	solve ()
	Kicks off computation for one time step.

Static Public Attributes
static const unsigned int	IMAX_W = 24

Private Member Functions
lerr_t	event_irrigate ()
	considers water input due to irrigation and/or liquid manure evetns
lerr_t	event_flood ()
	sets hydrologic conditions during flooding events, e.g., More...
double	balance_check (double *)
	checks each time step for correct water balance More...
void	CalcRaining (unsigned int)
	apportion daily rainfall to subdaily rainfall
double	rainfall_intensity ()
void	CalcIrrigation (unsigned int)
	apportion daily irrigation to subdaily irrigation More...
lerr_t	CalcIceContent ()
	Call of ice-content related processes.
lerr_t	CalcPotEvapoTranspiration ()
	Potential evaporation. More...
void	CalcInterception ()
	Calculates water quantity that is retained on leaf and stem bark surfaces. More...
double	get_leaf_water ()
	Collects leaf water from all canopy layers. More...
double	get_impedance_factor (size_t)
	Reduces water fluxes out of a layer if ice lenses have formed.
double	get_clay_limitation (size_t, double, double)
double	rootdensityfortranspiration_sl (double, double, double, bool)
lerr_t	CalcTranspiration ()
lerr_t	CalcTranspirationCouvreur (double const &hr_potTransinm)
void	CalcSnowEvaporation ()
void	CalcSurfaceWaterEvaporation ()
void	CalcSurfaceFlux ()
void	WatercycleDNDC_preferential_flow ()
void	CalcSoilEvapoPercolation ()
void	SoilWaterEvaporation (size_t)
double	WaterFlowCapillaryRise (size_t, double const &)
	Similar to WaterFlowBelowFieldCapacity(...) but restricts water flow by water availability from soil layer below.
double	WaterFlowAboveFieldCapacity (size_t, double const &)
double	WaterFlowSaturated (size_t, double const &)
double	get_fsl (size_t)
	Returns factor accounting for soil layer depth. Originally, a constant soil layer depth of 0.02 m has been assumed. Due to flexible layer depth, an additional factor has been introduced that accounts for a higher resistance of larger layers.
double	get_time_step_duration ()
	Time fraction with regard to daily rates.

Private Attributes
double	thornthwaite_heat_index
	heat index for potential evaporation using approach after Thornthwaite

Detailed Description

Watercycle model WatercycleDNDC.

Member Function Documentation

balance_check()

double ldndc::WatercycleDNDC::balance_check ( double *  _balance )

private

checks each time step for correct water balance

Balance check of water.

{
    double  balance = get_leaf_water() + wc_.surface_ice + wc_.surface_water + plant_waterdeficit_compensation;
    for ( size_t sl = 0; sl < m_soillayers->soil_layer_cnt(); ++sl)
    {
        balance += (wc_.wc_sl[sl] + wc_.ice_sl[sl]) * sc_.h_sl[sl];
    }
    balance += (   - m_wc.accumulated_irrigation_event_executed
                   - wc_.accumulated_irrigation_automatic
                   - wc_.accumulated_irrigation_ggcmi
                   - wc_.accumulated_irrigation_reservoir_withdrawal
                   - wc_.accumulated_precipitation
                   - wc_.accumulated_groundwater_access
                   + wc_.accumulated_percolation
                   + wc_.accumulated_runoff
                   + wc_.accumulated_wateruptake_sl.sum()
                   + wc_.accumulated_interceptionevaporation
                   + wc_.accumulated_soilevaporation
                   + wc_.accumulated_surfacewaterevaporation);

    if ( _balance)
    {
        double const  balance_delta( std::abs( *_balance - balance));
        if ( balance_delta > cbm::DIFFMAX)
        {
            KLOGWARN( "Water leakage:  difference is ", balance - *_balance);
        }
    }

    for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
    {
        if ( cbm::flt_less( wc_.wc_sl[sl], 0.0))
        {
            KLOGWARN( "Negative water content of: ", wc_.wc_sl[sl], " in layer: ", sl, ". Water content set to zero.");
            wc_.wc_sl[sl] = 0.0;
        }
        if ( cbm::flt_less( wc_.ice_sl[sl], 0.0))
        {
            KLOGWARN( "Negative ice content of: ", wc_.ice_sl[sl], " in layer: ", sl, ". Ice content set to zero.");
            wc_.ice_sl[sl] = 0.0;
        }
    }

    return  balance;
}

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

{
    //reset irrigation
    m_wc.irrigation = 0.0;

    if ( cbm::flt_in_range_lu(0.0, m_cfg.automaticirrigation, 1.0))
    {
        for ( size_t l = 0;  l < m_soillayers->soil_layer_cnt();  ++l)
        {
            // trigger irrigation events as soon as drought stress appears
            double const wc_target( sc_.wc_wp_sl[l] + m_cfg.automaticirrigation * (sc_.wc_fc_sl[l] -  sc_.wc_wp_sl[l]));
            if ( cbm::flt_less( wc_.wc_sl[l], wc_target))
            {
                // irrigation assumed until target water content
                double const irrigate_layer( (wc_target - wc_.wc_sl[l]) * sc_.h_sl[l]);
                m_wc.irrigation += irrigate_layer;
                wc_.accumulated_irrigation_automatic += irrigate_layer;
            }
        }
    }

    if ( cbm::flt_in_range_lu( 0.0, m_eventflood.get_saturation_level(), 1.0) &&
         cbm::is_valid( m_eventflood.get_irrigation_height()) &&
         cbm::is_valid( m_eventflood.get_soil_depth()))
    {
        double wc_sum( 0.0);
        double wc_fc_sum( 0.0);
        double wc_wp_sum( 0.0);
        for ( size_t l = 0;  l < m_soillayers->soil_layer_cnt();  ++l)
        {
            wc_sum += wc_.wc_sl[l] * sc_.h_sl[l];
            wc_fc_sum += sc_.wc_fc_sl[l] * sc_.h_sl[l];
            wc_wp_sum += sc_.wc_wp_sl[l] * sc_.h_sl[l];

            if ( cbm::flt_greater_equal( sc_.depth_sl[l], m_eventflood.get_soil_depth()))
            {
                break;
            }
        }

        // trigger irrigation events as soon as drought stress appears
        double const wc_target( wc_wp_sum + m_eventflood.get_saturation_level() * (wc_fc_sum -  wc_wp_sum));
        if ( cbm::flt_less( wc_sum, wc_target))
        {
            m_wc.irrigation += m_eventflood.get_irrigation_height();
            wc_.accumulated_irrigation_automatic += m_eventflood.get_irrigation_height();
        }
    }

    if ( m_cfg.ggcmi_irrigation)
    {
        if( (lclock()->subday() == 1) && (_i == 0))
        {
            if( m_veg->size() > 0 ) //check crop is present
            {
                double water_supply( 0.0);
                for ( size_t l = 0;  l < m_soillayers->soil_layer_cnt();  ++l)
                {
                    if ( cbm::flt_greater_zero( m_veg->mfrt_sl( l))) //only irrigate soil layers that contain roots
                    {
                        if( cbm::flt_less( wc_.wc_sl[l], sc_.wc_fc_sl[l]) &&
                            cbm::flt_equal_zero( wc_.ice_sl[l])) //only add water if no ice is present
                        {
                            water_supply += (sc_.wc_fc_sl[l] - wc_.wc_sl[l]) * sc_.h_sl[l]; //amount of water added (m)
                            wc_.wc_sl[l] = sc_.wc_fc_sl[l];
                        }
                    }
                }
                wc_.accumulated_irrigation_ggcmi += water_supply;
            }
        }
    }

    if ( cbm::is_valid( m_eventflood.get_water_table()))
    {
        /* Dynamic flooding */
        if ( cbm::is_valid( m_eventflood.get_irrigation_height()))
        {
            bool irrigate( false);

            /* Irrigate after surface water table drop */
            if ( cbm::flt_greater_equal_zero( m_eventflood.get_water_table()))
            {
                if ( cbm::flt_greater( m_eventflood.get_water_table(), wc_.surface_water))
                {
                    irrigate = true;
                }
            }
            else
            {
                if ( !cbm::flt_greater_zero( wc_.surface_water))
                {
                    for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
                    {
                        if ( cbm::flt_greater_equal( sc_.depth_sl[sl], -m_eventflood.get_water_table()))
                        {
                            //check if water content below AWD_DEPTH is saturated
                            if ( cbm::flt_less( wc_.wc_sl[sl], 0.9 * sc_.poro_sl[sl]))
                            {
                                irrigate = true;
                            }
                            break;
                        }
                    }
                }
            }
            // trigger AWD irrigation event
            if ( irrigate)
            {
                double water_supply( cbm::bound_min( 0.0, m_eventflood.get_irrigation_height() - wc_.surface_water));
                for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
                {
                    // set soil fully water saturated
                    if ( cbm::flt_greater( sc_.poro_sl[sl], wc_.wc_sl[sl]))
                    {
                        water_supply += ((sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl]);
                    }
                }

                //remove rainfall (less irrigation water is supplied if it if raining)
                water_supply = cbm::bound_min( 0.0, water_supply - m_mc.rainfall);

                if ( m_eventflood.have_unlimited_water())
                {
                    wc_.accumulated_irrigation_automatic += water_supply;
                }
                else
                {
                    if ( cbm::flt_greater( wc_.irrigation_reservoir, water_supply))
                    {
                        wc_.irrigation_reservoir -= water_supply;
                    }
                    else
                    {
                        water_supply = wc_.irrigation_reservoir;
                        wc_.irrigation_reservoir = 0.0;
                    }
                    wc_.accumulated_irrigation_reservoir_withdrawal += water_supply;
                }
                m_wc.irrigation += water_supply;
            }
        }
        /* Static flooding */
        else
        {
            if ( !cbm::flt_greater_zero( m_eventflood.get_water_table()))
            {
                double water_supply( 0.0);
                for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
                {
                    if ( cbm::flt_greater( sc_.depth_sl[sl], -m_eventflood.get_water_table()) &&
                         cbm::flt_greater( sc_.poro_sl[sl], wc_.wc_sl[sl]))
                    {
                        water_supply += (sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl];
                    }
                }

                //remove rainfall (less irrigation water is supplied if it if raining)
                water_supply = cbm::bound_min( 0.0, water_supply - m_mc.rainfall);

                /* Surface water will be actively drained */
                if ( cbm::flt_greater( wc_.surface_water, water_supply))
                {
                    wc_.accumulated_runoff += wc_.surface_water - water_supply;
                    wc_.surface_water = water_supply;
                }

                if ( m_eventflood.have_unlimited_water())
                {
                    wc_.accumulated_irrigation_automatic += water_supply;
                }
                else
                {
                    if ( cbm::flt_greater( wc_.irrigation_reservoir, water_supply))
                    {
                        wc_.irrigation_reservoir -= water_supply;
                    }
                    else
                    {
                        water_supply = wc_.irrigation_reservoir;
                        wc_.irrigation_reservoir = 0.0;
                    }
                    wc_.accumulated_irrigation_reservoir_withdrawal += water_supply;
                }
                m_wc.irrigation += water_supply;
            }
            else if ( cbm::flt_greater( m_eventflood.get_water_table(), wc_.surface_water))
            {
                //adjust surface and soil water in case of flooding event
                //water required to fill up surface water
                double water_supply( m_eventflood.get_water_table() - wc_.surface_water);
                
                //add water required to fill top x layers up to saturation
                unsigned long no_layers(3);
                if( no_layers > m_soillayers->soil_layer_cnt())
                {
                    no_layers = m_soillayers->soil_layer_cnt();
                }
                
                for ( size_t sl = 0;  sl < no_layers;  ++sl)
                {
                    water_supply += ((sc_.poro_sl[sl] - wc_.wc_sl[sl]) * sc_.h_sl[sl]);
                }
                
                //remove rainfall (less irrigation water is supplied if it if raining)
                water_supply = cbm::bound_min( 0.0, water_supply - m_mc.rainfall);
                
                if ( m_eventflood.have_unlimited_water())
                {
                    wc_.accumulated_irrigation_automatic += water_supply;
                }
                else
                {
                    if ( cbm::flt_greater( wc_.irrigation_reservoir, water_supply))
                    {
                        wc_.irrigation_reservoir -= water_supply;
                    }
                    else
                    {
                        water_supply = wc_.irrigation_reservoir;
                        wc_.irrigation_reservoir = 0.0;
                    }
                    wc_.accumulated_irrigation_reservoir_withdrawal += water_supply;
                }
                m_wc.irrigation += water_supply;
            }
        }
    }

    double const irrigation_scheduled( cbm::bound_min( 0.0, wc_.accumulated_irrigation - m_wc.accumulated_irrigation_event_executed));
    
    if ( cbm::flt_greater( irrigation_scheduled, 24.0 * rainfall_intensity()))
    {
        m_wc.irrigation += irrigation_scheduled / 24.0;
        m_wc.accumulated_irrigation_event_executed += irrigation_scheduled / 24.0;
    }
    else if ( cbm::flt_greater( irrigation_scheduled, (double)hours_per_time_step() * rainfall_intensity()))
    {
        m_wc.irrigation += (double)hours_per_time_step() * rainfall_intensity();
        m_wc.accumulated_irrigation_event_executed += (double)hours_per_time_step() * rainfall_intensity();
    }
    else
    {
        m_wc.irrigation += irrigation_scheduled;
        m_wc.accumulated_irrigation_event_executed += irrigation_scheduled;
    }
}

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

{
    day_potevapo = 0.0;
    if ( m_cfg.potentialevapotranspiration_method == "thornthwaite")
    {
        if ( lclock()->is_position( TMODE_PRE_YEARLY))
        {
            thornthwaite_heat_index = ldndc::thornthwaite_heat_index( m_setup->latitude(),
                                                                      lclock()->days_in_year(),
                                                                      m_climate->annual_temperature_average(),
                                                                      m_climate->annual_temperature_amplitude());
        }

        double const daylength(ldndc::meteo::daylength(m_setup->latitude(), lclock()->yearday()));
        double const avg_dayl(12.0);

        day_potevapo = m_param->WCDNDC_INCREASE_POT_EVAPOTRANS()
                            * ldndc::thornthwaite(
                                           mc_.nd_airtemperature,
                                           daylength, avg_dayl, lclock()->days_in_month(),
                                           thornthwaite_heat_index);
    }
    else if ( (m_cfg.potentialevapotranspiration_method == "penman") ||
              (m_cfg.potentialevapotranspiration_method == "priestleytaylor"))
    {
        /* clock */
        cbm::sclock_t const &  clk( this->lclock_ref());

        /* leaf area index */
        double lai( 0.0);
        for( PlantIterator  vt = m_veg->begin(); vt != m_veg->end(); ++vt)
        {
            lai += (*vt)->lai();
        }

        /* vapour pressure [10^3:Pa] */
        double vps( 0.0);
        ldndc::meteo::vps( mc_.nd_airtemperature, &vps);
        double vp( vps - m_climate->vpd_day( clk));

        if( m_cfg.potentialevapotranspiration_method == "priestleytaylor")
        {
            day_potevapo = m_param->WCDNDC_INCREASE_POT_EVAPOTRANS()
                                * ldndc::priestleytaylor( clk.yearday(),
                                                         clk.days_in_year(),
                                                         mc_.albedo,
                                                         m_setup->latitude(),
                                                         mc_.nd_airtemperature,
                                                         mc_.nd_shortwaveradiation_in,
                                                         vp,
                                                         m_param->PT_ALPHA());
        }
        else if( m_cfg.potentialevapotranspiration_method == "penman")
        {
            ldndc::surface_type surface( cbm::flt_greater_zero( lai) ? short_grass : bare_soil);

            day_potevapo = m_param->WCDNDC_INCREASE_POT_EVAPOTRANS()
                                * ldndc::penman(  surface,
                                                clk.yearday(),
                                                clk.days_in_year(),
                                                mc_.albedo,
                                                m_setup->latitude(),
                                                mc_.nd_shortwaveradiation_in * cbm::SEC_IN_DAY,
                                                mc_.nd_airtemperature,
                                                mc_.nd_windspeed,
                                                vp);
        }
    }
    else
    {
        KLOGERROR("[BUG] huh :o how did you get here? ",
                  "tell the maintainers refering to this error message");
        return  LDNDC_ERR_FAIL;
    }

    if ( m_cfg.evapotranspiration_method == "uniform")
    {
        // hourly evaporation demand is calculated from daily demand equally throughout the day
        hr_pot_evapotrans = day_potevapo * nhr / cbm::HR_IN_DAY;
    }
    else if ( m_cfg.evapotranspiration_method == "previouspattern")
    {

        hr_pot_evapotrans = day_potevapo * pot_evapotranspiration_ts[lclock()->subday()-1] / nwc;
    }
    else //if ( m_cfg.evapotranspiration_method === "nice name") // could be for example 'daytime_only'
    {
        double const hour(lclock()->subday());
        double const day(lclock()->yearday());
        double const hsun(0.833 * cbm::PI / 180.0);
        double const lat(m_setup->latitude() * cbm::PI / 180.0);
        double const dec(-23.44 * (cbm::PI / 180.0) * std::cos(2.0 * cbm::PI * (day - 172.0) / lclock()->days_in_year()));
        double const help((std::sin(hsun) - std::sin(lat) * std::sin(dec)) / (std::cos(lat) * std::cos(dec)));
        double dayl(0.0);
        if (help < -0.999) dayl = 0.0;
        else if (help > 0.999) dayl = 24.0;
        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)));
        if (dayl < 1.0) dayl = 0.0;
        double const rise(12.0 - dayl * 0.5);

        double factor(0.0);
        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);

        double fday(0.0);
        if (hour > rise && hour < (24.0 - rise)) fday = factor;

        hr_pot_evapotrans = day_potevapo * fday;
    }
    // the evaporation limit for all timesteps is the daily demand
    wc_.accumulated_potentialevaporation += hr_pot_evapotrans;

    return  LDNDC_ERR_OK;
}

Here is the call graph for this function:

CalcSnowEvaporation()

void ldndc::WatercycleDNDC::CalcSnowEvaporation ( )

private

Parameters

[in]	None
[out]	None

Returns

None

{
    //Average limit of snowEvaporation is 0.00013m day-1 (Granger and Male 1978)
    double hr_potESnow( cbm::bound(0.0, hr_pot_evapotrans, 0.00013 * nhr/24.0));

    if ( cbm::flt_greater_zero( wc_.surface_ice) &&
         cbm::flt_greater_zero( hr_potESnow))
    {
        //If there is a snowpack and there is a potential to evaporate
        if ( wc_.surface_ice > hr_potESnow)
        {
            //Partial evaporation from the snow
            hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - hr_potESnow);
            wc_.accumulated_surfacewaterevaporation += hr_potESnow;
            wc_.surface_ice -= hr_potESnow;
        }
        else
        {
            //Total evaporation of the snowpack
            hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - wc_.surface_ice);;
            wc_.accumulated_surfacewaterevaporation += wc_.surface_ice;
            wc_.surface_ice = 0.0;
        }
    }
}

CalcSoilEvapoPercolation()

void ldndc::WatercycleDNDC::CalcSoilEvapoPercolation ( )

private

Parameters

[in]	None
[out]	None

Returns

None

{
    /* Reset water fluxes */
    for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
    {
        m_wc.percolation_sl[sl] = 0.0;
    }

    /* Fill groundwater layer with water */
    for ( int sl = m_soillayers->soil_layer_cnt()-1;  sl >= 0; --sl)
    {
        if ( cbm::flt_greater_equal( sc_.depth_sl[sl], ground_water_table))
        {
            /* ensure that ice volume is less than porosity */
            double const ice_vol( wc_.ice_sl[sl] / cbm::DICE);
            if ( cbm::flt_greater( sc_.poro_sl[sl], ice_vol))
            {
                //air filled porosity is filled with groundwater
                double const gw_access( cbm::bound_min( 0.0,
                                                        (sc_.poro_sl[sl] - ice_vol - wc_.wc_sl[sl]) * sc_.h_sl[sl]));
                wc_.wc_sl[sl] += gw_access / sc_.h_sl[sl];
                wc_.accumulated_groundwater_access += gw_access;
            }
        }
        else
        {
            break;
        }
    }

    for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
    {
        if ( sl == 0)
        {
            wc_.accumulated_infiltration += wc_.surface_water;
            wc_.wc_sl[0] += wc_.surface_water / sc_.h_sl[0];
            wc_.surface_water = 0.0;
        }
        else
        {
            wc_.wc_sl[sl] += m_wc.percolation_sl[sl-1] / sc_.h_sl[sl];
        }

        if ( cbm::flt_greater( sc_.depth_sl[sl], ground_water_table))
        {
            //set constant flux for all ground water layers depending on last non-groundwater layer
            if ( sl > 0)
            {
                m_wc.percolation_sl[sl] = m_wc.percolation_sl[sl-1];
            }

            //in case of groundwater until soil surface first outflux is estimated by saturated water flow
            else
            {
                m_wc.percolation_sl[sl] = WaterFlowSaturated( sl, get_impedance_factor( sl));
            }
        }
        else if ( sl + 1 == m_soillayers->soil_layer_cnt())
        {
            m_wc.percolation_sl[sl] = 0.0;
            if ( cbm::flt_greater_zero( wc_.sks_sl[sl]))
            {
                if ( cbm::flt_greater( wc_.wc_sl[sl], sc_.poro_sl[sl]) &&
                     cbm::flt_greater( wc_.wc_sl[sl], sc_.wc_fc_sl[sl])) //security, should be always true if first condition is met
                {
                    //maximum allowed flow until field capacity
                    double const max_flow( (wc_.wc_sl[sl] - sc_.wc_fc_sl[sl]) * sc_.h_sl[sl]);
                    m_wc.percolation_sl[sl] = cbm::bound_max( WaterFlowSaturated( sl, get_impedance_factor( sl)),
                                                             max_flow);
                }

                //water flux above field capacity
                else if ( cbm::flt_greater_equal( wc_.wc_sl[sl] + wc_.ice_sl[sl], sc_.wc_fc_sl[sl]))
                {
                    m_wc.percolation_sl[sl] = WaterFlowAboveFieldCapacity( sl, get_impedance_factor( sl));
                }

                //water flux below field capacity and above wilting point
                else if (   cbm::flt_greater_zero( m_wc.percolation_sl[sl-1])
                         && cbm::flt_greater( wc_.wc_sl[sl], sc_.wc_wp_sl[sl]))
                {
                    m_wc.percolation_sl[sl] = m_param->FPERCOL() * m_wc.percolation_sl[sl-1];
                }

                //bound water flux by defined maximal percolation rate (e.g., during flooding events)
                m_wc.percolation_sl[sl] = cbm::bound_max( m_wc.percolation_sl[sl], max_percolation);
            }
        }

        //infiltrating water in current time step is exceeding available pore space of last time step
        //second condition: security, should be always true if first condition is met
        else if ( cbm::flt_greater( wc_.wc_sl[sl], sc_.poro_sl[sl]) &&
                  cbm::flt_greater( wc_.wc_sl[sl], sc_.wc_fc_sl[sl]))
        {
            //maximum allowed flow until field capacity
            double const max_flow( (wc_.wc_sl[sl] - sc_.wc_fc_sl[sl]) * sc_.h_sl[sl]);
            m_wc.percolation_sl[sl] = cbm::bound_max( WaterFlowSaturated( sl, get_impedance_factor( sl)),
                                                      max_flow);
        }
        else if ( cbm::flt_greater_equal( wc_.wc_sl[sl] + wc_.ice_sl[sl], sc_.wc_fc_sl[sl]))
        {
            m_wc.percolation_sl[sl] = WaterFlowAboveFieldCapacity( sl, get_impedance_factor( sl));
        }
        else if ( cbm::flt_greater( wc_.wc_sl[sl], sc_.wc_wp_sl[sl]))
        {
            m_wc.percolation_sl[sl] = WaterFlowCapillaryRise( sl, get_impedance_factor( sl));
        }
        else
        {
            m_wc.percolation_sl[sl] = 0.0;
        }


        // water content must be greater or equal wilting point
        if ( cbm::flt_less( (wc_.wc_sl[sl] + wc_.ice_sl[sl]) * sc_.h_sl[sl] - m_wc.percolation_sl[sl],
                             sc_.wc_wp_sl[sl] * sc_.h_sl[sl]) &&
             cbm::flt_greater_zero( m_wc.percolation_sl[sl]))
        {
            m_wc.percolation_sl[sl] = cbm::bound( 0.0,
                                                  (wc_.wc_sl[sl] + wc_.ice_sl[sl] - sc_.wc_wp_sl[sl]) * sc_.h_sl[sl],
                                                  m_wc.percolation_sl[sl]);
        }

        // Update of the current soil layer water content
        wc_.wc_sl[sl] -= m_wc.percolation_sl[sl] / sc_.h_sl[sl];

        SoilWaterEvaporation( sl);
    }

    /* Correct water content and percolation if water content is above porosity
     * If the layer below cannot take up all of the water assigned for percolation,
     * so if the air filled pore space is not sufficient,
     * the water remains in the upper layer.
     */
    for ( size_t sl = m_soillayers->soil_layer_cnt()-1;  sl > 0; sl--)
    {
        double const ice_vol( wc_.ice_sl[sl] / cbm::DICE);
        if ( cbm::flt_greater( wc_.wc_sl[sl] + ice_vol,
                               sc_.poro_sl[sl]))
        {
            double delta_wc( wc_.wc_sl[sl] + ice_vol - sc_.poro_sl[sl]);
            if ( cbm::flt_greater_equal( ice_vol, sc_.poro_sl[sl]))
            {
                delta_wc = wc_.wc_sl[sl];
                wc_.wc_sl[sl] = 0.0;
            }
            else
            {
                wc_.wc_sl[sl] = sc_.poro_sl[sl] - ice_vol;
            }
            m_wc.percolation_sl[sl-1] -= delta_wc * sc_.h_sl[sl];
            wc_.wc_sl[sl-1] += delta_wc * sc_.h_sl[sl] / sc_.h_sl[sl-1];
        }
    }

    double const ice_vol( wc_.ice_sl[0] / cbm::DICE);
    if ( cbm::flt_greater( wc_.wc_sl[0] + ice_vol,
                           sc_.poro_sl[0]))
    {

        double delta_wc( wc_.wc_sl[0] + ice_vol - sc_.poro_sl[0]);
        if ( cbm::flt_greater_equal( ice_vol, sc_.poro_sl[0]))
        {
            delta_wc = wc_.wc_sl[0];
            wc_.wc_sl[0] = 0.0;
        }
        else
        {
            wc_.wc_sl[0] = sc_.poro_sl[0] - ice_vol;
        }
        wc_.accumulated_infiltration -= delta_wc * sc_.h_sl[0];
        wc_.surface_water += delta_wc * sc_.h_sl[0];
    }


    CalcSurfaceFlux();

    wc_.accumulated_waterflux_sl += m_wc.percolation_sl;
    wc_.accumulated_percolation += m_wc.percolation_sl[ m_soillayers->soil_layer_cnt()-1];
}

CalcSurfaceFlux()

void ldndc::WatercycleDNDC::CalcSurfaceFlux ( )

private

Parameters

[in]	None
[out]	None

Returns

None

{
    if ( cbm::flt_greater_zero( m_eventflood.get_bund_height()))
    {
        if ( cbm::flt_greater( wc_.surface_water, m_eventflood.get_bund_height()))
        {
            wc_.accumulated_runoff += wc_.surface_water - m_eventflood.get_bund_height();
            wc_.surface_water = m_eventflood.get_bund_height();
        }
    }
    else
    {
        /* runoff fraction per time step equals RO fraction per hour * length of a time step */
        double const runoff_fraction( m_param->FRUNOFF() * get_time_step_duration());

        // if frunoff_ts < 1
        if ( cbm::flt_less( runoff_fraction, 1.0))
        {
            double const runoff( runoff_fraction * wc_.surface_water);
            wc_.surface_water -= runoff;
            wc_.accumulated_runoff += runoff;
        }
        else
        {
            wc_.accumulated_runoff += wc_.surface_water;
            wc_.surface_water = 0.0;
        }
    }
}

CalcSurfaceWaterEvaporation()

void ldndc::WatercycleDNDC::CalcSurfaceWaterEvaporation ( )

private

Parameters

[in]	None
[out]	None

Returns

None

{
    if ( cbm::flt_greater_zero( wc_.surface_water) &&
         cbm::flt_greater_zero( hr_pot_evapotrans))
    {
        if ( hr_pot_evapotrans > wc_.surface_water)
        {
            wc_.accumulated_surfacewaterevaporation += wc_.surface_water;
            hr_pot_evapotrans -= wc_.surface_water;
            wc_.surface_water = 0.0;
        }
        else
        {
            wc_.surface_water -= hr_pot_evapotrans;
            wc_.accumulated_surfacewaterevaporation += hr_pot_evapotrans;
            hr_pot_evapotrans = 0.0;
        }
    }
}

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

{
        // hourly potential transpiration
        double const hr_potTrans( hr_potTransinm * cbm::CM_IN_M);
        // hourly transpiration
        double hr_uptake( 0.0);

        // Calculate the root water uptake and hence the transpiration for every plant individually
        for( PlantIterator vt = m_veg->begin(); vt != m_veg->end(); ++vt)
        {
                // split up the potential transpiration equally to the various plants
                // frt mass of the current species
                double const mfrt_species( (*vt)->mFrt);
                // total frt mass
                double const mfrt_sum( m_veg->dw_frt());
                double const hr_potTransspecies = hr_potTrans * mfrt_species/mfrt_sum;  // in cm

                // If some soil layers are frozen and contain only ice, no fluid water is present and the water potential diverges.
                // The uptake from these layers is prevented by considering only roots in layers with more water than ice.
                // In this case the root distribution needs to be rescaled (otherwise the non-contributing layers would effectively contribute a potential = 0).
                double fFrt_unfrozen( 1.0);

                /* TODO: include icetowatercrit as proper parameter */
                double icetowatercrit( 0.1);
                double effsoilwatpot( 0.0);
                for( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
                {
                        // fine root conditions to prevent transpiration without uptake organs, wateruptake from frozen ground
                        if( cbm::flt_greater_zero( (*vt)->fFrt_sl[sl] * mfrt_species))
                        {
                                if(cbm::flt_less( wc_.ice_sl[sl]/wc_.wc_sl[sl], icetowatercrit))
                                {
                                        // overall water potential in this layer
                                        double const watpotl( Soillayerwaterhead( sl)); // negative
                                        effsoilwatpot += watpotl * (*vt)->fFrt_sl[sl];  // negative
                                }
                                else
                                {
                                        fFrt_unfrozen -= (*vt)->fFrt_sl[sl];
                                }
                        }
                        else
                        {
                                // roots are not wireless ;-)
                                break;
                        }
                }
                // normalize the root distribution in the potential, if some soil is frozen
                effsoilwatpot = effsoilwatpot / fFrt_unfrozen;

                // hydraulic conductivities of plant system
                double const Krs(vt->KRSINIT() * mfrt_species / vt->FRTMASSINIT() * nhr);               // cm(water)/cm(from pressure)/timestep
                double const Kcomp(vt->KCOMPINIT() * mfrt_species / vt->FRTMASSINIT() * nhr);   // cm(water)/cm(from pressure)/timestep

                // water potential in the leafs, depending on the potential transpiration
                double const leafwatpot( cbm::bound_min(vt->HLEAFCRIT(), effsoilwatpot - (hr_potTransspecies / Krs)));          // cm

                if( cbm::flt_greater_zero( leafwatpot))
                {
                        KLOGERROR( "leafwatpot is positive");
            return  LDNDC_ERR_FAIL;
                }
                if( !cbm::flt_greater_zero( effsoilwatpot - leafwatpot))
                {
                        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.");
            return  LDNDC_ERR_FAIL;
                }

                double const effuptake( Krs * (effsoilwatpot - leafwatpot));            // cm/timestep
                double uptake_tot( 0.0);
                double uptakecomp_tot( 0.0);
//              double uptakecompabs_tot( 0.0);

                for( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
                {
                        // fine root condition to prevent transpiration without uptake organs
                        if( cbm::flt_greater_zero( (*vt)->fFrt_sl[sl] * mfrt_species))
                        {
                                if(cbm::flt_less( wc_.ice_sl[sl]/wc_.wc_sl[sl], icetowatercrit))
                                {
                                        double const watpotl( Soillayerwaterhead( sl));
                                        double const compuptake( Kcomp * (watpotl - effsoilwatpot));
                                        double uptake_layer( (effuptake + compuptake) * (*vt)->fFrt_sl[sl]);    // absolute value in cm

                    wc_.accumulated_wateruptake_sl[sl] += uptake_layer * cbm::M_IN_CM;
                                        uptake_tot += uptake_layer;                                                                                             // absolute value in cm
                                        uptakecomp_tot += compuptake * (*vt)->fFrt_sl[sl];                                              // absolute value in cm
//                                      uptakecompabs_tot += fabs(compuptake * (*vt)->fFrt_sl[sl]);                             // absolute value in cm

                                        wc_.wc_sl[sl] -= uptake_layer * cbm::M_IN_CM /sc_.h_sl[sl];             // fraction of water in this layer, in m

                                        if( !cbm::flt_greater_zero( wc_.wc_sl[sl]))
                                        {
                                                KLOGERROR( "In soil layer ", sl, " the water content is negative: ", wc_.wc_sl[sl]);
                        return  LDNDC_ERR_FAIL;
                                        }
                                }
                        }
                        else
                        {
                                // roots are not wireless ;-)
                                break;
                        }
                }

        hr_uptake += uptake_tot * cbm::M_IN_CM;                                                 // absolute value in m for all species


                if( cbm::flt_greater_zero( fabs(uptakecomp_tot))){
                        KLOGWARN( "Compensatory water uptake ", uptakecomp_tot, " > 0.");
                }
                if( cbm::flt_greater( uptake_tot, hr_potTransspecies, 0.0000000001))
                {
                        KLOGERROR( "Total water uptake ", uptake_tot, " larger than potential transpiration ", hr_potTransspecies);
                }
                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))
                {
                        KLOGERROR( "T ", uptake_tot, "cm < Tpot ", hr_potTransspecies, "cm even though leafwatpot ", leafwatpot, " > HLEAFCRIT", vt->HLEAFCRIT());
                }
        }

    // total evaporation demand reduced to apply on other evaporating sources (ground, snow, surface water)
    hr_pot_evapotrans = cbm::bound_min(0.0, hr_pot_evapotrans - hr_uptake);

    return  LDNDC_ERR_OK;
}

event_flood()

lerr_t ldndc::WatercycleDNDC::event_flood ( )

private

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

• surface water table
• bund height

{
    lerr_t  rc = m_eventflood.solve();
    if ( rc != LDNDC_ERR_OK)
    {
        KLOGERROR("Irrigation event not successful!");
        return  rc;
    }

    /* max_percolation is gradually decreased after end of flooding event to inital value */
    double const maximum_percolation = m_eventflood.get_maximum_percolation();
    if ( cbm::is_valid( maximum_percolation))
    {
        if ( cbm::flt_greater_zero( maximum_percolation))
        {
            max_percolation = maximum_percolation * nhr / 24.0;
        }
    }
    else
    {
        /* time rate for gradual recovery of max_percolation */
        double const time_rate( get_time_step_duration() * 0.1);
        max_percolation -= (max_percolation - WaterFlowSaturated( m_soillayers->soil_layer_cnt()-1, 1.0)) * time_rate;
    }

    return  LDNDC_ERR_OK;
}

get_clay_limitation()

double ldndc::WatercycleDNDC::get_clay_limitation ( size_t  _sl, double  _water_avail, double  _water_ref  )

private

Parameters

[in]	_sl	Soil layer
[in]	_water_avail	Available water content
[in]	_water_ref	Reference water content
[out]	None

Returns

Clay limitation

{
    // layer specific water availability
    // organic fraction  [%]
    double const orgf( sc_.fcorg_sl[_sl] / cbm::CCORG);
    // clay factor "1"  [-]
    double const clay_f1( pow(10.0, (-sc_.clay_sl[_sl] - m_param->RCLAY())));
    // clay factor "2"  [-]
    double const clay_f2( std::max(0.0, 1.0 + std::min(1.0 - orgf, sc_.clay_sl[_sl])));

    // fraction of max uptake that can be used according to soil conditions  [%]
    return cbm::bound(0.0, clay_f1 + pow( _water_avail / _water_ref, clay_f2), 1.0);
}

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

{
    if ( m_veg->size() > 0u)
    {
        return wc_.wc_fl.sum();
    }
    return 0.0;
}

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

{
    if ( m_iokcomm->get_input_class< input_class_climate_t >())
    {
        return m_iokcomm->get_input_class< climate::input_class_climate_t >()->station_info()->rainfall_intensity * cbm::M_IN_MM;
    }
    else
    {
        return ldndc::climate::climate_info_defaults.rainfall_intensity * cbm::M_IN_MM;
    }
}

Here is the call graph for this function:

rootdensityfortranspiration_sl()

double ldndc::WatercycleDNDC::rootdensityfortranspiration_sl ( double  m_sl, double  l_sl, double  h_sl, bool  root_length_h2o_up  )

private

Parameters

[in]	m_sl	root mass in soil layer
[in]	l_sl	root length in soil layer
[in]	h_sl	height of soil layer
[in]	root_length_h2o_up,determines	if root mass or length is used

Returns

root density in this soil layer, either by mass or length

{
    double rootdensity_sl = -99.99;
    if( root_length_h2o_up == true)
    {
        rootdensity_sl = l_sl / h_sl;
    }
    else
    {
        rootdensity_sl = m_sl / h_sl;
    }
        
    return rootdensity_sl;
}

SoilWaterEvaporation()

void ldndc::WatercycleDNDC::SoilWaterEvaporation ( size_t  _sl )

private

Parameters

[in]	_sl	Soil layer
[out]	None

Returns

None

{
    //percolation and direct evaporation from the (upper) soil
    //if there is snow cover, evaporation from soil will be set to 0, and only sublimation is allowed.
    double const wc_min_evaporation( m_param->WCDNDC_EVALIM_FRAC_WCMIN() * sc_.wc_wp_sl[_sl]);
    if ( cbm::flt_less( sc_.depth_sl[_sl] - sc_.h_sl[_sl], m_param->EVALIM()) &&
         cbm::flt_greater( wc_.wc_sl[_sl], wc_min_evaporation) &&
         cbm::flt_greater_zero( hr_pot_evapotrans))
    {
        //limiting factor for soil water infiltration
        double const f_clay( 1.0 + m_param->SLOPE_CLAYF() * std::min(1.0 - sc_.fcorg_sl[_sl] / cbm::CCORG, sc_.clay_sl[_sl]));
        double const f_water( cbm::bound( 0.0,
                                          pow( (wc_.wc_sl[_sl] - wc_min_evaporation) / sc_.wc_fc_sl[_sl], f_clay),
                                          1.0));
        double const h_rest( std::min( sc_.h_sl[_sl], m_param->EVALIM() - ( sc_.depth_sl[_sl] - sc_.h_sl[_sl])));
        double const f_depth( h_rest / m_param->EVALIM() * std::max(0.0,
                                       1.0 - std::min((double)m_param->EVALIM(),
                                                      sc_.depth_sl[_sl] - ( 0.5 * sc_.h_sl[_sl]))
                                       / m_param->EVALIM()));
        double const soilevaporation_max( cbm::bound_max( f_water * wc_.wc_sl[_sl] * sc_.h_sl[_sl], hr_pot_evapotrans));

        double const soilevaporation( cbm::bound( 0.0,
                                      hr_pot_evapotrans * f_depth * f_water,
                                      soilevaporation_max));

        //update of the current soil layer water content
        wc_.wc_sl[_sl] -= soilevaporation / sc_.h_sl[_sl];
        wc_.accumulated_soilevaporation += soilevaporation;
        hr_pot_evapotrans = cbm::bound_min( 0.0, hr_pot_evapotrans - soilevaporation);
    }
}

WatercycleDNDC_preferential_flow()

void ldndc::WatercycleDNDC::WatercycleDNDC_preferential_flow ( )

private

Parameters

[in]	None
[out]	None

Returns

None

{
    //BY_PASSF is defined as the fraction of surface water that passes per hour (0.0 by default)
    if ( cbm::flt_greater_zero( m_param->BY_PASSF()) &&
         cbm::flt_greater_zero( wc_.surface_water))
    {

        double water_def_total( 0.0);
        for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
        {
            //water flux that flow through macro pores can maximally store
            water_def_total += cbm::bound_min( 0.0, (sc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]);
        }

        //duration of timestep (hr-1)
        double const fhr( get_time_step_duration());


        //fraction of surface water that is put into bypass flow (m timestep-1)
        double const bypass_fraction( cbm::bound_max( m_param->BY_PASSF() * fhr, 0.99));
        double const bypass_water( wc_.surface_water * bypass_fraction);
        wc_.surface_water -= bypass_water;
        wc_.accumulated_infiltration += bypass_water;

        if ( cbm::flt_greater_equal( water_def_total, bypass_water))
        {

            for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
            {
                double const water_def( cbm::bound_min( 0.0, (sc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]));
                double const water_add( water_def / water_def_total * bypass_water);

                wc_.wc_sl[sl] += water_add / sc_.h_sl[sl];
            }
        }
        else
        {
            for ( size_t sl = 0;  sl < m_soillayers->soil_layer_cnt();  ++sl)
            {
                double const water_def( cbm::bound_min( 0.0, (sc_.wc_fc_sl[sl] - wc_.wc_sl[sl] - wc_.ice_sl[sl]) * sc_.h_sl[sl]));
                wc_.wc_sl[sl] += water_def / sc_.h_sl[sl];
            }
            /* add surplus to percolation out of last soil layer */
            wc_.accumulated_percolation += bypass_water - water_def_total;
            wc_.accumulated_waterflux_sl[m_soillayers->soil_layer_cnt()-1] += bypass_water - water_def_total;
        }
    }
}

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

{
    double const sks( WaterFlowSaturated( _sl, _fact_impedance));
    double const fsl( get_fsl( _sl));
    double const fhr( get_time_step_duration());

    double slope(( _sl < m_soillayers->soil_layers_in_litter_cnt()) ? m_param->SLOPE_FF() : m_param->SLOPE_MS());

    //max 0.005 means that no sks values up to 9.88 cm min-1 are allowed !!
    double const travelTime( std::max( 0.005, 1.0 - log10( wc_.sks_sl[_sl] * cbm::M_IN_CM * cbm::MIN_IN_DAY * fhr)));
    return( cbm::bound_max( cbm::sqr(1.0 - (sc_.wc_fc_sl[_sl] / ((wc_.wc_sl[_sl] + wc_.ice_sl[_sl]) * slope)))
                            * 0.5 * wc_.wc_sl[_sl] * sc_.h_sl[_sl] * fsl * (1.0 - exp(-1.0 / travelTime)) * _fact_impedance,
                            sks));
}

WaterFlowSaturated()

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

private

Returns saturated water flow per time step accounting for impedance

{
    double const sks( wc_.sks_sl[_sl] * cbm::M_IN_CM * cbm::MIN_IN_DAY * get_time_step_duration() * _fact_impedance);
    return sks;
}

Member Data Documentation

IMAX_W

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

static

number of iteration steps within a day