Downscaling global ocean climate models improves estimates of exposure regimes in coastal environments

Study region

Monterey Bay is the largest bay on the California coast, and has a total area of approximately 1,162 km² (Fig. 1). Due to upwelling, the bay is highly productive and supports dense kelp forests, dominated by Macrocystis pyrifera, on rocky reefs that extend to approximately 20 m depth²⁷. These forests have the capacity to regulate pH^28,29,30 and sustain over 200 different species, from phytoplankton to marine mammals³¹.

Temperatures are the coolest in Monterey Bay during spring and early summer due to wind-driven upwelling³², driven by equatorward winds blowing over the California coast causing offshore Ekman transport³³. During the same months, filaments of water originating from Point Año Nuevo to the north are trapped within the Bay, where they warm due to solar radiation, forming a lens of warm water close to shore^34,35 (Fig. 1). Furthermore, a weak cyclonic eddy is observed within the bay due to the coastal geometry^32,34. Inside the bay, sea breezes and tides drive diurnal and semi-diurnal currents that can lead to significant variability in environmental conditions^17,36. During the upwelling period, salinity is approximately 34, temperature ranges from 9 to 13 °C at 17 m depth⁷, DO varies from as low as 100 μmol kg⁻¹ (3.2 mg L⁻¹ for T = 13 °C, S = 34 at 17 m) to as high as 300 μmol kg⁻¹ (9.62 mg L⁻¹ for T = 13 °C, S = 34 at 17 m), and pH varies from 7.7 to 8.1³ (Supplementary Fig. 1).

The three primary barotropic tidal constituents in the region (M2, K1, and S2) are responsible for over 80% of the tidal amplitude observed¹⁸. In the southern region of the bay, tides are mainly responsible for the cross-shelf velocity³⁶ and the interaction of these surface tidal currents with the steep topography create internal tides comprised of internal waves and bores occurring at tidal frequencies^37,38. Internal tides are formed when currents move over steep slopes, dense waters are forced into shallower regions, and, as these waters sink back to depth, an internal wave is generated³⁹. These internal waves have speeds⁸ on the order of 0.05–0.2 m s⁻¹ and, as the coast steepens, can break forming internal bores that move upslope and bring cold, low DO, low pH waters into nearshore kelp forest ecosystems⁴⁰. During upwelling, these processes drive variability in temperature, pH and DO over semi-diurnal and diurnal periods that can exceed the predicated changes in mean conditions predicted by global climate models for year 2100⁴¹.

Model description

To better understand how tides and winds affect exposure of nearshore organisms to variability in temperature, DO, and pH under current and RCP climate scenarios, we used dynamic downscaling and developed a 2D coupled biogeochemical hydrodynamic model using ROMS⁴². The model domain was created based on the Monterey Bay continental shelf described in Walter et al.⁶ with a maximum offshore depth of approximately 80 m (Fig. 2). The biogeochemical model used is described in Fennel et al.^43,44. We forced the model with representative wind (diurnal sea breeze accounting for regional winds), solar radiation, and tidal currents for the Monterey Bay region (Supplementary Figs. 2–6; Supplementary Table 1). Detail of the full model structure, including boundary and initial conditions, are included in the “Supplementary Information”.

We estimated depth profiles for temperature, DO, and pH for the downscaled coupled hydrodynamic-biogeochemical model. We considered a homogeneous salinity (S = 34) for the entire domain since the variation in salinity is less than 0.4 over an entire year and less than 0.05 during the upwelling season in the Monterey Bay region²⁷. Thus, initial and boundary stratification are assumed to be controlled by temperature alone. We represent initial (IC) and boundary conditions (BC) for temperature, DO, and pH using:

where var is the identified variable (e.g. temperature), z is the depth, α is a fit coefficient for each variable determined from a least-squares fit to observational data, and D_pyc is the depth of the pycnocline. We used this method to estimate profiles of temperature, dissolved inorganic carbon (DIC), and DO for present and all downscaled future scenarios setting D_pyc to 17.5 m. The data used for the initial and boundary conditions (BC) of phytoplankton and chlorophyll profiles were taken from Schuckmann et al.⁴⁵. The initial chlorophyll concentration was converted to zooplankton concentration (zoop = 0.34 × 10^–3 mmol m⁻³) using Eq. 3 from Wiebe⁴⁶, a method that has been used in other studies^47,48. Detritus was initially set to zero in the entire domain. All biogeochemical variables were forced hourly at the southern boundary.

For each scenario, we estimated depth profiles for temperature, dissolved oxygen, and DIC using Eq. (1). In order to fit Eq. (1), we obtained our estimated values of each variable near the surface, at 80 m depth on the shelf during upwelling for Present and 200 m depth for Future, respectively. We used 200 m depth from the future data set as this is the most common depth of source waters for upwelling in the region⁴⁹. Present surface and bottom values of temperature, DO, total alkalinity (TA), and DIC for present scenario were based on Koweek et al.²⁷. The mean of the 3-month period of strong upwelling (May, June, and July) was used to obtain values of temperature (surface and depth), oxygen (surface), and DIC (surface) from Representative Concentration Pathway (RCP) for the year of 2100 from the 4th report of the IPCC⁵⁰. Since only surface values for DO and DIC and no values for TA were available, we estimated values for these parameters at depth.

For DIC, we assumed that the ratio between surface and bottom values (80 m for Present and 200 m for Future) in present conditions will not change for future scenarios:

Therefore, to find the bottom values for future scenarios we divided RCP surface values by this ratio:

Surface and bottom TA values were kept the same as present conditions, following Feely et al.³. For DO, bottom values for RCP 2.6 and 8.5 were approximated from Figs. 5 and 6 of Bopp et al.⁵¹. We calculated the ratios between surface and bottom DO for Present, RCP 2.6, and RCP 8.5. The ratio of surface:bottom DO was not constant across scenarios, so we approximated the ratios for RCPs 4.5 and 6.0 using linear least squares fit (Table 1; example calculation in the “Supplementary Information”). We, then applied the values from Table 1 in order to calculate α assuming D_pyc = 17.5 m. We then used Eq. (1) to generate the initial and boundary conditions (Fig. 3).

We calculated pH and Ω_ar using the CO2SYS⁵² package in MATLAB using temperature, salinity, DIC, and Total Alkalinity (TA) from the simulations at the offshore location where the bottom depth was 15 m. We assumed concentrations of phosphate and silica based on Koweek et al.²⁷. We used dissociation constants for H₂CO₃ and HCO₃ from Dickson and Millero⁵³ and hydrogen sulfate ion constant (HSO^-4) from Dickson⁵⁴. All surface oxygen values were shifted positively 60 mmol m⁻³ (1.87 mg L⁻¹ for T = 13 °C and S = 34) in order to simulate the high primary production due to kelp forests²⁹, which is not specifically accounted for in the model. Overall, temperature at the bottom remained constant across all scenarios, with exception of RCP6.0 where temperature increased 0.2 °C (Supplementary Table 2).

Integrated exposure

Field observations^24,55 and laboratory experiments^5,56 have shown that below sub-lethal thresholds marine organisms inhabiting nearshore marine habitats in upwelling systems exhibit signs of physiological stress when exposed to elevated temperatures, low oxygen levels, or low pH waters, which is especially detrimental especially for calcifying species. Exposure of organisms to stressful temperature, DO, and pH conditions (φ_th where φ refers to temperature, DO, or pH) was done by subtracting the threshold value for a given organism and life stage from the model or observational data at 15 m water depth, then setting all positive values to zero for pH and O₂, and all negative values to zero for temperature. Next, we estimated integrated exposure (E_int) by integrating absolute exposure over a period of a week with a window interval of 1 h:

Thresholds of temperature, dissolved oxygen and pH (16 °C⁵⁷, 4.85²⁵ mg L⁻¹, and pH of 7.5²³ respectively) representing non-interactive negative impacts on juvenile red abalone growth were based on literature values^5,23,57. Integrated exposure quantified the time and degree of stress an organism experiences, similar to the degree heating week with units of ^oC w, or day (^oC d) measure used to estimate thermal stress on coral reefs⁵⁸, and has been previously used to understand the exposure of juvenile abalone populations to similar stressors in an empirical field study⁵⁷. Overall, red abalone threshold values chosen had a strong negative effect on the species^5,23,57, therefore, we used E_int as a proxy for estimating the potential impact of future conditions on abalone growth and survival.

Fertilization response

We estimated fertilization success using results of Boch et al.²⁵ where the fertilization response of red abalone (Haliotis rufescens) was quantified in response to multiple stressor climate conditions (high temperature, low DO, and low pH). Fertilization in abalone occurs over relatively short periods, therefore E_int would not provide an appropriate estimate in such cases. While the process of fertilization occurs over short periods, adult red abalone exhibit an extended spawning season, over which environmental conditions may vary greatly based on our modeling results. Thus, we used the equations from Boch et al.²⁵ to examine how fertilization success over a one-month period might be affected by environmental variability, specifically the interactive effects of ph and temperature. Changes in DO did not show a strong effect on fertilization in their experiments (Fig. 4; see Supplementary Table 3 for parameter values):

For temperatures = 13 °C:

For Temperatures = 18 °C:

where ({upbeta }_{0}) and ({upbeta }_{mathrm{pH}}) are intercepts, (upbeta ) and ({upbeta }_{2}) are slope segments, ({upbeta }_{mathrm{A}}) is slope of the pH-temperature interaction (pH × Temperature Group), ({upbeta }_{mathrm{B}}) is accounts for high temperature effects, and BP is the curve breaking point. Since only curves for 13 °C and 18 °C were available, we used 15.5 °C as a transition where Eq. (5) was applied for temperature less than 15.5 °C and Eq. (6) was applied for temperature greater than 15.5 °C. Since more complex interpolation schemes yielded similar results, we used this straightforward method for clarity.

Model evaluation

We first assessed whether the oceanographic model was able to reproduce current (observed) oceanographic conditions in Monterey Bay. The model was expected to reproduce the main dominant semi-diurnal and diurnal periods of oscillations observed in the region as well as primary production regulation of DO and pH (excluding kelp) in the biogeochemical model. We had two requisites to consider the model performance satisfactory. First, we required that the model was capable of reproducing local minima for DO. Second, the model needed to be capable of reproducing the mean and extremes for temperature and pH. In addition, we anticipated observing lower values of DO in surface waters since we did not account for the higher primary productivity observed in kelp forests²⁰, and therefore, DO oversaturation. To evaluate the biogeochemical model, we compared chlorophyll concentration in the model to satellite-derived chlorophyll estimates for the region. Monthly mean depth averaged chlorophyll by area (mg Chl m⁻¹) was calculated for the model and for the months of May, June, and July from Sea-Viewing Wide Field-of-View Sensor (SeaWiFS)⁵⁹ for the period of 2010–2017 to the closest region with data available next to our model simulation (cross section region in Fig. 1). Chlorophyll concentrations in the model were 2.12 mg Chl m⁻³ compared to 4.02 mg Chl m⁻³ estimated from SeaWifs. Thus, modeled values were within the range observed in the satellite data over this period.

In order to validate temporal variability in the model results, we estimated power spectra using the Thomson Multi-taper method (MTM)⁶⁰. This method was chosen due to its robustness for stationary data with low variance. Power spectra allowed us to quantify the variability by frequency and confirm that the model was reproducing variability at dominant periods (M2, K1) observed in the region, as initial and boundary conditions in the model were based on regional observations. We applied the analysis over a 3-week window of upwelling for temperature, DO, and pH and compared with Booth et al.⁷ data (Fig. 5). Spectral analysis was used to ascertain the dominate temporal constituents of temperature, DO, and pH for all scenarios (present and RCP) between 12- and 24-h periods. The spectra were then integrated numerically to calculate the variability of temperature, DO, and pH, and the 95% confidence interval at each frequency band. Bootstrap analysis was applied to the integrated exposure calculation for each variable. In order to assure convergence in our estimates, 1,000 iterations with replacement were done using a sample size of 50% of the data. In the end, the confidence interval (CI) was calculated based on 2.5% percentile and 97.5% percentile of the distribution of the estimated means.

Our model results reflect current variability in temperature, DO, and pH in southern Monterey Bay (Fig. 5). In addition, mean temperature and pH values were not significantly different from in situ data. However, mean DO in our model was lower than present day averages, likely due to the lack of oxygen-producing kelp in our model⁶¹. Thus, the model accurately simulates diurnal, semi-diurnal, and higher order tidal components⁶². For Monterey Bay, diurnal (cycles per day (cpd) = 1) and semi-diurnal (cpd = 1.93) are the main temporal components of variability in temperature, oxygen, and pH and are well represented by our model with overlapping 95% confidence intervals (Fig. 5b,d,f). Higher frequency variability (cpd = 3) is also within the 95% confidence interval when comparing model and observations. However, frequencies occurring between peaks are not well resolved (Fig. 5b,d,f). We expect this is because we used a 2D model, and therefore, the model was unable to resolve all the physical processes occurring in the Monterey Bay. The oscillations between peaks in the observed data were likely due other coastal ocean processes such 3D circulation and ocean surface waves³⁶ (as well as noise in the instruments). However, the variability at these periods did not have an appreciable effect on exposure calculations. Importantly, the model preserved the observed diurnal and semi-diurnal variability that is not present in global and regional scale climate models for future RCP scenarios (Fig. 5a,c,e).

Variability at 15 m

The model was able to simulate the main drivers of variability in temperature, DO, and pH (predominantly internal waves) for the region in study. Internal waves in the domain were seen as vertical changes in the u-component of the velocity (Fig. 2). Cross-shelf velocities ranged from − 0.05 to 0.1 m s⁻¹, and were within the range found in other studies^8,18,27. Before the arrival of the internal wave crest, isotherms were tilted downwards indicating previous downwelling. Flow in opposite directions between crests was observed during retreating of internal waves in the domain, as it has been observed in studies on internal waves with in situ data^8,9.

High variability in temperature, DO, and pH was observed in all model runs (Supplementary Figs. 7–9). Mean temperatures were 10.63 °C (SD = 0.39), 13.81 °C (SD = 0.46), 15.46 °C (SD = 0.50), 15.45 °C (SD = 0.49), 16.96 °C (SD = 0.95) for present, RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5, respectively (Supplementary Fig. 7). Overall, mean temperature increased as expected and exhibited similar variability (standard deviation [SD]) across all RCPs except for RCP 8.5 scenario. This was likely due to increased temperature stratification where surface waters warm faster than deeper waters resulting in higher temperature variability.

Mean dissolved oxygen values were 5.15 mg L⁻¹ (SD = 0.74) for present, 4.80 mg L⁻¹ (SD = 0.88) for RCP 2.6, 5.32 mg L⁻¹ (SD = 0.77) for RCP 4.5, 5.02 mg L⁻¹ (SD = 0.81) for RCP 6.0, and 4.93 mg L⁻¹ (SD = 1.21) for RCP 8.5 (Supplementary Fig. 8). The SD for RCP 8.5 was almost twice that of other scenarios. This was likely due to a stronger gradient in oxygen (surface remains saturated while the values at depth are lower). RCP 2.6 had the lowest mean DO but not different than the rest and RCP 4.5 was 0.52 mg L⁻¹ higher, though not significantly different than the present scenario. The consequence for the lowest mean DO in RCP 2.6 was related to the strength of density stratification⁶¹, and therefore, low oxygen at 15 m. Another scenario (not shown) was used where RCP 4.5 oxygen profile was applied using the density stratification from RCP 2.6 and the same low oxygen values found previously were also observed in the alternate run, supporting this mechanism.

pH variability was also high across all model runs (Supplementary Fig. 9). The mean value for pH decreased from 7.73 (SD = 0.07) for present to 7.44 (SD = 0.12) for the RCP 8.5 scenario. RCP 8.5 again had the highest variability among all scenarios. Otherwise, the pH range was ~ 0.075 for all other scenarios. Lower pH for the most extreme scenarios (RCP 6.0 and 8.5) has also been observed in large scale models¹².

Source: Ecology - nature.com