We are grateful to Reviewer 1 for their thoughtful and detailed comments on this manuscript, and we have included responses below on a comment-by-comment basis. Our overall takeways from these comments are the following: (1) the spectral analysis is noisy and the features of the spectra discussed in the manuscript are difficult for the reader to see, (2) the correlation between the along-axis wind stress and the PC loadings is poorly presented, and (3) the indirect connection between the two surface tracers and the upwelling process weakens the conclusions drawn about upwelling. To address these issues, we propose to (1) apply the multitaper method to our spectral analysis to reduce variance, (2) rewrite our discussion of the correlation between PC loadings and along-axis wind stress to improve clarity and include significance testing and spectral coherence, and (3) use the improved comparison between the PC loadings and the wind stress record to build a stronger connection between surface nitrate and upwelling. We are confident that these revisions along with the proposed changes below will satisfy the concerns raised by Reviewer 1 and improve the overall quality of the manuscript.

Major Comments:

Major Comment 1 (Strength of surface tracers as proxies for upwelling rather than vertical density structure variables):

We agree that an analysis of the density structure is ultimately needed to better understand the mechanism of upwelling in the system. However the correlation between PC loading and wind stress for nitrate PCA modes 1 and 3 is statistically significant, and we find this correlation to be compelling and sufficient evidence that wind-driven upwelling is producing these coastal surface nitrate anomalies even though the vertical density structure remains unexplored. Since this study focuses on the influence of upwelling on surface nutrients, we believe that nitrate is the most relevant variable to present. We include temperature to demonstrate how the different surface processes between the two tracers manifest in their upwelling responses to wind, but the focus is nitrate. To better support our claim that coastal nitrate anomalies are upwelling-sourced, we will make the following revisions to the correlation analysis between the PC loadings and the wind stress record: (1) more clearly describe the direction, or sign, of the wind stress record being included in the correlation calculation, (2) quantify the significance of the correlations using a p-value test or equivalent, and (3) present the spectral coherence between the PC loadings and the wind record to show the dominant frequency bands where wind influence is strong. Finally, with respect to the original request for analysis of the density fields, we will be submitting a more physics-oriented article in the coming year that diagnoses the density structure more explicitly.

Major Comment 2 (Noisy spectra, consider block averaging):

We will use the multi-taper method to reduce the variance of the spectra without compromising low-frequency resolution, and we will report the 95% confidence intervals across the tapers.

Major Comment 3 (Poor presentation of wind stress correlation, consider bandpass filtering):

We will explore the prospect of bandpass filtering to remove high frequency data from our correlation analysis, and we will clarify our reporting of the correlation between wind stress and different modes to better separate the positive wind stress cases and the negative wind stress cases. We will also move the temperature modes II and III to the supplement to reduce clutter. In preparing this response, we have also found that the spectral coherence between wind speed/stress and the nitrate modes (especially mode 3 along the eastern shore) is consistently strong in the fortnightly band. We will include these coherence plots to further support the relationship between wind and nitrate.

Major Comment 4 (Difficulty seeing relationship between PC loadings and wind in Figure 7):

This figure is overly burdened by the presence of temperature modes II and III. We will move these modes to the supplement and minimize their discussion in the main text. This action will help declutter this figure, but we will also clarify the text with respect to the nitrate modes.

Major Comment 5 (Wind time scales poorly defined):

We will quantify the dominant wind time scales from the ECCC and HRDPS records and compare to recent atmospheric studies (primarily Bakri et al. 2017, Int. J. Climatol. and Thompson et al. 2020, Int. J. Climatol.).

Major Comment 6 (Significance of open end channels in 2-layer model):

We will add more emphasis to this aspect of the basin geometry as a caveat to the hypothetical along-axis pycnocline tilt. We elaborate below in our response to the comment at line 414-415.

Minor comments:

Line 3 (Attribute model tuning to previous works):

Will clarify.

Line 7 ("climatology" is confusing - rephrase to "predominant wind pattern" or "alongaxis winds steered"):

Will rephrase.

Line 30 ("Basin scale" here should be replaced by "dynamical width" or similar):

Will rephrase.

Line 92-94 (Please provide a reference for the estuarine circulation/exchange):

Will expand Pawlowicz et al. 2007 reference to include this statement.

Line 139 ("partial steps at the bottom boundary" - the meaning of this is not clear to me):

The "partial steps" aspect of the model is not directly relevant to this study, so we will remove the statement.

Line 146 (How are the auxiliary biogeochemical parameters like silica, plankton species, etc. relevant?):

Biological consumption indeed controls the residence time of nitrate in the surface, and the undiscussed tracers are important for resolving this sink accurately. This biological sink is fundamental to our choice of nitrate as a tracer. We will improve the emphasis and clarity of this point.

Line 181 (Please indicate, roughly, the timing of the freshet):

We will add a freshet timing definition to Section 2.1.

Line 260 ("Provides significant physical driver" - this statement should be qualified to be less strong):

Will rephrase.

Line 263 (Please explain why these particular locations were chosen):

Since the upwelling response is coherent along most of the coastline, the choice of these locations is somewhat arbitrary as they are intended to simply orient the reader to the upwelling behavior of the system, before continuing on to the more quantitative PCA. To make these choices less arbitrary, we will use spatial averages over 4 regions instead of 4 isolated locations.

Line 272 ("Averaged over the SoG region": Please be a little more specific):

We will include polygons on the map (Figure 1) to show the spatial averaging regions for wind, nitrate and temperature.

Line 297 (I find the use of "low-frequency" here unclear - especially since it seems to include the diurnal band):

Will rephrase as "sub-tidal"

Line 303 ("which represent" should be qualified, e.g. "which we interpret to represent"):

Will rephrase.

Line 306/328 (Temperature mode III seems to have a strong N-S structure, but it is referred to as cross-axis):

We will move temperature modes II and III to the supplement and trim the presentation in the main text to remove clutter and confusion.

Line 330 ("time-averaged" here is confusing - reads as an average across all time points):

This wind stress metric is central to our analysis of the EOF results, so we will define it more clearly in a separate sentence.

Line 334 ("small amount of correlation" is confusing, please rephrase):

We will define this correlation more clearly as the "correlation to the reverse wind stress direction" in a separate sentence and in our overall discussion of Figure 8.

Line 341 ("Visibly correlate" - avoid if no significant correlation was found in the quantitative analysis):

We will rephrase all discussion of correlation in quantitative terms.

Line 367 (Please explain briefly which assumptions have gone into transforming Ri(U) to Ri(tau)):

The primary assumption in this transformation is that the bulk vertical shear stress is described by the surface layer friction velocity squared (u*^2) which is equal to the kinematic wind stress (tau/rho). This assumption is common in simplified 2-layer models of lakes (e.g. Spiegel and Imberger 1986 J. Fluid Mech.). Regardless, we use this quantity as a scale parameter only and not an indicator of turbulence. Nevertheless, we will add a clarifying statement describing the origin of our bulk Richardson number and our specific use case.

Line 374 & 376 (Surely "the coasts" are always important?):

We will rephrase this statement to clarify that we are talking about timescales before the coastal pressure gradient forces have become established.

Line 394-396 (There is an apparent contradiction here):

This is a typo. Zeta_{end} does not contribute significantly to upwelling on typical wind time scales.

Line 414-415 ("If the cross-axis ... fluxes": This is not self-evident to me):

We are describing the behavior of the infinite coast scenarios of the two layer model at the corners of the rectangular basin at (0, 0) and (L_C, L_A). An analytical solution is not easily found, but qualitatively, the upwind lower layer transports along the basin sides combined with the cross-wind lower layer transports at the basin ends result in upwelling at the upwind corner (0, 0) and downwelling at the downwind corner (L_C, L_A). These corner signals will grow in time, and propagate anticlockwise along the basin end walls, setting up cross-axis pressure gradient forces at the basin ends that directly oppose the cross-shore Ekman fluxes. The flow near the ends that was cross-shore will respond to the new force balance by becoming more along-axis, eventually leading to an along-axis, lake-like pycnocline tilt. The open channels at the end of the basin mentioned by Reviewer 1 in Major Issue number 5 are a major caveat to this setup. We will rewrite the opening sentence of this paragraph to more clearly describe this scenario. We also propose adding a diagram of this qualitative solution to aid the reader.

Figure 1 (Please add a scale bar and improve contrast of the Texada star marker):

Will add/modify.

Figure 4 (Difficult to see the wind time series. Much of what is in the figure caption belongs in the text):

We will add a zoomed panel to emphasize the period between the cutoffs.

Figure 4 (Please indicate the timing of the snapshots shown in Figure 2):

Will add.

Figure 5 (Spectra should be computed based on the productive seasons - as for the profiles):

We will make this change in addition to using variance-reduction methods.

Technical corrections:

Figure 6 (Needs a length scale):

Will modify.

Line 239 - 240 (I recommend using a standard date format):

Will modify.

Line 551 ("right" -> left?):

This is a typo, will correct.

Figure 2 (Include the predominant wind direction):

Will add.

Figure 5 (Improve contrast between median profile and IQ range):

Will enhance.

Figures 5, 8 (Please add units to PSD y-scales):

Will add.

Figures 7, right (The sharp red color makes it difficult to discern the other time series):

Since temperature modes II and III account for the least variance and are not central to our findings, we will exclude them from all figures.

Figures 8 (bottom) (Change color of the horizontal line in case of difficulties for colorblind readers):

Will modify.

Title and elsewhere (Hyphenate "wind-driven" since it is a compound adjective):

Will change throughout.

We are grateful to Reviewer 2, Jennifer Jackson, for her thoughtful and detailed comments on this manuscript, and we have included responses below on a comment-by-comment basis. Our overall takeways from these comments are the following: (1) the criteria for identifying particular PCA modes as representing a given physical process such as upwelling or mixing is poorly stated or missing, (2) the correlation between the along-axis wind stress and the PC loadings is poorly presented, and (3) the applicability of the 2-layer model presented in 4.1 to the SoG and the improvement it brings to the paper are unclear. To address these issues, we propose to (1) clearly state the criteria for diagnosing the physical phenomenon represented by a given PCA mode as stated below in the comment-by-comment responses, (2) rewrite our discussion of the correlation between PC loadings and along-axis wind stress to improve clarity and include significance testing and spectral coherence, and (3) clarify the applicability of the 2-layer model for describing upwelling in the SoG and include additional figures to aid the reader in interpreting the model solutions in the context of the PCA results. We are confident that these revisions along with the proposed changes below will satisfy the concerns raised by Dr. Jackson and improve the overall quality of the manuscript.

Major Comments:

Major Comment 1 (Important information is missing from the description on the EOF results):

The percent variance is displayed in the bottom-left corner of each EOF panel, but we understand that it was difficult to see and will add a reference to the location of the label in the figure caption. The modes are diagnosed based on the following three criteria: (1) the spatial pattern, i.e., coastal anomalies in the case of upwelling, spatially uniform anomalies in the case of mixing or diurnal heating, anomalies in the tidal mixing zones for tidal mixing, (2) spectral energy distribution: i.e., weak tidal peaks and broadly distributed subtidal energy for upwelling, prominant diurnal peak for solar heating, prominant tidal and subtidal peaks for tidal mixing, (3) the PC correlation with a given segment of the along-axis wind stress record, i.e., correlation with positive wind stress indicates western shore upwelling, correlation with negative wind stress indicates eastern shore upwelling, correlation with both signs of wind stress indicates wind mixing, and no correlation indicates tidal mixing. These criteria were all satisfied for nitrate modes 1-3 and temperature mode I and we will emphasize and clarify this point in our revisions. The remaining temperature modes are inconclusive and account for the least variance of all modes explored, we will move them to the supplement. The mixing-heating pattern is diagnosed based on the three criteria described above, and we will clarify this diagnosis in our revisions. We acknowledge that our description of the wind stress correlations as "positive" and "negative" is unclear and will improve the clarity to be consistent with our EOF diagnosis criteria (3) above.

Major Comment 2 (Lack of discussion of stochastic events in the manuscript, and absence of spectral peaks):

All spectral analysis of wind records that we have performed yield broadly distributed energy at subtidal frequencies rather than narrow band peaks. We also see this broad frequency distribution in the analysis of Sand Heads wind velocity by Halverson and Pawlowicz 2016, Atmos. Ocean. We expect these spectra given the stochastic nature of storm intervals. However, based on the request of Reviewer 1, we have proposed to use variance reduction methods and spectral coherence in our revisions to improve our analysis of the wind, nitrate and temperature spectra. The coherence especially will help identify the frequency bands where the wind influence on raw tracer records and the PC loadings is significant.

Major Comment 3 (I suggest adding a table that details all of the mathematical symbols):

We will add a table of symbols.

Major Comment 4 (Emphasize why case studies in 4.1 are needed and how they influence the model results. Include figures):

We included section 4.1 to provide physical context for why we primarily observe surface tracer upwelling signals as coastal bands along the eastern and western shores instead of blobs at the upwind ends of the SoG. While we recognize that the vertical structure of the SoG is not well-represented by a 2-layer model at rest, a 2-layer assumption becomes increasingly appropriate during upwelling events because the wind rapidly mixes the stratified surface layer and the depth of upwelling never penetrates below the intermediate layer. Additionally, the upwelling solutions that we describe for the 2-layer model retain their basic structure as the number of layers is increased (e.g., Csanady 1982). We will clarify these points in our revisions and add additional figures to guide the reader.

Major Comment 5 (Some key references are missing):

We appreciate having these references brought to our attention. We will reference the 2021 papers in our discussion of the effect of surface nutrients on phytoplankton, sections 2.1 and 4.4. We will reference Johannessen et al. 2014 in our discussion of the effects of the tidal mixing regions on surface nitrate.

Minor comments:

Line 28 (Add references to previous research on upwelling in enclosed basins):

References to upwelling in enclosed basins are included throughout the introduction between lines 25 and 57, however we will add/emphasize a few key references here as well.

Lines 45 to 57 (I found this paragraph confusing and it was difficult to understand the point):

The purpose of this paragraph is to present the secondary features of a basin that may affect upwelling after the Rossby radius is considered. The concepts of wind stress curl, spatial salinity gradients and wave damping are all revisted later in the manuscript, but topographic waves are not relevant to this study. We will rewrite this paragraph to improve the linkages to the rest of the paper.

Lines 151 to 152 (Please add a reference here):

Will add a reference to Hansen et al. 2013, Harmful Algae

Line 165 (How realistic are the 2.5 km winds in narrow channels? Do they impact the results?):

Our windrose comparisons between the HRDPS and observed records at the four open water stations suggest that HRDPS is sufficiently realistic to resolve upwelling forcing. From other analyses unrelated to this paper, we have found HRDPS skill to be strong in all open water locations and channels, and to weaken only in narrow inlets such as Howe Sound but certainly the Discovery Island channels as well. Since these inlets are already isolated from the primary upwelling areas of the SoG, we find this level of HRDPS skill satisfactory. We will add a clarifying statement in our revisions.

Line 208 (Figure 1 includes Juan de Fuca Strait yet only the region to the tidal mixing area is considered):

Will add a box to Figure 1 indicating the subregion used for PCA.

Lines 209 to 213 (Difficult to interpret what the authors are stating here. I suggest possibly adding this information to a figure):

We agree that this paragraph introducing the PCA methods is poorly written. We will rewrite for improved clarity and possibly add a figure.

Lines 214 to 226 (Are the references at the end of this paragraph for the whole PC and EOF equations?):

We built our PCA/varimax Python code using algorithms from Preisendorfer 1988 and Horst 1965. Although the algorithms are essentially identical to established tools available in Matlab and SPSS, they are poorly documented in these code libraries. The goal of these references is to aid the reader in locating the exact source of each algorithm. However, we acknowledge that the way we have made these references is unconventional. We will clarify these references in our revisions.

Lines 257 to 258 (I don't understand the sentence starting with "There is also a tendency..."):

Southeasterly bins greater than 10 m/s contain a higher proportion of the total data at the northern stations relative to the southern stations in winter. We will clarify in the text.

Figure 2 (The letters in the figure to identify the panels do not match the description in the caption):

This is a typo. Will correct.

Lines 274 to 275 (I can't see this result in the figures):

The fluctuations in surface nitrate and temperature are a result of the study and visible in Figure 4, but the seasonal formation of the vertical gradients is background information and is not shown. We will clarify this distinction and add a reference.

Lines 291 to 292 (It is difficult to see the correlation between winds and temperature/nitrate at individual locations in Figure 4):

We will modify figure 4 or add a companion figure that zooms in on the region of interest so this correlation is more visible. We hesitate to add a scatterplot here since this figure is more of a lead in to the more robust PCA component of the study.

Lines 306 to 310 (As mentioned above, it is not clear to me how these interpretations were made):

The mixing-heating mode is diagnosed by the uniform spatial EOF pattern, the prominant diurnal energy peak, and the correlation to both positive and negative along-axis wind stress. Discussion of temperature modes II and III will be removed. We will clarify the diagnosis of temperature mode I.

Figure 7 (What do positive and negative winds and PC amplitude mean?):

Positive and negative refers to the direction of the along-axis wind stress. PC amplitude is equivalent to PC loading, but we will choose a single term to use throughout our revisions.

Line 336 (How does the averaging window of 54 hours impact the storm data?):

The PC loading is not really a function of instantaneous wind as much as the time-integrated wind. This dependence is demonstrated in the infinite coast solution in 4.1, zeta_{side}. The averaging process is analogous to time-integration since the average is just the discrete integral divided by the window length. In this sense, the averaging process should not effect the energy imparted from a storm. As the text says, the window length is simply chosen to maximize the correlation coefficient, and this process can be interpreted as finding the critical storm duration required to produce a PC anomaly in the given mode. We will clarify this distinction in our revisions.

Figure 8 (What do positive and negative PC amplitudes mean? Figure 8b shows significant energy at fortnightly and monthly frequencies. Why does tidal mixing have significant energy here?)

As previously stated, we agree that the presentation of correlation used poor language and our revisions will clarify this section. We interpret the fortnightly and monthly peaks to result from the fortnightly tidal cycle and the resulting change in mixing strength in the tidal mixing regions. We will clarify this interpretation in the text.

Lines 571 (Much of the observational data used are available on CIOOS. Consider an acknowledgement):

We will add this acknowledgement.