Version 10.05:

• Optional inclusion of the method to account for the effect of ambient currents, with or without vertical shear, on the wave dynamics. The current field (e.g. tidal and wind-driven) can be taken from a circulation model or observations. This method is described in Rijnsdorp et al. (2024).
• Both first and second order Stokes theory, including transfer functions of sub- and super-harmonics, can be applied at the wavemaker. See for details the paper of Vasarmidis et al. (2024).
• The user can explicitly specify a bichromatic wave at the wavemaker.
• Merging or not merging of the individual output files in case of parallel runs is included as an option.

Version 10.01:

• the staggered Arakawa C-grid finite difference method with built-in mimetic properties¹ on Delaunay-Voronoi meshes - the covolume method - is employed for solving the nonlinear two-dimensional (2DH) and three-dimensional (3D) shallow water equations (SWEs) while second order spatial accuracy is obtained by means of high-resolution TVD schemes;
• a second order cell-centered finite volume discretization for 2DH/3D transport equations on unstructured triangular grid is developed that satisfies mass conservation and the maximum principle;
• a first order unstructured finite difference scheme using triangular cells for solving the vertical k-ε turbulence equations including wall functions is implemented;
• semi-implicit time integration, i.e. the θ-method, for SWEs with Newton-type iteration is the default in case of unstructured meshes;
• linear solvers (CG and BiCGSTAB) and ILU preconditioner are included for unstructured mesh matrices; the use of SPARSKIT (sparse matrix computations library) is now redundant;
• the following forces are added to the momentum balance: wind stress, baroclinic forcing and atmospheric pressure gradient.

Version 9.01:

• Optional inclusion of floating rigid bodies with six DoF motions and fluid-structure interaction (FSI) coupling. The associated method is described in Rijnsdorp et al. (2022). See also for details the User Manual. (Scroll down and search for command BODY DIM.)
• Time series of translation/rotation of rigid bodies and PTO power as output quantities are added.
• Supports internal wave generation for rectilinear non-uniform grid.

Version 8.01:

Version 7.01:

• Triangle
• Easymesh

Version 6.01:

• Internal wave generation is included as an option. The associated method is described in Vasarmidis et al. (2019). See also for details the User Manual. (Scroll down and search for command SOURCE.)
• Simulation of a full three-dimensional turbulent flow with anisotropic Reynolds stresses is included as an option. See for details the User Manual. (Scroll down and search for command VISCOSITY.)

Version 5.01:

• The simultaneous occurrence of free-surface flow and the pressurised flow is included as an option. The piezometric head is then calculated instead of the free surface elevation. This option accounts, for instance, for the presence of a non-moving floating body in the simulation. This enables to simulate the evolution of a complex wave field and its impact on a ship that is moored in coastal waters. The associated numerical approach that correctly describes global mass conservation in both the free surface and the pressurised region is extensively outlined in this paper.
• Output of hydrodynamic loads acting on a moored ship (forces and moments) is optionally included as well.
• Vegetation modelling is extended for the case of densely spaced, vertical cylinders (e.g. dense mangrove fields, porous brushwood groins), whereby the effect of (relatively low) porosity is optionally included.
• Wave damping due to vegetation schematised as horizontal cylinders is included as an option. The associated dissipation is based exclusively on the drag force.
• Time series of wave runup as output quantity is added.

Version 4.01:

• Two anti-creep methods are included as an option. Since SWASH makes use of the terrain-following coordinates it necessitates the transformation of the transport equation for salt, heat or sediment due to the curved z-planes. This includes some curvature terms, which may complicate the computation. However, when these curvature terms are neglected, this may lead to a false generation of vertical mixing. This effect, known as the artificial creeping, becomes evident when the bottom slope is relatively large in regions of strong stable stratification. Hence, in such as case, inclusion of the curvature terms, known as the anti-creepage terms, may reduce significantly the artificial creeping. The standard method is based on the local transformation. An alternative is the method of Stelling and Van Kester (1994), which computes the horizontal diffusion along strictly horizontal planes.
• A subgrid approach is included as an option to resolve the ambient flow and pressure field with a different resolution in the vertical. In this way, the solution for the non-hydrostatic pressure and vertical acceleration can be obtained on a relatively coarse vertical grid, while on a subgrid with a high vertical resolution the vertical structure of turbulent flow is resolved. For details, see this paper.
• Bottom roughness based on Colebrook-White method is included as an option.
• Inertia force due to vegetation is included as an option.
• Cohesive suspended sediment transport is included.
• Depth-limited wave growth by wind is included.
• A new output parameter is added: vorticity of the depth-averarged flow.

Version 3.14:

• The reduced pressure equation method is included as an option. Since most of the computational effort is devoted to inverting the Poisson pressure matrix, an effective way to minimize this effort is by reducing the dimension of the Poisson matrix. This leads to less computational cost and memory. In spite of this reduction it can still provide an accurate description of dispersive waves.
• An option to specify the bottom level in the centers of computational cells is added. To determine the total water depth at water level points, a bottom level in the cell center is required. As a rule, we take the average bottom level from the surrounding bottom corners to determine the total depth at cell center. However, in the vicinity of steep bottom slopes, use of the average bottom level may lead to an inaccurate representation of the flow. In that case it is more appropriate to take the bottom level at cell center being equal to the minimum of the surrounding bottom levels. Other alternatives (a maximum instead of a minimum and shift from a corner to cell center) are also included.
• Mass exchange of suspended sediment between the bed and the flow is included. Only suspended load for noncohesive sediment (e.g. sand) is modelled. The pickup function approach is employed. This function models the amount of sediment eroded from the bed into the fluid, i.e. upward sediment flux, and is related to the bottom shear stress. The downward sediment flux depends on the settling velocity only.
• The linearized Riemann invariants to be used as boundary conditions are included. This option is applicable for a subcritical flow in relative deep waters. Examples are tidal flows in a continental shelf or in a harbour.
• Output of the physical vertical velocity in case of a hydrostatic computation is included.

Version 2.00:

• Wave damping induced by aquatic vegetation (mangroves, salt marshes, etc.) at variable depths is included as an option. The calculation of this type of dissipation is specified by the drag coefficient, the stem diameter of plant (schematised as a vertical cylinder), the number of plants per square meter and the vegetation height. In addition to the vertical variation, the possibility of horizontal variation of the vegetation characteristics is included as well. This inclusion enables the vegetation in a given region to be varied so as to reflect real density variations in the field.
• A number of transport equations is added. Specifically, the user has the option to include the equations for conservative transport of salinity, heat and/or suspended sediment. These transport equations are coupled with the momentum equations through the baroclinic forcing term, whereas the equation of state is employed that relates density to salinity, temperature and sediment. This enables the user to simulate small-scale density driven flows (e.g. unstable stratified flows, turbidity flows, internal waves) in coastal seas, estuaries, lakes, and rivers.
• The equation for passive tracer is also included. This means that the density of water remains unchanged, whereas the transport is only influenced by the flow and turbulent dispersion.
• Output of transport constituents is included as well (depth-averaged, per layer, and/or time-averaged).
• Both vegetation and buoyancy production rates of turbulent kinetic energy are taken into account in the standard k-ε model.
• The air pressure at sea surface is optionally included. So, with SWASH the simulation of storm surges based on space varying wind and atmospheric pressure can be carried out as well. This underlies the capability of SWASH of simulating processes at scales that range from centimeters to tens of kilometers and from seconds to hours.
• An alternative to the many well-known wind drag parameterizations (Wu, Garratt, Smith and Banke, etc.) is added. It is based on a review of a large number of more recent observations, and will gives lower drag values for relatively high wind speeds. This parameterization, based on a second order polynomial fit, has been published recently in the next article.
• The use of a semi-implicit time integration as an alternative to the explicit time stepping is added. This time integration is based on the θ-method for both the surface level gradient and the continuity equation. As a consequence, unconditional stability is achieved with respect to the celerity of gravity waves. The enhanced stability of this time stepping allows larger time steps, by a factor of five to ten, compared to the explicit one, and may thus be beneficial to large-scale applications such as tidal flow, and wind and density driven circulation.

Version 1.20:

• If the user wants so, the computational grid can be repeated in one specific direction. It means that information leaving at one end of the domain enters at the opposite end. So, the wave field is periodic in one direction with the length of the domain in that direction.
• Second order bound long waves can be added to the first order irregular waves at the wave maker.
• New output parameter BRKP which indicates areas where wave breaking occur. In this way white spots can be visualised.

Version 1.10:

• Partial reflection and transmission through porous structures, like rubble mound breakwaters in harbors, can be specified. Both submerged and emerged structures can be schematized. This also includes the interaction of the flow with the rocks in the porous structure. The implemented approach is based on the well-known Darcy and Forchheimer formulation. For details, see the User Manual.
• Logarithmic wall-law for the standard k-ε model can be employed.
• Time-averaged velocities, like undertow, and time-averaged turbulence quantities can be outputted. Also significant or RMS wave height and wave-induced setup can be plotted as well.
• A seed of the random number generator used to select phases of the Fourier components can be specified. Hence, different seeds or realisations of the same boundary wave spectrum can be obtained.
• An alternative preconditioner for solving the pressure Poisson equation is included, the so-called ILUD preconditioner. This preconditioner scales much better than the classical ILU one when running parallel.

Version 1.05:

• Effects of vertical turbulent mixing are included. These are modelled by means of the standard k-ε model.
• Output of turbulence quantities, k, ε and ν_t, is included as well.
• Oblique waves can now be imposed at the wavemaker.
• Irregular, multidirectional waves described by any 2D spectrum can now be imposed. They may be vary along the boundary as well.
• SWASH allows nesting in SWAN.

Version 1.02:

Version 10.05:

• The MUSCL scheme for vertical advective term of u/v- and w-momentum equations is the default instead of first order upwind.

Version 10.01:

Version 9.01:

• In case of bound long waves to be specified at incoming boundaries, the frequencies will be sorted (if necessary).
• An error will be produced when specifying time-averaged quantities (e.g. wave heights) at a time step different than the last one.

Version 8.01:

• Second order Hancock scheme for global continuity equation.
• Sommerfeld condition now applied to water level in case of unstructured mesh.
• Temporarily dry points not excluded from output interpolation.
• Highest file reference number is set to 99999 (see initialisation file).

Version 7.01:

Version 6.01:

• Small change to halo data in connection with dry cells.
• Master computes random wave phase and broadcasts to other processors.

Version 5.01:

Version 4.01:

• The modelling of porosity layers has been changed, as follows:
  • empirical formula's of Engelund (1953) for frictional forces are replaced by Van Gent (1995)
  • values of dimensionless friction parameters adapted accordingly
  • effect of added mass included
  • extension w-momentum equation with frictional forces
  • extension k-ε model with closures due to porosity
• The CFL criterion for explicit layer-averaged models is adapted by replacing the depth-averaged velocity by the maximum of the velocity in the water column.
• Use double precision for time coding.

Version 3.14:

• The split preconditioner ILUD is replaced by the right preconditioner ILUD.
• The weighting between ILU(D) and MILU(D) has been adapted for a better convergence.
• Small updates in the parallel code.
• Crop output of PRINT file to minimum when requested.

Version 2.00:

Version 1.20:

• In case of three equidistant layers, a numerical approximation of dispersion relation is now used.
• More flexible buffering of data exchange in case of parallel runs.
• Non-hydrostatic pressure gradient in depth-averaged mode is made conservative at discrete level.
• Pressure gradient normal to the bottom is set to zero.

Version 1.10:

Version 1.05:

Version 10.05:

Version 10.01:

Version 9.01:

Version 8.01:

Version 7.01:

Version 6.01:

Version 5.01:

Version 4.01:

Version 3.14:

Version 2.00:

Version 1.20:

Version 1.10:

Version 1.05:

Version 8.01:

• Source distribution now also hosted on GitLab repository.
• Supports CMake for building SWASH (as part of the distribution on GitLab).
• Two nonstandard Fortran 90 statements added: stream I/O (Fortran 2003 standard) and execute_command_line (Fortran 2008 standard).
• Perl script platform.pl is modified in relation to oneAPI version of Intel Fortran compiler and gfortran version 10.

Version 1.05:

The purpose of describing the bug fixes in terms of problems solved, is to enable the user to identify previous SWASH runs that may have encountered these problems (noticed at the time of running or in hindsight with this (new) information). All the bug fixes have been implemented in the current version 10.05.

Solved in version 10.05:

• also Fourier series interpreted as short waves, including sub- and super-harmonic transfer functions
• correct vertical dimension of non-hydrostatic pressure in case of subgrid approach for output
• small bugfix in imposing wave conditions at open boundaries of rotated domain
• fix of combining a sponge layer on one side and at the same time wave conditions on the other adjacent right-angled side

Solved in version 10.01:

• small bugfix in VTK format related to halo regimes

Solved in version 9.01:

• small fix to set reference water density correctly

Solved in version 8.01:

• Sommerfeld radiation condition improved: applied beforehand instead of after solving the shallow water equations
• bug fix for specification of each interpolation output quantities
• grid-oriented velocities not outputted in case of unstructured mesh
• small bug fix to cos m-model
• fix with processing of spectral boundary conditions
• correction grid partitioning with permanently dry cells
• small correction collecting post processing data of parallel run
• exhaustive memory clean-up after terminating run

Solved in version 7.01:

• small fixes to sponge layers (spherical coordinates, only wet cells)

Solved in version 6.01:

• small fix to fixed layers
• small fix with respect to update Fourier series at different boundary points
• correction semi-implicit continuity equation in case of permanently dry cells
• correction imposing two-dimensional wave spectrum on east and north boundaries with repeating grids
• inclusion orientation of non-uniform / curvilinear grids in computing wave direction of imposed 2D spectrum
• small correction computing layer interfaces in case of pressurized flow
• small correction computing Ursell number for check incoming waves

Solved in version 5.01:

• correction transformation in output wind vector
• small correction outputting friction velocity as vector
• small correction in spectrum boundary condition in case of right-handed grid
• changed value minimum depth in k-ε turbulence model
• small fix lateral advection near cycling boundaries
• small fix in log-law based wall functions

Solved in version 4.01:

• change in flux limiter function in case of linear κ-scheme
• small fix transverse velocity gradient near closed boundaries
• fix in sign for weakly-reflective boundary conditions at east and south boundaries
• in case exception value is zero, not write to Matlab as NaN
• fix in solute transport at dry points
• take into account heading per frequency in case of wave spectra at boundary
• two small fixes in time-averaging the surface elevation meant for output
• fix in computing directional spreading to preserve symmetry at boundaries
• fix in implementation of the MacCormack scheme applied to momentum equations
• small fix in log-law boundary condition

Solved in version 3.14:

• bug fix in sponge layers as applied at west and south borders
• bug fix with respect to time stepping of output files for long runs
• bug fix in sponge layers for 2D simulations
• small bug in b.c. for turbulent quantities at free surface
• small bug in outputting friction velocity
• Neumann bc for transport constituents at closed boundaries
• correction read input files in case of instationary computation or multi-layered
• correction time stamp for transport constituents
• correction Gregorian date of December 31 for years 1599, 1999, 2399, etc.
• small correction wind direction
• small correction check wetting and drying

Solved in version 2.00:

• Bug fix in energy head computation.
• Small correction in porosity in case of submerged breakwater.
• Small correction in update boundary condition in case of spherical coordinates.
• Bug fix in calculating sponge layers in case of parallel runs.
• Bug fix of computing porosity in the continuity equation.
• Bug fix in periodic boundary conditions to prevent boundary layer effect.
• Bug fix in dimension of output quantity pressure.
• Small correction in computing second order bound long waves at wavemaker boundary.
• Small correction in computing horizontal eddy viscosity.
• Small correction in searching for point in curvilinear grid.

Solved in version 1.20:

• Transformation of velocities to contravariant physical components is corrected.
• Definition of mean wave period at open boundary is corrected.
• Distinction between depth- and layer-averaged output quantities for tables is corrected.
• Two small fixes on use of random seed.
• Two small fixes on computing phase for synthesizing time series of surface elevation.
• A small correction in the interpolation procedure as applied to computational grid.
• ILU preconditioning is corrected when the box scheme is activated
• merging FRAME output data in case of parallel computing is corrected
• unexpected behaviour removed when first output time is before start of computation
• open statement for existing Matlab file which should be replaced
• Julian date conversion is corrected (in particular for years 101, 102 and 103)
• several uninitializations removed
• some other bug fixes not relevant for SWASH (mainly unstructured mesh implementation)
• No effect of weakly reflection in case of spectrum on some segments of boundary. This is corrected.
• Timestamps for stationary output parameters have been removed.
• For porosity and structure height, SWASH takes 1 and 99999 (i.e. emerged), respectively, for points outside the input grid.

Solved in version 1.10:

• A small bug fix on processing boundary conditions when running parallel.
• Some bug fixes based on bound checks, floating point exceptions, etc.

Solved in version 1.05:

• Use of random seed in case of imposing wave spectra is corrected.
• Exchanging very large data at subdomains for very large parallel runs with multiple layers is now possible.
• Handling exception values for curvi-linear grids is corrected.
• Weighting of boundary values in case of interpolation is corrected.