This very good paper has two main threads: (1) presentation of an efficient way to implement effects from sea ice in an operational ocean model and (2) an improved physical understanding of why so may tide gauges in the Arctic display seasonal perturbations in their tidal harmonics.  Both justify publication, so I urge the editor's acceptance.

For the main results, shown in Figure 4 (total water level), I was somewhat disappointed to see how little improvement was obtained. But my enthusiasm returned when the tide results (Figure 5 and after) were presented. These results are most encouraging.

I have only a few minor comments/suggestions.

Figure 3 caption. It would be useful to have more information about the source of the landfast frequency. Simply citing the Nat. Snow and Ice Data Center is not very informative.

Figure 4 caption (very minor). It might help readers (or browsers) to note that the missing stations are those for which only historical observations are available, and those data do not overlap with the model timeframe.

Line 19: I think a better and more comprehensive reference for Arctic sea ice changes is a recent paper by Parkinson: https://doi.org/10.3389/frsen.2022.1021781

The Ray(2022) paper was published by Ocean Sciences: http://dx.doi.org/10.5194/os-18-1073-2022

For Eqns (1-4), if you are using LaTeX, you can obtain better sized parentheses and brackets by using \left and \right.  See the manual.

The discussion in the Supplement of poor results at Nome, Alaska, is interesting, in part because NOAA has noted its water level predictions for Nome are the worst of any U.S. tide gauge.  See:

https://tidesandcurrents.noaa.gov/pdf/Tide_Prediction_Error_for_the_United_States_Coastline.pdf.

The tides at Nome are small, and non-tidal variability is large, according to a spectrum of tidal residuals (http://dx.doi.org/10.1357/002224017821836761  Figure 3), but that does not explain why O1,K1 display such large annual perturbations while the semidiurnal tides do not.  A mystery, although it could be related to discharge, as the authors speculate.

Summary:

This study demonstrates a methodology to efficiently incorporate sea ice effects into a global operational total water level forecast model and analyzes the impact to tides, surge, and total water levels in the Arctic region, including seasonal modulation and attenuation of water levels.

The study and manuscript is of high quality and is a useful contribution to the field and the journal. My comments below are mainly minor ones asking for clarification of methodologies and ideas.

Specific comments:

1. Page 3 L59-60: “However, such parameterizations depend solely on ice concentration, which is insufficient to represent the ice strength and its impact on the air-sea momentum flux transfer” . Suggest that the authors either provide reference(s) to support this assertion or more information on the mechanism by which the ice concentration-only parameterization must be insufficient.
2. Page 3 L61-62: “In Canada, sea ice effects on TWL forecasts are a major concern, particularly in the Canadian Arctic and possibly on the east coast of Canada.” Why particularly Canadian Arctic and “possibly” the east coast? Is this because that’s where the most ice/ice change is or where the winter storms mostly impact? Please provide more detail. 
3. Page 6, L125, Eq (5): Suggest to add the condition for τ_s = τ_ao, i.e., when ice concentration, α = 0. 
4. Page 7 L136: Any reason why newer versions of GEBCO (2019-2022) bathymetry are not used?
5. Page 9 L208-210: The ice-ocean tide rotation angle values are described but not shown. Suggest to add row to figure 3 or add new figure to supplementary for rotation angle, φ.  
6. Page 10 L230: I’m not sure I quite understand what the scale factor, a^S is, how this factor is determined, and why it is a spatial constant. On L182 it says it is used to account for differences in valid depths between the external model and the TWL model. Does this mean the levels at which the model surface currents are valid or the seabed depths? But the text also seems to say it is actually determined from observations (which kinds? Section 2 does not describe any velocity data) not by comparing the models. I think this scale factor needs more and clearer explanation.
7. Page 11 L240. It is not clear to me how the following “processes" are removed:
   
   Run 2: TSI due to τ_io and Run 3: TSI due to τ_b. I think I can guess that in Run 4, u_ice^T - u_surf^T used in Eq. (9) is set to zero? But how is the tide-surge interaction removed from τ_io and τ_b? Please explain and/or write the mathematical conditions.
8. Page 13/14 L295-300. The authors mention the anomalously high positive M2 amplitude change for station 18 was not captured by the Run_AIO (actually predicts a negative change), and associates it with displacement of local amphidromes. However, Run_AIO captures the phase change very well which would seem to contradict this explanation. More evidence to explain this discrepancy is warranted. 
9. Page 21 L384: How did you remove, η_∈ “contributions from other process such as seiches”  to get the surge, η_S? Or is η_S really just the non-tidal residual here?
10. Page 21 L395-398 and Page 25 L455-459: Any ideas how to adjust for the underestimation of ice and current velocities on the northern Alaska coast due to insufficient ice break-up? The ice concentration-only based parameterization used in Joyce et al. (2019) gives the maximum ice-ocean drag coefficient at 50% ice concentration due to the “interplay between the increased number of ice floe edges and additional sheltering of the atmospheric flow downstream of ice floes with increasing sea ice concentration”. Although the ice velocities are underestimated in GIOPS, what about the ice concentrations? 

Technical corrections:

• Page 12 L265: “constitute” to “constituent”.