Variably Saturated flow

Numerous numerical models have been developed for simulating water movements in materials with various degrees of saturation. Clement et al (1994) reviewed many of them underlining the most common shortcomings. In most applications, the pressure-based form of the variably saturated flow equation through homogeneous and isotropic porous media is used (Cooley, 1983; Huyakorn et al., 1984):
                                                                                        (1)
where:
S_s is the specific storage of the medium [L^-1]
Ψ is the pressure head [L]
θ is the water content
η is the porosity
K(Ψ) is the hydraulic conductivity [L T^-1]
t is time [T]
This equation describes the movement of water in soils under the following assumptions:
1. The dynamics of the air phase do not affect those of the water phase;
2. The density of water is only a function of pressure;
3. The spatial gradient of the water density is negligible.
Finite elements solutions to this approach often incur mass-balance problems in unsaturated media (Celia et al., 1990).
The finite-difference algorithm developed by Clement et al. has the advantages of being computationally simple while remaining capable of modelling a wide variety of problems, including infiltration into dry soils. In this approach, a mixed form of Equation 1 is solved by calculating the pressure heads Ψ at each node for a given time step (n+1) using a modified Picard iteration, where at each iteration (m+1) we solve a system of linear equations defined as:
                                                                                                                                                                                                                                   (2)
where
 
 
 
                                                                                                              (3)
where
Θ is the moisture content
K is the hydraulic conductivity
C is the water capacity
Δx is the horizontal step
Δz is the vertical step
Δt is the time step
S_s is the specific storage
η is the porosity
As can be seen, all the coefficients used for the current iteration (m+1) are calculated from the results of the previous iteration (m). Θ, K and C vary for each node and are functions of Ψ, therefore need to be re-calculated after each iteration.
As discussed earlier, the link between Θ and Ψ is given by the Soil-Water Characteristic Curve (SWCC). The hydraulic conductivity of the materials can be expressed as a function of the pressure heads (Van Genuchten, 1980) by:
                                                                                               (4)
Finally, by definition, the water capacity represents the slope of the SWCC and therefore can be expressed as:
                                                                                                                                               (5)
Equation 2 applies to all interior nodes, while at boundary nodes it needs to be modified to consider the appropriate boundary conditions.
In the examples shown hereafter we consider a pavement structure like the one in Figure 1, where the water table depth is 2m from the top of the unbound layer. In normal conditions this system is in equilibrium and the water content does not vary with time.
Figure 1: Variably saturated flow – pavement structure
The boundaries considered are as follows (see also Figure 2):
► Top and right-hand side nodes: Neumann boundaries (no flow);
► Bottom nodes: Dirichlet boundaries (constant moisture content);
► Left-hand side nodes: symmetry.
Figure 2: Variably saturated flow – boundary conditions

Read More

fenv calculation 2

The inputs are organised in three hierarchical levels of detail, level 1 being the most accurate where most parameters are directly measured and level 3 being the simplest where only few parameters need to be known and are used to estimate everything else.
As noted before, the relationship between suction and degree of saturation for a soil is described by the Soil-Water Characteristic Curve (SWCC). Figure 1 shows various SWCCs for soils with different characteristics.
Figure 1: Soil-Water Characteristic Curves
The minimum inputs required (input level 3) to calculate the SWCC for a given soil are the plasticity index, PI, the percentage by weight passing the sieve #200, P200, and the effective grain size corresponding to 60% passing by weight, D60. As can be seen, these are all granulometric parameters of the material. P200 and PI are used for soils that contain a certain fraction of fines and have, therefore, a plasticity index larger than zero. They deliver the family of curves on the right hand side of Figure 1, where it is possible to see that thanks to the fine fraction the degree of saturation varies gradually with suction.
D60, instead, is used for coarser materials that have PI=0. In this case the SWCCs show a sharp change in degree of saturation with suction due to the lack of capillarity.
These curves and the determination of the parameters that define them are thoroughly described in the ME-PDG, Part 2 - Chapter 3, therefore won’t be reported in detail here.
The depth of the water table, which can be considered to vary seasonally or to be constant throughout the year, determines the values of suction at any point in the pavement, which are then used to calculate the degree of saturation by means of the SWCC. From the degree of saturation it is then possible to calculate F_u, the environmental factor for unfrozen or fully recovered material.
If the temperatures are below freezing point, the environmental factor takes the form of F_f, which has a constant value much larger than 1 since frozen materials are much stiffer.
Finally, when the soil is thawing the environmental factor becomes F_r (for recovering material) and is given a value lower than F_u. This value is at a minimum at the beginning of the recovering period, when the material is weakest, and tends gradually to F_u as the material returns to its normal, unfrozen state.
An example of how these factors appear with time for a number of nodes in a pavement is given in Figure 2.
Figure 2: Evolution of F_env with time
In Figure 2, it can be seed that on day 1 some nodes are frozen, some are thawing and the rest are unfrozen. As time passes, the coefficients for the frozen material remain constant for as long as the material stays frozen, while for the recovering material they increase up to the normal unfrozen value. Node 12, for instance, stays frozen until day 3 with F_f = 75, then starts thawing and the value of the environmental coefficient drops to 0.6. As the thawing goes on this value starts to increase and will eventually reach the value of F_u. Similarly, node 17 is already thawing at day 1 and at day 6 the thawing ends. At that point the material has completely recovered and the coefficient has reached a value of F_u of 0.9.

Read More

validation 2

This model has been validated against temperature data from the LTPP database for two test sites with different environmental conditions, one of them (section 1-0101) in Alabama and the other (section 31-3018) in Montana.
For both cases the structure analysed is that described in Table 1.

Layer	Thickness [m]	Specific Heat [J kg-1 K-1]	Thermal Cond. [J m-1 s-1 K-1]	Density [kg m-3]
HMA	0.2	1670	1.40	2400
Granular	5.2	840	0.76	2100

Table 1: Validation – pavement structure
The input parameters used are:

Cloud base factor	0.9
Vapour pressure	1.33mbar
S. s. w. absorptivity	0.98
S. emissivity	0.93

Table 2: Validation – inputs
The highest value for surface absorptivity was chosen as it corresponds to a newly constructed (or re-surfaced) road, which is the case of these test sites.
1.1.1.1 Alabama – latitude of 32.61°
The Figure 1 to Figure 5 show how the predicted monthly average temperatures compare to the measured daily ones at different depths for a test section in Alabama. In this case, measured temperatures were not available after the 12/10/1996.
Figure 1: Alabama – Prediction of surface temperature
Figure 2: Alabama – Prediction of temperature at 0.2m of depth
Figure 3: Alabama – Prediction of temperature at 0.4m of depth
Figure 4: Alabama – Prediction of temperature at 1m of depth
Figure 5: Alabama – Prediction of temperature at 2m of depth
1.1.1.2 Montana – latitude of 47.41°
Figure 6 to Figure 10 refer to a test site in Montana, with temperatures noticeably lower than our previous example.
Figure 6: Montana – Prediction of surface temperature
Figure 7: Montana – Prediction of temperature at 0.2m of depth
Figure 8: Montana – Prediction of temperature at 0.4m of depth
Figure 9: Montana – Prediction of temperature at 1m of depth
Figure 10: Montana – Prediction of temperature at 2m of depth

Read More

sensitivity analysis 2

Calculating the energy absorbed by the ground at a given time using the energy balance and the equations discussed in the previous sections requires a substantial number of inputs. It can be useful, therefore, to know which ones are those that can have a larger effect on the final results and which ones might, instead, require less accuracy.
Table 1 summarises the inputs needed in this part of the analysis. Considering that some of them are variables that are updated during the analysis and some others are meteorological data from external sources, the sensitivity analysis has focussed on the inputs that the user is required to input directly.

Fixed	During analysis	From Met. Station
Latitude	Time of day	Air temperature
Cloud base factor	Day of year	% of cloud cover
Vapour pressure	Surface temperature	Wind speed
S. s. w. absorptivity
S. emissivity

Table 1: Inputs to the temperature model
The structure analysed is given in Table 2. The meteorological data used in this case were extracted from an ME-PDG analysis for a location with the characteristics specified in Table 3.

Layer	Thickness [m]	Specific Heat [J kg-1 K-1]	Thermal Cond. [J m-1 s-1 K-1]	Density [kg m-3]
HMA	0.2	921	1.16	3200
Granular	0.4	754	0.38	2700
Granular	0.4	754	0.38	2700
Granular	5.0	754	0.38	2700

Table 2: Sensitivity analysis – pavement structure

Latitude	33.26º
Cloud base factor	0.9
Vapour pressure	1.33mbar
S. s. w. absorptivity	0.85
S. emissivity	0.93

Table 3: Sensitivity analysis – initial inputs
These values have been kept as the basis for the sensitivity analysis, varying them one at the time while keeping the others constant. Their effect on the simulated surface temperature is shown in Figure 1to Figure 5and is summarised in Table 4.
Figure 1: Sensitivity to latitude
Figure 2: Sensitivity to cloud base factor
Figure 3: Sensitivity to vapour pressure of air
Figure 4: Sensitivity to short wave absorptivity
Figure 5: Sensitivity to surface emissivity

Parameter	Range	Change [ºC] (Summer)	Change [ºC] (Winter)
Latitude	0 – 70 [º]	2.5	9.0
Cloud base factor	0.8 – 0.9	0.9	0.4
Vapour pressure	1.33 – 13.3 [mbar]	7.0	6.0
S. s. w. absorptivity	0.8 – 0.98	2.5	1.0
S. emissivity	0.8 - 1	2.5	2.0

Table 4: Sensitivity analysis – results (surface)
The effect of these parameters on the temperatures calculated at the bottom of the asphalt layer (depth of 0.2m) is, as expected, very similar, but these temperatures also affected by the material parameters such as thermal conductivity k and specific heat c. As shown in Table 5, the effect of varying k and c within the range suggested in the ME-PDG is quite small.

Parameter	Range	Change [ºC] (Summer)	Change [ºC] (Winter)
Specific heat	0.10 – 0.65 [J kg-1 K-1]	0.4	0.6
Thermal conductivity	0.44 – 0.81 [J m-1 s-1 K-1]	0.6	0.7

Table 5: Sensitivity analysis – results (bottom of asphalt layer)

Read More

Finite Differences method 2

The finite difference method is a commonly used iterative technique that allows numerical solution of complex equations. In particular, forward difference methods approximate the solutions to differential equations by replacing derivative expressions with approximately equivalent difference quotients (Dusinberre, 1961).
In the case of heat transfer through the layers of a road pavement, it is assumed that each layer has uniform properties and is of infinite horizontal extent in order to reduce the problem to the one-dimensional case of determining temperatures at any particular depth and time (see Figure 1).
Figure 1: Finite differenced method
As shown in Figure 111, the domain is portioned in space using a uniform vertical mesh of size Δz and in time using a uniform mesh of size Δt. For each layer, the values of material parameters such as density ρ [kg m^-3], thermal conductivity k [J m^-1 s^-1 K^-1] and specific heat c [J kg^-1 K^-1] are known. These values define the values of conductance K_ij (relative to the material between two points of the mesh) and heat capacity C_i (relative to the material around each mesh point) that are then used in the following equations to calculate the new temperatures T’_i from the current temperatures T_i:
                                                                                          (1)
where
with j ≠ i                                                                                                                (2)
with j ≠ i                                                                                                                   (3)
                                                                                                                               (4)
Q_gi is the energy that enters the system at the node i. In our case, only the surface node receives the energy Q_g (that has been calculated in the previous sections), therefore this term is zero for all the internal nodes of the mesh.
A consideration for equations 1-4 defines a criterion for the stability of this method. It can be seen, in fact, that if Δt is chosen large enough, the coefficients F_ii could become negative. This would be absurd because it would mean that the higher the temperature T_i is at the current time step the lower it is going to be at the next. In order to avoid this situation, therefore, the size of the time step has to comply with the following rule in every region of the system:
                                                                                                                                         (5)
which effectively states that the maximum acceptable value for Δt is directly proportional to the square of Δz.

Read More

Energy balance

Figure 1shows the energy sources that are considered when estimating the energy balance at the pavement surface. The energy that enters (or leaves) the ground comes from different sources such as the sun, the objects that surround the pavement and the air. This energy balance can be expressed as (Dempsey, 1985):
Q_g = Q_s + Q_l + Q_c                                                                                                                                            (1)
where
Qg is the energy absorbed by the ground [W m-2]
Qs is the net short wave radiation [W m-2]
Ql is the net long wave radiation [W m-2]
Qc is the convective heat transfer [W m-2]
Figure 1: Energy exchange between pavement surface and environment
It has to be noted that this approach does not include the effects of transpiration, condensation, evaporation and sublimation.
Q_s – net short wave radiation
The short wave radiation can be identified with the high frequency energy radiated by the sun that reaches the pavement surface and part of which is visible to the human eye. It is, therefore, a function of the extraterrestrial radiation incident on the atmosphere at a particular latitude, day and time and of the filtering and reflecting effect of the clouds that may be present at that location (ME-PDG, 2004).
                                                                                                                        (2)
where
a_s is the surface short wave absorptivity (0.8 – 0.98)
R is the potential extraterrestrial radiation [W m^-2]
S_c is the percentage of sunshine
A and B are constants that account for diffuse scattering and absorption by the atmosphere. The ME-PDG suggests values of 0.202 and 0.539 respectively but they will change with location.
The potential extraterrestrial radiation R is a function of latitude Ф [rad], declination δ [rad] and time of day. Its average value during one day, R_a, can be calculated by:
                              (3)
where
Ф is the latitude expressed in radians
δ is the declination expressed in radians and calculated by:
                                                                                              (4)
where day is the day of the year expressed as an integer.
The average value of R_a does not serve the purpose of estimating temperatures on an hourly basis. Nonetheless, assuming a parabolic distribution of R during daytime that starts from sunrise and finishes at sunset (Dempsey, 1985), it is possible to use R_a to estimate R at any given time. The total energy radiated during daytime must be equal to R_a multiplied by 24h, i.e. the areas of Figure 2(a) and (b) must be the same.
Figure 2: Calculating R from R_a (a) (b)
The time of sunrise and sunset for one particular day and latitude can be calculated as follows (University of Oregon, 2009):
Sunrise_angle = [rad]                                                                               (5)
Sunset_angle = [rad]                                                                        (6)
Sunrise_time = [s]                                                              (7)
Sunset_time = [s]                                                               (8)
From what given above, the parabola’s coefficients can be calculated:
                                                                                       (9)
                                                                                         (10)
[s]                                                                                    (11)
for                   (12)
R=0 for  or (13)
Q_l – net long wave radiation
Any object that is at a temperature higher than the absolute zero radiates energy. Since we are dealing with objects at low temperatures (certainly not comparable to the sun’s temperature), they radiate at low frequencies / long wave lengths, well outside the visible spectrum.
The general formula to calculate this energy is given by:
[W m^-2]                                                                                                                         (14)
where:
ε is the object’s emissivity
σ is the Stefan-Boltzman constant [W m^-2 K^-4]
T is the object’s temperature [K]
The net long wave radiation that contributes to the energy balance at the pavement’s surface is calculated considering the incoming (positive) energy radiated by the air and the outgoing (negative) energy radiated by the pavement itself. Therefore:
[W m^-2]                                                                                                          (15)
where:
R_ldc is the clear sky downward radiation:
                                                                                                                                         (16)
where the air emissivity εa can be calculated as a function of vapour pressure and Ta is the air temperature (Jacobs, 2004).
Q_x is the outgoing radiation:
                                                                                                                                           (17)
where εs is the surface emissivity and Ts is the surface temperature.
C_cc is the cloud cover correction factor:
                                                                                                                                      (18)
where N is the cloud base factor (0.9-0.8 depending on height) and W is the percentage of clouds.
Q_c – convective heat transfer
The pavement’s surface exchanges heat with the surrounding air by convection. This energy transfer is function of air and surface temperature and wind speed. The ME-PDG suggests the following equation:
                                                                                    (19)
where:
T_m = average of T_a and T_s [K]
U = average daily wind speed [m s^-1]

Read More