From Bore–Soliton–Splash to a New Wave-to-Wire Wave-Energy Model

We explore extreme nonlinear water-wave amplification in a contraction or, analogously, wave amplification in crossing seas. The latter case can lead to extreme or rogue-wave formation at sea. First, amplification of a solitary-water-wave compound running into a contraction is disseminated experimentally in a wave tank. Maximum amplification in our bore–soliton–splash observed is circa tenfold. Subsequently, we summarise some nonlinear and numerical modelling approaches, validated for amplifying, contracting waves. These amplification phenomena observed have led us to develop a novel wave-energy device with wave amplification in a contraction used to enhance wave-activated buoy motion and magnetically induced energy generation. An experimental proof-of-principle shows that our wave-energy device works. Most importantly, we develop a novel wave-to-wire mathematical model of the combined wave hydrodynamics, wave-activated buoy motion and electric power generation by magnetic induction, from first principles, satisfying one grand variational principle in its conservative limit. Wave and buoy dynamics are coupled via a Lagrange multiplier, which boundary value at the waterline is in a subtle way solved explicitly by imposing incompressibility in a weak sense. Dissipative features, such as electrical wire resistance and nonlinear LED loads, are added a posteriori. New is also the intricate and compatible finite-element space–time discretisation of the linearised dynamics, guaranteeing numerical stability and the correct energy transfer between the three subsystems. Preliminary simulations of our simplified and linearised wave-energy model are encouraging and involve a first study of the resonant behaviour and parameter dependence of the device.


Introduction
Early September 2010, three applied mathematicians at the University of Twente made requests to create a soliton in a make-shift wave tank for a new "research plaza" opening festivity. Part of that plaza contains a water feature or channel approximately 45 m long, 2 m wide, and 1.2 m deep. Normally filled to its edge with water and harbouring water plants and fish, at the time, it was only partially filled. A soliton is a wave with nonlinearity and dispersion in balance, such that the wave stays coherent and neither disperses nor breaks [13,33]. Solitons or solitary waves can be generated at the beginning of a rectangular channel with vertical walls: using either a piston moving bespokely, a block lowered at a finite yet fast speed into the water or by a quick sluice-gate removal between a higher rest-water level (h 1 ) lock section and a lower rest-water level (h 0 ) main section. We have used the latter for solitary-wave generation with an extra channel feature, sketched in Fig. 1 with dimensions in Table 1: a V-shaped channel end with vertical walls.
While a soliton is a well-known mathematical and fluid-mechanical feature in nonlinear science, a general audience would be entertained more when such a travelling heap of water would lead to an extreme wave, to mark the plaza opening. We, therefore, added the V-shaped contraction to create the highest splash possible by only varying water levels h 0 and h 1 before and after the sluice gate, given a wave-channel geometry. Ideas for wave-height amplification arose from work on hydraulic flow stowage in contractions [2] and on wave impact against sea walls by Peregrine [6,44], cf. Fig. 2. After some trials in two wave tanks, we created a "bore-soliton-splash" 1 . The highest bore-soliton-splash-case 8 in Table 2-consists of time-dependent evolution summarised in Fig. 3: (a) water initially at rest with (b) an excavator ready to lift the sluice gate out of the channel; after sluice-gate removal, a coherent compound of "2.5" solitary waves travels towards the contraction and the highest solitary wave breaks shortly after its inception into (c) a so-called hydraulic bore or spilling breaker; the slightly lower non-breaking second wave is tailed by a lower third wave, with all three waves travelling at a slightly different speeds; during its propagation, the bore diminishes in height due to turbulent energy dissipation and (d) breaking ceases just before the contraction; the smoothened wave then amplifies and reflects in the V-shaped channel contraction and upon reflection draws a deep trough in which the second slightly lower wave crashes precisely, leading to (e, f) generation of a wave or splash circa ten times the incoming, first, and highest solitary-wave height. The bore-soliton-splash let to a variety of new ideas including its modelling and also the creation and testing of a proof-of-concept of an inspired wave-energy device. Our paper aims to highlight and partially investigate these ideas and to relate this splash or man-made rogue wave to similar rogue or monster waves in our oceans.
Rogue, monster, or extreme waves are anomalously high and rare waves with wave height H rw , generally considered at sea, defined relative to a significant or ambientwave height H s of surrounding, preceding, and following seas. A straightforward definition of rogue waves states that the abnormality index AI = H rw /H s > 2, i.e., the   [14] and Khariff et al. [31] provide more advanced and precise definitions of rogue waves, but this common definition given above suffices here. Rogue waves have a rare, extreme occurrence and are, therefore, difficult to predict, either statistically or deterministically. Understanding their wave height and occurrence is relevant to maritime and coastal engineering given their potential to damage ships, maritime, and coastal structures, including sinking and disappearance of ships; an overview of such disasters is found in Nikolkina and Didulenkova [41]. There are different types of rogue waves, involving linear and nonlinear wave focussing in one horizontal dimension, spatial wave focussing due to coastal or submarine convergences, episodic   Fig. 3b taught before gate removal at circa 2.5 m/s. Cases 3, 6, and 8 were nearly identical, as evidenced by video footage [50], with (slight) differences attributed to differences in rest levels h 0,1 and sluice-gate removal operation waves such as tsunamis generated elsewhere, and crossing seas with pyramidal waves [7,16]. These different rogue-wave types have been (partially) explained within a hierarchy of different models, including, e.g., incompressible Euler equations with a free surface and passive or limited air motion, its potential-flow restriction, and numerous asymptotic approximations of these classical potential-flow water-wave equations such as Benney-Luke equations, Kadomtsev-Petviasvili's (KP) equation, nonlinear Schrodinger equation(s), and the Korteweg-De-Vries equation [13,28,31,40,42,43].
The main results achieved in this article are: • a detailed description of a man-made bore-soliton-splash rogue wave with an abnormality index of AI ≈ 10; • establishing and employing mathematical and numerical models for experimental cases 8 and 9 (see Table 2), with improved simulations beyond the one in Bokhove and Kalogirou [4]; and, • inspired by the bore-soliton-splash configuration, we invented a novel roguewave-energy device, and built and tested a scaled-down version; a first nonlinear mathematical model is developed here, for which we show simulations of its linearised dynamics.
The outline of our paper is as follows. Soliton-splash and bore-soliton-splash experiments are analysed in Sect. 2. Some mathematical and numerical solutions of soliton splashes with Benney-Luke's model are found in Sect. 3. In Sect. 4, our wave-energy device is introduced with one comprehensive and novel, nonlinear mathematical wave-to-wire model of the hydrodynamics, buoy motion, and power generation. After developing an intricate and novel compatible discretisation of that linearised model, numerical modelling results are presented. We finish with a discus-  Table 2. a Channel overview before sluice gate is removed. b Sluice gate with excavator used to smoothly remove the sluice gate. After sluice-gate removal, c the highest solitary wave in the compound becomes a bore or spilling breaker, dissipating turbulent energy, and diminishing amplitude, while it propagates to become d smooth again before the contraction. After reflection, it draws a trough in which the second wave falls, thus forming e a jet f collapsing after reaching its apex sion of open questions and challenges in Sect. 5, also highlighting a splash-inspired artwork.

Experimental Set-Up and Results
Our goal in 2010 was to create both a travelling soliton by removing a sluice gate separating two different water levels, initially at rest, and a splash of the highest possible amplitude in a V-shaped contraction. Given time constraints, the only way to determine whether our goal was reachable in practice was to resort to experimentation in two make-shift wave channels: the first one where the Roombeek, a brook, flows onto the University of Twente campus and the second one on the above-mentioned Research Plaza, see Fig. 3. On 20-09-2010, we managed to obtain two soliton splashes with h 1 ≈ 2h 0 in the Roombeek convincing us that it was possible to make a larger and reproducible Plaza-channel soliton-splash. Subsequently, seven test cases were completed on 27-09-2010, including six with different rest-water levels h 0 and h 1 , and one repeated case with the highest splash to ensure reproducibility on the opening day of the Research Plaza. The optimal case involved h 0 = 0.41 ± 0.01 m and h 1 = 0.9 ± 0.01 m. These two repeat cases and the general outcomes on 27-09-2010 showed that our experiments to create a bore-soliton-splash were reproducible on the opening day (30-09-2010). All cases are summarised and dated in Table 2 with repeat cases underlined and numbered by 3, 6, and 8. On 30-09-2010, this "optimal" case was successfully repeated, as shown in Fig. 3, followed by a case numbered 9 with the higher water level of h 0 = 0.43 m set in the main channel by the addition of sluice compartment's extra water from optimal case 8, while keeping h 1 = 0.9 m; case 9 resulted in a smooth solitary-wave compound without wave breaking and a lower splash. Its evolution in time is displayed in Fig. 4 as bespoke simulations introduced and explained later. Videos of (nearly) all cases are found on Zweers' YouTube channel [50] and numbered accordingly. Inspection of videos of three repeat cases 3, 6 and 8 reveals that there are some/minor differences, partially commented on in Table 2. We attribute differences to the estimated error of ±0.01 m in initial water levels h 0,1 and the manner of sluice-gate removal by the excavator, despite training to be as consistent as possible. Case 9 underscores these sensitivities to the initial conditions, because a 0.02 m change from h 0 = 0.41 m (cases 3, 6, and 8) to 0.43 m, while keeping h 1 = 0.9 m, within measurement error, led to a quite different splash. Note, however, that the outcomes are not chaotic, as three reasonably repeatable cases demonstrate.
Several splash types were observed in the nine cases including minor/major reflections, resonances between waves (cases 3, 6, and 8), smooth waves, sheets (case 9), and pyramidal waves (cases 3, 6, and 8). Rogue-wave amplitudes found ranged between H rw ∈ [0.6, 3.5] ± 0.5 m, i.e., concerning errors in H rw of 14-83%. For the highest cases of the bore-soliton-splash, we found H s = 0.35 ± 0.05 m for the highest and first solitary wave and a rogue-wave height of H rw = [3.25, 3.5, 3.5] ± 0.5 m leading to a maximum abnormality index, cf. [10], in the range of:

Fig. 4
Numerical solution of soliton splash event case 9 with μ = 0.04 and = 0.55, see [4,21]. When taking h 0 = 0.43 m instead of h 0 = 0.41 m as in case 8, with h 1 = 0.9 m the same in both cases, no wave breaking occurs [5]. Photo times at t = 8, 14, 15, 15 ± 0.5 s (relative) of observations found in [4] can be compared with simulation times at t = 8, 14, 15, 15.34 s. Values displayed are in metres. The simulation involves N k = 8010 elements of which N x N y = 20 × 390 = 7800 elements lie in the regular part of the channel and N x (N x + 1)/2 = 210 elements in the triangular contraction. There are N n = 8431 nodes with (N x + 1)(N y + 1) = 8211 nodes in the regular part of the channel and N x (N x + 2)/2 = 220 nodes in the triangular contraction truly rogish compared to observations of a typical AI ∈ [2,3] in the oceans. Coastal rogue waves can have a larger abnormality index, cf. [41] and the Tohoku Tsunami of 2011 [35] for which AI ≈ 5.25. Finally, in November 2011, we created a portable bore-soliton-splash with also a circa tenfold amplification in a miniature wave channel of approximately 0.7 m long, 0.1 m in width and 0.065 m high, where we used h 0 = 0.02 ± 0.001 m and h 1 ≈ 2h 0 to find H rw = 0.2 m. It again involved a solitarywave compound of a highest solitary wave followed by a second and third one of lower amplitudes. Perhaps surprisingly, nonlinear and inertial effects still dominate over friction and also over surface tension. A mini-splash video of this table-top dissemination is found online [50].

Mathematical and Numerical Modelling of Soliton Splashes
The bore-soliton-splash involves a series of mathematical and fluid-mechanical ingredients: dispersion, nonlinearity, a turbulent spilling breaker, and collapsing splash. Assuming incompressible fluid flow with a free surface, dispersion in a solitary wave is balanced by nonlinearity due to advection, while the hydraulic bore or spilling breaker highlights that this balance is temporarily and locally broken till turbulent dissipation reduces wave amplitude sufficiently to restore that balance, as we saw in Fig. 3c, d. When the flow is in balance, the soliton compound and splash can be modelled with a single-valued free surface in a singly connected domain till the apex of the splash is reached. Both spilling breaker and collapse of the splash are seen to involve multiply connected domains with bubbles and droplets. We will start our modelling of the bore-soliton-splash cases for a smooth singlevalued free-surface and using potential-flow equations and approximations thereof. Approximations used include a Benney-Luke pair/system of equations. Alternatively, one can explore the single, unidirectional KP equation in two horizontal spatial dimensions. These approximations have the advantage that dispersion is anomalously high which prohibits wave breaking and is, therefore, robust with the disadvantage being that outcomes during wave breaking will be less realistic, as follows. Numerical solutions are required to solve potential-flow and Benney-Luke equations in the actual wave channel, while exact solutions are available for the KP equation in an idealised domain for idealised settings. We will use variational principles and asymptotic theory to enhance numerical stability and robustness: our (novel) numerical techniques are direct, compatible space-time discretisations of relevant variational principles.
For potential-flow water waves, consider a free surface at z = h(x, y, t) = H 0 + η(x, y, t) over a flat bottom at z = 0 with vertical coordinate z, and horizontal coordinates x and y as well as time t. Acceleration of gravity g acts in the negative z-direction. Water at rest sits at z = H 0 and η = η(x, y, t) is the deviation from this rest level H 0 . Three-dimensional velocity u is approximated using a velocity potential φ = φ(x, y, z, t) as u = ∇φ with the gradient ∇ = (∂ x , ∂ y , ∂ z ) T . The horizontal part of the domain Ω h is defined by a main channel of width L x , with x ∈ [0, L x ] and length l y (x), with y ∈ [R(t), l y (x)]. The contraction is defined by y = l y (x) for y ∈ [L y − L c , L y ], where and the length of the contraction measured in the y-direction is L c . The piston wavemaker R(t) will be used later, and for the solid wall at y = 0 considered hitherto, we take R(t) = 0. Our derivation starts from Luke's [37] variational principle: modified to include a potential η R = η R (x, t) modelling the removal of the sluice gate (cf. [4]) and with the horizontal domain extent Ω h of the wave channel. Collecting all variables into the vector of unknowns variations of (3) are defined as follows: For R(t) = 0, the potential-flow water-wave equations resulting from (3) read: withn the unit normal vector and ∂Ω h the boundaries of the horizontal domain, i.e., any vertical walls-cf. the derivation in [37] with a minor change involving the sluice gate. To obtain the Benney-Luke system and a compatible/geometric finite-element discretisation on space and time, one can directly transform, discretise, and simplify (3b); for which details, we refer to [4,21].

Simulations
The Benney-Luke equations have been used to simulate cases 8 and 9. The simulation of case 9 is seen to be surprisingly good from the visual comparison between the photographic images (cf. Fig. 11.5 in [4]) and snapshots of the simulation (see Fig. 4). The simulation of case 9 shown is an improvement of the simulation in [4], because we have used better meshing with a symmetric mesh before the contraction, as seen in Fig. 10. The simulation of case 8 does not compare well with the photographic images and video material. The primary reason is that simulations of cases 8 and 9 with the Benney-Luke equations are fairly similar, as the comparison between these cases along the centreline of the wave tank in Fig. 5 reveal. Consequently, the reduction of the wave amplitude and, hence, wave speed of the first soliton due to wave breaking in case 8 does not occur. The observed resonance between the first and second solitary waves in which the second wave exactly falls within the trough drawn by the reflection is, therefore, absent in the simulation of case 8. The wave dispersion in the Benney-Luke equations is too strong. Either a potential-flow type model with parametrised wave breaking is required or a model with the dynamics of a water-air mixture, including the localised wave breaking inherent in such a two-phase model.

Novel Rogue-Wave-Energy Device
We have created, designed, and tested a novel wave-energy device inspired by the bore-soliton-splash event. Our rogue-wave-energy device involves wave-activated buoy motion. It has elements, i.e., wave-focussing due to a tapered channel and a power take-off (PTO) mechanism for the electro-magnetic-induction generator, of three existing wave-energy devices [15]:  Meanwhile, air compresses and decompresses by the rising and sinking wave leading to rapid air flow through the blow hole in which a wells' turbine is situated and generates electrical power when air flows in either direction.
Our device consists of a contracting channel with a wave buoy constrained to move in only one dimension, either in the vertical by sliding along a guiding mast or along a slightly curved arc pivoting around a horizontal axel at the contraction entrance. Attached to the buoy is either another vertical mast or a curved mast, to which magnets are attached that can move through a series of coils when the buoy is heaving due to the wave motion. An artistic rendering of the second version of the wave-energy device is given in Fig. 6. The latter magnet-and-coil system comprises a magnetic-induction motor, cf., the one in the Faraday shaking light shown in Fig. 7b. Relative to the version with the two vertical masts, one moving and one fixed, it has the advantage that the buoy can be taken out of action in storms and that a rotating axel is mechanically more robust than mast-guide ball bearings. Our device is intended to be part of a breakwater or dock, since waves will be absorbed.
In 2013, the "Berkeley wedge" wave-energy device was patented [38]. It is a waveabsorbing wedge moving on rails against a vertical wall attached to an induction motor. It is similar to our wave-energy device, but has no wave amplitude enhancing contraction like in the TapChan and OWC devices. During storms, the Berkeley wedge is sunk off into the water to protect it from damage. The Berkeley wedge operates in essence like a wavemaker reversed in time, following the principle that a good waveenergy device can also be good wavemaker, when time is reversed or energy is put in rather than being generated.
Before we advanced to any mathematical modelling, in the summer of 2013, we built and tested a proof-of-principle of our device to assess its viability. Our experimental set-up consisted of a straightforward wave tank with a hand-driven wavemaker, Fig. 7a, a shaped foam buoy with a mast topped by two magnets moving through a tube with four fixed coils, Fig. 7c, the latter parts coming out of two deconstructed shake or Faraday flashlights [22,25,39], Fig. 7b, with the buoy constrained to glide along a fixed, guiding mast. The shape of the buoy is close to a simplex with a slightly rounded Fig. 7 Overview of the working proof-of-principle of our new wave-energy device, here powering one LED: a wave tank with wavemaker, powered by OB, and contraction; b two Faraday shaking lights, one entire, and one deconstructed with the magnets put onto the mast and the coils wrapped around a plastic tube; c the tube guiding the magnet with its surrounding coils and wires leading to the LED; and, d the unit of contraction, guiding mast, and buoy-mast unit, at rest and slanted bottom face and a flat-top face, Fig. 7d. To the induction motor, either one LED was connected to demonstrate the power output or an Arduinoscope (a handmade oscilloscope using Arduino technology) to measure power output. Photographs of the set-up are given in Figs. 7 and 8. We powered one LED, but, given the AC-power generated, two (sets of) LEDs will be used in the mathematical model derived in the next section (see also [29,30]). Two (sets of) LEDs, circuited in parallel yet operating for currents in opposite directions, harness twice the wave energy into light. 3

Wave-to-Wire Monolithic Wave-Energy Model
A comprehensive mathematical model of the new wave-energy device will be developed next within a domain constructed to reproduce an existing small-scale wave tank at the University of Leeds, which is a larger tank than the one used for the proof-ofprinciple. Both the wave-to-wire formulation of this wave-energy device is novel as well as the compatible and robust discretisation of the linearised model. Such a discretisation is important, because it guarantees numerical stability and the correct two-way feedback between pairs of the three subsystems. The numerical wave tank has a piston wavemaker on its left side, and consists of a channel with a flat bottom at z = 0 that ends in a V-shaped contraction at the right end of the channel, cf. Fig. 9, as described in Sect. 3 and Eq. (2). A wave-energy buoy, here constrained to move only in the vertical, resides in the corner of the contraction. The shape of the buoy is described next.
The buoy is a simplex with a flat-top triangular face and has a slanted front face converging into one point at the bottom. The two remaining faces of the simplex align with the vertical walls of the V-shaped contraction. The slanted face is tilted to the vertical, as shown in the cross section at x = L x /2 in Fig. 9a). The buoy motion can be described by its position Z = Z (t) (note that given the constrained motion, this does not need to be the centre of gravity) and corresponding velocity W = W (t) = dZ /dt. Given this simplex geometry and assuming its flat-top face stays dry, the wetted buoy's height above the flat bottom is located at Here, α is the angle between the hull bottom and the horizontal, L y is the maximum length along x = L x /2 of the domain with the wavemaker at its rest position y = 0, and H k > 0 is the vertical distance between the position Z (t) and the keel of the buoy (see Fig. 9a). The keel of this V-shaped buoy, therefore, lies at (x, y, z) = Fig. 9 is defined as the point where the water meets the buoy. The overall waterline of the buoy is denoted by y = y b (x, t) and parameterised by x. At this waterline y b (x, t), for every x within the contraction, which x-interval varies over time, the height h = h(x, y, t) of the free surface of the fluid equals the buoy's surface , with the appropriate limits indicated. Taking the variation of this expression, one finds that which implies that y b is not an independent variable in the problem. The rest position of the buoy is directly determined geometrically via Archimedes' principle given the mass M and rest positionZ of the buoy, as follows. The angle θ c , defined by tan θ c = 2L c /L x , with θ c = 68.26 • presently, is the angle between the opening of the contraction and the line across the contraction (as shown in the mag-  Fig. 9). A summary of the dimensions of the existing wave tank as well as other relevant physical parameters used in the numerical calculations can be found in Table 3. At rest, y b = L b and given that tan such that the submerged part of the buoy is a smaller simplex, isomorphic to the entire buoy simplex, defined by the following four points: The volume V b of the submerged part of the buoy is then the displaced water mass divided by the density of water, i.e., given that tan Consequently, for this three-dimensional tetrahedral buoy, the rest position of the buoy's centre of mass is thus found to be: Attached to the buoy is a vertical mast with two magnets on top, moving through a set of coils, whose constrained and vertical movement by induction comprises the actuator. Magnet and mast are all included in the overall weight M of the buoy. The current I = I (t) is the derivative of the electrical charge Q = Q(t), i.e., I =Q. The coils have an overall inductance of L i . Rather than using the current I (t), we use the conjugate momentum P Q = P Q (t) as primary variable, defined in Appendix A by γ = 2πa 2 μN /L, magnetic dipole momentum μ of the magnet, a the radius of the coils, N the number of coil windings per metre, and L the length of the coils as well as a function G(Z ) defined in (44f). This function G(Z ) depends on the length of the mast H m , the length L m , and radius of the cylindrical magnet and the placement of the coils at z A novel model of the magnetic-induction actuator is developed in Appendix A from first principles, using the Maxwell's equations in a thin-wire approximation and given the cylindrical symmetry of the induction coils.
Hence, we can now formulate a comprehensive, or monolithic, variational principle of the entire, coupled water-wave problem, wave-activated buoy motion, and magneticinduction actuator, as follows: with φ s = φ(x, y, h(x, y, t), t) the velocity potential evaluated at the free surface, ∇ the three-dimensional gradient, and H the Hamiltonian. Herein, a scaled pressure p = p(x, y, z, t) acts as Lagrange multiplier to impose the incompressibility constraint D − 1 = 0 of a scaled density D = D(x, y, z, t), such that the density ρ(x, y, z, t) = ρ 0 D(x, y, z, t), cf. [9]. Moreover, the single-valued free water surface at z = h(x, y, t) with water depth h = h(x, y, t) is constrained [underlined terms in (12)] to be the dynamic shape of the wave-buoy using the Lagrange multiplier λ(x, y, t) over the wetted part y > y b (x, t) of the wave-buoy hull-see also [30]. We have used the Heaviside function Θ(y − y b ), zero for y < y b and unity for y ≥ y b , to single out this wetted part of the hull. The key reason to include the scaled density D and impose incompressibility condition D − 1 = 0 weakly is that variational principle (12) mathematically yields the boundary condition on λ at the waterline, as a consistency component of the entire formulation. This condition follows from the system of equations emerging next. In contrast, when one imposes the constraint D = 1 strongly, the (interface) condition for λ at the waterline needs to be imposed a priori.
As before, we collect all variables into the vector of unknowns with variations of (12) defined as in (4). Most variations of (12) emerge in a straightforward manner; only the variations involving the terms D∂ t δφ, D∇φ · ∇(δφ), and a comprehensive term involving the variations of the free surface δh in the upper integration limit are more complicated and require integration by parts in time and the use of Gauss' law-see [9,18] for more details. Given these hints, while leaving further derivation details to the reader, variation of (12) yields the following, fully nonlinear, equations of motion: δh with φ R = φ(R(t), y, z, t) the velocity potential evaluated at the wavemaker. Evaluation of (13a) at the free surface z = h and subtraction of (13f) yields that such that y = y b (x, t) at the waterline on the wave-buoy, from which it necessarily follows that λ (x, y b (x, t), t) = 0.
A straightforward way to model the energy harvesting is the use of two (sets of) LED diodes in parallel but positioned in opposite directions, such that only one (set of) LED(s) is active at one time. Using the Shockley equation, as model for the LEDvoltage-current relationship for current flow in either direction, yields that the voltage V s (I ) = −sign(I )n q V T ln (|I |/I sat + 1) is really a function of I = I (t); the sign-function and absolute value |I | used ensure operation and damping for currents in either direction. Since such an LED model leads to damping, it is less common to include it a priori in the variational principle and damping/loading is added a posteriori to the model by substitution of (15) for V s (t).
Parameters in the Shockley model for LEDs include the saturation current I sat , the quality factor n q , and the thermal voltage V T . In addition, we added resistance terms −(R c + R i )I , with a resistance R i of the wiring to the LEDs as well as a resistance R c of the coils, in the equation for P Q , to model losses for the circuit and LED diodes combined. Instead ofṖ Q = 0 in (13), one then obtainṡ The total electrical power output P is then the time integral P = T 0 I (t)V s (t) dt, with V s (t) the voltage across the LEDs.
In the rest state, we saw that the straight waterline lies at y = L b with rest depth H (x, y) = H 0 for y < L b , while for y ≥ L b , a rest depth H (x, y) as well as (rest-state) Lagrange multiplier Λ(x, y) are defined by: with the rest positionZ of the buoy. Hence To linearise the equations of motion (13), we consider the following decomposition of variables into rest-state and (small-amplitude) perturbations: φ(x, y, z, t) =φ(x, y, z, t), D(x, y, z, t) = 1 +D(x, y, z, t), the latter relation taken, such that the rest currentĪ = 0, since P Q = L i I − K (Z ). Upon linearising, the moving domain becomes a fixed domain y ∈ [0, l y (x)] with a fixed waterline at y = L b . After a Taylor expansion, the waterline condition forλ at y = L b becomes: where we used that Λ(x, L b ) = 0 by definition and omitted quadratic and higher order terms in the perturbation variables. Expansion of the waterline condition (7) leads to an explicit expression forỹ b (x, t) in terms ofZ (t) and the free surface at the linearised waterline η(x, L − b , t); to wit By combining (20) and (21), we obtain the desired boundary condition at the linearised waterline, that isλ It means thatλ(x, L b , t) is not an independent variable at y = L b . We highlight that the derivation of this subtle (weak) condition (22) onλ has only been possible, because we deferred the strong imposition of the incompressibility condition D = 1 orD = 0.
The final simplifications are that we consider the system in both the linear and shallow-water limits, yielding that the velocity potentialφ =φ(x, y, t) is a function of only the horizontal coordinates and time, now with ∇ = (∂ x , ∂ y ) T . The equations of motion and induction in these linear, shallow-water limits become: where we have added the effective circuit and coil resistances R i and R c and linearised (around I = 0 and usingZ ) model of the LED light, as well as a consistency equation, with (22) andn · ∇φ = 0 at the fixed, vertical walls. Note that from (23f),Ĩ = (P Q + γ G(Z )Z )/L i . When we rework (23g) and takeİ to be negligible, then the linearisation is seen to lead to a damping term proportional toŻ =W in the momentum equation, cf. [11,12]. In addition, the linearisation of the LED model (16) is seen to lead to an effective resistance n q V T /I sat , with linearised LED voltageṼ s (Ĩ ) = n q V T /I satĨ . The (linearised) electrical power outputP(t) is calculated as a function of time and is given byP = (Ṽ s +Ṽ r )Ĩ , withṼ r = (R c + R i )Ĩ the lost voltage due to circuit resistance.
A consistency equation arises by taking the time derivative of the primary constraint η −Z = 0, which leads to a secondary constraint, given by upon using two of the equations of motion to eliminate time derivatives. Subsequently taking the time derivative of this secondary constraint (24), while using two other equations of motion to eliminate the emerging time derivatives, yields the elliptic equation (23h) forλ.

Time Discretisation of the Linearised System
The time discretisation of the eight main equations in (23), i.e., in a count excluding boundary conditions, needs to be such that a time discretisation of the first seven equations is equivalent to, and consistent with, a time discretisation of the last seven equations, given that there are seven unknowns. The chosen time discretisation is a symplectic Euler one [19,34], given the conjugate pairs of variables {φ, η}, {W ,Z }, {P Q ,Q}. Hence, the time-discrete system with discrete time levels t n and t n+1 = t n + Δt reads: , and C 2 = (R c + R i + n q V T /I sat )/L i . By subtracting theZ -equation from the η-equation, we obtain: Given that η n −Z n = 0 andW n + ∇ · (H ∇φ n ) = 0, both imposed at time level n = 0 initially, ensuring that η n+1 −Z n+1 = 0 implies that we have to show that Using the equation forW n+1 and operating ∇ · (H (x)(·)) on Eq. (25b) forφ n+1 , the consistency condition is seen to be (25h). Hence, (25) is consistent.

Space Finite-Element Discretisation of the Linearised System
The model (23) is discretised in a few steps. The first step is to multiply the field equations in (23) by C 0 -test functions, and then integrate over space and by parts. The second step is to expand the fields using (special) C 0 -continuous and compact finite-element basis functions. We will use standard linear and compact Galerkin basis and test functions, which are unity at their home node and zero at neighbouring nodes of the elements connected to the home node. The result is a space-discrete system which will be revealed to be only consistent for certain special choices of the function spaces and expansions. Vice versa, we can first discretise time in a consistent manner, such that again the equations remain consistent as we showed already. Finally, by either discretising the space-discrete system properly in time or the proper time-discrete system in space, we obtain an internally consistent overall space-time discretisation fit for numerical implementation. The mesh is assembled using uniform quadrilateral elements up to the entrance of the contraction, with N x elements in the x-direction and N y elements in the y-direction, i.e., in the uniform section of the domain, the mesh comprises N x N y elements and (N x + 1)(N y + 1) nodes. In the contraction, the chosen mesh is still formed by quadrilateral elements, but now nodes are only aligned in every other line. A sample mesh is shown in Fig. 10. The tessellation of the contraction region increases the total number of elements by N x (N x + 1)/2 and the number of nodes by N x (N x + 2)/2. It is thus clear that the way the mesh is constructed provides a restriction in the choice of N x , which needs to be even. The total number of elements is N el = N x N y + N x (N x +1)/2 and the total number of nodes is N n = (N x + 1)(N y + 1) + N x (N x + 2)/2. The reststate waterline is such that it is aligned with one of the nodal lines parallel to the x-direction.
The nodes in the entire mesh are denoted by k, l = 1, 2, . . . , N n and the N n − N p +1 nodes under the buoy are in the ordered case denoted byk,l = N p , . . . , N n . The latter include the nodes on the waterline, which in the current linearised case lies on the line y = L b ; the N b nodes on the waterline are a subset thereof, in the ordered case denoted byb = N p , . . . , N p + N b − 1. When the N b nodes on the waterline are excluded in the latter, we use indexk. We multiply both wave equations (23b)-(23c) by the C 0 -test function ϕ k (x, y) and the elliptic equation (23h) as well as the constraint (23a) by ϕk(x, y) to obtain the weak forms after integration by parts, upon using the Neumann/Dirichlet conditions at y = 0, x = 0, L x for y ∈ [0, L y − L c ] and y = l y (x) for x ∈ [0, L x ] and y ∈ [L y − L c , L y ]. Also, λ h (x, y, t) = λk(t)ϕk(x, y)+ λb(t)φb(x, y), withφb being the part of the test function ϕb ≥ L b under the buoy and λb = g(ηb −Z ). Unfortunately, the consistency required cannot be shown for this choice of test functionφb at the waterline. We, therefore, made some adjustments: we take λ h (x, y, t) = λk(t)ϕk(x, y) to be the normal test function spanning across the waterline withk,l = N p , . . . , N n and thus smooth out the Heaviside function to allow inclusion of the full basis function ϕb at the waterline y = L b . Hence, the Heaviside function Θ(y − L b ) is removed. This changes de facto only the vector and matrix definitions ofSkl ,Qb and N kb below. The corresponding finite-element discretisation then becomes Skl + CQkQl λl = −gSk l η l − CQkQbλb −Skbλb with several mass and "Laplace" matrices defined by: Nonzero contributions only exist for the tilded matrices and vectors for certain index ranges.
The key consistency check is to ensure that the first seven equations in (28) are consistent with the last seven equations in (28). Consider the first seven equations. Take the time derivative of the primary constraint and eliminate the time derivatives using two of the other seven equations, to obtain the secondary constraint: Now, take the time derivative of this secondary constraint above and again eliminate the time derivatives using two different equations of these seven equations, to obtain the consistency equation

This consistency equation matches the last equation (28h) if and only if the following relations hold:Skl
These relations have been verified to hold up to machine precision. To date, we have not been able to verify these relations analytically. Finally, we find the consistent space-time discretisation by logically combining the time-discrete and space-discrete approaches derived in (25) and (28).

Numerical Results
We have set up a numerical code which simulates the full system as it evolves in time, including the generation/propagation of waves, their impact on the wave-energy buoy, the response of the buoy, and the power output. The numerical results presented next have been obtained using a mesh resolution of N x = 10 and N y = 50, i.e., the total number of elements in the calculations is N el = 555 and the total number of nodes is N n = 621. The time step used is Δt = 0.0028 s. At the start of the simulation the system is at rest and the water depth in the main wave tank is H 0 = 0.1 m. For t > 0, waves are generated from the left wall of the tank by a piston wavemaker that follows a periodic motion in time according to R(t) = A ω (1 − cos(ωt)), with amplitude A = 0.0653 m and frequency ω = 6π L y √ g H 0 = 9.3348 s −1 (which corresponds to a physical frequency of ω/2π = 1.4857 Hz). Therefore, on the left wall, ∂ yφ =Ṙ(t) = A sin(ωt).
The total energy of the system E(t) is computed numerically as a function of time. The respective energies of the water E w (kinetic and potential), the buoy E b , and the electro-magnetic system E i are defined by (the continuum forms are given below): Applying appropriate manipulations of equations in (23), by integrating in space the respective equations valid within the fluid and adding all of the resulting expressions, results in an equation for the time derivative of the total energy in time. This is not expected to be zero as in conserved systems but negative, due to the added damping in the system. In particular, we find: with E = E w + E b + E i the total energy andP =P l +P g the loss, separated into the resistive dissipationP l as well as the power gainedP g via the LEDs, defined by: The results of a simulation with the parameters described earlier can be seen in Figs. 11 and 12. The response of the buoy is shown in the top panel of Fig. 11, while the electrical current is shown in the middle panel and the power generated in the  Fig. 12.
In what follows, we consider a smoothed-out motion of the piston wavemaker which only attains its maximum amplitude A (as described in the first paragraph of this section) after two periods of oscillation. The wavemaker motion comprises two signals R 1 (t) and R 2 (t), which are out-of-phase initially, become in-phase after the first two oscillations, and then return back to be out-of-phase after N oscillations. The wavemaker hence stops operating after N oscillations at t = T w = 2π N /ω and the simulation continues until T = 2T w . Consequently N wave packets are sent into the wave tank. The motion of the wavemaker (and the respective signals R 1 , R 2 ) can be seen in Fig. 13, and the resulting wave deviation from the still water level H 0 = 0.1 m at the left wall is shown in Fig. 14.
Before performing any further simulations, we tested the accuracy of the numerical scheme both in time and space. The convergence of the Symplectic Euler method for the time-integration results is investigated by plotting the energy signal E(t) for two time steps Δt and Δt/2 (Fig. 15). The bounded oscillations in the energy signal (shown in Fig. 15 only after the wavemaker stops adding energy in the system, i.e., for t ≥ T w ), are seen to be approximately reduced by 41% when computed with the half time step indicative of an accuracy of ∼ (Δt) 3/4 . We note that the energy calculation using half the spatial resolution gives an accuracy of √ Δt, which suggests that the desired convergence of Δt may be achieved, albeit slowly, for a sufficiently resolved mesh. The spatial convergence rate evaluated using the formula, cf. [47]:  Fig. 13 The motion of the piston wavemaker in time, defined by the composition of two signals. The second signal is constructed to have a time-dependent phase difference from the first signal, so that the wavemaker reaches its maximum amplitude smoothly, and returns to the stationary position in a similar manner. The wavemaker operates for 0 ≤ t ≤ T w , with T w = 6.7309 s, and the simulation runs until T = 2T w = 13.4618 s where the norm is taken to be either the L 1 , L 2 or L ∞ norm. The subscripts 1, 2, and 3 correspond to the different mesh sizes in Table 4, obtained by halving the mesh spacing uniformly in each direction, or equivalently by doubling the number of elements. The convergence of the spatial discretisation, while expected to be of order 2, is seen to be circa 1.7.   Table 4 Convergence rates n using three different norms (L 1 , L 2 , L ∞ ) evaluated using the value of the velocity potentialφ at the final time of the simulation, i.e., at t = T . The time step used is such that the CFL condition is satisfied in similar fashion for each mesh resolution We next investigate the influence of various system parameters such as the induction coils, buoy mass, and wavemaker frequency on the amount of power generated (and used to illuminate the LEDs) or lost due to circuit resistance. Figure 16 Panel (a) shows that the total power generated increases linearly with C s , while the power lost decreases and is seen to asymptote for larger values of C s . We find that for larger values of C s , the currentĨ remains unaffected as a consequence of the fact that the wave state in the tank and the buoy response do not change with C s . The linear dependence of [P g ] on C s is due to the linearisation of the model. Furthermore, the increased power generation seen in panel (b) for smaller values of N is anticipated as it corresponds to lower resistance and inductance in the system (recall the definitions of R c and L i -see Table 5-which are both proportional to parameter N ). The curves in panel (b) are both seen to be proportional to 1/N 2 ; hence, taking even lower values or N (or equivalently a shorter coil) yields a much higher production of energy. The variation of the power with the mass of the buoy is seen to be minimal (panel c), but follows a monotonic increase. Higher values of M were not considered here to avoid having a very low centre of gravity relevant to the still water level. Different buoy, contraction, and wave tank geometries/sizes (and hence buoy masses) will be considered in future work. Finally, it is observed that an optimal power generation can be achieved at certain wavemaker frequency, in this particular case at ω = 10.89 s −1 (see also [3] for related results for the dependence on wavemaker amplitude and frequency for a slightly different set-up). This critical value appears to be the resonant frequency of the system considered. The energy of the system around this optimal frequency is illustrated in the top panel of Fig. 17, shown for times T w ≤ t ≤ T with T w = 6.7309 s and T = 2T w . The respective energies of the subsystems (water, buoy, and electro-magnetic system) are also shown in the bottom panels of Fig. 17. Clearly, the total energy shown in the top panel remains approximately constant (the small oscillations observed remain bounded in time), while the energy lost from the water around t = 7 s and t = 11.5 s is converted into kinetic energy of the buoy and electrical energy in the inductor, even though the latter is much smaller in this case. One way to improve the efficiency of the wave-energy converter is by adding more LEDs in serial and/or reducing the length of the induction coils (cf., Fig. 16a, b).

Summary and Discussion
In summary, we have reported in detail on the creation of the bore-soliton-splash, summarised modelling of this hydrodynamic splash, and showed how it inspired a novel wave-energy device. We will next provide further context of our work.

Relation to rogue waves at sea
The bore-soliton-splash is a nonlinear wave-resonance phenomenon in which a series of travelling solitons reflect in a V-shaped contraction leading to a tenfold resonant amplification of the initial main wave height. Once created in 2010, the phenomenon caught attention of the rogue-wave community. Rogue waves are extreme and rare waves, generally but not exclusively sea waves, at least twice as high as the wave height of the ambient sea. While the original bore-soliton-splash was engineered, it relates to several rogue-wave phenomena at sea, far from and near the coastline. Rogue waves can have several causes and emerge in different situations: rogue-wave emergence in crossing seas, either due to seas with two main directions, e.g., high seas generated by two hurricanes, or one hurricane changing direction. More rare are seas with waves and swell from three different main directions. Our V-shaped contraction walls can, therefore, be re-interpreted as virtual walls concerning two (virtual) waves travelling under two angles ±ϕ from the main wave's direction with two (virtual) planes of no-normal flow leading to a converging point. One difference is that the virtual case supports wave propagation with one splash at one space-time point, cf. [20], while our engineered set-up with solid walls necessarily leads to reflections.
Further modelling of the bore-soliton-splash and related rogue waves We reported (bore-)soliton-splashes including one with smooth solitary waves in which nonlinearity and dispersion are balanced without any wave breaking. The smooth soliton-splash of case 9 was successfully simulated using a compatible, geometric finite-element discretisation of a Benney-Luke model, a bidirectional simplification of the classic potential-flow model for water waves. Case 8 for the maximum bore-soliton-splash was not simulated correctly by the Benney-Luke model. Due to the lack of wave breaking, the case 8 simulation deviated significantly from reality in that the observed resonant interaction was absent. On one hand, it is possible to further explore simplified modelling of rogue waves in crossing seas using the Kadomtsev-Petviashvili (KP) equation, cf. [1,20,21,28,32]. On the other hand, more advanced modelling of the bore-soliton-splash will either require a full potential-flow model with localised and parameterised wave breaking or the use of models with actual, localised wave breaking while maintaining good dispersion properties. Single-phase or two-phase , and electro-magnetic system (bottom right). The overbar notation denotes the difference between total energy and the energy when the wavemaker no longer moves, i.e., at T w = 5.77 s. Here, the "optimal" wavemaker frequency ω = 10.89 s −1 is used mixture-theory models, including ones with a Van-der-Waals-type equation of state, may also be good candidates [23,46]. Alternatively, we have explored use of Smoothed Particle Hydrodynamics (SPH) 4 ; such a numerical method [24] can handle wave breaking and jet collapse into bubbles and droplets, and can be used to simulate the bore-soliton-splash of case 8. SPH lies at the other end of the spectrum of numerical techniques, compared to the geometric numerical techniques used by us. To date, SPH does, however, turn out to be too dissipative. It requires too many degrees of freedom and thus too much computational power to avoid detrimental numerical wave-amplitude dissipation. Consequently, dispersive wave propagation over longer distances remains relatively poor in SPH. In our attempts, to date, the splash amplitude simulated by SPH was too low, even though we had replaced wave propagation in the channel prior to the contraction by several laterally periodic slices, copies of one another, to significantly reduce computational resources, before waves in these slices were fed into the contraction.
Optimisation of the wave-energy device Our splash inspired the creation of a novel wave-energy device and we showed a working, experimental proof-of-principle model, Fig. 18 The steel-soliton-splash is an artistic rendering in stainless steel of two snapshots of the bore-soliton-splash of case 8 [48,49] but also developed and derived a new and fully nonlinear mathematical model of the combined water-wave dynamics, the wave-activated buoy, and the magnetic-induction power generator. Essential ingredients of this comprehensive model have been captured in one variational principle to which we a posteriori added dissipative effects of the electrical circuit, coils of the actuator, and LEDs used as the loads. The overall model was subsequently linearised and discretised using a finite-element method in space and time. This (linear) algebraic model was made fully compatible with the variational structure in the conservative and continuum limits. Its compatible, novel, and nontrivial discretisation was augmented with the resistances of the electrical circuit and coils of the induction motor as well as the LED loads. Preliminary simulations of the linear model showed promising results including (suboptimal) convergence and energy transfer between the three components. Finally, we investigated the resonant behaviour of the system as function of wave-frequency and load for a long wave-packet of harmonic waves. Nonlinear modelling, optimisation, and control of the wave-energy device require further exploration, and both the geometry of contraction, mass, and wave-buoy shape could be optimised for a given wave climate. We also aim to explore feedback control as function of contraction geometry, the number of coils of the induction motor, and the total load. In addition, higher order and more accurate spatial and temporal discretisation schemes require exploring.
Steel-soliton-splash artwork We finish on an artistic note. Our bore-soliton-splash inspired an artwork, the steel-soliton-splash [49], created by WZ. Snapshots of the video of case 8 were first outlined as silhouettes, two of which formed the basis for a three-dimensional artwork, scaled down by about a factor of three, and welded in stainless steel, see Fig. 18. In 2013, we donated the artwork to the Isaac Newton Institute of Mathematical Sciences in Cambridge, UK.
In spherical coordinates (r , ϕ, θ) and with the centre of the coil placed at the origin, the approximate expressions [27,36] of the two nonzero components of the magnetic flux, for distances r A m larger than the radius A m of the magnet, are given by the spherical-coordinate components B r = 2μ cos θ/r 3 , B θ = μ sin θ/r 3 of the magnetic flux B, or rewritten in cylindrical coordinates (ρ, φ, z) for the radial cylindrical-coordinate component of the magnetic flux as: The magnetic field B is expressed in Teslas, T, the permeability μ 0 of free space in Henry's per metre, H/m = Tm 2 /A, the dipole moment m in Ampere-square-metre, Am 2 , and the radius r in metres. The above holds for an infinitesimally short magnet and far away from the magnet, relative to the distance at which the magnetic flux is measured and relative to the length of the coils, cf. [11] [their equations (1) and (2)]. The magnetic field induced by a magnet of finite length L m is considered next and reduces in the far field to the above expression. The magnetic flux B is related to the magnetic field intensity H and magnetisation density M as follows (page 395 of [36]): In the absence of a current density J, the magnetic field density is irrotational, such that H = −∇Ψ outside the magnet, with magnetic scalar potential Ψ . Since ∇ · B = 0, we find from (37) that ∇ · H = −∇ · M and, therefore, that with magnetic charge density ρ m . The integral solution of (38) reads: For a magnet uniformly magnetised in the z-direction, we have a magnetisation density M = M 0ẑ or, alternatively, the magnetic density ρ m = 0 throughout the magnet except at the end surfaces where ρ m = ± μ 0 M 0 δ(z ∓ L m /2). The magnetic scalar potential Ψ is given by the triple integral (39) over the length L m of the magnet of radius A m < a with a the radius of the inductor's coils, cf. [36] (pp. 395) and [26]: in which, based on symmetry, it suffices to takeφ = 0 with x = ρ cosφ etc., such that |r − r | 2 = (R cos φ − ρ) 2 + R 2 sin 2 φ + (z − z ) 2 . The expression of the relevant component of the magnetic field evaluated for a magnet of length L m and radius A m < a then becomes: Far away from the magnet, this radial component can be approximated using Taylor expansions, such that which equals (36) with (35) when we identify m = M 0 π A 2 m L m as the magnetic/magnet's dipole moment with π A 2 m L m the volume taken up by the magnet. In (41), the magnet is placed in z ∈ [−L m /2, L m /2], while, in our wave-energy device, the magnet resides between z ∈ [Z (t) Consequently, we need to adapt expression (41) as follows: Second step Placing the magnet in a moving range z ∈ [Z (t) + H m − L m /2, Z (t) + H m + L m /2] with length H m above the centre of mass Z (t) and for fixed coils with N turns and a radius ρ = a over a coil length L over a fixed range z ∈ [Z + (1 + α h )H m − L/2,Z + (1 + α h )H m + L/2] with 0 < α h < 1, the expression (43a) can be integrated in z over the coil to obtain the coil density [cf. [11] their (6)-(8)] where we have used that at the faces of the moving magnet  Fig. 19 Sketch of the magnetic-induction motor set-up: on the left the system when the magnet passes its rest position and on the right away from the rest position with the latter approximation holding in the far field, starting from (44d) and (43b). Note that only when α h = 0 and the water and the buoy-mast-magnet system are at rest, the magnet sits in the middle of the coils. A sketch of the various coordinate systems is given in Fig. 19. A graph of γ G(Z ) versus Z and its approximation shows that the two functions lie close together, cf. [3]. Ohm's law for a circuit subject to a magnetic flux, moving with a speedŻ in the vertical, reads: with current density J, electrical field E, conductivity σ > 0 for a conductor, and unit vectorẑ in the vertical. Integration of Ohm's law (45) around the coils at ρ = a yields (extending [8,11]): where dl = a dθ(− sin θ, cos θ, 0) T and B =B ρ (cos θ, sin θ, 0) T with θ ∈ [0, 2π ], inductance L i of the inductor and the voltage drop V c . Note that the resistance of the coils is defined by R c = N (2πa)/(σ π D 2 /4) and the current magnitude by J = I /(π D 2 /4) with D the cross-sectional diameter of the coil. Hence, the circuit equation becomes: where we have added the resistance V s (I ) (15) of the two (sets of) LEDs modelled using combined Shockley equations, placed in parallel, as well as the resistance R i of the remaining wires to and from these LEDs. It should be possible to add I R i and the Shockley voltage V s directly in Ohm's law, but we simply added the two terms heuristically. Faraday's induction law (chapters 7 and 8 [36]) used above follows from one of the Maxwell's equations, ∇ × E = −∂ t B, for one circuit as: with induction Φ, surface S bounded by the path of the line integral, and surface element da on the coils of the inductor. The magnetic induction inside a long solenoid is B = μ 0 (N /L)I (cf. (8.59) in [36]), such that Φ = μ 0 (N /L)I πa 2 (cf. (8.60) in [36]) and, hence, for N windings L i = K μ 0 (N 2 /L)πa 2 (cf. (8.61) in [36]). For a long solenoid, K = 1, and for a short solenoid, K < 1 with K = 0.53 for a/L = 1 ( Table 8-1 in [36]). Third step Faraday's expression for the magnetic force F, the force on charged particles in the coils as used in the momentum equations of magneto-hydrodynamics, is: in which we have replaced the magnitude of J integrated across a string of coil with area ΔS by I /ΔS multiplied times ΔSδ(ρ − a) to focus its averaged effect at ρ = a exclusively. Hence, the vertical component F = F ·ẑ of the force between the conducting coils and the magnet is: To facilitate further analysis, we momentarily simplify the hydrodynamic force on the buoy, used in the main text, by a mass-spring component with spring constant k s . Combining this simplification with the vertical momentum equation of the simplified buoy-mast system, we then arrive at the following coupled mechanical and magnetical system:Q L iİ = γ G(Z )Ż −(R c + R i )I − sign(I )n q V T ln |I | with constant γ ≡ 2πa 2 μN /L and underlined dissipative terms. When we ignore the self-induction term L iİ and the Shockley expression for the LEDs in (51), we note that (R c + R i )I = γ G(Z )Ż ; elimination of I then shows that the magnetic force in the vertical momentum equation (51d) acts as a (nonlinear) drag, proportional tȯ Z or W , cf. the linear analogue in [11]. In the absence of the underlined, linear, and nonlinear dissipative terms in (51), the system (51) should be conservative, which will be explored next. In this conservative limit, we first rewrite (51) as: The system (52) has a Lagrangian L mm , defined in and satisfying the following variational principle: A Legendre transform of (54) yields the following conjugate momenta: The Hamiltonian H (Z , W , Q, P Q ) is then the Legendre transform of L(Q,Q, Z ,Ż ) as follows: H (Z , W , Q, P Q ) = P QQ + M WŻ − L mm (Q,Q, Z ,Ż ) The variational principle (54) in terms of these new variables then becomes: yielding, cf. (52), the system of equations: The modification (12) to the variational principle (57) without the potential energy of the spring and, instead, including the hydrodynamics was used in the main text. A summary of the physical parameters introduced in this appendix and some representative values can be found in Table 5.