A Regular Integral Equation Formalism for Solving the Standard Boussinesq’s Equations for Variable Water Depth

This paper begins with a question of existence of a regular integral equation formalism, but different from the existing usual ones, for solving the standard Boussinesq’s equations for variable water depth (or Peregrine’s model). For the question, a pseudo-water depth parameter, suggested by Jang (Commun Nonlinear Sci Numer Simul 43:118–138, 2017), is introduced to alter the standard Boussinesq’s equations into an integral formalism. This enables us to construct a regular (nonlinear) integral equations of second kind (as required), being equivalent to the standard Boussinesq’s equations (of Peregrine’s model). The (constructed) integral equations are, of course, inherently different from the usual integral equation formalisms. For solving them, the successive approximation (or the fixed point iteration) is applied (Jang 2017), whereby a new iterative formula is immediately derived, in this paper, for numerical solutions of the standard Boussinesq’s equations for variable water depth. The formula, semi-analytic and derivative-free, is shown to be useful to observe especially the nonlinear wave phenomena of shallow water waves on a beach. In fact, a numerical experiment is performed on a solitary wave approaching a sloping beach. It shows clearly the main feature of nonlinear wave characteristics, which has reached good agreement with the known (numerical) solutions. Hence, while being theoretical but fundamental in nonlinear computational partial differential equations, the question raised in the study may be solved.


Introduction
This paper concerns a theoretical question, which may be of a fundamental area of mathematical and computational methods applied to problems arising in the mathematical sciences and engineering. One may wonder if there exists a regular integral equation formalism, as radically distinct from existing ones, for integrating the standard Boussinesq's equations for variable water depth h (at still water), which were derived by Peregrine [1] as (1.2) subject to the initial conditions, Here, η, u and g denote free surface elevation, depth averaged velocity and acceleration of gravity, respectively. The nonlinear partial differential equations of (1.1) and (1.2), also known as Peregrine's model, account for the mass and momentum conservations, respectively. They physically describe shallow water waves (or long waves) on an uneven (or variable) sea bottom topography, the effect of which on propagating water waves is an issue of scientific and engineering importance (e.g., long waves on a beach, wave shoaling, and so on).
Especially when a special case of constant water depth is considered, Peregrine's model, (1.1) and (1.2), reduces to the classical Boussinesq's equations or the (single) classical Boussinesq's equation [2], whose regular integral equation formalism has recently been constructed by Jang [3]. Thus, the present study can be viewed as an extension of the previous work [3] to include the effect of a variable bottom topography on the wave propagation.
The above question of existence will be resolved, in this study, by explicitly constructing a regular integral equation formalism for the initial value problem of (1.1) and (1.2) together with (1.3). To this end, based on the basic idea of the procedures suggested by Jang [3], we first introduce a pseudo-water depth parameter in the present study to build a nonlinear integral formalism. This makes it possible to construct nonlinear integral equations of second kind, which are equivalent to the (present) initial value problem. Here, it should be emphasized that the integral equations constructed here are formulated from the two coupled nonlinear partial differential equations of (1.1) and (1.2), while the previous integral equation formulation by Jang [3] resulted from just a single (not coupled) nonlinear partial differential equation. Furthermore, it should be noted that the integral equations formulated here are not (inherently) related with, for example, the weak formulations of the finite element approach and the singular boundary integral equations often formulated by the boundary element method.
Applying the successive approximation to the present integral equation formalism (or the nonlinear integral equations of second kind constructed above), we derive a new formula, which is semi-analytic and derivative-free. It is a recurrence relation generating a sequence of improving approximate solutions of the initial value problem for Peregrine's model (or the standard Boussinesq's equations for variable water depth). It is straightforward, in the sense that the only thing required for iterating the recurrence relation is just the conventional numerical regular (not singular) integration in each iteration step. And, it is stable, because the integral equations constructed here are considered to be nonlinear integral equations of Hammerstein type, which are known to have the usual stability properties [4]. Further, the formula derived here gives an integral framework, which could be applicable to other more complicated models in a systematic way, and provides a high convergence rate for achieving an accurate solution, as demonstrated in the numerical experiment on a solitary wave.

Modifying the Boussinesq's Equations for Variable Water Depth
In this section, employing the concept of a pseudo-parameter proposed by Jang [3], we shall present a modified form of the standard Boussinesq's equations for variable water depth, (1.1) and (1.2), which is essential for constructing an integral equation formalism for solving the governing equations, as will be shown later.
We start by introducing a pseudo-parameter h 0 > 0, which will be referred to as a pseudo-water depth parameter in the present study. It is to be incorporated into (1.1) and (1.2), whereby we suggest a (h 0 -parameterized) modified form of (1.1) and (1.2), as in the following proposition.
Proposition 1 (modified system) Let h 0 > 0 be a pseudo-water depth parameter. Then, the standard Boussinesq's equations for variable water depth, (1.1) and (1.2), can alternatively be written in terms of h 0 as Here, ϕ and ψ represent the forcing terms in the form Proof With a pseudo-water depth parameter h 0 > 0, we add the term h 0 u x to both sides of (1.1) and the term − h 2 0 u x xt /3 to both sides of (1.2). This results in immediately the two equations of (2.1) and (2.2), with ϕ and ψ explicitly given in (2.3) and (2.4), completing the proof.

Remark 1
The pseudo-water depth parameter, h 0 > 0, multiplies (partial) derivatives of depth averaged velocity u, in the above Proposition. This is in sharp contrast to the case of Jang [3], where a pseudo-parameter multiplied wave elevation itself (not its derivative).
It should be noticed that the modified system of (2.1) and (2.2), incorporated with the pseudo-water depth h 0 > 0, is mathematically still equivalent to the original system of the governing equations of (1.1) and (1.2). Furthermore, the homogeneous version of the modified system of (2.1) and (2.2) has the form of the linearized classical Boussinesq equations, whose dispersion relation is characterized by (e.g., [5]) in which ω and k denote wave frequency and wave number, respectively, and c 0 characteristic velocity It is convenient to define a frequency function from (2.5) as the below.

Definition 1 (frequency function)
where R indicates the set of real numbers and R + the set of positive reals, is a (frequency) function such that

Building Integral Formalism
This section is devoted to building an integral formalism for the standard Boussinesq's equations for variable water depth, which is necessary to construct integral equations as will be discussed in the next section. We start by applying the integral transform techniques [3], to the modified system of (2.1) and (2.2) discussed in Proposition 1.

Laplace and Fourier Transforms
Let us take Laplace transforms L of wave elevation η and depth averaged velocity u with respect to time t, respectively, i.e., for a parameter s. We then obtain a system of ordinary differential equations from (2.1) and (2.2), as in the following.
Lemma 1 (ordinary differential equations) Let Laplace transforms L of η and u with respect to time t be denoted by η * and u * , respectively. Then, the coupled partial differential equations of (2.1) and (2.2) can be converted into ordinary differential equations for η * and u * , which have the form Proof We take Laplace transform L on (2.1) and (2.2), respectively, utilizing the identities from Laplace transform table. This, combined with the initial conditions (1.3), immediately leads to (3.3) and (3.4), which completes the proof.
Next, we denote Fourier transforms F of η * (x, s) and u * (x, s) with respect to space x by and respectively, for a parameter k. Then, their inverse images, via the inverse Fourier transform F −1 , may be written in the symmetric form (3.9) and respectively. We restrict our attention to the case where the wave motion of (1.1) and (1.2) vanishes in the far field, which can be described below.

Remark 2
The initial conditions (1.3) should be consistent with Assumption 1.
Under Assumption 1, carrying out F on the ordinary differential equations of (3.3) and (3.4) makes them reduce further to algebraic equations, being exhibited by the next Lemma.
The next Lemma deals with solving the algebraic Eqs. (3.12) and (3.13) in Lemma 2, where the calculation required for the solution may be long (and tedious).
Lemma 3η * andū * in the system of (3.12) and (3.13) can be represented as

17)
where ω B is the frequency function defined in Definition 1.
Proof (3.12) and (3.13) can be arranged as or alternatively, (3.18) in matrix notation. By applying inverse operation on (3.18), we are led to Here,η * is calculated as where (2.7) has been used. Similarly,ū * is determined from (3.19) as follows u * = 1 in which the frequency squared ω 2 B comes from (2.7). This completes the proof.

Inverse Transforms
In this subsection, we will recover wave elevation η and depth averaged velocity u from the results of the above Lemma, that is, we look for the inverse images forη * andū * . We first, for convenience, denote each term in the right hand side of (3.16) by and respectively, so thatη * in (3.16) can be simply represented as the sum In the similar way, (3.17) is expressed as by denoting each term in the right hand side of (3.17) by and respectively.

Lemma 4
Wave elevation η is recovered by taking the two successive inverse Fourier and Laplace transforms on both sides of (3.24). i.e., since Laplace transform formula L[cos(at)] = s/(s 2 +a 2 ) for a ∈ R from Laplace transform table. Thus, Note that the real part of the integrand in the last in (3.35) is even in k, while the imaginary part is odd in k. Therefore, (3.35) reduces to This completes the proof of (3.31).
(b) Derivation (3.32) Inverse Laplace transform of (3.21) becomes from the formula of Laplace transform a/(s 2 +a 2 ) = L[sin(at)] for a ∈ R. With the notation * (Fourier convolution), inverse Fourier transform of (3.36) can be calculated as Here, we notice that the resultant is real (i.e., its imaginary part vanishes) and the integrand of the last double integral in the above is even in k. This completes the proof of (3.32).
(c) Derivation (3.33) Denoting by the symbol • the Laplace convolution, we have the inverse Laplace transform of (3.22) written as a Duhamel (or convolution) integral where the property s/(s 2 + a 2 ) = L[cos(at)] for a ∈ R and the convolution theorem have been utilized. We further perform F −1 on the above, leading to where the imaginary part of the integral cancels out because its integrand is odd in k. This completes the proof of (3.33).
(d) Derivation (3.34) Inverse Laplace transform L −1 of f 4 in (3.23) would be a Laplace convolution (with the convolution notation •), in which Laplace transform formula a/(s 2 + a 2 ) = L[sin(at)] for a ∈ R has been used. Returning to definition of (Laplace) convolution •, we have a Duhamel integral Taking the inverse Fourier transform of (3.38) gives Notice that the integrand of the resulting integral is real and even in k. This completes the proof of (3.34).
In the same way as the above, depth averaged velocity u can also be uncovered, which is discussed in the next Lemma.

Lemma 5 The inverse image ofū
where each term can be found explicitly as Proof (a) Derivation (3.40) Let us take the inverse Laplace transform on (3.26) by using the formula a/(s 2 + a 2 ) = L[sin(at)] for a ∈ R. Then, F −1 (L −1 g 1 ) would be This completes the proof of (3.40).
The above can be inverse-Fourier transformed as which completes the proof of (3.41).
for a ∈ R, we perform the inverse Laplace transform on (3.28) Noting that L −1 L = I, where I denotes the identity, we have where the notation • indicates Laplace convolution. We next take further the inverse Fourier Transform F −1 on the above Duhamel integral, providing This completes the proof of (3.42).
(d) Derivation (3.43) With the use of s/(s 2 + a 2 ) = L[cos(at)] for a ∈ R and the notation • for Laplace convolution, the inverse Laplace transform of (3.29) would be a Laplace convolution where we use the fact that L −1 L equals the identity. By definition of Laplace convolution, The above result is a Duhamel integral, which can be inverse-Fourier transformed into This completes the proof of (3.43).

Building Integral Formalism
In this subsection, a system of nonlinear integral equations is constructed from the Lemmas developed in the above. We start with some preliminaries of background material of function space.
Proposition 2 (Banach space). Let B T be the linear space of continuous and bounded real- is a Banach space.
With the basis of Banach spaces in the above proposition, relevant operators may be defined in association with Lemma 4, as following.
for b ∈ B T , which are also linear.
Lemma 5 also motivates us to define more operators, as follows.

Definition 3 (linear operator)
The two transformations, U ( j) : B 0 → B T , for j = 1, 2, such that for b 0 ∈ B 0 , are linear integral operators. In an analogous manner, for b ∈ B T , U ( j) : are linear transformations as well.
Owing to Definitions 2 and 3, the nonlinear partial differential equations, (1.1) and (1.2), may be recast in the form of two coupled integral operator equations via (3.30) and (3.39). They can be, in the current study, regarded as an integral formalism corresponding to the standard Boussinesq's equations for variable water depth of (1.1) and (1.2).

Theorem 1 (integral formalism)
Let h 0 > 0 denote a pseudo-water depth parameter and let ϕ and ψ be the forcing functions of (2.3) and (2.4), respectively. Then, given h 0 > 0, wave elevation η and depth averaged velocity u satisfy

Constructing Integral Operator Equations
As mentioned above, Theorem 1, presented in Sect. 3.3, provides an integral formalism, which has a form of the coupled operator equations. It is regarded as a theoretical framework which plays fundamental role in constructing integral (operator) equations corresponding to (1.1) and (1.2). To this end, we need more operators which involve partial derivatives.

Definition 4 The transformations, denoted by
for first partials and for mixed partials. (For their explicit estimates, see Appendix A.1) Notice that the domains of the operators defined in the above Definition are B 0 , whereas those of operators in the next Definition will be B T .

Definition 5 By
we denote the transformations from B T to itself such that  We next, for systematic analysis, introduce new variables χ j , j = 1, 2, . . . , 7 such that The following properties are then immediately checked: from (4.5).

Proposition 3
The forcing terms ϕ and ψ in (2.3) and (2.4), respectively, can be expressed in terms of the new variable χ in (4.5) (4.10) (4.9) and (4.10) indicate that ϕ depends on wave elevation η, average velocity u and their spatial derivatives η x , u x , while ψ depends on average velocity u and their mixed derivatives u x , u t , u xt , u x xt . This implies and proves both (4.7) and (4.8) from (4.5) and (4.6).
In the next Definition, a h 0 -parameterized nonlinear operator shall be defined with the combination of the linear operators defined in the above Definitions, where the nonlinearity comes from the nonlinear functional forms of ϕ and ψ in Proposition 3.
Definition 6 (nonlinear operator) Given a pseudo-water depth parameter h 0 > 0 and initial wave and velocity profiles η 1 , u 1 ∈ B 0 , T h 0 is defined as a nonlinear transformation (or mapping) from B 7 T to itself, such that Here, ϕ and ψ are nonlinear scalar functions on the variable χ in (4.5), as represented in (4.9) and (4.10), respectively.
We are now in a position to arrive at a system of integral operator equations, which is equivalent to (1.1) and (1.2). , (1.1) and (1.2). Then, the χ ∈ B 7 T in (4.5), uniquely determined according to (4.6), is invariant under the (nonlinear) mapping T h 0 in (4.11) for a pseudo-water depth parameter h 0 > 0, i.e.,

Remark 3 Being referred to as a fixed point of the mapping T h 0 , χ ∈ B 7
T in (4.5) is invariant under T h 0 , no matter what the pseudo-parameters, h 0 > 0, are.

Remark 4
As a consequence of Remark 3, the solution of (4.12) does not depend on the choice of the pseudo-water depth parameters h 0 > 0. Moreover, (4.12) is known to have the stability properties, which follows from the fact that (4.12) is considered to be nonlinear integral equations of Hammerstein type [4].

Application of Integral Equation Formalism
The purpose of this section is to examine the classical problem (in shallow water hydrodynamics) of time-evolution of a solitary wave on a sloping beach, being an application of the integral equation formalism constructed in the previous section. The problem, in practice, is usually observed; if relatively gentle waves approach a shore, their crests would behave in some ways like separate waves so that they sometimes look like solitary waves [1].

A Solitary Wave on a Sloping Beach
This subsection gives a brief review of a solitary water wave and its interaction with a beach, together with discussion of their non-dimensionalization.
A solitary wave, also called a soliton, maintains its shape while it propagates at a constant speed, which is caused by a balance of nonlinearity and dispersion. It is of the solution of a widespread class of (weakly) nonlinear dispersive partial differential equations, which often arise in a wide variety of real physical systems (e.g., water waves, elastic waves and so on). In water waves, according to Lamb (1932, section 252) [6], a solitary wave of amplitude a 0 in water of uniform depth d 0 , which moves to left at speed c, is expressed as The propagation speed c in (5.1) is expressed as where the binomial expansion has been used to approximate the speed with assumption of a 0 /d 0 1. As a representative (or characteristic) length in the problem, henceforth, we select the water depth d 0 (without loss of generality), which is a well-defined proper reference length to measure each of the variables in (5.1) and (5.2). Then, the dimensionless form of the governing equations, (1.1) and (1.2), become [1] (see also Appendix B) with appropriate dimensionless (flow) variables, denoted bŷ Moreover, (5.1) and (5.2) may be expressed in dimensionless form respectively, withâ 0 andĉ dimensionless quantities [consistent with (5.5)], The solitary wave solution mentioned above enables us to find the initial conditions for wave propagation of a solitary wave on a gently sloping beach (See Fig. 1). Here, the variable water depth is simply specified as h(x) = αx for small α > 0 (5.9) with water in the region x > 0: i.e., a sloping beach of uniform small slope α is assumed for two-dimensional water waves which have the crests parallel to the shore line. From (5.1), we can impose the initial wave profile η solit of the solitary wave over the uniform sloping beach Fig. 1 A solitary wave over a gently uniform sloping beach of slope α of (5.9), in the form (5.10) so that it can describe a solitary wave in water of depth d 0 = αx (5.11) with its crest at x = d 0 /α. And, associated with η solit in (5.10), the initial depth averaged velocity u solit (x) is found from the continuity Eq. (1.1) to be [1] The non-dimensionalization of (5.10) and (5.12) leads tô respectively [1]. Here,η solit andû solit denote the initial dimensionless profiles of free surface elevation and depth averaged velocity, respectively, i.e.,η solit = η solit /d 0 ,û solit = u solit /(gd 0 ) 1/2 .

A Solitary Wave in Water of Constant Depth
Before detailed discussion of a solitary wave on a sloping beach, it seems instructive to investigate the (original) standard Boussinesq's equations (1.1) and (1.2) when water depth is constant, say d 0 ; (1.1) and (1.2) reduce to and respectively. The above equations can admit the soliton solution (5.1) if they are slightly modified as discussed the following Lemma.

Lemma 6 Let the momentum equation (5.16) be changed into
with a residual function R expressed in terms of the hyperbolic tangent function where the argument X is denoted by Then, the pair (5.15) and (5.17) allows exact solutions for a solitary wave in the form where x 0 is constant and c is given by (5.2).
Proof We begin with showing that (5.20) and (5.21) satisfy the mass conservation (5.15). Inserting η exact in (5.20) into (5.15) and subsequently integrating the resultant with respect to x, we have Re-arranging the above yields where the last is identical with the u exact (x, t) in (5.21), i.e., This implies that (5.20) and (5.21) are solutions for the mass conservation (5.15).
We next show that (5.20) and (5.21) satisfy the modified momentum conservation (5.17). To this end, we need some calculations; Using the above results, we substitute (5.20) and (5.21) into (5.17), resulting in (5.18). Therefore, conclusively, we have proved the Lemma. Figure 2 shows a typical plot of the residual function R depending upon amplitude a 0 , which indicates that the residual reduces as a 0 decreases. In what follows, gravity acceleration will be approximated as g ≈ 9.80665 m/s 2 .

Remark 5
The forcing terms ϕ in (2.3) and ψ in (2.4) for the standard Boussinesq's equations may be changed as according to the modified system of (5.15) and (5.17).
We would like to simulate the exact solitary wave solutions (5.20) and (5.21). The simulation can be accomplished by solving the integral operator Eq. (4.12) but with ϕ and ψ in (5.22) and (5.23), respectively, instead of solving the (modified) partial differential equations (5.15) and (5.17).
The numerical solution for (4.12) may be found by the successive approximation [3] (or Banach fixed point theorem, see [7]). That is, bearing in mind that (4.12) is of a fixed form and the product space B 7 T is a Banach space, we seek a fixed point solution χ in the solution The output of (5.24) gives not only both of free surface elevation η and depth averaged velocity u but also their (high) derivatives as a bonus via (4.5), if it converges.
Regarding the numerical calculation of (5.24), it should be pointed out that all the integrals involving T h 0 are not singular: e.g., the kernel of the operator Y (4) in (3.48) has two singularities of simple poles at k = ± i √ 3/ h 0 , which are off the real line R, because a pseudo-parameter h 0 , artificially introduced in Sect. 2, is real positive. Due to this fact, it suffices to perform conventional regular numerical integrations (such as the trapezoidal or Simpson's rule) for calculating the integrals (and thus the operator T h 0 ), which is demanded for running the recurrence relation (5.24). Figure 3 represents a flowchart for the implementation steps that are taken to implement the iterative algorithm (5.24).

Remark 6
Returning to the numerical simulation, we will iterate (5.24), whereby the initial value problem consisting of the two governing equations of (5.15) and (5.17) and the two initial conditions is to be solved; a 0 = 0.1m and d 0 = 1.0 m are taken. Here, the integrals appearing in the operator T h 0 in (4.11) are evaluated with the use of the usual Simpson rule for space and trapezoidal numerical integration for time, respectively, where the pseudo-parameter is selected as h 0 = 1.5 m. We discretize the time interval 0 sec < t < 1 · √ d 0 /g sec ≈ 0.3193 sec (or 0 <t < 1) into 40 equal length-parts so that its increment can be equal to t ≈ 0.008 sec (or t = 0.025). In an analogous way, we divide the spatial interval 10 m < x < 50 m (or 10 <x < 50) into 800 equally spaced panels to produce x = 0.05 m (or x = 0.05).  with L2 norm taken over the x only at end time. It is worth discussing that the smaller h 0 /d 0 would contribute to quicker convergence and more accuracy in both cases of Err η and Err u ; numerical convergence is almost attained within n ≈ 10. Thus, the case h 0 /d 0 = 1 corresponds to the most accurate result, where we find it interesting that the highest-derivative term in (5.28) vanishes, because d 2 0 −h 2 0 = 0; the more accurate results can be found in Table 1 with a finer mesh. Mathematically, the (numerically) converged solutions can be considered as    [7]. It would be instructive to calculate the rate of convergence of the iteration (5.24). Table 1 lists values for the rates of convergence [8] , (5.29) with the symbol | x, t , which means evaluation with the use of increments of x and t.
Here, the results are obtained for various x but with fixed t = 0.0125, showing that the convergence rates diminish with a finer step size. Analogously, the two rates 30) are plotted in Table 2, where t varies while x is fixed (i.e., x = 0.02). Here again, we observe that the two rates become smaller with a finer time step size.

Simulating the Solitary Wave Transformation
In this subsection, we shall simulate the propagating solitary wave over the sloping beach (governed by the original partial differential equations (1.1) and (1.2)) and investigate how it interacts with a beach. The interaction can be determined by solving the integral operator equation (4.12) by using the recurrence relation (5.24). Here, the effect of the pseudo-water depth parameter on the numerical solutions will be examined as well.  Figure 9 shows the numerical convergence behavior of (5.24) for α = 1/30, where a 0 = 0.1 m, d 0 = 1.0 m and h 0 = 1.5m are selected; the time interval is taken as 0 sec < t < 15 · √ d 0 /g sec ≈ 4.79 sec (or 0 <t < 15) with its increment equal to t ≈ 0.04 sec (or t = 1/8; 120 equal length-parts) and, similarly, the spatial interval as − 20 m < x < 60 m (or − 20 <x < 60) with its increment x = 1/4 m (or x = 1/4; 320 equally spaced panels). We can see a reasonable solution calculated just by the first iteration for small t and an almost converged solution (valid for the whole dimensionless time interval) found at n = 5. Especially, Fig. 10 depicts, at fixed dimensionless time t ≈ 4.79 sec (or t = 15), the convergence behavior for the dimensionless solitary wave transformation.  Figure 11 illustrates the time evolution of the dimensionless solitary wave transformation resulting from the interaction between wave and beach, which is in a good agreement with the Peregrine's result [1]. It is possible to observe that the wave steepens and gradually grows in amplitude. Further, the wave profile at initial stage looks symmetric, however, the wave profile tends to become asymmetric. In fact, the nonlinear (interaction) feature of solitary wave transformation was firstly explored by Peregrine [1], however, who made a finite difference approximation of his model of (1.1) and (1.2). Madsen & Mei [9] examined the same problem by using a system of partial differential equations (which is essentially equivalent to Peregrine's model) derived by Mei & Le Méhauté [10]. On the other hand, Chan & Street [11] directly approximated the Navier-Stokes equations (instead of using Peregrine's model) with finite difference equations to observe the solitary wave transformation.
We will examine one more example, which corresponds to a beach of different uniform slope, α = 1/20 (greater than the previous one). The convergence behavior of (5.24) is exhibited in Fig. 12, where the spatial discretization condition is the same as that of Fig. 9, but the time discretization condition is as follows; 0 sec < t < 10 · √ d 0 /g sec ≈ 3.19 sec (or 0 <t < 10) and its increment t ≈ 0.04 sec (or t = 1/8; 80 equal length-parts). It is seen that the convergence nature of the iterative solutions in Fig. 12 appears to be essentially similar to that of Fig. 9; so, Fig. 13, illustrating the convergence nature at fixed dimensionless time, looks alike Fig. 10. The time evolution of the deforming solitary wave is pictured in Fig. 14, presenting the same nonlinear features as that of Fig. 11 (e.g., the wave steepens and becomes asymmetric, growing in amplitude).  Figure 15 represents the wave (maximum) height for three cases of h 0 /d 0 (i.e., h 0 /d 0 = 1.0, 1.5, 2.0), compared with Peregrine's result for α = 1/20 [1]. Here, the full line indicates Green's law, derived from the linearized shallow water wave equation (which appears to present a reasonable approximation). It is found that there is an excellent agreement between ours and Peregrine's [1] especially when h 0 /d 0 = 1.0 (this phenomenon looks similar to Fig. 15 Wave-amplitude variation for α = 1/20;â means the (local) maximum height of the wave of Fig. 7 that of Sect. 5.2, where the most accurate result was found when h 0 /d 0 = 1.0 as shown in Fig. 7). This implies that the accuracy of (5.24) may depend on the artificially introduced dimensionless pseudo-parameter h 0 /d 0 (as is clearly confirmed in Fig. 15 as well as Fig. 7); in other words, it may be a crucial parameter which controls the solution accuracy [3].

Concluding Remarks
The main purpose of this paper is to construct new regular integral equations corresponding to the standard Boussinesq's equations for variable water depth (or Peregrine's model), but which should be different from existing integral equation formulations such as the singular integral equations from the boundary element method.
For that, we have introduced a pseudo-water depth parameter to convert the nonlinear partial differential equations of Peregrine's model into a one-parameter family of integral formalisms (in the form of operator equations) by using appropriate (Laplace and Fourier) integral transform techniques. This results in formulating a one-parameter family of the coupled nonlinear integral equations of second kind, being different from existing integral equation formulations. This implies that we have constructed new regular integral equations which correspond to Peregrine's model.
As an application, the regular integral equations constructed above immediately yield a new iterative formula for integrating the standard Boussinesq's equations for variable water depth through the fixed point iteration. The formula makes it possible to simulate a solitary wave moving toward a beach. The simulated results have been compared with the previously known (numerical) solutions and the agreement is excellent, where various pseudo-water depth parameters have been employed. An interesting nonlinear feature is observed on the transformation of a solitary wave over variable depth; the solitary wave gradually steepens and its amplitude grows with height, as expected.
It may be worth mentioning that there are a few roles of the pseudo-water depth parameter merged into the formula. First, it acts as a bridge between the original partial differential equations and a one-parameter family of the coupled nonlinear integral equations. Second, it may affect the solution accuracy so that it can be regarded as an accuracy control parameter, as was emphasized in the discussion in the preceding section where the effect of the pseudowater depth parameter on wave-amplitude variation was compared to both Peregrine's result and Green's law. Finally, in the presence of the pseudo-water depth parameter, off the real axis is the singularities of simple poles of the kernels of the integral operators used for the formula. As a consequence, the derived formula becomes singularity-free.
t (b 0 ) and U (2) t (b 0 ) can be calculated as The results of U x xt (b 0 ), j = 1, 2, in (4.2) are given as follows Here, the rule of Leibniz differentiation has been applied.