From Offshore to Harbor: An Integrated Multi-Scale Framework for Wave Transformation at Larache (Morocco)

Article information

J. Ocean Eng. Technol. 2026;40(3):301-324

Publication date (electronic) : 2026 June 24

doi : https://doi.org/10.26748/KSOE.2026.010

Sara Chagdali¹, Laila Mouakkir², Hassan Laarabi³, Mostafa Rachik⁴

¹Graduate Student, Laboratory of Analysis, Modeling and Simulation (LAMS); Faculty of Sciences Ben M’Sik, Hassan II University of Casablanca, Casablanca, Morocco

²Professor, Laboratory of Engineering and Materials (LIMAT), Faculty of Sciences Ben M’Sik, Hassan II University of Casablanca, Casablanca, Morocco

³Professor, Laboratory of Analysis, Modeling and Simulation (LAMS); Faculty of Sciences Ben M’Sik, Hassan II University of Casablanca, Casablanca, Morocco

⁴Professor, Moroccan School of Engineering Sciences (EMSI), Casablanca, Morocco

Corresponding author Sara Chagdali: +212522987470, sara.chagdali-etu@etu.univh2c.ma

Received 2026 February 2; Revised 2026 April 9; Accepted 2026 April 28.

Abstract

1. Introduction

Understanding and predicting the wave dynamics in coastal and estuarine environments is fundamental for safe navigation, port management, and the sustainable design of maritime infrastructure. Accurate modeling of wave transformation and harbor agitation is critical in regions where swell, tides, and currents interact, such as the northern Atlantic coast of Morocco. Over the past few decades, numerical wave models (Hervouet, 2007; Booij et al., 1999; Benoit et al., 1996; Komen et al., 1994) have transformed the ability of engineers and scientists to simulate these processes across a broad range of meteorological and oceanographic conditions. Phase-averaged (spectral) and phase-resolving approaches, ranging from third-generation wave models to mild-slope formulations, enable (1) wave propagation from offshore to coast (Booij et al., 1999; Benoit et al., 1996; Komen et al., 1994), (2) quantitative assessment of harbor agitation and operability (Berkhoff, 1972; Aelbrecht, 1997; Benoit, 2006), and (3) optimization of coastal and port protection strategies through scenario testing (Guillou & Chapalain, 2012; Chini et al., 2010; Do et al., 2022; Tan, 2023).

Despite their broad applicability, the predictive performance of these models remains strongly dependent on how numerical and physical processes are represented. Recent studies have examined the sensitivity of spectral wave models to their main input parameters. Park (2009), through applications of SWAN and TOMAWAC along the East Anglian coast, showed that frequency resolution has only a limited influence on model predictions, even though the coarsest frequency and directional discretization may lead to an underestimation of the significant wave height and mean wave period. By contrast, source and sink parameterizations associated with bottom friction and depth-induced breaking exert much stronger control over the simulated wave field. The same study also showed that whitecapping and quadruplet interactions have minor effects, whereas triad interactions transfer energy to higher frequencies and tend to produce shorter and steeper waves, especially in shallow-water environments.

These results show that model selection cannot be based solely on computational convenience, but must instead reflect the dominant physical processes to be resolved at each spatial scale. While phase-averaged models are generally well suited for offshore and nearshore wave propagation, phase-resolving models become necessary in harbor environments where diffraction, reflection, and wave–structure interactions strongly govern the local wave field. In this regard, Do et al. (2022) compared TOMAWAC, a phase-averaged spectral wave model, with ARTEMIS, a phase-resolving model within the TELEMAC–MASCARET system (Hervouet, 2007), to assess their respective abilities to reproduce wave diffraction in the presence of protective structures. They reported that ARTEMIS provided higher accuracy than TOMAWAC at the innermost locations for all tested wave heights and incident directions, whereas no substantial improvement in wave-height prediction was observed outside the shadow zone. This confirms that the added value of phase-resolving approaches is greatest in sheltered harbor areas where diffraction-dominated processes become significant.

Beyond the relative performance of individual models, the accurate simulation of wave transformation from offshore waters to inner-harbor basins requires a consistent multi-scale coupling strategy. Offshore propagation, nearshore transformation, and harbor agitation are governed by different physical processes and spatial scales. Hence, no single model is generally sufficient to represent the entire wave pathway with equal accuracy. Therefore, robust coastal and port modeling increasingly relies on coupled frameworks that combine the respective strengths of phase-averaged and phase-resolving approaches. In this context, Guillou and Chapalain (2012) emphasized the importance of nested computational domains for predicting wave heights within the inner basin of Le Havre harbor. Their study showed that coupling a regional phase-averaged model with a harbor-scale phase-resolving model provides an efficient and physically consistent framework for assessing harbor agitation, with the large-scale solution supplying boundary conditions to the inner high-resolution domain.

Although these studies have provided important insights, several gaps remain. Park (2009) treated sensitivity analysis mainly as a standalone exercise rather than as one component of a complete offshore-to-harbor modeling chain. In addition, the assessment was based primarily on comparisons at discrete observation points, which, although informative for local model behavior, provide only a limited view of the spatial structure of sensitivity effects. Guillou and Chapalain (2012) revealed the relevance of model coupling for tide-modulated wave penetration, but did not systematically address parameter sensitivity or explicit chaining criteria. Do et al. (2022) clarified the respective capabilities of TOMAWAC and ARTEMIS in a diffraction-dominated harbor environment, but remained focused on the comparative model performance at a specific site. The present study addresses these limitations by extending the sensitivity perspective of Park (2009) through an analysis of two-dimensional spatial distributions of significant wave height (H_m0), by integrating tide- and current-induced wave modulation explicitly through TELEMAC-2D, and by introducing physically based nesting criteria between TOMAWAC and ARTEMIS using the Miche breaking-steepness criterion and the Ursell number.

Beyond its scientific interest, the development of a consistent offshore-to-harbor modeling framework has strong practical implications for harbor safety, navigation risk assessment, and design optimization. Accurate prediction of wave penetration and agitation is necessary to evaluate vessel accessibility at harbor entrances, define safe operational windows, and identify zones where excessive wave disturbance may compromise mooring stability, loading operations, or basin operability. Such a framework is equally important for testing and optimizing engineering solutions, including breakwater configuration, harbor layout, and coastal protection strategies. In a broader perspective, integrated modeling chains of this kind are also fully aligned with the emergence of digital coastal management systems, in which interoperable numerical tools support scenario analysis, operational decision-making, and adaptive management of port and coastal environments.

Unlike validation-oriented investigations, this study does not aim to revalidate the numerical models themselves, which have already been extensively tested and verified in a wide range of configurations (Chini et al., 2010; Tan, 2023; Zhang et al., 2024; Machado et al., 2021; Scaldis Coast Project Team, 2023). No long-term in situ wave record suitable for direct validation was available for the Larache study area, and no direct field measurements were available for the present study. Instead, this study evaluated a comprehensive modeling methodology that ensures physical consistency, numerical robustness, and interoperability across TELEMAC–MASCARET modules. This approach reflects a methodological, rather than purely predictive, objective, aiming to establish reproducible procedures for coastal and port studies where direct observations may be scarce. Under these conditions, confidence in the results is supported by long-term quality-controlled offshore forcing data, by consistency with the expected physics of wave transformation in shallow coastal environments, and by comparisons with the conclusions from previous studies on spectral sensitivity, harbor-wave penetration, and tide-modulated wave propagation. This methodology can support future operational modeling, design optimization, and coastal-risk mitigation in Moroccan and Atlantic port environments.

Within this context, this study developed and applied an integrated, multi-scale numerical framework for the Port of Larache (Morocco), a dynamic estuarine site influenced by Atlantic swells and tidal currents. The modeling chain combines: (1) the spectral wave model TOMAWAC, systematically assessed in terms of mesh resolution, spectral/directional discretization, and source-term parameterization; (2) the hydrodynamic model TELEMAC-2D, resolving the tidal and regional current fields that modulate nearshore transformation (Hervouet, 2007); and (3) the phase-resolving model ARTEMIS, applied within the harbor to simulate diffraction, reflection and port agitation. The transition between the TOMAWAC and ARTEMIS domains is guided by the Miche breaking-steepness criterion and the Ursell number to ensure physically consistent nesting between computational scales (Miche, 1944; Ursell, 1953).

Accordingly, the objectives of this research are to

2. Study Area

The maritime front of Larache (Fig. 1) is a coastal zone in northwest Morocco, along the Atlantic Ocean. The data collected in this region can contribute to sustainable management strategies, essential for protecting coastal ecosystems and promoting local economic activities. This site is positioned at geographic coordinates 35°11’ North latitude and 06°10’ West longitude. The underwater beach has a gentle slope extending to the −10-meter isobath, where depressions form to depths of −15 meters. It is characterized by rich marine biodiversity and diverse coastal ecosystems (Saddik et al., 2024).

Fig. 1

Location of the Larache port site: (a) geographical position along the Moroccan Atlantic coast; (b) detailed satellite view of the port area and Loukkos River estuary

The wave forcing dataset used in this study was extracted from the Puertos del Estado SIMAR archive at the selected offshore point and spans a long-term period from 1958 to 2026, providing a robust basis for characterizing the long-term wave climate and distinguishing recurrent background conditions from low-frequency energetic events relevant to coastal and harbor applications (Puertos del Estado, 2024a, 2024b). After extraction, the dataset was quality-controlled by removing invalid or missing values and then statistically processed to derive a set of complementary descriptors, including long-term summary statistics, monthly climatology, annual maxima, the joint occurrence distribution of H_m0 and T_p, and directional wave roses.

Fig. 2 presents the temporal evolution of H_m0 over the full record, together with the monthly mean and the annual maxima. The time series reveals a persistent background regime dominated by moderate wave activity, with episodic high-energy events superimposed. Physically, this suggests that the site is most of the time subjected to relatively stable swell-dominated forcing, while energetic storms intermittently inject much larger wave energy into the system. The annual-maxima series confirms that energetic storms recur throughout the record and that extreme conditions are not isolated anomalies but an intrinsic part of the long-term wave climate. The highlighted event of 10 February 2014 is representative of this upper-energy regime, with H_m0 = 6.16 m, T_m = 9.51 s, T_p = 14.3 s, and a mean direction of 241°. This event illustrates the type of severe sea state that governs the upper tail of the forcing conditions relevant for design and modeling. Therefore, it was retained as the reference high-energy wave condition for the numerical simulations performed in this study.

Fig. 2

Long-term temporal variability of the significant wave height with the highlighted extreme event

Fig. 3 presents the joint distribution of H_m0 and T_p in terms of the occurrence frequency. The figure shows that the most recurrent sea states are associated with moderately significant wave heights, typically around 0.5–1.5 m, and peak periods mainly ranging from approximately 5 to 12 s, with a modal concentration around 9–11 s. This confirms that the long-term wave climate is dominated by persistent moderate swell conditions rather than by short-period locally generated wind waves alone. From a physical perspective, the prevalence of these intermediate-to-long periods suggests that offshore wave forcing is largely governed by remotely generated Atlantic swell, which can propagate over long distances with limited dissipation before reaching the study area. At the same time, the distribution exhibits a marked upper tail extending toward larger H_m0 and longer T_p, highlighting the presence of rarer but physically important storm events. These extreme states occur infrequently in percentage terms, but they are essential for defining the severe hydrodynamic forcing conditions relevant to exposed coastal and harbor structures. Long-period energetic waves are expected to induce stronger shoaling, enhanced wave setup, and potentially deeper penetration into semi-sheltered basins compared to shorter locally generated sea states.

Fig. 3

Joint distribution of significant wave height H_m0 and peak period T_p

Fig. 4 shows the overall wave rose derived from the same long-term record. The wave climate is strongly directional, with a dominant contribution from the western to northwestern sector, particularly from the WNW–NW quadrant. This suggests that the site is primarily exposed to incoming Atlantic wave propagating from the offshore western sector. From a physical perspective, this directional concentration suggests that the regional coastline orientation and offshore exposure favor a preferential incident-wave corridor, such that most of the wave energy reaching the area is delivered within a relatively narrow directional band. In terms of wave-height classes, the most frequent conditions fall within the 0–1.2 m and 1.2–2.5 m ranges, whereas more energetic waves, especially those exceeding 4 m, occur much less frequently. Therefore, the wave rose confirms the strong anisotropy of the wave forcing and the predominance of moderate sea states, while also identifying the directional window associated with the most energetic events.

Fig. 4

Wave rose showing the occurrence of significant wave height classes

Hydrographic zero is 2.14 meters below the General Leveling of Morocco (NGM). Table 1 lists the tidal amplitudes in Larache.

Table 1

Tidal amplitudes in the Larache region from Directorate of Maritime Public Domain, DPDPM (2024)

Bathymetry of the coastal area of Larache was developed based on two sources of nautical charts. The first chart was provided by SHOM (2012) (Fig. 5) (a), and the second, with greater accuracy near the coast, was provided by DPDPM (2024) (Fig. 5 (b)). These two charts offer a precise representation of maritime depths and seabed characteristics. Fig. 3 presents both charts, showing the depths of the estuarine area in Lambert 93 coordinates.

Fig. 5

Bathymetry for modeling: (a) SHOM (2012), (b) Bathymetry from DPDPM (2024)

Fig. 6 presents current fields derived from the Copernicus Marine Environment Monitoring Service (CMEMS) global reanalysis product for the Larache shelf on 16 June 2024 (Copernicus Marine Service, 2024), resolved at six upper-layer levels spanning 0.5–6.4 m depth. Surface velocities exceed 0.40 m s⁻¹ in the upper 2 m, and the jet progressively weakens with depth while remaining detectable down to 6.4 m. Offshore, current speeds are much weaker (<0.15 m/s) and a broad anticyclonic recirculation cell is apparent. These hydrodynamic features have direct implications for harbor engineering at Larache; the along-shelf jet modulates the incident wave field through current–wave interaction.

Fig. 6

Ocean circulation at different depths off the coast of Larache (Copernicus Marine Service, 2024)

3. Numerical Models

The TOMAWAC, TELEMAC-2D, and ARTEMIS modules from the TELEMAC-MASCARET system (TMs) (www.opentelemac.org) are used in this study. TMs is a software suite designed for the numerical modeling of free-surface hydraulics, sediment transport, and wave dynamics in 1D, 2D, or 3D. TMs uses the finite element method with unstructured grids. This section provides a brief description of the TOMAWAC, TELEMAC-2D, and ARTEMIS modules.

TOMAWAC is a third-generation spectral wave model phase-averaged wave propagation model that accounts for various phenomena, including wave generation, refraction, shoaling, nonlinear wave interactions, and energy dissipation due to bottom friction and rapid depth variation. The TOMAWAC model solves the wave energy density conservation equations, which can be expressed in the general form.

where N is the action density spectrum of the wave; σ is the intrinsic frequency; θ is the wave direction; C_g is the group velocity in space (x, y, σ, θ); and S_tot represents the total source terms. The TOMAWAC code considers various processes as source terms in wave dynamics. These source terms include Q₂₅ (wind-driven wave generation), Q_ds (whitecapping-induced energy dissipation), Q_nl (nonlinear quadruplet interactions), Q_bf (bottom friction-induced energy dissipation), Q_br (bathymetric breaking-induced energy dissipation), Q_tr (nonlinear triad interactions), and Q_ds;cur (enhanced dissipation due to waves interacting with currents). These terms govern the generation, redistribution, and dissipation of wave energy in offshore and coastal environments. The availability of different models for each source term allows for flexibility in simulating these processes under various conditions. Nevertheless, the application domain of TOMAWAC excludes harbor areas and, more generally, situations where the effects of reflection on structures and/or diffraction cannot be neglected, such as diffraction by breakwaters, piers, or shoals, and reflection, whether partial or total, from structures or significant depth irregularities. The following formulations, as implemented in the TOMAWAC user version 8.4 (TELEMAC-MASCARET system), were used in the present study:

TELEMAC-2D solves the shallow-water equations, derived by depth-integrating the Navier–Stokes equations. The equations below are expressed in Cartesian coordinates. The first equation represents mass conservation (continuity equation), while equations 2 and 3 correspond to the momentum equations in the x and y dimensions, respectively.

where t is time; H is the flow depth; u and v are depth-averaged velocities in x and y directions, respectively; u is the depth-averaged velocity vector (u, v); Z_b is the bed elevation; g is the acceleration due to gravity; v_t is momentum diffusion coefficient; ∇ is the gradient operator (∂/∂x, ∂/∂y), and ∇· is the divergence operator. S_h is a fluid source or sink, and S_x and S_y are momentum source terms in the x and y directions, respectively, encompassing factors such as wind, Coriolis force, bottom and lateral friction, wind stress, atmospheric pressure gradients, and wave forces. Considering only the bottom friction, S_x and S_y are expressed as follows, in which C_f is a dimensionless friction coefficient that can be calculated according to the Strickler equation:

For a direct coupling between TELEMAC-2D and TOMAWAC, both models must share the same horizontal mesh and be run over identical simulation durations, defined by the time step and the number of iterations. The time steps must either be identical or set as integer multiples of one another to ensure consistent data exchange. In addition, TOMAWAC cannot be forced with external current or water level files. Instead, it must receive hydrodynamic inputs directly from the coupled TELEMAC simulation. Finally, the driving forces along the x and y directions (FX, FY) must be included among the TELEMAC output variables specified in the steering file to allow proper transfer of momentum fields to TOMAWAC.

ARTEMIS is another wave module within the TMs (Hervouet, 2007). The main application of ARTEMIS concerns wave agitation in ports or small bays. The ARTEMIS module is a phase-resolving model that solves the mild-slope elliptic equation (Berkhoff equation) using the finite element method. This equation is derived from the Navier–Stokes equations under a set of assumptions and can be expressed in the general form.

where ø is the potential velocity; C and C_g denote the phase and group velocities, respectively; k is the wave number: k=2π/L; and ∇ and ∇ denote the gradient and divergence operators, respectively, as restricted to the horizontal coordinates. Nevertheless, the mild-slope model is linear, highlighting its limitations in shallow water. When the seabed rises, wave breaking amplifies the wave heights, leading to inconsistencies with the varying water depth. The original mild-slope equation was modified to account for bathymetric wave breaking and bottom friction, and improve the validity of the model and expand the applicability of ARTEMIS.

4. Methodology

A key component of the proposed methodology involves a multi-model coupling approach to accurately simulate wave propagation and its interactions within the port environment.

In the first step, a standalone TOMAWAC simulation is carried out over a 5,800 s period with a 1 s time step to assess how key numerical choices affect wave propagation. All sensitivity begins from a common reference configuration: a coastal mesh refined to a uniform 150 m resolution in the nearshore zone, a spectral discretization of 15 logarithmically spaced frequency bins, and 15 equally spaced directional bins. The model is initialized with still water at +2.15 m (H_s = 0 m) and forced at the offshore boundary by a Jonswap directional spectrum (H_s = 6.16 m, peak frequency f_p = 0.07 Hz, principal direction 240°). The significant wave height Hm0=4m0, derived from the spectral moment of order zero, is equivalent to the observed significant wave height H_s under the Rayleigh distribution assumption, and is used throughout this study as the representative wave height. Energy losses are represented by bottom friction (coefficient C_f = 0.038) and depth-induced breaking according to the Battjes and Janssen model (Battjes & Janssen, 1978) –Formula 1, while wind input is neglected. This baseline run provides the foundation against which spatial-mesh, frequency-bin, and directional-bin refinements, as well as source-term variants, are compared. After defining this baseline, each element is systematically varied as summarized in Table 2, which lists the analysis elements considered and the corresponding variants tested.

Table 2

Sensitivity analysis elements and variants tested

All runs remain 5,800 s long, and in each case, spatial ΔH_m0 maps (difference from the reference mesh) were generated, and the steady-state H_m0 and mean wave direction time series were extracted at four control nodes. The influence of each numerical choice was quantified by calculating three complementary statistical indicators against the reference configuration: the root-mean-square error (RMSE), the mean bias (Bias), and the 95^th percentile of the absolute error (P95).

The root-mean-square error (RMSE) between a test configuration and the reference configuration is defined, for a given variable X (e.g., H_m0 or mean wave direction), as follows:

where N is the number of points considered; Xitest is the value at point i for the test configuration; and Xiref is the corresponding value for the reference configuration.

The mean bias, which measures the systematic overestimation or underestimation of the reference solution by the test configuration, is given by

A positive value indicates an overall overprediction, whereas a negative value indicates an overall underprediction.

The 95^th percentile of the absolute error was also used to characterize the upper bound of the discrepancy distribution in a robust manner:

This indicator represents the error threshold below which 95% of the absolute differences fall, providing a practical measure of near-extreme mismatch that is not overly influenced by isolated outliers. For the mean wave direction, the error was calculated using the minimum angular difference to account for the circular nature of directional data.

An updated reference case was set once the mesh, frequency, and directional discretization analyses were complete. Fig. 7 shows the unstructured triangular mesh used for all TOMAWAC runs (31,708 triangles, 16,389 nodes) together with the four control points selected for time-series extraction. These nodes represent (1) deep-water conditions, (2) intermediate depth on the mid-shelf, (3) the nearshore platform, and (4) the harbor entrance. For each sensitivity case, the steady-state values of H_m0 and the mean wave direction were sampled at these locations to calculate the pointwise RMSE against the reference configuration. Therefore, a final configuration was retained for the remainder of the study.

Fig. 7

2D mesh with observation nodes

Once the sensitivity analysis of the TOMAWAC configuration was completed, an additional sensitivity test was performed on the offshore wave boundary conditions. This analysis considered perturbations of ±10% in H_m0 and ±5% in the wave direction to assess the robustness of the model response to plausible uncertainties in the incident forcing.

In the second step, the optimal TOMAWAC configuration was coupled with TELEMAC-2D to simulate wave propagation under tidal forcing and regional currents. This coupling incorporated water-level variations and tidal currents, ensuring that large-scale wave dynamics and localized tidal effects had been captured with sufficient realism before defining the transition interface to ARTEMIS for harbor-scale diffraction and reflection studies.

Thereafter, the focus shifted to defining the transition interface between the TOMAWAC and ARTEMIS simulations. TOMAWAC was applied from the offshore domain to the coastal region, where bottom-interaction processes and wave-breaking phenomena become significant. Its formulation accounted for key source terms, including whitecapping dissipation (Q_ds), nonlinear quadruplet interactions (Q_nl), and triad interactions (Q_tr). The relative importance of these mechanisms varied with depth. Whitecapping and quadruplet interactions dominated in deep water, while triad interactions governed energy redistribution in shallow coastal areas, driving shoaling, breaking, and refraction processes. A physically consistent interface between the two models must be defined to ensure that ARTEMIS is applied in the nearshore zone where phase-resolving processes such as diffraction and reflection prevail. Two complementary indicators were used for this purpose. The Miche criterion (Miche, 1944) was then applied to determine the limiting condition for wave-breaking. In practice, the criterion is simplified to the widely adopted formulation H_b/h=0.78, where H_b denotes the breaking wave height and h the local water depth. Thereafter, the Ursell number (Ursell, 1953), expressed as U_r=HL²/H³ (where H is the wave height, L the wavelength, and h the local water depth), was used to characterize the degree of nonlinearity in shallow-water wave transformation, identifying areas where dispersive effects become less significant relative to nonlinear processes. The linear dispersion relation ω² = g k tanh(kh) was solved for the wavenumber k (Newton iteration) and to obtain L = 2π/k using H_m0 as the operational wave height H and the peak frequency f_p to calculate the angular frequency ω = 2π f_p. The joint use of these two criteria provides a robust physical basis for positioning the ARTEMIS seaward boundary just shoreward of the zone where nonlinear transformations and wave breaking become dominant. In this way, the ARTEMIS computational domain is explicitly constrained to the region where diffraction, reflection, and port-induced agitation processes control the hydrodynamic response.

The modeling framework adopts a hierarchical multi-scale approach, combining a two-way coupling between TELEMAC-2D and TOMAWAC across the regional and intermediate domains with a one-way spectral nesting from TOMAWAC to ARTEMIS in the harbor domain. In the coupled TELEMAC-2D–TOMAWAC configuration, TELEMAC-2D updates the free-surface elevation, local water depth, and depth-averaged current field at each 1 s coupling time step. These hydrodynamic variables were then transferred to TOMAWAC, which solved the wave action balance under non-stationary water-level and current conditions. In return, TOMAWAC calculates the wave-induced momentum forcing components (F_x and F_y), which are fed back into TELEMAC-2D, ensuring a physically consistent two-way representation of wave–current interaction. At the inner scale, the directional wave spectrum obtained from TOMAWAC was extracted along the offshore boundary of the ARTEMIS mesh and imposed as the incident forcing. This nesting strategy preserved the large-scale transformation processes already resolved by TOMAWAC, including shoaling, refraction, dissipation, and current-induced modulation, while exploiting the phase-resolving capabilities of ARTEMIS to simulate diffraction, reflection, and harbor agitation. The robustness of the coupling procedure was ensured using dynamically compatible exchanged variables, consistent time discretization in the TELEMAC-2D–TOMAWAC system, and careful positioning of the TOMAWAC–ARTEMIS interface in a region where the incident wave field remains regular and physically interpretable as boundary forcing. The Miche and Ursell criteria were also used to confirm that the nesting boundary lay within an appropriate transition zone, reducing the risk of artificial discontinuities between the spectral and phase-resolving modeling frameworks.

5. Results

5.1. TOMAWAC results

5.1.1. Reference case

Under the reference-case conditions outlined in the methodology, an initial simulation showed that the wave train propagates shoreward while gradually losing height as the water depth decreases (Fig. 8). Wave rays refract until they become nearly orthogonal to the shoreline. The refraction coefficient between the offshore section and the estuary mouth is approximately 0.92, indicating only a modest change in direction. Inside the estuary, the wave heights drop sharply due to the shallow bathymetry (Fig. 6) and the sheltering effect of maritime structures so that only weak agitation reaches the harbor basin.

Fig. 8

Reference Case - Distributions of wave heights H_m0, for H_s = 6.16 m, T = 14 s, Dir = 240°

5.1.2. Mesh analysis

Fig. 9 compares the differential significant wave height fields (Δ H_m0) obtained with coastal grid resolutions of 80 m, 100 m, and 200 m with the reference 150 m mesh. In all cases, the color scale was dominated by dark blue tones, indicating point-by-point deviations below 0.10 m, i.e., less than 2 % of the imposed offshore significant wave height (H_s = 6.16 m). Such small discrepancies highlight the robustness of the numerical solution with respect to the spatial resolution. Slightly brighter patches appeared in the shoaling zone seaward of the breakwater headland, where the bathymetric gradient enhanced the local sensitivity to grid discretization. On the other hand, no systematic trend was observed. The refined 80 m and 100 m meshes did not artificially increase the crest amplitudes, while the coarser 200 m grid did not induce additional damping along the coastal slope. The overall wave propagation pattern remained unaffected, preserving the locations of refraction cells, interference features, and the lee-side energy shadow. These findings showed that the reference discretization Δx = 150 m already provided sufficient accuracy for the regional-scale objectives of this study, ensuring a faithful representation of the dominant processes without requiring computationally expensive refinement. Therefore, the 150 m mesh was retained as the optimal compromise between numerical accuracy and computational efficiency for the subsequent simulations, a choice that will be further corroborated through pointwise comparisons at selected control stations.

Fig. 9

Spatial difference maps of the significant wave height H_m0 for H_s = 6.16 m, T = 14s, Dir = 240°. Comparison between the reference mesh (Δx = 150 m) and three alternative resolutions: 80 m, 100 m, and 200 m

Fig. 10 compares the significant wave height (H_m0) and mean wave direction at four control points for four mesh sizes: 80 m, 100 m, 150 m (reference), and 200 m. The RMSE values remained generally low for both variables, even though their magnitudes varied with the local hydrodynamic conditions. For Hm0, the RMSE ranged from 0.018 to 0.028 m at Node 1, 0.029 to 0.049 m at Node 2, 0.038 to 0.095 m at Node 3, and 0.036 to 0.059 m at Node 4. The largest deviations occurred at Node 3, on the nearshore platform, indicating greater sensitivity to the mesh resolution in shallow-water transformation zones. For mean wave direction, circular RMSE remained below 1° at all nodes, varying from 0.22° to 0.38° at Node 1, 0.28° to 0.60° at Node 2, 0.51° to 0.74° at Node 3, and 0.05° to 0.18° at Node 4. These consistently small directional errors suggest that the predicted propagation direction is only marginally influenced by grid resolution, whereas the wave height is more responsive in shallow and morphologically complex areas. Overall, these results support the use of the 150 m mesh as a reasonable choice for numerical robustness in the present configuration.

Fig. 10

Comparison of the wave height H_m0 and mean direction at four control points for four different mesh sizes: 80 m, 100 m, 150 m (reference), and 200 m

5.1.3. Number of frequencies discretization

Although the reference configuration was originally defined with N_f = 15, the results are presented relative to the higher-resolution case N_f = 30 for clarity. Fig. 11 shows spatial ΔH_m0 maps for N_f = 15, 20, and 25 referenced to N_f = 30. All three maps were almost uniform, with very low values, indicating that the spatial response of Hm0 is only weakly sensitive to frequency discretization. No pronounced band, coherent error structure, or systematic increase in differences was observed across the domain, including in the nearshore zone. Any deviations remained highly localized and negligible in magnitude. This confirmed that frequency resolution has a very limited influence on the simulated wave-height field in the present setup, and that spectral convergence is effectively achieved with relatively moderate values of N_f.

Fig. 11

Spatial difference maps of significant wave height H_m0. Comparison between the reference configuration (N_f = 30) and three alternative frequency discretization: N_f = 15, 20, and 25

The discrepancies were negligible for all tested directional discretizations relative to the reference case Δθ = 15°, as shown in Fig. 12. The RMSE in Hm0 remained between 0.008 and 0.036 m across all nodes for Δθ = 20°, 25°, and 30°, i.e., below approximately 0.6% of the propagated incident wave height (H_m0 = 6.16 m). Similarly, the circular RMSE in mean wave direction remained below 0.25° at all locations. These results show that directional discretization has no significant influence on the simulated wave field in the present setup, and that Δθ = 15° is already sufficient to achieve a converged solution.

Fig. 12

Comparison of the wave height H_m0 and mean direction at four control points for four different frequency discretizations: N_f = 15 (reference), 20, 25, and 30

5.1.4. Number of directions discretization

For clarity, the spatial sensitivity to directional discretization is presented relative to the finest tested case (N_d = 35). Fig. 13 shows that the largest ΔH_m0 patterns occur for N_d = 15, mainly within the shoaling region and in the vicinity of the harbor entrance and breakwater, where directional spreading, refraction, and local energy redistribution are most sensitive to the angular resolution. The extent and magnitude of these discrepancies decrease progressively for N_d = 20 and N_d = 25, indicating a clear convergence of the simulated wave field as the directional spectrum is refined. Although the N_d = 30 solution is already very close to the N_d = 35 reference over most of the domain, residual differences remain visible in the most dynamically sensitive nearshore and harbor areas. In view of this behavior, N_d = 35 was retained as the final directional discretization to ensure the most robust and physically consistent wave field for the subsequent modeling steps.

Fig. 13

Spatial difference maps of significant wave height H_m0 for H_s = 6 m, T = 14 s, and wave direction 240°. Comparison between the reference configuration (N_d = 35) and three alternative frequency discretization: N_d = 15, 20, 25 and 30

Fig. 14 shows that directional discretization affects both significant wave height and mean wave direction, with a clear overall convergence trend as the number of directional bins increases toward the reference case (N_d = 35). For H_m0, the RMSE decreased from 0.179 to 0.006 m at Node 1, from 0.062 to 0.011 m at Node 2, from 0.077 to 0.013 m at Node 3, and from 0.038 to 0.018 m at Node 4 as N_d increased from 15 to 30. A similar reduction was observed for the mean wave direction, with circular RMSE decreasing from 2.42° to 0.18° at Node 1, from 0.69° to 0.27° at Node 2, from 2.82° to 0.48° at Node 3, and from 1.05° to 0.34° at Node 4. Although minor non-monotonic variations remained at some intermediate nodes, the overall behavior showed that coarse directional discretizations are less able to reproduce the reference solution, whereas higher N_d values substantially improve the wave-height and directional predictions. This pointwise convergence was consistent with the spatial ΔH_m0 patterns, which showed progressively weaker residuals as the directional resolution was refined. On this basis, N_d = 35 was retained as the reference directional discretization for the subsequent modeling steps.

Fig. 14

Comparison of wave height H_m0 and mean direction at four control points for four different directions of discretization: N_d = 35 (reference), 15, 20, 25, and 30

5.1.5. Source terms

Fig. 15 shows the sensitivity of the simulated significant wave height to different physical source-term parameterizations, expressed as differences in ΔH_m0 relative to the reference configuration in which wave breaking dissipation is represented by the Battjes and Janssen (1978) model (WBD – Formula 1). The introduction of triad interactions produced the most pronounced deviations, with local increases in wave height exceeding 0.40 m along the shoaling zone seaward of the harbor. This highlights the strong role of nonlinear triad interactions in redistributing the spectral energy in shallow water, leading to locally enhanced crest heights compared to the reference solution. In contrast, alternative formulations of whitecapping dissipation (F1 vs. F2) yielded negligible differences. The ΔH_m0 field remained almost uniformly close to zero, suggesting that deep-water whitecapping played only a marginal role in this coastal configuration. Similarly, the two tested quadruplet interaction formulations produce almost indistinguishable results, with pointwise differences well below 0.05 m, confirming that deep-water quadruplet interactions are not a controlling process at this nearshore scale. The largest discrepancies were observed for the wave-breaking dissipation tests. When bathymetric breaking-induced dissipation was switched off (Formula 0), ΔH_m0 exceeded 0.4 m in the surf zone and near the breakwater, reflecting a clear overestimation of wave energy in the absence of depth-induced breaking. Alternative dissipation formulations (F2, F3, and F4) lead to small deviations relative to the Battjes and Janssen (1978) model, mostly confined to the surf zone and within the estuary, with the reference formulation maintaining the most balanced dissipation pattern. Overall, these results confirmed that (1) nonlinear triad interactions can significantly modify the nearshore energy distribution; (2) alternative whitecapping and quadruplet schemes have a negligible impact in this context; (3) bathymetric breaking dissipation is essential, with the Battjes and Janssen (1978) model parameterization providing a physically consistent reference for this site.

Fig. 15

Two-dimensional ΔH_m0 fields for each source-term configuration

Fig. 16 summarizes the pointwise sensitivity to source-term parameterizations for (left) significant wave height H_m0 and (right) mean direction at four control nodes. Deep water (Node 1) was insensitive, where H_m0 remained ~5.25 m across all tests, and the mean direction varied by <1°, confirming that offshore conditions are dictated by boundary forcing rather than by nearshore source-term choices. At intermediate depth (Node 2), the changes are still minor (≈5.1 m; direction ~77°), indicating weak sensitivity before significant shoaling. The sensitivity increased on the nearshore platform (Node 3): introducing triad interactions (triad1) raised H_m0 to ~4.46 m compared with ~4.0–4.1 m in the other runs, which was consistent with the role of nonlinear triads in shallow-water energy redistribution and local crest enhancement. At the harbor entrance (Node 4), the wave response was controlled by depth-induced breaking. The reference simulation, using the Battjes and Janssen (1978) model (Formula 1), yields H_m0≈ 2.97 m. Suppressing bathymetric breaking (WBD0) led to a non-physical increase in the local energy level (H_m0 ≈ 4.1 m), highlighting the necessity of including depth-limited dissipation in this zone. All alternative breaking configurations (WBD2–WBD4) produced solutions that were more dissipative than the Battjes & Janssen reference. Of the tested options, WBD2 implementing the Thornton and Guza (1983) formulation (Formula 2) produced the strongest dissipation (H_m0 = 2.19 m). WBD4, corresponding to Izumiya and Horikawa (1984) (Formula 4), also substantially lowered the energy (H_m0 ≈ 2.48 m). By contrast, WBD3, representing Roelvink-type breaking (Roelvink, 1993) (Formula 3), delivered the solution that departs the least from the Battjes and Janssen reference, with H_m0 = 2.83 m, even though it remains marginally over-dissipative. In contrast, the mean wave direction was effectively invariant at all four control points (≈71°, 77°, 93°, and 113° for Nodes 1–4, respectively, with departures <1–2°), confirming that, for this configuration, the choice of source-term parameterization primarily affected the energy balance rather than the directional field. Overall, the pointwise diagnostics corroborate the map-based findings: (1) triads materially affect the shallow-water energy; (2) whitecapping and quadruplet schemes have negligible impact in this setting; and (3) bathymetric breaking dissipation is essential, with the Battjes and Janssen (1978) model formulation providing the most balanced nearshore/entrance energy levels.

Fig. 16

Comparison of the wave height H_m0 and mean direction at four control points, four control nodes, and associated RMSE diagnostics

Fig. 17 shows that the three configurations with triad interactions (“Quadruplet + Triad”, “Whitecapping + Triad”, and “Whitecapping + Triad + Quadruplet”) all produced a similar nearshore response: a narrow, shore-parallel high-energy core along the shoaling zone and around the harbor, with ΔH_m0 > 0.40 m. This suggests that triads alone are sufficient to generate coastal amplification. “Quadruplet + Triad” and “Whitecapping + Triad” were almost indistinguishable, so adding Quadruplet or whitecapping on top of triads does not remove the triad-driven peak. When all three processes were activated, the core pattern was unchanged, but the signal outside the red band was slightly weaker, suggesting that the full combination (triad + quadruplet + whitecapping) limits the lateral spread of the anomaly rather than its maximum. In contrast, the only triad-off run (“Whitecapping + Quadruplet”) showed only modest increases (ΔH_m0≈ 0.10–0.15 m) and no harbor hotspot, confirming that triads are the decisive term, while whitecapping and quadruplet act more on the spread than on the peak.

Fig. 17

Two-dimensional ΔH_m0 fields for mixing terms

Fig. 18 confirms the dominant role of triad interactions. At Nodes 1 and 2 (offshore and pre-shoaling), H_m0 remained relatively unchanged (≈5.3–5.4 m and ≈5.1–5.2 m), suggesting that deep-water conditions were insensitive to the tested source-term combinations. The contrast was explicit in shallow water (Nodes 3 and 4). At Node 3, H_m0 ≈ 4.0–4.1 m in the reference and whitecapping–quadruplet cases, but increased to ≈4.6 m when triads are activated (Q+T, W+T). At Node 4, the case without triads drops from approximately 4.0 m at Node 3 to 2.97 m, while the cases with triads dropped from approximately 4.6 m to 2.8 m, so the ordering observed at Node 3 was reversed by depth-induced breaking. The mean wave direction remained effectively invariant at all nodes (<1–2°). The pointwise results were consistent with the H_m0 maps, highlighting the effects of triads more clearly, but they did not reflect the slight whitecapping–quadruplet effects visible on Fig. 15.

Fig. 18

Comparison of the wave height H_m0 and the mean direction at four control points, four control nodes, and associated RMSE diagnostics

Accordingly, a baseline setup consisting of a coastal mesh refined to a uniform 150 m resolution in the nearshore zone was retained for the remainder of the study: a spectral discretization of 20 logarithmically spaced frequency bins and 25 equally spaced directional bins; initialization with still water at +2.15 m (H_s = 0 m); and offshore forcing by a JONSWAP directional spectrum (H_s = 6.16 m, peak frequency f_p = 0.07 Hz, principal direction 60°). Energy losses are represented by bottom friction with C_f = 0.038 and depth-induced breaking. Nonlinear triads, whitecapping, and quadruplets are all included using their Formula 1 implementations.

5.1.4. Quantitative synthesis of sensitivities

Fig. 19 presents a strict pointwise sensitivity summary based on three complementary indicators: RMSE, Bias, and P95. In this representation, the bubble label gives the RMSE, the color indicates the sign and magnitude of the bias, and the bubble size reflects the 95^th percentile of the absolute discrepancy. Overall, these metrics provide a compact but robust interpretation of model sensitivity, distinguishing between the average deviation, systematic over- or underestimation, and the magnitude of the largest local departures.

Fig. 19

Comparison of the wave height H_m0 and mean direction at four control points, four control nodes, and associated RMSE diagnostics

Fig. 19 shows the clear hierarchy of sensitivities when RMSE, Bias, and P95 were analyzed jointly along the offshore-to-harbor sequence defined by the four control points. Nodes 1, 2, 3, and 4 correspond to the most offshore location, an intermediate propagation sector, the nearshore transformation zone, and the innermost point closest to the harbor entrance, respectively.

Frequency discretization remained weakly influential at all four nodes, indicating that the spectral representation in frequency had already converged before the waves entered the coastal transformation zone. The mesh sensitivity was also limited overall, but it became more visible at Node 3, where the larger RMSE, Bias, and P95 values suggest that bathymetric representation mainly affects the local shoaling and refraction pattern once the wave field reaches shallow water. By contrast, directional discretization produced a much stronger and spatially structured response from Nodes 1 to 3, indicating that angular sampling of the spectrum controls how the incident wave energy is redistributed during refraction and directional spreading. The strongest sensitivities were associated with the source terms. Depth-induced breaking primarily controlled the wave conditions at Node 4, where suppressing breaking produced a very large positive bias and extreme P95, evidencing unrealistically deep wave penetration toward the harbor. On the other hand, overly dissipative breaking formulations removed too much energy before it reached the inner domain. Triad interactions exhibited a different signature, with their strongest effect at Node 3 compared to Node 4. Hence, triads primarily reorganize wave energy within the shoaling zone, increasing or decreasing the local wave height before the final dissipation stage. The contrast between the strong response of Node 3 to triad-containing combinations and the dominant sensitivity of Node 4 to breaking formulations showed that the nearshore wave field was governed by a two-step process: first, a redistribution of energy within shallow water, and second, a depth-limited filtering that controls the wave energy ultimately able to propagate toward the harbor entrance. Overall, once a reasonable numerical resolution is achieved, the dominant source of uncertainty is no longer numerical but physical, and is governed mainly by the balance between the triad-induced energy redistribution at Node 3 and depth-induced breaking controlling transmission toward Node 4.

5.1.5. Uncertainty Propagation

Fig. 20 compares the sensitivity of nearshore H_m0 to perturbations in the offshore wave height and mean wave direction. The upper panels show that a ±10% change in offshore H_s produces the largest response, with |ΔH_m0| reaching approximately 0.30–0.40 m over a wide part of the offshore and intermediate domains. This suggests that the model response is strongly controlled by the incident wave-energy level. By contrast, the lower panels show that a ±5° directional perturbation has a weaker but still clearly identifiable effect, with |ΔH_m0| generally remaining lower and being concentrated mainly in the transformation zone before wave breaking and around the harbor-approach area.

Fig. 20

Spatial response of significant wave height H_M0 to perturbations in the offshore wave height and incident wave direction

Perturbations in the offshore wave height directly modified the incident energy flux and propagated through the model domain as a first-order amplitude effect, accounting for the broad, relatively uniform response visible in the upper panels. In contrast, directional perturbations mainly acted by altering the angle of incidence, and hence the local balance of refraction, shoaling, and energy focusing as waves propagate across irregular bathymetry. Therefore, their effect was weaker offshore and became most visible where wave rays are progressively refracted and redistributed in shallow water. The reduced sensitivity observed in the innermost zone further suggests that, once the waves enter the depth-limited coastal sector, part of the directional imprint of the offshore forcing was damped by shallow-water dissipation and refraction, which rotated the wave field toward directions more consistent with the local bathymetric control. Overall, offshore wave-height uncertainty exerted a stronger, more spatially persistent influence on H_m0 than moderate directional perturbations, whose effects remained more localized within the nearshore transformation zone.

5.2. TOMAWAC-TELEMAC-2D

5.2.1. Currents

Fig. 21 shows the velocity fields under tidal forcing alone (top) and under combined tide and large-scale imposed currents (bottom) at low tide (frame 4) and high tide (frame 11). In the tide-only runs, the flow at low tide (frame 4) was essentially one-directional along the coast: a strong coast-parallel jet (>3.0 m/s) formed in the southern nearshore and runs northward toward the harbor. This jet was mostly set by wave-induced longshore currents (north-coming waves) that act in the same direction as the tidal flow in that sector, so the two forcings reinforce each other. At high tide (frame 11), the water depth increased; the radiation-stress gradients weakened; the velocities dropped below ~2 m/s, and the pattern became more diffuse, with the first signs of a wider circulation cell offshore. The picture changed when regional currents were imposed at the northern boundary (tide + current cases) because the wave-driven littoral current interacted with the imposed large-scale flow. At low tide, the nearshore jet was still present, but it was wider and diverted offshore, where it interacted with the incoming background current. This interaction was clearer at high tide; the combined tide + wave-driven current + imposed regional currents produced well-defined recirculation cells over the inner shelf and lowered the velocity peaks compared with low tide. Overall, the figure showed that the tidal phase controls the strength (low tide implies faster currents); wave-induced longshore currents controlled the nearshore branch, and the imposed regional flow mainly reshaped the circulation by creating or strengthening the recirculation zones.

Fig. 21

Spatial distribution of the velocity magnitude and vectors under tidal and tidal-current scenarios

Therefore, the next step was to examine how these modifications in current dynamics influence the wave field itself.

5.2.2. Wave height

Fig. 22 shows the spatial pattern of ΔH_m0, relative to the reference run, for two forcings (tide only vs. tide + background current) at low tide (frame 4) and high tide (frame 11). Fig. 18 shows that the water depth controls wave propagation. At low tide, the water was shallow, so waves interacted strongly with the bottom and with the nearshore current. The wave-driven longshore current and the tidal current ran in the same direction along the southern coast, which pulled wave energy along the shore and toward the harbor entrance. Hence, ΔH_m0 stayed close to the coast and entered the inlet. The water was deeper at high tide. The wave–bottom interaction was weaker, and the tidal flow was more diffuse (with recirculation), so it no longer guided waves along the coast. The wave energy then remained offshore on the inner shelf, despite ΔH_m being larger. The behavior was the same when background (regional) currents were added. Two-way coupling between TELEMAC-2D and TOMAWAC was needed because these depth, current, and wave effects act together: TOMAWAC must receive the time-varying water levels and currents to propagate and refract waves correctly, and TELEMAC-2D must receive the wave radiation stresses to reproduce wave-driven currents and setup. Without full coupling, nearshore energy would be underestimated, and the focusing area near the harbor would be mislocated.

Fig. 22

Spatial variation of significant wave height H_m0 during low and high tide conditions with and without current interactions

Fig. 23 shows that the tidal and current forcings affect the four control points differently. At Node 1 (offshore), the tidal signal was clearly visible in both runs (tide-only and tide + current). Both rose well above the reference and follow the same tidal phasing. After approximately 35,000 s, however, the tide + current case became higher than the tide-only case and stays higher for the rest of the record. This indicates a different response when currents are present. The recirculating flow field maintained wave energy slightly longer than the tide alone, despite both runs being driven by the same tidal elevation. At Node 2 (intermediate depth), the same ordering was observed but with smaller differences. The reference stayed close to 5.0 m. The tide-only run rose to approximately 5.8–5.9 m, and the tide + current run was marginally higher over part of the tidal cycle. At Node 3, only the two tidal runs were close to each other (tide ≈ 5.7–5.8 m and tide + current ≈ 5.8–5.9 m), and both stayed above the reference case (≈ 5.1 m), showing that local breaking and dissipation largely mask the additional effect of the background current. At Node 4, the behavior was the opposite of that at the offshore points. The two tidal runs were slightly above the reference at low tide, but slightly below it at high tide. This matches the maps shown in Fig. 18. At high tide, the water was deeper; the wave–bottom interaction was weaker, and the tidal flow became more diffuse (with recirculation). Therefore, it no longer guided wave energy toward the inlet, and the waves stayed offshore on the inner shelf, even though ΔH_m0 was larger there. The same pattern persisted when the background (regional) currents were added.

Fig. 23

Temporal evolution of significant wave height H_m0 at selected control nodes under tidal forcing and combined tidal-current conditions

5.3. Defining Interface

TOMAWAC predicts energetic offshore seas in the fully activated configuration (triads, whitecapping, and quadruplets), with H_m0 peaks ≈ 6 m. Two complementary indicators were then mapped: the Ursell number U_r = H L²/ h³, which diagnoses nonlinearity in shallow water, and the Miche index H/h (with the breaking threshold γ = 0.78) (Fig. 24). This reference case was benchmarked against two hydrodynamic forcings: tide only and tide + regional current, each analyzed at low tide and high tide. The results revealed a coherent, phase-controlled migration of the thresholds. The U_r = 26 and H/h = γ isolines shifted seaward at low tide (shallower water increases nonlinearity and steepness) and landward at high tide (greater depth relaxes bottom control) because U_r ∝ h⁻³ and H/h ∝ h⁻¹. This tide-driven behavior was systematic and physically consistent, not a numerical artifact. Practically, the joint behavior of Ursell and Miche defines a mobile breaking envelope governed by the tidal level. Among all scenarios, the “tide + regional current, high tide” case generates the most seaward Miche limit (H/h = 0.78). Therefore, it defines the controlling envelope, and the ARTEMIS interface was placed seaward of this contour for stable harbor-scale nesting. This ensures that incoming seas are non-breaking and only weakly nonlinear, avoids contamination by partially broken waves, and preserves the correct phasing for diffraction inside the port.

Fig. 24

Definition of interface based on the criteria of Miche (1944) and Ursell (1953)

Accordingly, the black contour shown in Fig. 25 corresponds to the selected ARTEMIS computational domain. This boundary was deliberately chosen upstream of the nearshore breaking and nonlinear triad interactions (Q_tr), which TOMAWAC cannot fully resolve in shallow waters. The model avoided spurious breaking inside the spectral wave simulations by positioning the ARTEMIS domain seaward of these phenomena, ensuring that offshore processes are captured by TOMAWAC, while nearshore effects such as diffraction, reflection, and local transformations were explicitly handled by ARTEMIS. This choice guaranteed robust numerical coupling and ensured that the simulated agitation within the estuary and port areas was physically consistent and numerically stable.

Fig. 25

Nested meshes: global (TOMAWAC) and local refinement (ARTEMIS) at the harbor

5.4. ARTEMIS Results

Fig. 26 shows the spatial distribution of significant wave height extracted from the ARTEMIS simulations. A monodirectional random wave condition was imposed and extracted from TOMAWAC results, with an incident wave height of 1.2 m and a peak period of 14 s. The partial reflection coefficients of R = 0.45 for rubble-mound breakwaters and R = 0.90 for vertical quay walls were applied to reproduce the realistic reflective behavior within the port basin. The bathymetry provided by DPDPM (Fig. 3(b)) showed that the main channel leading to the estuary was largely filled with sediment, which reduced the local depths and promoted additional wave energy dissipation. The results revealed a clear reduction in wave height at the port entrance, driven by the combined effects of the port layout, sediment deposits, and variable bathymetry. As the waves propagated inland, variations in depth enhanced dissipation, while refraction, diffraction, and bottom friction contributed to the progressive attenuation of wave energy. Offshore high-wave zones (red areas) gradually transitioned to moderate and low values (green and blue areas) inside the port. This pattern confirmed the roles of bathymetry and harbor structures in attenuating wave energy and showed how sediment accumulation within the access channel acts as a natural buffer, reducing the wave impact within the harbor.

Fig. 26

Distributions of wave heights modeled with ARTEMIS inside the estuary and the Larache port

6. Discussion

6.5. Computational Cost

Table 3 provides a summary of the computational cost of the main model configurations used in this study. For the standalone TOMAWAC simulation, 5800 s of simulated time (approximately 1.61 h) required 42 min 54 s of wall-clock time on 32 CPU cores, corresponding to about 0.44 h of computation per simulated hour. For the coupled TELEMAC-2D–TOMAWAC configuration, a fully transient 24 h simulation required 18 h 31 min on 32 CPU cores, corresponding to approximately 0.77 h of computation per simulated hour. In this coupled configuration, the hydrodynamic computation uses a time step of 1 s, and the coupling between TELEMAC-2D and TOMAWAC is applied at every time step. When normalized by simulated physical time, the coupled configuration is approximately 1.75 times more computationally expensive than standalone TOMAWAC, while remaining within the same order of magnitude.

Table 3

Summary of computational cost metrics for the different model configurations used in this study

The ARTEMIS runtime must be interpreted differently. In contrast to TOMAWAC and TELEMAC-2D–TOMAWAC, this study did not use ARTEMIS as a transient model over a prescribed physical duration. Instead, it calculates a steady-state, phase-resolving solution for a given incident wave condition and directly provides the final spatial distribution of harbor agitation. Hence, its runtime of 8 min 24 s represents a convergence time for a given wave case rather than a cost per simulated hour. Overall, these results show that the computational efficiency of the proposed framework relies on a scale-adapted distribution of modeling tasks: TOMAWAC for large-scale spectral propagation, TELEMAC-2D–TOMAWAC for non-stationary wave–current interaction, and ARTEMIS for detailed harbor-scale agitation. This organization is also beneficial for scalability because the most computationally intensive processes are only executed where they are physically necessary. In practice, this allows using standalone TOMAWAC and ARTEMIS for repeated scenario testing at moderate cost, while reserving the coupled TELEMAC-2D–TOMAWAC configuration for selected transient conditions where wave–current interactions and tidal variability must be explicitly represented. Therefore, the framework is well-suited to pre-operational studies and scenario-based forecasting on modern parallel computing platforms, while leaving room for further optimization through adapted coupling frequencies, time-step strategies, and reduced simulation windows.

7. Conclusions

This study developed and applied a comprehensive multi-scale numerical modeling framework to investigate wave agitation at Larache Port, located on the northern Atlantic coast of Morocco. The proposed methodology captures the full cascade of wave transformation processes, from offshore propagation under regional forcing to diffraction- and reflection-dominated behavior within the harbor basin, by combining the spectral wave model TOMAWAC, the hydrodynamic model TELEMAC-2D, and the phase-resolving model ARTEMIS. The sensitivity analysis showed that mesh resolution and frequency discretization exert only a limited influence on the predicted nearshore wave field once a reasonable level of convergence is achieved, whereas directional discretization remains a more influential numerical parameter. By contrast, nonlinear triad interactions and depth-induced breaking parameterizations emerge as the dominant controls on wave shoaling, energy redistribution, and dissipation in shallow-water conditions. In addition, the dynamic coupling with TELEMAC-2D highlighted the major influence of tides and currents on wave transformation, showing that hydrodynamic modulation alters the incident wave heights, modifies the breaking thresholds, and reorganizes the spatial structure of nearshore wave processes. The transition from TOMAWAC to ARTEMIS was established using Miche and Ursell criteria, ensuring physically consistent nesting between the spectral and phase-resolving domains. This approach allows restricting ARTEMIS computations to zones where phase-resolving effects such as diffraction, reflection, and resonance become dominant, while avoiding applications in areas still controlled by active wave breaking. The results further show a marked attenuation of wave energy within the estuary and harbor basin, confirming the protective role of the existing port configuration and emphasizing the importance of physically based domain decomposition in complex maritime settings.

Beyond the specific case of Larache, this study contributed three major advances to coastal and harbor wave modeling. First, it provided clear evidence that triad interactions and breaker parameterizations govern shallow-water wave dynamics, whereas whitecapping and quadruplet formulations remained of secondary importance in the nearshore zone considered. Second, it shows that wave conditions in harbor approaches cannot be fully understood without accounting for tide- and current-induced modulations, which have direct consequences for harbor agitation, operability, and engineering design. Third, it establishes a robust and reproducible nesting methodology between phase-averaged and phase-resolving models for estuarine and harbor environments. The proposed framework is not limited to Larache and is intended to be transferable to other harbors exposed to complex interactions between offshore swell, tidal forcing, nearshore transformation, and basin-scale agitation. Its structured configuration procedure, physically based chaining criteria, and modular character make it particularly relevant for application to other Moroccan ports and, more broadly, to Atlantic and estuarine environments where similar hydrodynamic constraints are encountered. The methodology also has strong potential for storm applications. Such a framework can provide valuable insight into extreme-event behavior, navigation risk, and infrastructure vulnerability under energetic wave conditions, when increased water levels, intensified currents, and stronger wave forcing may enhance harbor penetration and agitation. In this respect, it constitutes a useful basis for routine coastal and port studies and for resilience assessment and climate-adaptation strategies under more severe hydrometeorological conditions.

Future work should focus on confronting the modeling chain with field measurements of waves, currents, and water levels in order to assess its predictive capability. Additional applications to storm scenarios and to other Moroccan and regional ports would help test the robustness and generality of the methodology across a wider range of environmental conditions. In the longer term, integrating this multiscale framework into real-time or near-real-time forecasting systems could support harbor operations, early-warning procedures, and adaptive coastal management. Such developments would represent a significant step toward digital coastal management systems in which interoperable numerical models serve as research tools and as operational decision-support systems for safer, more resilient, and more sustainable port infrastructure.

Article information Continued