The impact of sensitivity and uncertainty of soil physical parameters on the terms of the water balance_ Some case studies with default R packages. Part I_ Theory, methods and case descriptions

These papers (part I and part II) emphasize the need for sensitivity and uncertainty analyses. A number of techniques are applied, e.g. latin hypercube sampling, impact response surfaces and Sobol-analyses. Five examples are presented, four of them concerning the numerical model SWAP. The data generation and analysis is performed with standard R packages. Although the computations can be made on any computer, the most timeconsuming examples in this paper have been run on a High Performance Computer Cluster. With the relatively simple Impact Response Surface technique it is shown that variation of the saturated hydraulic conductivity has far less impact than changing the moisture content at saturation. Analyses according to the Sobol-Jansen method show that when the soil physical relationships are described according to Damiano, then the parameter b has a very large influence on the results. If the well-known Mualem Van Genuchten equations are applied, most variation can be explained by the parameter n.


Introduction
Climate change and a growing demand of water for agriculture and urban development cause an increasing pressure on fresh water resources. Therefore, it is important to use these scarce resources as effective as possible. To increase knowledge in agronomy and hydrology, models are used to analyse experimental data, whereas experiments are needed to parameterize models (De Jong Van Lier et al., 2015). During the past decades, these models have become more and more detailed. Though data processing capacities have grown exponentially lately, running a model is time-consuming. Sensitivity and uncertainty analyses should be part of the model development process, but are often neglected.
Efficient sensitivity analysis, particularly for a global sensitivity analysis (GSA) to identify the most important or sensitive parameters, is crucial for understanding complex hydrological models. The sensitivity of model outcomes to input parameters is a key issue in this context. Efficient parameter identification is an important issue for mechanistic agro-hydrological models with a complex and nonlinear property. According to Pianosi et al. (2015), GSA is a term describing a set of mathematical techniques to investigate how the variation in the output of a numerical model can be attributed to variations of its inputs. GSA can be applied for multiple purposes, including: (i) to apportion output uncertainty to the different sources of uncertainty of the model, e.g. unknown parameters, measurement errors in input forcing data, etc. and thus prioritize the efforts for uncertainty reduction; (ii) to investigate the relative influence of model parameters over the predictive accuracy and thus support model calibration, verification and simplification; to understand the dominant controls of a system (model) and (iii) to support model-based decision-making. Sarrazin et al. (2016) state that sensitivity analysis aims to characterize the impact that changes in the model input factors (e.g. parameters, initial states, input data, time/spatial resolution grid etc.) have on the model output (e.g. a statistic of the simulated time series, such as the average simulated streamflow, or an objective function, like the Root Mean Squared Error). Over-parameterization is a well-known and often described problem in numerical models, especially for distributed models. Therefore, methods to reduce the number of parameters via sensitivity analysis are important for the efficient use of these models (Van Griensven et al., 2006). Excellent overviews of sensitivity and uncertainty analysis methods have been presented by Norton (2015) and by Pianosi et al. (2016).
We choose to work with the numerical model SWAP (Soil-Water-Atmosphere-Plant) (Kroes et al., 2017), which simulates the soil moisture flow in the one-dimensional vadose zone and is capable of simulating detailed crop growth. The original version of this model was called SWATR (Feddes et al., 1978). The most-cited publication about the successor of this model, called SWATRE, was written by Belmans et al. (1983). Developments continued, and recently SWAP version 4.0.1 (Kroes et al., 2017), was released (swap.wur.nl), which includes the WOFOST 7.17 model for crop growth (Boogaard et al., 2014). Supit (1994) was the first to publish a detailed system description of the WOFOST model. Recent versions are published on the internet (wofost.wur.nl) and an overview of 25 years WOFOST modelling is given by De Wit et al. (2018).
A global sensitivity analysis of an earlier version of SWAP has been presented by Wesseling et al. (1997). These authors concluded that knowledge of the boundary conditions was more important than knowledge of the soil parameters. Shafiei et al. (2014) present parameter uncertainty of the SWAP model and its effect on model prediction within the generalized likelihood uncertainty estimation (GLUE) framework for two irrigated agricultural fields in a dry region of Iran. Parameter uncertainty analysis of soil hydraulic parameters showed that in spite of similarity of soil texture in both the considered fields, the estimated parameters (i.e. posterior distribution) exhibit different behaviors. This was caused by the dynamics of soil structure which varies considerably within cultivated fields during the growing season. Their study reveals the importance of uncertainty analysis to estimate the degree of reliability associated with model predictions as an important first step for providing decision makers with realistic information about the models outputs.
Although the developments of multi-dimensional models are promising, the data and computation demand of multi-dimensional models of the biosphere is huge. Therefore these models are still mainly used for detailed research applications and are less suitable for regional operational practice.
The objective of the research presented in this paper was threefold: (i) show the influence of the soil physical parameters on the output of unsaturated zone models; (ii) investigate the usability of standard R packages to perform sensitivity and uncertainty analyses on models and (iii) show the advantages of a High Performance Computer Cluster. We will show that a sensitivity or uncertainty analysis is rather simple if the correct hardware and software tools are used (HPC, R, Julia). First we'll briefly present the theory of soil moisture flow. Then the applied sensitivity and uncertainty methods are described, together with the applied tools. These are applied to 5 different cases from earlier projects. The output from a simple program that computes soil moisture storage in a soil is used in the first example. The next four examples are based upon earlier research projects in Argentina and Brazil where the SWAP model was applied. This paper describes the theory behind the simulations, the applied methods and it presents a description of the cases. In the accompanying paper (part II) the results are presented and discussed.

Moisture flow
The 1-dimensional (vertical) flow of water in the soil can be described by Darcy's law: The mass conservation law states: (2) where = moisture content at depth z [L 3 .L −3 ] (cm 3 cm −3 ). This means the change of moisture content with time is equal to the change of flux density with position. Introduction of the differential moisture capacity C [L −1 ](cm −1 ), which is defined as the derivative of moisture content with respect to the prevailing pressure head, or and combining the first two equations yields the partial differential equation for the description of transient (vertical) soil moisture flow: which has been first published by Richards (1931). The K h ( )-and h ( ) relationships are called the hydraulic conductivity function and the soil moisture retention function. The values and shapes of these functions vary with soil type and depth. They may be measured in the laboratory or in the field or derived by inverse modeling. See e.g. Wesseling et al. (2008) and Wesseling (2009)

Description of soil physical characteristics
A lot of attempts have been made over the year to describe the soil physical characteristics (the hydraulic conductivity function and soil moisture retention function) in such a way they can be applied in numerical studies. These methods vary from purely mathematical approximations to physically based equations. In the present paper we apply two different methods: the widely-used Mualem-Van Genuchten method and the one that was recently published by Damiano (2008).
One of the most widely spread methods to describe soil physical relationships has been introduced by Van Genuchten (1980). The K h ( )and h ( )-relationships are described by S-shaped curves that are defined completely with only 6 parameters: where Damiano (2008) presents the relationship between pressure head (= potential) and the moisture content in two parts, depending on the position of related to the value at the inflection point i .
In case The hydraulic conductivity K is computed as

Steady state pressure head profiles
Nowadays most soil water simulation models are transient models. Nevertheless, when characterizing and classifying soils, usually properties of these soils are considered that are obtained under steady state conditions. Examples of these properties are storage capacity and critical distance (see e.g. Wösten et al. (2012) or Wösten et al. (2013)). Therefore our first example will be an application of the steady-state moisture flow. Steady state usually implies that the situation of the soil has not changed for a long time. This can occur after a long period of drought or in case of long-lasting precipitation. In that case we can rewrite Eq. 1 to find a relation between the change in position and the change in pressure head: It is known that the pressure head is zero at the freatic surface. Assuming z = 0 at that position and defining z as positive upward, the relationship between pressure head h and vertical position z (the pressure head profile) can be computed by integrating the previous If the vertical flux value is positive, so directed upward (capillary rise), this integral can be approximated numerically by When computations are performed for a heterogeneous soil profile, the soil physical characteristics of the different soil layers should be considered.
For negative values of q (infiltration), instability may occur, especially when the value of K(h) is close to the (absolute value of) the flux q. Then the ratio is approximating −1 and the denominator in Eq. 14 will approach 0, causing dz to be . For that reason Eq. 14 is, in case of infiltration, rewritten as: Now the pressure head profile can be described by This can be computed numerically as In case of a multi-layer profile, care should be taken that K(h) depends on the position z as well.
The equations presented above are applied peviously for the computation of steady-state pressure head profiles (Wesseling and Brandyk, 1984;Wesseling, 1991). These pressure head profiles can be presented graphically as in Fig. 2a. From the known soil moisture retention curve (Fig. 1) the corresponding moisture contents can easily be found. These moisture contents are shown in Fig. 2b with a blue color. The solid fraction of the soil is presented in brown, clearly showing it is a 3-layer profile.
One of the properties used to classify a soil sample of profile is the storage capacity (also called saturation deficit) at 1 or at 2 mm d · 1 (Wösten et al., 2012). This value can easily be computed from the moisture content profile, as it is represented by the size of the noncolored area in Fig. 2b. Or, written mathematically:

Transient soil moisture flow
If more detailed data (in time and depth) is required, then Eq. 4 has to be solved numerically. Due to the increasing capability of digital computers, the number of numerical simulation models has increased strongly during the past decades. Most of these models simulate both saturated and unsaturated moisture flow in one or more dimensions and solve the Richards equation either with the finite difference method, the finite element method or the boundary element method. Many additional features are available in some of these models, e.g. heat transport, solute transport, frozen soils, simple and detailed crop growth models, etc. One of the most widely-used models is SWAP, which is described by Kroes et al. (2009) and Kroes et al. (2017). See Fig. 3 for a schematic overview of its features.

Uncertainty and sensitivity
It is generally known that uncertainty and sensitivity analysis should be an integral part of the modeling process (Saltelli et al., 2000). Model developers, however, often ignore this need. Main causes are a lack of time, and fear for the high amount of computer time required, because in general a lot of model runs have to be made.
Several definitions of 'uncertainty' and 'sensitivity' can be found in literature (Saltelli, 2017;Baroni and Tarantola, 2014). Here we'll use the definitions found in EPA (2009): • Uncertainty Analysis -Investigates the effects of lack of knowledge or potential errors of the model (e.g. the uncertainty associated with parameter values or model design and output). • Parameter uncertainty is discussed separately in some conventions.
This type of uncertainty is assigned to the data used to calibrate parameter values.
• Model uncertainty: simplification of real-world processes, misspecification of the model structure, use of inappropriate variable or parameter values, aggregation errors, application/scenario.
In general, the term 'data uncertainty' is used to refer to the uncertainty caused by measurement errors, analytical imprecision and limited sample sizes during data collection and treatment (EPA, 2009). In contrast to data uncertainty, variability results from the inherent randomness of certain parameters or measured data, which in turn results from the heterogeneity and diversity in environmental processes. Variability can be better characterized, but is hard to reduce, with further study. Separating variability and uncertainty is necessary to provide greater accountability and transparency. However, variability and uncertainty are inextricably intertwined and ever present in regulatory decision making (EPA, 2009).
According to Xu et al. (2016), parameter sensitivity analysis (SA) is a prerequisite step in the model-building process. The SA method identifies parameters that do or do not have a significant impact on model simulation of real world observations and is critical for reducing the number of parameters required in model validation. Generally, SA can be divided into two different schools: the local SA school and the global one. In the first approach, the local response of model output is obtained by varying the parameters one at a time while holding the others fixed to certain nominal values. This approach has been adopted by some studies because of its easy application. Yet, local SA methods have the known limitations of linearity and normality assumptions and local variations. For complex non-linear models, only global sensitivity analysis (GSA) methods are able to provide relevant information on the sensitivity of model outputs to the whole range of model parameters. An example of a GSA is presented by Wesseling et al. (1997). Huang et al. (2018) describe an integrated system for the dynamic prediction and assessment of agricultural yield using the Sunway TaihuLight supercomputer platform. This system enables parallelization and acceleration for the existing AquaCrop, DNDC (DeNitrification and DeComposition) and SWAP (Soil Water Atmosphere Plant) models, thus facilitating multi-model ensemble and parameter optimization and subsequent drought risk analysis in multiple regions and at multiple scales. By means of testing with varying core group numbers these authors show that computation time can be reduced by between 2.6 and 3.6 times. Based on the powerful computing capacity, a county-level modelparameter optimization (2043 counties for the period from 1996 to 2007) by Bayesian inference and multi-model ensemble using BMA (Bayesian Model Average) methods were performed, demonstrating the enhancements in predictive accuracy that can be achieved.
Model predictions with respect to yield are highly sensitive to soil hydraulic parameterization, even when this sensitivity does not show up strongly in predicted soil water content or pressurehead values (De Jong Van Lier et al., 2015).
Over-parameterisation is a well-known and often described problem with hydrological models, especially for distributed models. Therefore, sensitivity analysis methods that aim to reduce the number of  J. Wesseling, et al. Computers and Electronics in Agriculture xxx (xxxx) xxxx parameters that require fitting with input and output data are common. These methods identify parameters that do or do not have a significant influence on model simulations of real world observations for specific catchments (Van Griensven et al., 2006).

Materials and methods
Several techniques are applied during this research: Impact Response Surfaces (IRS), Sobol sensitivity analysis and Latin Hypercube Sampling (LHS). These methods will be described briefly here. As was mentioned in the introduction, the research presented in this paper had multiple goals. These goals will be reached by performing computations in five different case studies. These cases will be described here, the results of the simulations will be presented in the next section. A summary of the applied methods, input variables and output variables for these cases is presented in Table 1.

Impact response surfaces (IRS)
One of the methods widely applied in meteorological sciences and crop sciences is the Impact Response Surface. This is a two-dimensional chart showing the impact of the (combined) change of two parameters on some output value of a model (Fronzek et al., 2018).

Sobol sensitivity analysis
A frequently used method for sensitivity analysis is the variancebased one presented in Sobol (1993) and Sobol (2001). The total variance of the model output, produced by parameter variation, is decomposed into partial variances. These partial variances can be related to the model parameters. The parameter sensitivity is quantified by first-and total-order indices (S 1 and S Ti resp.). These are computed as where V stands for variance, E for expectation, Y for model output, X i is the i th input parameter and X~i the matrix of all model inputs except the i th input parameter. The first order index measures the contribution of an individual parameter x i to the total variance of the model outcome (main effect). The total-order index involves the main effect of x i and S interactions with other parameters. The difference between the indices gives the interaction effect that parameter i has with other parameters. High sensitivity is associated with high values of S i and S Ti . Additive models without existence of parameter interaction have the property that S i and S Ti are equal and the sum of all values S i and all values S Ti is equal to 1. If S Ti is greater than S i and the sum of all S i is less than 1, then the model is non-additive (Saltelli et al., 2008;Stahn et al., 2017). The required number of function evaluations N e can be computed from where N d is the number of draws from the random generator and N p is the number of parameters to be investigated. This method has been applied by several authors, see e.g. Stahn et al. (2017) and Wainwright et al. (2014). Mulder et al. (2016) applied this method in combination with the SWAP model to obtain the sensitivity of the parameters for evaporation. Baroni and Tarantola (2014) developed a framework based on this theory. They tested it with the SWAP model and concluded that the boundary conditions of the system had a large influence on the results of the simulations. This is in agreement with the conclusions of Wesseling et al. (1997). Pianosi et al. J. Wesseling, et al. Computers and Electronics in Agriculture xxx (xxxx) xxxx (2015) developed a similar toolbox in Matlab. This toolbox presents a generally applicable surrounding for both unexperienced users and experts. These authors advise the application of different GSA methods to the same problem for at least two reasons. Firstly, as methods differ in their ability to address specific questions (e.g. input ranking, screening, mapping, analysis of individual contributions or of interactions), the insights provided by several methods can complement each other so that a more complete picture of the problem at hand is obtained. Secondly, since methods rely on different assumptions (e.g. linear/non-linear input-output relationship, skewed/non-skewed output distribution) whose degree of validity is sometimes not clearly defined, the application of multiple methods is a practical way to validate, reject or reinforce the conclusions of GSA. Zhan et al. (2013) propose an efficient integrated approach that integrates a qualitative screening method (the Morris method) with a quantitative analysis method based on the statistical emulator (variance-based method with the response surface method, named the RSMSobol method) to reduce the computational burden of GSA for time-consuming models. Xu et al. (2016) developed two modules for parameter sensitivity analysis and inverse estimation, respectively. In addition, a new solute transport module with numerically stable schemes was developed for ensuring stability of the combination of the models SWAP for soil moisture flow and EPIC for crop growth. Their method was tested and validated with a two-year dataset in a wheat growing field. Fourteen parameters out of the forty-nine total input parameters were identified as the sensitive parameters. These parameters were first inversely calibrated by using a numerical case, and then the inverse calibration was performed for the real field experimental case. Their research indicates that the proposed global method performs successfully to find and constrain the highly sensitive parameters efficiently that can facilitate application of the SWAP-EPIC model. Several adaptations of the Sobol method have been published, see e.g. Jansen (1999), Saltelli (2002), Sobol et al. (2007), Saltelli et al. (2010), Janon et al. (2014), Le Gratiet et al. (2014).
In the present study we applied the method described by Jansen (1999) (hereafter referred to as Sobol-Jansen) because this estimator is good for large first-order indices, and for both large and small total indices of the Sobol variance-based method.
Despite it is a sound advise to use several methods on the same dataset, we decided to apply only the Sobol-Jansen method in our project because we just want to show that performing sensitivity analyses is simple when using the appropriate tools.

Latin hypercube sampling (LHS)
Latin hypercube sampling (LHS) is a statistical method for generating a near-random sample of parameter values from a multidimensional distribution. The sampling method is often used to construct computer experiments or for Monte Carlo integration. The LHS has been described in McKay et al. (1979). An independently equivalent technique was proposed by Eglajs and Audze (1977). It was further elaborated by Iman et al. (1980).
In the context of statistical sampling, a square grid containing sample positions is a Latin square if (and only if) there is only one sample in each row and each column. A Latin hypercube is the generalisation of this concept to an arbitrary number of dimensions, whereby each sample is the only one in each axis-aligned hyperplane containing it.
When sampling a function of N variables, the range of each variable is divided into M equally probable intervals. M sample points are then placed to satisfy the Latin hypercube requirements; note that this forces the number of divisions, M, to be equal for each variable. Also note that this sampling scheme does not require more samples for more dimensions (variables); this independence is one of the main advantages of this sampling scheme. Another advantage is that random samples can be taken one at a time, remembering which samples were taken so far.

Hardware
The majority of the SWAP computations described in this paper were performed on the High Performance Computer (HPC) of Wageningen University and Research. The other SWAP computations were performed on a Paradigit computer with an Intel i7-3770 K CPU, running at 3.5 GHz and with 8 cores. The operating system on this machine was Ubuntu Linux version 18.04. The program for the computation of the storage capacity of a soil profile was developed and applied on a Lenovo W530 laptop with an Intel i7-3740QM CPU, running at 2.70 GHz and having 8 cores as well. The laptop was running Fedora Linux 28 in an Oracle VirtualBox under Windows 7.

Software
The model SWAP (Kroes et al., 2017) has been written in Fortran. Source code and executable (Windows) are freely available from the internet (swap.wur.nl). The program for the storage capacity has been written in Julia (Version 0.6.3). The generation of input files, analysis of output and creation of charts has been done in R with the additional packages ggplot2, directlabels, latex2expp, lhs and sensitivity. Fig. 4 presents the data flow for the computations: a) in case of a single computer and b) in case of computations on the HPC and control and processing on the personal computer. The software for submitting a specified number of jobs on the HPC was written in Julia.

The considered cases
To show the power and some of the capabilities of the methods described above, we selected 5 cases from earlier studies. The advantage is that the parameterization was done already, so we only had to perform the sensitivity and uncertainty analyses. A short description of the cases is presented below. The soil physical properties of the soil are described according to Damiano (2008) as presented in Eqs. (9)-(13) with the parameters from the Arrecifes soil: a = 10.8 kPa, I = 17.5 kPa, b = 6.94, D = 0.025, K s = 0.045 md −1 and S = 0.451. These parameters are used to compute the storage capacity of the soil profile with a groundwater depth of 2.20 m and a steady flux of 2 mm d −1 . The influence of each parameter is now estimated with the Sobol-Jansen procedure (Jansen, 1999). The possible changes are up to 5% 1 . A simple R procedure was developed to read the default values of the soil physical parameters, generate new parameter sets, run the program for each parameter set, read the storage capacity and finally analyse the data.

The influence of
S and K S of a soil in the Argentinian pampa's To show the power of the relatively simple and straightforward Impact Response surfaces (IRS), we consider the work presented by Kroes (2018) and Kroes et al. (2019). According to these authors, groundwater recharge results from the dynamic interaction between climate, land use and soil hydrology as it occurs in the critical zone, a thin portion of the biosphere connecting the lithosphere, atmosphere and hydrosphere. The variation of this interaction in space and time can be very large and is influenced by human and natural activities. Especially in flat, poorly drained sub-humid plains, such as the Argentine Pampas. Local observations of crop and soil parameters are used with special attention to a detailed determination of soil hydraulic properties. The soil parameters were obtained from Lozano et al. (2014) who published hydraulic properties based on observations in dominant soils in the Pampas. These input data were read into the SWAP model and simulations were performed with 26 years of local meteorological data. For this case, we investigated the influence of the saturated conductivity K s and the moisture content at saturation s in Eqs. (5) and (7) on the output of the model. For these computations we developed three scripts in R. The first script generates the input files for the SWAP runs by generating a regular grid with K s and s combinations. It was assumed that the groundwater level was deep and constant. Because we knew from experience that a fixed groundwaterlevel may yield unrealistic fluxes through the bottom of the profile, we assumed a deep aquifer with a constant head. Water flow between the profile and the aquifer was assumed to take place over a resistance of 100 d. To see the influence of the (constant) head in the aquifer on the output as well, data was prepared for several head values. These computations were performed on the HPC. After performing the computations, the output files were read by the second R script and the output values under considerations were stored in a file on the local computer where they were read by a third script that created the desired plots using the R packages ggplot2 and directlabels,

The influence of some parameters of an Argentinian soil
In this example we took the data presented by Kroes (2018) and Kroes et al. (2019) again. Soybeans were assumed to grow each year from Nov. 14th until May 5th. A heterogeneous soil with 7 soil layers was assumed. We were interested in the influence of the parameters of the top soil layer (20 cm thick) on the terms of the water balance. Therefore, we applied the Sobol-Jansen procedure as described by Jansen (1999) on the parameters under consideration. These parameters were n K , , , s s and l (see Eqs. 5 and 7). Beside these parameters, we included the parameter c , which is a critical stress index for compensation of root uptake. This parameter is crop-dependent and is read from the crop file. To analyze the influence of the groundwater level on the water balance terms, the groundwater level was made one of the input parameters to vary. The considered simulation period was from Jan. 1st, 19891st, to Dec. 31st, 2015 Because it was foreseen this case would require a lot of simulation runs, it was decided to do the computations on the HPC (method b in Fig. 4).

The influence of the uncertainty of soil parameters in a Brazilian soil
In Brazil, moisture contents were measured at 36 locations in an experimental area in Piracicaba during one year (Campos Oliveira et al., submitted for publication). The locations were chosen in such a way it was assumed all locations had similar soil properties with a 20 cm toplayer on a thick homogeneous subsoil. Meteorological data were obtained from the local station. The Mualem -Van Genuchten parameters n , , r s and were obtained from pressure plate measurements and K s was obtained from inverse modelling using water content measurement in the field with the Hydrus-1D model. The lvalue was fixed at 0.5. The soil moisture retention curves and hydraulic conductivity functions obtained this way are presented in Fig. 5. This figure shows there is a large variation of values for the considered soils. Specially the soil moisture retention curves may differ considerably.
It was expected there would be a correlation between the values of the parameters. The computed correlations are presented in Table 2. As expected, there are both positive and negative correlations. The most striking one is the correlation between r and n. This correlation coefficient is 0.8908 for the top soil and 0.945 for the subsoil.
In the original research we assumed a free outflow as the bottom boundary condition of the simulated profile. Because we wanted to see the influence of the uncertainty of the soil physical data on the calculated groundwater level as well, in the present study a deep aquifer was assumed to be at the bottom boundary of the soil profile. The head in the aquifer was −500 cm and there was a connection with the groundwater over a resistance of 250 d.
To show the influence of these parameter distributions, we applied the Latin Hypercube Sampling technique.

The contribution of soil parameters on the model output for a Brazilian soil
To find out which of the 10 parameters considered in the previous paragraph had the largest influence on the considered water balance terms, we applied the Sobol-Jansen technique once more. Contrary to the previous paragraph, where we considered the growing season only, now we considered the entire year.

Results and discussion
The results of the analyses described above will be presented in part II of this paper.

Software availability
All software applied to obtain the results presented in this paper, is available on Github (https://git.wur.nl/wesse016/unsens/).

Declaration of Competing Interest
The authors declare that there is no conflict of interest.