Nomenclature

1 Introduction

The consolidation of soils due to the dissipation of excess pore pressures is a common issue in geotechnical engineering. A considerable number of researches on the consolidation theory for unsaturated soils have been conducted, and achieved great progress. Scott^[1] estimated the consolidation of unsaturated soils with occluded air bubbles. Biot^[2] proposed a general consolidation theory that is suitable for analyzing unsaturated soil with occluded air bubbles. Barden^[3] presented an analysis of one-dimensional (1D) consolidation of compacted unsaturated clay. Assuming that the air and water phases are continuous, Fredlund and Hasan^[4] proposed a 1D consolidation theory, in which two partial differential equations (PDEs) were used to describe the dissipation processes of pore pressures in unsaturated soils. This theory is now widely accepted. Later, Dakshanamurthy et al.^[5] extended the 1D consolidation theory proposed by Fredlund and Hasan^[4] to analyze the three-dimensional (3D) case. By assuming that all the soil parameters remain constant during consolidation, Fredlund et al.^[6] presented a simplified form of 1D consolidation equations for unsaturated soils.

Several analytical approaches have been developed since the inception of the 1D consolidation theory of unsaturated soils. Using the Laplace transformation and Cayley-Hamilton techniques, Qin et al.^[7-8] presented an analytical solution to 1D consolidation for a single-layer unsaturated soil subjected to an instantaneous loading and an exponentially time-dependent loading when considering a permeable top and impermeable base of the soil layer. On the other hand, Qin et al.^[9] also gave semi-analytical solutions under different boundary conditions. Zhou et al.^[10] used two alternative terms, and , which consist of pore-air and pore-water pressures, to convert the nonlinear inhomogeneous PDEs into traditional homogeneous PDEs, and then the solutions for single and double drainage conditions were obtained using the separation of variables method. Shan et al.^[11] also deduced an exact solution of 1D consolidation for unsaturated soils with different boundary and loading conditions by adopting the separation of variables method. However, the final equations have been left undisclosed as a result of cumbersome derivation, and it should be noticed that the solutions given by Shan et al.^[11] are difficult to be followed by engineers. By adopting the eigen-function expansion method and Laplace transform techniques, Ho et al.^[12] discussed a simple yet precise analytical solution for 1D consolidation of an unsaturated soil deposit under homogeneous boundary conditions subjected to an instantaneous loading. Ho and Fatahi^[13] gave an analytical solution for the two-dimensional (2D) plane strain consolidation of an unsaturated soil stratum subjected to time-dependent loadings using the same method as Ho et al.^[12].

This paper presents semi-analytical solutions to predict the dissipation of pore-air and pore-water pressures and the settlement for unsaturated soils deposit using the 1D consolidation theory proposed by Fredlund and Hasan^[4] under different initial conditions, homogenous boundary conditions, and time-dependent loadings. To obtain final solutions, inhomogeneous governing equations for unsaturated soils are first transformed into homogeneous equations, and then the equations are solved using the Laplace transform. It is found that the current solutions are more general and in good agreement with the existing solutions in the literature, and they also can be degenerated into that of Terzaghi's consolidation^[14] for fully saturated soils. Finally, examples are given to illustrate consolidation behavior of unsaturated soils. Changes in the pore-air and pore-water pressures and settlement under varying air to water permeability ratios, time-dependent loadings, loading parameters, initial and boundary conditions are sufficiently demonstrated in this paper.

2 Mathematical model 2.1 Governing equations

Fredlund and Hasan^[4] proposed 1D consolidation theory for an unsaturated soil stratum, whose width is considered to be infinite, and the thickness is measured and denoted by h. The soil system of interest indicates that the dissipations of air and water and settlement only occur along the vertical direction (z-direction), as shown in Fig. 1.

The main assumptions for 1D consolidation of unsaturated soils are made as follows:

(ⅰ) Solid particles and water phase are incompressible.

(ⅱ) Air and water phases are assumed to be continuous.

(ⅲ) The effects of air diffusing through water, air dissolving in the water, and the movement of water vapor are ignored.

(ⅳ) The coefficients of permeability with respect to air and water phases and volume change for the soil remain constant throughout the consolidation process.

(ⅴ) The loading and deformation take place only along the vertical direction.

The governing equations for water and air phases after the application of the load q are written as follows^[11]:

(1)

(2)

where u_a and u_w are the pore-water and pore-air pressures (kPa), C_a and C_w are interactive constants with respect to the air and water phases, respectively, C_σ^a and C_σ^w are the consolidation coefficients for the air and water phases, respectively, and C_v^a and C_v^w are coefficients of volume change with respect to the air and water phases, respectively. The consolidation parameters can be expressed as follows:

(3a)

(3)

(3c)

(3d)

(3e)

(3f)

where m_1k^w and m₂^w are the coefficients of water volume change with respect to a change in the net normal stress σ-u_a and matric suction u_a-u_w, respectively, m_1k^a and m₂^a are the coefficients of air volume change with respect to a change in σ-u_a and u_a-u_w, respectively, and the subscript ‘k' stands for K₀-loading. k_w and k_a are the coefficients of water and air permeability, respectively, and the unit weight of water γ_w = 9.8 kN/m³. The acceleration of gravity g=9.8 m/s², and S_r0 and n₀ are the initial degree of saturation and initial porosity, respectively. The molecular mass of air M=0.029 kg/mol, and the absolute pore-air pressure u_a⁰ = u_a⁰ + u_atm, where u_atm is the atmospheric pressure. The universal gas constant R=8.314 J/(mol·K), and T is the absolute temperature.

2.2 Initial and boundary conditions

The initial conditions are as follows:

(4a)

(4b)

where a_a and a_w are the changing rates of initial pore-air and pore-water pressures along the depth, respectively, b_a and b_w are initial pore-air and pore-water pressures at the top surface, respectively.

The top and bottom boundaries are considered to be permeable or impermeable to air and water phases, and then the possible boundaries are as follows:

The top boundaries are

(5a)

(5b)

The bottom boundaries are

(6a)

(6b)

As shown later, the method in this paper can solve Eqs. (1) and (2) under any combination of the above top and bottom boundary conditions.

3 Derivation of semi-analytical solutions 3.1 General semi-analytical solutions

Equations (1) and (2) can be rewritten as follows:

(7)

(8)

where

Equations (7) and (8) can be transformed into an equivalent set of PDEs of and ,

(9)

(10)

Details of the derivation and the variables of , , Q₁, Q₂, β₁, and β₂ can be found in Appendix A.

It is shown that if Eqs. (9) and (10) are solved by the separation of variables method, it only can be solved under the conditions of the identical permeability of water and air, such as the boundary conditions of single drainage or double drainage^[10]. In the paper, the Laplace transform is used to solve Eqs. (9) and (10) for overcoming the above limitation. By applying the Laplace transform to Eqs. (9) and (10), respectively, the following equations can be obtained:

(11)

(12)

where Q(s) is the result of the Laplace transform of upon the time t, , and .

General solutions of two-order ordinary differential equations (11) and (12) are

(13)

(14)

where and are special solutions to Eqs. (11) and (12), χ₁² = s/Q₁, χ₂² = s/Q₂, and C₁, C₂, D₁, and D₂ are arbitrary functions to be determined from boundary conditions.

From Eqs. (A10) and (A11) shown in Appendix A, we have

(15)

(16)

Combining Eqs. (13) and (14) with Eqs. (15) and (16) gives

(17)

(18)

Taking the derivatives with respect to Eqs. (17) and (18) leads to

(19)

(20)

3.2 Solutions of initial boundary value problems

Applying the Laplace transform to the initial condition (see Eq. (4)) gives

(21a)

(21b)

Since the pore-air pressure dissipates more easily than the pore-water pressure, it may be impractical to consider both the top and bottom of the soils permeable to water but impermeable to air. Therefore, the other possible boundary conditions can only be the following four cases.

Case 1 The top boundary is permeable to water and air, and the bottom boundary is impermeable to water and air. Applying the Laplace transform to the boundary condition for Case 1 leads to

(22)

Substituting Eq. (22) into Eqs. (17) -(20) yields

(23)

(24)

Case 2 The top boundary is permeable to water and air, and the bottom boundary is permeable to water and air. Applying the Laplace transform to the boundary condition for Case 2 leads to

(25)

Substituting Eq. (25) into Eqs. (17) and (18) yields

(26)

(27)

with and .

Case 3 The top boundary is impermeable to water and permeable to air, and the bottom boundary is impermeable to water and air. Applying the Laplace transform to the boundary condition for Case 3 leads to

(28)

Substituting Eq. (28) into Eqs. (17), (19), and (20) yields

(29)

(30)

with δ₁ = sinh(χ₁h) cosh(χ₂h), δ₂ = sinh(χ₂h) cosh(χ₁h), δ₃ = sinh(χ₁h) cosh(χ₂(h − z)), and δ₄ = sinh(χ₂h) cosh(χ₁(h − z)).

Case 4 The top boundary is permeable to water and air, and the bottom boundary is impermeable to water and permeable to air. Applying the Laplace transform to the boundary condition for Case 4 leads to

(31)

Substituting Eq. (31) into Eqs. (17), (18), and (20) yields

(32)

(33)

with δ₅ = cosh(χ₁h) cosh(χ₂z), δ₆ = cosh(χ₂h) sinh(χ₁z), δ₇ = cosh(χ₁h) cosh(χ₂(h − z)), and δ₈ = cosh(χ₂h) sinh(χ₁(h − z)).

According to the approach with two stress-state variables for unsaturated soils^[6],

(34)

where m_1k^s = m_1k^a + m_1k^w, and m₂^s = m₂^a + m₂^w.

By applying the Laplace transform to Eq. (34), we have

(35)

The settlement of unsaturated soil layer in the Laplace transformed domain can be calculated by

(36)

In conclusion, the expressions of , and are obtained for 1D consolidation in unsaturated soils under any boundary condition except for the boundary impermeable to air and permeable to water, which is not practical.

By adopting the Crump's method^[15] to perform the Laplace inversion, the semi-analytical solutions of pore-air pressure, pore-water pressure, and soil layer settlement in the time domain can be obtained. Details of the Crump's method can be found in Appendix B.

3.3 Verification

With the consideration of the instantaneous loading, the parameters of Case 1 in Eqs. (23) and (24) become A_σ = 0, W_σ = 0, β₁ = 0, and β₂ = 0. Substituting the expressions of Q₁, Q₂, c₁₂, and c₂₁ into Eqs. (23) and (24), the following expressions can be obtained:

(37)

(38)

where

Equations (37) and (38) are similar to the existing solutions provided by Qin et al.^[7]. Furthermore, Eq. (24) can be degenerated into the solution to Terzaghi's consolidation equation for saturated soils^[17].

4 Examples and discussion

In the case study, the parameters are assumed as follows: h = 10 m, n₀ = 50%, S_r0 = 80%, k_w = 10⁻¹⁰ m/s, m_1k^s = −2.5 × 10⁻⁴ kPa⁻¹, m_1k^w = −0.5 × 10⁻⁴ kPa⁻¹, m₂^s = −1.0 ×10⁻⁴ kPa⁻¹, m₂^w = −2 × 10⁻⁴ kPa⁻¹, u_a⁰ = 5.0 kPa, u_w⁰ = 40 kPa, u_atm = 101.3 kPa, q₀ = 100 kPa. Since the pore-air pressure dissipates more easily than the pore-water pressure, the air permeability k_a varies from 10^-10 m/s to 10^-8 m/s, while the water permeability k_w is kept constant at 10^-10 m/s.

Four different types of loadings, i.e., instantaneous, ramp, exponential, and sinusoidal loadings, are simulated and discussed in separate examples. Furthermore, by adopting the boundary condition for Case 1, which is expressed by Eq. (22), the investigation is carried out to study the effects of the air to water permeability ratio (k_a/k_w) and different depths on changes of the pore-air pressure, pore-water pressure, and settlement under the applied loadings. In addition, parametric studies are then conducted to investigate the effects of the loading parameters on the changes of the pore pressures and the settlement.

4.1 Consolidation under instantaneous loading

In order to compare with the results from Qin et al.^[7], the variations of pore-water and pore-air pressures with the time at z=8 m are computed for the single drainage with the constant initial pore-air and pore-air pressures distributions. Figure 2 presents the results under the instantaneous loading for different k_a/k_w values. It is found that results of the two solutions are identical. Moreover, the proposed solutions are useful for general applications and can be used under different time-dependent loadings, boundary and initial conditions.

4.2 Consolidation under ramp loading

The ramp loading, which can simulate the construction loading, can be expressed as

(39)

where q₀ is the initial surcharge, t is the time, and a is the linear loading rate. In this section, an initial surcharge q₀ = 100 kPa and a linear loading rate a=10^-4 kPa/s are adopted. Substituting the Laplace transform of Eq. (39) into Eqs. (23) -(24), (26) -(27), (29) -(30), (32) -(33), and (36), respectively, the pore-air pressure, pore-water pressure, and settlement for different boundary conditions can be obtained.

Figure 3 illustrates the dissipation patterns of pore-air and pore-water pressures varying with k_a/k_w. It is found that the values of k_a/k_w have significant effects on the dissipation rates of pore-air and pore-water pressures. The larger the value of k_a/k_w is, the more quickly the pore-air pressure dissipates at the later stage of consolidation and the pore-water pressure dissipates at the intermediate stage of consolidation. At the later stage, the dissipation patterns of the pore-water pressure for different values of k_a/k_w converge into a single curve because the pore-air pressure has dissipated completely. As the loading increases linearly with the time, the pore-air and pore-water pressures increase exponentially and reach the highest value at about 10⁶ s. Furthermore, it can be observed that when the pore-air pressure dissipates to zero, it is about the same time that the plateau period occurs, and this indicates the end of the early stage.

For the parametric study, the linear loading rate a ranges from 10^-3 kPa/s to 10^-5 kPa/s. Figure 4 represents the dissipation patterns of pore-air and pore-water pressures for different values of the linear loading rate a. It is found that the increasing loading rate a results in a quicker increase in the ramp loading, which causes higher pore pressures after 10⁵ s, 10⁶ s, or 10⁷ s. The dissipation patterns of the pore-water pressure for different values of the linear loading rate a converge into almost a single curve. The plateau then appears when the pore-air dissipates completely, which is similar to Fig. 3.

Figure 5 shows the dissipation patterns of pore-air and pore-water pressures for different depths with k_a/k_w = 10 and the linear loading rate a=10^-4 kPa/s. The closer to the top boundary the depth is, the earlier the pore-air pressure dissipates at the beginning. Pore-air and pore-water pressures increase exponentially and reach the highest value at about 10⁶ s. The plateau period at the later stage is also observed as that in Figs. 3(b) and 4(b). It is noteworthy that plateau gets longer with the depth.

Figure 6(a) illustrates the change of the normalized settlement (w^∗ = w/(m_1k^s q₀h)) for different values of k_a/k_w for the linear loading rate a=10^-4 kPa/s. Figure 6(b) shows the change of the normalized settlement for different values of the linear loading rate a with k_a/k_w =10. At the beginning, the settlement rates for different values of k_a/k_w are very similar, as shown in Fig. 6(a). After 10⁶ s, as k_a/k_w increases, the settlement increases, too. On the other hand, the higher linear loading rate a leads to an increase in the settlement rate during the later stages of consolidation as shown in Fig. 6(b). Also these settlement curves tend to be the same at the later stage for different values of k_a/k_w or different linear loading rates of a.

4.3 Consolidation under exponential loading

Being a simplified mathematical form of the construction loading, the exponential loading is expressed as

(40)

where q₀ is the initial surcharge (kPa), A is the parameter influencing the loading magnitude, and b is the loading parameter controlling the rate of exponential loading. In this study, an initial surcharge q₀ = 100 kPa, A=1, and the exponential loading rate b=5 × 10^-4s^-1 are adopted. Substituting the Laplace transform of Eq. (40) into Eqs. (23) -(24), (26) -(27), (29) -(30), (32) -(33), and (36), the pore-air pressure, pore-water pressure, and settlements for different boundary conditions can be obtained.

Figure 7 illustrates the pore-air and pore-water pressure dissipation patterns at z=8 m for different k_a/k_w. The application of loading results in the exponential increase in both pore-air and pore-water pressures at the beginning of consolidation. At the later stage, when the external loading becomes constant, pore pressures no longer increase but gradually decrease with the time, and eventually return to zero. It is found that the larger k_a/k_w leads to an increase in the dissipations of the pore-air pressure at the later stage of consolidation and the pore-water pressure at the intermediate stage of consolidation. It should be noted that, during the later stage, the dissipation patterns of the pore-water pressure for different values of k_a/k_w converge into a single curve similar to Fig. 3(b).

For the parametric study, the exponential loading rate b ranges from 5× 10^-3s^-1 to 5× 10^-5 s^-1. Figure 8 plots the dissipation patterns of pore-air and pore-water pressures for different values of the exponential loading rate b. It is found that the increasing loading rate b may result in an earlier increase in the exponential loading. Then, when the pore pressures reach the maximum value, the first plateaus appear as shown in Fig. 8, and the larger the exponential loading rate b is, the longer the plateau period gets. It should be noted that only pore-water pressure experiences the second plateau period observed in Fig. 8 (b). During this time, the dissipation patterns of pore-water pressures for different values of the loading rate b begin to converge into a single curve at the later stage.

Figure 9 depicts the dissipation patterns of pore-air and pore-water pressures for different depths with k_a/k_w = 10 and the exponential loading rate b=5× 10^-4s^-1. At the beginning, the pore-air and pore-water pressures increase at the same rate, and it is found that the pore-air pressure dissipates more slowly with the depth, but it is prone to dissipate completely at almost the same time for different depths as shown in Fig. 9(a). A similar dissipation behavior can be observed for the pore-water pressure as shown in Fig. 9(b). Similar to Figs. 7(b) and 8(b), the plateau period becomes longer as the point of interest is further away from the top surface.

Figure 10(a) shows the change of the normalized settlement for different values of k_a/k_w with the exponential loading rate b=5× 10^-4 kPa/s. Figure 10(b) shows the change of the normalized settlement for different values of the exponential loading rate b with k_a/k_w = 10. It is found that the settlement increases with the external loading. At the beginning, the settlement with different values of k_a/k_w proceeds at the same rate in Fig. 10(a). After 10⁴ s, the rate of settlement increases with the increase in the value of k_a/k_w. On the other hand, the higher exponential loading rate b increases the rate of settlement during the early stage of consolidation as shown in Fig. 10(b). However, these settlement curves tend to converge into a single curve at the later stage regardless of different values of k_a/k_w or different exponential loading rates b.

4.4 Consolidation under sinusoidal loading

In this section, the application of sinusoidal loading expressed in Eq. (2) is investigated,

(41)

where q₀ is the initial surcharge (kPa), c is the parameter that influences the loading amplitude, and θ is the angular frequency (rad/s) for the sinusoidal loading function. In this study, an initial surcharge q₀ = 100 kPa, the parameter c=1, and the angular frequency θ =2π × 10^-7 rad/s are adopted, for simulating the ground surface loading induced by pumping in and out oil in the oil storage tank, and stacking and removing the container in the port yard, etc. Substituting the Laplace transform of Eq. (41) into Eqs. (23) -(24), (26) -(27), (29) -(30), (32) -(33), and (36), the pore-air pressure, pore-water pressure, and settlement for different boundary conditions are obtained.

Figure 11 illustrates the dissipation patterns of pore-air and pore-water pressures varying with k_a/k_w. After applying the sinusoidal loading, the pore-air pressure does not fully dissipate after a long period of time but rather oscillates continuously with constant amplitudes in Fig. 11(a). At the beginning of the periodical oscillations of pore-air pressure, the curve of pore-water pressure also appears to be oscillating continuously, but the mean value of pore-air pressure decreases in Fig. 11(b). This phenomenon occurs because the simulated loading function consists of loading-unloading curves, which influence the changes of the pore pressures. On the other hand, the larger k_a/k_w is, the smaller the amplitudes of the pore pressures oscillations are.

For the parametric study, the angular frequency θ includes 2π × 10^-7 rad/s and 2π × 10^-8 rad/s. Figure 12 presents the dissipation patterns of pore-air and pore-water pressures for different values of the angular frequency θ. It is found that the increasing angular frequency θ results in a large amplitude of pore pressures curves. It can also be observed that the smaller the angular frequency θ is, the later the patterns of the pore pressures start to oscillate. The amplitudes of pore-water pressure curve decrease at the beginning of the periodical oscillations of pore-air pressure curve, which is the same as that observed in Fig. 11(b).

Figure 13 shows the dissipation patterns of pore-air and pore-water pressures for different depths with k_a/k_w = 10 and the angular frequency θ =2π × 10^-7 rad/s. The closer to the top boundary the depth is, the earlier the pore pressures dissipate and the smaller the amplitudes of pore-air and pore-water pressures are.

Figure 14(a) illustrates the change of settlement for different values of k_a/k_w with the angular frequency θ =2π × 10^-7 rad/s. Figure 14(b) shows the change of the normalized settlement for different values of the angular frequency θ with k_a/k_w = 10. It is found that the larger k_a/k_w will result in the higher amplitude of the settlement curve. On the other hand, the higher angular frequency θ increases the frequency of periodical oscillations for settlement at the later stages of consolidation compared with the smaller angular frequency θ in Fig. 14(b). The smaller the angular frequency θ is, the later the periodical oscillations of settlement start.

4.5 Parametric study on isochrones of pore pressures under different initial conditions

From Eq. (4), different initial conditions can be outlined, when parameters of a_a, b_a, a_w, and b_w take different values. There are three kinds of initial conditions, as shown in Fig. 15.

(a) Triangle , and u_w⁰(z) = 4z kPa;

(b) Rectangle u_a⁰(z)=5 kPa, and u_w⁰(z) =40 kPa;

(c) Trapezium , and u_w⁰ (z) = 40 + 4z kPa.

Figures 16 and 17 present the isochrones of pore-air and pore-water pressures against the time t under three different initial pore pressure distributions. k_a/k_w = 10 is adopted for analysis. It can be observed that the initial condition has a noticeable effect on the pore pressure distribution along the depth. At the early stage of consolidation, the pore pressure distribution along the depth is similar with initial conditions. As the pore pressures dissipate, the pore pressure isochrones tend to be consistent.

5 Concluding remarks

In this paper, general semi-analytical solutions to 1D consolidation equations for unsaturated soils are proposed by introducing two related variables and and the Laplace transformation techniques. The proposed semi-analytical solutions are more general for practical applications and simpler compared with the previous available solutions in the literature. The verification for instantaneous loading shows good agreement with the existing analytical results in the literature. The proposed solutions can also be degenerated into that of Terzaghi's consolidation for fully saturated soils. In the case studies, the consolidation behavior of unsaturated soils is investigated under three different initial conditions, homogeneous boundary conditions and four time-dependent loadings. It is found that the smaller k_a/k_w is, the more slowly the pore pressures dissipate at the beginning. The pore-air pressure dissipates completely in a relatively short time, and the dissipation curves of the pore-water pressure tend to be the same after the pore-air pressure has completely dissipated. The loading parameters a, b, and θ have significant effects on the dissipation of pore-air pressure at the early stage, while the effects on the dissipation of pore-water pressure are observed at the later stage. For the parameter depth z, it is shown that the pore-air pressure tends to dissipate faster at the depth close to the top surface. The initial condition has distinctly affected the pore pressure distributions upon the depth. The pore pressure distribution along the depth is similar to the initial condition at the early stage of consolidation, and with the dissipation of pore pressures, the isochrones of pore pressures tend to be consistent.

Appendix A

Multiplying Eqs. (7) and (8) by arbitrary constants c₁ and c₂, respectively, and adding these two equations together lead to the following equation:

(A1)

Equation (A1) could be transformed into a conventional diffusion equation with the variable by introducing a constant Q that satisfies the following relationships:

(A2)

(A3)

In order to make Eqs. (A2) and (A3) hold, the constant Q must satisfy the following condition:

(A4)

Equation (A4) is a quadratic equation in terms of Q that has the following two roots Q₁ and Q₂:

(A5)

when Q = Q₁, the solutions to Eqs. (A2) and (A3) for variables c₁ and c₂ are c₁₁ and c₂₁, whereas when Q = Q₂, the solutions for variables c₁ and c₂ are c₁₂ and c₂₂.

Without loss of generality, it is possible to assume that c₁₁ = c₂₂ = 1, and c₁₂ and c₂₁ can be expressed as follows:

(A6)

(A7)

Equation (A1) can therefore be rewritten in the following forms:

(A8)

(A9)

where

(A10)

(A11)

(A12)

(A13)

The initial conditions for and are as follows:

(A14)

(A15)

Applying the Laplace transform to Eqs. (A14) and (A15), respectively, we have

(A16)

(A17)

Appendix B

Given a function F(s) defined for complex values of s, the routine estimates values of its inverse Laplace transform by the Crump's method^[15]. The Crump's method applies the epsilon algorithm to the summation in Durbin's Fourier series approximation^[16],

(B1)

In this study, the method for taking values of parameters is explained as follows:

(ⅰ) τ = t_fac × max(0.01, t_j), where t_fac = 0.8.

(ⅱ) , where α_b should be specified equal to or slightly larger than the value of α, and α has two alternative interpretations, i.e., α is the smallest number such that |f(t)| ≤ m × exp(αt) for large t, and α is the real part of the singularity of F(s) with the largest real part. E_r is the required relative error in the values of the inverse Laplace transform, and thus E_r must be in the range 0 ≤ Er < 1.0.

(ⅲ) The values of t_j, j=1, 2, …, n, must be supplied in the monotonically increasing order.

The routine calculates the values of the inverse function f(t_j) in the decreasing order of j.