A User-Friendly Tool to Characterize the Moisture Transfer in Porous Building Materials: FLoW1D

This paper presents a user-friendly tool—FLoW1D (One-Dimensional Water Flow)—for the estimation of parameters that characterize the unsaturated moisture transfer in porous building materials. FLoW1D has been developed in Visual Basic for Applications and implemented as a function of the well-known Microsoft Excel© spreadsheet application. The aim of our work is to provide a simple and useful tool to improve the analysis and interpretation of conventional tests for the characterization of the hygric behavior of porous building materials. FLoW1D embraces the conceptual model described in EN 15026 for moisture transfer in building elements, and its implementation has been verified and validated correctly. In order to show the scope of the code, an example of an application has been presented. The hygric characterization of the limestone that is mostly employed in the Cathedral of Santa Maria and San Julian in Cuenca (Spain) was conducted based on an analysis of the conventional water absorption by capillarity tests (EN 15801).


Introduction
The conservation of architectural heritage is a social and economic responsibility, both for its cultural value and its value as a source of wealth for society [1]. The conservation of architectural heritage mainly focuses on the maintenance and rehabilitation of the materials used in the construction of these heritage elements-generally, porous building materials (PBMs). The degradation of these materials is usually caused by the transport of moisture through their pores in many environments [2,3]. Among others, different pathologies associated with water transport can be highlighted, such as the appearance of saline efflorescences and the development of alveolarization processes, as well as the fragmentation of the stone material due to the generation of stress associated with freezing processes or volume changes [4].
For this reason, it is necessary to determine the hygric properties of PBMs before conducting any action aimed at their conservation [5], as shown in standard UNE 41810 [6]. This standard shows the criteria and methodologies of intervention for stone materials and is focused on the conservation of cultural heritage; this is mainly defined as those activities framed within the processes of the stabilization of pathologies and preventive conservation strategies. At the European level, the EN 16515 [7] standard regulates the guidelines for the characterization of natural stone used in cultural heritage. The characterization of these materials is an essential stage in the correct definition of a conservation strategy that includes the definition of possible interventions, whether repairing or even replacing damaged stones. The standard EN 16515 [7] includes methodologies to conduct mineralogical, chemical, physical and mechanical characterization. The execution and interpretation of the tests proposed by EN 16515 [7] to specifically define the hygric behavior of PBMs is not easy; the use of numerical tools that are capable of simulating the moisture transfer in unsaturated porous media has become necessary for a satisfactory interpretation of the experimental tests.
A wide range of numerical tools is currently available to simulate the hygrothermal behavior of porous materials. Revisions by Delgado et al. [8] and Hens [9] or the multiphysics models developed in recent years [4,[10][11][12] highlight this fact. These tools adopt a well-structured conceptual framework that is described, for example, in EN 15026 [13], which incorporates the conceptual model of the main physical processes that condition the moisture and energy transfer in building elements. However, regardless of the conceptual level of the model, its practical scope is usually conditioned by the material parameters (the parameters of the state functions used in the modeling of the transfer of moisture and energy). It is therefore of the utmost importance to have a solid experimental basis for its estimation. In the current state-of-the-art approaches, the standard experimental tests focused on obtaining the hygric properties of the porous building materials used in architectural heritage are not aimed at triggering processes characterized by a single material parameter; in contrast, they are generally tests in which several physical processes are coupled together. This is the case, for example, for the water absorption by capillarity (WAC) test, described in EN 15801 [14], in which the coupled flow of both liquid water and vapor is present. In addition, the magnitudes that the international standards indicate as parameters for estimation (in the case of the WAC test, the water absorption by capillarity coefficient) are not usually material parameters but behavioral indices of these materials.
However, it does not seem reasonable to redefine a experimental framework that is deeply rooted in the technical community. It seems more logical to obtain the experimental data and to propose new strategies in order to estimate the parameters of the state functions considered in standard EN 15026.
There is extensive experience in the field of inverse problems and parameter estimation. The works of Carmeliet and Roels [15], Roels et al. [16], Gómez et al. [17], Galván et al. [18], Amirkhanov [19] and Rouchier et al. [20] are examples of parameter estimation from laboratory tests. One of these approaches could be adopted together with the freeware codes included in the review by Delgado et al. [8], or any other numerical solver, for the solution of the direct problem (simulation). However, a numerical tool should be as friendly as possible for the technicians in charge of the characterization of the PBM or even other partially saturated porous materials such as concrete, mortar, rammed earth or soils. For this reason, the program FLoW1D (One-Dimensional Water Flow) has been developed in the environment of the well-known Microsoft Excel© worksheet application using the Visual Basic for Applications (VBA) language. This allows users to modify the code according to their needs. However, if they do not wish to do so, the developed program can be used with the same simplicity as a worksheet function. In this way, the developed code can be used in conjunction with one of the many available optimization modules in Microsoft Excel© (such as Solver), facilitating the parameter estimation process.
In this work, the fundamentals of the conceptual model adopted for the analysis of isothermal tests are described. Additionally, we have also defined the numerical implementation in VBA and the interface in the Microsoft Excel© environment. The tool has been completely checked with three verification exercises. In order to illustrate the scope of the code, an example of an application has been included; the proposed case features an estimating of the unsaturated moisture transfer parameters (intrinsic permeability and water retention curve) of limestone used in the construction of the Cathedral of Santa María and San Julián in Cuenca, Spain.
Finally, it is important to note that the numerical tool has been made freely available to the scientific and technical community (Supplementary Material). It can be used in its current state or, if desired, it can be modified to suit other needs such as different sample sizes, discretization schemes or even for the application of alternative boundary conditions (as described below in Section 3.2).

Experimental Test
The WAC test has been selected for the parameter estimation exercises, as other authors have done for cement-lime mortars containing phase change materials [21]. In this well-known test, there are processes of the coupled transport of water in both the liquid phase and vapor in a partially saturated porous medium. The test consists of placing a sample of the selected PBM on a saturated permeable bed so that the water rises by capillarity. The increments of the water mass absorbed in each sampling time interval are measured, and the WAC coefficient-a performance index of the tested materialis calculated. However, the same data can be used to estimate the material parameters that characterize the unsaturated moisture transfer (intrinsic permeability and the water retention curve parameters) by the analysis of the inverse problem.

Conceptual Model
The formulation used for balance equations, field equations and the structure of state functions is based on the conceptual model described in EN 15026. For the simulation of WAC tests, an isothermal and one-dimensional flow has been assumed. In addition, the PBM is considered non-deformable.
The water mass balance, without considering water sources or sinks in the domain, is defined by the expression ∂w ∂t where w is the water content per total volume unit (kg m −3 ), g is the water flow rate (kg m −2 s −1 ) and ∇· is the divergence differential operator. The matric suction (defined as the difference between gas pressure and liquid pressure, s = P G − P L ) has been selected as the state variable. In the water content term, the contributions of the liquid, w w , and vapor, w v , phases have been considered according to where φ is the porosity (m 3 m −3 ), Sr is the degree of saturation (m 3 m −3 ), θ is the volumetric water content (m 3 m −3 ), θ sat is the volumetric water content in saturated conditions (m 3 m −3 ), ρ w is the liquid water density (kg m −3 ) and ρ v is the water vapor density (kg m −3 ), estimated by the psychrometric equation (Equation (A2)). The water flow rate, g, is obtained from the contributions of the liquid, g w , and vapor, g v , phases: The liquid and vapor fluxes were defined by adopting the approaches from Cabrera et al. [22]. The term corresponding to the liquid phase is defined as where k is the hydraulic conductivity (s) and ∇s is the suction gradient (Pa m −1 ). The hydraulic conductivity has been calculated as where γ w is the specific weight of the water (N m −3 ), κ is the relative permeability (dimensionless) and K sat is the saturated conductivity (m s −1 ), defined as where K is the intrinsic permeability (m 2 ). The relative permeability has been determined by the van Genuchten-Mualem [23] approach: where m (dimensionless) is a fitting parameter used in the van Genuchten [24] water retention curve: where α (Pa −1 ) and n (dimensionless) are also van Genuchten model parameters.
The water vapor flow rate is calculated as where δ p is the vapor permeability (kg Pa −1 m −1 s −1 ) and ∇P v is the water vapor pressure gradient (Pa m −1 ). The vapor permeability is calculated by where δ 0 (kg Pa −1 m −1 s −1 ) is the vapor permeability of the still air and µ (dimensionless) is the diffusional resistance factor. The vapor permeability of the still air, δ 0 , is calculated by the expression where MM w is the molar mass of the water (kg mol −1 ), R is the universal gas constant (Pa m 3 mol −1 K −1 ), T is the absolute temperature (K) and D v is the binary diffusion coefficient of the water vapor (m 2 s −1 ), which is calculated as in [25] D v = 5.9 × 10 −6 T 2.3 P G where P G is the gas pressure (Pa), which is assumed to be constant and equal to the atmospheric pressure. The diffusional resistance factor, µ, can be calculated as in [25] where τ v (dimensionless) is the tortuosity.

Numerical Model and System Abstraction
To solve the boundary value problem, an explicit finite difference scheme has been used. A detailed description of the spatial and temporal discretization applied is included in Appendix C, but Figure 1 provides a scheme of the discretization of the domain used. The lower grid point (i = 1) corresponds to the base of the sample in contact with the saturated permeable layer. Consequently, a Dirichlet condition of suction equal to 0 Pa has been applied in this boundary. The boundary condition applied to the upper grid point (i = n) is estimated using the formulation described in EN 15026 for moisture transport across the interfaces between the material and environment.
In the WAC tests, the specimens are kept in a desiccator with a controlled relative humidity (RH) of 32% before being tested. Therefore, the initial condition set is the suction corresponding to this RH, which is calculated using the psychrometric law (Equation B2). The numerical scheme described in Appendix C has been programmed using the Visual Basic for Application (VBA) language and implemented as a Microsoft Excel© function called FLoW1D. Figure 2 shows a snapshot of the tool interface. The code offers the possibility of calculating the relative permeability by means of the van Genuchten [24] or the Brooks and Corey [26] approaches. The van Genuchten model is the default option. If the user wants to use the Brooks and Corey model, the value of λ (empirical parameter of the model) must be included in the input of variables (cell C3 in Figure 2a).
The arguments of FLoW1D (Par, Tini, Tfin, So) are entered in the same way as in a predefined Excel function, as shown in Figure 2b. The Par argument is the vector (vectors are identified by using bold font throughout the paper) of material parameters and geometric dimensions (cells C1 to C8 in Figure 2a). The "Tini" and "Tfin" arguments correspond to the initial and final time of the time increment modeled. The So argument is the vector corresponding to the suction at the "Tini" time. Given that a match is made between the spreadsheet cells and the grid points of the discretization, the first value for all components of So (the initial condition of the boundary value problem) must be defined in cells H3 to H103 in Figure 2b (if the domain is discretized in 101 grid points). By default, the value computed in cell C12, obtained from the initial RH after applying the psychrometric law (Equation B2), is assumed as the first value of So. For the following computational times (defined in cells I2, J2, K2 and so on), So is the result of the previous time step. The boundary conditions are defined as indicated in Figure 1. For all computational times, a Dirichlet boundary condition is assumed at the bottom grid point. Therefore, cells H3, J3, K3 and so on are equal to the value "s1" defined in cell C14 (Figure 2a). If the psychrometric law (Equation B2) is applied, s1 is a function of the RH value consigned in cell C13 and consequently is sufficient to change that value to modify the suction value at grid point 1 throughout the computational process. Although a Dirichlet condition has been assumed by default in grid point 1, FLoW1D has the option to activate a von Neumann Sample Sat urated permeable layer RH initial = 32% Figure 1. System abstraction discretization, boundary and initial conditions. RH: relative humidity.
In the WAC tests, the specimens are kept in a desiccator with a controlled relative humidity (RH) of 32% before being tested. Therefore, the initial condition set is the suction corresponding to this RH, which is calculated using the psychrometric law (Equation (A2)).
The numerical scheme described in Appendix C has been programmed using the Visual Basic for Application (VBA) language and implemented as a Microsoft Excel© function called FLoW1D. Figure 2 shows a snapshot of the tool interface. The code offers the possibility of calculating the relative permeability by means of the van Genuchten [24] or the Brooks and Corey [26] approaches. The van Genuchten model is the default option. If the user wants to use the Brooks and Corey model, the value of λ (empirical parameter of the model) must be included in the input of variables (cell C3 in Figure 2a).
The arguments of FLoW1D (Par, Tini, Tfin, So) are entered in the same way as in a predefined Excel function, as shown in Figure 2b. The Par argument is the vector (vectors are identified by using bold font throughout the paper) of material parameters and geometric dimensions (cells C1 to C8 in Figure 2a). The "Tini" and "Tfin" arguments correspond to the initial and final time of the time increment modeled. The So argument is the vector corresponding to the suction at the "Tini" time. Given that a match is made between the spreadsheet cells and the grid points of the discretization, the first value for all components of So (the initial condition of the boundary value problem) must be defined in cells H3 to H103 in Figure 2b (if the domain is discretized in 101 grid points). By default, the value computed in cell C12, obtained from the initial RH after applying the psychrometric law (Equation (A2)), is assumed as the first value of So. For the following computational times (defined in cells I2, J2, K2 and so on), So is the result of the previous time step. The boundary conditions are defined as indicated in Figure 1. For all computational times, a Dirichlet boundary condition is assumed at the bottom grid point. Therefore, cells H3, J3, K3 and so on are equal to the value "s1" defined in cell C14 (Figure 2a). If the psychrometric law (Equation (A2)) is applied, s1 is a function of the RH value consigned in cell C13 and consequently is sufficient to change that value to modify the suction value at grid point 1 throughout the computational process. Although a Dirichlet condition has been assumed by default in grid point 1, FLoW1D has the option to activate a von Neumann boundary condition, such as that indicated for the grid point n in Figure 1. If the user activates this optionand since the term corresponding to the liquid phase is assumed to equal to zero-the flow g B0 applied to grid point 1 is given, according to the EN 15026 standard, by the expression where s d,0 is the equivalent vapor diffusion thickness of bottom interfaces, P v,0 is the vapor pressure obtained by the relative humidity at the environment in the bottom interface and P v,1 is the vapor pressure obtained by the psychrometric law (Equation (A2)) at grid point 1. By default, at the top of the domain (grid point n in Figure 1), the von Neumann-type condition applies. In this case, s d,n , P v,n and P v,an are used instead of s d,0 , P v,0 and P v,1 , respectively. FLoW1D has also been prepared to change this top boundary condition, and it is able to shift to a Dirichlet condition. To do this, cells H103, I103, J103 and so on must be equal to "sn" (cell C16), as a function of the RH value consigned in cell C15.

Verification Process
A comprehensive verification of the tool has been carried out; here, we present only three of the performed exercises. In Table 1, the parameters and initial and boundary conditions adopted in each exercise are shown. Therefore, in the default option (bottom, Dirichlet condition; top, von Neumann condition), FLoW1D will be activated for the cell range 3 to 103. According to the initial calculation proposal, 101 grid points have been adopted for the spatial discretization. The user can easily increase or reduce this number. The activation of FLoW1D for a column of cells, as an "array formula" (or "CSE formula"), requires the selection of the whole range of solution vectors and for the user to simultaneously press the keys Control, Shift and Enter [27]. Thus, the results of the FLoW1D function are the suction in each of the grid points in which the sample has been discretized for the desired computational times (see Figure 2b).
Although the state variable used is suction, the magnitude measured directly in the WAC test is the water content. To obtain the value of w associated with the value of s in each grid point, the function ws was implemented. In this function, besides s, the arguments are porosity φ (cell C2), the van Genuchten parameters n, m and α (cells C5, C6 and C7, respectively), and the temperaturê T (cell C8). The results of this function are shown in Figure 2c. Finally, it should be noted that the total mass of water contained in the sample is w tot,cal ; this variable represents the variation of the accumulated water mass (g) in the tested material specimen, m w (Figure 2d). This calculated variable is compared with the experimentally obtained value of w tot,obs . In summary, the main equation of FloW1D is Equation (1), which represents the water mass balance in the material. The state variable of this balance is w (water content per total volume unit). From the values of w at each grid point, it is possible to obtain the total mass of water accumulated in the tested sample (m w ) considering the total volume of the sample and the observation times.

Verification Process
A comprehensive verification of the tool has been carried out; here, we present only three of the performed exercises. In Table 1, the parameters and initial and boundary conditions adopted in each exercise are shown.
3. In the first two verification exercises (QE1 and QE2), a flux was imposed at the top of the domain, maintaining perfect drainage at the bottom. In the first case, a linear retention curve and constant hydraulic conductivity were assumed. These simplifications allowed the use of an analytical solution [28]. Figure 3a shows the satisfactory fit obtained with FLoW1D (with a mean absolute relative error (MARE) of 0.004%). For the second example presented, the same physical problem has been simulated, but the van Genuchten water retention curve and the non-constant hydraulic conductivity model defined in Table 1 were adopted. As a reference, this physical problem has been simulated using the well-known software SEEP/W-Geoslope ® [29]. Figure 3b shows the good fit obtained (MARE = 0.22%).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 18 In addition to the liquid flow, the implementation of vapor transport has also been verified. Thus, among other exercises, benchmark 2 [30] defined in the benchmark report generated in the international project HAMSTAD [31] (a project focused on the standardization of formulation of heat, air and moisture transfer models) has been simulated. This benchmark-named exercise QE3 in Table  1-studies the 1D isothermal water flow in a reference material caused by a sudden change of the environmental RH. The results in Figure 4 show the good fit obtained (MARE = 1.77 %), giving confidence regarding the use of FLoW1D to characterize the water transport in porous building materials.

Model Application
In order to illustrate the coupled use of FLoW1D with optimization modules in the Microsoft In addition to the liquid flow, the implementation of vapor transport has also been verified. Thus, among other exercises, benchmark 2 [30] defined in the benchmark report generated in the international project HAMSTAD [31] (a project focused on the standardization of formulation of heat, air and moisture transfer models) has been simulated. This benchmark-named exercise QE3 in Table 1-studies the 1D isothermal water flow in a reference material caused by a sudden change of the environmental RH. The results in Figure 4 show the good fit obtained (MARE = 1.77%), giving confidence regarding the use of FLoW1D to characterize the water transport in porous building materials. In addition to the liquid flow, the implementation of vapor transport has also been verified. Thus, among other exercises, benchmark 2 [30] defined in the benchmark report generated in the international project HAMSTAD [31] (a project focused on the standardization of formulation of heat, air and moisture transfer models) has been simulated. This benchmark-named exercise QE3 in Table  1-studies the 1D isothermal water flow in a reference material caused by a sudden change of the environmental RH. The results in Figure 4 show the good fit obtained (MARE = 1.77 %), giving confidence regarding the use of FLoW1D to characterize the water transport in porous building materials.

Model Application
In order to illustrate the coupled use of

Model Application
In order to illustrate the coupled use of FLoW1D with optimization modules in the Microsoft Excel© environment, WAC tests of the limestone used in the construction of the Cathedral of Santa María and San Julián in Cuenca, Spain were analyzed. Two different lithotypes of this limestone have been identified: the first lithotype is used in ornamental elements and named ornamental porous stone (OPS), and the second lithotype is a stone used in structural elements, denominated as structural porous stone (SPS). The chemical analysis (see Table S1) performed by X-ray diffraction (see Figure S1) shows that calcite is the major material in both lithotypes, along with a significant fraction of quartz (Q peak in Figure S1b) in SPS. The real and apparent densities [32] are 2746 and 1498 kg m −3 , respectively, for OPS. In the case of SPS, these densities are 2623 and 1886 kg m −3 , respectively. Figure 5(a1,a2) shows samples of the two tested materials. In the image of the optical microscope ( Figure 5(b2)), ooliths of an approximate size of 10 µm are observed with a rounded shape and low crystallinity index, with small regular calcite envelopes whose nuclei are usually formed by quartz clasts, characterizing SPS. Furthermore, the distribution of pore sizes ( Figure S2) obtained in the Mercury Intrusion Porosimetry (MIP) test shows that the mean pore size ranges between 1 and 20 µm for OPS and 0.1 and 10 µm for SPS [33]. Finally, regarding experimental characterization, Table 2 includes the WAC coefficients obtained according to EN 15801 [14].
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 18 fraction of quartz (Q peak in Figure S1b) in SPS. The real and apparent densities [32] are 2746 and 1498 kg m -3 , respectively, for OPS. In the case of SPS, these densities are 2623 and 1886 kg m -3 , respectively. Figures 5(a1) and 5(a2) shows samples of the two tested materials. In the image of the optical microscope ( Figure 5(b.2)), ooliths of an approximate size of 10 μm are observed with a rounded shape and low crystallinity index, with small regular calcite envelopes whose nuclei are usually formed by quartz clasts, characterizing SPS. Furthermore, the distribution of pore sizes ( Figure S2) obtained in the Mercury Intrusion Porosimetry (MIP) test shows that the mean pore size ranges between 1 and 20 μm for OPS and 0.1 and 10 μm for SPS [33]. Finally, regarding experimental characterization, Table 2 includes the WAC coefficients obtained according to EN 15801 [14].  In order to carry out the parameter estimation, FLoW1D was used with the predefined complement Solver of Microsoft Excel©. The Generalized Reduced Gradient (GRG) was the selected optimization option of Solver. This algorithm, developed by Abadie and Carpentier [34], is an extension of the Fank and Wolfe [35] reduced gradient method. Its mathematical structure has undergone continuous development [36][37][38][39], but the GRG2 version of this algorithm is the one implemented in Solver [40,41].
A least squares objective function was adopted. The intrinsic permeability, K, and the parameters of the van Genuchten retention curve, n and m, were the estimated material parameters. The porosity of the two samples tested was determined by ISO 17892-1 [42], obtaining results of 0.38 and 0.23 for OPS and SPS, respectively. A tortuosity value equal to 1 was assumed, in agreement with the work of Philip and de Vries [43].
The parameter α of the van Genuchten retention curve can be obtained by the relationship  In order to carry out the parameter estimation, FLoW1D was used with the predefined complement Solver of Microsoft Excel©. The Generalized Reduced Gradient (GRG) was the selected optimization option of Solver. This algorithm, developed by Abadie and Carpentier [34], is an extension of the Fank and Wolfe [35] reduced gradient method. Its mathematical structure has undergone continuous development [36][37][38][39], but the GRG2 version of this algorithm is the one implemented in Solver [40,41].
A least squares objective function was adopted. The intrinsic permeability, K, and the parameters of the van Genuchten retention curve, n and m, were the estimated material parameters. The porosity of the two samples tested was determined by ISO 17892-1 [42], obtaining results of 0.38 and 0.23 for OPS and SPS, respectively. A tortuosity value equal to 1 was assumed, in agreement with the work of Philip and de Vries [43].
The parameter α of the van Genuchten retention curve can be obtained by the relationship where the s o is the initial suction and Sr o is the degree of saturation at the beginning of the WAC test obtained from the water mass in Equation (2): where w o is the water content at the begging of the WAC test. Figure 6 shows the results of two WAC tests carried out with OPS and SPS, respectively, and the simulation developed with the tool FLoW1D.
where the o s is the initial suction and o Sr is the degree of saturation at the beginning of the WAC test obtained from the water mass in Equation (2): where o w is the water content at the begging of the WAC test. Figure 6 shows the results of two WAC tests carried out with OPS and SPS, respectively, and the simulation developed with the tool FLoW1D. As can be observed in Figure 6, a satisfactory fit was obtained for both tests. The estimated values of K, n and m, are indicated in Table 3. This table also presents the parameter α, which is calculated using Equation (15) using estimated values of n and m.

Conclusions
This work presents FLoW1D, a numerical tool designed as a support for the analysis of the conventional experimental tests for the hygric characterization of porous building materials. The code simulates the one-dimensional moisture transfer in a partially saturated porous media. It has been verified against other widely used commercial codes and contrasting benchmarks, using other porous materials as a reference.
To increase its usability, FLoW1D has been developed as a function in the environment of the As can be observed in Figure 6, a satisfactory fit was obtained for both tests. The estimated values of K, n and m, are indicated in Table 3. This table also presents the parameter α, which is calculated using Equation (15) using estimated values of n and m.

Conclusions
This work presents FLoW1D, a numerical tool designed as a support for the analysis of the conventional experimental tests for the hygric characterization of porous building materials. The code simulates the one-dimensional moisture transfer in a partially saturated porous media. It has been verified against other widely used commercial codes and contrasting benchmarks, using other porous materials as a reference.
To increase its usability, FLoW1D has been developed as a function in the environment of the well-known Microsoft Excel© spreadsheet application using the Visual Basic for Applications language. This aspect provides flexibility to the numerical tool, which could be modified easily for the analysis of the hygric characterization of other partially saturated porous materials whose hygric behavior can be reproduced by the implemented state functions. Additionally, the computational method used is relatively simple (an explicit finite difference scheme has been used); therefore, no in-depth knowledge is required to make changes to the FLoW1D code. Finally, it is important to highlight the fact that FLoW1D is presented as an application developed in an Microsoft Office platform, providing greater accessibility to the scientific and technical community as the Microsoft Excel© spreadsheet application is a very common working environment used in laboratories.
The results presented in this work show that FLoW1D allows both the simulation of two WAC tests and the estimation of hygric parameters (the intrinsic permeability and associated parameters of the water retention curve) with robustness and ease when it is used in conjunction with the well-known Microsoft Excel© complement Solver. Therefore, FloW1D allows test analysis and parameter estimation to be performed in the same environment.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix B. Water Properties
The density of water vapor, ρ v , is calculated by the law of ideal gases: where MM w is the relative mass of water (kg mol −1 ), R is the ideal gas constant (Pa m 3 mol −1 K −1 ), T is the absolute temperature (K) and P v is the water vapor pressure (Pa), obtained from the psychrometric equation [44] as a function of suction, s.
where, if the vapor is an ideal gas, the vapor pressure in saturated conditions is calculated as where ρ v,sat is the density of water vapor under saturated conditions. It is obtained by the expression [45] ρ whereT is the temperature ( • C). Finally, the dynamic viscosity of the water (N s m −2 ) is defined as [45] µ w = 0.6612 (T − 229) −1.562 (A5) where T is the absolute temperature (K).

Appendix C. Numerical Model
To solve the boundary value problem, an explicit centered finite difference scheme has been developed. The time derivative of the water content is expressed as On the other hand, the water flow rate (a scalar magnitude in this 1D model) can be written as The water flow rate is calculated at the center of two consecutive grid points; i.e., in i + 1 2 and i − 1 2 (see Figure A1). Therefore, it can be assumed that g j−1 The described formulation constitutes a forward-time central-space (FTCS) scheme (with an evaluation of the spatial differences at the midpoints of the grid for a better evaluation of the non-linear storage term).
The computational time interval is defined by the initial and final time; i.e., "Tini" and "Tfin". The time step of the interval, ∆t, was defined in the code, which was an arbitrarily chosen value and remained constant throughout the calculation process. The time was updated after each time step until the final value of the interval ("Tfin") was reached, when the calculation process ended.
The described formulation constitutes a forward-time central-space (FTCS) scheme (with an evaluation of the spatial differences at the midpoints of the grid for a better evaluation of the nonlinear storage term).
The computational time interval is defined by the initial and final time; i.e., "Tini" and "Tfin".
The time step of the interval, Δt , was defined in the code, which was an arbitrarily chosen value and remained constant throughout the calculation process. The time was updated after each time step until the final value of the interval ("Tfin") was reached, when the calculation process ended. Figure C1. Definition of the points at which the water flow rate was calculated.

9.
Hens  Figure A1. Definition of the points at which the water flow rate was calculated.