Alongshore Propagation of Subtidal Sea Level Fluctuations Around the Korean Peninsula Over Varying Stratification and Shelf Topography

Understanding subtidal (a period band of 3–20 days) coastal sea level fluctuations is of increasing importance for the sustainable use of coastal areas under challenging conditions resulting from regional and global sea-level rise. The wind-forced coastal-trapped waves (CTW) theory often accounts for subtidal sea level fluctuations observed off global coasts with a considerable range of propagation speeds. Here, the propagation speeds of subtidal sea level fluctuations observed at seven coastal tide-gauge stations around the Korean Peninsula (KP) from 1997 to 2017 were compared with phase speeds modeled for 10 segments of realistic bottom topography around the KP and four seasons of density stratification using a wind-forced CTW model. Alongshore variations in the modeled phase speed (2.0–10.0 m s–1, increasing with the bottom depth, shelf width, and vertical density difference) were consistent with those of the observed propagation speed (4.7–18.9 m s–1), despite systematic underestimation by 50%, i.e., the mean speed of 6.0 vs. 11.8 m s–1. This underestimation is discussed considering subtidal sea level fluctuations of non-CTW origin that were not incorporated into the CTW model, such as (1) Barometric sea level response to atmospheric pressure disturbances, and (2) Upwelling/downwelling response to local alongshore wind stress. This study suggests that subtidal coastal sea level fluctuates around the KP within and beyond CTW dynamics.


INTRODUCTION
Understanding subtidal (defined as a period from 3 to 20 days) coastal sea level variability is important for effectively adapting to rising global and regional sea levels (Stammer et al., 2013). Subtidal coastal sea level variability is often explained by alongshore propagating coastal-trapped waves (CTWs), which play an important role in coastal ocean circulation and turbulent mixing on continental slopes (Hughes et al., 2019;Schlosser et al., 2019). The CTWs generated by alongshore wind forcing propagate thousands of kilometers away along the coast with changing phase speeds due to local density stratification and shelf bottom topography (Wang and Mooers, 1976;Huthnance, 1978;Brink, 1982;Battisti and Hickey, 1984;Hughes et al., 2019;Schlosser et al., 2019). The subtidal sea level fluctuations in many coastal areas related to stratification and shelf topography have been studied using a classical CTW model; for example, phase speeds ranging from 2.72 to 14.10 m s −1 off the South African coast varying with shelf widths of 10-100 km (Schumann and Brink, 1990) and phase speeds ranging from 1.83 to 4.43 m s −1 off the west coast of India varying with shelf widths of 50-80 km (Amol et al., 2012).
However, owing to the assumptions of classical CTW theory, namely a straight coastline or limited consideration for external forcing (e.g., only alongshore wind forcing is considered), the CTW model often results in sea level fluctuations that are significantly different from the observations (Cho et al., 2014;Brunner et al., 2019). Recent studies have pointed out that the classic CTW model is not universally applicable to most coastlines that do not meet the required assumptions of a straight coastline with similar shelf topography (Brunner et al., 2019). The coastal areas around the Korean Peninsula (KP) are one such region where the assumptions for the CTW model cannot be optimally satisfied (Figure 1). The regions off the coasts around the KP consist of contrasting shelf topographies between the west and south coasts with a shallow and gentle (<0.003) slope connected to the Yellow and East China Sea and the east coast with a deep and steep (>0.01) slope connected to the East Sea (Japan Sea). Despite some previous studies explaining the subtidal coastal sea level fluctuations around the KP based on the classical CTW model under varying conditions of stratification and shelf topography (Lee and Chung, 1982;Lyu et al., 2002;Cho et al., 2014), the reasons for the significant difference between the observed (6.1-15.6 m s −1 ) and modeled (4.6-10.3 m s −1 ) propagation speeds are still poorly understood (Cho et al., 2014). In contrast, the subtidal sea level fluctuations off the east coast of the KP are also related to non-isostatic sea level responses to atmospheric pressure disturbances over a period of 2-20 days (Lyu et al., 2002;Nam et al., 2004;Jung et al., 2008), as well as a local upwelling/downwelling response to alongshore wind forcing beyond CTW dynamics (Park and Nam, 2018). Thus, our understanding of subtidal sea level fluctuations around the KP is far from complete; for example, how well the subtidal sea level fluctuations can be explained using the classic CTW and other dynamics beyond CTWs are poorly understood.
This study aims to (1) Characterize alongshore propagations of subtidal sea level fluctuations observed around the KP, (2) (Figure 1). The TG sea level data spans 21 years, from January 1997 to December 2017. Sea level fluctuations on the west and south coasts (BS, YS, MP, and GS) are generally much larger due to a higher tidal range of up to several meters, as compared with that along the east coast (SC, MH, and PH). Hourly sea level atmospheric pressures observed near the TG stations were used to correct the TG sea levels for the effect of atmospheric pressure loading required to extract sea level fluctuations dynamically linked to subsurface pressure perturbation.
Sea surface wind stress data from Copernicus Marine Environment Monitoring Service (CMEMS) blended satelliteobserved wind products with a time interval of 6 h and a spatial resolution of 0.25 • × 0.25 • were used to simulate the wind-driven subtidal coastal sea level fluctuations. Such blended wind products have been widely used, but low performances for coastal areas compared to offshore areas have been reported (Bentamy, 2019). Moreover, satellite-observed wind products may have an issue with quality, often yielding significant differences from in situ measured winds, as also reported for the study region, e.g., near the east coast of the KP due to the orographic effect (Nam et al., 2005). Although recent satellite wind products such as the CMEMS blended wind provide reasonable patterns of spatio-temporal wind variations , the data were compared with the wind products of Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA2) with a time interval of 6 h and a spatial resolution of 0.50 • × 0.33 • (Gelaro et al., 2017). The CMEMS and MERRA2 wind stress data at 74 grid points (W1-W74) from January 1997 to December 2017 were selected to use considering the typical wavelength of mode-1 CTWs (100-1,000 km) and the locations of the seven TG stations (Figure 1 and Supplementary Table 1).
Bottom topography data gridded with 2-arc minute grids within 100 km of the coastline were used to characterize the shelf topography. To characterize the seasonal density stratifications around the KP, vertical profiles of the water temperature (T) and practical salinity (S) were obtained from ship-based hydrographic data provided by the South Korea National Institute of Fisheries Science (NIFS). The NIFS hydrographic data collected at 15 standard depths (0,10,20,30,50,75,100,125,150,200,250,300,350,400, and 500 m; red horizontal grid lines in Figure 2k) have been provided on a bi-monthly basis at nominal stations around the KP since 1961; the data collected at selected stations around the KP (offshore region 50-100 km distant from the coast; pink squares in Figure 1) from 1997 to 2017 were used here. To supplement the hydrographic data, particularly for the northern part of the east coast of the KP where NIFS data are not available (Figure 1), World Ocean Atlas 2018 (WOA18) data were used, which are monthly averaged annual vertical climatology (12 datasets) profiles of temperature and salinity, constructed with a spatial resolution of 0.25 • × 0.25 • . The 82 grids from the WOA18 data around the northern part of the east coast of the KP and an additional 48 grids from the WOA18 data matching the NIFS stations were selected (Figure 1).

Data Processing
The coastline was defined as a location where the water depth is ≤10 m, matching grids for wind stress data (wind grids) in an alongshore interval of 14-40 km. Local coastline directions were determined as the direction connecting two neighboring coastline locations. A local coordinate system was used to determine the alongshore (y) and cross-shore (x) directions, where the alongshore direction was defined as a direction rotated clockwise from due north by θ, which is the same as the local coastline direction by definition (Supplementary Table 1 and Figure 1). The sea surface wind stress was decomposed into alongshore (τ y ) and cross-shore (τ x ) components as τ y = τ x 0 sin(θ) + τ y 0 cos(θ) and τ x = τ x 0 cos(θ) − τ y 0 sin(θ) with positive up-shelf (e.g., poleward on the east coast, eastward on the south coast, and equatorward on the west coast) and positive seaward, respectively, where τ x 0 and τ y 0 denote the eastward and northward wind stresses, respectively. The term "remote, " as opposed to "local, " location was defined here as a location more distant than one wind stress or coastline grid away from the reference location. To the north of the northernmost TG station (SC) along the east coast, three grids for wind stress (W9, W18, and W32; squares in Figure 1), which well-represent regional wind patterns, were selected for comparisons of regional wind and observed TG sea levels (η obs ). The 6-hourly alongshore and cross-shore wind stresses (τ y and τ x ) from both CMEMS and MERRA2 wind products were linearly interpolated to yield hourly time-series matching η obs . The ocean tides were removed from η adj using the tidal level prediction program known as the Unified Tidal Analysis and Prediction (UTide) developed for the MATLAB software (Codiga, 2011). A band-pass filter (5th order Butterworth) with high and low cut-off frequencies of 1/3 and 1/20 cpd was applied to all time-series of TG sea levels (η obs ), atmospheric pressures (P), and alongshore and cross-shore wind stresses (τ y and τ x ). The observed and modeled sea levels were normalized, e.g.,η = η−µ σ , for comparison among the various sea levels, where µ, σ, andη are the mean, standard deviation, and normalized form of η, respectively.
To directly link the ocean dynamics (subsurface ocean pressure) to sea level, the sea level response to atmospheric pressure should be removed from the observed sea level. As fast movement of the atmospheric pressure patterns over a large (e.g., synoptic) spatial scale causes nearly homogeneous sea levels observed at neighboring TG stations and minimization of the phase differences, fast propagations of the subtidal sea level fluctuations may be derived unless the sea level response to atmospheric pressure can be effectively removed from the observation. In general, the sea level response to atmospheric pressure is isostatic and is known as an inverted barometer (η 1 ). For example, the sea level drops by 1 cm in response to an atmospheric pressure rise of 1 hPa as follows (Wunsch and Stammer, 1997): where ρ 0 , g, and P 0 are the reference density (=1,000 kg m −3 ), gravity acceleration, and the reference pressure (=1013.25 hPa), respectively. However, the sea level response to atmospheric pressure loading becomes non-isostatic for a period of 2-20 days owing to the Helmholtz-like resonance occurring in a semiclosed deep basin connected to the outside through narrow and shallow straits, as in the case of the East Sea (Lyu et al., 2002;Nam et al., 2004); the effective removal of the sea level response to atmospheric pressure is not possible using η IB . The nonisostatic sea level response (η a ) considering the realistic bottom topography can be expressed as a function of time t, following Inazu et al. (2006), as: where the amplitude of non-isostatic sea level response is modified using G, the gain departure from the isostatic case as a function of time departure t − t, where t is time lag from the isostatic sea level response to atmospheric pressure (where G = 1.0 cm hPa −1 and t = 0) for each TG station (Supplementary Table 2). Thus, the adjusted sea levels (η adj ) were calculated from η adj =η obs − η a at the corresponding TG FIGURE 2 | (A) Cross-sectional profiles of the shelf topography within the coastal area around the KP [Segs. 1-10; from (a-j)] and (B) vertical profiles of the squared buoyancy frequency (N 2 ) averaged over four seasons (spring, summer, fall, and winter; pink, red, sky blue, and blue colors, respectively) for 10 segments [Segs. 1-10; from (k-t)] at the upper 300 m. In (a-j), water depth (y-axis) ranges from the surface to 3,500 m for Segs. 1-6 and 400 m for Segs. 7-10, respectively. In (a,k), the cross-sectional profiles of shelf topography and vertical profiles of N 2 for the East Australian coast (Church et al., 1986) and west coast of Peru (Brink, 1982;Illig et al., 2018) are overlaid with green solid and dotted lines, respectively. In (k), the 15 standard depths of the NIFS hydrographic observations are denoted by gray horizontal lines. In (t), standard deviations from the seasonal-mean and alongshore-mean N 2 (seasonal std. and alongshore std.) are shown with black solid and dashed lines, respectively, and zoomed in at depths 50-100 m. stations to extract the sea level components dynamically linked to subsurface pressure.
To estimate the observed propagation speeds (OPS) of subtidal sea level fluctuations and the modeled phase speed (MPS) of the mode-1 CTW around the KP over four seasons (spring; MAM, summer; JJA, fall; SON, winter; DJF), the coastal area around the KP was divided into 10 segments (Segs. 1-10) based on the alongshore changes in the coastline curvature and shelf topography of Segs. 1-4 and the locations of the TG stations for Segs. 5-10 (Figure 1). Four coastal zones were defined based on the maximum water depth (North Zone: Segs. 1-4; East Zone: Segs. 5-7; South Zone: Segs. 8-9; and West Zone: Seg. 10). A total of 84 periods (21 years × 4 seasons) of 3-month-long timeseries of sea levels were constructed for each TG station. A few cross-shore profiles of the bottom topography were averaged to estimate a representative shelf topography across the segment (Figures 2a-j). The total water depth offshore (seaward) of the slope was assumed to be constant. The alongshore distances between two TG stations were calculated by the summation of the distance of coastline grids between the two TG stations ( Table 1). The mean time lag in each segment was calculated by averaging the time lags with the maximum correlation between η adj at two TG stations. The OPS in each segment was calculated by dividing the alongshore distances by the mean time lag ( Table 1).
The vertical profiles of NIFS T and S were quality controlled and quality assured following (Bushnell and Worthington, 2016). First, outliers exceeding the global ranges of T (from −2.5 to 40 • C) and S (from 2 to 41) were removed. Next, outliers exceeding three standard deviations from the monthly climatological mean of T (ranging from 0.2 to 33 • C) and S (ranging from 27.1 to 35.7 • C) at each station and depth were removed. Conservative temperature (T C ), absolute salinity (S A ), and the density (ρ) at the surface pressure (zero dbar) were calculated from NIFS T and S using the Gibbs Seawater Oceanographic Toolbox (McDougall and Barker, 2011). To prevent vertical density inversion, data with higher density at upper than lower depth levels were removed from the analysis. A vertical profile of the buoyancy frequency squared (N 2 ) at each station was estimated using: All N 2 profiles were vertically interpolated at intervals of 1 m and were averaged over each season and segment. The N 2 derived from WOA18 data (N 2 W ) were used as reference and compared with the N 2 derived from NIFS data (N 2 obs ) for the same season and segment (Segs. 5-10); the difference was within 3% and quantified as |1 − N 2 W /N 2 obs |. A few N 2 profiles averaged over each segment were used to represent N 2 profiles in the segment (Figures 2k-t).
To analyze the sea level response to (P − P 0 ), periods of strong and weak (τ y ) fluctuations were selected based on the variance averaged over the total period for the subtidal time scale using the MATLAB toolbox for wavelet analysis (Torrence and Compo, 1998). Periods of "weak" and "strong" fluctuations of atmospheric pressure were defined as the periods when the subtidal (P − P 0 ) variance was lower and higher than the temporal mean ( = 0.25), respectively. In addition, the time-series of (P − P 0 ), η obs , η obs − η IB , and η obs − η a (η adj ) were compared for 2-day-long periods for both weak (13-15 July 2007) and strong (22-24 March 2012) (P − P 0 ) fluctuations.

Wind-Forced CTW Model and Its Application
The cross-sectional modal structure (F), frictional decaying coefficient (a), wind coupling coefficient (b), and MPS (c) were computed using the wind-forced CTW model. For linear, hydrostatic, Boussinesq, f-plane, free-surface, and long wave approximations, the equation of motion can be solved by expanding the subsurface oceanic perturbation pressure, p, as a sum of modes. Here, the subsurface perturbation pressure corresponds to the pressure due to the adjusted (not observed) sea level (η adj ), i.e., p =ρ 0 gη adj = ρ 0 g(η obs − η a ). To incorporate the non-isostatic sea level response to atmospheric pressure in the region, Eq. (2) instead of Eq. (1) was used in this study. Each mode contribution is a multiplication of the cross-sectional modal structure (F) and time-and alongshore varying amplitude (φ) as follows: where z is the vertical coordinate, t is the time, and F n is the normalized F of the n-th mode, which is weighted by the amplitude of the n-th mode CTWs, such that φ n is as follows (Battisti and Hickey, 1984;Clarke and Van Gorder, 1986): where a nm is the frictional decay coefficient of the n-th mode to m-th mode, which is calculated as (Jordi et al., 2005;Cho et al., 2014): where h(x) is the cross-shore profile of water depth across the segment, f is the Coriolis parameter at each latitude, and r is the bottom friction coefficient (Clarke and Van Gorder, 1986;Jordi et al., 2005;Cho et al., 2014). In Eq. (5), the b n and c n are the wind coupling coefficient and MPS of n-th mode CTWs, respectively, where b n is calculated as follows (Jordi et al., 2005;Cho et al., 2014): The frictional decay coefficient a and wind coupling coefficient b represent the bottom stress of friction and surface wind stress associated with bottom topography and density stratification affecting the CTW propagation in a certain mode (Brink, 1982(Brink, , 1991Battisti and Hickey, 1984). The amplitude of mode-1 CTWs (φ 1 ) can be calculated, in an integrating form, as a function of the distance (y) from the origin (y = 0) and time (t) as follows (Battisti and Hickey, 1984;Clarke and Van Gorder, 1986): where the φ 1 at the origin (y = 0) was set to zero following previous studies on CTWs in the region (Cho et al., 2014; TABLE 1 | C z and C s of the observed propagation speeds (OPS, m s −1 ) in two coastal zones (East Zone, West and South zones collectively), for six segments (Segs. 5-10), and the C of the OPS for four seasons (spring, summer, fall, and winter) and six segments.

(Time lag in h)
Winter (December-February) 11.0 ± 1.6 (2.2) 18.9 ± 2.8 (2.5) 14.9 ± 1.9 (2.3) 7.6 ± 1.2 (5.2) 10.9 ± 2.5 (5.8) 8.6 ± 2.2 (4.3) Seasonal standard deviations The C z , C s , and C are denoted with their uncertainties. Time lag corresponding to the OPS in hours is denoted with numbers in parentheses. The seasonal standard deviations (m s −1 ) for the C of OPS are denoted with italic letters below the C of each segment. Park and Nam, 2018). This boundary condition is justified as the wind stresses at up-shelf grids around the origin was not large compared to those at down-shelf grids (e.g., W9, W18, and W32) so that inclusion of wind stresses at the up-shelf grids does not change the results significantly. The MATLAB code package "bigr.m" was used to calculate the CTW parameters (Brink, 2018). Multiple cross-shore profiles of bottom topography across the shelf within each segment and vertical N 2 profiles averaged over each segment area were used as the inputs for the CTW model to represent the topographic and stratification condition of the segment (Figure 2). Considering the assumption of alongshore uniformity for the CTW model, each segment was divided considering alongshore changes (in order of magnitude comparable to wavelength of the CTWs) in coastline direction, cross-shore shelf topography, and stratification. The applicability of the CTW theory assuming straight coastline and uniform along-shelf conditions was tested by comparing the results with observations. The f and r were fixed to constants of 8.803 × 10 −5 s −1 and 0.005 cm s −1 , respectively (Cho et al., 2014). Results for the MPS were not sensitive (<1% change) to the selection of r in a reasonable range (from 0.005 to 0.1). The grid size was set to 100 (cross-shore) by 30 (vertical), and the width of the domain (maximum x) was set to 200 km, where the water depth was set to be constant for x > 100 km (seaward of the slope) (Cho et al., 2014). Both OPS and MPS were averaged over each season, segment, and coastal zone (Tables 1, 2). The season-mean OPS and MPS ( C) were calculated based on data from the corresponding season over the total period of 21 years, namely 84 separate periods (21 years and four seasons). The segment-mean OPS and MPS (C s ) and coastal zone-mean speeds (C z ) were calculated based on data for corresponding seasons and segments, i.e., for 252 separate periods (21 years × 4 seasons × 3 segments) for each coastal zone. The mean speeds were obtained by weightedaveraging the OPS and MPS over the season, segment, and coastal zone, where the weights are derived from the maximum correlation coefficients among the TG sea level fluctuations. The C variability were quantified from interannual spreads using standard errors calculated as σ √ Y (σ: standard deviation of season-mean OPS and MPS in each year and Y is the number of samples set to 21). Similarly, variability in C s and C z was estimated based on the standard errors as σ √ Y [σ : standard deviation of the segment-mean or coastal zone-mean OPS and MPS in each segment or coastal zone and Y is the number of samples set to 84 (segment-mean) or 252 (coastal zonemean)]. The CTW model results for four seasons and 10 segments were compared in terms of four parameters, estimated separately from the model: (1) Total water depth (H i ); (2) Shelf width (W i ) or the distance between the coastal wall and shelf break in segment C s ; (3) Vertical density difference (D i,j ) between the lower boundary of the seasonal thermocline (seasonal thermocline set to 70 m) and the bottom; and (4) Depthaveraged buoyancy frequency squared (M i,j ) for segment i and season j. These parameters were estimated to quantify the bottom topography and density stratification used as the model inputs, and to discuss the model results in this region in comparison to those in other regions. Here, the shelf break was defined as the location closest to the coastal wall where the bottom slope exceeds the mean bottom slope over the segment. The TABLE 2 | C z and C s of the modeled phase speeds (MPS, c 1 : m s −1 ) of the first CTW mode (n = 1) in three coastal zones (North Zone, East Zone, West and South zones collectively) and for ten segments (Segs. 1-10), and the C of the MPS, frictional decay coefficients (a 11 : 10 −9 cm −1 ), and wind coupling coefficients (b 1 : s 1/2 cm −1/2 ) of the first CTW mode for four seasons (spring, summer, fall, and winter) and ten segments.

Coastal Zone
North Zone East Zone West and South Zones The C z and C s for all segments and C for Segs. 5-10 (East Zone, West and South zones) are shown with their uncertainties. The seasonal standard deviations (m s −1 ) of C for the MPS are denoted with italic letters below the C for each segment. Burger number (Bu) in the ten segments for four seasons. The frictional decay coefficients and Burger numbers are highlighted by underlining and italic letters, respectively. The parameters for Segs. 5-10 (East Zone, West and South zones) are highlighted by gray shading.
H i and W i are parameters characterizing shelf topography, whereas D i,j and M i,j are the stratification conditions, defined as: where N 2 i,j is the buoyancy frequency squared (N 2 ) averaged over segment i and season j. Note that D i,j was calculated below d th , i.e., the depth where the standard deviation of the alongshore-mean exceeds that of the season-mean N 2 (Figure 2t). The Burger number, , at the shelf break for each season and segment was estimated to diagnose the relative importance between stratification and bottom relief, e.g., CTWs are more baroclinic (barotropic) for Bu 1 (Bu 1).

Decomposition of Wind-Driven Coastal Sea Level
Subtidal coastal sea level variations modeled using the classic wind-forced CTW theory and local and remote wind forcing beyond CTW dynamics were assumed to be a linear sum of three components: where η 1 , η 2 , and η 3 denote sea level fluctuations due to (1) The alongshore propagation of wind-forced CTWs, (2) Local upwelling/downwelling response to local alongshore wind stress, τ y , and (3) Sea level set-up/set-down by cross-shore wind stress, τ x , respectively. The CTW dynamics underlying η 1 are associated with the balances among the pressure gradient, bottom frictional stress, and surface wind stress in the alongshore direction, e.g., where τ y s and τ y b are the alongshore surface wind stress and bottom frictional stress, respectively, and geostrophic balance in the cross-shore direction. The η 2 value was derived from dynamics distinctly different from the CTW dynamics where the alongshore balance between the alongshore Coriolis force (due to non-zero crossshore currents) and alongshore surface wind stress via crossshore Ekman transport, e.g., fu = 1 ρ ∂τ y s ∂z , yielding alongshore current v as a time-integrated form of alongshore wind stress, along with the thermal wind relationship in the cross-shore direction. Details on the dynamics deriving η 2 are available from supporting information in Park and Nam (2018). The main difference between the dynamics of η 1 and η 2 is whether to allow the pressure gradient in the alongshore direction (η 1 ) or not (η 2 ), e.g., uniform pressure in the alongshore direction. Unlike the geostrophic balance (thermal-wind relationship) in cross-shore directions in the other two dynamics, the last term, η 3 , is not derived from the cross-shore geostrophic balance but non-negligible cross-shore wind stress where coastal setup (set-down) is associated with the balance among the crossshore pressure gradient, cross-shore Coriolis force due to the alongshore current, and surface wind stress in the cross-shore direction, e.g., −fv = − 1 ρ ∂p ∂x + τ x s ρ (Liu and Weisberg, 2007). The η 3 dynamics yield subtidal fluctuations of barotropic and baroclinic currents and sea level set-up/set-down primarily by the cross-shore wind stress, and may play a role in certain conditions, as in Liu and Weisberg (2007). Alongshore propagations of the coastal set-up (set-down) sea level were assumed to be governed by the balance between the alongshore pressure gradient and bottom friction, following mode-1 CTW dynamics.
The η 1 was directly calculated from the output of the CTW model (mode-1 CTWs), i.e., the perturbation pressure, P 1 (0) = F 1 (0, z) φ 1 y, t , at the coastal wall (x = 0), as follows: The η 2 can be estimated using the alongshore current (v) driven by the alongshore wind stress, based on Park and Nam (2018), expressed in the following form: and using the geostrophic balance in the cross-shore direction: where δ s , x, and T d are the maximum thickness of the upper layer in the two-layered system, the first baroclinic Rossby radius of deformation, and the time scale of the exponentially decaying wind stress, respectively. Here, x was defined as 1 f 0 ∫ δ s Ndz, and T d was set to 24 h following Park and Nam (2018). Subtidal sea level fluctuations induced by alongshore propagations of both windforced mode-1 CTWs and local upwelling/downwelling response to alongshore winds are expressed as follows: The propagation speed of subtidal η 1,2 fluctuations (MPS by η 1,2 ) was calculated using the same method applied to the OPS (those of η adj ) described in the previous section as the speed can be departed from c 1 due to the effects of local upwelling/downwelling dynamics (η 2 ) in response to alongshore wind varying along the coast, although the η 2 itself does not propagate along the coast (no alongshore pressure gradient). The local set-up/set-down of sea level due to cross-shore wind stress (η 3 ) is defined as follows Liu and Weisberg (2007): where the cross-shore distance, X, for integrating the cross-shore wind stress was set to 100 km, matching the maximum distance within the segment; for example, between the coastal wall and offshore boundary. Then, η 3 was obtained in the same manner as Eq. (8) was applied to η 1 , from the assumption that the alongshore propagation speed was the same as c 1 , yielding: The results of the comparisons between the remote alongshore wind stress (τ y ) and local sea levels (η adj and η 1 ) and between the modeled (η 1 , η 2 , η 3 , and η mod ) and observed sea levels (η adj ) were quantified using cross-correlation coefficients (R) and linear regression coefficients: slope (β) and intercept (α).

Characteristics of Observed Sea Level Propagation
Subtidal η adj fluctuations were found to propagate in the downshelf direction from SC to GS [equatorward along the east coast of the East Zone (Segs. 5-7), westward along the south coast of the South Zone (Segs. 8-9), and poleward along the west coast of the West Zone (Seg. 10)], with speeds varying over seasons and segments around the KP (Figures 3, 4). For example, on April 4-7, 1999 (corresponding to spring), the subtidal η adj fluctuations observed in SC on April 4 gradually propagated to GS (dotted line in Figure 3A).    Figures 4A,B). The η adj propagated faster along the east coast (East Zone; Segs. 5-7) than along the west and south coasts (West and South zones; Segs. 8-10), yielding significantly higher C z for the OPS of 13.1 ± 0.6 m s −1 along the East Zone compared to that at 8.0 ± 0.5 m s −1 along the West and South zones ( Table 1). The C s of the OPS reached a maximum speed of 16.2 ± 1.2 m s −1 in Seg. 6, which was 41-46% higher than that in Seg. 5 (11.6 ± 0.8 m s −1 ) and Seg. 7 (11.3 ± 0.9 m s −1 ) and was statistically significant (Table 1 and Figure 5A). The C s in Seg. 8 (6.3 ± 0.5 m s −1 ) and 9 (8.3 ± 0.7 m s −1 ) was significantly lower than those in other segments, and the minimum C s among all of the segments occurred in Seg. 8. The C s in Segs. 5, 7, and 10 (12.3 ± 1.4 m s −1 ) were not significantly different from each other.

h in spring and 25 h in summer in 1999 (dashed lines in
Seasonal variations in the C of the OPS differed among the segments (Table 1 and Figure 5A). The C in spring was considerably higher in Segs. 5 (17.4 ± 1.5 m s −1 ) and 10 (17.8 ± 2.4 m s −1 ) by 57-98% and 42-164%, respectively, than the C in other seasons in each segment. In contrast, C in Segs. 6 (13.5 ± 1.4 m s −1 ) and 7 (7.6 ± 0.8 m s −1 ) showed the lowest speed in spring compared to those in other seasons. Seasonal variation in C in Seg. 8 was not significant. In Seg. 9, a higher C (by 79-137%) occurred in less stratified seasons (i.e., spring and winter) than in more stratified seasons (i.e., summer and fall).

Modeled Characteristics of Mode-1 CTW Propagation
A more baroclinic modal structure, F 1 , smaller frictional decaying coefficient, a 11 , smaller wind coupling coefficient, b 1 , and higher MPS, c 1 , were obtained along the east coast (East Zone) than along the west and south coasts of the KP (West and South zones), indicating significant alongshore variations in the CTW characteristics (Table 2 and Figure 5C). The C z of the MPS was higher in the East Zone (6.8 ± 0.6 m s −1 ) than in the West and South zones (4.1 ± 0.4 m s −1 ) in a pattern similar to that in the observations ( Table 2). The C s of the MPS was highest in Seg. 6 (9.1 ± 0.3 m s −1 ) and lowest in Seg. 8 (3.7 ± 0.7 m s −1 ), which was also consistent with the observations (Table 2 and Figure 5C). In contrast, the C s in Seg. 7 (6.6 ± 0.4 m s −1 ) was 43% higher than in Seg. 5 (4.6 ± 0.0 m s −1 ), unlike in the observations; this difference was statistically significant. The a 11 and b 1 were approximately 10-and 2-fold higher in the West and South zones than in the East Zone, respectively ( Table 2). The F 1 was nearly baroclinic in the East Zone, yielding a Bu that generally exceeded unity, and more barotropic at the West and South zones with a relatively small Bu, particularly in spring and winter (Table 2 and Figure 6).
Seasonally, F 1 became more baroclinic with a smaller frictional decaying coefficient, a 11 and higher MPS, c 1 in summer and fall than in spring and winter, particularly in the West and South zones (Table 2 and Figure 5C). In contrast, wind coupling coefficient, b 1 was nearly constant (within 10%) over the seasons in the West and South zones. In the West and South zones, c 1 in summer and fall was 60-120% higher than that in spring and winter, indicating faster CTW propagation in the stratified seasons (Table 2 and Figure 5C). The seasonal standard deviations of c 1 were larger in the West and South zones (1.1-1.6 m s −1 ) than in the East Zone (<0.9 m s −1 ), implying that along the southern and western coasts of the KP, the CTWs were more sensitive to the seasonal variations of stratification ( Table 2). The a 11 was also more sensitive to seasonal changes in stratification in the West and South zones (1.29-8.49 × 10 −9 cm −1 ), where a 11 in summer and fall (1.29-4.49 × 10 −9 cm −1 ) was half of that in spring and winter (3.84-8.49 × 10 −9 cm −1 ), as compared to the East Zone (0.22-0.50 × 10 −9 cm −1 within 16% variation; Table 2). The F 1 in the West and South zones demonstrated a more baroclinic structure in summer and fall with a relatively large Bu (0.02-0.57), whereas in spring and winter, the F 1 demonstrated a more barotropic structure with a small Bu (<0.04); however, such seasonal variation was not clear in the East Zone, which was nearly baroclinic in all seasons (Figure 6).
The modeled sea level in association with the mode-1 CTWs (η 1 ) propagated in the down-shelf direction from SC to GS around the KP. The modeled sea level η 1 had the highest correlation with the remote alongshore wind, τ y , at W18, where the variance in the subtidal τ y fluctuation was the largest among the three grids for wind stress (W9, W18, and W32; Table 3 and Figure 7). The RMS difference in alongshore wind at the three grids between CMEMS and MERRA2 for the total period from January 1997 to December 2017 is 3.4 m s −1 , which is less than the standard deviations of 4.8 m s −1 (CMEMS) and 5.8 m s −1 (MERRA2) for the same period. The correlation patterns between the modeled sea level η 1 at the TG stations and remote alongshore wind τ y from the CMEMS wind products at the three grids were similar to those from the MERRA2 products ( Table 3). Regardless of kinds of wind products, the modeled sea level was significantly correlated with remote alongshore winds. For example, the equatorward τ y (downwelling favorable wind stress) at W18 on 20 and 24 September and 2 October led to drops in η 1 propagating from SC to GS around the KP (Figures 7A,B). From September to October 1999, η 1 at each TG station showed maximum correlation coefficients with τ y , with the greatest among the three grids at W18 (0.47-0.87). The τ y at W9 had a higher correlation coefficient (0.21-0.66) with η 1 at each TG station than τ y at W32 (<0.22), despite the farther distance (Table 3).

Comparison Between Observed and Modeled Sea Level Propagations
The CTW model systematically underestimated the propagation speed (i.e., MPS < OPS; Figures 5A,C). However, despite this systematic bias, common features were found in both OPS and MPS (c 1 ): C z of OPS and MPS was consistently higher in the East Zone (Segs. 5-7) than in the West and South zones (Segs. 8-10). Overall, except for Seg. 10, the alongshore structures of C s were consistent between the OPS and MPS. For example, both the OPS and MPS were the highest in Seg. 6, the second highest in Seg. 7, and the lowest in Seg. 8. However, the relatively high OPS in Seg. 10, which was similar to Segs. 5 and 7 in the East Zone, was not simulated by the CTW model. In addition, the C s of OPS (6.3-16.2 m/s) was double that of the MPS (3.7-9.1 m/s) in all segments.
The seasonal variations in C of the MPS (c 1 ) showed higher speeds in summer and fall than in spring and winter, but were not clear in the C of the OPS (Figures 5A,C). Most seasonal variations in the C of the OPS and MPS were not consistent, except for the common low speed in spring in Segs. 6-7. In particular, the higher C of the OPS in spring (compared to other seasons) in Segs. 5 and 10 was not found in the MPS. The seasonal standard deviations TABLE 3 | Cross-correlation coefficients (Left: CMEMS satellite wind, Right in italics: MERRA2 reanalysis wind) between τ y at three wind grids (τ y at W9, W18, and W32) and η 1 , and between τ y at three wind grids and η adj at each TG station from September 1 to October 30, 1999. Time lags in hours corresponding to the maximum correlation coefficients between τ y and η 1 , and between τ y and η adj are denoted by underlined numbers in parentheses. Regression coefficients [e.g., the slope (β), intercept (α), and R] between η adj and η 1 , η 2 , η 3 , and η mod ( = η 1 + η 2 + η 3 ) in the same period are denoted below. Note that the bold numbers indicate significance at the 95% confidence level (p < 0.05).

Stations
Frontiers in Marine Science | www.frontiersin.org FIGURE 7 | Time-series of (A) τ y (colored solid lines for CMEMS blended satellite wind and colored dashed lines for the MERRA2 reanalysis wind) at three wind grids (W9, W18, and W32) located north of SC (colored squares in Figure 1), (B) η adj (solid line) and η 1 (dotted line) at the 7 TG stations 1 September to 30 October 1999. In (A,B), the propagations of η adj and η 1 , as well as the related τ y described in section "Comparison Between Observed and Modeled Sea Level Propagations," are highlighted by colored arrows. The time-series of η adj and η mod ( = η 1 + η 2 + η 3 ) are shown at (C) MH and (D) GS. In the upper panel of (C,D), the η adj and η mod are represented by black solid and black-dashed lines, respectively. The η 2 (red-dashed line) and η 1 (red solid line), which are the most dominant at each station, were overlapped in the upper panels of (C,D), respectively. In the lower panels of (C,D), the red solid, red-dashed, and blue solid lines denote η 1 , η 2 , and η 3 , respectively. of the C of the OPS (1.3-4.3 m/s) were substantially larger than those of the MPS (<1.6 m/s; Tables 1, 2).
Subtidal fluctuations in η 1 in each segment reconstructed by the mode-1 CTWs showed a significant correlation (0.39-0.69, p < 0.013) with those of η adj observed at the corresponding TGs (Table 3). However, some sea level rises or drops in η adj were not consistent with those in η 1 , indicating that η adj could not be solely explained by the mode-1 CTWs (Figure 7). For example, in both η adj and η 1 , the sea level rise on 20-22 September, which was led by downwelling favorable (equatorward) τ y at W18, propagated sequentially from SC to GS (red arrows in Figures 7A,B). The sea level rise on October 23-24, which was led by downwelling favorable τ y at W18 simulated by the model (η 1 ), was not observed in η adj (blue arrows in Figures 7A,B). The sea level rise on September 24-26 in η 1 , led by the τ y at W18, was not observed at SC and MH, but at downstream TG stations from PH to GS (green arrows in Figures 7A,B). In addition, the subtidal fluctuations in η 1 were more correlated with η adj at PH-GS (R between η 1 and η adj : 0.511-0.688) than with those at SC (R between η 1 and η adj : 0.406) and MH (R between η 1 and η adj : 0.389).

Effects of Mode-1 CTWs on Alongshore Propagation of η adj
The shelf topography in each segment around the KP, which is characterized by water depth H and shelf width W, affected the propagation speed of the subtidal coastal sea level fluctuations, η adj , around the KP. We observed higher speeds for wider W and deeper H, as is common in many coastal areas, e.g., areas off the west coast of the United States (Battisti and Hickey, 1984), west coast of Peru (Brink, 1982), and the South African coast (Schumann and Brink, 1990). The effects of the water depth on the MPS (c 1 ) were confirmed from the results, i.e., increasing with an increasing H and decreasing bottom frictional decay (smaller a 11 ) (Table 1 and Figures 5C, 8A). Interestingly, the OPS had a similar pattern to the MPS at higher propagation speeds at deeper depths (Tables 1, 2 and Figures 5A, 8A). In Segs. 5-7 (East Zone), where H was deeper than Segs. 8-10 (West and South zones), the overall OPS and MPS were relatively high (Tables 1, 2 and Figures 5A,C, 8A). The effects of the shelf width became clearer when comparing Seg. 5 (narrower shelf) and 6 (wider shelf) for a high H, and between Seg. 8 (narrower shelf) and 9 (wider shelf) for a low H. The c 1 increased with W, accounting for the higher OPS and MPS (c 1 ) in Seg. 6 and 9 when compared to Seg. 5 and 8, respectively (Tables 1, 2 and Figures 5A,C,  8B). In summary, despite the systematic underestimation (up to 60%) in the CTW model, the alongshore variations in the OPS were explained by the CTW dynamics affected by the shelf topography (H and W). However, since the OPS was also affected by the other kind of wave propagation, e.g., fast barotropic Kelvin waves yielding much higher propagation speed  (Battisti and Hickey, 1984), South African coast (blue) (Schumann and Brink, 1990;Illig et al., 2018), west coast of Peru (cyan) (Brink, 1982;Illig et al., 2018), north-western Mediterranean Sea coast (green) (Jordi et al., 2005), and Eastern Australian coast (Church et al., 1986), are shown for comparison.
or nearly synchronous sea level response to atmospheric forcing (2nd event in Figure 3B), signal mixing with such processes beyond the CTW dynamics that is not incorporated into the CTW model may cause the model's systematic underestimation (higher OPS than MPS). Density stratification, characterized by the vertical mean buoyancy frequency squared, M, and the density difference between the bottom and the lower boundary of the seasonal thermocline (d th , 70 m), D, affects the MPS following mode-1 CTW dynamics; in other words, higher speeds with the stronger density stratification with larger D and M, but their effect on OPS was partial. This faster propagation of the CTWs by stronger seasonal stratification has been simulated in many studies, but has not yet been well observed (Wang and Mooers, 1976;Battisti and Hickey, 1984;Schumann and Brink, 1990;Brunner et al., 2019). Increasing stratification in summer and fall with the faster propagation of the CTWs is related to decreasing a 11 , because the increased stratification mitigates the effect of bottom friction and wave damping (Brink, 1982(Brink, , 1991. Although higher MPSs (c 1 ) were found with more stratified conditions, yielding a large D and M and small a 11 , as suggested by previous studies (Wang and Mooers, 1976;Brink, 1991;Brunner et al., 2019), the greatest seasonal variations in the OPS could not be explained by the CTW dynamics (Tables 1, 2 and Figures 5A,C). However, the low C of the OPS in Segs. 6 and 7 in spring, consistent with the MPS, suggested that the seasonal variations in the OPS could also be partly affected by seasonal stratification; in other words, the smaller the D and M, the lower the OPS (Figures 5A, 8C,D). In contrast, the effect of vertical density differences between the bottom and d th became clearer based on the comparison between the C s in Segs. 5 and 7. The higher or similar C s values of the OPS and MPS in Seg. 7, in comparison with Seg. 5, were due to the larger D in Seg. 7, despite the deeper water in Seg. 5 (Tables 1, 2 and Figures 5A,C, 8D). However, although a higher c 1 could be mainly explained by the stronger stratification with a larger M in Segs. 8-10 in summer and fall, the seasonal variations in the C of the OPS remain unclear. In summary, the CTW dynamics affected by the density stratification (D and M) partly accounted for both seasonal and alongshore variations in the OPS around the KP.
The MPSs around the KP reported in this study were not aligned with those found in other coastal regions, which is possibly due to distinctly different regimes in the shelf topography and density stratification, as compared with those of other studies. Relatively shallow water depth (small H) in the study area compared to those off the west coast of Peru (Brink, 1982;Illig et al., 2018) and the East Australian coast (Church et al., 1986; Figure 9) with no significant difference (slightly narrower) in shelf width W from those in other coastal areas (Figure 2a), cannot explain the faster propagating CTWs in the study region (with small H and W) than other regions (Figure 9). Note that the MPS generally decreased toward the left-bottom corner in Figure 9A, except for those found in this study. Thus, a striking difference in the stratification regime is responsible for the abnormal MPS in the study region. As the seawater that occupies below a depth of ∼100 m in the study region is cold (water temperature is less than 1 • C) and nearly homogeneous , the N 2 is extremely small, as confirmed by the comparison of the vertical N 2 profiles in Seg. 1 to those off the west coast of Peru (Brink, 1982;Illig et al., 2018) and the East Australian coast (Church et al., 1986); however, strong seasonal variation in N 2 in the upper 100 m allows even higher stratification in the study region (Figure 2k). The stratification characterized by N 2 averaged over the upper 100 m yielded the distinctly different regimes from the other regions, accounting for the high MPS in the study region (Figures 9B,C).
The correlation between η adj and τ y was higher at W18 than at W9 or W32, which demonstrates the importance of remote wind forcing and alongshore/temporal scales, relevant to the mode-1 CTWs, on the subtidal sea level fluctuations propagating around the KP (Table 3 and Figures 7A,B). Subtidal sea level fluctuations associated with remote wind forcing at locations separated by hundreds of kilometers, rather than local winds, following the dynamics of mode-1 wind-forced CTW, have been reported in this (Cho et al., 2014;Park and Nam, 2018) and other coastal areas (Battisti and Hickey, 1984;Clarke and Van Gorder, 1986). The subtidal sea level fluctuations along the entire coastline of the KP (1,830 km from Segs. 1 to 10) can be affected by alongshore wind forcing off the north coast (W9 and W18) following the dynamics of the mode-1 wind-forced CTW, considering the sufficiently long frictional damping length scales (Csanady, 1978;Park and Nam, 2018). The frictional damping length scales corresponding to frictional decay coefficients of 0.33-5.46 × 10 −9 cm −1 for the North Zone ( Table 2) is 1.8-30.3 × 10 3 km, supporting the sensitivity of subtidal sea level fluctuations along the entire KP coast to alongshore wind off the north coast. The signal of subtidal sea level fluctuations generated by the alongshore wind forcing off the north coast would arrive in the west coast within a few days following the propagation speed of mode-1 CTWs. For example, for the distance from W18 to GS (1,330 km) and the C s of c 1 in Segs. 3-10 (3.7-9.1 m s −1 ), the sea level fluctuation led by τ y at W18 arrived at GS within 1.7-4.2 days, which was confirmed

Effects of Wind and Atmospheric Forcing on Alongshore Propagation of η adj
The subtidal fluctuations in the atmospheric pressure loading to the sea level over a large area with a minimized time lag could induce higher OPS beyond CTW dynamics, resulting in a systematically higher OPS than MPS. The subtidal fluctuations in the atmospheric pressure pattern around the KP are known to be mainly induced by the eastward movement of low pressure developed inland in Asia in spring and summer (Jung et al., 2008). As shown in the period of strong fluctuations in atmospheric pressure (March 22-24, 2012; Figures 10B-E), with variance higher ( = 0.50) than the mean ( = 0.25; Figure 10A), the eastward propagation of sea level fluctuations (from MP to PH; Figure 10C), induced by the pattern of eastward movement of subtidal fluctuations in atmospheric pressure (Figure 10B), remained stable after two different corrections (Figures 10D,E). In contrast, in the period of weak fluctuations in atmospheric pressure (July 13-15, 2007; Figures 10F-I), with variance lower than the mean ( = 0.15; Figure 10A), the sea level fluctuations induced by the pattern of subtidal fluctuations in atmospheric pressure ( Figure 10G) were removed, and the westward propagation of the sea level relevant to the propagation of the CTWs became distinct after corrections (Figures 10H,I). The C s of the OPS based on η adj , chosen only in periods of weak fluctuations in atmospheric pressure (OPS-LP), was smaller by 15 ± 13% than the C s of the OPS based on the total η adj (Figures 5A,B, 11A). In addition, the significantly higher C of the OPS at Segs. 5 and 10 in spring, which could not be accounted for by the MPS, disappeared in the C of the OPS-LP ( Figure 5B). When the subtidal fluctuations in atmospheric pressure were strong, their loading effects on sea level would not optimally be removed by the corrections, such that they would increase the correlations between the sea level at adjacent TG stations by a similar pressure loading effect, yielding a larger OPS than the theoretical CTW phase speed.
The upwelling/downwelling response of η adj to the local τ y could cause the systematically higher OPS (than MPS) by inducing a higher propagation speed with a minimized time lag due to the sufficiently large spatial (alongshore) scale of subtidal τ y fluctuations (longer than the distance between neighboring TGs and comparable to the wavelength of mode-1 CTWs). The sea levels between neighboring TG stations responded almost simultaneously to the local τ y ; thus, inducing minimized time lags and higher propagation speeds. When considering both the upwelling/downwelling response to the local τ y and mode-1 CTWs (MPS by η 1,2 ), the C s of the MPS was 82 ± 50% larger than that of the MPS predicted by only η 1 (black-dashed and solid lines in Figure 11A). The upwelling/downwelling response of the sea level to the local τ y was more significant in Segs. 5-8 than in Segs. 9-10, yielding a larger difference between the MPS by η 1,2 and MPS by η 1 in Segs. 5-8 due to the wider alongshore scales of the local τ y , as confirmed by the longer decorrelation length scales of τ y in Segs. 5-8 ( Figure 11B).
Subtidal sea level fluctuations around the KP primarily induced by mode-1 CTWs (η 1 ) could be modified by the , the causative forcing factors affecting the subtidal sea level fluctuation discussed in this study are denoted with gray (P, atmospheric pressure loading) and sky-blue arrows (τ y and τ x , alongshore and cross-shore wind forcing). The alongshore variations in the observed propagation speed of the subtidal sea level fluctuations are also represented with red arrows around the coastline of KP. In the right panels, the conditions for (B) the high MPS and (C) the low MPS following mode-1 CTWs dynamics are shown, i.e., water depth (H), shelf width (W), and the vertical density difference (D) between the lower boundary of the seasonal thermocline (d th ) and the bottom.
upwelling/downwelling response to local alongshore wind forcing (η 2 ) or local set-up/set-down of sea level due to crossshore wind forcing (η 3 ), but the impacts varied among TG stations. The predominant η 2 was in SC-MH, the predominant η 1 was in BS-GS, and η 3 made the least contribution at all TG stations. The η 1 made the largest contribution to η adj in BS-GS, as confirmed by the largest slopes ranging from 0.192 to 0.584 among all factors, as well as larger correlation coefficients between η 1 and η adj in PH-GS when compared to those in SC-MH (Table 3 and Figures 7C,D). In contrast, η 2 had the greatest contribution to η adj in SC-MH, with slopes ranging from 0.207 to 0.357 (Table 3 and Figure 7C). This result is consistent with the fact that the upwelling/downwelling response to local alongshore wind primarily explains the alongshore current variability at time scales shorter than approximately 16 days in the mid-east coast of South Korea (Park and Nam, 2018). The individual η 3 had the least contribution to η adj (slope in the linear regression β of −0.162 to 0.140) at all TG stations among all components, along with a very low R (0.00-0.33); however, the sum of all of the components (η mod ), including η 3 , better accounted for the subtidal sea level fluctuations in SC-MH (slope β of 0.251-0.682) than did any individual components, yielding a larger R (0.442-0.447) than any individual components (Table 3 and Figures 7C,D). The results reveal that the upwelling/downwelling response to local wind forcing plays a significant role in shaping subtidal sea level fluctuations around the KP beyond the wind-forced CTW dynamics where the time-varying wind forcing is considered with constant wind-coupling coefficient b.

CONCLUSION
Alongshore propagation in the down-shelf direction from SC to GS (equatorward along the east coast, westward along the south coast, and poleward along the west coast) around the KP of the coastal subtidal (3-20 days) sea level fluctuations were examined and modeled for mode-1 CTWs under varying shelf topographies and stratifications. In particular, alongshore and seasonal variations in the propagation speed were presented and discussed in terms of the dynamic response of the coastal sea level to atmospheric pressure and wind forcing ( Figure 12A). The observed subtidal sea level fluctuations propagated with a higher speed of 13.1 ± 0.6 m s −1 along the east coast (Segs. 5-7) and lower speed of 8.0 ± 0.5 m s −1 along the west and south coasts (Segs. 8-10). This was similar to the alongshore variations of mode-1 CTWs (6.8 ± 0.6 m s −1 along the east coast and 4.1 ± 0.4 m s −1 along the west and south coasts), despite systematic underestimation by the CTW model that may be linked to dynamics beyond the wind-forced CTWs. The phase speeds of mode-1 CTWs increased with a deeper water depth, wider shelf width, and larger vertical density difference; these also accounted for most alongshore variations in the observed subtidal sea level propagation speeds around the KP (Figures 12B,C).
The barometric sea level response to the rapidly moving synoptic weather system and upwelling/downwelling response to local alongshore wind could account for the systematically higher propagation speeds of the observed sea level fluctuations. Subtidal sea level fluctuations around the KP were still affected by atmospheric pressure loading at the coastal sea level, which would not be fully corrected by known methods that consider the non-isostatic sea level response to atmospheric pressure in the semi-enclosed basin. Moreover, subtidal sea level fluctuations off the east coast beyond the CTW dynamics are partly explained by the upwelling/downwelling response to local (not remote) alongshore wind forcing, rather than the set-up/set-down response due to cross-shore wind forcing. The sufficiently larger alongshore scales of both the atmospheric pressure and local alongshore wind relative to distances among TG stations induced almost simultaneous sea level fluctuations at neighboring TG stations, thereby decreasing the time lags and deriving systematically higher propagation speeds. The atmospheric pressure loading, which would not be fully removed, may induce a 15% higher propagation speed, and the upwelling/downwelling response of coastal sea level to local alongshore wind may induce an 82% higher propagation speed for subtidal sea level fluctuations than that predicted by the CTW model. In addition to the previous works about the missing points of the classical CTW theory, including the complex coastline that may impact scattering of the CTWs (Huthnance et al., 1986;Jordi et al., 2005;Brunner et al., 2019), and non-linear advection that may hamper the model prediction (Cho et al., 2014), this study suggests the new point of subtidal coastal sea level fluctuations, particularly within semi-enclosed marginal seas: (1) Nonisostatic sea level response to atmospheric pressure loading and (2) Upwelling/downwelling response to local wind forcing play an important role in shaping the subtidal sea level fluctuations beyond the CTW theory.

DATA AVAILABILITY STATEMENT
The TG station data were provided by the Korea Hydrographic and Oceanographic Administration (KHOA) and are available at https://www.khoa.go.kr/koofs. The shipbased hydrographic data were provided by the National Institute of Fisheries Science (NIFS) and are available at https://www.nifs.go.kr/kodc/. The sea surface pressure data were provided by the Korea Meteorological Administration (KMA) and are available at https://data.kma.go.kr/cmmn/main.do. The bottom topography data were from ETOPO2 (Earth topography 2 arc-minute) data provided by the US National Geophysical Data Center and are available at https://www.ngdc.noaa.gov/mgg/global/etopo2.html.
The monthly climatology data were provided by WOA18 (World Ocean Atlas 2018) and are available at https: //www.nodc.noaa.gov/OC5/woa18/. The blended wind stress data were provided by the Copernicus Marine Environment Monitoring Service (CMEMS) and are available at https: //marine.copernicus.eu/services-portfolio/access-to-products. The reanalysis wind stress data were provided from the Modern-Era Retrospective Analysis, version-2 (MERRA-2) and are available at https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/.

AUTHOR CONTRIBUTIONS
KL and SN contributed to conception of the study. KL organized the methodology, software, performed the formal analysis, visualized all of the figures, and wrote the original draft of the manuscript. KL, SN, and J-HP investigated the formally analyzed data and contributed to writing-review and editing. SN supervised all processes. All authors contributed to the manuscript in multiple ways.

FUNDING
This study was also a part of the project titled, "Deep Water Circulation and Material Cycling in the East Sea (20160040), " funded by the Ministry of Oceans and Fisheries (MOF), South Korea. This study was funded by the KHOA, MOF, through a project named titled, "Analysis and Prediction of Sea Level Change in Response to Climate Change around the Korean Peninsula."