Solitary-wave loads on a three-dimensional submerged horizontal plate: Numerical computations and comparison with experiments

A parallelized three-dimensional (3D) boundary element method is used to simulate the interaction between an incoming solitary wave and a 3D submerged horizontal plate under the assumption of potential flow. The numerical setup follows closely the setup of laboratory experiments recently performed at Shanghai Jiao Tong University. The numerical results are compared with the experimental results. An overall good agreement is found for the two-dimensional wave elevation, the horizontal force and the vertical force exerted on the plate, and the pitching moment. Even though there are some discrepancies, the comparison shows that a model solving the fully nonlinear potential flow equations with a free surface using a 3D boundary element method can satisfactorily capture the main features of the interaction between nonlinear waves and a submerged horizontal plate.


I. INTRODUCTION
Submerged horizontal plates are common coastal engineering structures that are used for several purposes. When located close to or at the mean water surface, they can act as effective breakwaters for offshore wave control and harbour protection as discussed by Yu [1]. Meanwhile, a pulsating reverse flow may occur below these breakwaters under certain circumstances. Turbines can be put under the plate (Graw [2]; Carter [3]) to convert wave energy. A lot of other coastal structures, such as bridges, docks and very large floating structures, can be modeled as submerged plates in order to study the effects of storm surges, tsunamis and other extreme wave events.
Periodic wave scattering by submerged plates has been widely studied. Siew and Hurley [4] provided first-order reflection and transmission coefficients for the scattering of long waves by a submerged plate. Patarapanich [5] studied forces and moments exerted on plates both experimentally and numerically using a finite element method. The effects of various parameters such as the ratio of plate length to wavelength and the submergence depth were investigated as well. Cheong et al. [6] extended the eigenfunction expansion method to the complete range of water depths and compared reflection and transmission results with finite-element simulations. Dong et al. [7] used a modified matched eigenfunction expan-sion method to analyse wave scattering on a submerged horizontal plate over variable topography.
Due to their simplicity for experimental studies, solitary waves have been used extensively by researchers. Lo and Liu [8] conducted shallow water experiments of solitary waves incident on a submerged plate. Experimental results were compared with numerical simulations and analytical solutions based on linear long wave theory. Strong vortices were observed near the trailing edges using particle image velocimetry (PIV). Seiffert et al. [9] conducted a series of laboratory experiments to investigate the forces exerted on a submerged plate by solitary waves. Hayatdavoodi and Ertekin [10] studied waveinduced loads due to solitary and cnoidal waves using Green-Naghdi theory and the influence of several parameters was discussed. Dong et al. [11] conducted experiments and simulated solitary-wave interactions with a submerged horizontal plate both on a flat bottom and on a sloping beach. Christou et al. [12] studied the influence of the angle of attack when a solitary wave propagates over a thin finite square plate. They used Hydro3D, an open source Large Eddy Simulation code. Xie et al. [13] used a multiphase flow model combined with the largeeddy simulation approach to investigate the interaction of a solitary wave with a thin submerged plate. Wang et al. [14] performed three-dimensional (3D) experiments and measured the spatial and temporal variation of the twodimensional (2D) free-surface deformation using a multilens stereo reconstruction system. The hydrodynamic loads were measured by underwater load cells. Wave focusing induced by the plate led to an increased maximum elevation along the streamwise centerline of the plate. A 6-stage loading process based on the maxima of the vertical wave force and the pitching moment was proposed.
One of the conclusions is that the vertical wave force on the plate is reduced compared to that obtained in previous 2D experiments. Although strong vortices were observed at the trailing edges of the plate, it is legitimate to ask the following question: can this problem only be solved by Computational Fluid Dynamics (CFD) or can fully nonlinear potential flow theory still be applied to this problem? Of course, it depends on the Reynolds number. In Xie et al. [13], the Reynolds number based on the wave speed and plate thickness is approximately equal to 10 4 , which justifies the use of CFD methods. In the experiments of Wang et al. [14], the Reynolds number is one order of magnitude larger (≈ 10 5 ) because of a thicker plate and a larger water depth.
In the present paper, we first describe the numerical method in Section 2. The laboratory experiment is reviewed in Section 3. Numerical results for the wave elevation, the vertical force and the moment are provided in Section 4. They are compared with experimental results. Velocity fields are also shown. The effect of vortices that is neglected in the numerical simulations is discussed. In Section 5, we investigate the influence of the thickness of the plate. Additional results are provided as Supplementary Material.

II. NUMERICAL METHOD
The fully nonlinear potential flow model with a free surface is used to solve the problem of a solitary wave impacting on a submerged horizontal plate. The fluid domain is denoted by Ω, with boundary Γ. The boundary includes the free surface, the wavemaker, the bottom, the submerged plate and a vertical wall far downstream of the plate.

A. Mathematical formulation
The velocity potential φ(x, t), where x = (x, y, z) is the vector of spatial coordinates with z the vertical coordinate and t is the time, is used to represent inviscid irrotational flows. The continuity equation in the fluid domain is Laplace's equation for φ: We follow the approach described in Grilli et al.'s [15]. The three-dimensional free space Green's function is defined as where r = x − x l with r = |r| being the distance from the source point x to the collocation point x l , and n representing the normal unit vector pointing out of the domain at point x.
Green's second identity transforms Laplace's equation (1) into a integral equation on the boundary: where α(x l ) is proportional to the exterior solid angle made by the boundary at the collocation point x l .
On the free surface, φ satisfies the nonlinear kinematic and dynamic boundary conditions, written in a mixed Eulerian-Lagrangian form, with the material derivative D/Dt ≡ ∂/∂t + ∇φ · ∇: where x is the position vector of a free-surface fluid particle and g the acceleration due to gravity. The atmospheric pressure has been set equal to 0. In the case of wave generation by a wavemaker moving with velocity U , the normal velocity is continuous over the surface of the wavemaker: At the bottom Γ b (t) and along other fixed parts of the boundary, the no-flow condition ∂φ/∂n = 0 is prescribed.

B. Time integration
Following the method implemented in Grilli et al.'s [16] 3D model, second-order explicit Taylor series expansions are used to express both the new position and potential on the free surface. First-order coefficients are given by the boundary conditions (4) and (5). The pairs ∂φ/∂t, ∂ 2 φ/∂t∂n that are needed to obtain secondorder coefficients are computed by solving another integral equation similar to equation (3). For the evaluation of the tangential derivatives, a fourth-order interpolation scheme is employed.
The time-step is adapted by finding the minimum distance between two nodes on the free surface. Grilli et al. [17] found an optimal value for the constant Courant number C 0 of roughly 0.4. In order to maintain the stability when strong nonlinear free surface deformations occur, an equally-spaced regridding method is adopted every 10 time steps, starting when the crest of the solitary wave arrives at the front edge of the plate. Lagrangian points would otherwise concentrate and eventually lead to a crash of the computations. In the literature, researchers use similar smoothing techniques to remove instabilities. For instance, Longuet-Higgins et al. [18] used a filter every 5, 10 or 20 time-steps. Ming Xue et al. [19] applied a similar technique every N s (N s typically 3 or 6) steps. Grilli et al. [16] and Fochesato and Dias [20] also used a free surface node regridding method.

C. Spatial discretization
In this subsection, we follow closely Fochesato and Dias [20]. The boundary is discretized into N collocation nodes and M high-order elements are used to interpolate in between q of these nodes. Within each element, the boundary geometry and the field variables are discretized using polynomial shape functions N j (ξ, η), where ξ denotes the intrinsic coordinates of the reference element: where the indices j = 1, . . . , q locally number the nodes within each element. We choose cubic reference elements (q = 16), which provides C 2 continuity in between elements.
The integral equation (3) is transformed into a sum of integrals over the boundary elements, each one being calculated on the reference element. The change of variables is described by the Jacobian of the transformation J i for the i th element. Thus, the discretized form of the integrals can be written as The associated discretized boundary integral equation leads to a sum on the N boundary nodes, where l = 1, . . . , N and K D lj , resp. K N lj , are Dirichlet, resp. Neumann, global matrices.
When the collocation node l doesn't belong to the integrated element, a standard Gauss-Legendre quadrature is used. When it does belong to the element, r becomes zero at one of the nodes and a self-adaptive singular quadrature [21] is implemented to handle the presence of the singularity on the boundary.
Instead of computing the diagonal coefficient of the Neumann matrix K N ll , the rigid mode technique is used: which provides the diagonal term of a row by minus the sum of its off-diagonal coefficients. In terms of the intersecting parts of different boundaries, such as between the free surface and the lateral boundaries, the boundary conditions and the normal directions are generally different. Therefore, double-nodes are used to represent these corners and the continuity of the potential is imposed.

D. Domain decomposition
In order for the simulations to mimic as closely as possible the experiments, domain decomposition [22] [23] [24] is used to boost the efficiency. The fundamental idea of domain decomposition is to divide the computational region into two or more sub-domains. The boundary integral equation is solved independently in each sub-domain. One major issue in domain decomposition is to satisfy continuity between adjacent sub-domains.
Following the so-called D/D-N/N scheme introduced by De Haas et al. [23], the computational domain is decomposed into the sub-domains Ω 1 and Ω 2 which are separated by an interface Γ . On the interface, the potential and its derivative are unknown, and an initial guess needs to be imposed. An iterative procedure is then used to get the exact potential or its derivative on the interface. The properties of Laplace's equation lead to continuous potential values and of their derivatives on the interface. The scheme can be extended straightforwardly to the case of more sub-domains.
The only issue in this iterative scheme is to deal with the case when the interface has a Dirichlet boundary. Because the double nodes share the same geometry and boundary condition between the free surface and the interface, their coefficients in the matrix will be exactly the same resulting in singular algebraic equations. To deal with this difficulty, the so-called semi-discontinuous elements were proposed in [24]. The idea is to use different geometrical points to do the integration based on a discontinuity coefficient γ. It brings extra error in the sim-ulation. As stated in [24], this method provides slightly better conditioned matrix equations, but the convergence is still slow when using an iterative solver. Here, we propose another way to resolve this difficulty.
The singularity occurs because the double nodes share the same Dirichlet boundary condition. Thus if we can somehow change one of the double nodes into a Neumann boundary condition, then the continuity of potential can be imposed again. Although the normal derivatives are unknown on the free surface, in the iterative procedure we can always get an inaccurate solution from the previous step. Therefore, when the interface is of Dirichlet boundary type, we can change one of the double nodes on the free surface into a Neumann boundary condition using the normal derivative from the previous step. This slightly modified scheme has been satisfactorily implemented in the numerical code.
In order to illustrate this new scheme, we first recapitulate the original algorithm [24]: (0) Choose an initial guess φ k on the interface Γ, (k = 0); (i) Solve Laplace's equation in each sub-domain to get ∂φ k 1 /∂n 1 and ∂φ k 2 /∂n 2 on Γ; (ii) Take an average of the solutions (the normal vectors n 1 and n 2 are opposite): (iv) Take an average of the solutions: (v) Calculate the maximum error k+2 on Γ: (vi) If k+2 satisfies a prescribed criterion, then exit the iteration. Otherwise, go to step (1) and repeat with k = k + 2.
The problem mentioned above occurs at step (i) since the interface and the free surface both share Dirichlet boundary conditions. As shown in Fig 1, double nodes are used on the intersection of the free surface and the interface. We can change those double nodes on the free surface into a Neumann boundary condition. In the initial step (0), the normal derivative we need can be obtained from the solution of the previous time step. In step (iii), when we solve Laplace's equation with a Neumann boundary condition on the interface, we get normal derivatives on the free surface. Those values can be used during the iteration.
The main advantage of this domain decomposition method is that it is superlinear -see Table I and Fig  2. The reason is that the assembly of the full matrix is O(N 2 ) and here we use a direct solver which is O(N 3 ). We need to mention that in this scheme a few extra nodes are added on the interface. Hollow circles denote the nodes on the free surface and filled circles represent the nodes on the interface. Here a Dirichlet boundary condition is applied on the interface. Therefore, the double nodes on both surfaces have the same Dirichlet boundary condition, which leads to a singularity. The new scheme resolves this issue by imposing a Neumann boundary condition on the double nodes on the free surface, using the values from the previous time step. The experiments were conducted in the Tsunami Basin for Offshore Regions in Shanghai Jiao Tong University [14]. The wave flume is 42 m long and 4 m wide with a piston type wavemaker installed at one end. The platetype structure is made of organic glass. It is 200 cm long, 78 cm wide, 10 cm thick and mounted on the bottom in the middle region of the flume (Fig 3 and Fig 4). Four piezoelectric force balance units are installed inside the plate. A dynamometric system is used to measure wave loads on the structure. The constant water depth is denoted by h.
The free-surface elevation is measured by resistancetype wave gauges. Twenty wave gauges are spread around the plate. Their exact locations and labels are shown in Fig 5. Taking the center of the plate as the origin (0,0) of the (x, y)−plane, we list the (x, y)−coordinates of all the wave gauges (WG) in Table 1. s

B. Wave generation
Solitary waves are generated by a piston-type wavemaker on the left side of the wave tank. A third-order solitary wave profile [25] is defined by (note that the wave profile is uniform across the flume, so it is given as a func-  tion of x only) where η is the wave elevation, H the wave height, h the still water depth, α = H/h, g the acceleration due to gravity and c the wave celerity. The symbols sech, resp. tanh, denote sech(k(ct − x)), resp. tanh(k(ct − x)).
In reference to the improved Goring et al. [26] wave generation method introduced by Malek-Mohammadi et al. [27], the horizontal velocity of the piston is determined by where 3 . This wave generation method combined with the thirdorder solitary wave profile has the capability of generating steady solitary waves of dimensionless wave amplitude up to H/h = 0.5 [28]. It is implemented in the numerical code.

IV. RESULTS
The size of the computational domain is exactly the same as that of the wave flume: 42 m long and 4 m wide with a moving boundary condition at the left end. At the bottom and along other fixed boundaries, the impermeable condition is applied. The right end is far way from the plate so that the reflected wave has no influence. The discretization used in the simulations is 300 × 30 × 10 elements on the tank boundaries and 25×10×5 on the plate boundaries along the x, y and z directions, respectively (i.e. 25450 elements and 26982 nodes). The numbers of elements on the tank boundaries are chosen based on the finest mesh used in [16]. Additional convergence tests, shown in Fig 6, have been performed to check the level of refinement needed along the plate. As can be seen, refining the mesh further has a negligible effect on the hydrodynamic load. In Xie et al. [13], where a CFD code was used, the computational domain is discretised by a uniform mesh 1600 × 96 × 160.
As shown in Fig 4, G is the submerged depth and B the distance from the bottom to the lower surface of the plate. For the comparison between experiments and numerical simulations, 41 cases have been investigated. They correspond to various combinations of water depth, wave height and submerged depth that are listed in Table 3. In five of these combinations, breaking waves were observed in the experiments. We first compare the numerical results for the surface elevations, horizontal and vertical forces, and pitching moments with the experimental results. Then we present results for the velocity fields. Finally, the basic hydrodynamic loading process that occurs when the solitary wave propagates past the plate is analyzed. Numerical results for all cases can be found in the Supplementary Material.   Therefore they are not shown here.
To show all the cases at all the wave gauges would take too much space. Instead, we decided to show first the differences between the wave gauges for a given experimental setup -see Fig 7 -and then the differences between the various cases at a given wave gauge -see Fig 8. Throughout the discussion, the numerical and experimental results are synchronized using the data from WG3, which is located on the middle line upstream of the plate.
The experimental case we selected for the comparison at all wave gauges has the following characteristics: B = 30 cm, h = 50 cm, H/h = 0.5. It corresponds to a large amplitude non-breaking wave (H/h = 0.5) with a plate which is relatively close to the free surface (G = 10 cm). Overall, the numerical results shown in Fig 7 compare well with the experimental data, both for the peak and for the oscillations that follow (the oscillations are clearly visible at WG11 for example). The wave height increases from the lateral side to the middle line. The highest peak occurs at WG15 when the solitary wave leaves the plate. This is due to wave shoaling.
By looking at the first three wave gauges WG1, WG2 and WG3, we notice that the wave reflected by the leading edge of the plate is quite small compared to the large deformations at WG11 and WG15. This is different from what happens in 2D experiments [8], where the reflected wave is clearly observed.
Checking the three wave gauges WG7, WG10 and WG15 that surround WG11 in Fig 7, it can be seen that the oscillations are noticeable at all three gauges. This indicates that the large deformation above the plate propagates in all directions. Fig 8 shows the free surface elevation for all selected cases at WG11, which is above the center of the plate. It can be seen that when the plate is closer to the free surface, the free surface deforms more. For example, after the solitary wave passes over the plate, the oscillations that appear are larger when B = 40 cm than when B = 20 cm.
Generally speaking, as the wave amplitude becomes larger, the agreement between experimental and numerical results becomes better. This is partly due to the electrical noise in the experiments. For large amplitude waves, the influence of electrical noise becomes negligible. For small amplitude waves, the numerical peak is smaller than the experimental peak. Fig 9 shows the overall free-surface deformation above the plate for the case B = 40 cm, h = 60 cm and H/h = 0.3. As the solitary wave passes above the plate, the wave amplitude becomes larger due to shoaling. Once the large amplitude wave reaches its maximum amplitude, it propagates in all directions, causing reflection from the trailing edge. Meanwhile, the reflection also occurs at the leading edge. These two effects compete simultaneously and leave a pitfall just above the plate. The large amplitude wave keeps radiating and overcomes the pitfall leaving a bulge above the plate.

B. Horizontal force, vertical force and moment
As anticipated because of the symmetries involved in the experimental and numerical setups, it has been confirmed both by [14] and in the present work that the lateral force F y , the yaw moment M z and the roll moment M x are negligible. Therefore we concentrate on the horizontal force F x , the vertical force F z and the (pitching) moment M y for comparisons between the various cases. The basic behavior of F x , F z and M y was explained in [14], where a 6−stage process was introduced to highlight the various peaks in the vertical force and moment.
As said above, we ran a lot of cases. In Fig 10, we show the results for two cases which we believe are representative of all cases. In these two cases, all parameters are identical except the depth of submergence of the plate. A third case is shown in Fig 11. In terms of horizontal and vertical forces, the numerical results agree well with the experimental data, even if the first peak is systematically under predicted by the numerical code. The agreement for the moment is not as good: after some time, the numerical values deviate from the experimental values. The greater the free-surface deformation, the longer the moment agrees (this is shown in the Supplementary Material). The moment is dominated by the pressure distribution on the upper and lower surfaces of the plate. What is intriguing in the experimental data is that the moment decreases first and then amplifies as seen for example in Fig 10. In the subplot at the lower left, the experimental moment is very small for t between 12 and 14 s but oscillations of increasing amplitude appear for t > 14 s. In the subplot at the lower right, the experimental moment develops oscillations of increasing amplitude for t > 15 s. The match between experimental and numerical values is poor and remains unexplained at this stage. For the horizontal force, as observed in both cases, the plate first experiences a positive force and then a negative force. This sequence occurs when the solitary wave hits the leading edge and then leaves the trailing edge. The basic structure of the vertical force has been explained in [8]. It is of interest to explore the oscillations observed in the hydrodynamic loads as well as in the wave gauge data. Their connection will be explained in the next section.
As shown in Fig 11, where results for a different set of parameter values are presented, we can define the positive maximum horizontal force as f x + and the negative minimum force as f x − . The first peak values of the vertical force and moment are defined as f z + 1 and M y + 1 . The second peaks are defined as f z + 2 and M y + 2 . The negative minima are defined as f z − and M y − . As can be seen in Figs 12, 13 and 14, the discrepancy between experimental and numerical values increases as the wave amplitude increases. As explained in [8,29], the discrepancy is most likely due to the presence of a boundary layer and vortex shedding along the plate, which alter the pressure distribution on the plate. In general, the numerical results capture the trend of these extreme values. Note that the purpose of the experiments of [14] was not to detect boundary layer effects. The flow near the four edges is obviously more complicated than the rest of the flow. Strong vortices can form near these edges.
In the following figures (Figs 15 and 16), we show the 2D free-surface elevation at four selected times. forms above the plate due to shoaling, which generates great dynamic pressure on the upper surface at time B, and leads to the first negative vertical force. Then at time C the solitary wave just passes over the plate. The pressure at the trailing edge is larger than that at the front edge. The large negative dynamic pressure on the upper surface leads to the second positive vertical force. Finally, at time D, the second bulge forms above the plate and generates positive dynamic pressure on the upper surface. Though the dynamic pressure is now much smaller than that at time B, it still leads to the second negative vertical force because the dynamic pressure on the lower surface becomes negligible. A discussion on the pressure distribution on a submerged circular plate due to a solitary wave can be found in Wu et al. [30]. We have partially explained the process of oscillations in section 4.1 that covers the above figures. Another important factor for these sustained oscillations is the reflection from the lateral walls. To explore that, we changed the width of the wave tank. The results for the vertical force are shown in Fig 19. When the flume width is multiplied by 2, the oscillations have a smaller amplitude. It is surprising that the frequency does not change much. With four times the width, a similar structure persists up to the third peak of the vertical force f z + 5 . It looks like the magnitude of f z + 5 is not affected by an even larger width, since with doubling the width from 8 to 16, the value only decreases a little bit.

C. Velocity field
The numerical code used in the present study can also be used to compute velocities inside the fluid domain.  The flow from the left is due to the reflection from the upper surface of the plate. Meanwhile, the flow from the right is caused by the propagation of the bulge mentioned above. Once the focused wave reaches its peak, it spreads again.

V. INFLUENCE OF THE THICKNESS OF THE PLATE
Having shown the validity of the numerical method, several purely numerical experiments are conducted to study the effects of the thickness of the plate. As shown in Fig 4, once we change the thickness we may also change G and B, which are the submerged depth and the distance between the plate and the bottom, respectively. In the following discussion, we decided to fix the water depth once for all to h = 60 cm and the dimensionless wave height to H/h = 0.3. We first take G = 20 cm and consider three different cases for the thickness: δ = 10 cm, δ = 20 cm and δ = 30 cm.
The horizontal and vertical forces are shown in Fig  22. The vertical force only increases slightly with the thickness. The horizontal force increases with the thickness, which is not surprising. We can normalize the horizontal force with the thickness as shown in Fig 23. We also checked that the free-surface elevations at three wave gauges along the middle line are the same within graphical accuracy. Therefore we can conclude that under the same submerged depth G, the thickness of the plate has negligible effects on both the hydrodynamic loads and the free-surface elevation around the plate. For the latter one, Lo and Liu [8] have shown that in the two-dimensional case varying the thickness changes the shapes of both the reflected and transmitted waves.      the difference in depth of submergence, we computed the vertical force with G = 10 cm, B = 40 cm, δ = 10 cm and G = 10 cm, B = 20 cm, δ = 30 cm. We found that the main feature of the vertical force is maintained with different thicknesses under the same depth of submergence. Therefore, the characteristics of the hydrodynamic loads and wave elevations around the plate are dominated by the depth of submergence.

VI. CONCLUDING REMARKS
The interaction between a solitary wave and a fully submerged three-dimensional horizontal plate is investigated numerically and the results are compared with those of three-dimensional laboratory experiments described in [14]. The wave focusing phenomenon that was observed experimentally is also observed numerically. The free-surface elevation increases gradually along the center line of the plate and exceeds the incoming wave amplitude. The presence of the three-dimensional plate modifies the shape of the solitary wave. A wave amplification is produced by the local shoaling effect of the plate and the shoaling-induced wave refraction. A larger amplitude of the solitary wave leads to a stronger focusing up to the end of the plate.
The horizontal wave force is characterized by a peak followed by a trough. The vertical wave force is characterized by a series of peaks and troughs. The loading process is described based on the peaks of the vertical force and the pitching moment. The process is linked with the oscillations of the free-surface elevation. When the wave approaches or leaves the plate, the channel flow under the plate contributes to the positive peaks of the vertical force. As the wave crest approaches the center of the plate, the negative force caused by the dynamic pressure on the top side of the plate dominates. The strong focusing of large-amplitude waves reduces the second positive peak. The pitching moment is mainly generated by the vertical force, and the observation of surface elevation indicates that the time for the peak of pitching moment depends on the occurrence of the maximum surface ele-

vation.
Then the velocity field is illustrated. Unfortunately, there are no experimental results to compare it with. Finally, the influence of the plate thickness is found to be negligible.
In the introduction, we asked the question: can the problem of a solitary wave impacting on a submerged horizontal plate only be solved by Computational Fluid Dynamics (CFD) or can fully nonlinear potential flow theory still be applied to this problem? In the experiments of Wang et al. [14] that we used for the comparisons, the Reynolds number was of the order of 10 5 . Overall, the potential flow model gives good results. However, there are some discrepancies, especially for the pitching  moment. In the future, it is suggested to compare the experimental results also with a Navier-Stokes solver, especially when breaking waves have been observed in the experiments. The present problem could become an excellent benchmark to study the performance of codes that can handle wave breaking.

SUPPLEMENTARY MATERIAL
Results for all the cases are available in the supplementary material.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.