#include "DrainageToSaturationZone.h"

//Soil water is represented as three components, with average soil moisture deficit representing total storage between soil surface and groundwater table: 
//Note: 1) evapotranspiration zone or StorageDeficit_FieldCapacity_SoilEvapZone_Max_m with capillarity held water for ET, releasing gravitational water into unsaturated zone; it can saturate from below
//Note: 2) gravitational zone or AveSMD representing gravitational water, draining into saturated zone; it can saturate from below into evapotranspiration zone
//Note: 3) saturated zone for subsurface runoff; when it extends into evapotranspiration zone it fills the gravitational voids
//Note: The evapotranspiration and unsaturated zones are conceptually situated alongside each other, each filling a different volumetric component
//Note: The evapotranspiration zone storage depth extends to the limit of evapotranspiration accessing capillarity held water, a constant for all TI bins
//Note: The unsaturated zone extends from the groundwater table to the ground, and accounts for all gravitational held water, including evapotranspiration zone
//Note: Storage_VadoseZone_TI_m (m) is unsaturated zone storage, initialized to 0; represents only gravitational water
//Note: Drainage_VadoseZone_TI_m (m) is unsaturated zone runoff, initialized to 0; represents release of gravitational water

//Note: Catchment average soil moisture deficit (StorageDeficit_VadoseZone_m, AveSMD) determines water storage in each Topographic Index bin
//Note: Saturation excess runoff is an outcome of AveSMD and local SMD values, based on variable source area and topographic index (TI) theory
//Note: The catchment Bulk Area uses local SMD values to determine saturation excess in each TI bin, StorageDeficit_VadoseZone_TI_m
//Note: The Green Infrastructure (GI) Area uses AveSMD via Depth_to_Groundwater_Table_m = StorageDeficit_VadoseZone_m / Porosity_Drainable_m3pm3
//Note: Green Infrastructure (GI) uses Depth_to_Groundwater_Table_m in Groundwater_StormwaterDevice.cpp to determine saturation by groundwater


void DrainageToSaturationZone::calculate(Inputs *input, CompactRagged* beC, int TI_bin, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//ia = TI_bin as current TI bin to get appropriate topographic index properties
	int ia = TI_bin;

	//Storage routine to determine saturation excess runoff and drainage
	//Runoff_SatExcess_TI_m (m) is reset to 0 at each time step
	beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_TI_m")[ia] = 0.0;
	//Discharge_Groundwater_TI_m (m) is reset to 0 at each time step
	beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Flux_SatExGW_to_SoilFieldCapacity_TI_m")[ia] = 0.0;

	//Storage_VadoseZone_TI_m (m) is increased by Drainage_SoilEvapZone_TI_m (m) flowing out of ET zone, including Drainage_macroPore_m (m), created in DrainageToEvapotranspirationZone.cpp
	// TOPMODEL code uses: IF(SRZ(IA).LT.0.)THEN SUZ(IA) = SUZ(IA) - SRZ(IA) and SRZ(IA) = 0; note SRZ(IA) is a deficit, negative deficits are flow into SUZ 
	beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_TI_m")[ia] + beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_TI_m")[ia];

	//if Storage_VadoseZone_TI_m (m) > StorageDeficit_VadoseZone_TI_m (m), saturation excess runoff occurs; note Storage_VadoseZone_TI_m is a storage and StorageDeficit_VadoseZone_TI_m is a deficit (or capacity)
	// TOPMODEL code uses: IF(SUZ(IA).GT.SD(IA)) THEN EX(IA) = SUZ(IA) - SD(IA) and SUZ(IA)=SD(IA)
	if (Inputs::isGreaterThan(beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_TI_m")[ia], beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia])) {

		//Runoff_SatExcess_TI_m (m) equals Storage_VadoseZone_TI_m (m) - StorageDeficit_VadoseZone_TI_m (m); note Storage_VadoseZone_TI_m is a storage and StorageDeficit_VadoseZone_TI_m is a deficit (or capacity)
		//Note: Consider Refactor to check if Runoff_SatExcess_TI_m * TI_binArea_frac fits within available pervious depression storage
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_TI_m")[ia] - beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia];

		//Storage_VadoseZone_TI_m (m) reset to equal the StorageDeficit_VadoseZone_TI_m (m) deficit (or capacity)
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia];

		//Flux_SatExGW_to_SoilFieldCapacity_TI_m (m) is StorageDeficit_FieldCapacity_SoilEvapZone_TI_m (m)
		//Note: This is the amount to fill soil voids to field capacity in the evapotranspiration zone given it is saturated
		//Note: The Flux_SatExGW_to_SoilFieldCapacity_TI_m may have been activated in DrainageToEvapotranspirationZone::calculate
		//Note: Evapotranspiration zone manages the moisture between Soil_FieldCapacity_m3pm3 and Soil_WiltingPoint_m3pm3, which is filled from below when Storage_VadoseZone_TI_m > StorageDeficit_VadoseZone_TI_m
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Flux_SatExGW_to_SoilFieldCapacity_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia];
		//StorageDeficit_FieldCapacity_SoilEvapZone_TI_m (m) is a deficit, aka capacity, which is reduced to zero saturation by groundwater
		//Note: Setting StorageDeficit_FieldCapacity_SoilEvapZone_TI_m (m) to zero will enable maximum rates for EvapoTranspiration_SoilEvapZone_TI_m (m)
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia] = 0.0;

		//Drainage_SoilEvapZone_MacroPore_frac (fraction) and Drainage_SoilEvapZone_Infil_frac (fraction) defined as fractions that sum to one, based on relative contribution to Drainage_SoilEvapZone_TI_m (m)
		//Note: Drainage_macroPore_PreRunoffSatExcess_m (m) is a constant used in quotient, not Drainage_macroPore_m (m) which may be reduced during TI loop
		if (Inputs::isLessThanOrAlmostEqual(beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_TI_m")[ia], 0.0)) {
			beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_MacroPore_frac") = 0;
		}
		//Else Drainage_SoilEvapZone_TI_m > 0 then enter into division by Drainage_SoilEvapZone_TI_m
		else {
			//Drainage_SoilEvapZone_MacroPore_frac (fraction) is ratio of Drainage_macroPore_PreRunoffSatExcess_m and Drainage_SoilEvapZone_TI_m, updated for each TI bin
			beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_MacroPore_frac") = Inputs::SafeDivide(beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_macroPore_PreRunoffSatExcess_m"), beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_TI_m")[ia]);
		}
		//Drainage_SoilEvapZone_Infil_frac (fraction) is 1 minus Drainage_SoilEvapZone_MacroPore_frac, so the two terms sum to one
		beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_Infil_frac") = 1 - beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_MacroPore_frac");

		//Remove Runoff_SatExcess_TI_m[ia] within TI loop from Drainage_SoilEvapZone_m (m), Infiltration_viaInfilEx_m (m), Infiltration_viaSatEx_m (m), and Drainage_macroPore_m (m)
		//Fraction of Runoff_SatExcess_TI_m removed is proportional to Drainage_SoilEvapZone_MacroPore_frac, Drainage_SoilEvapZone_Infil_frac, InfiltExcessGovernedArea_frac, and TI_binArea_frac
		//Note: Not updating: a) Storage_VadoseZone_TI_m[ia] (m) which was updated when it is reset to StorageDeficit_VadoseZone_TI_m; 
		//Note: Not updating: b) StorageDeficit_FieldCapacity_SoilEvapZone_TI_m[ia] (m) which was updated in DrainageToEvapotranspirationZone.cpp and SoilEvapZoneETCalc.cpp
		//Note: Not updating: c) Infiltration_m (m), sum of Infiltration_viaSatEx_m and Infiltration_viaInfilEx_m, which will be updated at end of TI for loop
		//Drainage_SoilEvapZone_TI_m (m) is reduced by Runoff_SatExcess_TI_m (m)
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_TI_m")[ia] - beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_TI_m")[ia];
		//Infiltration_TI_m (m) is reduced by Runoff_SatExcess_TI_m (m)
		//Note: Runoff_SatExcess_TI_m is most efficiently removed from Infiltration_TI_m rather than from Infiltration_viaInfilEx_m and Infiltration_viaSatEx_m
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Infiltration_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Infiltration_TI_m")[ia] - beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_TI_m")[ia] * beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_Infil_frac");
		//Drainage_macropore_TI_m (m) is reduced by Runoff_SatExcess_TI_m (m) scaled to Drainage_SoilEvapZone_Infil_frac
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_macropore_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_macropore_TI_m")[ia] - beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_TI_m")[ia] * beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_MacroPore_frac");
	}

	//else if Evapotranspiration_Depth_m > (StorageDeficit_VadoseZone_TI_m / Porosity_Drainable_m3pm3) then Evapotranspiration satisfied by groundwater
	//Note: Recharging StorageDeficit_FieldCapacity_SoilEvapZone_TI_m with groundwater 1-time step behind RootZoneEvapoTranspirationCalc::Calculate
	else if (Inputs::isGreaterThan(beC->by_key(MapPixel_ID, DataFolder_ID, "Evapotranspiration_Depth_m"), (beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia] / input->InputXml["Porosity_Drainable_m3pm3"]))) {

		//Evapotranspiration_Depth_Saturated_frac (frac) = (Evapotranspiration_Depth_m - (StorageDeficit_VadoseZone_TI_m / Porosity_Drainable_m3pm3)) / Evapotranspiration_Depth_m
		//Note: Computing the fraction of evapotranspiration depth soil profile that is saturated
		double Evapotranspiration_Depth_Saturated_frac = Inputs::SafeDivide((beC->by_key(MapPixel_ID, DataFolder_ID, "Evapotranspiration_Depth_m") - (beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia] / input->InputXml["Porosity_Drainable_m3pm3"])), beC->by_key(MapPixel_ID, DataFolder_ID, "Evapotranspiration_Depth_m"));

		//Evapotranspiration_Depth_UnSaturated_frac (frac) = 1 - Evapotranspiration_Depth_Saturated_frac
		//Note: Computing the fraction of evapotranspiration depth soil profile that is unsaturated
		double Evapotranspiration_Depth_UnSaturated_frac = 1 - Evapotranspiration_Depth_Saturated_frac;

		//Flux_SatExGW_to_SoilFieldCapacity_TI_m (m) = StorageDeficit_FieldCapacity_SoilEvapZone_TI_m * Evapotranspiration_Depth_Saturated_frac
		//Note: Evapotranspiration zone manages the moisture between Soil_FieldCapacity_m3pm3 and Soil_WiltingPoint_m3pm3, which is filled from below when Evapotranspiration_Depth_m > (StorageDeficit_VadoseZone_TI_m / Porosity_Drainable_m3pm3) to depth of Evapotranspiration_Depth_Saturated_frac
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Flux_SatExGW_to_SoilFieldCapacity_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia] * Evapotranspiration_Depth_Saturated_frac;
		//StorageDeficit_FieldCapacity_SoilEvapZone_TI_m (m) = 0.0 * Evapotranspiration_Depth_Saturated_frac + StorageDeficit_FieldCapacity_SoilEvapZone_TI_m * Evapotranspiration_Depth_UnSaturated_frac
		//Note: StorageDeficit_FieldCapacity_SoilEvapZone_TI_m is a deficit, aka capacity, which is reduced by groundwater
		//Note: Setting StorageDeficit_FieldCapacity_SoilEvapZone_TI_m (m) to zero will enable maximum rates for EvapoTranspiration_SoilEvapZone_TI_m (m)
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia] = 0.0 * Evapotranspiration_Depth_Saturated_frac + beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia] * Evapotranspiration_Depth_UnSaturated_frac;
	}

	//Updating scalar values of infiltration and drainage fluxes within TI loop by adding TI bin vectors, updated above, now scaled to the TI bin area
	//Infiltration_viaInfilEx_m (m) is incrementally increased by product of Infiltration_TI_m and InfiltExcessGovernedArea_frac, scaled to TI_binArea_frac, during the TI bin loop
	beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_m") + beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Infiltration_TI_m")[ia] * beC->by_key(MapPixel_ID, DataFolder_ID, "InfiltExcessGovernedArea_frac") * beC->by_key(MapPixel_ID, DataFolder_ID, "TI_binArea_frac");
	//Infiltration_viaSatEx_m (m) is incrementally increased by product of Infiltration_TI_m and (1-InfiltExcessGovernedArea_frac), scaled to TI_binArea_frac, during the TI bin loop
	beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaSatEx_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaSatEx_m") + beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Infiltration_TI_m")[ia] * (1 - beC->by_key(MapPixel_ID, DataFolder_ID, "InfiltExcessGovernedArea_frac")) * beC->by_key(MapPixel_ID, DataFolder_ID, "TI_binArea_frac");
	//Drainage_macroPore_m (m) is incrementally increased by Drainage_macropore_TI_m * TI_binArea_frac during the TI bin loop
	beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_macroPore_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_macroPore_m") + beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_macropore_TI_m")[ia] * beC->by_key(MapPixel_ID, DataFolder_ID, "TI_binArea_frac");
	//Drainage_SoilEvapZone_m (m) is incrementally increased by Drainage_SoilEvapZone_TI_m * TI_binArea_frac during the TI bin loop
	beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_m") + beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_TI_m")[ia] * beC->by_key(MapPixel_ID, DataFolder_ID, "TI_binArea_frac");
	//Flux_SatExGW_to_SoilFieldCapacity_m (m) is incrementally increased by Flux_SatExGW_to_SoilFieldCapacity_TI_m * TI_binArea_frac during the TI bin loop
	beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_SatExGW_to_SoilFieldCapacity_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_SatExGW_to_SoilFieldCapacity_m") + 
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Flux_SatExGW_to_SoilFieldCapacity_TI_m")[ia] * beC->by_key(MapPixel_ID, DataFolder_ID, "TI_binArea_frac");

	//Drainage routine, treating all water as gravitational other than water in soil evapotranspiration zone water held by capillarity 
	//Drainage_VadoseZone_TI_m (m) reset to zero at start of each TI value
	beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_VadoseZone_TI_m")[ia] = 0.0;
	
	//if StorageDeficit_VadoseZone_TI_m (m), the local soil moisture deficit, > 0 then there is an unsaturated zone, and it can drain to lower level groundwater
	if (Inputs::isGreaterThan(beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia], 0.0 )) {
		//Drainage_VadoseZone_TI_m (m), unsaturated zone flow to groundwater = (unsaturated zone storage / (local soil moisture deficit * Unsaturated zone time delay (Td_dt, adjusted to time step) constant)
		//Eq 16 of Beven and Wood (1983) for time delay of water release from unsaturated zone, for water above field capacity
		//Beven and Wood (1983). Catchment Geomorphology and the Dynamics of Runoff Contributing Areas, Journal of Hydrology. 65:139-158. In Special Issue: Scale Problems in Hydrology, Editors: Rodriguez-Iturbe and V.K. Gupta.

		//if Td_h (h) greater than 0 then there is a delay in the drainage of Storage_VadoseZone_TI_m (m) and generation of Drainage_VadoseZone_TI_m (m), adjusted for SimulationTimeStep_Duration_sec
		if (Inputs::isGreaterThan(beC->by_key(MapPixel_ID, DataFolder_ID, "VadoseZone_Drainage_Delay_h"), 0.0)) {
			beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_VadoseZone_TI_m")[ia] = Inputs::SafeDivide(beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_TI_m")[ia], (beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia] * beC->by_key(MapPixel_ID, DataFolder_ID, "VadoseZone_Drainage_Delay_h") * Ratio_s_to_hr / input->SimulationTimeStep_Duration_sec[timeStep]));
		}
		//if Td_h (h) <= 0 there is no delay in the drainage of Storage_VadoseZone_TI_m (m) and generation of Drainage_VadoseZone_TI_m (m)
		else {
			beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_VadoseZone_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_TI_m")[ia];
		}

		//if (unsaturated zone flow > unsaturated zone storage), set it equal to the unsaturated zone storage
		if (beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_VadoseZone_TI_m")[ia] > beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_TI_m")[ia]) {
			beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_VadoseZone_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_TI_m")[ia];
		}		
		//Storage_VadoseZone_TI_m (m) storage decreased by flow out of the unsaturated zone and into groundwater 
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_TI_m")[ia] - beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_VadoseZone_TI_m")[ia];		
	}

	//Control for negative values with if conditional
	//Note: Each of these variables had subtracted Runoff_SatExcess_TI_m, which may have led to negative values
	if (Inputs::isLessThan(beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_m"), 0)) { beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_m") = 0; }
	if (Inputs::isLessThan(beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaSatEx_m"), 0)) { beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaSatEx_m") = 0; }
	if (Inputs::isLessThan(beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_macroPore_m"), 0)) { beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_macroPore_m") = 0; }
	if (Inputs::isLessThan(beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_m"), 0)) { beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_m") = 0; }

	//Drainage_VadoseZone_m (m) = sum Drainage_VadoseZone_TI_m for each TI bin for time step, weighted by watershed fraction; pervious cover area depth
	//Note: HydroPlus weights discrete TI parameter values, i.e, SoilParameter[ia], with TI_binArea_frac, catchment area bounding discrete value
	beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_VadoseZone_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_VadoseZone_m") + beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_VadoseZone_TI_m")[ia] * beC->by_key(MapPixel_ID, DataFolder_ID, "TI_binArea_frac");

	//Storage_VadoseZone_m = sum Storage_VadoseZone_TI_m for each TI bin for time step, weighted by watershed fraction
	//Note: HydroPlus weights discrete TI parameter values, i.e, SoilParameter[ia], with TI_binArea_frac, catchment area bounding discrete value
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_m") + beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_TI_m")[ia] * beC->by_key(MapPixel_ID, DataFolder_ID, "TI_binArea_frac");
}