*2.1. Model Equation with Two-Bump Configurations*

The forced KdV equation is written as follows

$$\begin{aligned} \eta\_t(\mathbf{x},t) + 2\lambda \eta\_x(\mathbf{x},t) - \frac{3}{2} (\eta^2)\_x(\mathbf{x},t) - \frac{1}{3} \eta\_{xxx}(\mathbf{x},t) &= f\_x(\mathbf{x}), \quad -\infty < \mathbf{x} < \infty \\ \eta(\mathbf{x},0) = \eta\_0(\mathbf{x}) \\ \eta(\mathbf{x},t) = \eta\_x(\mathbf{x},t) = \eta\_{xx}(\mathbf{x},t) &= 0 \text{ at } \mathbf{x} \pm \infty, \end{aligned} \tag{1}$$

where *η*(*<sup>x</sup>*, *t*) is the free surface elevation measured from an undisturbed water level, and *λ* is the measurement of the perturbation of the upstream uniform flow velocity from its critical value. The external forcing *f*(*x*) is given by the topography of the rigid bottom. Our main interests in this work are numerical stability when stationary wave solutions evolve in time. Hence, the problem (1) defined in (−∞, ∞) can be restricted to a finite interval [−*b*, *b*] for a positive constant *b*, where this interval is large enough to use a finite difference scheme. First, at these finite boundaries, we employ the artificial boundary condition, which has been proposed in [21]. Then, the time-independent nonlinear boundary value problem is efficiently solved using the Newton method with the artificial boundary condition which is described as

$$\begin{cases} 2\lambda \eta(\mathbf{x}) - \frac{3}{2}\eta^2(\mathbf{x}) - \frac{1}{3}\eta\_{xx}(\mathbf{x}) = f(\mathbf{x}), \ -b \le \mathbf{x} \le b\\ \eta\_{\mathbf{x}}(\mathbf{x}) = \mp\sqrt{6\lambda}\eta(\mathbf{x}) \text{ at } \mathbf{x} = \pm b \end{cases} \tag{2}$$

As shown in [21], the method provides a good alternative to find various stationary solutions under arbitrary two-bump configurations. Once the stationary wave solutions have been found, we take these stationary solutions as initial conditions in (1). Next, the finite difference method described in the next subsection is performed to solve the time-dependent forced KdV equation in a finite interval *x* ∈ [−*b*, *b*] for *t* > 0.

Moreover, the bottom topography *f*(*x*) is modeled by the following functions for two bumps

$$h(\mathbf{x}; a, P\_1, P\_2) = \begin{cases} P\_1 \cos^4(\frac{\pi}{2}(\mathbf{x} + a)) & \text{if } |\mathbf{x} + a| \le 1, \\\ P\_2 \cos^4(\frac{\pi}{2}(\mathbf{x} - a)) & \text{if } |\mathbf{x} - a| \le 1, \\\ 0 & \text{elsewhere.} \end{cases}$$

Note that *h*(*x*; *a*, *P*1, *<sup>P</sup>*2) consists of two bumps (centered at *a* and −*<sup>a</sup>*) with the distance between the boundaries of two bumps, *d* = <sup>2</sup>(*a* − <sup>1</sup>).
