<p>A new tool is developed for simulating three-dimensional (3D) water-wave motion using the first fully variational 3D discretisation in space and time of Luke’s variational principle (VP), with an additional focus on the formation and analysis of extreme waves generated within in-house experimental wavetanks. The resulting “numerical wavetank” is able to emulate laboratory sea states in which complex wave-wave and shoaling interactions occur. After first transforming the time-dependent free surface and oscillatory wavemaker into a static rectilinear domain with fixed boundaries using a <InlineEquation ID="IEq1"> <EquationSource Format="TEX">\(\sigma \)</EquationSource> <EquationSource Format="MATHML"><math> <mi>σ</mi> </math></EquationSource> </InlineEquation>-coordinate transformation, a fully variational approach is used to derive a system of weak formulations that leads to a non-autonomous space-discrete Hamiltonian system to which robust (stable and mass-conserving) temporal integrators are applied. Specifically, time-discrete VPs with second-order Störmer–Verlet and modified-midpoint time-integration have been derived (and directly implemented through automation) in the finite-element environment <i>Firedrake</i> rather than their explicit weak forms, with spectrally accurate higher-order finite elements. Verification and validation of the new tool are demonstrated via a novel convergence analysis and comparisons of its numerical results with a new analytical solution of three-dimensional two-soliton interaction as well as data post-processed from wavetank experiments.</p>