Pilcher et al. used an ocean-biogeochemical model to investigate carbonate system variability in the Bering Sea, with a focus on the long-term trends of ocean acidification (OA) variables (Ω, pH, H+, pCO2), motivated by the need of relevant OA indices for marine resources management. They extended the temporal coverage of previous modeling studies, quantifying spatiotemporal trends during 1970-2022. Their model results showed a significant acceleration of the OA trend over the last 25-years. The simulated bottom trends are greater than the surface trends, presumably associated with an increased respiration/remineralization in response to enhanced phytoplankton production. This is an interesting and valuable study that contributes to better understand changes in the carbonate system in the Bering Sea, but I think additional work is required to clarify/improve the model settings and further explain the model results. 

Main comments

1) Model forcing

The main results in this study relate to the long-term patterns in OA variables. A relevant question is therefore how robust the derived 1970-2022 trends are, and I have some concern about this. The author recognized potential issues associated with the use of multiple products to derive the surface model forcing (CORE for 1970-1994, NCEP-CFSR for  2011-2021, and NCEP-CFSv2 for 2011-2022) and boundary conditions of physical variables (Northeast Pacific model (NEP) for 1970-1994 and CFS for 1995-2022). Model results show that the CORE-to-CFS transition most likely altered the long-term patterns in the OA progression. Maybe I am missing something, but there are available  atmospheric reanalysis products that cover the entire study period (e.g., ERA5: 1950-2022), so I am not clear why the authors decided to use those three different products. Could you clarify this?

An additional CORE-forced hindcast ending in 2003 was intended to clarify the CORE-CFS transition impact on the trends, but I do not think this extra analysis significantly helped to that goal. Instead of comparing if the CORE-CFS trends for 1970-2022 are like the CORE trends for 1970-2003, I would compare if trends derived from CORE and CFS forced experiments are consistent during the overlapping period of these two products. A comparison over the overlapping period could also help to clarify potential impacts on seasonal and interannual variability.

Maybe you could re-run the full hindcast using ERA5 as the only atmospheric forcing, so that the “forcing issue” would be limited to the NEP-to-CFS shift in the boundary conditions. In that case, you can run additional experiments to compare trends over the overlapping periods of NEP and CFS.

2) Salinity trend shift

Figure S1 show a strong trend shift in salinity associated with the change in the model forcing products (CORE to CFS).  I suggest reporting the mean alkalinity series at surface and bottom as supplementary figure, since salinity and alkalinity are usually strongly correlated. Did you get a similar trend change in alkalinity? Since salinity and alkalinity are drivers of Ω, pH, and pCO2, then that shift could significantly impact all the reported OA trends. Could you discuss about it? 

3) pCO2 patterns

There are not many observations during fall-winter, but the M2 records in 2021 suggests that the model  is overestimating pCO2 during those seasons. If that is correct, then you could conclude that the model has an overall positive bias in pCO2 (and negative bias in pH), which could have a strong impact on the air-sea CO2 fluxes. You could discuss about it and tone done your results related to CO2 fluxes.

4) Underlying drivers

A Taylor series decomposition could provide valuable insights about the underlying drivers of OA variability, helping to identify the causes for the carbonate system trend changes, support the hypothesis of a biological driven increase in the bottom OA trends, and identify if salinity and alkalinity play any role on the OA progression changes. Although the authors mentioned  that diagnostic mechanisms are beyond the scope of this study, I strongly recommend adding a Taylor decomposition analysis, especially considering that a paper describing historical OA pattern was already published (Pilcher et al., 2019).

Minor comments:

I suggest tone done all the modeling results, considering that an ocean-BGC model rather suggest than demonstrate the OA patterns. For example, in the abstract: “surface Ωarag decreases by -0.043 decade-1 and surface pH by -0.014 decade-1” => “model results suggest that Ωarag decreases by 0.043 decade-1 and surface pH by -0.014 decade-1”

168-171: Could you provide the spatial resolution of the Northeast Pacific model and the CFS products used for the boundary conditions?

174: How did you estimate the empirical climatological profiles for iron?

393: “variables are comparable” what do you mean?

415: “seasonally occurring” => provide the season name.

455: I suggest using the same x-axis range (1970-2022) for the two panels, independently that the M8 station data extend only until the 80s.

This study present a long term hindcast (1970 to 2022) of the carbonate system on the Bering Sea shelf. Such long term hindcast are very useful to understand the impact of environmental changes on the ecosystem and are very useful for fisheries management. The model results show an acidification trend, which is more important at the bottom than at the surface. The trends of the different OA variables are reasonable and within the range observed elsewhere. However, there is no atmospheric forcing data that covers the whole period of the hindcast and thus 3 different products were used. This approach seems to introduce an artificial trend when the forcing is switched from one product to the other. This issue is well discussed in the manuscript but unfortunately it creates uncertainty in the results. I think a little more could be done to understand what is going on when moving from the CORE forcing to the CFSR forcing. It would be very useful to know how the temperature and salinity change at the open boundaries of the model, is it why a general decrease in salinity is observed? How is the salinity-DIC relationship modified, does this lead to more or less DIC flowing into the system. Water from which depth at the boundary will affect the Bering Sea shelf? What is the impact of not modifying the TA relationship at the boundary? Some processes affecting DIC (not the air-sea flux), water mass changes, can also slightly affect TA and thus the relationship. Modifying one variable but not the other could introduce a bias in pH and omega. Was this explored with the Earth System Models or only dismissed? And finally, the large trend towards the end of the hindcast was attributed to an increase in primary production. Could you explain why PP is increasing? The salinity gradient between the surface and the bottom seems smaller with the CFSR forcing, is there a decrease in stratification that leads to a greater nutrient availability? Or is there more light available? How does the PAR forcing (SW, clouds) compare between the different products? Is it a temperature effect (no info is given on the temperature change)?

To summarize, I think it is a very useful study but giving a little more information on the differences between the different forcing, and some explanations on the causes of the observed changes (i.e., salinity, primary production) would help build confidence in the results. The manuscript is very well written and clear. I think once some of the above questions are answered and delt with if necessary, the manuscript would meet all the criteria for publication in this journal.

Specific and technical comments :

Line 58 : Could refer to Figure 1

Line 107: lower pH values -> not clear lower than what. It would be better to actually put the pH values.

Line 110: Seung and Punt not in the reference list.

Line 111: on the Bering …

Line 112: some discussion on the role of OA: who had these discussions, any ref to a workshop or something like that?

Line 125: could you specify what was the length of the previous hindcast?

Line 147: add parenthesis (2020).

Line 160: How many rivers are there in the model? How many rivers have discharge and DIC, TA values? The same values are added to all the rivers that do not have measurements (if that is the case)? What do the TA and DIC values look like (i.e. range of values). What are the largest rivers in term of runoff. I think some more info would be useful especially that the river data seem to have a large impact in some regions of the Bering Sea.

Line 166: 1995-2010 instead of 2011? Why didn’t you use the CFSv2 reforecast for your whole 1995-2021 period? These 3 different forcing actually creates 2 transitions.

Lines 175 and following: The adjustments made to the relationship and the resulting DIC over the full salinity range is not clear to me. I think a figure would help. Also could you name the 3 ESMs used? In the earth system models the DIC changes include more than the results of atmospheric CO2 influx and other processes could impact TA as well. So at depth it might be relevant to include the changes in TA as well (from the ESMs) to preserve the TA-DIC balance.

Line 120. Wd should be in italics

Line 265: How are the off shelf break water affecting shelf water? A related question is how are the model boundary conditions affecting the shelf water, i.e. only the properties at similar depths (down to 200 m) or deeper?

Line 282: Why is omega more variable?

Figure 4: There is a smaller peak that precedes the larger peak in the model. What is it caused by? In the observation we rather see a slowdown of the increase.

Line 330: Would the TA and DIC bias originate from your treatment at the open boundary (i.e, increasing DIC and not TA).

Line 356: This decrease in salinity will lead to important changes in the DIC at your open boundary (?) and thus to the important change in H⁺, pH and omega on the bottom? This comes back to my previous comments about the usefulness of actually displaying the changes in the S-DIC relationship and mentioning the importance of the OBC on the shelf conditions.

Line 397: Could you add a reference?

Line 417: Wrong x axis at M8

Line 460-461: Could you say why? What are the changes in temperature between the two forcings?

Figure 10: p missing in DpCO2 panel title.

Line 510: Why is primary production increasing with the CFSR forcing? Is there a change in stratification (the difference is surface and bottom salinities seems to suggests that) that would provide more nutrients to the upper layer? Any changes in PAR from atmospheric forcing?

Line 511: Any oxygen observations that would support the remineralization rates at the bottom?

Figure 15: replace are by area in legend