In-Stream Tidal Energy Resources in Macrotidal Non-Cohesive Sediment Environments: Effect of Morphodynamic Changes at Two Bays in the Upper Gulf of California

Because of the need to diversify the renewable energy matrix and because hydrokinetic tidal energy technologies are mature, many in-stream tidal energy resource studies are available globally. Still, there are many questions regarding the effect of seabed changes on tidal energy resources. For coastal regions in particular, where the seabed is generally more mobile than in deep waters, bathymetric evolution could significantly affect tidal energy production. Here, two models are used to analyse the potential effect of natural morphodynamic change on tidal energy resources at two macro-tidal sandy bays, Adaír Bay and San Jorge Bay, in the Upper Gulf of California, Mexico. One of the models is (purely) hydrodynamic, and the other is a morphodynamic model (with hydrodynamic–morphodynamic coupling). The models are validated against tidal current observations obtained with acoustic Doppler current profilers in the region of interest, using three different error statistics, which showed good agreement between models and observations. The results also showed that the most significant bed changes and the largest renewable energy resources are located near the shore. Moreover, there was a good correlation between (a) regions with the most significant depth changes and (b) the areas where the difference in annual energy production with and without depth change was largest. Finally, a two-year simulation with the morphodynamic model permitted to analyse the seabed evolution of a zonal profile off Punta Choya, the headland between the two bays. This profile evolved towards a featureless equilibrium, as expected from the morphological classification for macro-tidal sandy environments under a dominant tidal forcing. However, most importantly, this natural evolution would not be detrimental to tidal energy exploitation at the site.


Introduction
The Gulf of California (GC), located in the northwest Mexican coast (See Figure 1), is the only marginal sea in the Eastern Pacific. Many studies (both numerical and observational) have described different oceanographic phenomena in the GC and their dynamics, for example, trapped coastal waves, coastal upwelling, gyre evolution and mixing, or tides and residual tidal currents, to name a few [1][2][3][4].
Some authors have suggested that tides in the GC may be a possible source of renewable energy [4,5]. In the upper GC, located north of the Great Island Region, and in particular close to the head, tides are dominated by the semi-diurnal M 2 principal component, with tidal ranges near the head varying between six and eight meters [3]. Spring-neap cycles are also evident throughout the GC, as observed in tidal gauge data obtained from http://predmar.cicese.mx/ (last accessed on 9 April 2021). Tidal ranges and spring-neap cycle characteristics both affect tidal energy production for tidal lagoon developments. Other studies have also shown that tidal currents are commercially viable energy resources in the Great Island Region and the upper GC [6].
The first analysis assessing the effects of morphology on renewable energy resources used the Infiernillo Channel as a case study [7]. The Infiernillo Channel is a strait between Tiburón Island and Mexico's mainland; it is around 40 km long and 8 km wide in its widest section. Magar [7] showed that variations of bathymetric conditions had remarkable effects on the monthly tidal power density (TPD) predicted by the models, specifically on the TPD magnitude and the locations of maximum TPD. This dependence of the TPD on the bathymetric input conditions may cause substantial technical impacts on Renewable Energy device deployments.
Coastal headlands are also potential in-stream tidal energy sites because of the local flow accelerations caused by the headland [8]. Therefore, one may envisage assessing the tidal energy resources around headlands within the GC. One such headland is Punta Choya, shown in Figure 1. However, there are no islands in front of it that constrain the flow, as in the Infiernillo Channel, the two headlands delimiting the entrance to Adaír Bay constrict the flow within the Bay. Near Punta Choya, the tidal range is more than six meters, and flow speeds between slack waters could be sufficiently large for in-stream tidal array developments. The energy demand of the region depends on commercial activities and the electricity and illumination needs of nearby towns, and therefore the power transmission lines are well developed [9]. While new renewable energy developments would require an increase in transmission line capacity, having the transmission line infrastructure nearby is advantageous.
Regional models may be used to evaluate either the effect of tidal turbines on the environment or the impact of the environment on the tidal turbines' efficiency. Chatzirodou et al. [10] and Robins et al. [11] have studied the former in the context of tidal sandbanks and formulated a level of "acceptable change". This research focuses on the latter by evaluating the effect of morphological change on energy production. Such an effect may be crucial for in-stream tidal plant planning if minor changes in morphology lead to significant changes in possible energy production. This study aims to assess tidal currents and morphological evolution under climatologically mild conditions in the Upper GC region and their effect on the tidal energy resources. We set up two models with high resolution in the area of interest. The models have the same open-boundary forcings at the mouth of the GC and same initial bathymetry, but one of them is only hydrodynamic, and the other is a morphodynamic model. The high-resolution part of the domain covers two macro-tidal bays on the Upper GC's northeastern side: Adaír Bay and San Jorge Bay (See Figure 1). The models are validated against acoustic Doppler current profiler (ADCP) observations. A climatologically mild year is assumed, so the models are forced only with tides at the open boundary. The flow velocity fields, the mean annual tidal power density (TPD) and the annual energy production (AEP) are determined for one simulation year. The morphodynamic model is run for a second year to assess how the seabed is evolving in the longer term.
The paper is organised as follows. In Section 2, the model equations and the initial and boundary conditions are described, together with the calibration parameters and statistical measures used for the calibration and validation of the models against the ADCP observations. In Section 3, the mean peak flow during spring tides, U STM , the annual average of the tidal power density (TPD), and some considerations on the annual energy production (AEP) are presented and discussed. Then, the implications of natural seabed changes for tidal energy resource exploitation in the area are analysed, particularly in the regions with the largest U STM , TPD, and AEP values. A two-year simulation is undertaken with the morphodynamic model to assess the seabed evolution over more extended periods. It is worth recalling that the model assumes mild weather conditions and that it is unlikely that this is valid over two years in a row. Still, the study allows us to understand the system's behaviour when tidal forcings prevail. Figure 1. Bathymetry of the study area. The coastline shown corresponds to the mean sea level. The water depths are higher than expected in some nearshore regions; these regions were masked in order to avoid inaccurate flow speed computations. This is not expected to have a significant effect on the computations in the rest of the domain. SJI1, SJI2, and AB are the locations of the acoustic Doppler current profilers (ADCPs) used for calibration and verification of the numerical model. The indices a, b, c and d indicate locations used for the in-depth energy resource analysis.

Hydrodynamics Formulation
The numerical simulations of the hydrodynamics for both GC_MD and GC_MD mor models were carried out with the hydrodynamics modules of Delft3D-FLOW. Delft3D-FLOW solves the incompressible Navier-Stokes equations under the shallow water and the Boussinesq approximations to simulate shallow water flows [12]. In 3D simulations, the default version of Delft3D uses curvilinear coordinates in the horizontal dimensions and σ terrain-following coordinates in the vertical dimension [13]. In this work, we implemented the two-dimensional, depth-averaged model with the default coordinate settings. The depth-averaged system of equations consists of two momentum conservation equations for the velocity components (U, V), a conservation equation for the hydrostatic pressure P, and the continuity (mass conservation) equation for the sea surface height ζ: ∂ζ ∂t where ζ is measured from the reference level corresponding to mean sea level (MSL) and d is the water depth from MSL, defined as positive downwards. The terms (F x , F y ) are the components of a vector force representing the unbalance of the horizontal Reynolds stresses; the Reynolds stresses are modelled using the eddy viscosity concept introduced by [14]. The forces (F x , F y ) are necessary because, although the model is capable of resolving the turbulent scales, usually the hydrodynamic grids are too coarse to resolve the fluctuations. Therefore the momentum equations are Reynolds-averaged, introducing socalled Reynolds stresses. The Reynolds stresses are related to the Reynolds-averaged flow quantities through a turbulence closure formulation. (M x , M y ) are the external forcing and source/sink terms due, for instance, to discharge or withdrawal of water, water-structure interactions, or wave stresses; D H is the depth-averaged horizontal turbulent eddy viscosity, and τ s − τ b represents the difference between the surface and the bottom shear stresses, respectively [15].

Sediment Transport and Morphodynamics Formulations
The sediment transport and morphodynamics are solved with an advection-diffusion equation and with Exner's equation, respectively. For the GC_MD mor model setup, we assumed that the seabed was mainly composed of non-cohesive sediments with a median sediment diameter D 50 of 200 µm, corresponding to a sandy seafloor; a specific density of 2650 kg m −3 ; a dry bed density of 1600 kg m −3 ; and an initial sediment layer thickness of 5 m. Note that in the northeastern side of the Upper GC, the sediment on the seafloor varies from medium to fine sand, as one can see from Figure 2 (modified from [16]). Figure 2 shows that in the high-resolution sub-domain (see Figure 3) D 50 is larger than 7.8 µm, and generally between 62.5 µm and 0.25 mm close to the coast, which corresponds to the range for non-cohesive sediments. This sediment composition agrees with the qualitative sediment sample observations we performed during the fieldwork campaigns of 2017 and 2018. Note that Carriquiry et al. [16] had very few sediment sampling stations in Adaír Bay and San Jorge Bay. However, since we have no riverine sources of sediment and no mangrove or saltmarsh fields with high organic content in Adaír Bay and San Jorge Bay, it is justifiable to assume that the sediments here are mostly non-cohesive. Thus, it is reasonable to take a uniform median grain diameter of the order of 0.2 mm throughout the study region as a first modelling approximation.  [16]). The white marks correspond to sediment sampling stations. The hashes correspond to the following sediment classes: Vertical, widely spaced (between the coast and contour 1) for 0.125 mm ≤ D 50 < 0.25 mm; horizontal, thinly spaced (between contours 1 and 2 from coast) for 62.5 µm ≤ D 50 < 0.125 mm; vertical, thinly spaced (between contours 2 and 3 from coast) for 31.25 µm ≤ D 50 < 62.5 µm; 45 • incline, thinly spaced (between contours 3 and 4 from coast) for 15.6 µm ≤ D 50 < 31.25 µm; minus 45 • incline (thinly spaced between contours 4 and 5, and very thinly spaced inside contour 4) for 7.8 µm ≤ D 50 < 15.6 µm; horizontal, very thinly spaced (inside contour 5) for D 50 ≤ 7.8 µm.
Van Rijn's model for sediment bedload and suspended load transport is used in the computations [17]. In this approach, the suspended sediment concentration C of the single sediment fraction [in kg m −3 ], satisfies in the depth-average formulation a two-dimensional advection-diffusion equation of the form: where (U, V) are the depth-averaged flow velocity components (m s −1 ), D H the sediment diffusivity (m 2 s −1 ), and S the sediment source/sink terms (kg m −3 s −1 ). The bottom boundary condition for non-cohesive sediments is a balance between the erosion (E) and deposition (D) flux terms, which in turn depend on the settling processes at a bottom reference bed level z = z b : with w s the (terminal) settling velocity (m/s). As there are neither estuaries nor other sources of sediment discharge in the area of interest, the source/sink terms in Equation (6) are directly linked to E − D, as defined in Equation (6). The model also includes a bedload transport formulation that is considered in the bathymetric evolution, the reader is referred to the Delft3D manual for further information [15]. As mentioned above, the same timestep is used for the momentum and the advection-diffusion equations. Note that model outputs, both at the observation points and in map format, are stored at 30-min intervals only due to data storage constraints.

Tidal Forcing Assumptions
Gross and Magar [18] ran a simulation that included the potential tidal forcing over the whole domain and the tidal forcing at the mouth of the GC. They concluded that the potential tidal forcing effect was negligible compared to that caused by the co-oscillation with the Pacific Ocean at the mouth. Therefore, here both the GC_MD and GC_MD mor models were forced at the mouth of the GC with 13 principal tidal harmonics (M2, S2, N2, K2, K1, O1, P1, Q1, MF, MM, M4, MS4, MN4) from the TPXO8.1 model, a Global Inverse Tidal Model [19,20]. The tides were the only forcing implemented in the model because we assumed climatologically mild conditions, and thus the tides dominate over other forcing mechanisms. This assumption is justified because the renewable energy resources generated by tidal currents in this region are significantly larger than those generated by wind-driven currents, even during a climatologically average year (see [6] for more details).

Bathymetry
The bathymetry chosen for the simulations was the General Bathymetric Chart of the Oceans (GEBCO2014). GEBCO2014 has a spatial resolution of approximately 900 m (see Figure 1). Because the GEBCO data was unreliable in the coastal lagoons, some regions were masked. Since this work focuses on the offshore regions around the coastal shoals of Adaír Bay and the areas in the two bays with water depths of 10 m or more, this masking does not affect the results presented here. Depth outliers in other parts of the domain were also removed. It is worth highlighting that although the Navy and other institutions provided some gridded bathymetry data obtained with in-situ surveys, their coverage was very small compared to the region of interest, and the resolution was also 900 m.

Sub-Domain Selection
The morphodynamic model GC_MD mor is equal to the hydrodynamic model GC_MD extended with sediment transport and morphological updating of the bathymetry. In both models, the computational domain is split into three sub-domains: A coarse domain d01 with a spatial resolution of 6480 m covering the western approach (where it connects to the Pacific Ocean) and the deep parts of the GC; an intermediate domain d02 with a spatial resolution of 2160 m in the coastal regions and the north part of the Gulf; and an inner domain d03 with a resolution of 240 m covering Adaír Bay and San Jorge Bay, located in the Northeastern part of the Gulf. The area covered by each of the three sub-domains is shown in Figure 3, with particular emphasis on the sub-domain with the highest resolution, namely d03, with grid cell size ∆x = 240 m. The grids are regular with square cells and are orientated roughly along the longitudinal axis of the GC. The communication between sub-domains occurs along the d01-d02 and the d02-d03 boundaries. All the sub-domains need to have the same parameter settings and timestep intervals, which increases computational costs. However, these costs can be compensated by carrying out the computations separately (using P-threading) at each sub-domain [15]. In fact, the domain-decomposition approach is more flexible, accurate, and generally faster than other nested or unstructured grid techniques.

Timestep Selection
The timestep used for all sub-domain computations is determined using the Courant-Friedrichs-Lewy (CFL) parameter, establishing a theoretical stability condition for subdomain d03. In general terms, the CFL parameter is defined as: where u is a characteristic speed, ∆t the model timestep, and ∆x = 240 m the grid size of the highest-resolution sub-domain, d03. For Delft3D-FLOW, in general, CFL < 10 is sufficient to achieve model stability. For this condition to hold in sub-domain d03, it was sufficient to take ∆t = 60 s. The same timestep is used in both the GC_MD and GC_MD mor models. Some authors use a morphodynamic acceleration factor, MORFAC [15], where ∆t mor is the morphological timestep, and ∆t the hydrodynamical timestep already discussed above. MORFAC generally is 3 to 10, or even larger, but the morphodynamic model still needs to be calibrated and tested with MORFAC = 1. Therefore, here we set ∆t mor equal to the ∆t satisfying the CFL condition in the sub-domain d03.

Calibration and Validation of the Models
The models were calibrated and validated using the tidal signal extracted from the model and from in-situ measurements. Because the hydrodynamic assumptions in both models are identical, we used GC_MD for the calibration and validation. The in-situ water levels and flow velocities were obtained from bottom-mounted acoustic Doppler current profilers (ADCPs) that we installed between 2017 and 2019 at three different locations (see in Figure 1). SJI1 and SJI2 were installed near San Jorge Island, and AB was installed west of Punta Choya in Adaír Bay. The average water depth at SJI1 was 16 m, and the average tidal range was 6.5 m. The average water depth at SJI2 was 22 m, and the average tidal range was 6.3 m. The average water depth at AB was 17 m, and the average tidal range was 7.3 m. The ADCPs were used either for calibration or validation according to their location and period of operation (the dates are in dd/mm/yyyy format): The ADCPs sampled at a frequency of one Hertz, but only 30-min or 60-min averages were stored in memory. Measurements were taken at one-meter intervals and 2.10 m upwards from the seafloor for SJI1 17 , and 1.89 m upwards for the other monitoring campaigns. The measurement cells closest to the sea surface were removed due to the signal's contamination by sound reflections from the water-air interface. The remaining cells' data was depth-averaged, and we computed the tidal signal using the t_tide program [21]. t_tide was also used to determine the tidal signal from the depth-averaged velocities predicted with the model. An example of the reconstructed tidal signals is shown in Figure 4, for the GC_MD mor model, and the data obtained with the ADCP in Adaír Bay. It is important to note that the tidal speed peak signal's timings obtained with the models and those in the ADCP observations are almost the same, as one can see from Figures 4 and 5. This is expected because of the predictability of tides. Therefore, only the velocity magnitudes are used for calibration and validation.
The model skill was assessed by determining the agreement between model X = U mod and ADCP data Y = U obs , for the reconstructed tidal signal. This agreement was analysed by comparing the Pearson correlation coefficient, ρ X,Y [-] (non-dimensional), defined as [22]: where cov is the covariance between the two time series, and σ X and σ Y the standard deviations of X and Y, respectively; the root mean square error, RMSE U , [m s −1 ]: and at times the absolute relative error, RE U [%], The time interval used for the model-measurement intercomparisons was 30 min at the SJI1 17 , SJI1 18 and SJI2 moorings, and 60 min at the AB mooring.

Calibration Results
We compared the tidal signal predicted by the model against the tidal signal observed at SJI1 17 over part of the sampling period to calibrate the model parameters. The parameters we varied were the reflection parameter (dimensionless), the horizontal viscosity (m 2 s −1 ), and the friction coefficient C f (m 1/2 s −1 ). The reflection coefficient was introduced by Stelling [12] in the boundary conditions to make the boundaries less reflective for disturbances with the eigenfrequency of the model area; in ocean and sea models, it is recommended to set it to zero ( [15], p. 203). The horizontal viscosity, on the other hand, is a parameter usually tuned according to the grid cell size, as it may affect model stability ( [15], p. 371). For models with cell sizes between 20 and 40 m, a horizontal viscosity of 0.5 m 2 s −1 works well. For grids with cell sizes in the order of 100 m, the default value of 1 m 2 s −1 is reasonable. Preliminary results showed that the model predictions were not very sensitive to slight variations of the reflection parameter or the horizontal viscosity [23], so we set them to the default values. Thus, we only tested in detail the model sensitivity to variations of the friction coefficient C f , based here on the Chézy formulation. Table 1 shows Corr U and RMSE U for seven values of C f (m 1/2 s −1 ), for simulations that were run from StartTime = 01/06/2017 to EndTime = 04/06/2017. The first three days were the model's spin-up time and were discarded in the computation of the errors. The final C f was taken as 45 m 1/2 s −1 , when RMSE U is smallest. Note that we expected that we would need to reduce the default value of C f to values around 45 m 1/2 s −1 , based on results by Baston et al. [24]. It was also for this value of C f when the plots of modelled and observed tidal signals, such as the one shown in Figure 4 for AB, looked best.

Validation Results
We used the ADCP moorings SJI1, SJI2 and AB, introduced in the previous section, to validate the results obtained with the models. Table 2 summarises the results of these validations for the yearly mean speed, U, obtained from the depth-averaged model predSictions (mod) and ADCP observations (adcp) of the flow speed. Overall, the correlation coefficient, the RE and the RMSE values obtained in the modelled-observed speed comparisons were acceptable, even though one may consider that a relative error of around 20% is at the limit of the acceptability threshold. SJI2 was the location and monitoring period with the smallest relative error, but a relatively large RMSE of 0.115 m/s. This large RMSE for the flow speeds causes significant errors for the TPD and the AEP validations at this location. We obtain the best agreement between the GC_MD model and the ADCP reconstructed tidal speed time series at the AB mooring, shown in Figure 4. However, in some tidal cycles, GC_MD seems to overestimate the flow speed.
To validate the morphodynamic model GC_MD mor , we used the depth-averaged reconstructed tidal ADCP measurements taken at the SJI1 17 mooring, and the corresponding depth-averaged, reconstructed tidal signal predicted by GC_MD mor . We found a good correlation of 0.88 (see Table 2) and an excellent relative error of 9% between model and measured data, but in Figure 5 we see that the GC_MD mor model overestimates the flow speeds at the ADCP moorings. Table 2 shows the value of the RMSE U , which is quite high, considering the mean values of the velocity field at this location, also shown in Table 2.

Results and Discussion
We ran a one-year simulation of the GC_MD and GC_MD mor with the calibrated parameters, to map and analyse the yearly averaged spring tidal speed maxima, U STM , the yearly averaged Tidal Power Density (TPD), and the Annual Energy Production (AEP). The instantaneous TPD, TPD i , is defined as: with U i being the instantaneous speed. Its mean is defined simply as and the AEP is defined as: Here N is the total number of hourly data in the simulated year, or the number of times when TPD i is above the threshold value of 50 W m −2 , depending on the variable used. This technical TPD threshold corresponds to a low-speed threshold of 0.46 m s −1 , typical of most tidal energy devices, as discussed in [6].
The mean spring tidal speed maxima, U STM , is computed by identifying the maximum flow speed in spring tides, which for the year of simulation gives 26 maxima, and then computing the mean of those values. Figure 6 shows the yearly mean of the spring tidal maxima flows, U STM , for the GC_MD (Subplot 1) and the GC_MD mor models (Subplot 2). The largest U STM values in San Jorge Bay were also on or near regions with bottom sand ridges, but the effect of the bottom bathymetry is more noticeable in Adaír Bay. The largest values of U STM were observed as well in the nearshore of Adaír Bay, where we have a series of shoals, with the site just North of Punta Choya (see Figure 1) being the most important region for tidal energy resources. It was expected that Punta Choya would play an essential role for tidal energy resources because it is a headland, and other researchers have observed that coastal flows tend to accelerate around headlands [8]. It is worth mentioning that these possible tidal energy exploitation sites had not been identified with low-resolution models [6], which highlights the limitations of those models for nearshore tidal energy characterisation. Figure 6 also shows the effect of coastal morphodynamics on the spring tidal maxima flows, for example, that the areas where the yearly averaged U STM is 1 m s −1 or above expand when the bottom is allowed to evolve under tidal action. This effect will be explored in more detail using the yearly tidal power density. Figure 7 shows the yearly mean TPD obtained for the GC_MD model (Subplot 1) and the morphodynamic model GC_MD mor model (Subplot 2) over the simulation year. We found areas in the nearshore region of Adaír Bay, where the maximum TPD values are above 200 W/m 2 , but in general below 400 W/m 2 . We observe TPD values that are significantly lower than those reported for the seven US sites analysed in the Tidal In-Stream Energy Conservation (TISEC) program [25]. The TISEC site with the lowest power density was Head Harbour Passage (New Brunswick), with a TPD of 940 W/m 2 . In contrast, the Minas Passage in Nova Scotia has a TPD of 4.5 kW/m 2 , i.e., more than 20 times larger than the largest TPD observed in our study site. For sites with lower TPD, the return on investment will take longer. Small scale also increases the costs by a factor of three [25], and this applies to pilot sites. Despite these drawbacks, this site could still be a useful test site for some devices. As mentioned in the previous section, the best tidal current energy resources are located off the headland separating Adaír Bay and San Jorge Bay, but other tidal current energy hot spots can also be observed along the eastern coast of Adaír Bay. These regions are, consistently, in the same locations of the domain in the GC_MD and the GC_MD mor models. However, in GC_MD mor , the regions with high TPD are spread over larger areas. We observe, in general, an increase in yearly averaged TPD throughout Adaír Bay in the  Table 3). However, initially the TPD at (a) is much larger than at (b), so even if percentually the yearly TPD does not increase as much in (a) as in (b), the tidal current energy resource is still better at (a). At location (c) in the nearshore of San Jorge Bay, in contrast, the tidal current energy resource decreases in GC_MD mor in comparison to GC_MD, by 16.9%. Finally, at location (d), %(TPD) = 91.26%, but initially the tidal current energy resources are the lowest here, and therefore it remains low in GC_MD mor . Please note that %(AEP) = %(TPD).  Figure 8 shows the AEP maps, with contours of 500 kWh m −2 delimiting high-energy regions in Adaír bay; these contours can be seen in the results for both the GC_MD (subplot 1a) and the GC_MD mor models (subplot 2a). When using the GC_MD mor model (subplot 2a), i.e., when the seabed evolves in time, these high-energy regions cover a much larger area. In particular, the region with maxima AEP located at the southeastern boundary of Adaír Bay extends further to the west in the model GC_MD mor , in comparison to GC_MD. In the plots 1b and 2b of Figure 8, we also observe this spreading of what we may call good sites for tidal energy exploitation, as measured by the spread of the 200 kWh m −2 contour, extending in San Jorge Bay in GC_MD mor , and an increase from 200 to 400 kWh m −2 in the locations of large AEP. However, in San Jorge Bay these regions are significantly smaller than in Adaír Bay, for both the GC_MD (plot 1a vs plot 1b in Figure 8) and the GC_MD mor models (plot 2a vs. plot 2b in Figure 8). We can identify the sites with the best renewable energy resources from the AEP maps. Most importantly, however, is the fact that the plots in Figure 8 show the impact of natural seabed evolution on areas with exploitable renewable energy resources. Indeed, within a year the size of such regions has more than doubled in Adaír Bay, a reassuring result for any would-be tidal energy developer interested in implementing a tidal energy demonstration plant in the region.

Technical Constraints
Some devices, such as Flumill [26] for deployments in waters as shallow as 10 m and up to 50 m deep, or Minesto [27] for deep water deployments, can operate at the low tidal stream speeds of around 1.0 m s −1 that are found in some locations in the Adaír and San Jorge Bays (as shown in Figure 6). Another technical constraint to consider is the cut-in speed. Most commercially viable technologies require cut-in tidal speeds of around 0.5 m s −1 before they start generating any power from the flow; this cut-in speed corresponds to TPD = 50 W m −2 . In Figures 9 and 10, we only consider TPD values of 50 W m −2 or above, defined as TPD = TPD >50 , and use this for the corresponding computations of AEP = AEP >50 in Figures 10. The mesh resolution of 240 m, combined with the 30-arcsecond resolution (equivalent to around 900 m) for the bathymetry, gives a better representation of the tidal energy resources in shallow regions in comparison to the global model used by [6]. However, we observe that the global model ( Figure 5 of [6]) and the GC_MD model (left subplot, or subplot 1, in Figure 9 in this paper) both show that the northern side of San Jorge Bay, and the nearshore areas in Adaír Bay, do not provide commercially viable conditions for in-stream tidal energy exploitation. Despite this, we can also see some areas with TPD >50 over 200 W m −2 in Adaír Bay, such as the nearshore area marked as (a) in Figures 9 and 10, or the area at the southern edge of Adaír Bay marked as (b). Unfortunately, at the observation point AB inside Adaír Bay, we obtained low TPD values and a correspondingly low AEP, but this does not affect the use of AB for model validation purposes. The area marked (c) in San Jorge Bay is the only area in this Bay with TPD >50 above 200 W m −2 . As can be seen on the left subplot of Figures 10, the area (a) has the largest AEP >50 of the four marked regions. In fact, at (a) the AEP obtained after one year with the model GC_MD is actually 1122.6 kWh m −2 , as shown in Table 3. The only place where such high values of AEP are observed is in the nearshore region of Adaír Bay.  (subplots 1a and 2a) and San Jorge Bay (plots 1b and 2b) obtained with the GC_MD (subplots 1a and 1b) and with the GC_MD mor models (subplots 2a and 2b).   Figures 9 and 10 show the predicted TPD >50 and AEP >50 for the GC_MD (on the left) and the GC_MD mor (on the right) models. Firstly, we observe that the best locations for in-stream tidal energy exploitation do not shift in space after one year when we have seabed changes (model GC_MD mor results on the right) in comparison to when we have a fixed bathymetry (model GC_MD results on the left). In Adaír Bay, the area covered by these locations increases in the GC_MD mor model, which is reassuring as the seabed in this region is naturally mobile, and therefore the GC_MD mor predictions are expected to be closer to reality. In contrast, in San Jorge Bay, although we see an increase in in-stream tidal energy resources predicted by GC_MD mor at the locations where the resources were good in GC_MD, the area they cover does not increase in comparison to the area they cover in the GC_MD simulation. While the resources remain marginally exploitable in both simulations, these analyses indicate that Adaír Bay is a better location than San Jorge Bay for a demonstration in-stream tidal energy site, and its suitability improves when the tides force the seabed to evolve, at least in the first year.

Effects of Morphological Evolution on Tidal Energy Resources
We consider now the percentage of time, %T, when the TPD is equal or above the cut-in TPD of 50 W m −2 , and evaluate the impact of morphodynamic evolution. Subplots 1a and 1b in Figure 11 show %T when TPD is above 50 W m −2 , for the GC_MD model, and subplots 2a and 2b in Figure 11 show %T when TPD is above 50 W m −2 , for the GC_MD mor model. The maps in subplots 1a and 2a show Adaír Bay, where we can identify regions where TPD is above 50 W m −2 for more than 40% of the time, for both for GC_MD and GC_MD mor . There are also vast regions in Adaír Bay where %T is close to 60%. Figure 11 shows results that are consistent with those in Figures 9 and 10, obtaining good agreement between areas with the largest TPD >50 and AEP >50 and areas with the largest %T. We also note that the region along and near the same longitude as that of location (b) (see Figure 9), is the region with the largest morphological changes near the coast, in particular close to the eastern shore of the bay, but this impact is positive from an in-stream renewable resource perspective. On the other hand, in San Jorge Bay we found that the percentage of time when TPD is above the cut-in value is generally below 40% (Subplots 1b and 2b in Figure 11), and in the GC_MD model (no morphodynamics), the regions satisfying this constraint are quite small (Figure 11, 1b). These results show that the optimal location for tidal energy resource exploitation is near the southeastern corner of Adaír Bay, west of the peninsular hook known locally as "Punta Choya". The results show that when the seabed evolves with time, the regions with high tidal energy potential remain centred in the same locations, but the area they cover spreads and becomes larger within a year. In these regions %T is between 40% and 60%, for both GC_MD and GC_MD mor simulations.
We now investigate in more detail the bathymetric changes over one year, and the changes they cause on tidal energy resources. Figure 12 shows the difference between final and initial bathymetry for GC_MD mor after one year of simulation, Figure 13 shows the difference in AEP predicted by GC_MD and GC_MD mor for that same year, and the subplots in Figure 14 show local depth change for that same year at the locations marked as (a), (b), (c), and (d) in Figure 12. We first describe and discuss the depth and AEP change pattern maps (Figures 12 and 13), after one year of simulation, over Adaír Bay and San Jorge Bay. In general, we observe a close correlation between regions with large depth changes ( Figure 12) and regions with large AEP changes (Figure 13), in particular in areas where the AEP predicted by GC_MD was already large. The largest depth changes are in the nearshore region of Adaír Bay (see Figure 12a), particularly off Punta Choya (marked on the map of Figure 13). Figure 11. Percentage time, %T, when TPD is above 50 W m −2 , for Adaír Bay (subplots 1a and 2a) and San Jorge Bay (plots 1b and 2b) obtained with the GC_MD (subplots 1a and 1b) and with the GC_MD mor models (subplots 2a and 2b).  Table 3 shows d i (m), the water depth for the GC_MD model, which remains unchanged; d f (m), the final water depth at these locations in GC_MD mor , with the initial water depth being the one reported for GC_MD; AEP (kWh m −2 ), the AEP obtained with the GC_MD model; AEP mor (kWh m −2 ), the AEP obtained with the GC_MD mor model; and the relative change in AEP relative to the purely hydrodynamic model.
In the nearshore of Adaír Bay, marked as (a), the change in water depth can be larger than one meter, as can be seen in Figure 12a. Since the initial water depth here is of 5.7 m, with a relative change in water depth is of 22.8%. This is within the range of 15 to 25% of relative bathymetric change that can be qualitatively observed in this region (see also Figure 14a). The bathymetry in the nearshore of Adaír Bay consists of a series of sandy shoals, which are flattened over time by the tides over the simulation year. This shoal flattening leads to an increase in water depth (yellow/light grey regions in the nearshore) at the crests of the shoals, and a decrease in water depth on the slopes of the shoals (blue/dark grey regions, adjacent to the yellow/light grey regions previously mentioned). That is, the tides tend to flatten the shoals in the nearshore region, but do not make the shoals migrate. However, there is a net transport of sediment towards the sea. This region is very shallow and perhaps some tidal fences or horizontal Flumill devices could be used for energy exploitation, however they may be exposed at low tide, so the tidal energy plant would have to be an exclusion zone to ship traffic. Implementing an exclusion zone could cause security problems and would reduce the temporal windows for energy generation. We observe an increase of 34.5% in AEP in GC_MD mor relative to GC_MD at location (a)-see Table 3. In this area we have relatively large changes in depth in one year, but a relative increase in AEP that is modest in comparison to the change in AEP at other areas of the study region.  The southern edge of Adaír Bay is the second area that we investigated, marked as (b) in Figure 12. In this region the relative change in water depth is significantly smaller than in the nearshore area of Adaír Bay, because, as we can see in Table 3 and in Figure 14b, the initial depth was of 11 m, and the change in water depth after one year was a bit larger than 3 cm. Therefore, the relative water depth change is really small, around 0.27%. Despite this very small relative water depth change, we observe an increase of 69.3% in AEP between GC_MD and GC_MD mor . The local change in water depth, being almost negligible at location (b), cannot alone explain such a large increase in AEP when the morphology is allowed to evolve. Therefore, the local change in AEP is likely driven by a local water speed increase due to bathymetric changes in nearby locations. It is important to note that location (b) is still quite shallow, even if it is not as shallow as location (a), and that this will also constrain the size of the devices that can be installed here.
The third area investigated, marked as (c) in Figure 12, is the nearshore area of San Jorge Bay. Here we have a small, curved sand ridge that initially increases in size due to sediment feedbacks from the nearshore, but after that the water depth increases as the ridge is flattened by the action of the tides. The initial increase in sand ridge size is of only 20 cm, and the ridge flattening induced a final bathymetry change from the initial depth of around 20 cm as well, after on year of simulation. As shown in Table 3, the initial water depth is 2.7 m, so here the relative change in bathymetry between the initial and final bathymetries is of 7.4%. Here the AEP in the GC_MD mor model is smaller than that predicted by the GC_MD model, corresponding to a decrease of 16.9% in AEP. Finally, it is important to note that this area is possibly too shallow for tidal energy exploitation, except for devices developed specifically for such environments.
Finally, we consider the least-shallow region at the southern edge of San Jorge Bay, marked as (d) in Figure 12. As shown in subplot d) in Figure 14, the initial water depth here is of 30.32 m, and the bathymetry only changed around 1 cm, so the relative change in water depth is only 0.03%. Despite this, we have an increase in AEP of around 202 kWh m −2 between the GC_MD mor and the GC_MD models, corresponding to an increase of 91.3% in AEP, between GC_MD mor and GC_MD. The AEP predicted by both models, however, is much smaller in this area than in the nearshore of Adaír Bay. As with location (b), it is likely that the observed difference in AEP at location (d), where the change in water depth is almost negligible, is due to a local increase in water speed that is driven by larger water depth variations at nearby locations within the domain.
From this analysis of changes in annual energy production at these two macro-tidal bays, based on a coupled hydrodynamic-morphodynamic model relative to a purely hydrodynamic model, we may infer the following. In both models the largest AEP is located near the shore and in the area with the largest tidal range of the region, these two conditions are satisfied in the northeastern side of Adaír Bay, and specifically in area just north of "Punta Choya" headland (location a). Other shallow areas also undergo bathymetric changes but they are not as important, and the AEP in these regions, although still considerable, may not be exploitable, due to the extreme shallowness in the area (location c).
Large changes in AEP between the two models are observed even in areas where the bathymetry changes very little, such as the deeper regions in the mouth of Adaír Bay (location b), or south of San Jorge Bay (location d). These large AEP changes are likely to be due to local velocity accelerations associated with bed changes at other locations. However the AEP in these regions were smaller in absolute terms than in the shallow regions, and even if relatively they increase more in the hydrodynamic-morphodynamic model, in comparison to the hydrodynamic model, than at other locations, they are still less attractive in AEP terms, than the shallow regions, for tidal current energy exploitation. This implies that those shallow regions are to be preferred under mild-weather conditions, and this choice is still best even when we take into account the change in AEP caused by morphological changes over a year.
Finally, we now assess the depth changes along a transect across the region with largest AEP changes, which corresponds to a transect off Punta Choya and passing through the area (b), as shown in Figure 13. Figure 15 shows depth changes as a function of distance to shore after one year and after two years of simulation. The region between 2 and 5 km from the shoreline is an erosion zone, while the region between 5 and 8 km from the shoreline is an accretion zone. We also observe a small accretion zone between 11 and 13 km from the shoreline. The figure shows there is a net seawards migration of sediment over time. The profile is likely to be converging towards an equilibrium profile under the action of the tidal forcing. The equilibrium profile has a slope that is much smaller than that for the equilibrium beach profiles defined by Dean [28] for wave-driven sandy beaches, in line with macrotidal environments such as this one [29]. However, the tendency towards an equilibrium is not observed at other locations further offshore, as Figure 14 demonstrates. Figure 15. Bathymetric profile along a transect off Punta Choya: Initial profile, and evolution after one year and after two years. We observe how the sandy shoal located between two and six kilometers from the shore gets flattened out by the tides, and the sand moves offhore; this means this part of the profile is tending towards a featureless (equilibrium) profile. Further offshore, the profile changes very little after two years of simulation.
In this paper, we have demonstrated the importance of including the hydrodynamicmorphodynamic coupling in numerical simulations, but also that bathymetric evolution is not detrimental, at most locations, to tidal current energy production. Finally, it is also important to note that when the devices are present, they extract energy from the system. In that case, other forcings' effects may become more important than we have assumed in this work and should not be neglected.

Conclusions
In this work, we developed two regional domain-decomposition models with three sub-domains, covering the Gulf of California (GC), Mexico, and with a high resolution sub-domain on the northeastern side of the Upper GC, in order to assess the likely effects of morphological evolution on tidal energy resources at San Jorge Bay and Adaír Bay, two macrotidal sandy bays located in this region. The three sub-domains, gridded models with one-way coupling (from coarse grid to fine grid only), had regular grids with cell sizes of 6480 m, 2140 m, and 240 m, respectively. The GC_MD and the GC_MD mor models were set up not only with the same grids, but also the same boundary forcings and same parameter settings. The main difference between the two was that one of them, the GC_MD model, was purely hydrodynamic, while the GC_MD mor model, was a coupled hydrodynamicmorphodynamic model. We validated both models using in-situ measurements at San Jorge Bay and Adaír Bay obtained with three ADCPs, with very acceptable results. We used three renewable energy indicators, namely the mean tidal spring flow speed, the tidal power density, and the annual energy production, to characterise the tidal energy resources in the region, and to analyse the effects of depth changes on the annual power production, AEP, after one year of simulation. We computed the AEP with the two models, and we found that the regions with the most substantial depth changes in the GC_MD mor simulations are well correlated with the regions where the AEP difference between the GC_MD and the GC_MD mor model predictions was largest. Most of these regions are found near the shore, with the area with most significant changes being west of Punta Choya, that is, west of the headland joining San Jorge Bay and Adaír Bay. A closer analysis of the seabed evolution along a transect off Punta Choya, based on a two-year simulation, showed that the seabed profile here is likely to evolve towards an equilibrium profile, with a smaller slope than that expected from Dean's for wave-driven beaches. Based on the morphological drivers and characteristics of macrotidal, sandy beaches, such result is expected in some regions near the shore. However, this tendency towards an equilibrium profile is not observed at other locations further offshore. Finally, future research should explore the impact of other forcings and the tidal energy devices on the environment and the available tidal energy resources.  Data Availability Statement: Some of the data for model calibration, model validation, and visualization of results will be publicly available in the ZENODO platform, under the same title as this paper.