A methodology for estimating the response of the coastal ocean to meteorological forcing : A case study in the Bohai Bay

The sea level (SL) variations at the coastal ocean result from multiscale processes and are substantially contributed by the SL changes due to the meteorological forcing. In this study, a new methodology, named as IBR, is developed to estimate the response of the coastal ocean to meteorological forcing. The response is taken as the combination of the static ocean response calculated using the inverted barometer formula and the dynamic ocean response estimated using the multivariable linear regression involving atmospheric pressure and wind component at the dominant wind orientation. 15 The dominant wind orientation is determined based on the averaged values of the magnitude squared coherences between the adjusted SL and wind at every wind orientation. The IBR is implemented to estimate the response of the coastal ocean at two stations, E1 and E2 in the Bohai Bay, China. The analysed results indicate that at both E1 and E2, the adjusted SLs are related more to the regional wind, which is the averaged value in the Bohai Bay of the 10 m wind in the ERA-Interim data, than to the local wind; the dominant regional 20 wind orientation is 75°. The estimated response using IBR with the regional meteorological forcing is much closer to the observed values than other methods, including the classical inverted barometer correction, the dynamic atmospheric correction, the multivariable linear regression and the IBR with local forcing, demonstrating that IBR with regional forcing have the best skill in estimating the response. The large deviations between the observed values and the estimated values using IBR with the regional meteorological forcing are mainly due to the remote wind, which is not considered in the IBR. 25 This case study indicates that the IBR is a feasible and relatively effective method to estimate the response of the coastal ocean to the meteorological forcing.


Introduction
The elevated sea level (SL, see Appendix A) at the coast is dangerous for the nearby city, while the depressed sea level can render navigation of coastal bays and harbors difficult and hazardous (Tilburg and Garvine, 2004).The SL variations at the coast result from multiscale processes, with the superposition of global mean SL, regional SL and local SL changes, as shown in Melet et al. (2016) and Melet et al. (2018).SL changes due to the meteorological forcing, including sea surface atmospheric pressure (AP) and wind, make substantial contributions to total SL changes in coastal ocean (Melet et al., 2018;Zhang et al., 2019).Traditionally, the ocean response to meteorological forcing is estimated simply using the inverted barometer (IB) correction, which makes the response be poorly accounted for (Wunsch and Stammer, 1997).The classical IB formulates the static response of the ocean to AP forcing, and the wind effects are totally ignored (Carrè re and Lyard, 2003).However, the coastal ocean response is different from the IB response in shallow water and the influence of the wind forcing cannot be ignored (Carrè re and Lyard, 2003;Lv et al., 2018).
The dynamical atmospheric correction (DAC), which is a combination of the high frequencies of simulated response forced by AP and wind using numerical model and the low frequencies of the IB correction (Carrere et al., 2016), has been widely implemented to estimate the barotropic response of the global ocean to meteorological forcing, especially in the altimetry SL estimations.Tierney et al. (2000) pointed out that the simulated SL signals at periods less than 20 days, which were not uniform in space but can be enhanced in some semi-enclosed regions (Fukumori et al., 1998), were more consistent with observations when using AP plus wind forcing than using wind forcing alone.Hirose et al. (2001) used a barotropic, shallow water model to simulate the response of the ocean to atmospheric disturbances and found that the AP-driven results accounted for a small portion of the observed SL variance while the wind-driven results explained a large amount of the variability in the observations at middle and high latitudes.Carrè re and Lyard (2003) simulated the global ocean response to atmospheric wind and pressure forcing using the MOG2D-G model, and found the model correction could reduce the SL variance at high latitudes, continental shelf areas and shallow waters compared to the classical IB correction.Melet et al. (2016) evaluated the relative importance of processes causing coastal SL variability at different time-scales, in which the SL induced by AP and wind was estimated using the DAC.
The regression analysis was also widely implemented to investigate the response of oceans to the meteorological forcing.Tilburg and Garvine (2004) developed a simple linear-regression model based on modest wind forecast capability and records of local coastal sea level, wind and pressure, and found that this empirical model was adequate for general use.
Based on the results of regression analysis, Andres et al. (2013) hypothesized that the annual SL changes along the shoreline were significantly influenced by the local winds and barotropic response.Calafat and Chambers (2013) found that a multivariable linear regression (MLR) involving local wind and AP can account for a substantial fraction of the annual SL variation at the Boston and New York tide gauges.When the MLR involving the local wind and AP was implemented, Lv et al. (2018) found that the response of the coastal ocean to AP had the similar form to the IB correction and the estimated results were much closer to the observations than those using only wind or IB correction or combination of them.
In the DAC, the high frequencies of simulated response forced by AP and wind and the low frequencies of the IB correction represent the dynamic ocean response and the static ocean response, respectively (Piecuch et al., 2016).It is noted that the fine and accurate topography data in the coastal and shallow area is difficult to be obtained, which strongly influences the wind-driven barotropic SL variability (Fukumori et al., 1998), so the dynamic ocean response simulated by numerical model may be not accurate in the coastal and shallow area.Moreover, the static ocean response is always ignored in the regression analysis.Besides, the SL change may be not determined by the local meteorological forcing, as shown in Amin (1982), Thompson et al. (2014) and Piecuch et al. (2016).
The purpose of this paper is to present a new methodology, in which the response of the coastal ocean to the meteorological forcing is estimated using the combination of the static ocean response, which is calculated using IB correction, and the dynamic ocean response, which is estimated using regression analysis, and to demonstrate the feasibility and effectivity of this new method.Besides the local meteorological forcing as shown in Lv et al. (2018), the influence of the regional meteorological forcing will also be considered in this study; in addition, the low-pass wind and AP are used in the regression analysis in this study, rather than the original forcing as used in Lv et al. (2018).The details of the rest of the paper are as follows: the new methodology and data are shown in section 2; the estimated response of the coastal ocean to the meteorological forcing in the Bohai Bay is shown in section 3; the comparisons with the other methods and the reason for the large deviations in the estimated results are discussed in section 4; and the conclusions can be found in section 5.

Methodology for estimating the response of the coastal ocean to meteorological forcing
As shown in Sandstrom (1980), Hsueh and Romea (1983), Castro and Lee (1995) and Lv et al. (2018), the response of the coastal ocean to meteorological forcing is mainly at low frequencies and is the combination of the following two main components (Piecuch et al., 2016): where, hL is the low-frequency SL, and hstatic and hdynamic indicate the static ocean response and the dynamic response, respectively.
According to the DAC, the static ocean response can be estimated using IB correction.However, the dynamic ocean response in the DAC, estimated using the simulated results by numerical model, may not be accurate in the coastal and shallow area, as the resolution of topography data is not always high enough; on the contrary, the regression analysis is not affected by the resolution of topography data.The response of the ocean to meteorological forcing, including AP and wind, is approximately linear (Li and Yang, 1983), as demonstrated in Calafat and Chambers (2013) and Lv et al. (2018), so the dynamic ocean response can be estimated using the linear combination of AP and wind.As the static and dynamic ocean responses are linear in the meteorological forcing, the low-frequency AP and low-frequency wind will be used to estimate the ocean response to avoid introducing the high-frequency signals into the estimated low-frequency SL.Therefore, the response of the coastal ocean to meteorological forcing is estimated as follows: where, LowAP is the low-frequency AP; LowWind is the low-frequency wind component at the dominant wind orientation, which can be determined based on the result of magnitude squared coherence, as shown in Lv et al. (2018); hIB is the SL estimated using the IB formula with LowAP; α0 is a constant, α1 and α2 are coefficients.
As the IB formula is explicit, as described in Paraso and Valle-Levinson (1996) and Kurapov et al. (2017), the linear combination of low-frequency AP and low-frequency wind is the adjusted SL (ASL) in fact and is equal to the difference between low-frequency SL and IB correction.α0, α1 and α2, can be solved in a least squares sense by regressing the ASL onto the low-frequency wind component at the dominant wind orientation and low-frequency AP: When α0, α1 and α2 are determined, the response of the coastal ocean to the meteorological forcing can be estimated using Eq.
(2).As both the classical IB correction and regression analysis are implemented, this new methodology is labelled as IBR hereafter.

In-situ observations
In-situ observations were obtained from two stations, E1 and E2, in the Bohai Bay (Fig. 1).At each station, the total SL was measured using the moored pressure gauge, which was accurate to within 5 cm (Lv et al., 2018); in addition, the meteorological observations, including 10 m wind and sea level AP, were measured using the XZY3-type automatic observing system produced by the National Ocean Technology Center, SOA of China, of which the 10 m wind was measured by the XFY3-1 propeller anemograph that had been widely used in most coastal station systems in China (Wang et al., 2015).The XFY3-1 propeller anemograph was accurate to within 5° for wind direction and 0.5 m/s for wind speed (Lv et al., 2018).The hourly in-situ observations at E1 and E2 started at 0000 UTC 19 August and ended at 0000 UTC 18 November 2014.These meteorological observations were taken as the local forcing in this study.
Based on the IB correction (Paraso and Valle-Levinson, 1996), adjustment was made to the hourly total SL using the observed AP, and ASL was obtained.Similar to Lv et al. (2018), the ASLs were filtered using a cosine-lanczos filter (Duchon, 1979) with a high frequency cut-off of 0.8 cpd, through which the tidal, local inertial and higher frequency signals were eliminated, and the low-pass result was labelled as LASL.The wind components were similarly filtered, and the lowpass u wind, v wind and wind speed were labelled as LUW, LVW and LW, respectively.The similarly filtered total SL and AP were labelled as LSL and LAP, respectively.Similar to Sandstrom (1980), Hsueh and Romea (1983), Castro and Lee (1995) and Lv et al. (2018), all the hourly low-pass data, LASL, LUW, LVW, LW, LSL and LAP, were sub-sampled at 6hourly intervals, and the results were labelled as SLASL, SLUW, SLVW, SLW, SLSL and SLAP, respectively

ERA-Interim data
ERA-Interim data is the global atmospheric reanalysis data produced by the European Centre for Medium-Range Weather Forecasts (Dee et al., 2011).The ERA-Interim data is from 1 January 1979 to present and the minimum time interval is 6 hours.The details about the ERA-Interim data can be found in Dee et al. (2011).Six-hourly ERA-Interim analysis data of AP and 10 m wind, with a horizontal resolution of 0.125° latitude and longitude, were used to calculate the regional meteorological forcing in this study.
To evaluate the quality of the ERA-Interim data, the six-hourly AP and wind components in the ERA-Interim data and the sub-sampled in-situ observations, at E1 and E2, are compared in Fig. 3.It can be seen that the magnitudes and trends were similar to each other.The mean absolute errors (MAEs) and the correlation coefficients between the meteorological forcing in the sub-sampled in-situ observations and those in the ERA-Interim data are listed in Table 1.The MAEs for AP at both E1 and E2 were less than 0.85 mbar, indicating that AP in ERA-Interim data was close to the observed values in the insitu observations, as shown in Fig. 3a and Fig. 3d.The MAEs for u wind and v wind at E1 were about 1.6 m/s, but the correlation coefficients were larger than 0.8, indicating that both the u wind and v wind in the ERA-Interim data at E1 were in good agreement with the observed values in the in-situ observations, as shown in Fig. 3b and Fig. 3c.The same conclusions can be obtained at E2. Overall, the meteorological forcing data in ERA-interim data agreed well with the in-situ observations and can be used to calculate the regional forcing.
The regional meteorological forcing in the Bohai Bay, including AP and 10 m wind, was defined as the averaged value over the region spanning 117.5°-119°E and 37.9°-39.3°Nas shown by the area R1 in Fig. 1b.The regional meteorology was then similarly filtered with a high frequency cut-off of 0.8 cpd.The low-pass results were labelled as ERA-LAP and ERA-LW, respectively.The difference between the sub-sampled in-situ observed SL and the IB correction calculated using ERA-LAP was taken as the ERA-ASL, which was also filtered with a high frequency cut-off of 0.8 cpd, and the result was labelled as ERA-LASL.

Relationship between adjusted sea level and local wind
Following Lv et al. (2018), the magnitude squared coherences between SLASL and SLW were calculated, as functions of frequencies and wind orientation measured clockwise from north at an interval of 5 degrees, to find the relationship between the local ASL (i.e., SLASL) and local wind (i.e., SLW).The calculated results at E1 and E2 are shown in Fig. 4a and Fig. 4d, respectively.The corresponding averaged values of the magnitude squared coherences are shown in Fig. 4b and Fig. 4e for E1 and E2, respectively.It can be seen that the averaged magnitude squared coherences reached the maximum when the wind orientation was 80° at E1 and 10° at E2, demonstrating that the dominant wind orientations were 80° and 10° at E1 and E2, respectively, which were different from those at which the maximum magnitude squared coherences were obtained (Fig. 4a and Fig. 4d).The wavelet coherence and phase (Grinsted et al., 2004) between SLASL and SLW at the dominant wind orientation at E1 are shown in Fig. 4c, and those at E2 are shown in Fig. 4f.It was indicated that at E1, the SLASL and the SLW at 80° were antiphase in most of the periods at the 5% significance level, especially at 8-16 days.At E2, the lag time between SLW at 10° and SLASL was much larger than that at E1 (Fig. 4f).
As shown in Fig. 2, the SLSLs at E1 and E2 were almost equal, but the SLUW and SLVW at E1 were far from those at E2. Besides, the correlation coefficient between SLASL at E1 and the SLW at the dominant wind orientation (80°) was -0.56, while it was just -0.17 at E2.It was guessed that the SLASL at E2 may be more related to the SLW at E1 than that at E2.The magnitude squared coherences between SLASL at E2 and SLW at E1 were calculated and are shown in Fig. 4g; besides, the averaged values at each wind orientation are shown in Fig. 4h.The patterns of Fig. 4g and Fig. 4h are similar to Fig. 4a and Fig. 4b, respectively, indicating that the E2 SLASL response to E1 SLW was similar to the E1 SLASL response to E1 SLW.
The maximum value of the averaged magnitude squared coherences between SLASL at E2 and SLW at E1 was 0.41, which was much larger than that between SLASL at E2 and SLW at E2; moreover, the dominant wind orientation was 85°, which is similar to that shown in Fig. 4b; besides, the relationship between SLASL at E2 and SLW at E1 was negative and the lag time was much smaller than that between SLASL at E2 and SLW at E2 (Fig. 4i).It was hypothesized from the aforementioned results that the wind influencing the SLASL at E2 may be not the local wind, but the regional wind.

Relationship between adjusted sea level and regional wind
The magnitude squared coherences between the regional ASL (i.e., ERA-LASL) and regional wind (i.e., ERA-LW), at E1 and E2, were calculated and are shown in Fig. 5a and Fig. 5d, respectively.At both E1 and E2, the periods at which the coherences reached the maximum value were much less than those between SLASL and SLW, while the variations of the averaged values (Fig. 5b for E1 and Fig. 5e for E2) were similar to those shown in Fig. 4. Based on the averaged value of the magnitude squared coherences, the dominant wind orientation influencing ERA-LASL, at both E1 and E2, was 75°, further indicating that the regional wind affected the ASL in the Bohai Bay.Besides, at both E1 and E2, the ERA-LASL and ERA-LW at 75° were nearly antiphase in most of the periods at the 5% significance level (Fig. 5c and Fig. 5f), showing that the relationship between them was negative.
The correlation coefficient between the ERA-LASL at E1 and ERA-LW at 75° was -0.62, which was much larger than that between SLASL at E1 and SLW at 80° (-0.56).The correlation coefficient between SLASL at E2 and SLW at 10° was just -0.17 and that between SLASL at E2 and SLW at 85° at E1 was -0.54, which were less than that between ERA-LASL at E2 and ERA-LW at 75° (-0.58).The aforementioned results indicated that the ASL, at both E1 and E2, had much stronger relationship with the regional wind than local wind, providing guidance for estimating the response of the coastal ocean to the meteorological forcing in the Bohai Bay.

Response of the coastal ocean to meteorological forcing in the Bohai Bay
Based on the aforementioned conclusions, the regional meteorological forcing was the dominant driver of lowfrequency SL in the Bohai Bay, and will be used in the IBR to estimate the response of the coastal ocean to meteorological forcing.The ERA-LASLs at E1 and E2 were regressed onto the ERA-LW at 75° and ERA-LAP, as follows: The regression coefficients in Eq. ( 4) are listed in Table 2.It was noted that the MLR with Eq. ( 4) at both E1 and E2 were significant (p<0.001).It was apparent from the regression coefficients in Table 2 that ERA-LASLs at E1 and E2 were negatively correlated with both ERA-LAP and ERA-LW, as both α1 and α2 were negative; besides, the dynamic responses of the ocean at E1 to both ERA-LAP and ERA-LW at 75° were larger than those at E2, as α1 and α2 at E1 were larger than those at E2, which may be because the water depth at E1 was much smaller than that at E2, as shown in Fig. 1b.When ERA-LAP increases by 1 mbar, the SLSL will decrease by 1.94×10 -2 m at E1 and 1.78×10 -2 m at E2, with both the static ocean response and the dynamic ocean response considered, which were of the same order as those in Lv et al. (2018), suggesting that MLR in IBR was reasonable.
Combined the IB correction using ERA-LAP and the result of MLR with Eq. ( 4), the SLSL at E1 was estimated using Eq. ( 2) and is shown in Fig. 6a, while the result at E2 is shown in Fig. 6b.At both E1 and E2, the estimated SLSLs reproduced most of the local maxima of the observed values, although the observations were slightly underestimated in some cases; however, the estimated SLSLs deviated greatly from the observed values at most of the local minima.In addition, the SLSLs at both E1 and E2 were mainly attributed to the SL variabilities resulting from regional wind (i.e., ERA-LW), rather than AP (i.e., ERA-LAP), as shown in Fig. 6.The correlation coefficient and MAE between the observed and estimated SLSL at E1 (E2) were 0.70 (0.65) and 13.13 cm (13.44 cm), respectively, as listed in Table 3.Following MWRPRC (2014), Kurapov et al. (2017) and Lv et al. (2018), the frequency of occurrences when the estimated SLSL was within 0.15 m from the observed SLSL, labelled as FO, was taken as another metrics for the evaluation.FOs at E1 and E2 were 72.33% and 71.78%, respectively, as listed in Table 3.The results indicated that the estimated responses of the coastal ocean to the meteorological forcing, using IBR with regional forcing, could account for a significant fraction of the observed SLSLs.

Comparison with the classical IB correction
Traditionally, the IB correction was the classical method to estimate the response of the coastal ocean to the meteorological forcing (Wunsch and Stammer, 1997).As listed in Table 3, when the IB correction with local SLAP was used to estimate the response at E1 (E2), the MAE between the estimated SLSL and the observed values was 20.13 cm (19.28 cm), which was much larger than that when IBR with regional forcing was implemented; besides, FO was only 48.22% (47.40%), showing that the estimated SLSLs was far from the observed values on the whole, as shown in Fig. 7a (Fig. 8a).Comparison of the SLSLs at E1 (E2) estimated using the IB correction and IBR with regional forcing showed that IBR with regional forcing can improve FO by as much as additional 24.11% (24.38%).The aforementioned results indicated that the classical IB correction cannot accurately estimate the SLSL and the dynamic ocean response cannot be ignored.

Comparison with the dynamic atmospheric correction
As the ocean has dynamic response to meteorological forcing at high frequencies when considering large scales (Vinogradova et al., 2007), the DAC can significantly improve the altimetry product compared with the classical IB correction.To compare the estimated results using the IBR with regional forcing and the DAC in the Bohai Bay, MITgcm (Marshall et al., 1997) was used to simulate the dynamic response of the ocean to the meteorological forcing at E1 and E2.
The computing area was the whole Bohai Sea, as shown in Oregon State University Tidal Inversion Software (Egbert and Erofeeva, 2002).For the sake of simplicity, the model settings were similar to those in Exp1-Exp4, in which the model was two dimensional and the constant TS1 profile in Fig. 9 was used.
The horizontal resolution was 7.5′×7.5′ in Exp1-tide and 2′×2′ in Exp2-tide, respectively.The harmonic constants of each constituent, analysed from the simulated water level of last 30 days from 0000 UTC 1 August 2014, were compared with the observations at 14 tide gauge stations, whose locations are shown in Fig. 1b.Whether in Exp1-tide or Exp2-tide, the MAEs between the simulated and observed amplitudes of all the constituents were less than 10 cm, as shown in Fig. 10, while the MAEs for the phase lag were less than 16°, showing that the numerical model was validated and acceptable.
Following Carrè re and Lyard ( 2003) ,the residual signal variance and the ratio of the variation reduction compared with the IB correction, at E1 and E2, were calculated in all the numerical experiments and are listed in Table 5.When only AP was taken as the meteorological forcing, the residual signal variations at E1 (E2) were reduced less than 1.31% (1.84%), whatever the dimension, horizontal resolution and vertical stratification were assigned, showing that the DAC with only AP led to better estimated result than that using the IB correction, but the improvement was slight, similar to the conclusions in Hirose et al. (2001).On the contrary, the AP plus wind forcing reduced the variance at E1 (E2) by 23.12% (19.79%) to 24.57% (20.71%), further showing that the wind was an important driver of SLSL and cannot be ignored.when the simulated results in Exp8 and Exp10 were compared, it can be seen that the simulated response of sea level to AP and wind at E2 was improved when the topography was finer, but the results at E1 were not, which may be because the horizontal resolution of the ERA-Interim meteorological forcing was just 7.5′.Tierney et al. (2000) concluded that the density stratification did not make much difference in modelling SL high frequency signals (periods shorter than several weeks) because the response was essentially barotropic, which was also shown in Vinogradova et al. (2007).The results in Exp5 and Exp6 indicated that the near real density stratification TS2 did not improve the performance of DAC, when AP was taken as the meteorological forcing.Improvement was not observed either, when comparing Exp7 and Exp8 where both AP and wind were considered.However, when the density stratification TS1 was used, the residual variance can be further decreased if the model was changed to three dimensional from two dimensional in vertical direction, as can be seen when the results in Exp2 and Exp7 (Exp1 and Exp5, or Exp3 and Exp9) were compared.
Considering the sum of the ratio of the variation reduction at E1 and E2, the best performance of DAC was obtained in Exp7, in which both AP and wind were included in the forcing and the variation reduction reached 24.57% at E1 and 20.30% at E2 when compared with the IB correction, while the estimated results using the IBR with regional forcing induced a greater reduction of 42.27% at E1 and 36.79% at E2; besides, the IBR with regional forcing can improve FO by as much as additional 20.82% at E1 and 17.53% at E2, when compared to the results of DAC in Exp7, as listed in Table 3.The results indicated that the IBR with regional forcing had much better performance than the DAC, as shown in Fig. 7b and Fig. 8b.

Comparison with the multivariable linear regression
Based on the same observations used in this study, Lv et al. (2018) compared several regression methods and found that the best estimated results were obtained when the MLR involving wind component at dominant orientation and AP was used, as follows: where SAP and SW are the sub-sampled AP and wind component at the dominant wind orientation, respectively; γ0 is a constant, γ1 and γ2 are coefficients.
As listed in Table 3, the MAEs between the estimated SLSL in Lv et al. (2018) and the observed values were 14.51 cm at E1 and 16.03 cm at E2, which were much less than those when IB correction or DAC were used; besides, FOs at E1 and E2 were also larger than those with IB correction or DAC, showing that MLR with Eq. ( 5) was a useful method to estimate the response.
However, it was obvious that the estimated results using MLR with Eq. ( 5) at E1 and E2 were much farther from the observations than those using IBR with regional forcing, as shown in Table 3, Table 5,  some high frequency signals in the original wind, which will introduce significant harmonic motion into water bodies (Militello and Kraus, 2001), the estimated SLSLs using MLR with Eq. ( 5) included much more high frequency signals than those using IBR with low-pass regional forcing, at E1 (Fig. 7a) and E2 (Fig. 8a).Overall, the IBR with regional forcing had much better performance than MLR with Eq. ( 5) presented in Lv et al. (2018).

Comparison with the IBR using local meteorological forcing
For IBR with local meteorological forcing implemented, the MLR is as follows: The regression coefficients are listed in Table 2.When the regression coefficients were determined, the estimated SLSLs at E1 (E2) can be obtained using Eq. ( 2) and are shown in Fig. 7b (Fig. 8b).The estimated SLSL at E2 was slightly farther from the observed values than that using MLR with Eq. ( 5) (Table 3), as the correlation coefficient between SLASL and local wind at E2 was too small; however, the estimated SLSL at E1 was slightly closer to the observed values than that using MLR with Eq. ( 5), showing that the IBR with local forcing was not always worse than MLR with Eq. ( 5).As shown in Table 3, at both E1 and E2, the estimated SLSLs using IBR with regional meteorological forcing were much closer to the observed values than those using all the other methods, including IBR with the local forcing, demonstrating that IBR with regional forcing had the best skill in estimating the response of the coastal ocean to the meteorological forcing in the Bohai Bay.

Discussion about the response using IBR with regional meteorological forcing in the Bohai Bay
As the dominant wind orientation was 75° and the relationship between ERA-LASL and ERA-LW at 75° was negative, at both E1 and E2, the regional across-shore wind was an important factor influencing the ASL, which may be related with the fact that the Bohai Bay is a part of the Bohai Sea, a semi-enclosed coastal sea; in detail, the onshore wind will increase the ASL and the offshore wind will decrease the ASL, the same as in Lv et al. (2018) and different from continental shelves (Andres et al., 2013;Hsueh and Romea, 1983;Zhao and Cao, 1987).Besides, as shown in this study, the regional meteorological forcing was the dominant driver of SLSL in the Bohai Bay, rather than the local forcing, indicating that the dominant factors should be determined carefully at a specific location.
Although the estimated results using IBR with regional forcing were much closer to the observed response than those using other methods, including IB, DAC, MLR with Eq. ( 5) and IBR with local forcing, at both E1 (Fig. 7) and E2 (Fig. 8), it should be noted that the estimated results were far from the observations during some extreme events, as shown in Fig. 6. extreme at yearday 306 (316) due to swell swash (Melet et al., 2018).As shown in Fig. 12, before the SLSL reached the local maximum at yeardays 277.25 and 284, the wind was not strong; however, the wind in the Bohai Sea at yearday 284 was much stronger than 10 m/s, while the ERA-LW at 75° was just about -5 m/s, indicating that the wind set-up in the Bohai Sea may be the cause of this extreme event at E1 and E2 and a smaller lag time between SLSL and wind than that when the SLSL reached it local minima at yeardays 306 and 316.On the contrary, the wind during yearday 276.5 to 277.25 was not stronger in the Bohai Sea and north Yellow Sea, as shown in Fig. 12, so extreme SLSL at yearday 277.25 was not related with the wind and the cause was not clear.Overall, except the extreme maximum at yearday 277.25, the other three extreme events at yeardays 284, 306 and 316 were influenced by the remote wind that was not considered in the IBR, and therefore the estimated results using IBR with regional forcing during these extreme events deviated greatly from the observed SLSLs at both E1 and E2.The aforementioned results indicated that besides the regional meteorological forcing, there were other known and unknown factors influencing the SLSL in the Bohai Bay, which should be further investigated in the future.

Conclusions
The response of the coastal ocean to the meteorological forcing was a significant part of the total SL variations and cannot be estimated accurately using the traditional IB correction, as shown in Carrè re and Lyard ( 2003) and Lv et al. (2018).
Building on the description of the static ocean response and dynamic ocean response in Piecuch et al. (2016), the DAC in Carrere et al. (2016) and the MLR with Eq. ( 5) in Lv et al. (2018), a new methodology, named as IBR, was developed in this study to estimate the response of the coastal ocean to the meteorological forcing.In IBR, the response was a combination of the static ocean response and the dynamic ocean response.The former component was calculated using IB formula and the latter component was estimated using MLR with Eq. ( 3) involving low-pass AP and wind component at the dominant wind orientation.The observed SL from two stations, E1 and E2 located in the Bohai Bay, China, were taken as a case study to evaluate the feasibility and effectivity of IBR and compare it with other methods.
The magnitude squared coherences between SLASL and SLW were calculated and the averaged values at every orientation were taken as the indicator to find the dominant wind orientation.It was found that the correlation coefficient between SLASL at E2 and SLW at the dominant wind orientation at E1 (R=-0.54)was much larger than that between SLASL at E2 and SLW at the dominant wind orientation at E2 (R=-0.17),indicating that the SLASL at E2 was more related to SLW at E1 than to SLW at E2 and the response in the Bohai Bay was not attributed to the local meteorological forcing.
The AP and 10 m wind in the ERA-Interim data were spatially averaged in the Bohai Bay (R1 area in Fig. 1b) and the results were taken as the regional meteorological forcing.The ERA-LASLs, at both E1 and E2, were mainly forced by the ERA-LW at 75°; besides, the correlation coefficient between ERA-LASL at both E1 and E2 and ERA-LW at 75° was much larger than that between local SLASL and local wind.
Based on the regional meteorological forcing, including ERA-LAP and ERA-LW at 75°, the IBR was implemented to estimate the response of the ocean to the forcing at E1 and E2.The MAE between the estimated and observed response was 13.13 cm (13.44 cm) and FO was 72.33% (71.78%) at E1 (E2), indicating that the estimated response was much closer to the observations than those obtained using the other methods, including IB correction, DAC, MLR with Eq. ( 5) and IBR with local forcing, as shown in Table 3, Table 5 and Fig. 7 (Fig. 8).The aforementioned results demonstrated that the developed new methodology IBR was a feasible and relatively effective method to estimate the response of the coastal ocean to the meteorological forcing.Besides, most of the extreme events were influenced by the remote wind that was not considered in the IBR, so the estimated results using IBR with regional meteorological forcing deviated greatly from the observed values.The frequency of occurrences when the estimated SLSL within 0.15 m from the observed SLSL.
2 The best performance of DAC was obtained in Exp7 on the whole, so the estimated results in Exp7 were taken as the results 5 of DAC.

Fig. 1 .
The six-hourly AP and 10 m wind, extracted from the ERA-Interim data, were interpolated spatially and temporally to obtain the surface forcing.Flather radiation condition was used at the east open boundary, which is shown in Fig. 1b.As the meteorological forcing, model dimension, horizontal resolution and vertical stratification may affect the simulated results, several experiments (Exp1-Exp10) were carried out to discuss the factors.The detailed model settings of the numerical experiments are listed in Table 4.In all the above experiments, the model began with a cold start from 0000 UTC 1 June 2014, and continued running with a time step of 60s until 0000 UTC 20 November 2014.The initial 78 days were used for spin-up.For validation of the numerical model, ocean tides, which were similar to atmospheric forced ocean waves (Carrè re and Lyard, 2003), were simulated using MITgcm with the similar model settings.Four dominating constituents, including M2, S2, K1 and O1, were implemented as tidal forcing at the east open boundary (Fig. 1b); the tidal information was extracted from Ocean Sci.Discuss., https://doi.org/10.5194/os-2019-32Manuscript under review for journal Ocean Sci. Discussion started: 17 May 2019 c Author(s) 2019.CC BY 4.0 License.Hirose et al. (2001) concluded that the fine topography was preferable, but they only used a two-dimensional barotropic ocean model.The same conclusion can be drawn from the comparison of the simulated results in Exp2 and Exp4.However, Fig. 7a and Fig. 8a.As there were Ocean Sci.Discuss., https://doi.org/10.5194/os-2019-32Manuscript under review for journal Ocean Sci. Discussion started: 17 May 2019 c Author(s) 2019.CC BY 4.0 License.
Spitz and Klinck (1998) pointed out that the response of the sea water to wind forcing in the Bay was complex, depending on local forcing and nonlocal forcing.As shown in Fig. 11, before the SLSL reached the local minimum at yearday 306 (316), the wind was much stronger than 10 m/s and in NW-SE direction in the north Yellow Sea and Bohai Strait from yearday 305.25 (315.25) to 305.5 (315.5),reducing the SLSL locally due to wind set-up and may make the SLSLs at E1 and E2 Ocean Sci.Discuss., https://doi.org/10.5194/os-2019-32Manuscript under review for journal Ocean Sci. Discussion started: 17 May 2019 c Author(s) 2019.CC BY 4.0 License.

Figure 1 .
Figure 1.(a) Map showing the general location of the Bohai Sea (rectangle with dotted lines).(b) Map showing the locations of observation stations (black stars), E1 and E2, in the Bohai Bay; the location of the tide gauge stations (red circles) and the east open boundary (magenta triangles); the area R1 (rectangle with dotted lines), where the ERAwind is spatially averaged; and bathymetry of the Bohai Sea (colors).5

Figure 3 .
Figure 3. (a) Time series of the sub-sampled in-situ observations of AP at E1 (red line) and the six-hourly AP at E1 in ERA-Interim data (blue line); (b) same as (a) but for u wind component; (c) same as (a) but for v wind component; (d-f) same as (a-c) but for those at E2.

Figure 6 .Figure 7 .
Figure 6.The observed values of SLSL (gray line), the estimated SLSL using IBR with the regional meteorological forcing (black line) and the components, at (a) E1 and (b) E2.The black dotted lines show the extreme events.

Figure 10 .
Figure 10.Comparison of simulated and observed amplitude for (a) M2, (b) S2, (c) K1 and (d) O1 in numerical experiment Exp1-tide, in which the horizontal resolution is 7.5′.(e-h) Same as (a-d), but for those in numerical experiment Exp2-tide, in which the horizontal resolution is 2′.

Table 4 . The detailed model settings of the numerical experiments, when the DAC was used
Ocean Sci.Discuss., https://doi.org/10.5194/os-2019-32Manuscript under review for journal Ocean Sci. Discussion started: 17 May 2019 c Author(s) 2019.CC BY 4.0 License.