1. Introduction
When a clean rigid rectangular tank with a finite liquid depth is longitudinally excited with the forcing frequency close to the lowest natural sloshing frequency, the free-surface nonlinearity causes an energy flow (Faltinsen & Timokha Reference Faltinsen and Timokha2009; Pilipchuk Reference Pilipchuk2013) from the lowest (primary-excited) to higher natural sloshing modes so that the nonlinearity becomes a prevailing mechanism preventing an infinite resonant wave-amplitude response. Viscous damping plays then a secondary role (Faltinsen & Timokha Reference Faltinsen and Timokha2001). This explains why the single-dominant (applicable for smaller excitation amplitudes, Faltinsen et al. Reference Faltinsen, Rognebakke, Lukovsky and Timokha2000) and adaptive (it should be adopted with increasing the forcing amplitude and at the critical liquid depth, Faltinsen & Timokha Reference Faltinsen and Timokha2001) ‘inviscid’ weakly nonlinear modal theories/systems provide satisfactory agreement with the measured steady-state wave-amplitude characteristics. However, these theories fail to predict the phase-lag response curve, which is, theoretically, piecewise (the phase lag possesses the values 0 and 
$\pm {\rm \pi}$) in the inviscid approximation while accounting for damping makes it continuous. The phase-lag response is most sensitive to the damping in those mechanical systems. Hence, having known a set of measured phase lags and introducing physically and mathematically consistent unknown a priori damping terms in the modal systems, one could implement (according to the so-called hidden-physics concept by Raissi & Karniadakis Reference Raissi and Karniadakis2018; Raissi, Perdikaris & Karniadakis Reference Raissi, Perdikaris and Karniadakis2019; Ahmed et al. Reference Ahmed, Pawar, San, Rasheed, Iliescu and Noack2021), as discussed by Ahmed et al. (Reference Ahmed, Pawar, San, Rasheed, Iliescu and Noack2021) for the reduced order modelling as well as by Saltari et al. (Reference Saltari, Pizzoli, Coppotelli, Gambioli, Cooper and Mastroddi2022) for violent liquid sloshing due to vertical excitations, a machine learning procedure to implicitly deduce the introduced viscous damping terms.
Appropriate measurements of the phase lag were recently reported by Bäuerlein & Avila (Reference Bäuerlein and Avila2021) who utilised them to evaluate a constant damping ratio in the Duffing equation, which was adopted by the authors as a phenomenological mathematical model of horizontal motions of the liquid-mass centre (spring and pendulum-type phenomenological models are often employed for sloshing: see Miles Reference Miles1962; NASA 1968; Aliabadi, Johnson & Abedi Reference Aliabadi, Johnson and Abedi2003; Godderidge, Turnock & Tan Reference Godderidge, Turnock and Tan2012, and references therein). This damping ratio, 
$\xi$, had to equal 
$15\times 10^{-3}$ to more or less satisfactory approximate the experimental data. Using the same experimental set-up, Bäuerlein & Avila (Reference Bäuerlein and Avila2021) analysed logarithmic decrements of the first (lowest) natural sloshing mode and, thereby, computed its damping ratio 
$\xi _1^{exp}=8.4\times 10^{-3}$. This experimental value is larger than its theoretical prediction (
$\xi _1^{(0)}=5.7\times 10^{-3}$) by Keulegan (Reference Keulegan1959). Because the first differential equation of the single-dominant nonlinear modal system by Faltinsen et al. (Reference Faltinsen, Rognebakke, Lukovsky and Timokha2000) governs perturbations of the first natural sloshing mode, Bäuerlein & Avila (Reference Bäuerlein and Avila2021) inserted the linear damping term with 
$\xi _1^{exp}$ into this equation, derived a steady-state wave (periodic) solution, and showed that even in that case, the single-dominant modal theory is not able to satisfactory fit the measured steady-state resonant response.
A way of getting an excellent agreement is a purely data-driven (phenomenological) modelling dealing with dynamical systems with respect to amplitude and phase-lag parameters as it is demonstrated by Cenedese et al. (Reference Cenedese, Axås, Bäuerlein, Avila and Haller2022) exploiting the Stuart–Landau-type equation. Another way consists of following the hidden-physics concept suggesting a machine learning of physics-based dynamical models (ROMs) on sloshing derived from the original free-surface boundary problem. In the present paper, we step the second way. The single-dominant inviscid modal system by Faltinsen et al. (Reference Faltinsen, Rognebakke, Lukovsky and Timokha2000) is modified by incorporating the theoretically consistent damping terms. Implementing Moiseev's asymptotic scheme, a periodic (steady-state) solution of the modified system is derived, which, as the derivations show, depends on three unknown a priori parameters associated with the introduced damping. One parameter is responsible for the constant damping ratio of the second (-order) generalised coordinate (mode) but the other two are coefficients in the linear regression (by the dominant non-dimensional steady-state wave amplitude) chosen as the amplitude-dependent damping ratio for the first (lowest-order) generalised coordinate (mode). After defining the loss function as an integral distance between experimental (measured) and theoretical phase-lag curves, the gradient descent is applied to compute the three damping-related parameters. How it works is demonstrated with experimental phase lags by Bäuerlein & Avila (Reference Bäuerlein and Avila2021). An excellent agreement is shown except for the cases when, according to observations by Bäuerlein & Avila (Reference Bäuerlein and Avila2021), applicability of the single-dominant modal theory is questionable. Discrepancy between the linear damping quantities predicted by Keulegan (Reference Keulegan1959), coming from direct measurements by Bäuerlein & Avila (Reference Bäuerlein and Avila2021), and following from the present machine learning procedure is extensively discussed.
2. Single-dominant modal theory
We consider two-dimensional sloshing of an ideal liquid with irrotational flows in a rectangular tank of the width 
$L$ and breadth 
$B$. The non-dimensional hydrodynamic statement suggests that 
$L$ is the characteristic size and 
$T=2{\rm \pi} /\sigma _1$ is the characteristic time, where 
$\sigma _1$ is the lowest natural sloshing frequency. The tank is filled to a finite liquid depth which means that the non-dimensional (
$L$-scaled) depth 
$h \gtrsim 0.4$ (Faltinsen & Timokha Reference Faltinsen and Timokha2009, chap. 8). The tank oscillates harmonically with the forcing frequency 
$\sigma$ close to 
$\sigma _1$. Based on rigorous mathematical results by Faltinsen (Reference Faltinsen1974) and Ockendon & Ockendon (Reference Ockendon and Ockendon1973) concerning the steady-state resonant sloshing in a rectangular tank with a finite liquid depth (§ 8.2.2 in the textbook by Faltinsen & Timokha Reference Faltinsen and Timokha2009), Faltinsen et al. (Reference Faltinsen, Rognebakke, Lukovsky and Timokha2000) proved that, if the non-dimensional forcing amplitude is relatively small, the nonlinear resonant sloshing, both steady-state and transient, can be described by employing the so-called single-dominant nonlinear modal system whose hydrodynamic generalised coordinates obligatory satisfy specific asymptotic relations in terms of the forcing amplitude. After inserting the framed viscous damping terms, the system takes the following non-dimensional form:
\begin{align} &\ddot{\beta}_1 +\beta_1 +\fbox{$2\xi_1[\dot{\beta}_1+ \varXi_1 (\beta_m,\dot{\beta}_m |_{m=1,2})]$}+d_1(\ddot{\beta}_1\beta_2 +\dot{\beta}_1\beta_2) + d_2(\ddot{\beta}_1 \beta_1^2 +\dot{\beta}_1^2\beta_1)\nonumber\\ &\quad +d_3\ddot{\beta}_1\beta_1=P_1\eta_{2a}\bar{\sigma}^2\cos(\bar{\sigma} t -\theta), \end{align}
\begin{align} &\ddot\beta_2 +\bar{\sigma}_2^2\beta_2 +\fbox{$2\bar{\sigma}_2 \xi_2 [\dot\beta_2+\varXi_2(\beta_m,\dot\beta_m |_{m=1,2})]$} +d_4\ddot\beta_1\beta_1 +d_5\dot\beta_1^2 =0, \end{align}
\begin{align} &\ddot{\beta}_3 +\bar{\sigma}_2^2\beta_3 +\fbox{$2\bar{\sigma}_3\xi_3 [\dot{\beta}_3+ \varXi_3(\beta_m,\dot{\beta}_m |_{m=1,2})]$}+ q_1\ddot{\beta}_1\beta_2+ q_2 \ddot{\beta}_1\beta_1^2 + q_3\ddot{\beta}_2\beta_1 + q_4\dot{\beta}_1^2\beta_1\nonumber\\ &\quad+q_5\dot{\beta}_1\dot{\beta}_2 =P_3\eta_{2a}\bar{\sigma}^2\cos(\bar{\sigma}t -\theta), \end{align}
\begin{align} &\ddot\beta_m +\bar{\sigma}_m^2\beta_m +\fbox{$2\bar{\sigma}_m \xi_m \dot\beta_m$} =P_m\eta_{2a}\bar{\sigma}^2\cos(\bar{\sigma} t -\theta),\quad m\geqslant 4, \end{align}
where the hydrodynamic generalised coordinates 
$\beta _m(t)$ come from the functional (modal) representation of the free surface and necessarily fulfil the aforementioned Moiseev's asymptotic relations proven by Faltinsen (Reference Faltinsen1974) and Ockendon & Ockendon (Reference Ockendon and Ockendon1973) (see also § 8.2.2 by Faltinsen & Timokha Reference Faltinsen and Timokha2009),
$$\begin{gather} z=\zeta(y,t)=\sum_{m=1}^\infty \beta_m(t) \cos\left({\rm \pi} m \left(y+\frac{1}{2}\right)\right), \end{gather}$$
$$\begin{gather}\beta_1=O(\epsilon^{1/3}); \quad \beta_2=O(\epsilon^{2/3});\quad \beta_3=O(\epsilon); \quad \beta_n\lesssim O(\epsilon),\quad n\geqslant 4; \quad P_1\eta_{2a}=\epsilon \ll 1, \end{gather}$$
respectively. Here, 
$\sigma _m$ are the natural sloshing frequencies, 
$\bar {\sigma } =\sigma /\sigma _1$, 
$P_m=2/({\rm \pi} m) \tanh ({\rm \pi} m h) ((-1)^m-1), \bar {\sigma }_m^2=\sigma _m^2/\sigma _1^2=m \tanh ({\rm \pi} m h)/\tanh ({\rm \pi} h)$, the hydro- dynamic coefficients 
$d_n, q_n$ depend on the non-dimensional liquid depth 
$h$, 
$\eta _{2a}=O(\epsilon )\ll 1$ is the non-dimensional forcing amplitude, and 
$\theta$ is the phase lag in the external horizontal harmonic forcing. The computed values of the hydrodynamic coefficients are tabled by Faltinsen & Timokha (Reference Faltinsen and Timokha2009, chap. 9). A novelty with respect to the original work by Faltinsen et al. (Reference Faltinsen, Rognebakke, Lukovsky and Timokha2000) is the framed damping terms which we do not know a priori. As notified by Faltinsen & Timokha (Reference Faltinsen and Timokha2009, § 7.4), applicability of modal theories implicitly assumes a low-viscous liquid, i.e. the incorporated damping terms should be relatively small. What does it mean in the present paper is discussed below.
Because (2.1) is based on the Moiseev-type asymptotic relations (2.2b) with neglecting the 
$o(\epsilon )$-order quantities, only 
$\beta _1$ and 
$\beta _2$ are nonlinearly coupled (there is an energy flow between these two lower natural sloshing modes); 
$\beta _3$ is the ‘driven’ generalised hydrodynamic coordinate of the highest asymptotic order 
$O(\epsilon )$ but the natural sloshing modes with 
$m\geqslant 4$ are considered within the framework of the linear sloshing theory.
The same asymptotic relations should be true for the incorporated damping terms, that is, 
$\varXi _i$ are only nonlinear functions of 
$\beta _1$ and 
$\beta _2$ and their first derivatives. Furthermore, in the linear limit with 
$\eta _{2a}=0$ (free oscillations), the linear damping ratios 
$\xi _i$ can be associated with logarithmic decrements of the natural sloshing modes 
$\cos ({\rm \pi} i (y+\tfrac {1}{2}))$ in (2.2a). These ratios should be small for low-viscous liquids and, because our nonlinear analysis will centre around the two nonlinearly coupled generalised coordinates, 
$\beta _1$ and 
$\beta _2$, one can assume, as the ‘most worse’ case, that 
$\xi _1\sim \xi _2$ possess the lowest asymptotic order 
$O(\epsilon ^{1/3})$. Finally, because the asymptotic modal theory neglects the 
$o(\epsilon )$-order contribution in (2.1), one concludes 
$\varXi _3=\varXi _2=0$ but 
$\varXi _1=\varXi _1(\beta _1,\dot {\beta }_1)$ is a quadratic function, which we do not know a priori as well as 
$\xi _i$.
The steady-state resonant sloshing is associated with periodic solutions of (2.1). Following Moiseev’s asymptotic scheme (Faltinsen & Timokha Reference Faltinsen and Timokha2009, chap. 8), one can find an asymptotic approximation of these solutions even if (2.1) is equipped with the framed damping terms. In Moiseev's approximation, the lowest-order generalised hydrodynamic coordinates 
$\beta _1$ and 
$\beta _2$ take the form
\begin{equation} \beta_1(t)=a\cos(\bar{\sigma} t) + O(a^3), \quad \beta_2(t)=a^2[l_0+l_1\cos(2\bar{\sigma} t)+l_2\sin(2\bar{\sigma} t)]+O(a^4), \end{equation}
where 
$a=O(\epsilon ^{1/3})>0$ is the dominant (lowest-order) wave amplitude,
\begin{equation} l_0=\frac{d_4-d_5}{2\bar{\sigma}_2^2},\quad l_1(\xi_2)=\frac{(d_4+d_5)(\bar{\sigma}_2-4)}{2((\bar{\sigma}_2^2-4)^ 2+16\bar{\sigma}_2^2\xi_2^2)},\quad l_2(\xi_2)=\frac{2(d_4+d_5)\bar{\sigma}_2\xi_2}{(\bar{\sigma}_2^2-4)^2+16\bar{\sigma}_2^2\xi_2^2}; \end{equation}
the amplitude 
$a>0$ and the phase lag 
$\theta$ come from the so-called secular system
\begin{equation} a(m_1(\xi_2) a^2+\hat\sigma^2-1)=\epsilon\cos\theta, \quad a(m_2(\xi_2)a^2-2\xi_1[1 + \xi a])=\epsilon\sin\theta, \end{equation}
which is derived after substituting (2.3a,b) into (2.1a) and collecting all quantities at 
$\cos (\bar {\sigma } t)$ and 
$\sin (\bar {\sigma } t)$ after applying the cosine angle-difference identity in the right-hand side and dividing by 
$\bar {\sigma }^2$; here,
\begin{equation} m_1(\xi_2)={-}\tfrac 12 d_2 -2 d_3l_1+d_1({-}l_0+\tfrac 12 l_1), \quad m_2(\xi_2)=\tfrac 12 l_2(d_1-4d_3),\quad \hat\sigma=\bar{\sigma}^{{-}1} \end{equation}and
\begin{equation} \xi={-}\frac{\bar{\sigma}}{\rm \pi} \int_{0}^{2{\rm \pi}/\bar{\sigma}} \varXi_1(\cos(\bar{\sigma} t), -\sin (\bar{\sigma} t))\sin(\bar{\sigma} t)\,{\rm d} t \geqslant 0, \end{equation}remembering that, from physical reasons,
\begin{equation} \int_{0}^{2{\rm \pi}/\bar{\sigma}} \varXi_1(\cos(\bar{\sigma} t), -\sin (\bar{\sigma} t))\cos(\bar{\sigma} t)\,{\rm d} t =0. \end{equation}
The steady-state analysis does not require to explicitly know 
$\varXi _1(\beta _1,\dot \beta _1)$. Any quadratic function 
$\varXi _1$ satisfying (2.7) and (2.8) is allowed. An appropriate expression of 
$\varXi _1$ for transient wave computations by (2.1) is discussed in § 4 Conclusions.
Taking the sum of squares in (2.5a,b) makes it possible to rewrite the secular system to the form
$$\begin{gather} a^2[ (m_1(\xi_2)a^2 +(\hat\sigma^2-1))^2 + (m_2(\xi_2)a^2 -2\xi_1(1+\xi a))^2]=\epsilon^2, \end{gather}$$
$$\begin{gather}\theta=\mbox{atan2} ( [m_2(\xi_2)a^2-2\xi_1(1+\xi a)], \quad [m_1(\xi_2)a^2+(\hat\sigma^2-1)] ) \end{gather}$$
(
${\rm atan}2(y,x)$ is the 2-argument arctangent) that can be solved to get an analytical prediction of the resonant response curves. Indeed, (2.9a) gives
\begin{align} \bar{\sigma}^{{-}2}(a;\xi_1,\xi_2,\xi,\epsilon)&= \hat\sigma^2(a;\xi_1,\xi_2,\xi,\epsilon)\nonumber\\ &=1-m_1(\xi_2)a^2\pm \sqrt{\frac{\epsilon^2}{a^2}-[m_2(\xi_2)a^2-2\xi_1(1+\xi a)]^2}, \end{align}which has the physical meaning when both the right-hand side and expression under the square root are non-negative that yields the left and right bounds in
\begin{equation} a_{min}(\xi_1,\xi_2,\xi,\epsilon)\leqslant a\leqslant a_{max}(\xi_1,\xi_2,\xi,\epsilon). \end{equation}Furthermore, inserting (2.10) into (2.9b) derives
\begin{equation}
\theta(a;\xi_1,\xi_2,\xi,\epsilon) = \mbox{atan2} (
m_2(\xi_2)a^2-2\xi_1(1+\xi a),
m_1(\xi_2)a^2+(\hat\sigma^2(a;\xi_1,\xi_2,\xi,\epsilon)-1
),
\end{equation}
where 
$\hat \sigma ^2$ should be taken from (2.10).
Hence, when varying 
$a$ in the interval (2.11) for a fixed set of 
$\xi _1,\xi _2,\xi$ and 
$\epsilon$, (2.10) and (2.12) parametrically define the phase-lag response curve in the 
$(\bar {\sigma },\theta )$-plane but (2.10) as function of 
$a$ determines the wave-amplitude response curve in the 
$(\bar {\sigma },a)$-plane.
3. Learning 
$\xi _1$, $\xi _2$ and $\xi$ from measurements of the phase lag
The unknowns 
$\xi _1$, 
$\xi _2$ and 
$\xi$ can be ‘learnt’ by using a set of the measurements of 
$\theta _i^{(n)}$ with the given 
$\bar {\sigma }_1^{(i,n)}$ and 
$\epsilon _n=P_1\eta _{2a}^{(n)}, n=1,\ldots, N_{\eta _{2a}}$. For this purpose, we introduce the distance function 
$D$ between the phase-lag response curve by (2.10)–(2.12) and 
$(\bar {\sigma }_1^{(i,n)},\theta _i^{(n)})$ as
\begin{align} &D(n,i;\xi_1,\xi_2,\xi)\nonumber\\ &\quad =\sqrt{\min_{a_{min}\leqslant a\leqslant a_{max}} [ (\bar{\sigma}(a;\xi_1,\xi_2,\xi,\epsilon_n)-\bar{\sigma}_1^{(i,n)})^2+ (\theta(a;\xi_1,\xi_2,\xi,\epsilon_n) -\theta_i^{(n)})^2 ]}; \end{align}
it depends on the indexes 
$n$ and 
$i$ enumerating 
$\bar {\sigma }_1^{(i,n)}$, 
$\theta _i^{(n)}$ and 
$\epsilon _n$. Therefore, the integral distance
\begin{equation} C(\xi_1,\xi_2,\xi)=\sum_{n=1}^{N_{\eta_{2a}}}\sum_i D(n,i;\xi_1,\xi_2,\xi)>0 \end{equation}
between all the adopted measurement points and asymptotic solution (2.10)–(2.12) may be taken as the loss function in the machine learning procedure. The wanted damping-related parameters 
$\xi _1,\xi _2$ and 
$\xi$ realise the absolute minimum of the mean-square error-function 
$C$. Minimisation of the loss function (3.2) can be done by the gradient descent. The absolute minimum of 
$C$ exists and, from geometrical and physical point of view, is unique. Hypothetically, if the absolute minimum is realised on a manifold in the 
$(\xi _1,\xi _2,\xi )$-space, standard gradient-descent algorithms should effectively detect it.
3.1. Employing the measurements by Bäuerlein & Avila (2021)
Experimental data on the liquid-mass centre (horizontal steady-state wave amplitude and phase lag) were reported by Bäuerlein & Avila (Reference Bäuerlein and Avila2021). Horizontal position of the liquid-mass centre is, according to Faltinsen & Timokha (Reference Faltinsen and Timokha2009, (8.73)), described by
\begin{equation} y_C(t)={-}\frac{2}{{\rm \pi}^2 h}\beta_1(t) +O(a^3)={-}\frac{2}{{\rm \pi}^2 h} a\cos\bar{\sigma} t +O(a^3). \end{equation}
This means that the theoretical steady-state phase lag and wave amplitude in experiments by Bäuerlein & Avila (Reference Bäuerlein and Avila2021) could in the lowest-order approximation be associated with 
$\theta$ and 
$2a/({\rm \pi} ^2 h)$ coming from the asymptotic solution (2.10)–(2.12).
The experimental model tests by Bäuerlein & Avila (Reference Bäuerlein and Avila2021) were done with the non-dimensional liquid depth 
$h=0.8$ that implies (Faltinsen & Timokha Reference Faltinsen and Timokha2009, table 9.1) 
$d_2=3.142, d_2=2.533, d_3=-0.021, d_4=-0.042, d_5=-3.225$. The experimental forcing amplitudes were 
$\eta _{2a}= 0.0009$, 0.0017, 0.0032 and 0.0064. The low-viscous (fresh water) liquid with the kinematic viscosity 
$\nu = 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$ is used. At the same time, visual observations with 
$\eta _{2a}=0.0064$ discovered wave breaking, overturning and serious contribution of higher modes/harmonics that indicates inapplicability of the single-dominant modal theory for the experimental series with the largest excitation amplitude.
Utilising experimental series with 
$\eta _{2a}=0.009$ (I), 0.0017 (II) and 0.0032 (III) in the above-proposed machine learning procedure makes it possible to compute 
$\xi _1$, 
$\xi _2$ and 
$\xi$. The theoretical accuracy should statistically depend on the accuracy of the measured data. Ideally, if all the measurements 
${\rm I}+{\rm II}+{\rm III}$ are used for training, that computes 
$\xi _1=\xi _1^{l}=0.007167965$, 
$\xi _2=\xi _2^{l}=0.01002107$ and 
$\xi =\xi ^{l}=0.2636$. Other acceptable trainings were: 
${\rm I}+{\rm III} \rightarrow 0.00710169946$, 0.0095860284, 0.193322; 
${\rm I}+{\rm II} \rightarrow 0.00733898346$, 0.010288244, 0.0931; 
${\rm III} \rightarrow 0.0073849521$, 0.005228571, 0.37915. Inserting these values into expressions (2.10)–(2.12) makes it possible to draw the theoretical phase-lag response curves to be compared with measurements adopted for both training and testing. Our numerical experiments showed that the training strategies III, 
${\rm I}+{\rm II}$, 
${\rm I}+{\rm III}$ and 
${\rm I}+{\rm II}+{\rm III}$ lead to visually almost identical curves. The triad 
$\xi _1=\xi _1^{l}=0.007167965$, 
$\xi _2=\xi _2^{l}=0.01002107$, 
$\xi =\xi ^{l}=0.2636$ with 
${\rm I}+{\rm II}+{\rm III}$ is therefore adopted in the graphical analysis.
Figure 1(a) confirms applicability of the single-dominant modal theory (2.1). Agreement is excellent, especially, in contrast to Bäuerlein & Avila (Reference Bäuerlein and Avila2021, figure 14) who adopted 
$\xi _1=\xi _1^{exp}= 8.4\times 10^{-3}$ (coming from the measured logarithmic decrements) and 
$\xi _2=\xi =0$ in their computations by the single-dominant modal system. Agreement in figure 1(a) is also much better than with using the Duffing mathematical model by Bäuerlein & Avila (Reference Bäuerlein and Avila2021, figure 13). As matter of the fact, we showed that damping of the second-order mode (non-zero 
$\xi _2$ in (2.1b)) and the quadratic damping 
$\varXi _1$ in the first modal equation (2.1a) matter.
The theoretical (lines) and experimental (symbols) phase-lag (in grads) response curves for the experimental set-up by Bäuerlein & Avila (Reference Bäuerlein and Avila2021). The measured phase lags with the forcing amplitudes 
$\eta _{2a}=0.0009$, 0.0017 and 0.0032 from the (a) were used in the machine learning procedure from § 3 to compute 
$\xi _1=\xi _1^{l}=0.007167965$, 
$\xi _2=\xi _2^{l}=0.01002107$ and 
$\xi =\xi ^{l}=0.2636$ and substitute them into (2.10)–(2.12). Panel (b) is drawn for the largest experimental forcing amplitude 
$\eta _{2a}=0.0064$ when the single-dominant modal system is, generally, not applicable; it detects a zone where this happens.
Even though, as we remarked above, the single-dominant modal theory is not applicable for the largest experimental forcing amplitude 
$\eta _{2a}= 0.0064$, the computed 
$\xi _1=\xi _1^l$, 
$\xi _2=\xi _2^l$, 
$\xi =\xi ^l$ and 
$\epsilon =P_1\eta _{2a}=P_1\times 0.0064$ were used in (2.10)–(2.12) to draw the theoretical phase-lag response curve and compare it with the corresponding measurements by Bäuerlein & Avila (Reference Bäuerlein and Avila2021). The numerical results are shown in figure 1(b). The figure detects a zone where the single-dominant modal theory is clearly non-applicable (wave breaking and other free-surface phenomena indeed matter!). Far from this zone, agreement between the learnt solution and experiments is also excellent.
We used asymptotic solution (2.10)–(2.12) with 
$\xi _1=\xi _1^{l}=0.007167965$, 
$\xi _2=\xi _2^{l}=0.01002107$ and 
$\xi =\xi ^{l}=0.2636$ to compare the theoretical mass-centre amplitude by (3.3) with measurements by Bäuerlein & Avila (Reference Bäuerlein and Avila2021). The results are presented in figure 2. Figure 2(a) is drawn for the cases from figure 1(a) but figure 2(b) corresponds to 
$\eta _{2a}=0.0064$. The dashed lines in figure 2(a) imply the inviscid single-dominant theory for 
$\eta _{2a}=0.0009$ and 0.0032.
The same as in figure 1 but for the liquid-mass centre amplitude. The fully inviscid branching by the single-dominant theory is shown in (a) by the dashed lines for 
$\eta _{2a}=0.0009$ and 0.0032 to confirm that it cannot predict horizontal position of the jump-down bifurcation point 
$P$, which is a function of viscous damping. For the smaller forcing, 
$\eta _{2a}=0.0009$, discrepancy between viscous and inviscid theoretical amplitude near the amplitude peak is huge while, for the larger one, 
$\eta _{2a}=0.0032$, it does not look dramatical.
Damping is responsible for the horizontal position of the jump-down bifurcation point 
$P$ (demonstrated for the Duffing-type systems, e.g. by Faltinsen & Timokha Reference Faltinsen and Timokha2009, § 8.2.1.3). As for the horizontal position of 
$P$, agreement between the ‘learnt viscous’ solution and measurements by Bäuerlein & Avila (Reference Bäuerlein and Avila2021) is rather good in figure 2(a) but the single-dominant theory is not applicable in the case in figure 2(b). Except in the left of 
$P$ on the upper subbranch, discrepancy between theoretical viscous and inviscid amplitudes does not look dramatic for the larger forcing amplitude 
$\eta _{2a}=0.0032$, whereas both theories are rather inaccurate to fit the measurements. The situation is opposite for the lowest tested forcing amplitude 
$\eta _{2a}=0.0009$. The latter facts confirm earlier conclusions by Faltinsen & Timokha (Reference Faltinsen and Timokha2001) who showed for clean tanks that viscous damping may cause almost negligible effect on the wave amplitude with increasing excitation amplitude when, contrarily, the free-surface nonlinearity and associated energy flow between the generalised hydrodynamic coordinates (modes) play dominant roles. Faltinsen & Timokha (Reference Faltinsen and Timokha2001) showed that agreement for the wave-amplitude response curves with increasing forcing amplitude can be improved when using an adaptive modal theory with zero damping or small damping following from the single-dominant prediction.
3.2. Linear damping ratios $\xi _1$ and $\xi _2$
Neglecting nonlinear terms in (2.1) leads to an infinite set of linear oscillators, equations of the so-called linear modal theory by Faltinsen & Timokha (Reference Faltinsen and Timokha2009, chap. 5), which is applicable to non-resonant sloshing of a low-viscous liquid. When posing 
$\eta _{2a}=0$, its solution determines a superposition of decaying standing Stokes waves whose logarithmic decrements are associated with the linear damping ratios 
$\xi _i$.
Rigorous linear mathematical theory of viscous unforced sloshing was created by Selim Krein (Reference Krein1964). Krein's theorem states that only a finite number of oscillatory sloshing modes exists, i.e. interpreting his result in terms of the linear modal theory with 
$\eta _{2a}=0$, there exists 
$N$ such that 
$\xi _i>\bar {\sigma }_i$ for 
$i\geqslant N$. Numerical analysis by Barnyak & Barnyak (Reference Barnyak and Barnyak1996) also shows that only lower sloshing modes and frequencies of viscous liquids can be approximated within the framework of the inviscid hydrodynamic theory. The inequality and computations are in conflict with assumptions of the linear damped modal theory, which requires small 
$\xi _i$ to weakly affect the natural sloshing frequencies and modes, namely, the theory requires 
$\xi _i\ll \bar {\sigma }_i$ on an asymptotic scale. The latter can only be fulfilled for a few lower natural sloshing modes. Now, readers may understand why we postulated 
$O(\epsilon ^{1/3})=\xi _i\ll \bar {\sigma }_i=O(1)$ only for 
$i=1$, 
$2$ and 3.
When associating 
$\xi _i\ll 1$ with logarithmic decrements of lower natural sloshing modes, these can be estimated from below (Faltinsen & Timokha Reference Faltinsen and Timokha2009, chap. 6) by using asymptotic formulas, which are based on the laminar boundary-layer theory for low-viscous liquids and accounting for the bulk damping effects, 
$\xi _i^{(0)}=\xi _i^{layer}+\xi _i^{bulk}\lesssim \xi _i$. Miles (Reference Miles1967), Miles & Henderson (Reference Miles and Henderson1998) and Faltinsen & Timokha (Reference Faltinsen and Timokha2009, chap. 6) discussed why and when the bulk damping may matter as well as how other potentially valuable viscous mechanisms including the dynamic contact angle, wave breaking, free-surface contamination and roof impact can significantly increase the damping. Taking Keulegan's estimate for 
$\xi _i^{layer}$ and (6.139) by Faltinsen & Timokha (Reference Faltinsen and Timokha2009) for 
$\xi _i^{bulk}$ gives
\begin{equation} \xi_i^{(0)}=\sqrt{\frac{\nu}{2\sigma_i}}\frac 1 {B}\left(1+ \frac{B}{L}\left[1 +{\rm \pi} i\frac{1-2h}{\sinh (2{\rm \pi} i h)} \right] \right)+ 2\nu \frac{{\rm \pi}^2i^2}{L^2\sigma_i^2}, \end{equation}
where 
$\nu$ is the kinematic viscosity. Specifically, the formula provides 
$\xi _i^0/\bar {\sigma }_i\to \infty$ as 
$i\to \infty$ but, as we stated above, (3.4) is valid only for lower indices 
$i$ and strongly requires 
$\xi _i^{(0)}\ll O(1)$.
For the experimental set-up by Bäuerlein & Avila (Reference Bäuerlein and Avila2021, 
$L=0.5$ m, 
$B=0.05$ m and 
$h=0.8$) with fresh water, (3.4) outputs 
$\xi _1^{(0)}=0.0057$ and 
$\xi _2^{(0)}=0.0047$; the computed values are practically not affected by the bulk damping and, therefore, one can say that these imply the Keulegan's approximation. The Keulegan's damping ratios are definitely lower than 
$\xi _1^l=0.007167965$ and 
$\xi _2^l=0.01002107$ coming from the machine learning procedure with 
$\eta _{2a}= 0.0009$, 
$0.0017$ and 0.0032. The inequality is correct because (3.4) estimates from below.
Bäuerlein & Avila (Reference Bäuerlein and Avila2021) conducted special model tests attempting to estimate 
$\xi _1$. They interrupted harmonic excitations of the steady-state wave with 
$\eta _{2a}= 0.0064$ at an instant 
$t=t_0$ and, thereafter, measured logarithmic decrements of the wave decay. The averaged experimental value of 
$\xi _1$ coming from the measured decrements was estimated at 
$\xi _1=\xi _1^{exp}=0.0084$. This value is larger than 
$\xi _1^l=0.007167965$ because of, most probably, very specific free-surface phenomena discovered for model tests with 
$\eta _{2a}= 0.0064$ that could affect logarithmic decrements even after interrupting the harmonic forcing.
On the other hand, the theoretical linear damping rate of the second natural sloshing mode, 
$\xi _2^{(0)}=0.0047$, is too low with respect to 
$\xi _2^l=0.01002107$. The discrepancy between 
$\xi _2^{(0)}$ and 
$\xi _2^l$ may be explained by an ‘inaccurate’ modelling of energy flow between the second and higher modes. The single-dominant modal theory simply ignores this. Using an adaptive modal system by Faltinsen & Timokha (Reference Faltinsen and Timokha2001) would lead to a more accurate distribution of 
$\xi _i$, 
$i\geqslant 2$ and improve agreement in figure 2(a). However, the adaptive modal systems do not have analytical solutions, that is, the machine learning procedure becomes then much more complicated.
4. Conclusions
Estimates of viscous damping in liquid sloshing dynamics can be done by applying a machine learning procedure to multimodal systems with unknown a priori viscous damping terms. Such a procedure was developed for the single-dominant modal system/theory by Faltinsen et al. (Reference Faltinsen, Rognebakke, Lukovsky and Timokha2000) when the steady-state wave (periodic) solution can analytically be derived by using Moiseev's asymptotic scheme. The solution becomes a function of three unknown parameters, 
$\xi _1$, 
$\xi _2$ and 
$\xi$, responsible for damping in the hydrodynamic system. Having known measurements of the phase lags makes it possible to compute (learn) these parameters. The measured phase lags by Bäuerlein & Avila (Reference Bäuerlein and Avila2021) are employed to evaluate abilities of the proposed procedure and find 
$\xi _1=\xi _1^l$, 
$\xi _2=\xi _2^l$ and 
$\xi =\xi ^l$. An excellent agreement is demonstrated except when, according to visual observations in Bäuerlein & Avila (Reference Bäuerlein and Avila2021), the single-dominant system is not applicable. Following the physics-driven machine learning instead of the data-driven approach by Cenedese et al. (Reference Cenedese, Axås, Bäuerlein, Avila and Haller2022) did not provide good agreement for all the steady-state amplitude parameters but made it possible to make several physical conclusions on nonlinear sloshing. The obtained results can also be generalised to transient waves and coupling with tank motions, including in ship dynamics.
One conclusion is that effect of higher (here, second) natural sloshing modes (nonlinear energy transfer from the primary-excited to the higher modes) and viscous damping of these modes cannot be neglected. Furthermore, viscous damping of the primary-excited (first) natural sloshing mode should, generally speaking, be a function of the wave amplitude and, moreover, for the single-dominant system, the function is just the linear regression. Finally, the proposed procedure can be an efficient tool for estimating the viscous damping of the lowest (dominant) natural sloshing mode, but it fails for the higher (order) modes. To get more accurate 
$\xi _2$, one should use an adaptive multimodal system by Faltinsen & Timokha (Reference Faltinsen and Timokha2001) where 
$\beta _1\sim \beta _2$. Such a system is also required for an accurate prediction of the wave-amplitude response curves.
Derivations of the steady-state wave (periodic) solution does not require an exact analytical expression for 
$\varXi _1(\beta _1,\dot {\beta }_1)$. Any analytical quadratic form of 
$\varXi _1$ is allowed if it deduces the two integrals (2.7) and (2.8). An example could be
\begin{equation} \varXi_1(\beta_1,\dot{\beta}_1)=\xi \frac{3{\rm \pi}}{8} |\dot{\beta}_1|\, \dot{\beta}_1, \end{equation}which physically implies a drag force. Based on model tests by Faltinsen et al. (Reference Faltinsen, Rognebakke, Lukovsky and Timokha2000) and other authors with transient waves, one should check whether (4.1) is indeed applicable.
Funding
A.T. acknowledges financial support of the National Research Foundation of Ukraine (project number 2020.02/0089).
Declaration of interests
The authors report no conflict of interest.
