Received: 4 January 2024

Accepted: 22 February 2024

DOI: 10.1002/vzj2.20326
Vadose Zone Journal

TECHNICAL NOTE
S p e c i a l S e c t i o n : Tr i b u t e t o R i e n v a n G e n u c h t e n , R e c i p i e n t o f t h e 2 0 2 3 Wo l f P r i z e f o r A g r i c u l t u r e

Implementation of the Brunswick model system into the Hydrus
software suite
Efstathios Diamantopoulos1
1 Chair of Soil Physics, University of

Bayreuth, Bayreuth, Germany
2 Department of Environmental Sciences,

University of California, Riverside,
California, USA
3 Soil Science Section, University of Kassel,

Witzenhausen, Germany
Correspondence
Efstathios Diamantopoulos, University of
Bayreuth, Bayreuth, Germany. Email:efstathios.diamantopoulos@uni-bayreuth.de

Jirka Simunek2

Tobias K. D. Weber3

Abstract
The Brunswick modular framework for modeling unsaturated soil hydraulic properties (SHP) over the full moisture range was implemented in the Hydrus suite. Users
can now additionally choose between four different variants of the Brunswick model:(i) van Genuchten–Mualem (VGM), (ii) Brooks–Corey, (iii) Kosugi, and (iv) modified van Genuchten. For demonstration purposes, simulation results for two different
setups, (i) bare soil evaporation and (ii) root water uptake, are presented, along with acomparison of the original VGM model and its Brunswick variant. Results show thatthe original VGM model underestimates the simulated cumulative evaporation andcumulative transpiration due to the inconsistent representation of the SHP in the dry
moisture range. We also implemented a two-step hydro-PTF (pedotransfer function)
into the Hydrus suite that converts the parameters of the original VGM model (fromRosetta) to the corresponding Brunswick variant. In that way, physically comprehensive simulations are ensured if no data on SHP are directly available, but information
on physical soil properties (e.g., texture and bulk density) exists.

1

INTRODUCTION

Accurate simulations of water dynamics in soils with theRichards–Richardson equation (Richards, 1931; Richardson,
1922) require the correct representation of the soil hydraulicproperties (SHP) from saturation to very dry conditions (Idenet al., 2021). According to Weber et al. (2023), we can easily obtain 200 different models of SHPs if we combine the
22 water retention models listed in Du et al. (2016) with the 9models of relative conductivity presented by Assouline and Or(2013). Among all those models, the most widely used model
for representing SHP is by far the van Genuchten–Mualem

Abbreviations: BC, Brooks-Corey; BW, Brunswick; Cum. Transp.,
cumulative transpiration; PTF, pedotransfer function; KS, Kosugi;
modVGM, modified van Genuchten; RWU, root water uptake; SHP, soil
hydraulic properties; VGM, van Genuchten-Mualem.

model (VGM) (Mualem, 1976; van Genuchten, 1980). During the last years, many different models have been proposedin the literature with the aim to extend or replace the VGM
model, with the focus either on the wet range or the dry
range (Fayer & Simmons, 1995; Lebeau & Konrad, 2010;
Peters, 2013; Vogel et al., 2000), among hundreds of others.
Such a system for modeling SHPs was presented by Weberet al. (2019), which, for the most popular model configuration, comes with a pedotransfer function to predict the model
parameters (Weber et al., 2020).
In this technical note, we present implementations of theBrunswick model variants (Streck & Weber, 2020) into theHydrus suite, Version 5. These models are: (i) van GenuchtenMualem (VGM) (Mualem, 1976; van Genuchten, 1980), (ii)
Brooks–Corey (BC) (Brooks & Corey, 1964), (iii) Kosugi(KS) (Kosugi, 1996a), and (iv) modified van Genuchten(modVGM) (Vogel & Cislerova, 1988). As an example, weThis is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original
work is properly cited.
© 2024 The Authors. Vadose Zone Journal published by Wiley Periodicals LLC on behalf of Soil Science Society of America.
Vadose Zone J. 2024;e20326.
https://doi.org/10.1002/vzj2.20326

wileyonlinelibrary.com/journal/vzj2

1 of 10

compared the performance of the classical VGM model tothe Brunswick variant (VGM-BW). For this, we first fitted
the two models to water retention and unsaturated hydraulicconductivity data for three soils with contrasting properties. Second, for the three selected soil types, we tested thedifference between the two SHPs models by conducting two
different simulations: (i) bare soil evaporation and (ii) root
water uptake (RWU). For the RWU simulations, we used two
different models: the Feddes model (Feddes et al., 1978) andthe Nimah–Hanks model (Nimah & Hanks, 1973). The maindifference between the two models is that in the Nimah–
Hanks model RWU is calculated using the unsaturated
hydraulic conductivity of the local soil and the pressure headgradient between the soil and the roots, whereas, in the Feddes model, the local pressure head determines actual uptake.
Third, we discuss the systematic differences in the simulation
results, identifying the relevance of including a comprehensive representation of the SHP in the dry range. Finally, we
present the equations for all BW variants implemented in
Hydrus in the Appendix.

2

MATERIAL AND METHODSIn the following paragraphs, we present the governing equation of water flow in soils, the different models for SHP, andthe simulation setups in Hydrus.

2.1

DIAMANTOPOULOS ET AL.

Vadose Zone Journal

Water flow in soils

Water fluxes in unsaturated isotropic soils are described by
the RE (Richards, 1931; Richardson, 1922) assuming RWU
is represented as a sink term:
𝜕𝐾(ℎ)
𝜕𝜃(ℎ)
= −∇ ⋅ 𝐪 − Γroots = ∇ ⋅ [𝐾(ℎ)∇ℎ] +− Γroots ,
𝜕𝑡
𝜕𝑧
(1)
where 𝑧 is the vertical coordinate [L], ℎ is the pressure head[L], 𝜃 is the volumetric water content [L3 L‒3 ], 𝐪 is the Darcy
flux [L T−1 ], K [L T−1 ] is the saturated/unsaturated hydraulicconductivity as a function of 𝜃 or ℎ, and Γroots is the sinkterm representing RWU [T−1 ]. Equation (1) requires that the
air pressure in the soil is equal to the atmospheric pressure(single-phase flow) and that a local equilibrium between watercontent and pressure head is always valid (Diamantopoulos &Durner, 2012) and is described by the water retention curve
𝜃(ℎ).

Core Ideas∙ A modular framework was implemented in Hydrus
to describe soil hydraulic properties.
∙ The new model is more robust in describing soil
hydraulic properties in the dry moisture range.
∙ Hydrus simulations with the new model predicthigh bare soil evaporation and transpiration fluxes.

2.2.1

The classical approach

In the classical approaches, the water retention curve is
described as a function of the pressure head:𝜃(ℎ) = 𝜃r + (𝜃s − 𝜃r )Φe (ℎ),

where Φe (ℎ) is the effective water saturation function, withvalues between 1.0 and 0.0, and 𝜃s [L3 L‒3 ] and 𝜃r [L3 L−3 ]are the saturated and residual water contents, respectively.
The unsaturated hydraulic conductivity is given as afunction of the pressure head by:𝐾(ℎ) = 𝐾s 𝐾r (ℎ),

Models for SHPs

In the following, we limit the description of the SHP to a minimum. The full mathematical expressions for all original andextended models are provided in Appendix A.1–A.4.

(3)

where 𝐾r [-] is the relative unsaturated hydraulic conductivity
function, which is scaled with the saturated hydraulic conductivity, 𝐾s [L T−1 ]. Full expressions for Φe (ℎ) and 𝐾r (ℎ) are
given in Appendix A.2–A.5 for all four models presented in
the introduction.

2.2.2

The VGM-BW model

According to Weber et al. (2019), the SHP in Equations (2)and (3) may be extended using the Brunswick model framework. Therein, the soil water retention and the unsaturated
hydraulic conductivity are partitioned into two domains. One
represents water stored and conducted in completely filledcapillaries, and another one water stored and conducted inpartially full pores. These two domains are, respectively,
named the capillary (suffix c) and the non-capillary (suffix nc) pore domains. The water retention function in BW isgiven as:𝜃(ℎ) = 𝜃sc 𝑆c (ℎ) + 𝜃snc 𝑆nc (ℎ),

2.2

(2)

(4)

where 𝜃sc [L3 L−3 ] and 𝜃snc [L3 L−3 ] are the saturated watercontent values of the respective domains. 𝑆c (ℎ) [-] and 𝑆nc (ℎ)[-] describe the effective water saturation as a function of thepressured head, respectively. The function 𝑆c (ℎ) is a rescaled

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20326, Wiley Online Library on [01/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

2 of 10

Vadose Zone JournalΦe (ℎ) by𝑆c (ℎ) =Φe (ℎ) − Φe (ℎ0 ).
1 − Φe (ℎ0 )(5)This ensures a water content of 0 at a fixed very low pressure
head, here ℎ0 = −106.8 cm, as it is described in Appendix A.1.
The unsaturated hydraulic conductivity function is given by
𝐾(ℎ) = 𝐾sc 𝐾rc + 𝐾snc 𝐾rnc .

Modeling soil evaporation

In Hydrus, evaporation dynamics at the soil surface are simulated by assuming a Neumann-type boundary condition. The
water flux at the soil surface, 𝑞top [L T-1 ], is set equal to thepotential evaporation flux (𝑞Pot. E. ) [L T-1 ], provided by theuser. However, if the ℎsoil (𝑧 = 0) [L] becomes equal or lessthan a critical pressure head ℎcrit [L], the boundary condition is switched to a Dirichlet-type boundary condition withℎsoil (𝑧 = 0) = ℎcrit [L]. In that way, the potential evaporation
rate is reduced to the actual evaporation rate, and it is calculated by applying the Darcy law to the soil surface. For more
information, we refer to Šimunek et al. (2006).

2.4

at each depth, calculated from Equations (3) and (6) for thetwo SHP models, respectively. The main difference between
the two models is that for the same values of ℎsoil (𝑧) and forthe same plant (the same stress function, 𝛼(ℎ)), the modelof Feddes et al. (1978) will calculate the same Γroots forany soil. This is different in the model of Nimah and Hanks(1973), in which the conductivity curve controls the RWU
calculations.

(6)

Again, c and nc stand for the capillary and non-capillary poredomains, 𝐾sc [L T-1 ] and 𝐾rc [-] are the saturated and relative hydraulic conductivities of the capillary pore domain, andlikewise, 𝐾snc [L T-1 ] and 𝐾rnc [-] the saturated and relative
hydraulic conductivities of the non-capillary pore domain.
Note that the full BW model framework includes a term for
isothermal vapor conductivity. In Hydrus, the option to simulate isothermal vapor conductivity is a standard feature thatcan be selected and, therefore, not included here.

2.3

3 of 10

2.5

Experimental data

We fitted the VGM and VGM-BW models to experimental
water retention and conductivity data for three different soils
from the recently published database of Hohenbrink et al.
(2023). We selected a very sandy soil (Sand) (99.8%, sampleID: “KIT_063”), a clay loam soil (CLoam) (37.4% sand,
34.1% loam, and 28.5% clay, sampleID: “TUBS_106”), and aclay soil (Clay) (66.4% clay, sampleID: “KIT_140”). Details
about the experimental protocol can be found in Hohenbrink
et al. (2023). We followed the fitting protocol for tabulated
SHPs data as in Diamantopoulos and Durner (2015) andPeters and Durner (2008). The parameters for the three soils
and two models are presented in Tables 1 and 2.

2.6

Simulation setup

2.6.1

Evaporation

We assumed a 10-cm-long soil profile with a uniform initialpressure head distribution, a zero flux at the bottom of thecolumn, and a potential evaporation rate of 0.5 cm day−1 at𝑧 = 0 cm. In total, we ran six simulations for three different soil types, parameterized with the VGM and VGM-BW
models.

Modeling RWU
2.6.2

We calculated Γroots in Equation (1) based on the two macroscopic approaches of Feddes et al. (1978) and Nimah andHanks (1973). In this technical note, we do not present
the models systematically but focus on the basic differences
between them. For full details about the two models, we referto the original publications of Feddes et al. (1978) and Nimahand Hanks (1973), respectively. The model of Feddes et al.
(1978) calculates Γroots , based on potential transpiration andon the soil pressure head ℎsoil(𝑧) at each numerical node, byapplying a stress function 𝛼(ℎ), which varies between 0 and 1(1: no water stress, 0: no RWU). In the model of Nimah andHanks (1973), Γroots is calculated by the difference between
the “effective” pressure head at the root surface, ℎroot (𝑧) andℎsoil (𝑧), and the unsaturated hydraulic conductivity of the soil

Root water uptake

We assumed a 150-cm long soil profile with a uniform initial pressure head distribution. We assumed a seepage faceboundary condition at the bottom of the profile and a linearly decreasing RWU distribution with depth. Similarly, as
for evaporation experiments, we ran simulations for three different soils and two RWU models. Finally, we assumed thepotential transpiration rate equal to 0.5 cm day−1 . The criticalpressure heads for the Feddes model were set to: 𝑃 0 = −10cm, 𝑃 𝑂𝑝𝑡 = −25 cm, 𝑃 2𝐻 = −320 cm, 𝑃 2𝐿 = −600 cm,
and 𝑃 3 = −16, 000 cm (Sugar Beat; Wesseling et al., 1991).
For the Nimah and Hanks model, the maximum value of thepressure head that the roots can apply was assumed to be equalto -16,000 cm or pF = 4.2 (wilting point).

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20326, Wiley Online Library on [01/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

DIAMANTOPOULOS ET AL.

DIAMANTOPOULOS ET AL.

Vadose Zone Journal

TA B L E 1Parameters for the van Genuchten–Muealem model for three different soils.

Soil type𝜽r [-]

𝜽𝒔 [-]

𝜶 [cm−1 ]𝒏 [-]𝑲s [cm day−1 ]

𝝀 [-]

Sand

0.059

0.306

0.0213

3.4

74.9

−0.4

CLoam

0.217

0.679

0.149

1.1

1000

−0.5

Clay

0.059

0.608

0.015

1.1

8.6

4.6

T A B L E 2 Parameters for the van Genuchten–Muealem model to the Brunswick variant model for three different soils. For all the simulations
with this model, we fixed 𝑎 = 1.5 and ℎ0 = 106.8 .
Soil type𝜽snc [-]

𝜽sc [-]

𝜶 [cm−1 ]𝒏 [-]𝑲sc [cm day−1 ]𝑲snc [cm day−1 ]

𝝀 [-]

Sand

0.063

0.242

0.021

3.5

74.3

0.0

−0.4

CLoam

0.400

0.257

0.111

2.8

71.0

0.17

−1.7

Clay

0.000

0.611

0.020

1.1

80.0

1e-3

10.0

2.7

Initial conditions

As the primary emphasis of this study lies in the representation of the dry range of SHPs, we deliberately selected distinct
initial pressure heads for the three soil types (green lines in
Figure 1). This ensures that any variations in the simulation
results arise exclusively from the representation of the dry
range in the two models and not from the wet range.

et al., 2006). For the conductivity data, both models perform
almost equally, with a distinct change in the slope around pF= 2.5 for the Sand, which is typical for corner and/film flow
(Diamantopoulos & Durner, 2015; Tuller & Or, 2001) components. Only the VGM-BW model can describe this change
in slope in the logK versus pF data, while the VGM modelcannot (Figure 1).

3.2
3
3.1

Simulations of bare soil evaporation

RESULTS AND DISCUSSION
Fitting of VGM and VGM-BW to data

Both the VGM and VGM-BW models describe the water
retention and conductivity data in the wet range very well
(Figure 1, Tables 1 and 2). Any disparities visible in the wet
range between the two models can be attributed to the absenceof conductivity data in the wet range, along with compensation made by the fitting algorithm to characterize data in thedry range (mostly retention data), accounting for limitations
inherent in the models, especially in the VGM model. However, the difference between the VGM and VGM-BW models
becomes apparent in describing the dry range of both curves,
since the VGM model predicts an almost constant and nonzero water content for pF >2.5 (h in cm) for the Sand, pF > 4
for the CLoam, and pF > 6.8 for the Clay, respectively. On theother hand, the VGM-BW model describes the retention datain the dry range very well, showing an almost linear decreaseof the water content versus pF in that range, guaranteeing bydefinition that water content will be zero around a value of pF= 6.8 (in cm) (“oven dryness,” (Blöcher et al., 2023)). It has
also been shown before that this dry component in the water
retention curves is required for simulating bare soil evaporation at the continuum scale (Iden et al., 2021, 2021; Saito

For all three soils, the simulations with the VGM-BW modelyielded a higher Cumulative Evaporation value (Cum. Evap.)
at the end of the simulation (Figure 2). The difference waslow (around 0.2 cm) for the Sand and CLoam soil and slightly
higher for the Clay soil (around 0.8 cm). For the VGM model,
the simulation for the Sand did not finish because Hydrus
experienced numerical problems due to the constant watercontent value after pF > 2.5 (Figure 2, top, red x sign).
Although the differences in Cum. Evap. were low, more
important implications is that for the Sand and CLoam soil,
the simulations predict a later time for stage 2 evaporation
to be reached when the water content at the top of the soil
equals 0.069 cm3 cm‒3 and 0.217 cm3 cm‒3 , respectively.
This is physically implausible since, under prolonged highatmospheric demand, we expect a near zero water content
at the top of the soil surface, with corresponding effects on
the water content profiles. We note that Iden et al. (2021)
have illustrated that a correct description of bare soil evaporation experiments of different soils with isothermal Richards
equation requires a water retention model that ensures a zerowater content at oven dryness (Saito et al., 2006), as well as afilm component in the conductivity function, and an isothermal vapor conductivity component which is not included inour simulations.

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20326, Wiley Online Library on [01/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

4 of 10

Vadose Zone Journal

5 of 10F I G U R E 1 Water retention curves (left column) and hydraulicconductivity curves right column) for the three soils and for the van
Genuchten–Mualem model (VGM) and VGM model to the Brunswick
variant models. The green line shows the initial value of the pressure
head for each soil type and for all simulations. The red line shows thepressure head at the wilting (-16,000 cm or pF=4.2) point for bothRWU models.

F I G U R E 2 Cumulative actual evaporation for the three soil types(rows) and the two soil hydraulic property models (blue line: van
Genuchten-Mualem model [VGM], orange line: VGM model to theBrunswick variant). For the VGM model and for the Sandy soil, thesimulation could not be completed because once the soil reaches theresidual water content, no water can be removed by evaporation andnumerical problems occur (red X sign in the top figure).

3.3

from the cumulative potential transpiration line (black line in
Figure 3) always starts later than for VGM. This differenceis very small for the Clay soil since, until pFwp = 4.2, bothretention curves are similar (Figure 1, bottom left).
Similar to the Feddes model, the VGM-BW model yieldedhigher values of Cum. Transp. at 100 days (Figure 3, rightpanel) for the Nimah and Hanks RWU model. For the Sand,
both the VGM and VGM-BW simulations show a longer
period where actual transpiration is equal to potential, whichindicates an inherent compensation built in the model of
Nimah and Hanks, where the roots can apply higher suction
to satisfy the demand. However, after a critical time (different
for the two SHP models), actual transpiration is reduced due tothe lower conductivity of the dryer soil. We note that at the endof the simulation, the VGM model simulates a pressure headof –3100 cm, residual water content, and very low conductivity on the order of 10−9 cm day−1 at the soil surface (dryer
node due to higher root length density). On the contrary, for
the VGM-BW model, these values were -11,000 cm, 0.04 cm3
cm‒3 , and 10−6 cm day−1 , respectively. For the CLoam soil,

Simulations of root water uptake

For all three soils and the two RWU models, the simulations
with the VGM-BW model yielded a higher Cumulative Transpiration value (Cum. Transp.) at the end of the simulation(Figure 3, left panel). For the Feddes model, the differencein the final value of Cum. Transp. for the Sand and CLoam
soils was around 4 and 6 cm, respectively. This is due to theshape of the water retention curve in the dry range: from thewet range until the permanent wilting point (in our simulationℎwp = −16, 000 cm or pFwp = 4.2) and the concept of residual water content. For the VGM model, the retention curveshows a change in slope at around the residual water content, indicating less water capacity for RWU. This is especiallyvisible for the Sandy and CLoam soils. Due to this change
in slope in the dry range, the VGM simulation reaches thewilting point pressure head faster than the VGM-BW model
and, therefore, transpiration stops. This is not the case for the
VGM-BW model due to the linear part of the water retention
curve (Figure 1) in the dry range. For VGM-BW, deviation

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20326, Wiley Online Library on [01/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

DIAMANTOPOULOS ET AL.

DIAMANTOPOULOS ET AL.

Vadose Zone Journal

with the use of a two-step hydro-PTF functions (Weber et al.,2020), that allow the conversion of the VGM parameters to
the VGM-BW parameters. This guarantees more physicallycoherent simulations if SHP data are missing.

4

CONCLUSIONSWe implemented the Brunswick modular SHP model and thetwo-step hydro-pedotransfer function of Weber et al. (2020)in the Hydrus suite. The model allows for a more physicallycoherent representation of the water retention and hydraulicconductivity curves in the dry range. We also show that thewidely used VGM model may underestimate bare soil evaporation and transpiration when compared to the Brunswick
model. Given the simulation results of this study, we suggest replacing the VGM model with the VGM-BW model,
especially for simulations under arid conditions. Future work
will target the implementation of the bimodal version of the
VGM-BW model into the Hydrus suite.

F I G U R E 3 Cumulative actual transpiration for the three soil
types (rows), the two root water uptake models (left: Feddes, right:
Nimah & Hanks), and the two soil hydraulic property models (blueline: van Genuchten-Mualem model [VGM], orange line: VGM modelto the Brunswick variant).

due to the lower ℎinit , the drop in the hydraulic conductivity is
the limiting factor for reducing actual transpiration and overall lower uptake values in comparison with the Feddes model.
For the clay soil, differences between VGM and VGM-BW
are larger due to the larger differences in the unsaturated conductivity for the two models (Figure 1, right panel). We notethat the change of the slope in the conductivity data in thedry range can be a result of corner flow (Diamantopoulos &Durner, 2015) and/or film flow (Tuller & Or, 2001). Convincing experimental evidence of this change in slope exists forsandy soils, but not for heavier soil classes, meaning that our
results, especially for the CLoam and Clay soil and the Nimahand Hanks model, should be interpreted with caution, since nochange in the slope in the conductivity data can be observedor hypothesized.
Overall, there are significant differences between thesimulation results with the VGM and VGM-BW models concerning bare soil evaporation and RWU with the Hydrus suite.
This indicates the high need for water retention and unsaturated conductivity data of high quality, especially in the dry
range, and SHP models that can describe those data accurately (Peters, 2013; Weber et al., 2019). This can now bedone in Hydrus by the use of the Brunswick model as well as

AU T H O R C O N T R I B U T I O N S
Efstathios Diamantopoulos: Conceptualization; data curation; formal analysis; investigation; methodology; software; visualization; writing—original draft; editing. Jirka
Simunek: Software; writing—original draft; editing. Tobias
K. D. Weber: Writing—original draft; editing.
ORCID
Efstathios Diamantopoulos https://orcid.org/0000-00017870-0291
Jirka Simunek https://orcid.org/0000-0002-0166-6563
Tobias K. D. Weber https://orcid.org/0000-0002-34485208
REFERENCES
Assouline, S., & Or, D. (2013). Conceptual and parametric representation of soil hydraulic properties: A review. Vadose Zone Journal,
12(4), 1–20. https://doi.org/10.2136/vzj2013.07.0121
Blöcher, J. R., Diamantopoulos, E., Durner, W., & Iden, S. C. (2023).
Validating coupled flow theory for bare-soil evaporation under different boundary conditions. Vadose Zone Journal, 22(6), e20277.
https://doi.org/10.1002/vzj2.20277
Brooks, R. H., & Corey, A. T. (1964). Hydraulic properties of porous
media (Hydrology Papers No 3). Colorado State University.
Diamantopoulos, E., & Durner, W. (2012). Dynamic nonequilibrium ofwater flow in porous media: A review. Vadose Zone Journal, 11(3),0197. https://doi.org/10.2136/vzj2011.0197
Diamantopoulos, E., & Durner, W. (2015). Closed-form model forhydraulic properties based on angular pores with lognormal size distribution. Vadose Zone Journal, 14(2), 1–7. https://doi.org/10.2136/
vzj2014.07.0096
Du, E., Di Vittorio, A., & Collins, W. D. (2016). Evaluation of hydrologic
components of community land model 4 and bias identification. International Journal of Applied Earth Observation and Geoinformation,
48, 5–16.

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20326, Wiley Online Library on [01/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

6 of 10

Fayer, M. J., & Simmons, C. S. (1995). Modified soil water retentionfunctions for all matric suctions. Water Resources Research, 31(5),1233–1238.
Feddes, R. A., Kowalik, P. J., & Zaradny, H. (1978). Water Uptake byplants. In Simulation of field water use and crop yield (pp. 16–30).
Wiley.
Hohenbrink, T. L., Jackisch, C., Durner, W., Germer, K., Iden, S. C.,
Kreiselmeier, J., Leuther, F., Metzger, J. C., Naseri, M., & Peters, A.
(2023). Soil water retention and hydraulic conductivity measured in awide saturation range. Earth System Science Data Discussions, 2023,
1–20.
Iden, S. C., Blöcher, J. R., Diamantopoulos, E., & Durner, W. (2021).
Capillary, film, and vapor flow in transient bare soil evaporation (1):
Identifiability analysis of hydraulic conductivity in the medium to dry
moisture range. Water Resources Research, 57(5), e2020WR028513.
Iden, S. C., Diamantopoulos, E., & Durner, W. (2021). Capillary, film,
and vapor flow in transient bare soil evaporation (2): Experimental
identification of hydraulic conductivity in the medium to dry moisture
range. Water Resources Research, 57(5), e2020WR028514.
Kosugi, K. (1996). Lognormal distribution model for unsaturated
soil hydraulic properties. Water Resources Research, 32(9), 2697–
2703.
Lebeau, M., & Konrad, J.-M. (2010). A new capillary and thin film flow
model for predicting the hydraulic conductivity of unsaturated porous
media. Water Resources Research, 46(12), W12554.
Mualem, Y. (1976). A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research, 12(3),513–522.
Nimah, M., & Hanks, R. (1973). Model for estimating soil water, plant,
and atmospheric interrelations: II. Field test of model. Soil Science
Society of America Journal, 37(4), 528–532. https://doi.org/10.2136/
sssaj1973.03615995003700040019x
Peters, A. (2013). Simple consistent models for water retention andhydraulic conductivity in the complete moisture range. Water
Resources Research, 49(10), 6765–6780.
Peters, A., & Durner, W. (2008). Simplified evaporation methodfor determining soil hydraulic properties. Journal of Hydrology,
356(1-2), 147–162.
Richards, L. A. (1931). Capillary conduction of liquids through porous
mediums. Journal of Applied Physics, 1(5), 318–333.
Richardson, L. F. (1922). Weather prediction by numerical process.
Cambridge University Press.
Saito, H., Šimůnek, J., & Mohanty, B. P. (2006). Numerical analysis ofcoupled water, vapor, and heat transport in the vadose zone. Vadose
Zone Journal, 5(2), 784–800. https://doi.org/10.2136/vzj2006.0007
Šimůnek, J. (2006). Models of water flow and solute transport in theunsaturated zone, Encyclop. American Cancer Society.
Šimunek, J., Van Genuchten, M. Th., & Šejna, M. (2006). The HYDRUS
software package for simulating two-and three-dimensional movement of water, heat, and multiple solutes in variably-saturated media(Technical Manual Version 1.0). PC Progress.
Streck, T., & Weber, T. K. D. (2020). Analytical expressions for noncapillary soil water retention based on popular capillary retention
models. Vadose Zone Journal, 19(1), e20042. https://doi.org/10.1002/
vzj2.20042
Tokunaga, T. K. (2009). Hydraulic properties of adsorbed water films in
unsaturated porous media. Water Resources Research, 45(6), 460.
Tuller, M., & Or, D. (2001). Hydraulic conductivity of variably saturated porous media: Film and corner flow in angular pore space. Water
Resources Research, 37(5), 1257–1276.

Vadose Zone Journal

7 of 10van Genuchten, M. Th. (1980). Closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Science Society ofAmerica Journal, 44(5), 892–898. https://doi.org/10.2136/sssaj1980.
03615995004400050002x
Vogel, T., & Cislerova, M. (1988). On the reliability of unsaturated
hydraulic conductivity calculated from the moisture retention curve.
Transport in Porous Media, 3, 1–15.
Vogel, T., Gerke, H. H., Zhang, R., & Van Genuchten, M. Th.
(2000). Modeling flow and transport in a two-dimensional dualpermeability system with spatially variable hydraulic properties.
Journal of Hydrology, 238(1-2), 78–89.
Weber, T. K. D., Durner, W., Streck, T., & Diamantopoulos, E. (2019).
A modular framework for modeling unsaturated soil hydraulic properties over the full moisture range. Water Resources Research, 55(6),4994–5011.
Weber, T. K. D., Finkel, M., Conceição Gonçalves, M., Vereecken,
H., & Diamantopoulos, E. (2020). Pedotransfer function for theBrunswick soil hydraulic property model and comparison to thevan Genuchten–Mualem model. Water Resources Research, 56(9),e2019WR026820.
Weber, T. K. D., Weihermüller, L., Nemes, A., Bechtold, M., Degré,A., Diamantopoulos, E., Fatichi, S., Filipović, V., Gupta, S., &Hohenbrink, T. L. (2023). Hydro-pedotransfer functions: A roadmap
for future development. Egusphere. (preprint). https://doi.org/10.
5194/egusphere-2023-1860
Wesseling, J., Elbers, J., Kabat, P., & Van Den Broek, B. (1991).
Swatre: Instructions for input. internal note, Winand Staring Centre, Wageningen, The Netherlands. International Waterlogging andSalinity Research Institute, Lahore, Pakistan.

How to cite this article: Diamantopoulos, E.,
Simunek, J., & Weber, T. K. D. (2024).
Implementation of the Brunswick model system into
the Hydrus software suite. Vadose Zone Journal,
e20326. https://doi.org/10.1002/vzj2.20326

APPENDIX A: SOIL HYDRAULIC FUNCTIONS ANDANALYTICAL SOLUTIONSIn the following subsections, we first define the noncapillary saturation and conductivity functions, followed by
the description of the VGM, BC, and KS models, modified
VGM, in the original and BW variants. In each model case,the analytical solutions for the non-capillary saturation model
to Equation (A.1) are given. Therefore, we will use the equality 𝑥 = log10 (|ℎ|) as in the original publication (Streck andWeber, 2020).
The BW non-capillary saturation and conductivity functions
To model the non-capillary saturation 𝑆nc as a function ofpressure head ℎ, (Weber et al., 2019) developed a universal
and flexible integral approach. The main advantage is thatit can, in principle, be used for any model of Φe (ℎ). In the
A.1

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20326, Wiley Online Library on [01/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

DIAMANTOPOULOS ET AL.

DIAMANTOPOULOS ET AL.

Vadose Zone Journal

following, we adopt the notation by (Streck and Weber, 2020).
For convenience, we first introduce 𝑥 = log10 (|ℎ|). Then we
write 𝑆nc as a function of 𝑥 which is given by
𝑥 (

∗pF

𝑆nc (𝑥) =

∫−∞

) ′
pF ′
Φe (𝑥 ) − 1 𝑑𝑥

(A.1)

′

pF

where Φe (𝑥) [-] is Φ𝑒 (ℎ) in log10 , and 𝑥 is the dummy
variable of integration. Streck and Weber (2020) explicitly
explain that the integral equals the area between full saturation′(𝑥 = −∞, corresponding to ℎ = 0 cm) and the capillary saturation at a given pF value. Analytic solutions for special casesof 𝑆c (ℎ) were derived by Streck and Weber (2020). Numericalintegration from −∞ is not tractable. Thus, numerical integration is done by replacing −∞ with −3, which corresponds toa pressure head of 10−3 cm. Weber et al. (2019) transform
Equation (A.1) to normal space by substituting 𝑥 by ℎ. Using
the relationship1𝑑𝑥= − log10 (𝑒) ,
𝑑ℎ
ℎ

(A.2)

then yields
∗
𝑆nc
(ℎ) =
∫

−10𝜖

(

′

′

Φe (ℎ ) − 1

′

𝑑ℎ .

(A.3)

∗ (ℎ)
𝑆nc
∗ (ℎ )
𝑆nc
0

(A.4)

𝐾rnc [-] is a function of the non-capillary saturation (Peters,
2013; Tokunaga, 2009)(
𝐾rnc (ℎ) =

|ℎ0 |
ℎr

pF
𝑆nc
(𝑥) =

[
(
1
1 )
B(1 + (𝛼10𝑥 )𝑛 ; 1, 0) − B 1 + (𝛼10𝑥 )𝑛 ; , 0𝑛 ln(10)𝑛+𝑛 −where B is the incomplete beta function.
The unsaturated hydraulic conductivity of the capillary part
(Equation 6) for the VGM variant of the Brunswick model
framework is given by𝐾rc (ℎ) = 𝑆c (ℎ)

𝜆A.2
The VGM and VGM-BW models
According to van Genuchten (1980), the effective water saturation, Φe is given as a function of pressure head by:(A.6)

where 𝛼 [L−1 ] and 𝑛 [-] are shape parameters and 𝑚 = 1 −1∕𝑛. The relative unsaturated hydraulic conductivity is given[

1−1 − Φe (ℎ)1∕𝑚]𝑚 )2

1 − Φe (ℎ0 )1∕𝑚

(A.9)

Equations (VGM-BW): 4, 5, A.6, A.8, 6, A.9, A.5
Parameters (VGM-BW): 𝜃sc , 𝜃snc , 𝛼, 𝑛, 𝐾sc , 𝐾snc , 𝜆, 𝛼s , ℎ0
A.3
The BC and BC-BW models
The Brooks–Corey (BC) (Brooks & Corey, 1964) effective
saturation model is given as a function of pressure head byΦe (ℎ) ={|𝛼ℎ|−𝑛ℎ < −1∕𝛼1ℎ ≥ −1∕𝛼,(A.10)

with 𝑛 [-] and 𝛼 [L-1 ] as flexible parameters.
The unsaturated hydraulic conductivity function is given by

where ℎr [L] ensures matching dimensions in 𝐾rnc (ℎ) and is
set to 1 cm. 𝛼s controls the slope of decrease in 𝐾rnc (ℎ).

Φ𝑒 (ℎ) = [1 + [|𝛼ℎ|]𝑛 ]−𝑚 ,

]
∞
(
)∑
1
1
1
,
+
−1
𝑛
𝑛
𝑘(𝑘 + 1)(𝑛𝑘 + 1)
𝑘=1

(A.8)

)−𝛼s (1−𝑆nc (ℎ)),
(A.5)

(A.7)

where 𝜆 [-] is a shape parameter of the conductivity
function.
Equations (VGM): 2, 3, A.6, A.7
Parameters (VGM): 𝜃s , 𝜃r , 𝛼, 𝑛, 𝐾s , 𝜆According to Streck and Weber (2020), the analytical solution of the non-capillary saturation function for the van
Genuchten case is given as(∗ (ℎ) neither reaches from 1 at full saturation to 0 at
Since 𝑆nc
complete dryness, shifting and re-scaling are required, whichis achieved by

𝑆nc (ℎ) = 1 −𝐾r (ℎ) = Φe (ℎ)𝜆 [1 − (1 − Φe 1∕𝑚 )𝑚 ]2 ,

)

ℎ′

ℎ

by:

𝐾𝑟 (ℎ) = Φe (ℎ)2∕𝑛+𝜆+2 .

(A.11)

Equations (BC): 2, 3, A.10, A.11
Parameters (BC): 𝜃s , 𝜃r , 𝛼, 𝑛, 𝐾s , 𝜆By substituting 𝑛 = 1∕𝛽 the solution to Equation (A.1) for
the BC-BW model is given as⎧ 𝛽+ln (1∕𝛼)⎪pF
𝑆nc (𝑥) = ⎨⎪0⎩ln(10)−1𝛽− ln(10)(𝛼10𝑥 ) 𝛽 − 𝑥

𝛼10𝑥 ≥ 1otherwise.
(A.12)

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20326, Wiley Online Library on [01/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

8 of 10

Vadose Zone Journal

The unsaturated hydraulic conductivity function for thecapillary part is given by:𝐾rc (ℎ) = Φe (ℎ)2∕𝑛+𝜆+2 .

(A.13)

Equations (BC-BW): 4, 5, A.10, A.12, 6, A.13, A.5
Parameters (BC-BW): 𝜃𝑠𝑐 , 𝜃𝑠𝑛𝑐 , 𝛼, 𝑛, 𝐾𝑠𝑐 , 𝐾𝑠𝑛𝑐 , 𝜆, 𝛼𝑠 , ℎ0
A.4
The Kosugi and Kosugi-BW models
The Kosugi (KS) (Kosugi, 1996a) effective saturation model
is given as a function of pressure head by⎧⎪0.5 erf c
Φe (ℎ) = ⎨⎪1⎩(ln(ℎ∕𝛼)√2𝑛(𝐾r (ℎ) = Φe (ℎ)ℎ<0𝜃(ℎ) =(){
𝜃𝛼 + 𝜃𝑚 − 𝜃𝛼 Φe (ℎ) ℎ < ℎs𝜃sℎ ≥ ℎs ,{(A.14)

Φe (ℎ) =otherwise.

])2[ln (ℎ∕𝛼)
𝑛
1
.
erf c
+√
√
2
2𝑛
2

(A.15)

Equations (KS): 2, 3, A.14, A.15
Parameters (KS): 𝜃s , 𝜃r , 𝛼, 𝑛, 𝐾s , 𝜆Note that in Hydrus, the symbol 𝛼 instead of ℎ0 and 𝑛
instead of 𝜎 are used, compared with the original notation in
Kosugi (1996a).
Substituting 𝑥𝑚 = log10 𝛼 the analytical of the integral
model using Equation (A.14) is
(
)
⎧ 𝑥−𝑥
ln(10)(𝑥−𝑥𝑚 )
𝑛
𝑚
√
√
−
⎪ 2 erf cln(10)
2𝜋
2
𝑛
⎪
pF
)
(
2
2
𝑆nc (𝑥) = ⎨ln (10)(𝑥−𝑥𝑚 )− (𝑥 − 𝑥𝑚 )ℎ<0
2 𝑛2
⎪ exp −
⎪0otherwise,⎩(A.16)and the hydraulic conductivity curve of the capillary poredomain is given as(
𝐾rc (ℎ) = 𝑆c (ℎ)

A.5
The modified VGM and VGM-BW models
A more flexible model of SHPs (modVGM), proposed from
Vogel and Cislerova (1988), is implemented in Hydrus. Vogel
and Cislerova (1988) modified the equation of van Genuchten(1980) to add flexibility in the near saturation SHPs. The soil
water retention in Equation (2) is given in an alternative form
by:

(A.18)

where Φe (ℎ) is given by:)

The relative hydraulic conductivity is given by
𝜆

9 of 10𝜆1(1+(𝛼∣ℎ∣)𝑛 )𝑚ℎ < ℎs1.0ℎ ≥ ℎs ,and 𝜃𝑚 [L3 L−3 ] is a fictitious parameter slightly higher than𝜃s , and 𝜃𝛼 [L3 L−3 ] is an extra parameter, with 𝜃𝛼 ≤ 𝜃r , that
adds flexibility in the dry range of the water retention.
The relative conductivity function for the modified van
Genuchten model is given by:
1⎧𝐾r (ℎ)⎪𝐾(ℎ−ℎ )(𝐾s −𝐾𝑘 )𝐾r (ℎ) = ⎨ 𝐾𝑘 + 𝐾 𝑘(ℎ −ℎ
s
s s
𝑘)
⎪1.0
⎩

ℎ ≤ ℎ𝑘ℎ𝑘 < ℎ ≤ ℎs

(A.20)ℎ ≥ ℎs ,where𝐾𝐾r1 (ℎ) = 𝑘
𝐾s

[

] [
]
Φe 𝜆 𝐹 (𝜃r ) − 𝐹 (𝜃) 2,
Φe𝑘
𝐹 (𝜃r ) − 𝐹 (𝜃𝑘 )(A.21)and[(𝐹 (𝜃) = 1 −𝜃 − 𝜃𝛼𝜃𝑚 − 𝜃𝛼

)1∕𝑚 ]𝑚√√ )2
erf [erf c−1 [2Φ𝑒 (ℎ0 )] + 𝑛∕ 2] − erf [erf c−1 [2Φ𝑒 (ℎ)] + 𝑛∕ 2].
√
1 + erf [erf c−1 [2Φe (ℎ0 )]] + 𝑛∕ 2Note that the rescaled 𝑆c (ℎ) is used as multiplier before thebracket and Kosugi’s original Φe (ℎ) from Equation (A.15) in
the fraction which is evaluated in two instances at Φe (ℎ0 =106.8 ).
Equations (KS-BW): 4, 5, A.14, A.16, 6, A.17, A.5
Parameters (KS-BW): 𝜃sc , 𝜃snc , 𝛼, 𝑛, 𝐾sc , 𝐾snc , 𝜆, 𝛼s , ℎ0

(A.19)

Φe𝑘 =𝜃𝑘 − 𝜃r.
𝜃s − 𝜃r(A.22)(A.17)(A.23)The modified VGM model allows for a non-zero minimum
capillary height ℎs [L], by adjusting the value of 𝜃𝑚 , to a valueslightly higher than 𝜃s . In that case, 𝜃 = 𝜃s (and 𝐾(ℎ) = 𝐾s )15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20326, Wiley Online Library on [01/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

DIAMANTOPOULOS ET AL.

DIAMANTOPOULOS ET AL.

Vadose Zone Journal

for ℎ ≥ ℎs . This change has no effect on the water retention
curve, but a big effect on the conductivity function, especiallyfor fine-textured soils, with 1.0 < 𝑛 ≤ 1.3 (Šimůnek, 2006).
A second conductivity point (less or equal to 𝐾s ) can be provided by the user at which 𝐾𝑘 = 𝐾(𝜃𝑘 (ℎ𝑘 )) [L T−1 ], where
the conductivity value is linearly interpolated between the twopoints ((ℎ𝑠 , 𝐾𝑠 ) and (ℎ𝑘 , 𝐾𝑘 )).
A simpler version of this model is also implemented in
Hydrus, where users can fix ℎs = −2 cm, as proposed by(Vogel and Cislerova, 1988), and recommended for finetextured soils. In that model ℎ𝑘 = ℎs , 𝜃𝑘 = 𝜃s , and 𝐾𝑘 =𝐾s .
Equations (modVGM): A.18, A.19, 3, A.20, A.21, A.22,
A.23
Parameters (modVGM): 𝜃s , 𝜃𝑚 , 𝜃r , 𝜃𝛼 , 𝛼, 𝑛, 𝐾s , 𝐾𝑘 , 𝜃𝑘 , 𝜆The Brunswick variant of the modified VGM model(modVGM-BW) is given by:
𝜃(ℎ) = 𝜃𝑚 𝑆c (ℎ) + 𝜃snc 𝑆nc (ℎ),(A.24)

where 𝜃𝑚 is higher than 𝜃sc and 𝑆c (ℎ) is calculated from
Equation (5) where Φe is given by Equation (A.19).
The relative conductivity function for the capillary part in
modVGM-BW is calculated by:1⎧𝐾rc (ℎ)⎪ 𝐾𝑘(ℎ−ℎ )(𝐾sc −𝐾𝑘 )𝐾rc (ℎ) = ⎨ 𝐾 + 𝐾 𝑘(ℎ −ℎ
sc
sc s
𝑘)
⎪1.0
⎩

ℎ ≤ ℎ𝑘ℎ𝑘 < ℎ ≤ ℎs

(A.25)ℎ ≥ ℎs ,1 (ℎ) is calculated by Equations (A.21–A.23),
where 𝐾rc
where 𝜃r = 𝜃𝛼 = 0. The simpler version with ℎs = −2 cm wasalso implemented. Note that for 𝑆nc the function derived for
the VGM model can be used directly since Φe is identical in
the VGM-BW and modVGM-BW case.
Equations (modVGM-BW): A.24, 5, A.19, A.8, 6, A.25,
A.21, A.22, A.23
Parameters (modVGM-BW): 𝜃sc , 𝜃𝑚 , 𝜃snc , 𝛼, 𝑛, 𝐾s , 𝐾𝑘 , 𝜃𝑘 ,𝜆, 𝑎s , ℎ015391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20326, Wiley Online Library on [01/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

10 of 10