1. Introduction
Nonlinear internal waves are found in many parts of the world’s oceans, for example, [Reference Helfrich and Melville32, Reference Jackson, Da Silva and Jeans33]. The phenomenon is commonly driven by tides or cramming of the pycnocline by winds. Waves moving in the shallow coastal ocean are susceptible to instability near the bottom, e.g., [Reference Boegman and Stastna8], within the pycnocline, e.g., [Reference Lamb36], near the surface, e.g., [Reference Vrećica, Pizzo and Lenain54, Reference Zhang and Alford57], and may contribute to vertical fluxes, e.g., [Reference Stastna and Legare51]. Regarding internal wave-driven instability in the bottom boundary layer, see also, e.g., [Reference Carr, Davies and Shivaram13, Reference Ellevold and Grue19, Reference Xu and Stastna56, Reference Zulberti, Jones and Ivey58].
Occurrence of nonlinear internal waves has been extensively studied in selected areas such as the South China Sea [Reference Chang, Lien, Lamb and Diamessis14], the New Jersey Shelf [Reference Shroyer, Moum and Nash47], the Strait of Gibraltar [Reference Alpers, Brandt and Rubino4], the Portuguese Shelf [Reference Quaresma, Vitorino, Oliveira and da Silva43] and the Australian North West Shelf [Reference Zulberti, Jones and Ivey58]. In contrast, studies of internal waves in arctic waters are few. Tidally forced internal wave-induced mixing has been documented on the shelf north of Svalbard [Reference Fer, Koenig, Kozlov, Ostrowski, Rippeth, Padman, Bosse and Kolås21]. A statistical study suggests several hot spots of internal waves on the Norwegian Continental Shelf [Reference Dokken, Olsen, Wahl and Tantillo17]. The present study site is located north of the Arctic Circle, at 69 degrees north, 14 degrees east, on the eastern side of the Norwegian Sea. The location hosts places for deep water formation where effects due to internal waves may contribute. Internal wave-induced instability may potentially stimulate local primary and secondary production. In our case, internal wave-driven instability may support nutrients to a cold water coral reef (Hola) located downstream of the internal wave groups in consideration.
Tidal generation of upstream nonlinear internal waves at an obstacle occurs by the lee-wave mechanism (e.g., [Reference Apel, Holbrook, Liu and Tsai5, Reference Grue28, Reference Liu, Holbrook and Apel38, Reference Maxworthy39, Reference Ramp, Yang and Bahr45]). The stratified flow effects at the Knight Inlet Sill are relevant for the case studied here (e.g., [Reference Farmer and Armi20]). Using field measurements and numerical simulations, [Reference Cummins, Vagle, Armi and Farmer16] studied upstream waves during the ebb tide, at the landward side of the sill. The depression of the pycnocline was found to steepen at its leading edge, evolving into an undular bore, eventually forming a sequence of solitary-like waves behind the forward front. Above and at the lee side of the sill, at the seaward side, a stagnant wedge and unstable jet over the steep downslope side were formed. Kelvin–Helmholtz instability occurred near the shear interface. Klymak and Gregg [Reference Klymak and Gregg34] ascribed asymmetry of the lee-wave response during flood and ebb tide to the density difference between the seaward and landward sides of the sill. The sudden widening, split, or bifurcation of the isopycnals was ascribed to breaking or mixing, and potentially both. Lamb [Reference Lamb35] pointed to the effect of viscous boundary-layer separation on the large lee wave and formation of downslope jet beneath the lee wave. Stastna and Peltier [Reference Stastna and Peltier52] investigated the analogy between topographically driven lee waves at a sill and the downslope wind storms associated with breaking internal mountain waves in the atmosphere.
The local lee wave formation at our study site shares similarity to that in Knight Inlet, in the sense that both have large vertical excursion, and the breaking processes may be similar. However, the internal flow conditions differ, where the barotropic tidal flow in our case is approximately one-tenth of that in Knight Inlet. The pycnocline is also much deeper (ca. 25 m) compared to the shallow stratification in Knight Inlet (ca. 4–6 m). The dominant shelf current in our case also differs from the conditions in Knight Inlet.
The internal wave generation on the seasonal pycnocline considered here occurs by interplay between: i) the boundary current along the shelf/slope, ii) the tide of the diurnal and semidiurnal modes and iii) a particular double ridge-valley shelf topography. The shelf/slope current is the dominant contributor to the generation. The tidal components are competing, where the diurnal component dominates at the position in question. We point to a resonant condition at the generation site on the shallow shelf in Vesterålen, where the internal diurnal wavelength corresponds exactly to the topographical ridge separation length. The correspondence enhances the excursion of two internal tidal troughs, which locally become strongly nonlinear. A related situation occurs in Luzon Strait in the Pacific Ocean, where a resonance between the semidiurnal tide and a double ridge is taking place, generating very large internal waves, moving along the comparatively deeper stratification of the South China Sea, see [Reference Alford, Peacock and MacKinnon2, Reference Buijsman, Legg and Klymak11, Reference Buijsman, Legg, Alford, Farmer, MacKinnon, Nash, Park, Pickering and Simmons12].
Satellite images exhibiting internal wave groups are used to identify their position, group separation distance and individual wavelengths. An ocean general circulation model (OGCM) is used to calculate the shelf/slope current, tide and seasonal pycnocline. The particular OGCM resolves the mesoscale circulation and hydrography of the Norwegian Coast at a horizontal resolution of 800 m [Reference Albretsen, Sperrevik, Staalstrøm, Sandvik, Vikebø and Asplin1, Reference Asplin, Albretsen, Johnsen and Sandvik6]. The methodology to combine data of nonlinear internal waves from satellite and ocean models at the mesoscale is tested out here. In this perspective, the recent Surface Water Ocean Topography (SWOT) satellite comes with new opportunities for obtaining accurate nonlinear internal wave data world wide.
Main research questions of the study include: 1. What are the properties of the internal wave groups that systematically appear on satellite images of the area in focus, and how is the generation mechanism? 2. Depending on the degree of nonlinearity, waves may potentially break and induce mixing of the water column, however, may breaking thresholds such as local amplitude, Richardson number (relevant to shear-instability) or other parameters, be calculated from the data sources (satellite images and ocean model)? 3. The wave groups occur on the steady boundary current along the Norwegian shelf/slope. A coral reef Hola is located downstream of the waves. Are the waves strong enough to induce breaking and mixing, possibly supporting nutrients to the reef?
Following the Introduction, the study area, shelf topography, current along the shelf and the modes of the barotropic tide are described (§ 2). The ocean model and a three-layer model for how to obtain the linear internal long-wave speed are described (§ 3). The internal wave groups are identified from satellite images (§ 4). The internal tide, current on the shelf and the generation period are given (§ 5). The propagation properties of the wave groups are discussed (§ 6). Finally, the discussion and summary are given (§ 7).
2. Study area, current and tide
2.1. Location and shelf topography
The Norwegian shelf and shelf slope are described in the MAREANO project [Reference Buhl-Mortensen, Hodnesdal and Thorsnes10]. Our study area is found in Vesterålen, at 69 degrees north and 14 degrees east, on eastern boundary of the Norwegian Sea (Figure 1a). Directed north east/south west, the Vesterålen shelf has a narrow width of 30–60 km. Lofoten-Vesterålen is known for its productive ecosystem (ref kilo rapport). The local shelf topography is characterized by crosswise ridge–valley–ridge–valley system (Figure 1b). The cold water coral reef Hola is located in the deepest marine valley (depth 250 m). A current in north east direction defines upstream and downstream positions. The upstream ridge Staurgrunnen has a depth of 80 m at its peak, while the shelf further upstream is 120 m deep. Staurgrunnen is marked by a cross on the 64 km long section a–b indicated in the figure. The downstream ridge Vesterålsgrunnen extends 20 km along the shelf direction. The distance between the downstream edge of the north-eastern side of Vesterålsgrunnen and the downstream edge of the north-eastern side of Staugrunnen shallow, marked by the two small circles on section a–b, is 36 km. (The distance corresponds exactly to the separation between internal troughs IT1 and IT2 discussed in Section 5.1.)

Figure 1. a) Location. b) Shelf topography viewed from the offshore. The coral reef Hola is located in the main marine valley. A trench follows the shoreline. The Vesterålen archipelago is indicated in dark red colour and the shelf slope in dark blue. Depth profile along section b–a and its elongation (black solid curve). Depths of −80 m and −250 m indicated. The position of the shallow Staurgrunnen is marked by a cross on section a–b. The shelf slope has NE/SW direction. Section a–b is 64 km long. Distance between the two circles on section a–b is 36 km. Made from kartkatalog.geonorge.no, DMT50, gis/paraview, by Håkon Dreyer and Lars Willas Dreyer. c) Shelf/slope current (colour scale). Time-average between 18 August 12:00 and September 12:00, 2019. d) Ellipses of the barotropic tide. Diurnal mode (black), semi-diurnal (green), 6.2 hr ×10 (red), 3.1 hr ×10 (magenta).
2.2. Current along the shelf
The North Atlantic Current (NAC) is a topographically steered western boundary current transporting warm and saline Atlantic water along the Norwegian Coast and into the Barents Sea, see [Reference Skagseth, Furevik, Ingvaldsen, Loeng, Mork, Orvik, Ozhigin, Dickson, Meincke and Rhines48]. The mean flow of the NAC is 40 km wide and of a mean velocity of 0.3 m s−1. The Norwegian Coastal Current (NCC) transports fresher water from the Baltic Sea and river runoff along the coast, see [Reference Christensen, Sperrevik and Broström15]. The NAC and NCC are modelled by the OGCM Norkyst (Section 3.1). The velocity vectors averaged vertically and in time, between 18 August and 2 September 2019, are in the range 0–0.6 m s−1 in the area of investigation (Figure 1c).
2.3. Diurnal and semidiurnal modes of the barotropic tide
Tides in Vesterålen are dominantly composed of the diurnal and semidiurnal modes. The authors [Reference Gjevik and Straume24], and more recently, [Reference Børve, Isachsen and Nøst9], have resolved the components of the barotropic tide in model-domains covering the North Sea, the Norwegian–Greenland–Barents Seas and the Atlantic Ocean, with favourable comparison to measurements. The OGCM computes the ellipses of the barotropic tide in the area of investigation, including the diurnal, semidiurnal, 6.2 and 3.1 hr modes. The dominant diurnal tidal amplitude (sum of K1 and O1 components) at the location in focus is 8.3 cm s−1 (Figure 1d).
3. Methods
3.1. The ocean model
The OGCM Norkyst is used to calculate the seasonal pycnocline, baroclinic (and barotropic) tide, current along shelfslope and on the shelf (the NAC and coastal current). Results of the calculations are used to analyse the signatures of the internal wave groups observed in three satellite images; see Section 4 below.
Norkyst covers the entire Norwegian coast, from Skagerak in the south to the Barents Sea in the north. The model resolves
• tides (forced barotropic),
• pressure-gradient driven velocity (2D plus 3D) through hydrography and air pressure,
• wind-driven ocean circulation with hourly forcing,
and does not resolve
• variability of sub-mesoscale dynamics below 800 m horizontal resolution.
Note that Norkyst cannot resolve the generation or upstream propagation of the groups of short, nonlinear and dispersive internal waves in question.
Norkyst is an operational setup of the Regional Ocean Modelling System (ROMS) at MET Norway, resolving the mesoscale circulation and hydrography of the Norwegian Coast on a horizontal resolution of 800 m with 35 vertical layers [Reference Albretsen, Sperrevik, Staalstrøm, Sandvik, Vikebø and Asplin1, Reference Asplin, Albretsen, Johnsen and Sandvik6]. ROMS solves the Reynolds averaged, hydrostatic primitive equations using a bottom-following coordinate system with free surface [Reference Shchepetkin and McWilliams46]. Accounting for the strongly variable bathymetry of the various small seamounts on the shelf, as well as the shelf-slope, the physical conditions of the coastal areas of Norway, including the coastal current, tide and wave effects, such as the Kelvin waves are reproduced [Reference Röhrs, Sutherland, Jeans, Bedington, Sperrevik, Dagestad, Gusdal, Mauritzen, Dale and LaCasce44].
Norkyst receives tidal forcing from a global inverse barotropic model of ocean tides, TPXO7.26 [Reference Egbert and Erofeeva18], through the prescription of amplitude and phase for sea surface elevation and currents for the major eight primary harmonic constituents (M2, S2, N2, K1, K2, O1, P1, Q1) of diurnal and semidiurnal frequencies. Atmospheric forcing is provided by surface fields from AROME-MetCoOp [Reference Müller, Homleid, Ivarsson, Køltzow, Lindskog, Midtbø, Andrae, Aspelien, Berggren, Bjørge and Dahlgren41], MET Norway’s regional weather forecast model. Lateral boundary conditions for Norkyst are provided by the Topaz basin scale model [Reference Xie, Bertino, Counillon, Lisæter and Sakov55]. The seasonal model pycnocline is a consequence of atmospheric forcing, river runoff along the coast of Norway, and the Baltic Sea brackish water reaching the model domain through its southern lateral boundary conditions. Baltic Sea and river freshwater is largely propagating northwards along the coast as the Norwegian Coastal Current [Reference Christensen, Sperrevik and Broström15].
At its given horizontal resolution, Norkyst explicitly resolves the horizontal diffusion of momentum and tracers due to chaotic mesoscale circulations features, such as eddies and larger fronts. The sub-mesoscale domain must be described as merely eddy-permitting, e.g., Norkyst exhibits fronts at the scale of a few kilometers but not to a sufficient detail to account for full diffusion due to such mechanisms. This lack of diffusive processes resolved is filled by implicit numerical diffusion and explicit horizontal diffusion. The setup of this model employs a 4th-order Akima scheme for advection of tracers (i.e., temperature and salinity) and an explicit Laplacian horizontal exchange coefficient for tracers and momentum of 10 m
$^2\,$s−1. Horizontal momentum advection is calculated using a 3rd order upstream scheme. The explicit diffusion is imposed along geopotential surfaces, while momentum diffusion is applied along model grid surfaces. Vertical diffusion of momentum and tracers is parameterized using a 2.5 order generic length scale mixing scheme [Reference Umlauf and Burchard53].
3.2. Internal wave speed. Three-layer model
The purpose of this section is the calculation of the linear internal long-wave speed c 0 of the density stratification. The density field is calculated by the OGCM along the computational section a–b for day 2019-08-27 at 13:00 (Figure 2a). Timing is two days before satellite image of signatures of internal wave groups in question (Figure 3a). Pycnocline, basically located between 15 and 50 m depth, has a density jump
$\Delta \sigma_T\simeq 1$ kg m−3, approximately. Density as a function of depth is detailed along section a–b, at 13.5 deg (leftmost position), 13.8 deg (above Staurgrunnen), 14.4 and 14.6 deg (both above Vesterålsgrunnen). An average of the four density profiles is indicated by the thick black solid line in Figure 2b,c.

Figure 2. a) Pycnocline along section a–b on 2019-08-27 at 13:00 hr computed by the OGCM Norkyst. Longitude position in degrees east (horizontal axis), depth in m (vertical axis). b) Density at position 13.4 deg (black dashed dot), 13.8 deg (black dashed), 14.4 deg (blue dashed dot), 14.6 deg (blue dashed), average over all four positions (black solid line). Fitted three-layer model to the average density (red solid line) with
$(N^2_1,N^2_2,N^2_3)=(0.059,0.435,0.017)/1000$ s−2,
$(h_1,h_2,h_3)=(13.5, 21.2,115.3)$ m. c) Same but fitted three-layer model (red dashed line) with
$(N^2_1,N^2_2,N^2_3)=(0,0.435,0)/1000$ s−2,
$(h_1,h_2,h_3)=(13.5, 23,113.5)$ m.

Figure 3. SAR images of internal wave groups in Vesterålen. a) Sentinel-1 29 Aug 2019. Indication of a series of groups 1–7 and another series of groups 8–10. b) Sentinel-2 15 Aug 2017. Groups 1*–6* indicated. c) (overleaf) ERS-1 16 Aug 2000. Groups 1**–6** indicated. Indication of the computation section a–b (64 km long). Location of the shallow Staurgrunnen (cross).
The phase speed c 0 due to the averaged density function is calculated by a three-layer model as follows. A pycnocline of constant buoyancy frequency
$N^2_2$ and depth h 2 is defined. The pycnocline separates an upper layer 1 of depth h 1 and constant
$N_1^2$ from a lower layer 3 of depth h 3 and constant
$N_3^2$. The lower layer 3 assumes a constant average depth along section a–b of 150 m. The three-layer model fitted to the averaged density function in Figure 2b has
$(N_1^2,N_2^2,N_3^2)=(0.059/1000, 0.435/1000, 0.017/1000)$ s−2 and
$(h_1,h_2,h_3)=(13.5,21.2,115.3)$ m. The resulting model density function is indicated by the red solid line in the plot. Following [Reference Fructus and Grue23], their section 4.1, the eigenvalue problem of the three-layer fluid is solved in the case of long waves, giving,

where
$T_j=\hat K_j\cot(\hat K_j h_j)$,
$\hat K_j=N_j/c_0$, (
$j=1,2,3$). This obtains
$c_0=0.422$ m s−1. A variant of the three-layer representation accounts for the density jump of
$\sigma_T=1$ kg m−3 across the pycnocline with zero N 2 in layers 1 and 3. The value of
$N_2^2$ is obtained by the maximum buoyancy frequency of the averaged density profile, denoted by
$N_0^2$. The thickness of the model pycnocline h 2 is then defined by
$h_2=g'/N_0^2$ where the reduced acceleration of gravity is
$g'=g\Delta \rho/\rho\simeq 10^{-2}$ m s−2 (
$\Delta \rho=1$ kg m−3, g the acceleration of gravity). This variant has
$N_0^2=N_2^2=0.435/1000$ s−2 and
$h_2=23$ m. The upper and lower layers have
$N_1^2=N_3^2=0$ and depths
$h_1=13.5$ m and
$h_3=113.5$ m. The resulting density function is shown by the red dashed line in Figure 2c. Calculation of the wave speed obtains
$c_0=0.417$ m s−1 which differs by 1 per cent from the calculation using the more accurate three-layer representation of the density.
Alternatively, c 0 is calculated from the maximum buoyancy frequency averaged along section a–b. This obtains
$N_0^2=0.60\cdot 10^{-3}$ s−2 and
$h_2=g'/N_0^2=16$ m on 27.08.2019 at 13.00. The average depth of the
$26.6\sigma_T$ isoline becomes
$H_0=25$ and is assumed to describe the middle depth of the model pycnocline. The three-layer model with
$N_1=N_3=0$ obtains
$c_0=0.419$ m s−1. The value fits very well with the predictions above. The pycnocline along a-b computed by OGCM for conditions three days later, on 30.08.2019 (at 13:00), obtains a somewhat smaller maximum buoyancy frequency of
$N_0^2=0.38\cdot 10^{-3}$ s−2 and depth
$h_2=g'/N_0^2=26$ m, suggesting a thicker pycnocline. The average depth of the
$26.6\sigma_T$ isoline becomes
$H_0= 25$ m, still. The three-layer model gives
$c_0=0.410$ m s−1. We may conclude from this subsection that the stratification along the computation section a–b has
• a linear internal long wave speed of
$c_0\simeq 0.42$ m s−1 on 27.08.2019 and
$c_0\simeq 0.41$ m s−1 on 30.08.2019;
• an average depth of the
$26.6\sigma_T$ isoline of
$H_0\simeq 25$ m on both days.
In Section 5.1, we estimate the speed
$c_0=0.422$ ms−1 directly from the separation distance between the two main crests IT1 and IT2, spanning across section a–b, corresponding to the internal diurnal tidal mode. This prediction results from the OGCM, where velocity, density and pressure variables are resolved simultaneously over the entire computational shelf area and over a substantial time window. The prediction fitting very well to the direct computation from the local stratification suggests that the evaluation of c 0 is robust. Moreover, the effect of a potential background velocity shear, which is accounted for in the calculation of the troughs IT1 and IT2, seems not to be present in our case. The effect of the earth’s rotation on c 0 has been evaluated by [Reference Grimshaw, Helfrich and Johnson25] and contributes to a correction smaller than 0.01 per cent for the conditions in our case, see Appendix A.1.
4. Observations from satellite images
4.1. Groups of internal waves
Three satellite images show several internal wave groups as identified by synthetic aperture radar (SAR) technique [Reference Alpers3] (Figure 3). Groups of the strongest visibility, localized above the shallow shelf next to the shelf break, are investigated. Wave groups appear in series and are found at positions south west (upstream) of the computational section a–b. Local depth where they are found is 120 m, approximately. The position of Staurgrunnen shallow is marked by red cross in the images. The distance between the groups is found to be approximately constant in each image.
Images were obtained in the month of August, by Sentinel-1 on 29 August 2019, Sentinel-2 on 15 August 2017 and ERS-1 on 16 August 2000. The 2019-image (Figure 3a) shows two series of wave groups, where series one includes groups indicated by numbers 1–7, while series two includes groups indicated by numbers 8–10. The two series may originate from two different sources. We are mostly concerned with groups 1–7. Groups 8–10 are briefly discussed below (Section 4.3). Groups 2A,B–6 are almost perfectly aligned. They are separated by an average distance of 12 km. The position of group 7 is found at the shelf break where the depth and the counter current both increase. We note that groups 1A,B have a shorter separation distance compared to groups 2–6. A reason for this is not known.
A similar series of wave groups is identified on satellite image taken on 15 August 2017 (Figure 3b). The image also shows the common separation distance. Groups 3–6 have an average separation of 8.7 km while the distance between groups 1–3 is somewhat smaller.
Satellite image from 16 August 2000 shows rather curved group fronts (Figure 3c). The groups are found both along the shelf break and off the shelf. Some wave trains are found well to the ocean side of the image. An average separation distance between groups 3–6 is 11 km, approximately.
Positions of the observed wave groups relative to Staurgrunnen shallow are summarized in Table 1. Wave group positions are also plotted as function of the group number (Figure 4a). Distance from Staurgrunnen grows linearly with group number. Linear fits to the positions of groups 2–6 in individual cases show all over separation distance of 10.1±1.7 km. This is an average obtained from the three images. The individual cases show a gap distance of 11.8 km (2019-image), 8.5 km (2017-image) and 9.8 km (2000-image).

Figure 4. a) Distance from Staurgrunnen to group front vs. group number. 29.08.2019 (
$\bullet$) with fitted line (solid), 15.08.2017 (+) with fitted line (- -), 16.08.2000 (×) with fitted line (-.) b) Nonlinear excess propagation speed,
$c/c_0-1$, KdV-model (solid line), measurement by Carr et al. (symbol), fit (
$A/H=(\Delta c/c_0)/0.331+\gamma (\Delta c/c_0)^7$, dashed).
Table 1. Distance in km between Staurgrunnen and the front of groups 1–6 obtained from satellite images

4.2. Individual waves
Wave groups include several individual waves. Distance between successive crests/troughs (or distance between the associated maximum slopes) indicated in the satellite images is denoted by wavelength λ. The leading wavelength λ 1 spans
$550-1830\,$m (2019-image),
$300-1050\,$m (2017-image) and
$270-830\,$m (2000-image) (Table 2). The average of λ 1 over the groups becomes 970 m (2019-image), 760 m (2017-image) and 530 m (2000-image). The number of waves identified in the groups is in the range
$n\sim 3-11$ (2019-image),
$n\sim 2-8$ (2017-image) and
$n\sim 4-19$ (2000-image). The number of waves averaged over the groups becomes n = 7 (2019-image), n = 4 (2017-image) and n = 9 (2000-image).
Table 2. First wave length λ 1 in meters and number of waves (in parenthesis) of wave groups 1–7 extracted from satellite image on 29.08.2019 (row 2), 15.08.2017 (row 3), 16.08.2000 (row 4). aSecond wave length λ 2

Ratio between the first and the second wavelength, and the second and the third wavelength, within each group, averaged over all groups presented in Table 2 becomes 1.3, approximately, and implies that longer waves travel ahead of shorter waves.
Nondimensional wavelengths of the 2019-image, of
$\lambda_1/H_0\sim 22-73$, suggest that the upstream waves are up to 3.5 times longer than young waves, which have a typical separation length of
$20H_0$. The wavelength range may be compared to similar published observations of
$\lambda_1/H_0\sim 16-87$, e.g., [Reference Halpern31, Reference Maxworthy39, Reference Melville and Helfrich40], [Reference Cummins, Vagle, Armi and Farmer16, Reference Ostrovsky and Grue42, Reference Stanton and Ostrovsky49], see Appendix A.2.
4.3. Wave groups 8-10 in the 2019-image
Wave groups may also radiate from a source at the coast. Signatures found on the satellite image from 2019 are indicated by groups 8–10 in Figure 3a. We are somewhat uncertain if the signal marked by 9 really corresponds to an internal wave group. Distance between groups 8 and 10 is estimated to be 18 km, approximately. The two groups may be generated by the internal Kelvin wave propagating along the coast, where the locally highly variable depth and the narrow trench may magnify the excursions of the pycnocline.
4.4. Prominent wave groups in the satellite images
Wave groups 1–6 in the satellite images are the most prominent ones. These groups
• are located south west of the end point ‘a’ of section a–b;
• are separated by a common distance of 10.1±1.7 km, which is an average obtained from the three satellite images;
• have a gap distance of 11.8 km in the 2019-image;
• have waves that are rank ordered.
5. Calculations by ocean model
5.1. Internal tide
Internal tide and motion of the pycnocline are calculated by the OGCM. Isosurface
$26.6\sigma_T$ exhibits two main internal tidal troughs indicated by IT1 and IT2 (Figure 5a). Troughs extending across the shelf are straight, aligned and orthogonal to the current. Timing of the plot on 27.08.2019 at 13:00 corresponds to the maximum of the depression IT1, evaluated at the position on the computational section a–b, during simulation through several days of August 2019. The plot in Figure 2a shows the density stratification for the same timing along section a–b. Both plots illustrate the excursions of the pycnocline two days prior to the satellite image of the internal wave groups, taken on 29.08.2019 at 16:51 hr.

Figure 5. a) Depth (m in colour scale) of the 26.6σT isosurface on 27.08.2019 at 13:00 hr. Internal tidal troughs IT1–2 and computation section a–b indicated. b) The excursion of the 26.6σT isosurface along section a–b in a Hovmöller diagram for the period 10.08–20.09.2019. Distance in km from the end point ‘a’ os section a–b (vertical axis). c) Vertically averaged current speed along (
$U^*$) and across (
$V^*$) section a–b, evaluated at the average position of the trough IT1 on 27.08.2019, for the period 10.08–20.09.2019. The magnitude of c 0 (red dashed line).
The separation
$\lambda_{1,2}$ between troughs IT1 and IT2, between 36 and
$37\,$km, corresponds very well to the distance of 36 km between the north-eastern side of Vesterålsgrunnen and the north-eastern edge of Staugrunnen shallow (Figure 1b). Resonant interaction between the coastal flow processes and the geometrical features of the double ridge-valley topography is evident, where the local amplitude of the involved modes of motion are enhanced. We point to the fact that the troughs IT1 and IT2, found to occur near each of the marked circles in Figure 1b, are both large and strongly nonlinear, where the pycnocline greatly exceeds its average level. Period
$P_{1,2}$ is connected to distance
$\lambda_{1,2}$ by the internal wave speed. We may assume that the main period is the diurnal one, in accordance with the dominant barotropic tidal mode at the location. This suggests that
$P_{1,2}=24\,$hr obtaining

This speed corresponds very well to the calculations of c 0 of the stratification of section a–b presented above (Section 3.2).
The dominant diurnal period during days 25–28 August is further documented by calculations of the internal tide near end-point ‘a’ of section a–b. We may compare plots b), h) and n) in Figure 6. They illustrate that positions of IT1 at 13:00 hr on the days 25, 26 and 27 August are similar. A second trough at a downstream distance of 11 km from IT1 has emerged in plot n) on 27 August at 13:00. The second trough is also seen in the sectional plot (Figure 2a). We may compare the sequence of plots at a later time 17:00. Plots i) and o) on 26–27 August show that IT1 becomes oriented south-north at this timing, where both plots are rather similar. The plots j) and p) at a later time 21:00 on 26–27 Aug are also rather similar. The calculations presented in the figure illustrate the dominance of the 24 hr mode of the internal tide near the position of end point ‘a’ of section a–b during the days 26–28 August.

Figure 6. The motion of the internal tidal trough IT1 near the end-point ‘a’ of the computational section a–b (indicated in plot a) through 25 August 09:00 to 28 August 05:00, 2019. Depth (m in colour scale) of the 26.6σT.
The IT1 computed by the OGCM never moves beyond Staurgrunnen shallow or end position ‘a’ of the computational section, however. The depth of IT1 (in colour scale) at the intersection with a–b is presented in the Hovmöller plot (Figure 5b), where the computed 24 hr oscillation period is pronounced during the days 26–28 August. The response is moderate earlier and later than the days 26–28 August. The computational result is in contrast to the observation in the satellite image on 29.08.2019 (Figure 3a). The signatures of the internal wave groups, exhibiting a constant gap distance between groups 2–6, suggest that the favourable conditions for generation at the 24 hr period occurred not only during 26–28 August, but also through the entire period 22–28 August. The discrepancy between observation and computation suggests that the OGCM has a slow response in representing the internal tidal component at the 24 hr mode, at the particular shelf topography in Vesterålen. The model shows a significant, large negative excursion at the 24 hr mode, only during parts of the observation on the satellite image, where constant separation (24 hr) between wave groups is found during several days.
5.2. Current on the shelf
Current is represented by a vertical average with components
$U^*$ and
$V^*$ along and across section a–b, respectively. Figure 5c shows computed
$(U^*,V^*)$ at the average position of IT1. Computed current at the end point ‘a’ and at the position of IT1 are found to be the same. The lateral velocity
$V^*$ is approximately zero during the actual days in August 2019. Time-average of
$U^*$ is approximately 0.42 m s−1 during 26–27 August, approximately 0.38 m s−1 during 28-29 August and approximately
$0.40\pm 0.03\,$m s−1 during 26–29 August. The computation suggests that the series of wave groups in Figure 3a is formed in near-critical condition with
$Fr=U^*/c_0$ close to unity.
5.3. Wave properties of IT1
Trough IT1 is dominant in the area of investigation. The sectional plot in Figure 2a shows a length of the depressed 26.2σT isoline in the along shelf direction, of approximately 3.6 km. The 26.6σT isopycnal drops from 17 m at upstream position to 70 m at the trough, corresponding to net excursion of 53 m.
5.3.1. Field measurements
Field measurements using chains of conductivity-temperature-density (CTD) sensors covering the water column have been carried out at Staurgrunnen shallow during July–September 2021 and 2022. The recordings, currently being analysed (unpublished), show some main trends. These are important to this paper and are presented here in brief:
• At downstream, north east position, a main depression is measured. Corresponding to the calculated IT1 of the internal tide, this has a wavelength of 3.6 km, obtained as an average over several events of the field measurements, corresponding very well to the computed wavelength of IT1. The measurements exhibit systematic excursions of the pycnocline between 25 and 50 m.
• The measured depression is found to propagate upstream in south-westerly direction. Elapsed time between two measurement positions obtains the local propagation speed relative to the ground. Adding the measured current velocity, the true propagation speed during several events is close to or somewhat smaller than the linear long-wave speed c 0.
• Measurement at Staurgrunnen peak exhibits a steep front and fission into a sequence of shorter waves of period between 20 and 40 minutes. The short waves are ordered according to the local excursion.
• A main wave event occurs each day corresponding to 24 hr period.
6. Propagation properties
6.1. Space-time relation of group fronts
Results presented above justify that wave groups observed in the satellite image from 2019 are generated at the diurnal period. (Note: the events from 2017 or 2000 cannot be simulated by the OGCM Norkyst.) Separation distance between the group fronts found to be 11.8 km in the 2019 case, combined with a time separation of 24 hr, makes a propagation speed relative to ground of
$\Delta c=0.137$ m s−1. Accounting for the effect of the counter current, of
$U^*=0.40 \pm 0.03$ m s−1 through days 26–29 August, the true frontal speed becomes

Wave groups locate upstream of end point ‘a’ of the computational section, where the local water depth is H = 120 m. This reduces the layer depth ratio to
$H_0/H\simeq 0.21$ and the reference speed to
$c_0=0.41$ m s−1. The dimensionless propagation speed becomes

The potential effect of the earth’s rotation on the wave group speed has been evaluated in Appendix A.1. The contribution is negative, of 0.018 times c 0, slightly increasing the nonlinear speed. However, the correction is within uncertainty bounds and disregarded.
6.2. Speed-amplitude relation of internal solitary waves
6.2.1. Weakly nonlinear theory
Weakly nonlinear Korteweg–de Vries (KdV) theory may estimate the nonlinear-dispersive propagation properties of the wave groups. Solutions of one-directional propagation are expressed by the stream function
$\psi(x,y,t)=a(x,t)\phi(y)$ where
$\phi(y)$ denotes eigenfunction and (x, y) horizontal and vertical coordinates. The amplitude function
$a(x,t)$ satisfies (e.g., [Reference Benney7, Reference Grue27, Reference Stastna and Lamb50])

where
$\epsilon=A/H_0$ (A amplitude) and
$\mu=\Lambda^2/H_0^2$ (Λ wavelength parameter). The nonlinear and dispersive coefficients read
$\alpha_0=(3c_0/2)\int_{-H}^0\rho \phi_y^3dy/\int_{-H}^0\rho \phi_y^2dy$ and
$\beta_0=(c_0/2)\int_{-H}^0\rho \phi^2dy/\int_{-H}^0\rho \phi_y^2dy$, respectively. Eigenfunction corresponding to waves of mode one satisfies
$\phi(y)=B_1 y/(-h_1)$ in layer 1,
$\phi=B_3(y+H)/h_3$ in layer 3 and
$\phi(y)=B_2^c \cos \hat K_2(y+H_0)+B_2^s \sin \hat K_2(y+H_0)$ within the model pycnocline (layer 2), using the three-layer formulation. Coefficients are determined by continuity of ϕ and ϕy at the upper and lower boundaries of the pycnocline, and
$\hat K_2$ is defined below equation (3.1). The ϕ is normalized by unity. Coefficients due to the three-layer model-pycnocline in Figure 2c (with
$N_2^2=0.435/1000$ s−2,
$N_1^2=N_3^2=0$,
$(h_1,h_2,h_3)=(13.5,23,83.5)$ m) are compared to two-layer (interfacial) case where
$\alpha_0H_0/c_0=-(3/2)(1-H_0/(H-H_0))$ and
$\beta_0/(c_0 H_0^2)=(1/2) (H/H_0-1)$ (Table 3). KdV-soliton of shape
$-A \mbox{sech}^2[(x-ct)/\Lambda]$ obtains dimensionless nonlinear excess speed of
$\Delta c^{KdV}/c_0=(1/3) A\,\alpha_0/c_0=0.331\, A/H_0$, where
$A/H_0\simeq 1$ gives speed corresponding to the observed speed of the group fronts.
Table 3. Nonlinear and dispersive coefficients of the KdV-equation. Nonlinear excess speed and wavelength of KdV-soliton

6.2.2. Nonlinear ISWs
Finite amplitude internal solitary waves (ISWs) have speed and wavelength exceeding predictions by KdV-theory. Experiments of ISWs by [Reference Carr, Davies and Shivaram13] have laboratory stratification similar to the field conditions where the wave groups in the 2019-satellite image are observed. The relative water depth of
$H_0/H=0.21$ is also the same in the laboratory and field conditions. This suggests that dimensionless wave speed and amplitude are also similar in the two situations. We refer to experimental laboratory runs 070207, 020307, 220207, 200207, 210207 in [Reference Carr, Davies and Shivaram13]. Measured speed of ISWs has an average of

(The value of c 0 was re-computed using the three-layer method.) The corresponding amplitude averaged over the five cases becomes

This indicates that amplitude of the strongly nonlinear ISW is somewhat higher than suggested by KdV-model (when speed is the same). Relation between amplitude and wave speed is illustrated in Figure 4b. The amplitude of ISWs is bounded by
$A+H_0=(1/2) H$ and implies that
$A/H_0$ is limited by 1.4 in the actual condition, where the field waves have amplitude close to the maximum.
7. Discussion and summary
7.1. Potential instability of the wave groups and IT1
The laboratory measurements by [Reference Carr, Davies and Shivaram13] exhibit instability and vortex formation in the bottom boundary layer, below the back shoulder of the ISWs. Regarding this kind of instability, see, also, e.g., [Reference Boegman and Stastna8], [Reference Ellevold and Grue19, Reference Xu and Stastna56, Reference Zulberti, Jones and Ivey58]. The dimensionless frontal speed of the wave groups in the 2019 field case is matching the
$c/c_0$ of the ISWs measured by Carr et al. The layer depths of the density stratification are also similar with
$H_0/H=0.21$ in the two cases. All of the runs in Carr et al. referred to above, carried out at a laminar stratified Reynolds number of
$Re_w=c_0H/\nu=6\cdot 10^4$ (ν kinematic viscosity), exhibit instability and vortex formation in the BBL below the wave. The threshold amplitude for instability in Carr et al. is
$A/H_0=1.06$, corresponding to 26.6 m in the field case. Simulations extending the instability boundaries by Carr et al. suggest that the threshold amplitude reduces when the scale increases [Reference Ellevold and Grue19]. The field condition of the 2019 case of
$Re_w=5\cdot 10^7$ suggests the threshold amplitude of
$A/H_0=0.24$ (6 m) for ISW-driven instability in the BBL to occur, using the laminar predictions by Ellevold and Grue. The similar turbulent instability process at the field scale may alter the stability threshold. However, the process does occur due to ISWs at the field scale [Reference Zulberti, Jones and Ivey58], as also found in the field measurements of IT1 (unpublished).
Convective breaking near the surface due to ISWs occurs when the amplitude is strong enough (e.g., [Reference Grue, Jensen, Rusås and Sveen30], and also in unpublished field measurements at the position of IT1). The strength of the wave groups in the 2019 case is in match with ISWs that exhibit this kind of breaking of the stratified surface layer. Regarding potential Kelvin–Helmholtz instability of the wave groups, we note that the laboratory ISWs in [Reference Carr, Davies and Shivaram13], and the calculations in [Reference Ellevold and Grue19] did not exhibit shear instability of the pycnocline similar to that in the field condition. Shear instability of the wave groups in the 2019 case is unlikely or may occur only occasionally, depending on the thickness of the pycnocline.
The trough IT1, corresponding to the generating source of the wave groups, as calculated in Figure 2a, is illustrated by the isolines of the local depression. Potential shear instability of the pycnocline may be estimated by the local Richardson number (Ri) and is evaluated by parameters along the isolines assuming ideal theory. The authors [Reference Fructus and Grue23] and [Reference Fructus, Carr, Grue, Jensen and Davies22] derived an exact relation between local Ri, upstream buoyancy frequency (N 2), wave speed (c), local horizontal velocity (u) and local excursion (A), obtaining,

where conservation of mass above the isoline of initial level
$\tilde h_0$ estimates
$1-u/c\simeq \tilde h_0/(\tilde h_0+A)$, and parameters in the evaluation are:
$N^2=0.8/1000$ s−2, A = 53 m,
$\tilde h_0=17$ m,
$c\simeq c_0$, relevant to the field condition. The low Ri-value suggests that depression IT1 is susceptible to Kelvin–Helmholtz instability within the pycnocline. The local water depth varies between
$H\sim 115-140$ m where IT1 is located. The large relative amplitude of
$A/H\sim 0.38-0.46$ suggests the occurrence of instability and vortex formation in the bottom boundary layer. The amplitude relative to the initial upstream level of
$A/\tilde h_0=53/17=3.1$ suggests the occurrence of convective breaking in the upper part of IT1, e.g., [Reference Grue, Jensen, Rusås and Sveen30, Reference Zhang and Alford57]. Indeed, field measurements of IT1 exhibit convective breaking near the surface, shear instability within the pycnocline and instability in the BBL (results under processing, unpublished).
The evaluations suggest that instability occurs both in the surface layer and BBL. At the bottom, resuspension of sediments or other particles located at the sea floor may occur. The instabilities occur where the wave groups and trough IT1 are located on the shelf, corresponding to a horizontal extension of approx. 60 km in the along shelf direction and approx. 25 km in the crosswise direction, see Figure 3a. The additional shear instability within the pycnocline of IT1 implies that both instability, mixing and vertical fluxes of mass, momentum and energy occurs throughout the entire water column, where IT1 is located. Mixed material is convected by the persistent current along the shelf/slope to downstream positions. This includes the location of the coral reef Hola.
7.2. Summary
Groups of internal waves on the shelf in Vesterålen, Norway, and the generating source, are analysed by data from a small set of satellite images and calculations using an ocean model, where the latter resolves the mesoscale circulation and hydrography of the Norwegian Coast.
The degree of nonlinearity of the groups is expressed by the dimensionless nonlinear excess propagation speed as determined by the analysis. The nonlinearity is pivotal to estimate the strength and potential breaking power of the waves. The observed field waves are connected to KdV-theory and laboratory models of ISWs of density stratification fitting to the field condition. In particular, a laboratory study [Reference Carr, Davies and Shivaram13], of stratification and wave nonlinearity similar to the field waves, exhibits ISW-driven instability and vortex formation in the bottom boundary layer, where instability occurs at even weaker nonlinearity when the scale is increased [Reference Ellevold and Grue19]. The instability likely occurs in the field condition, potentially inducing resuspension of particles in the bottom boundary layer, over the shelf area where the wave groups are located. Moreover, waves are strong enough for convective breaking of the stratified surface layer to occur, however, shear-instability of the pycnocline is unlikely.
The wave groups, propagating upstream, against the persistent current along the shelf/slope, oriented in north east direction, are generated by the approx. 3.6 km long, 50 m deep, and 20 km wide depression lee wave, driven by the internal tide at the particular double ridge-valley shelf topography. Evaluation of the breaking power of the depression lee wave suggests that shear-instability of the pycnocline, convective breaking in the surface layer and instability and vortex formation in the bottom boundary layer are all likely, suggesting that breaking and mixing occur throughout the entire water column (where unpublished field measurements confirm that the three different breaking mechanisms occur). Instability drives vertical fluxes of mass, momentum and energy. Mixed material is transported to downstream positions by the shelf/slope current. This includes the coral reef Hola located to the north east of the locally generated internal wave phenomenon which indeed may supply nutrients to the reef and should be further investigated by expertise within marine biology.
The generation occurs at the dominant diurnal period (24 hr) at the location, at the double ridge-valley topography extending across the narrow Vesterålen shelf, where the local wavelength of the diurnal internal tide fits exactly to the lee-slope to lee-slope distance of the topography. This makes two deep lee wave troughs of the internal tide extending across the shelf, oriented orthogonal to the shelf/slope current. The strongest trough, located upstream of the coral reef, generates the considered upstream wave groups. They systematically appear during the month of August, as observed on satellite images, where conditions of the pycnocline (the linear long wave speed), the internal tide and the interaction with the particular topography may be particularly favourable. The present analysis, supported by field measurements during July–August in 2021 and 2022 (unpublished), should encourage wider investigations of the phenomenon during extended periods of the year.
The OCCM in use predicts: a) the source(s) of internal wave groups that are emitted at the particular topography, b) the current on the shelf, c) the oscillation period of the wave generation, d) the characteristics of the seasonal pycnocline including the position of the isopycnals and the buoyancy frequency, whereby the internal long-wave speed and reference depth of the stratification are obtained (using three-layer model). Together, points b) and d) document that the wave groups in consideration are generated at the critical Froude number, in this case. From the satellite data, the number of wave groups is identified. This determines the duration of the generation and the spatial frequency of the wave groups.
Our methodology, to analyse data from satellite in combination with ocean model calculations at the mesoscale, for identification of nonlinear internal wave-generation and propagation, is general and may be applied also at other locations. The recently launched Surface Water and Ocean Topography (SWOT) satellite provides new opportunities for obtaining internal wave data of enhanced accuracy. Combined with ocean model calculations, properties of the generating source(s) and the local currents can be used to find the nonlinear propagation speed of the finite-amplitude internal waves. This is a measure that, so far, cannot be obtained from the satellite images. The degree of the nonlinearity of the waves and the consequences for instability may be estimated.
Data availability statement
The model data from Norkyst-800 is openly available from the Norwegian Meteorological Institute: Norkyst-800 ocean forecast archive, Norwegian Meteorological Institute [data set], https://thredds.met.no/thredds/fou-hi/fou-hi.html.
Acknowledgement
Acknowledgements We acknowledge with gratitude the comments and suggestions by one anonymous referee, which significantly helped improve the paper.
Author contributions
Conceptualization: J.G. Methodology: J.G; J.R. Data curation: J.G.; J.R. Data visualization: J.G.; J.R. Writing original draft: J.G. All authors approved the final submitted draft.
Funding statement
This research was supported by NFR300329, Ecopulse from the Research Council of Norway.
Ethical standards
The research meets all ethical guidelines, including adherence to the legal requirements of the study country.
Appendix A.
A.1. Effect of rotation
Grimshaw et al. [Reference Grimshaw, Helfrich and Johnson25] explored the effect of rotation on the linear internal long-wave speed. For the case of a two-layer fluid, which provides a relevant approximation to our case of a moderately thick pycnocline, they derived:
$c_0^2\simeq (g' h_1h_2)/[(h_1+h_2)(1+f^2(h_1+h_2)/3g')]$, where
$f=2\Omega \sin(69)=1.36\cdot 10^{-4}$ s−1 is the Coriolis parameter. This is given below their equation (A11). This shows that c 0 reduces because of the rotation. In our case, with
$h_1+h_2=150$ m and
$g'=10^{-2}$ ms−2, the modification becomes
$f^2(h_1+h_2)/3g'\simeq 0.000093$, which is very small compared to unity. The rotational effect on c 0 may be disregarded in our case.
The rotation effect on the nonlinear wave speed was evaluated for the Ostrovsky equation. The correction becomes
$c'=\gamma/k^2$ where
$\gamma=f^2/2c_0$ and k the wavenumber. The contribution to the group velocity is negative by the same amount,
$-\gamma/k^2$ ([Reference Grimshaw, Helfrich and Johnson25], their equation 8). In our case, the average length of the wave groups identified in the satellite image of 2019 is 3.6 km. A corresponding wavenumber is
$k=1.75\cdot 10^{-3}$ m−1. The negative contribution to the group velocity becomes
$-\gamma/k^2=-0.00719$ ms−1 and contributes to the third decimal place in the evaluation of the speed. The implication is that the observed group speed has a nonlinear contribution that is somewhat higher than the estimates given in the main section of the paper.
A.2. Wavelength
Shapes of exact internal solitary waves (ISWs) may be obtained by Dubreil-Jacotin-Long (DJL) theory or interfacial models (e.g., [Reference Grue, Jensen, Rusås and Sveen29, Reference Lamb and Wan37]). Ideal strongly nonlinear and fully nonlinear interfacial models were compared to field data of the COPE experiment [Reference Ostrovsky and Grue42, Reference Stanton and Ostrovsky49] for a wide range of layer depths and amplitudes. Reference calculations are relevant for interpretation of individual crests (troughs) of observed waves. However, the length of an ideal ISW generally differs from observed crest to crest separation, which evolves due to rank-ordered nonlinear and three-dimensional propagation properties.
The crest-to-crest separation distance is evaluated in existing publications. Short internal waves in the form of oscillations of the seasonal pycnocline at a sill in Massachusetts Bay had an observed period of 7.5 min and a calculated wave speed of 0.77 ms−1. Wavelength
$\lambda_1=344$ m and amplitude
$A_1\simeq 10$ m of the first wave and average pycnocline depth of
$H_0=16.7$ m were obtained [Reference Halpern31]. The experiments by [Reference Maxworthy39] exhibited a first wavelength of
$\lambda_1/H_0\simeq 19$ (experimental H 0 of 0.6 cm), and the three-dimensional calculations of Maxworthy’s experiments by [Reference Grue28] showed
$\lambda_1/H_0=21$ and amplitude of
$A_1/H_0=0.5$ after 70 % of the tidal period. The laboratory experiments by [Reference Melville and Helfrich40] and the nonlinear reproduction of those experiments by [Reference Grue, Friis, Palm and Rusås26] obtained
$\lambda_1/H_0=25$ and
$A_1/H_0=1.09$ (experiment) and 1.04 (computation). Cummins [Reference Cummins, Vagle, Armi and Farmer16] provided observations and calculations of the upstream ISWs at the sill crest in Knight Inlet in British Columbia. The observed leading wavelength was
$\lambda_1=62$ m at times 16:00:21 and 16:15:11, 83 m at 16:29:00 and 189 m at 18:27:24 (their figure 6), and the amplitude was
$A_1=9.6$ m at time 17:02:00 (their figure 7). Measured excursions of the thermocline in the COPE experiment exhibit
$\lambda_1/H_0=87$ and
$A_1/H_0=4$ [Reference Ostrovsky and Grue42, Reference Stanton and Ostrovsky49]. The results summarized in Table 4 suggest that the range of observed wavelengths in the satellite image from 2019 fit with data from existing publications.
Table 4. Wavelength (λ 1) and amplitude (A 1) of the first (or second) wave of a wave group divided by the average pycnocline depth H 0, as estimated from field data, laboratory tests or computation.
