4.1. Summation of double series with symmetric elements

The double series such as (3.49)–(3.51), as well as many others, in § 3 can be written in the form

(4.1) \begin{equation} \mathcal{S} = \sum _{n=-2}^{\infty } \sum _{m=-2}^{\infty } A_{n m}, \end{equation}

where $A_{n m} = A_{m n}$ for any $n, m \geqslant -2$ . We have

(4.2) \begin{align} \mathcal{S} = \sum _{n=-2}^{\infty } A_{n n} + \sum _{n=-1}^{\infty } \sum _{m=-2}^{n-1} A_{n m} + \sum _{m=-1}^{\infty } \sum _{n=-2}^{m-1} A_{n m} = 2 \sum _{n=-1}^{\infty } \sum _{m=-2}^{n-1} A_{n m} + \sum _{n=-2}^{\infty } A_{n n}. \end{align}

For example, the formula (3.51) for $x\gt \xi$ would read as

(4.3) \begin{align} W^{(2+)}_p(x) &= - \frac {1}{4 \omega ^2} \sum _{n=-2}^{\infty } \frac {\kappa _n^8}{(2 \kappa _n)^4 Q_n^2} \ \tilde {E}_{n n} \ \textrm{e}^{2 {i} \kappa _{n} (x-\xi )} \nonumber \\ & \quad - \frac {1}{2 \omega ^2} \sum _{n=-1}^{\infty } \frac {\kappa _n^4}{Q_n} \textrm{e}^{\textrm{i} \kappa _{n} (x-\xi )} \sum _{m=-2}^{n-1} \frac {\kappa _{m}^{4}}{Q_{m}\left (\kappa _{n}+\kappa _{m}\right )^{4}} \ \tilde {E}_{n m} \ \textrm{e}^{\textrm{i} \kappa _{m}(x-\xi )}. \end{align}

Equation (4.2) helps to accelerate convergence of double series with symmetric elements.

4.2. Convergence of the numerical results

The numerical results corresponding to the hydroelastic vibrations of an ice plate fixed to a wall and subjected to an oscillating point pressure at $x=\xi$ are obtained and investigated for ice density $\rho _i = 917 \,\textrm{kg m}^{-3}$ , Young’s modulus $E = 4.2 \times 10^9 \,\textrm{N m}^{-2}$ , Poisson’s ratio $\nu = 0.33$ and ice thickness $h_i = 1.5 \,\textrm{m}$ (Squire et al. Reference Squire, Robinson, Langhorne and Haskell1988). Thus, the mass of the ice plate per unit area is $M = 1375.5 \,\textrm{kg m}^{-2}$ and the flexural rigidity is $D = 1.32\times 10^9\,\textrm{N m}$ . Further, the water density is $\rho = 1026 \,\textrm{kg m}^{-3}$ and the physical water depth is $H^{\prime } = 15 \,\textrm{m}$ . Hence, the characteristic length of this ice plate is $L = 19.05 \,\textrm{m}$ and, correspondingly, $\alpha = M/(\rho L) = 0.0704$ . Similar characteristics of ice have been used in several earlier studies, see for example Brocklehurst et al. (Reference Brocklehurst, Korobkin and Părău2010) and Wu et al. (Reference Wu, Khabakhpasheva, Ni and Korobkin2023a ). The water depth is comparable with the characteristic length $L$ , which implies that the finite depth matters in our analysis. The frequency and position of the external point pressure will vary in the following calculations.

We obtain the roots $\kappa _n$ of the dispersion relation (3.2) using the Wolfram Mathematica software. To check the correctness of the numerical values of these roots, and decide upon a truncation number $N$ for the series in (3.6) and (3.7), we calculate the following series ((47) and (88) of Korobkin et al. Reference Korobkin, Malenica and Khabakhpasheva2019):

(4.4) \begin{align} l_0 (\omega ) = \sum _{n=-2}^{\infty } \frac {1}{Q_n}, \quad l_1(\omega ) - \omega ^2 = \sum _{n=-2}^{\infty } \frac {\kappa _n^2}{Q_n} - \omega ^2, \quad l_2(\omega ) = \sum _{n=-2}^{\infty } \frac {\kappa _n^4}{Q_n}, \end{align}

which must all be equal to zero for any dimensionless frequency $\omega$ .

After performing the numerical calculations for higher values of $\omega$ (up to 100), it was observed that the series $l_j(\omega ) \,(j=0,1,2)$ keep growing quadratically and their values remain no more close to zero. It is only within the interval $(0,3.77)$ that the accuracy of the series $l_j$ is in the range $10^{-9} {-} 10^{-20}$ , see figures 3 and 4. The decreasing accuracy with increasing $\omega$ can be attributed to numerical sensitivity. At high frequencies, the large and oscillatory nature of the wavenumbers $\kappa _n$ not only complicates root-finding, but also leads to catastrophic cancellation in the calculation of $Q_n$ , where the subtraction of large terms significantly amplifies round-off errors. Similar behaviour of the series was reported by Korobkin et al. (Reference Korobkin, Malenica and Khabakhpasheva2019), where the roots $\kappa _n$ were determined using another algorithm. Thus, for all the results hereafter, the dimensionless frequency $\omega$ is considered only in the interval $(0,3.77)$ , which is for long waves. For a characteristic length $L=19.05 \,\textrm{m}$ , this corresponds to a dimensional angular frequency $\omega ^\prime$ in the range $(0,2.70) \, \textrm{rad s}^{-1}$ , covering long flexural-gravity waves with wavelengths greater than approximately 70 m $(\kappa _0^\prime = \kappa _0/L = 1.7431/L)$ . All series in our calculations are truncated at 200 imaginary roots of the dispersion relation (3.2).

Figure 3.

Convergence of series $l_2$ with the number of imaginary roots of the dispersion relation (3.2) for $\omega = 2$ .

Figure 4.

The series $l_0$ , $l_1 - \omega ^2$ and $l_2$ against dimensionless frequency $\omega$ .

4.3. Limiting cases

4.3.1. Static deflection and the small-frequency limit

In the static regime $(\omega \to 0)$ , the fluid is at rest and the velocity potential vanishes, $\phi = 0$ . In our original formulation, the static problem reduces to the plate (2.3), where $\phi = 0$ , subject to the boundary conditions (2.4) and the far-field condition $w \to 0$ as $x \to \infty$ . This solution would be valid for arbitrary values of $\varepsilon$ . It should be noted that this is a nonlinear problem, which implies that correction to the linear solution is of order $\mathcal{O} (\varepsilon ^2 )$ and is due to the nonlinearity of the plate equation only. This nonlinear boundary value problem can be solved numerically. We limit ourselves to the case of small $\varepsilon$ .

In the leading order approximation as $\varepsilon \to 0,$ the static linear solution can be obtained analytically. Since there is no fluid motion, the right-hand side of (2.17) becomes zero and the leading-order problem reduces to the static fourth-order ordinary differential equation with a point load, whose solution is given by

(4.5) \begin{align} W^{(1)}_{\textit{static}}(x) = \begin{cases} \mathcal{C}_1 \ \textrm{e}^{(-1+\textrm{i}) \mu x} + \mathcal{C}_2 \ \textrm{e}^{(-1-\textrm{i}) \mu x} \quad (x\gt \xi ), \\ \mathcal{D}_1 \ \textrm{e}^{(1+\textrm{i}) \mu x} + \mathcal{D}_2 \ \textrm{e}^{(1-\textrm{i}) \mu x} + \mathcal{D}_3 \ \textrm{e}^{(-1+\textrm{i}) \mu x} + \mathcal{D}_4 \ \textrm{e}^{(-1-\textrm{i}) \mu x} \quad (0\lt x\lt \xi ), \end{cases} \end{align}

where $\mu = 1/\sqrt {2}$ , and the unknown coefficients $\mathcal{C}_1,$ $\mathcal{C}_2,$ $\mathcal{D}_1, \mathcal{D}_2, \mathcal{D}_3$ and $\mathcal{D}_4$ are determined using the edge conditions (2.18) at the wall, and the matching conditions (3.8) and (3.9) at $x=\xi$ .

The second-order problem (2.23)–(2.30) in the static case gives $w^{(2)} = 0$ . This is because in the static regime, all time-dependent terms vanish and the fluid potential is zero, reducing (2.25) to a homogeneous equation with homogeneous boundary conditions, see (2.26) and (2.30). Thus, within the second-order approximation, the static deflection is given by

(4.6) \begin{equation} w_{\textit{static}}(x) = {\textrm{Re}}\,\{W^{(1)}_{\textit{static}}(x)\} + \mathcal{O}\big (\varepsilon ^2\big ). \end{equation}

To validate our dynamic formulation, we compute the full second-order unsteady solution, derived in § 3, for a small frequency $\omega =0.01$ and compare it with the static solution (4.5). In this quasi-static limit, the inertial effects in both the fluid and the ice plate are expected to be negligible. Therefore, no propagating waves are present and the dynamic solution should converge to the static displacement profile. This is confirmed in figure 5, where the second-order dynamic solution for $\omega =0.01$ , shown by the solid black curve, matches closely the static solution (4.6) shown by the blue circles.

Figure 5.

Comparison of the static deflection profile (4.6) with the unsteady second-order solution for $\omega =0.01$ .

4.3.2. Oscillating point pressure at large distance from the wall

Within the leading-order approximation, the wave is generated by an oscillating point pressure applied far from the wall $(\xi \gg 1)$ , and propagating towards it can be regarded as an incident wave with frequency $\omega$ and dimensionless wavenumber $\kappa _0$ . The amplitude of this wave depends on the strength of the external pressure, but is independent of the presence of the wall. For the point pressure given in (2.50), the wave amplitude is equal to $|C_0|= 1/(2 \omega ^2 \kappa _0 Q_0)$ , see (3.7) and (3.14). Therefore, our leading-order results for large $\xi$ can be directly compared with those of Brocklehurst et al. (Reference Brocklehurst, Korobkin and Părău2010), who studied ice deflection near the wall due to an incident wave. All physical parameters in our calculations, see § 4.2, match those of Brocklehurst et al. (Reference Brocklehurst, Korobkin and Părău2010), except for the water depth and ice thickness, which were $H^\prime =350 \,\textrm{m}$ and $h_i = 1.6 \,\textrm{m}$ , respectively, in their study. Consequently, this requires updating $L = 19.8942 \,\textrm{m}$ and $\alpha = 0.0719$ .

We choose $\omega = 0.5295$ for the comparison and set $\xi =500$ . This frequency $\omega$ was selected so as to fit the characteristics of the incident wave used by Brocklehurst et al. (Reference Brocklehurst, Korobkin and Părău2010). Note that Brocklehurst et al. (Reference Brocklehurst, Korobkin and Părău2010) used the water depth as the characteristic length scale, so the dimensionless distance from the wall in their notation was $x^{*} = x^\prime /H^\prime$ . To enable a direct comparison, we rescale our results using the same scales. The ice deflection near the wall, $0\lt x^*\lt 3$ , is presented in figure 6, where the blue solid line represents our solution and the red dots correspond to data extracted from figure 9(a) of Brocklehurst et al. (Reference Brocklehurst, Korobkin and Părău2010). The comparison demonstrates good agreement between the two solutions.

While the presented framework can be extended to wave scattering problems, this study intentionally focuses on the foundational case of a localised force to systematically isolate and understand the key nonlinear mechanisms such as resonant second-harmonic generation and strain localisation without the added complexity of a full wave scattering problem. While the framework in this subsection demonstrates its capability to model an incident wave field, we chose not to pursue the second-order response for that specific case in this work to maintain a sharp focus on the novel physics of the forced response near a boundary.

Figure 6.

Comparison of ice plate deflection near the vertical wall $(x=0)$ for waves generated by a distant oscillating pressure $(\xi \gg 1)$ and for the incoming wave case studied by Brocklehurst et al. (Reference Brocklehurst, Korobkin and Părău2010).

4.4. Comparison of leading order and second-order strains

The relative elongations (or contractions) of the filaments (or fibres) of an elastic plate are referred to as strains. The elastic strains are maximal on the upper and lower surfaces of the ice plate. In dimensionless variables, they are given by

(4.7) \begin{equation} \epsilon (x, t) = \mp \epsilon _{\textit{sc}} \ \frac {\partial ^2 w}{\partial x^2}, \end{equation}

where $\epsilon _{\textit{sc}} = (h_i \varepsilon )/(2 L)$ is the strain scale. The negative sign applies to the upper surface of the plate, while the positive sign corresponds to the lower surface. Tensile strain is indicated by a positive value, whereas compressive strain is represented by a negative value. The yield strain for ice $\epsilon _Y$ is estimated as $8 \times 10^{-5}$ , see Brocklehurst et al. (Reference Brocklehurst, Korobkin and Părău2010) and Korobkin & Khabakhpasheva (Reference Korobkin and Khabakhpasheva2025) for discussions on this value. We assume that the elastic strains in a continuous ice caused by an external load cannot exceed the yield strain value. We use this value as a reference threshold to identify regions where the induced elastic strain is critically high. It is important to note that this value applies to the failure of a free-floating ice plate. The fracture process at the clamped wall is more complex, likely involving localised crushing, spalling or a different failure mode due to the extreme stress concentration and constraint. Therefore, while strains exceeding this yield value in the interior strongly indicate a high risk of crack initiation, the actual failure mechanism at the wall boundary itself may differ and require a separate, more localised analysis (Zou, Xiao & Jordaan Reference Zou, Xiao and Jordaan1996).

Using the asymptotic expansion (2.13) for the ice deflection $w(x,t)$ and the expressions (2.33), (2.38) for the leading order deflection $w^{(1)}(x,t)$ and the second-order deflection $w^{(2)}(x,t)$ , respectively, we obtain the first-order strain $\epsilon ^{(1)}(x,t)$ as

(4.8) \begin{equation} \frac {\epsilon ^{(1)}}{\epsilon _{\textit{sc}}} =\mp \textrm{Re}\, \left \{ W_{xx}^{(1)}(x) \ \textrm{e}^{-\textrm{i} \omega t}\right \} \end{equation}

and the full second-order strain $\epsilon ^{(2)}(x,t)$ as

(4.9) \begin{align} \frac {\epsilon ^{(2)}}{\epsilon _{\textit{sc}}} = \mp {Re}\left \{ W_{xx}^{(1)}(x) \textrm{e}^{-\textrm{i} \omega t} + \varepsilon \ W_{xx}^{(2)}(x) \ \textrm{e}^{-2 {i} \omega t} + \varepsilon \ W_{0xx}^{(2)}(x)\right \}. \end{align}

Both strains $\epsilon ^{(1)}(x,t)$ , which is predicted by the linear theory, and $\epsilon ^{(2)}(x,t)$ , which is the second-order strain determined in the present study, are periodic in time with the period $2\pi /\omega$ . In the following, we focus on the maximum in time amplitudes:

(4.10) \begin{equation} \epsilon _m^{(1)}(x)= \max _{0 \lt t \lt \frac {2 \pi }{\omega }} |\epsilon ^{(1)}(x,t)| = \epsilon _{\textit{sc}} \ |W^{(1)}_{xx} (x)|, \quad \epsilon _m^{(2)}(x)=\max _{0\lt t \lt \frac {2 \pi }{\omega }} |\epsilon ^{(2)}(x,t)|. \end{equation}

The linear strain amplitude $\epsilon _m^{(1)}(x)$ depends on the scale $\epsilon _{\textit{sc}}$ , which is proportional to the strength of the external load through the small parameter $\varepsilon$ , and the function $W^{(1)}(x)$ , which depends on the structural parameter $\alpha$ , the dimensionless frequency $\omega$ , the dimensionless distance of the application of the load $\xi$ from the wall and the dimensionless water depth $H$ . Correspondingly, the inequality $\epsilon _m^{(1)}(x)\lt \epsilon _Y$ can be treated as a restriction in terms of the parameter $\varepsilon$ and finally in terms of the strength of the external point load.

To contextualise the magnitude of the forcing, we define a critical value $\varepsilon _*$ that estimates when the linear strain would reach the yield threshold. This provides a reference scale to ensure our weakly nonlinear analysis remains in a physically plausible regime where the ice would not have already fractured according to a linear model. Thus, acceptable values of $\varepsilon$ lie in the range $0\lt \varepsilon \lt \varepsilon _*$ , where

(4.11) \begin{equation} \varepsilon _* = \frac {2 L}{h_i} \: \frac {\epsilon _Y}{\max _{x\geqslant 0} | W_{xx}^{(1)}(x)|}. \end{equation}

We restrict our analysis to $0\lt \varepsilon \lt \varepsilon _*$ , ensuring that the linear component of the strain does not surpass the yield strain $\epsilon _Y$ . This does not imply the plate model is invalid up to this point; rather, it is a conservative criterion to maintain the relevance of our elastic model before the onset of pervasive plastic deformation or fracture.

Figure 7.

Magnitude of $W^{(1)}_{xx}(x)$ in the ice plate as a function of $x$ for (a) varying values of $\xi$ and fixed dimensionless frequency $\omega =1.4$ and (b) varying values of $\omega$ for fixed $\xi = 10$ .

The function $|W^{(1)}_{xx}(x)|$ is shown in figure 7 for different $\xi$ and the frequency $\omega$ . Other parameter values are listed in § 4.2. The figure reveals how the load position $\xi$ and frequency $\omega$ govern the failure risk. The key finding is a fundamental shift in the location of maximum strain, which dictates where the ice is most likely to fracture. For low frequencies ( $\omega \lt 0.5$ ) and loads positioned away from the wall ( $\xi \gt 2$ ), the maximum strain occurs not at the wall, but directly beneath the load itself (see figure 7 b). Physically, the ice responds quasi-statically, bending locally under the oscillating load. As the frequency increases or the load moves closer to the wall, the behaviour of the system changes. Wave reflection from the rigid boundary becomes dominant, creating a standing wave pattern that pins the maximum strain at the clamped edge ( $x=0$ ), as seen in figure 7(a). The wall, not the load, becomes the critical point of failure.

To illustrate the physical implications, consider a point load $190.5\,\textrm{m}$ from the wall ( $\xi =10$ ) oscillating at $\omega ^\prime = 0.6\,\textrm{rad s}^{-1}$ ( $\omega =0.84$ ). Our analysis gives $\varepsilon _* \approx 0.0036$ , corresponding to a total load magnitude of approximately $13.15 \times 10^3 \, \textrm{N m}^{-1}$ , which provides a conservative threshold for critically high elastic strains in the ice plate. This threshold serves as a practical reference for assessing strain localisation under oscillatory loading, rather than a prediction of fracture at the wall. Table 1 summarises this for the parameters in figure 7. A practical problem addressed in this paper concerns engineering activities on floating ice. These include the use of vibrators placed on the ice to investigate ice properties and bottom topography, to detect subsurface flows and moving objects beneath the ice, and to image subsurface geological structures, particularly in the context of oil and gas exploration (Holdsworth & Glynn Reference Holdsworth and Glynn1978; Eisen et al. Reference Eisen, Hofstede, Diez, Kristoffersen, Lambrecht, Mayer, Blenkner and Hilmarsson2015).

Table 1.

Maximum values of the linear strain $|W^{(1)}_{xx}(x)|$ and corresponding critical strain amplitudes $\varepsilon _*$ for various parameter combinations. Results are extracted from figure 7.

Figure 8.

The maximum value $\varepsilon _*$ of the small parameter $\varepsilon$ as a function of frequency $\omega$ and load location $\xi$ : (a) full parameter range and (b) magnified view highlighting cases where the load is applied very close to the wall.

To further understand the dependence of ice plate strain response on the frequency and position of the point load, we compute $\varepsilon _*$ using the formula (4.11) over a range of non-dimensional frequencies $\omega \in (0,3.5]$ and loading locations $\xi \in (0,10]$ , see figure 8(a). We observe a minimum values of strains for $\xi \gt 3$ and $0.5 \lt \omega \lt 3.5$ , and maximum values of strains for small $\xi$ . The magnified view in figure 8(b) focuses on cases where the load is applied very close to the wall. To contextualise these values in a physical scenario, we relate the parameter $\varepsilon _*$ to a physical load via $F_{phys} = \varepsilon _* \rho g L^2$ N m⁻¹. For our parameters ( $\rho =1026 \,\textrm{kg m}^{-3}$ , $L=19.05$ m), the range $\varepsilon _* = 0.004-0.022$ investigated corresponds to physical loads of $14.6$ kN m⁻¹ to $80.4$ kN m⁻¹. Distributed over a characteristic width, this represents the load of a snowmobile (over 0.5 m) up to light aircraft or industrial equipment (over 3 m), confirming the analysis addresses realistic forcing. The figures show that $\varepsilon _*$ is more sensitive to the distance of the load from the wall than to the frequency of the load. Note that $\xi =0.2$ corresponds to $3.81 \,\textrm{m}$ in the dimensional variables, which is small compared with the water depth of $15\,\textrm{m}$ in our calculations. The locations, $x_{{max}}(\omega ,\xi )$ , where the maximum strains are achieved, are shown in figure 9. Only for frequencies smaller than 0.5 and dimensionless distances $\xi \gt 2$ , the maximum strain is achieved not at the wall. For the range $\omega \lt 0.5$ and $\xi \gt 2$ , the position $x_{{max}}$ of the maximum strain tends to coincide with the load position $\xi$ . Overall, these findings highlight the strong influence of both the loading frequency and position on the strain distribution. They provide valuable guidance for the design and safe operation of floating structures and point-load systems acting on ice plates.

Figure 9.

The $x$ -coordinate of the point on the ice plate, $x_{{max}}$ , where the maximum of $|W^{(1)}_{xx}(x)|$ is achieved.

Figure 10.

Strain components for $\xi = 0.2$ showing (a,b) magnitude of the leading order scaled strain $\epsilon ^{(1)}_m(x)/\epsilon _Y$ and (c,d) the difference between leading and second-order strains scaled by $\varepsilon _*$ , $(\epsilon _m^{(2)}(x) - \epsilon _m^{(1)}(x) )/\epsilon _Y/\varepsilon _*$ , both as functions of the distance $x$ from the wall. Results are shown for (a, c) $\omega = 0.0056$ , $\varepsilon _* = 0.0117$ and (b,d) $\omega = 2.8$ , $\varepsilon _* = 0.0132$ .

Figure 11.

Strain components for $\xi = 5$ showing (a,b) magnitude of the leading order scaled strain $\epsilon ^{(1)}_m(x)/\epsilon _Y$ and (c,d) the difference between leading and second-order strains scaled by $\varepsilon _*$ , $(\epsilon _m^{(2)}(x) - \epsilon _m^{(1)}(x) )/\epsilon _Y/\varepsilon _*$ , both as functions of the distance $x$ from the wall. Results are shown for (a, c) $\omega = 0.0056$ , $\varepsilon _* = 0.00575$ and (b,d) $\omega = 2.8$ , $\varepsilon _* = 0.00522$ .

Figure 12.

Strain components for $\xi = 10$ showing (a,b) magnitude of the leading order scaled strain $\epsilon ^{(1)}_m(x)/\epsilon _Y$ and (c,d) the difference between leading and second-order strains scaled by $\varepsilon _*$ , $(\epsilon _m^{(2)}(x) - \epsilon _m^{(1)}(x) )/\epsilon _Y/\varepsilon _*$ , both as functions of the distance $x$ from the wall. Results are shown for (a, c) $\omega = 0.0056$ , $\varepsilon _* = 0.00574$ and (b,d) $\omega = 2.8$ , $\varepsilon _* = 0.00522$ .

The scaled leading order strains $\epsilon _m^{(1)}(x)/\epsilon _Y$ are shown in panels (a,b) of figures 10, 11 and 12 for two dimensionless frequencies, $\omega = 0.0056$ and $\omega =2.8$ , and three load distances from the wall, $\xi = 0.2, 5, 10$ . Panels (c,d) in the same figures compare the magnitude of the leading order strains with the magnitude of the second-order strains $\epsilon _m^{(2)}(x)/\epsilon _Y$ (see (4.10)) for the same parameter combinations. In each case, $\varepsilon$ is taken as its upper limit $\varepsilon _*$ defined by (4.11). The second-order strain (nonlinear strain) in the ice plate is obtained from (4.9) using the expressions for $W^{(1)}(x)$ , $W^{(2)}(x)$ and $W_0^{(2)}(x)$ derived in §§ 3.1, 3.2 and 3.3, respectively. Note that here, we focus on the second-order corrections to the leading order strain.

The results illustrate the pronounced influence of second-order corrections on the strain field, particularly for frequencies close to the resonant frequencies. We did not calculate all resonant frequencies in the present study. One of such frequencies detected is $\omega = 0.0014$ . At this frequency, $2 \kappa _0$ is the real root of the dispersion relation (3.2) for the double frequency $2 \omega$ , where $\kappa _0$ is the real root of (3.2) for the original frequency $\omega$ . In the low-frequency case of $\omega = 0.0056$ , a marked deviation is observed between the nonlinear strain and the linear prediction. This deviation becomes especially significant due to the near-singular behaviour of the second-order coefficients $E_{nm}$ (see (3.45) and the discussion following the equation).

For a load placed very close to the wall at $\xi =0.2$ , the leading order strain attains its maximum at the wall and decays rapidly with distance, see figure 10(a), and approaches a constant value in the far field for moderate frequencies, see figure 10(b). The sharp changes of strains in figure 10(a,b) indicate that the strain derivative is not continuous at $x=\xi$ , see (3.8), for a point load. However, as the load is moved away from the wall to $\xi =5$ (figure 11 a), the peak of the linear strain shifts towards the load location and exhibits increased amplitude and spatial extent. Figure 10(c) shows that the second-order correction of the strains is small for the loads applied near the wall. The maximum deviation of the second-order strains from the leading order ones is achieved at the wall, where it is estimated as $-0.4 \varepsilon _* \epsilon _Y$ . The nonlinear effects reduce the strains at the wall, but increase them at the place where the load is applied, see figures 10(c) and 10(d). However, for $\xi =5$ , see figure 11(c), the normalised second-order correction reaches values as low as $-12 \varepsilon _* \epsilon _Y$ , highlighting the dominance of second-order effects at resonance and remote loading. The nonlinear effects reduce the maximum strains for small frequencies, see figures 10(c), 11(c), 12(c). Their contributions for moderate frequencies are less clear, see figures 10(d), 11(d), 12(d). A similar pattern is observed at $\xi =10$ , see figure 12, where the leading-order strain peaks near the load and the second-order correction remains significant, with a maximum deviation of approximately $-10 \varepsilon _* \epsilon _Y$ . The nonlinear correction becomes increasingly important for remote loading cases $(\xi =5,10)$ and low frequencies, with deviations from the linear prediction reaching up to 10 times $\varepsilon _* \epsilon _Y$ . This indicates that for accurate modelling, especially for frequencies close to the resonant frequencies, second-order effects must be retained, particularly when the load is placed far from the boundary.

In contrast, at a moderate frequency of $\omega = 2.8$ , the linear strain response behaves differently. For all values of the load position $\xi$ , the maximum linear strain consistently occurs at the wall, rather than at the location of the load. For instance, when $\xi = 5$ (see figure 11 b), the normalised linear strain $\epsilon _m^{(1)}(x)/\epsilon _Y$ at $x = 5$ drops below 1, indicating that the ice is not close to yielding at the load point. Additionally, the second-order correction across the domain remains small, with normalised deviations within $\pm 1$ , confirming the weak influence of nonlinear effects. These results highlight that in contrast to the low-frequency regime, the moderate-frequency case exhibits a weak strain response and minimal nonlinear effects. Nonlinear effects become significant primarily in the low-frequency regime, particularly when the applied load is positioned away from the wall. In these cases, the second-order contributions noticeably alter the strain field. In contrast, at relatively higher frequencies, the linear solution remains a reliable and accurate approximation across all load positions $\xi$ , indicating that nonlinear corrections are minimal in this regime. Importantly, the second-order corrections are seen to be both positive and negative, indicating that nonlinearity may either amplify or reduce the leading-order strain depending on the spatial location and loading configuration, an effect that could influence the onset and localisation of material failure in practice.

Figure 13.

Contour plot showing the scaled difference $\Delta \varepsilon = (\epsilon _m^{(2)}(x_{{max}}) - \epsilon _m^{(1)}(x_{{max}}) )/\epsilon _Y/\varepsilon _*$ across the parameter space of $\omega$ and $\xi$ . The quantity represents the difference between the maxima of second-order strain and maxima of leading order strain normalised by $\varepsilon _Y$ and the scale $\varepsilon _*$ .

Figure 13 presents the scaled difference between the maximum second-order strain and the maximum leading order strain, $\Delta \varepsilon = (\epsilon _m^{(2)}(x_{{max}}) - \epsilon _m^{(1)}(x_{{max}}) )/\epsilon _Y/\varepsilon _*$ , over the $(\omega , \xi )$ parameter space. Owing to the definition of the reference scale $\varepsilon _*$ , see (4.11), which is based on the maximum of the leading order strain, the scaled leading order maximum is identical for all $(\omega , \xi )$ , and the observed variation in $\Delta \varepsilon$ therefore reflects solely the behaviour of the second-order strain. A distinct transition is evident near $\xi =2$ . As shown in figure 9, this value of $\xi$ marks a change in the location of the maximum leading order strain, from a region adjacent to the wall for $\xi \lt 2$ to a region close to the oscillating point load for $\xi \gt 2$ . This spatial reorganisation of the dominant strain location leads to a marked change in the second-order strain. Consequently, the difference between the second- and first-order strain maxima exhibits a sharp variation across $\xi =2$ , separating a wall-dominated regime from a load-dominated regime in the weakly nonlinear response.

The finding that nonlinearities are most significant at lower frequencies, not relatively higher, stems from two key factors. First, resonant energy transfer is crucial. At low frequencies, the system is near resonance for second-harmonic generation, where the self-interaction of the primary wave efficiently pumps energy into a freely propagating second harmonic, leading to large, growing nonlinear strains. At moderate frequencies, the system is off-resonance, suppressing this mechanism. Second, the role of the wall evolves with frequency. At moderate frequencies, wave reflection from the wall creates a standing-wave pattern that localises the maximum strain at the boundary itself. This strong linear boundary-layer effect dominates the response, leaving little room for nonlinearities to manifest. At low frequencies, however, the influence of the wall is more geometric, setting up a domain where resonant, propagating nonlinear waves can develop and dominate the strain field away from the boundary.

The significant second-order effects reported here are fully consistent with, and confined to, the regime of validity of the present weakly nonlinear formulation. The perturbation expansion is constructed under the explicit assumptions that all dependent variables and their spatial and temporal derivatives remain order one, that the ice behaves as a thin elastic plate ( $\alpha \lt 1$ ), and that the excitation frequency $\omega$ lies within the long-wave regime. Together with the requirement $\varepsilon _* \ll 1$ , these conditions ensure that the series remains asymptotically ordered. For the results presented, the values of the critical parameter ( $\varepsilon _* \leqslant 0.022$ , see figure 8) satisfy this condition, confirming that the large second-order corrections constitute a meaningful correction without violating the underlying assumptions. A violation of this condition very close to resonance would signify a breakdown of the weakly nonlinear theory, necessitating a more advanced analysis to capture the fully saturated state.