| Title: | Simulation of Soil Water Flow Based on Darcy's Law |
|---|---|
| Description: | Provides a simple model to simulate soil water content and water flux changes at multiple soil depths. Calculations mainly follow Allen et al. (1998, ISBN:92-5-104219-5) and Teh (2017, ISBN:9789673447626). |
| Authors: | Minguo Liu [aut, cre] |
| Maintainer: | Minguo Liu <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.1.1 |
| Built: | 2026-10-24 18:25:45 UTC |
| Source: | https://github.com/cran/soilwatR |
In the absence of RH max and RHmin, it can be used to estimate actual vapour pressure.
cal_ActualVapourPressure_from_RHmean(RHmean, Tmax, Tmin)cal_ActualVapourPressure_from_RHmean(RHmean, Tmax, Tmin)
RHmean |
mean relative humidity(%). |
Tmax |
daily maximum temperature (degrees Celsius). |
Tmin |
daily minimum temperature (degrees Celsius). |
A vector of actual vapour pressure
FAO Irrigation and drainage paper 56 (P39)
The atmospheric pressure, P, is the pressure exerted by the weight of the earth's atmosphere.
cal_atmosphericPressure(elevation)cal_atmosphericPressure(elevation)
elevation |
elevation above sea level (m) |
Assuming 20°C for a standard atmosphere. Evaporation at high altitudes is promoted due to low atmospheric pressure as expressed in thepsychrometric constant. The effect is, however, small and in the calculation procedures, the average value for a location is sufficient.
A vector of atmospheric pressure (Kpa)
FAO Irrigation and drainage paper 56 (P31)
Based on lat, z, J, Tmax, Tmin, n, RHmax, RHmin, windSpeed parameters, reference evapotranspiration was calculated by Penman-Monteith.
cal_ET0(Latitude, Altitude, J, Tmax, Tmin, Rs, RHmean, Wind)cal_ET0(Latitude, Altitude, J, Tmax, Tmin, Rs, RHmean, Wind)
Latitude |
latitude (radian), positive for the northern hemisphere and negative for the southern hemisphere. |
Altitude |
station elevation above sea level [m]. |
J |
is the number of the day in the year between 1 (1 January) and 365 or 366 (31 December). |
Tmax |
daily maximum air temperature (degrees Celsius). |
Tmin |
daily minimum air temperature (degrees Celsius). |
Rs |
Solar radiation [MJ m-2 d-1]. |
RHmean |
daily mean relative humidity %. |
Wind |
wind speed at 2 m height [m s-1]. |
A vector of reference evapotranspiration (mm/day)
The FAO Penman-Monteith method is maintained as the sole standard method for the computation of ETo from meteorological data.
cal_ET0_from_PM(delta, Rn, G, gamma, Tem, u2, es, ea)cal_ET0_from_PM(delta, Rn, G, gamma, Tem, u2, es, ea)
delta |
slope vapour pressure curve (kPa °C). From cal_slopeOfSaturationVapourPressureCurve() |
Rn |
net Radiation at the crop surface [MJ m-2 day-1]. From cal_netRadiation() |
G |
soil heat flux density [MJ m-2 day-1]. |
gamma |
psychrometric constant (kPa °C). |
Tem |
air temperature at 2 m height [°C]. |
u2 |
wind speed at 2 m height [m s-1]. |
es |
saturation vapour pressure [kPa]. |
ea |
actual vapour pressure [kPa]. |
A vector of reference evapotranspiration [mm day-1].
Ten-day or monthly time step :
Notwithstanding the non-linearity in the Penman-Monteith equation and some weather parameter methods, mean ten-day or monthly weather data can be used to compute the mean ten-day or monthly values for the reference evapotranspiration. The value of the reference evapotranspiration calculated with mean monthly weather data is indeed very similar to the average of the daily ETo values calculated with daily average weather data for that month.
When the soil is warming (spring) or cooling (autumn), the soil heat flux (G) for monthly periods may become significant relative to the mean monthly Rn. In these cases G cannot be ignored and its value should be determined from the mean monthly air temperatures of the previous and next month.
Daily time step:
Calculation of ETo with the Penman-Monteith equation on 24-hour time scales will generally provide accurate results.
As the magnitude of daily soil heat flux (G) beneath the reference grass surface is relatively small, it may be ignored for 24-hour time steps.
FAO Irrigation and drainage paper 56 (P65)
The extraterrestrial radiation, Ra, for each day of the year and for different latitudes can be estimated from the solar constant, the solar declination and the time of the year.
cal_extraterrestrialRadiation_for_daily(J, lat)cal_extraterrestrialRadiation_for_daily(J, lat)
J |
is the number of the day in the year between 1 (1 January) and 365 or 366 (31 December). |
lat |
latitude (radian), positive for the northern hemisphere and negative for the southern hemisphere. |
A vector of extraterrestrial radiation for daily(MJ m-2 day-1)
FAO Irrigation and drainage paper 56 (P45-47)
calculating field capacity
cal_FC(Sand_frac, Clay_frac, OM)cal_FC(Sand_frac, Clay_frac, OM)
Sand_frac |
Sand fraction |
Clay_frac |
Clay fraction |
OM |
Organic matter content percentage |
A vector of field capacity
Calculating inverse relative distance Earth-sun
cal_inverseRelativeDistance_Earth_sun(J)cal_inverseRelativeDistance_Earth_sun(J)
J |
is the number of the day in the year between 1 (1 January) and 365 or 366 (31 December). |
A vector of inverse relative distance Earth-sun (radian)
FAO Irrigation and drainage paper 56 (P46)
Typical values for the crop coefficient at the end of the late season growth stage, Kc end, are listed in Table 12 for various agricultural crops.
cal_Kcend_for_singleKc(RHmine, u2e, Ktable, he)cal_Kcend_for_singleKc(RHmine, u2e, Ktable, he)
RHmine |
mean value for daily minimum relative humidity during the mid-season growth stage , for 20 <= RHmine<= 80 |
u2e |
mean value for daily wind speed at 2 m height over grass during the mid- season growth stage (m/s), for 1 <= u2e <= 6 |
Ktable |
value for Kc mid taken from Table 12 |
he |
mean plant height during the mid-season stage (m) for 0.1 m < h < 10 m |
A vector for the K value of the late season stage
only applied when the tabulated values for Kc end exceed 0.45. The equation reduces the Kc end with increasing RHmin. This reduction in Kc end is characteristic of crops that are harvested 'green' or before becoming completely dead and dry (Kc end >= 0.45).
No adjustment is made when Kc end (Table) < 0.45 (Kc end = Kc end (Tab)). When crops are allowed to senesce and dry in the field (as evidenced by Kc end < 0.45), u2 and RHmin have less effect on Kc end and no adjustment is necessary.
FAO Irrigation and drainage paper 56 (P125-126)
Calculating Kcini value
cal_Kcini_for_SingleKc(Pmean, ET0, tw, type, fw)cal_Kcini_for_SingleKc(Pmean, ET0, tw, type, fw)
Pmean |
is the average depth od infiltrated water per wetting events(mm) |
ET0 |
mean ET0 during initial period(mm/day) |
tw |
is the mean interval between wetting events(days) |
type |
soil type:"coarse soil textures" and "medium and fine soil textures" |
fw |
the fraction of surfaces wetted by irrigation or rain (0-1) |
A vector for the K value of initial stage
FAO Irrigation and drainage paper 56 (P263-269)
Typical values for the crop coefficient at the end of the late season growth stage, Kc end, are listed in Table 12 for various agricultural crops.For specific adjustment in climates where RHmin differs from 45 % or where u2 is larger or smaller than 2.0 m/s.
cal_Kcmid_for_singleKc(RHmine, u2e, Ktable, he)cal_Kcmid_for_singleKc(RHmine, u2e, Ktable, he)
RHmine |
mean value for daily minimum relative humidity during the mid-season growth stage , for 20 <= RHmine <= 80 |
u2e |
mean value for daily wind speed at 2 m height over grass during the mid- season growth stage (mls), for 1 <= u2e <= 6 |
Ktable |
value for Kc mid taken from Table 12 |
he |
mean plant height during the mid-season stage [m] for 0.1 m < h < 10 m |
A vector for the K value of the mid-season stage
FAO Irrigation and drainage paper 56 (P122-123)
Due to the non-linearity of the above equation, the mean saturation vapour pressure for a day, week, decade or month should be computed as the mean between the saturation vapour pressure at the mean daily maximum and minimum air temperatures for that period.
cal_meanSaturationVapourPressure(Tmax, Tmin)cal_meanSaturationVapourPressure(Tmax, Tmin)
Tmax |
the daily maximum air temperature(degrees Celsius). |
Tmin |
the daily minimum air temperature(degrees Celsius). |
Using mean air temperature instead of daily minimum and maximum temperatures results in lower estimates for the mean saturation vapour pressure. The corresponding vapour pressure deficit (a parameter expressing the evaporating power of the atmosphere) will also be smaller and the result will be some underestimation of the reference crop evapotranspiration. Therefore,the mean sauration vapour pressure should be calculated as the mean between the saturation vapour pressure at both the daily maximum and minimum air temperature.
A vector of mean saturation vapour pressure (es)
FAO Irrigation and drainage paper 56 (P36)
Calculating net longwave radiation Rnl
cal_netLongwaveRadiation(TKmax, TKmin, ea, Rs, Rso)cal_netLongwaveRadiation(TKmax, TKmin, ea, Rs, Rso)
TKmax |
maximum absolute temperature during the 24-hour period [K]. |
TKmin |
minimum absolute temperature during the 24-hour period [K]. |
ea |
actual vapour pressure [kPa]. |
Rs |
measured or calculated solar radiation [MJ m-2 day-1]. From cal_solarRadiation(). |
Rso |
calculated clear-sky radiation [MJ m-2 day-1]. From cal_skySolarRadiation_withas_bs() or cal_skySolarRadiation_withas_elevation(). |
A vector of net outgoing longwave radiation [MJ m-2 day-1]
The rate of longwave energy emission is proportional to the absolute temperature of the surface raised to the fourth power. This relation is expressed quantitatively by the Stefan-Boltzmann law. The net energy flux leaving the earth's surface is, however, less than that emitted and given by the Stefan-Boltzmann law due to the absorption and downward radiation from the sky. Water vapour, clouds, carbon dioxide and dust are absorbers and emitters of longwave radiation. Their concentrations should be known when assessing the net outgoing flux. As humidity and cloudiness play an important role, the Stefan-Boltzmann law is corrected by these two factors when estimating the net outgoing flux of longwave radiation. It is thereby assumed that the concentrations of the other absorbers are constant. An average of the maximum air temperature to the fourth power and the minimum air temperature to the fourth power is commonly used in the Stefan-Boltzmann equation for 24-hour time steps. The term (0.34-0.14*sqrt(ea)) expresses the correction for air humidity, and will be smaller if the humidity increases. The effect of cloudiness is expressed by (1.35 Rs/Rso - 0.35). The term becomes smaller if the cloudiness increases and hence Rs decreases. The smaller the correction terms, the smaller the net outgoing flux of longwave radiation. Note that the Rs/Rso term in Equation 39 must be limited so that Rs/Rso <= 1.0. Where measurements of incoming and outgoing short and longwave radiation during bright sunny and overcast hours are available, calibration of the coefficients in Equation 39 can be carried out.
FAO Irrigation and drainage paper 56 (P51-52)
The net radiation (Rn) is the difference between the incoming net shortwave radiation (Rns) and the outgoing net longwave radiation (Rnl).
cal_netRadiation(Rns, Rnl)cal_netRadiation(Rns, Rnl)
Rns |
incoming net shortwave radiation. From cal_netSolarRadiation(). |
Rnl |
outgoing net longwave radiation. From cal_netLongwaveRadiation(). |
A vector of net radiation
calculating the net rain
cal_netRain(lai, rain)cal_netRain(lai, rain)
lai |
leaf area index m2 m-2 |
rain |
rain mm |
A vector of net rain (mm)
The net shortwave radiation resulting from the balance between incoming and reflected solar radiation.
cal_netSolarRadiation(alpha, Rs)cal_netSolarRadiation(alpha, Rs)
alpha |
albedo or canopy reflection coefficient, which is 0.23 for the hypothetical grass reference crop [dimensionless]. |
Rs |
the incoming solar radiation [MJ m-2 day-1]. From cal_solarRadiation() |
A vector of net solar or shortwave radiation [MJ m-2 day-1].
FAO Irrigation and drainage paper 56 (P51)
This function calculates the potential evapotranspiration according to the single crop coefficient method, and then converts the potential evapotranspiration into potential evaporation and potential transpiration based on LAI
cal_PET(data, Kc_mid, Kc_end, fw, soil_type)cal_PET(data, Kc_mid, Kc_end, fw, soil_type)
data |
A data frame. |
Kc_mid |
value for Kc mid |
Kc_end |
value for Kc end |
fw |
the fraction of surfaces wetted by irrigation or rain (0-1) |
soil_type |
soil type:"coarse soil textures" and "medium and fine soil textures" |
A data frame including the potential evaporation and potential transpiration
The stages of data should include all four stages. If a crop has multiple growth cycles, each cycle should include all four stages.
Calculating psychrometric constant
cal_psychrometriCconstant(atmospheric_pressure)cal_psychrometriCconstant(atmospheric_pressure)
atmospheric_pressure |
atmospheric pressure (kPa). |
A vector of psychrometric constant (kPa/degree Celsius)
FAO Irrigation and drainage paper 56 (P32)
calculating reduction rate for Potential transpiration
cal_Red_PETc(C3, vwc_root, pwp_root, sat_root)cal_Red_PETc(C3, vwc_root, pwp_root, sat_root)
C3 |
Judge whether it is a C3 plant. If not, it is C4 |
vwc_root |
Current volume moisture content of root zone |
pwp_root |
Volume moisture content of wilting point in root zone |
sat_root |
Saturated volume moisture content of root zone |
A vector for the reduction rate of transpiration
calculating reduction rate for Potential evaporation
cal_Red_PETs(vwc1, sat1)cal_Red_PETs(vwc1, sat1)
vwc1 |
Volumetric moisture content of topsoil fraction |
sat1 |
Saturated volume moisture content of topsoil fraction |
A vector for the reduction rate of evaporation
calculating Soil saturated water content
cal_Sat(Sand_frac, Clay_frac, OM, THETA_33)cal_Sat(Sand_frac, Clay_frac, OM, THETA_33)
Sand_frac |
Sand fraction |
Clay_frac |
Clay fraction |
OM |
Organic matter content percentage |
THETA_33 |
wilting point from cal_FC function |
A vector of Soil saturated water content
saturation vapour pressure at the air temperature T.
cal_saturationVapourPressure(Tem)cal_saturationVapourPressure(Tem)
Tem |
air temperature (degrees Celsius). |
A vector of saturation vapour pressure at the air temperature T (kPa).
FAO Irrigation and drainage paper 56 (P36)
The calculation of the clear-sky radiation, Rso, when n = N, is required for computing net longwave radiation.
cal_skySolarRadiation_withas_elevation(z, Ra)cal_skySolarRadiation_withas_elevation(z, Ra)
z |
station elevation above sea level [m]. |
Ra |
extraterrestrial radiation [MJ m-2 day-1]. From cal_extraterrestrialRadiation_for_daily() |
A vector of clear-sky solar radiation [MJ m-2 day-1].
FAO Irrigation and drainage paper 56 (P51)
Calculating slope of saturation vapour pressure curve
cal_slopeOfSaturationVapourPressureCurve(Tem)cal_slopeOfSaturationVapourPressureCurve(Tem)
Tem |
is air temperature (degrees Celsius). |
In the FAO Penman-Monteith equation, where it occurs in the numerator and denominator, the slope of the vapour pressure curve is calculated using mean air temperature.
A vector for slope of saturation vapour pressure curve at air temperature T
FAO Irrigation and drainage paper 56 (P3)
Calculating solar declination with FAO56 method
cal_solarDeclination_in_FAO(J)cal_solarDeclination_in_FAO(J)
J |
is the number of the day in the year between 1 (1 January) and 365 or 366 (31 December) |
A vector of solar declination (radian)
FAO Irrigation and drainage paper 56 (P46)
Calculating sunset hour angle
cal_sunsetHourAngle(lat, solar_declination)cal_sunsetHourAngle(lat, solar_declination)
lat |
latitude (radian), positive for the northern hemisphere and negative for the southern hemisphere. |
solar_declination |
solar declination(radina). |
A vector of sunset Hour Angle(radian)
FAO Irrigation and drainage paper 56 (P46)
calculating wilting point
cal_WP(Sand_frac, Clay_frac, OM)cal_WP(Sand_frac, Clay_frac, OM)
Sand_frac |
Sand fraction |
Clay_frac |
Clay fraction |
OM |
Organic matter content percentage |
A vector of wilting point
Converting the unit of angle in longitude and latitude into the unit of radian.
convert_angert_to_radian(anger)convert_angert_to_radian(anger)
anger |
Longitude or latitude in angle. |
A vector of longitude or latitude in radian
Convert degrees Celsius to Fahrenheit
convert_degreesCelsius_to_Fahrenheit(degrees_Celsius)convert_degreesCelsius_to_Fahrenheit(degrees_Celsius)
degrees_Celsius |
temperature in degrees Celsius(°C). |
A vector of temperature in Fahrenheit(°F)
FAO56
Linear interpolation is performed by using the values on both sides of the missing values.
linear_interpolation(DataVector)linear_interpolation(DataVector)
DataVector |
data vector.Note that the starting value of vector needs to be no missing value. |
A vector
Graphically display the output results of the model
plot_result(Result)plot_result(Result)
Result |
The output result of the soilwatr function |
A combined plot if ggpubr is installed; otherwise a list of ggplot objects.
teh data including Input data, soilProfileSet data and obs data. Input data is a form of data frame includes daily upper boundary layer data of two irrigation treatments. The soilProfileSet is a matrix data that stores the soil profile setting information. Obs is the measured change of soil water content at different soil depths.
SIalfalfaSIalfalfa
A list with 3 data.
Simulation of soil water flow in different soil layers
soilwatr(data, param, soilProfileSet, verbose = FALSE)soilwatr(data, param, soilProfileSet, verbose = FALSE)
data |
A data frame. The model input data, including Day, rain_mm, lai, T_mm, E_mm. |
param |
A list. Input parameters of the model. |
soilProfileSet |
A matrix. Definition of soil profile and its physical characteristics. |
verbose |
Logical; if |
A matrix of simulated results.
data(SIalfalfa) param <- list( initialRootingDepth = 1.2, RootGrowthRate = 0.0035, has_watertable = FALSE, subintervals = 20, maxsteps = 5 ) SIalfalfa_T1 <- SIalfalfa$Input[SIalfalfa$Input$Tre2 == "T1", ] soilProfileSet_T1 <- SIalfalfa$soilProfileSet soilProfileSet_T1["THETA_v_curr", ] <- c(0.12, 0.13, 0.28, 0.31, NA) SI_T1_result <- soilwatr( data = SIalfalfa_T1, param = param, soilProfileSet = soilProfileSet_T1 ) head(SI_T1_result)data(SIalfalfa) param <- list( initialRootingDepth = 1.2, RootGrowthRate = 0.0035, has_watertable = FALSE, subintervals = 20, maxsteps = 5 ) SIalfalfa_T1 <- SIalfalfa$Input[SIalfalfa$Input$Tre2 == "T1", ] soilProfileSet_T1 <- SIalfalfa$soilProfileSet soilProfileSet_T1["THETA_v_curr", ] <- c(0.12, 0.13, 0.28, 0.31, NA) SI_T1_result <- soilwatr( data = SIalfalfa_T1, param = param, soilProfileSet = soilProfileSet_T1 ) head(SI_T1_result)