Abstract
We establish existence and uniqueness of solutions to the Cauchy problem associated with a new one-dimensional weakly-nonlinear, weakly-dispersive system which arises as an asymptotical approximation of the full potential theory equations for modelling propagation of small amplitude water waves on the surface of a shallow channel with variable depth, taking into account the effect of surface tension. Furthermore, numerical schemes of spectral type are introduced for approximating the evolution in time of solutions of this system and its travelling wave solutions, in both the periodic and nonperiodic case.
1. Introduction
In this paper we study the propagation of water waves on the surface of a shallow channel with variable depth, considering the effect of surface tension. To describe this phenomenon we will derive a new water wave model from Euler’s equations (in dimensionless variables) for an inviscid, incompressible liquid bounded above by a free surface and bounded below by an impermeable bottom topography [1]:with the nonlinear free surface conditionsat . Here denotes the potential velocity and the wave elevation measured with respect to the undisturbed free surface . The dimensionless parameters and are small positive real numbers which measure the strength of nonlinear and dispersive effects, respectively. The parameter measures the ratio inhomogeneities/wavelength and the parameter is associated with the surface tension. The Neumann condition at the impermeable bottom isThe bottom topography is described by , where We will see in Section 3 that system (1)–(3) is equivalent to leading order provided that , to the weakly-nonlinear, weakly-dispersive systemThe variables and are the space and time coordinates, respectively, with denoting the elevation of the free surface, and represents the potential fluid velocity measured at the channel’s bottom. The metric coefficient is related to the channel’s bottom and it is defined aswhere is the coordinate transformation used in the derivation of system (5) to map the original physical channel with variable depth onto a strip in the complex plane. This strategy of change of variable was introduced by Hamilton [2] and later used successfully by Nachbin in [3] to study wave propagation over a channel with a highly variable topography. Observe from (6) that the coefficient is infinitely differentiable although the function describing the channel’s bottom is not smooth. In case of constant depth, the coefficient is identically one and system (5) reduces to a system derived by Quintero and Montes in [4].
Wave-topography interaction has been the subject of considerable mathematical research [5–18]. The physical applications range from coastal surface waves [19] to atmospheric flows over mountain ranges [20, 21]. In particular, the interaction of waves with fine features of the topography is of great interest. As pointed out in the introduction to the Orography proceedings [21] of the European Centre for Medium-Range Weather Forecasts (ECMWF), “the representation … of subgrid-scale orographic processes is recognized as crucial to numerical weather prediction at all time ranges.” In the atmospheric literature, orography implies mountain ranges [20].
In previous works, some weakly-nonlinear, weakly-dispersive models have been developed to describe the interaction of a long pulse with small amplitude that propagates on the surface of a channel with a variable bottom [2, 3, 5, 7–9, 22–24]. However, these models either neglect the effect of surface tension on the free surface where the wave propagates or are not applicable to bottoms described by a discontinuous/nondifferentiable function or the wave elevation is removed in their physical derivation. For instance, Milewski [23] derived a bidirectional scalar Benney-Luke type equation in terms of the potential velocity which includes the surface tension effect and the influence of a variable bottom. However, the asymptotical derivation of this model implies eliminating the wave elevation and consequently neglecting several second-order terms in the parameters and .
Some of the features of new formulation (5) are the following:(i)It can be applied to study wave propagation over a shallow channel with a discontinuous or nondifferentiable bottom, provided that the bottom’s fluctuations satisfy . This is a consequence of introducing the conformal mapping mentioned above. Observe that all coefficients of the reduced equations (5) in the new coordinate system result in being infinitely differentiable since they depend on the smooth function .(ii)One-dimensional system (5) models bidirectional waves and it incorporates the simultaneous effects of surface tension and variable depth upon the shape of a water wave that propagates on the surface of an irregular shallow channel.(iii)Furthermore, in the derivation of system (5), we do not eliminate the wave elevation (which prevents us from neglecting additional terms of order ), as done, for example, in [23]. Therefore, the new formulation (5) is expected to be a more accurate approximation of the full potential theory equations.
In the present paper, we establish existence and uniqueness of a solution to Cauchy problem (5) using classical semigroup theory and Banach’s fixed point principle. In [4] the well-posedness of system (5) was analyzed but only in the case that . In second place, we formulate a Galerkin-spectral numerical scheme of spectral accuracy in space to approximate the solutions of system (5) based on the Fourier basis and using an implicit-explicit (IMEX) second-order strategy for time stepping. This type of temporal discretization is described, for instance, in [25], and it has been used in conjunction with spectral methods [26, 27] for the time integration of spatially discretized PDEs of diffusion-convection type. IMEX schemes have also been successfully applied to the incompressible Navier-Stokes equations [28] and in environmental modelling studies [29]. On the other hand, we develop a numerical solver to compute travelling wave solutions of system (5) by using a Fourier-collocation strategy combined with a Newton-type iterative procedure. We also indicate how to determine appropriate starting points for this iterative process in order to achieve convergence. Existence results on travelling wave solutions (in the periodic and nonperiodic case) of system (5) with have also been established in [4]. Travelling wave solutions exist as a consequence of a balance between nonlinear and dispersive effects present in a system; these waves travel with a constant speed without any temporal evolution in shape or size when the frame of reference moves with the same speed of the wave. In the last decades, the study of travelling waves has grown enormously because they appear in several and varied fields of application, such as fluid mechanics, optics, acoustics, oceanography, and astronomy. Thus, to determine existence and properties of such type of solutions is a fundamental problem in the theory of ordinary and partial differential equations of great interest for both pure and applied mathematicians.
The rest of this paper is organized as follows. In Section 2, we introduce the functional spaces and notation to be employed in the paper. In Section 3, we present in detail the derivation of model (5) starting from the potential theory equations including the surface tension effect. In Section 4, we discuss existence and uniqueness of a solution of the Cauchy problem associated with system (5) by using Banach’s fixed point principle and semigroup theory. In Section 5, we introduce the numerical schemes for approximating the evolution of a solution of system (5) and computing their travelling wave solutions (periodic and nonperiodic cases). Section 6 presents a set of numerical simulations to check the accuracy of the numerical schemes developed in the paper. Furthermore, some numerical experiments are included to illustrate the interaction between the wave-topography effects and surface tension. Finally in Section 7 are the conclusions of the paper.
2. Preliminaries
To analyze existence of solutions of problem (5), we will use the standard notation. For , we will denote by (or simply ) the Banach space of measurable functions in such that if and , if . We define the norm in for byand in by . is a Hilbert space for the scalar productWe set . For function , the Fourier transform is defined asand the inverse Fourier transform is defined byWe will also denote by (or ) and (or ) the extensions of these operators to . The convolution of two functions, , is defined asWe recall that . For , we define the Sobolev space (sometimes written for simplicity as ) as the completion of the Schwartz space (rapidly decaying functions) defined as with respect to the normFor simplicity, we also denote this norm by . The inner product in is defined asThe product norm in the space is defined by for . Sometimes, we will also use the equivalent product normFor , we will denote by the space of continuous functions , that is, the space of continuous functions , , with the supremum norm and the product normfor .
3. Governing Equations
We start by presenting the potential theory formulation for Euler’s equations (in dimensionless variables) for an inviscid, incompressible liquid bounded above by a free surface and bounded below by an impermeable bottom topography and including the effect of surface tension [1]:with the nonlinear free surface conditionsat . Here denotes the potential velocity and the wave elevation measured with respect to the undisturbed free surface . The dimensionless parameters and measure the strength of nonlinear and dispersive effects, respectively. The parameter measures the ratio inhomogeneities/wavelength and is related to the surface tension effects. The Neumann condition at the impermeable bottom isThe bottom topography is described by , whereThe bottom profile is described by the (possibly rapidly varying) function . The topography is rapidly varying when . The scale represents the total length of the irregular section of the coast. The undisturbed depth is given by and the topography can be of large amplitude provided that . The fluctuations are not assumed to be small, nor continuous, nor slowly varying.
Let us consider a symmetric flow domain by reflecting the original one about the undisturbed free surface (cf. Figure 1). This domain is denoted by where and can be considered as the conformal image of the strip where with . Then with and , a pair of harmonic functions on . Following the strategy suggested by Hamilton in [2] and Nachbin [3] within the weakly-nonlinear, weakly-dispersive regime (, ) and using the relationships whereand the variable free surface coefficient defined asthe potential theory equations can be approximated to order in the orthogonal curvilinear coordinates with by the equationwith conditions at the free surface ,and condition at the channel’s bottomObserve that the change of variables lets the origin of the curvilinear coordinate system at the bottom. The Jacobian for the coordinate transformation is represented by , and corresponds to the position of the free surface in the curvilinear coordinate system. By performing an asymptotic simplification as in [1] (page 464) through a power series expansion in terms of the dispersion parameter near the bottom of the channel in the forminto (26)–(28) for being small, we find that free surface conditions (27) can be approximated to order , by the equationsHere denotes the potential fluid velocity at the channel’s bottom . We point out that (30) and (31) with (no surface tension) correspond to those derived in [3] (bottom of page 915 and top of page 916). In the derivation of (30)-(31), we used the relationshipwhich means that, at leading order, the Jacobian of the conformal coordinate transformation is time independent.
Let us introduce the new variable . Thus observe that from (30) This relationship allows us to change the form of dispersive terms in the equations above. In particular, by using the decomposition in (30), we obtain that system (30)-(31) can be approximated to order by the following equations:We point out that the system above is applicable to modelling of propagation of water waves over an arbitrary rapidly varying depth. In case of a slowly varying channel’s topography and with being a stationary random process with standard deviation and correlation length of order one, the metric coefficient can be expanded as (see [9])Thusand system (35) can be approximated to order byNote that the coefficient is smooth even when the function describing the bottom is discontinuous or nondifferentiable. The function is time independent and becomes identically one for a channel with constant depth. In this case, system (38)-(39) reduces to that studied in [4]. Moreover, in applications, this coefficient is bounded and . These properties will be important to obtain existence and uniqueness results for system (38)-(39). We also point out that the function actually depends on the dispersion parameter . For this reason, it will be denoted by whenever we need to emphasize this dependence.
4. Existence and Uniqueness
System (38)-(39) can be written aswhereSystem (40) is supplemented with the initial conditionsTaking Fourier transform in the spatial variable in system (40), we havewhere
4.1. Analysis of the Linear Semigroup
Consider system (40) with ,which has unique solution in the formwhereBy using the inverse Fourier transform in (46), we have thatwherewith
Theorem 1. The family of linear operators is a -semigroup in . Furthermore is its infinitesimal generator in .
Proof. Let . ThenHereafter will denote a generic constant and we recall that . Therefore , , is a family of continuous linear operators in . On the other hand, it is easy to see that , , . Finally Observe that , as and , , as . By using Lebesgue’s dominated convergence theorem, we can conclude thatas . We conclude that , , is a strongly continuous semigroup in .
On the other hand, To see that is the infinitesimal generator of the semigroup , , observe that, for , But by virtue of and using again Lebesgue’s dominated convergence theorem, we get that
4.2. Analysis of the Nonlinear Term
Theorem 2. Let . The application maps on itself andFurthermore,where is a constant and , .
Proof. In the first place, let : Taking into account the inequalities valid for (applying Corollary in [31]), and using the fact that the function is bounded, we get that On the other hand, for , we have that
4.3. Local Existence and Uniqueness
Theorem 3. Let and . Then there exists and unique , which satisfies the integral equation
Proof. Let be fixed constants and consider the nonlinear operatordefined in the complete metric spaceLet us prove that if , then . Indeed if , then We remark that if , thenGiven that is a -semigroup, then as . To analyze the second expression, let us suppose that . We obtain that Observe that since as , , and using Lebesgue’s dominated convergence theorem, we get that as .
On the other hand, due to Theorems 1 and 2, we arrive at Therefore we get that as . This means that .
Let us prove existence of small time such that . Let . For any , for being small enough. Therefore .
Finally, let us show that time exists so that the operator is a contraction on . Let , . Then using Theorem 2, we arrive at We see that it is possible to select time such thatso that the operator is a contraction on the closed ball . The fixed point principle guarantees the existence of unique solution of integral equation (64). Finally, uniqueness of this solution in the metric space can be established via Gronwall’s lemma.
5. Numerical Schemes
In this section we describe the numerical scheme we propose to approximate solutions of system (38)-(39). In the first place, it is convenient to rewrite it assubject to the initial conditions , , where we introduced the new variable . In the numerical solver to be introduced, the computational domain is discretized by equidistant points, with spacing , and the unknowns and are expanded as truncated Fourier series in space with time-dependent coefficients:with The time-dependent coefficients , , are calculated by means of the equationand analogously for . Substituting these expressions into (75) and projecting the resulting equations with respect to the -orthonormal basis and the inner productit follows thatwhereand denotes the operatorWhen the period is taken large enough, this numerical scheme can be applied to approximate the rapidly decaying solutions of system (75) on the entire real line . This technique was used successfully in [9].
5.1. Temporal Discretization
Note that (80) can be seen as a system of ordinary differential equations where the unknowns are the time-dependent Fourier coefficients of the solutions. To solve it, we use an implicit-explicit scheme (IMEX) in the formHere , denote the approximations of the unknowns , , respectively, at time , where is the time step of the method, and denote the approximations of the functions evaluated at time . Similarly, denote the approximations to the Fourier transforms of the functions and , respectively, with respect to the variable , evaluated at time . The numerical approach adopted for solving system (75) ensures the scheme results to be linearly unconditionally stable which can be easily verified. Further, observe that the linear dispersive terms are approximated by using an implicit strategy, in contrast to the nonlinear terms and terms where the variable coefficient is present which are treated in explicit form. Implicit-explicit schemes (IMEX) were already applied in [25] for scalar dispersive evolution equations. The main advantage of the numerical scheme described here is that at each time step we can solve explicitly the approximations , , , of the Fourier coefficients of the unknowns and from (83), without using implicit Newton-type iterations. Thus, the scheme results in being cheap and its computer implementation is easier.
In the scheme, the spatial derivatives and are computed by exact differentiation of the truncated Fourier series. For instance,The numerical calculations presented in this paper were carried out in double precision by using MATLAB R2012b on a Mac platform. The Fourier-type integral appearing in the operator (see (82)) is approximated through the well-known Fast Fourier Transform (FFT) routine.
5.2. Approximating Travelling Wave Solutions
In this section we are interested in computing solutions of system (75) over a channel with flat bottom (i.e., when the metric coefficient is equivalent to 1) of the form which are named as travelling wave solutions. The parameter is called the wave velocity.
After dropping the tildes, travelling wave solution with speed of system (75) must satisfy the following equations:In first place, we are interested in finding approximations to even periodic travelling wave solutions with period , of system (75). Thus let us introduce truncated cosine expansions for and , whereand analogous expressions for . By substituting expressions (87) into (86), evaluating them at the collocation pointswe obtain a system of nonlinear equations in the formwhere the coefficients are the unknowns. Nonlinear system (90) can be solved by Newton’s iteration. Computation of the cosine series in (87) and the integrals in (88) is performed using the FFT (Fast Fourier Transform) algorithm. The Jacobian of the vector field is approximated by the second-order accurate formulawhere and . We stop Newton’s iteration when the relative error between two successive approximations and the value of the vector field are smaller than .
The starting point for Newton’s procedure in the periodic frame is taken asIn second place, we are also interested in approximating solitary wave solutions of system (75), that is, when the functions and their derivatives decay to zero at . In this case, we can also approximate by a truncated Fourier series as in (87) using a large enough length . To derive an appropriate initial point for Newton’s iteration, we compute an approximate solitary wave solution of system (75) with and being small. To accomplish this, we will use original system (38)-(39) in the variables . Let us observe that from (39)Substituting this into (38) and neglecting second-order terms in , we derive the Benney-Luke type equation for the potential :We look for a solitary wave solution of (94) in the formTherefore, abandoning the tildes, the function must satisfy the following equation:Integrating the equation above and using the fact that and their derivatives decay to zero at , we arrive atMultiplying the previous equation by and integrating again, we obtain thatFurthermore, letting , we find from (93)Assume the solution form of to bewhere are constants. Then, by replacing this into (98), we get that the constants are given byNow substituting into (99), we arrive at the following approximation of the wave elevation:
6. Description of the Numerical Experiments
In this section we compute some approximations to solutions of system (75) using the numerical schemes described in the previous section. In first place, some travelling wave solutions are approximated by using the numerical scheme discussed in Section 5.2 and then these solutions are checked through computer simulations conducted with numerical solver (83). Finally, we apply it to compute the evolution of a gaussian pulse over a channel with variable depth taking into account the surface tension effect.
6.1. Periodic Travelling Wave Solutions
Using Newton’s iteration as explained in Section 5.2, we develop some numerical simulations to compute periodic travelling wave solutions of system (75). The results are shown in Figures 2, 3, 4, and 5 for different values of the modelling parameters and wave speed . The initial step for Newton’s scheme is taken as in (92). The number of FFT point is and we recall that here the channel has constant depth; that is, . In Figure 6 we check the accuracy of the periodic travelling wave solutions computed by comparing the prediction of numerical scheme (83) at time (using and ), with the profile in Figure 2 shifted at a distance of to the right. Observe that the profiles coincide with good accuracy and they propagate with the expected wave velocity. The difference between the two profiles is of order . Similar results were obtained for the other periodic travelling wave solutions shown in Figures 3, 4, and 5.
6.2. Solitary Wave Solutions
Numerical experiments in the case of solitary wave solutions are shown in Figures 7, 8, and 9 for different values of the modelling parameters and wave speed . The initial step for Newton’s scheme is taken as in (100)-(102) centered at the position . The number of FFT point is and we recall that here the channel has constant depth; that is, . In this nonperiodic scenario, the spatial computational domain is the interval which is large enough so that the pulses do not reach the computational boundaries within the time interval modeled. In Figure 10 we check the accuracy of the periodic travelling wave solutions computed by comparing the prediction of the numerical scheme (83) at time (using and ), with the profile in Figure 2 shifted at a distance of to the right. As in the experiments in the previous section, we see that the profiles coincide with good accuracy and they propagate with the expected wave velocity. The difference between the two profiles is of order . Similar results were obtained for the other periodic travelling wave solutions shown in Figures 8 and 9.
6.3. Variable Depth
Finally in this section we include some numerical simulations to illustrate the dynamics of incoming Gaussian pulses in the form governed by (75), under the simultaneous effects of surface tension (), dispersion (), and channel’s topography (coefficient ). We further remark that the incoming pulses are located at two units to the left of the irregular part of the channel which covers the interval . In order to focus on these three phenomena, we only consider the linear case; that is, .
In Figure 11 is displayed the output of numerical scheme (83) for the value of the dispersion parameter and without surface tension effect, that is, , for a constant depth channel. The numerical parameters are and and the computational domain is the interval which again is large enough so that the pulses do not reach the computational boundaries. We point out that, as a consequence of dispersion (which produces a difference between the phase velocity of the Fourier modes contained in the pulse), the incident pulses progressively develop an oscillatory tail that propagates to the left, while the leading wave propagates to the right. We contrast this experiment to the dynamics of a pulse under the influence of surface tension which is presented in Figure 12. Note that, in this case, both the oscillatory tail and the wave front propagate to the right side of the computational domain. We conclude that the surface tension effect alters the dispersive characteristics of model (75).
We can anticipate the influence of the additional third-order term (due to the surface tension effects) on the solutions of system (75) (with and flat channel ) by analyzing its linear dispersion relation given byWe recall that this expression arises by studying the propagation of solutions of system (75) in the form , , where are constants.
In Figure 13, the phase velocities (with ) and (no surface tension) are compared as functions of ( is called the wavenumber). Observe that, in the presence of surface tension, the phase speed is larger than 1, and it approximates to the value as tends to infinity. In contrast, without surface tension effects, the phase speed is lower than 1 and Fourier components with high wavenumber travel more slowly than the Fourier modes with small wavenumber. These observations are in accordance with the results presented in Figures 11 and 12. We further point out that the fact that the phase velocity in system (75) is bounded by for all has an advantage in improving the stability of numerical solvers formulated to approximate solutions of this system. In particular, the properties of numerical solver (83) employed in the present paper, such as the convergence, stability, and error estimates, will be studied in a future research.
After analyzing the flat bottom case, we consider a numerical simulation in Figure 15 where the channel’s topography is variable. The topography considered is shown in Figure 14 together with the metric coefficient computed using the Schwarz-Christoffel Toolbox designed for MATLAB by Driscoll [30]. Note that the irregular part covers the interval (10 units long). A nonnegligible value of surface tension is incorporated into model equations (75). Again the level of dispersion is taken as . The numerical parameters are the same as in the previous computer simulation. In addition to the effect of the surface tension on the direction of propagation of the oscillatory tail, we now observe irregular fluctuations (coda) developing behind the wave front and propagating to the left end of the computational domain. This tail signal is due to the multiple reflections of the incoming pulses produced by every variation of the channel’s depth. These phenomena have been studied in previous works using other water wave models that neglect the surface tension effect [9–11]. In the present paper, our purpose was to initiate the study of these wave phenomena and numerical scheme (83) and system (75) have proved to be valuable tools in developing this research in a future work.
To see in more detail the effect of surface tension on the wave-topography interaction, in Figure 16, we show the evolution of the Gaussian pulse on the surface of a channel with a large irregular bottom covering the interval (10 times larger than in the previous numerical simulation) in the case of (no surface tension) and . As in the previous experiments, the physical depth is described by a linear piecewise function and the coefficient is computed again using Driscoll’s package for MATLAB. The numerical parameters for this computer simulation are , , , and the computational domain is the interval . Observe that the pulse without the presence of surface tension appears to be delayed with respect to the pulse propagating under the influence of a nonnegligible value of surface tension and the long channel’s topography introduced. These phenomena and further research using full model equations (35) for the case of a rapidly varying bottom (i.e., ) will be studied in detail in a future work. We remark that numerical scheme (83) could be adapted to solve full system (35).
7. Conclusions
In this paper we introduced a new model for describing the propagation of a small amplitude water wave over the surface of a shallow channel with a possibly varying depth. It was derived formally as an asymptotical approximation of the full potential theory equations in the weakly-nonlinear, weakly-dispersive regime. Existence and uniqueness of a solution to the initial value problem associated with (5) were also studied using Banach’s fixed point argument and the classical semigroup theory in appropriate Sobolev spaces. The main features of model (5) considered in this work are the following:(i)It can be applied to study wave propagation over a shallow channel with a discontinuous or nondifferentiable bottom provided that the bottom’s fluctuations satisfy . Observe that in the derivation of system (5) the function was not assumed to be small nor continuous since the conformal mapping has a smoothing effect on the topography. Note that the effective variable coefficient that appears in final system (5) is defined as a convolution between the function describing the bottom’s fluctuations and the function (see (25)). This fact implies that and all variable coefficients in system (5) are infinitely differentiable, even when the function is discontinuous or nondifferentiable. In Figure 14, the coefficient is displayed for a (nondifferentiable) piecewise linear topography. We see that the coefficient is a regularized version of the physical topography.(ii)One-dimensional system (5) models bidirectional waves and it incorporates the simultaneous effects of surface tension and variable depth upon the shape of a water wave that propagates on the surface of an irregular shallow channel.(iii)Furthermore, in the derivation of system (5), we do not eliminate the wave elevation (which prevents us from neglecting additional terms of order ), as done, for example, in the model derived by Milewski [23]. This means that we did not reduce (30)-(31) to only one equation for the potential velocity where does not appear. Therefore, new formulation (5) is expected to be a more accurate approximation of the original potential theory equations.
In second place, we introduced a numerical Galerkin-spectral discretization of second order in time and spectral accuracy in space to approximate solutions of system (5) and a Fourier-collocation strategy combined with a Newton-type iterative procedure to compute its travelling wave solutions. In all computer simulations, these solutions were observed to propagate approximately without changing their shape within a large time interval.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This research was supported by Universidad del Valle, Calle 13 No. 100-00, Cali, Colombia, under the Research Project C.I. 988.