Case-Study VII: Bayesian inference of HYDRUS-1D using soil moisture data

The seventh case study considers the modeling of the soil moisture regime of an agricultural field near Jülich, Germany. Soil moisture content was measured with Time Domain Reflectometry (TDR) probes at 6 cm deep at 61 locations in a 50 × 50 m plot. The TDR data were analysed using the algorithm described in Heimovaara et al. (1990) and the measured apparent dielectric permittivities were converted to soil moisture contents using the empirical relationship of Topp (1980). Measurements were taken on 29 days between 19 March and 14 October 2009, comprising a measurement campaign of 210 days. For the purpose of the present study, the observed soil moisture data at the 61 locations were averaged to obtain a single time series of water content. Precipitation and other meteorological variables were recorded at a meteorological station located 100 m west of the measurement site. Details of the site, soil properties, experimental design and measurements are given by Scharnagl et al. (2011) and interested readers are referred to this publication for further details.

The HYDRUS-1D model of Simunek et al. (2008) was used to simulate variably saturated water flow in the agricultural field (see Figure 7.01). This model solves Richards' equation for given (measured) initial and boundary conditions

where θ (cm3/cm3) denotes moisture content, t (days) denotes time, z (cm) is the vertical (depth) coordinate, h (cm) signifies the pressure head, and K(h) (cm day-1) is the unsaturated soil hydraulic conductivity.

Figure 7.01. Schematic representation of HYDRUS-1D model setup for a field plot near Jülich, Germany.

To solve Equation (7.01) numerically the soil hydraulic properties need to be defined. We use herein the van Genuchten-Mualem (VGM) model (van Genuchten, 1980)

where θs (cm3/cm3) and θr (cm3/cm3) signify the saturated and residual soil water content, respectively, α (cm-1), n (-) and m = 1 - 1/n (-) are shape parameters, Ks (cm day-1) denotes the saturated hydraulic conductivity, and λ represents a pore-connectivity parameter. The effective saturation, Se (-) is defined as

Observations of daily precipitation and potential evapotranspiration were used to define the upper boundary condition. In the absence of direct measurements, a constant head lower boundary condition was assumed, hbot (cm), whose value is subject to inference with DREAM.

Table 7.01 lists the parameters of the HYDRUS-1D model and their prior ranges which are subject to inference using the 210-day period of observed soil moisture measurements.

† is the normal probability density function with mean a and standard deviation b

‡ is the uniform probability density function with lower bound a and upper bound b

£ These parameters are sampled in the log space

Table 7.01: Parameters of the HYDRUS-1D model and their prior uncertainty ranges

The pedotransfer toolbox of ROSETTA was used with soil textural data from the field site to derive prior estimates of the MVG parameters (Scharnagl et al. 2011). These estimates were used to create an explicit prior distribution for each soil hydraulic parameter. These distributions are defined in the last column of Table 7.01, and make sure that the inference with DREAM honors the values derived separately from ROSETTA. Note that ROSETTA does not provide estimates of the parameter hbot, and this parameter was therefore assumed to have a uniform prior with ranges listed in the bottom right cell of the Table.

A Gaussian likelihood was assumed to characterize in probabilistic terms the distance between the HYDRUS-1D simulated soil moisture data and their observed values. This likelihood function is given by

where  and  are the observed and simulate soil moisture data, respectively, I denotes the number of measurements, and |·| denotes the modulus operator (absolute value). The initial population is drawn from the lower and upper ranges listed in Table 7.01 using Latin hypercube sampling. We use N = 10 chains with DREAM to generate samples from the posterior parameter distribution using standard settings for the algorithmic variables.

Implementation of plugin functions

The complete source code can be found in DREAM SDK - Examples\D3\Drm_Example07\Plugin\Src_Cpp

///////////////////////////////////////////////////////////////////////////////////////////////////

// EvaluatorModel.cpp

///////////////////////////////////////////////////////////////////////////////////////////////////

//-------------------------------------------------------------------------------------------------

 HRESULT STDMETHODCALLTYPE CEvaluatorExample07::GetInputParams

//-------------------------------------------------------------------------------------------------

(

 IDreamInputParams* pIDreamInputParams // Pointer to DREAM input data to be defined

)

///////////////////////////////////////////////////////////////////////////////////////////////////

/*

 DESCRIPTION:

 The purpose of this function is to define default DREAM parameters for this particular project.

 This function is called by the main program to initialize data of each newly created simulation

 case. This data can be later modified in the GUI (except parameters locked by flag

 DreamSolver::disable_edit) or can be overwritten while reading case data from a file.

 In the AppTest console application, there is no GUI and this functions actually defines all

 input parameters that are subsequently used by DREAM Solver.

 RETURN VALUE:

 S_OK    - Operation successful

 S_FALSE - Function not implemented (program will continue)

 E_FAIL  - Operation failed (program will stop)

*/

///////////////////////////////////////////////////////////////////////////////////////////////////

{

 // Always check pIDreamInputParams pointer

 if(pIDreamInputParams == nullptr)

 {

   return E_INVALIDARG;

 }

 // Set default input parameters

 pIDreamInputParams->SetProblemDimension(7, DreamSolver::disable_edit);

 pIDreamInputParams->SetNumberOfMarkovChains(10, DreamSolver::enable_edit);

 pIDreamInputParams->SetNumberOfGenerations(1000, DreamSolver::enable_edit);

 pIDreamInputParams->SetLikelihoodChoice(eLikelihood::LIK_102, DreamSolver::disable_edit);

 pIDreamInputParams->SetInitialDistribution(eInitDistrib::eInitUniform, DreamSolver::enable_edit);

 pIDreamInputParams->SetUsePriorDistribution(1, DreamSolver::enable_edit);

 pIDreamInputParams->SetPriorDistributionType(ePriorType::eUnivariate, DreamSolver::enable_edit);

 pIDreamInputParams->SetOutlierTest(eOutlierTest::eOutlierIqr, DreamSolver::enable_edit);

 pIDreamInputParams->SetBoundHandling(eBoundHandling::eBoundUnbounded, DreamSolver::enable_edit);

 // Example can run in parallel mode (recommended mode)

 pIDreamInputParams->SetParallelMode(eParallelModeYes, DreamSolver::enable_edit);

 return S_OK;

}

//-------------------------------------------------------------------------------------------------

 HRESULT STDMETHODCALLTYPE CEvaluatorExample07::InitData

//-------------------------------------------------------------------------------------------------

(

 IDreamPluginInfo* pIPluginInfo            // Plugin information

, IDreamInputParViewer* pIDreamInputParams  // Current DREAM input parameters

)

///////////////////////////////////////////////////////////////////////////////////////////////////

/*

 DESCRIPTION:

 This function is called to prepare:

 A/ DREAM input data before accessing them via the IDreamDataSource interface and

 B/ all other local data required by evaluator for the calculation.  

 It can be also used to read needed data from a file in directory pIPluginInfo->m_strModelDataDir

 (directory ...Drm_ExampleXX\Model\Data).

         

 Remark 1: This function prepares DREAM input data in the evaluator according to current

 DREAM parameters defined by pIPluginInfo interface. Note that some parameters can be

 different from default parameters defined in function CEvaluatorExample::GetInputParams().

 Remark 2: This function should be called before using functions

 CEvaluatorExample::GetMinMaxData, ::GetNormalData, ::GetPriorData, ::GetPriorDataCustom,

 ::GetMeasurementData and ::GetBayesData. However, function CEvaluatorExample::GetInputParams()

 can be called independently.

 RETURN VALUE:

 S_OK    - Operation successful

 S_FALSE - Function not implemented (program will continue)

 E_FAIL  - Operation failed (program will stop)

*/

///////////////////////////////////////////////////////////////////////////////////////////////////

{

 // Check parameters

 if(pIPluginInfo == nullptr || pIDreamInputParams == nullptr)

 {

   // Invalid argument

   return E_INVALIDARG;

 }

 // Load data only one time

 if(LoadData(pIPluginInfo->m_strModelDataDir))

 {

   return S_OK;

 }

 return E_FAIL;

}

//-------------------------------------------------------------------------------------------------

 HRESULT STDMETHODCALLTYPE CEvaluatorExample07::GetMinMaxData

//-------------------------------------------------------------------------------------------------

(

 IDreamInputParViewer* pInputPar // Interface to DREAM input parameters

, IDreamDataMinMax* pIMinMaxData  // Interface to Min/Max data

)

///////////////////////////////////////////////////////////////////////////////////////////////////

{

 // Always check pIMinMaxData pointer

 if(pIMinMaxData == nullptr)

 {

   return E_INVALIDARG;

 }

 // Fill min max data

 pIMinMaxData->AddItem( 0.0430,  0.0910);

 pIMinMaxData->AddItem( 0.4090,  0.4810);

 pIMinMaxData->AddItem(-2.5528, -2.0706);

 pIMinMaxData->AddItem( 0.1790,  0.2670);

 pIMinMaxData->AddItem(-2.2366, -0.0800);

 pIMinMaxData->AddItem(-5.4900,  6.2700);

 pIMinMaxData->AddItem(-250.0,   -50.0);

 // Data added successfully

 return S_OK;

}

//-------------------------------------------------------------------------------------------------

 HRESULT STDMETHODCALLTYPE CEvaluatorExample07::GetPriorData

//-------------------------------------------------------------------------------------------------

(

 IDreamInputParViewer* pInputPar // Interface to DREAM input parameters

, IDreamDataPrior* pIPriorData    // Interface to prior distribution data

)

///////////////////////////////////////////////////////////////////////////////////////////////////

{

 pIPriorData->AddUnivariateItem(0.067, 0.006, eDistributionUnivariate::eUnivariate_Normal);

 pIPriorData->AddUnivariateItem(0.445, 0.009, eDistributionUnivariate::eUnivariate_Normal);

 pIPriorData->AddUnivariateItem(-2.310, 0.060, eDistributionUnivariate::eUnivariate_Normal);

 pIPriorData->AddUnivariateItem(0.223, 0.011, eDistributionUnivariate::eUnivariate_Normal);

 pIPriorData->AddUnivariateItem(-1.160, 0.270, eDistributionUnivariate::eUnivariate_Normal);

 pIPriorData->AddUnivariateItem(0.390, 1.470, eDistributionUnivariate::eUnivariate_Normal);

 pIPriorData->AddUnivariateItem(-250.0, -50.0, eDistributionUnivariate::eUnivariate_Uniform);

 return S_OK;

}

//-------------------------------------------------------------------------------------------------

 HRESULT STDMETHODCALLTYPE CEvaluatorExample07::PrepareWorkDir

//-------------------------------------------------------------------------------------------------

(

 BSTR strSourceDir

, BSTR strDestinationDir

)

///////////////////////////////////////////////////////////////////////////////////////////////////

{

 // Copy all files to the new directory

 if(!DrxCopyFiles(strSourceDir, strDestinationDir))

 {

   return E_FAIL;

 }

 // Create new level_01.dir with the correct path

 std::wstring LevelDirFile(strDestinationDir);

 LevelDirFile += L"\\level_01.dir";

 std::wfstream ofs(LevelDirFile, std::ofstream::out | std::ofstream::trunc);

 if(!ofs)

 {

   std::wcerr << L"Unable to create " << LevelDirFile << L"\n";

   return E_FAIL;

 }

 

 std::wstring NewPath(strDestinationDir);

 ofs << NewPath;

 ofs.flush();

 ofs.close();

 return S_OK;

}

//-------------------------------------------------------------------------------------------------

 HRESULT STDMETHODCALLTYPE CEvaluatorExample07::EvaluateProposal

//-------------------------------------------------------------------------------------------------

(

 IDreamInputParViewer* pInputPar // Interface to DREAM input parameters

, int iSim              // input parameter (i-th simulation / Markov chain)

, IDreamMatrix* x       // input matrix (m_NumberOfMarkovChains x m_ProblemDimension)

, IDreamMatrix* res     // output matrix (Data_Dimension x m_NumberOfMarkovChains)

, BSTR strWorkingDir    // path to the current working directory (for current calculation thread)

, BSTR strModelDataDir  // path to the permanent directory with model data

, BSTR strModelBinDir   // path to the permanent directory with model executables

)

///////////////////////////////////////////////////////////////////////////////////////////////////

/*

 DESCRIPTION:

 The purpose of this function is to evaluate the current proposal (i.e. one single step

 of IDreamInputParams::GetNumberOfGenerations() steps) for iSim-th Markov chain. It computes

 the likelihood (or log-likelihood or the residual vector to be compared with the measurement

 - see the technical manual) for iSim-th row of input matrix "x" and returns the result

 in output matrix "res".  

 

 Parameter iSim is from interval <0, IDreamInputParams::GetNumberOfMarkovChains() - 1>.

 If pInputPar->GetVectorization() == eVectorization::eVectYes, then iSim = 0 and all proposals

 (for all Markov chains) should be evaluated in one step, i.e. during one call of this function.

 RETURN VALUE:

 S_OK    - Operation successful

 S_FALSE - Function not implemented (program will continue)

 E_FAIL  - Operation failed (program will stop)

*/

///////////////////////////////////////////////////////////////////////////////////////////////////

{

 vector<double> ix(x->_colCount);

 // Copy row data

 for(int iCol = 0; iCol < x->_colCount; iCol++)

 {

   ix(iCol) = x->_data[iSim * x->_colCount + iCol];

 }

 // Back-transform parameters

 matrix<double> xv(1, ix.size(), 0);

 row(xv, 0) = ix;

 xv(0, 2) = pow(10.0, xv(0, 2));

 xv(0, 3) = pow(10.0, xv(0, 3));

 xv(0, 4) = pow(10.0, xv(0, 4));

 // Working directory

 std::wstring path(strWorkingDir);

 // Modify SELECTOR.IN

 if(!CDreamMath<>::ReplaceData(path + L"\\SELECTOR.IN", "thr",

                               project(xv, range(0, 1),range(0, 6))))

 {

   return E_FAIL;

 }

 // Modify PROFILE.DAT

 matrix<double> initial(m_loadData.m_initial);

 column(initial, 2) = vector<double>(initial.size1(), xv(0, 6));

 if(!CDreamMath<>::ReplaceData(path + L"\\PROFILE.DAT", "Mat", initial))

 {

   return E_FAIL;

 }

 // Modify ATMOSPH.IN

 matrix<double> boundary(m_loadData.m_boundary);

 column(boundary, 6) = vector<double>(boundary.size1(), xv(0, 6));

 if(!CDreamMath<>::ReplaceData(path + L"\\ATMOSPH.IN", "tAtm", boundary))

 {

   return E_FAIL;

 }

 // Run HYDRUS-1D

 std::wstring strHydrusCalcExe = strModelBinDir;

 strHydrusCalcExe += L"\\H1D_CALC.EXE";

 if(!CDreamMath<>::RunEXE(strHydrusCalcExe, path))

 {

   return E_FAIL;

 }

 // Read text file containing the time series of simulated soil water contents

 matrix<double> outputData;

 if(!ReadDataHydrus1D(path + L"\\OBS_NODE.OUT", "time ", outputData))

 {

   return E_FAIL;

 }

 vector<double> hoy   = column(outputData, 0);

 vector<double> water = column(outputData, 2);

       // Filter sims for measurement dates

 size_t nData = m_loadData.m_hoy.size();

 vector<size_t> indeces(nData, 0);

 for(size_t i = 0; i < nData; i++)

 {

   size_t iIndex;

   if(!Find<double>(iIndex, hoy, m_loadData.m_hoy(i), [](double a, double b)

   {return abs(a - b) < FLT_EPSILON ? true : false; }))

   {

     res->_data[iSim] = -1.e100; // If HYDRUS-1D did not converge,

                                 // set the log-likelihood to some arbitrary low value

     return S_OK;

   }

   indeces(i) = iIndex;

 }

 // Compute residuals

 vector<double> waterModify;

 CDreamMath<>::ObtainValuesFromIndeces(waterModify, indeces, water);

 vector<double> epsilon = m_loadData.m_water - waterModify;

       

       // Compute log-likelihood

 ElementPow<double>(epsilon, 2);

 res->_data[iSim] = -(double)nData / 2.0 * log(sum(epsilon));

 return S_OK;

}

///////////////////////////////////////////////////////////////////////////////////////////////////

// Local functions

///////////////////////////////////////////////////////////////////////////////////////////////////

//-------------------------------------------------------------------------------------------------

 bool CEvaluatorExample07::ReadDataHydrus1D

//-------------------------------------------------------------------------------------------------

(

 const std::wstring& file

, const std::string& keyWord

, matrix<double>& loadData

)

///////////////////////////////////////////////////////////////////////////////////////////////////

{

 std::ifstream ifs(file, std::ifstream::in);

 if(!ifs)

 {

   std::wcerr << "Unable to open " << file <<"\n";

   return false;

 }

 typedef boost::tokenizer<boost::char_separator<char>> Tokenizer;

 boost::char_separator<char> sep(" ");

 std::string data;

 size_t iLine = 0;

 size_t startLine = 0;

 size_t endLine = 0;

 size_t nTok = 0;

 while(!ifs.eof())

 {

   getline(ifs, data);

   Tokenizer tok(data, sep);

   if(data.find(keyWord) != std::string::npos)

   {

     nTok = std::distance(tok.begin(), tok.end());

     startLine = iLine + 1;

   }

   else if(data.find("end") != std::string::npos)

   {

     endLine = iLine;

   }

   iLine++;

 }

 size_t nRow = endLine - startLine;

 loadData = matrix<double>(nRow, nTok);

 ifs.clear();

 ifs.seekg(0, std::ios::beg);

 iLine = 0;

 size_t iRow = 0;

 while(!ifs.eof())

 {

   getline(ifs, data);

   if(startLine <= iLine && iLine < endLine)

   {

     Tokenizer tok(data, sep);

     vector<double> values(nTok);

     size_t i = 0;

     for(Tokenizer::iterator iter = tok.begin(); iter != tok.end(); ++iter)

     {

       try

       {

         if(*iter == "T")

         {

           values(i) = 1;

         }

         else if(*iter == "F")

         {

           values(i) = 0;

         }

         else

         {

           values(i) = boost::lexical_cast<double>(*iter);

         }

       }

       catch(boost::bad_lexical_cast&)

       {

         std::wcerr << "warning: value on Line" << iLine << ", Column " << i << "from "

                    << file << " could not be converted to double = 0." << std::endl;

         values(i) = 0.;

       }

       i++;

     }

     row(loadData, iRow) = values;

     iRow++;

   }

   iLine++;

 }

 ifs.close();

 return true;

}

//-------------------------------------------------------------------------------------------------

 bool CEvaluatorExample07::LoadData(std::wstring path)

//-------------------------------------------------------------------------------------------------

{

 // Load structure with observation data

 if(!CDreamMath<bool>::ReadData(path + L"\\Meas_WaterInd.txt", m_loadData.m_waterInd))

   return false;

 if(!CDreamMath<>::ReadData(path + L"\\Meas_Water.txt", m_loadData.m_water))

   return false;

 vector<double> water;

 CDreamMath<double>::ObtainValuesFromBool(water, m_loadData.m_waterInd, m_loadData.m_water);

 m_loadData.m_water = water;

 if(!CDreamMath<>::ReadData(path + L"\\Meas_Hoy.txt", m_loadData.m_hoy))

   return false;

 vector<double> hoy;

 CDreamMath<double>::ObtainValuesFromBool(hoy, m_loadData.m_waterInd, m_loadData.m_hoy);

 m_loadData.m_hoy = hoy;

 m_loadData.m_obsNode = 1;

 // Provide data needed to modify initial condition

 if(!CDreamMath<>::ReadData(path + L"\\ProfileDat_InitialConditions.txt", m_loadData.m_initial))

   return false;

 // Provide data needed to modify the lower boundary condition

 if(!CDreamMath<>::ReadData(path + L"\\AtmosphIn_BoundConditions.txt", m_loadData.m_boundary))

   return false;

 return true;

}