Numerical Modeling of Submarine Pipeline Scouring under Tropical Storms

: Submarine pipelines are the lifelines of the national economy. Under the inﬂuence of typhoons, high-speed currents and waves continuously erode the seabed, leading to suspension or even rupture of pipelines. Therefore, it is of great importance to study the sediment transport under the action of waves and currents. A numerical model of sediment scouring and deposition combining wave and currents is established, which considered tidal current, storm surges, wind waves, and sediments in the East China Sea. Combining with the monitoring of the actual laying condition of pipelines, it is found that the area with the most serious scouring is around KP300. It is shown that the typhoon weather with high intensity and density will lead to the suspension of pipelines, which is noteworthy in the construction of marine engineering.


Introduction
Oil and gas pipelines are regarded as the lifelines of marine energy transmission. Since the first submarine pipeline was laid in the Mexico Gulf by the United States in 1954, the bottom of the sea has seen vigorous development of submarine pipelines all over the world [1]. Nowadays, more than 3000 km of submarine pipelines have been laid in China. Along the coast of Zhejiang Province, submarine pipelines of groups of Chunxiao Gas Field, Pinghu Oil and Gas Field and Hangzhou Bay Cross-sea Oil have been constructed [2].
The pipeline of Pinghu Oil and Gas Field landed in Zhoushan City, Zhejiang Province where two pipelines broke on 2000 (Chen et al. 2005) [3], leading to oil spills. As a result, hundreds of millions of dollars were spent, and serious environmental damage was caused by widespread oil pollution. Therefore, it is essential to ensure the safety of submarine pipelines. Figure 1 shows how accidents happened due to the local scour at pipelines. In 1982, six European leading gas storage and transport companies launched a campaign to collect data on accidents in pipeline systems, setting up the European Pipeline Accident Data Organization (EGIG), which generated statistics about the causes of the failure of submarine gas pipelines between 1970 and 2010 [4]. Now, EGIG is a co-operation the southeast of Shanghai city and 350 km to the south of Ningbo city, Zhejiang province. In 2007, after Typhoon Lily, pipeline suspension was detected near KP300, thus section KP287.665~KP301.906 was selected as the digital-analog area. Several exposed and suspended sections were found in the offshore section of Chunxiao gas field group subsea pipeline, which seriously threatens the safety of the pipeline.
such as complex terrain, geomorphologic region and hydrodynamic environment, and can also provide reference for the monitoring of submarine pipelines.

Study Area
As is shown in Figure 2, the Chunxiao Gas Field Group (28°10′ N~28°48′ N, 124°45′ E~125°15′ E) is located in the East China Sea continental shelf, which is about 450 km to the southeast of Shanghai city and 350 km to the south of Ningbo city, Zhejiang province. In 2007, after Typhoon Lily, pipeline suspension was detected near KP300, thus section KP287.665~KP301.906 was selected as the digital-analog area. Several exposed and suspended sections were found in the offshore section of Chunxiao gas field group subsea pipeline, which seriously threatens the safety of the pipeline.

Topography
The entire seafloor where the submarine pipeline is located gradually inclines in the SE direction with a maximum water depth of 110 m and an average depth of 90~110 m. The topography of the seabed in the whole area is basically floating without severe undulations or steep ridges, whose average slope is 0.3 × 10 −3 . There are tidal channels, underwater shore slopes, comb-shaped sand ridges and the central submarine plain areas of the continental shelf from the land to the sea.
The local topography has a great impact on the hydrodynamic environment. In this research, the KP287.66~KP301.906 section of the submarine pipeline is located on the outer

Topography
The entire seafloor where the submarine pipeline is located gradually inclines in the SE direction with a maximum water depth of 110 m and an average depth of 90~110 m. The topography of the seabed in the whole area is basically floating without severe undulations or steep ridges, whose average slope is 0.3 × 10 −3 . There are tidal channels, underwater shore slopes, comb-shaped sand ridges and the central submarine plain areas of the continental shelf from the land to the sea.
The local topography has a great impact on the hydrodynamic environment. In this research, the KP287.66~KP301.906 section of the submarine pipeline is located on the outer edge of the Mountain Niubi water channel, which belongs to the underwater bank slope section where water depth is between 11 and 15 m. Besides, the sedimentary landform is the main feature in this area, and the seabed sediments are mainly deposited by the coastal rivers since the Holocene.
In this paper, the 1:250000 coastline data from National Geophysical Data Center of National Oceanic and Atmospheric Administration (NOAA) is selected (https: //www.ngdc.noaa.gov/mgg/shorelines/shorelines.html accessed on 17 May 2021), whose resolution and precision can meet the requirements.
As is shown in Figure 3, the high-resolution water depth data accumulated by the Second Institute of Oceanography of the State Oceanic Administration for many years were used in the area with complex coastal topography, and 30" × 30" bathymetric data provided by the UK Ocean Data Center were used in offshore areas. (https://www.bodc.ac.uk/ data/hosted_data_systems/gebco_gridded_bathymetry_data/ accessed on 17 May 2021). whose resolution and precision can meet the requirements.
As is shown in Figure 3, the high-resolution water depth data accumulated by the Second Institute of Oceanography of the State Oceanic Administration for many years were used in the area with complex coastal topography, and 30" × 30" bathymetric data provided by the UK Ocean Data Center were used in offshore areas. (https://www.bodc.ac.uk/data/hosted_data_systems/gebco_gridded_bathymetry_data/ accessed on 17 May 2021).

Seafloor Sediments
The subsoil in the study area can be roughly divided into clayey soil, silty soil and sandy soil, whose spatial distribution is shown in Figure 4. The seabed is mainly clay soil from KP287.66 to KP301.906. The offshore sea is sandy soil, and the silts mainly appear as interlayers in the sea area near KP150. In 2007, the sea area with pipeline suspension was monitored. The sediment was mainly composed of silt and clay, of which clay was slightly dominant, accompanied by a small amount of shell fragments, and the sediment particle size was between 5 and 75 μm.

Seafloor Sediments
The subsoil in the study area can be roughly divided into clayey soil, silty soil and sandy soil, whose spatial distribution is shown in Figure 4. The seabed is mainly clay soil from KP287.66 to KP301.906. The offshore sea is sandy soil, and the silts mainly appear as interlayers in the sea area near KP150. In 2007, the sea area with pipeline suspension was monitored. The sediment was mainly composed of silt and clay, of which clay was slightly dominant, accompanied by a small amount of shell fragments, and the sediment particle size was between 5 and 75 µm.

Marine Hydrology
According to measured data, the tidal current in the pipeline area is typical irregular semi-diurnal. The measured maximum velocity of the bottom layer was 96 cm/s and the residual velocity was 2.2 cm/s. Moreover, the flow rate of high tide is slightly higher than that in low tide. The movement of the tidal current is typical reciprocating flow near the shore, showing a rotating flow in the outer sea area (outside the Mountain Niubi water channel). Generally, the coastal waters of Zhejiang are mainly wind-generated waves which are more sensitive to typhoons and not controlled by the incoming surges outside the wind-zone.

Wind Field
Statistics in recent decades show that the average number of tropical cyclones affecting the region is about four times per year. The years that were most affected by tropical cyclones were 1959 and 1960, respectively, while the year that was least affected was 1996.

Marine Hydrology
According to measured data, the tidal current in the pipeline area is typical irregular semi-diurnal. The measured maximum velocity of the bottom layer was 96 cm/s and the residual velocity was 2.2 cm/s. Moreover, the flow rate of high tide is slightly higher than that in low tide. The movement of the tidal current is typical reciprocating flow near the shore, showing a rotating flow in the outer sea area (outside the Mountain Niubi water channel). Generally, the coastal waters of Zhejiang are mainly wind-generated waves which are more sensitive to typhoons and not controlled by the incoming surges outside the wind-zone.

Wind Field
Statistics in recent decades show that the average number of tropical cyclones affecting the region is about four times per year. The years that were most affected by tropical cyclones were 1959 and 1960, respectively, while the year that was least affected was 1996. Tropical cyclones generally occur from May to November, which are mainly concentrated from July to September, accounting for about 80% of the total. August is the peak of tropical cyclone activity that accounts for 35% of the total.

Storm Surge Model
Under the assumption of incompressible liquid and shallow water, the Navier Stokes Equations were solved ignoring vertical acceleration. Orthogonal curvilinear grid was adopted for a horizontal plane (Arakawa C-type grid system is used for equation variable arrangement), while σ grid was applied for vertical direction which can fit the terrain. K-ε turbulence model was selected, and finite difference method was used to discrete the equation.
The vertical average mass conservation equation is as follows.
∂ζ ∂t Q represents the role of source and sink, such as drainage and precipitation. Momentum conservation equation in a direction ζ: Momentum conservation equation in a direction η: In σ coordinate, the vertical velocity component is solved by mass conservation equation: Here, ω refers to the vertical velocity perpendicular to the σ coordinate plane, which changes as the σ coordinate plane moves up and down. In the above formulas, H is the water depth, H = h + ζ is the water level, and is the water depth relative to the mean sea level; G ξξ and G ηη are the conversion coefficients of curvilinear coordinate system to rectangular coordinate system; U and V are the average velocity in ξ and η direction respectively. g is the acceleration of gravity; f is the parameter of Coriolis force; F ξ and F η are turbulent momentum fluxes in ξ and η directions, respectively. P ξ and P η represent the water pressure gradient in ξ and η directions; M ξ and M η represent the source or sink of momentum in the ξ and η directions, respectively; q in and q out represent source and sink items.
ADI scheme is a common method for solving parabolic and elliptic equations in a finite difference method. It is improved from crank Nicolson scheme. It has second-order accuracy in time and space and is unconditionally stable. It divides a time step into two steps and transforms a five diagonal linear equation system into a tridiagonal linear equation system, which greatly improves the computational efficiency and is widely used. In the horizontal direction, the model is discretized by the ADI method.
Step 1: Step 2: Among them: where λ is the linear bottom friction coefficient, → d is the external force, such as wind and atmospheric pressure. In the first step, the calculation process is from l to l + 1 2 , that is, the calculating time is from t = l∆t to t = l + 1 2 ∆t. In this process, V Equation (3), U Equation (2) and continuity Equation (1) coupled by free surface gradient are solved; In the second step, the calculation process is from l + 1 2 to l + 1, that is, the calculation time is from l + 1 2 ∆t to (l + 1)∆t. In this process, U Equation (2), V Equation (3) and continuity Equation (1) coupled by free surface gradient are solved.

Storm Wave Model
In SWAN, the wave spectrum is described by the spectral action balance equation given in Hesselmann [20]. The SWAN model does not use two-dimensional spectral density, but rather uses the two-dimensional dynamic spectral density to represent random waves. In the flow field, the dynamic spectral density is conserved, but the energy spectral density is not. The dynamic spectral density N(σ, θ) is the ratio of the energy spectral density E(σ, θ) to the relative frequency σ. In a Cartesian coordinate system, the dynamic spectrum balance equation is expressed as follows: In the formula, the first term on the left ∂ ∂t N is the change of wave action with time, the second term is ∂ ∂x c x N, the third term ∂ ∂y c y N is the propagation of wave action in geographical space, the fourth term ∂ ∂σ c σ N is Doppler shifting caused by topography and current, and the fifth term ∂ ∂θ c θ N is refraction caused by topography and current, Is the source term, including wind-induced wave, dissipation, nonlinear wave action and wave breaking. They are the propagation velocity in the direction. S(= S(σ, θ)) is the source term, including wind-induced wave, dissipation, nonlinear wave action and wave breaking. c x , c y , c σ , c θ are propagation velocities in (x, y, σ, θ) direction. Various physical processes Water 2021, 13, 1425 7 of 18 considered in the model are input through the right source term, including wind energy input, energy dissipation and nonlinear wave interaction, the specific forms are as follows.
Wind energy input: The input of wind energy is described by resonance mechanism feedback mechanism in SWAN model: The coefficients are determined by wave frequency and direction as well as wind speed and direction.
Energy dissipation: The dissipation of wave energy is mainly composed of three parts, including white wave dissipation S ds,w (σ, θ), bottom friction dissipation S ds,b (σ, θ) and wave breaking dissipation S ds,br (σ, θ).
The white wave is mainly controlled by the wave steepness. In the current third generation wave model, the white wave dissipation adopts the Hsselmann formula [20].
where Γ is a coefficient related to wave steepness, σ and k are average wave frequency and average wave number. The dissipation caused by topography is mainly caused by bottom friction, seafloor infiltration and reflection scattering of irregular topography. In the continental shelf of a sandy seabed, seabed dissipation is mainly caused by bottom friction, which can be expressed as follows: where C bottom is the bottom friction coefficient. For the dissipation caused by wave breaking, swan adopts the following formula: For the dissipation caused by wave breaking, the following formula is used in SWAN: E tot is the total energy and D tot < 0 is the total energy dissipation rate caused by wave breaking. It mainly depends on the wave breaking factor γ = H max , which can be a constant or a variable in SWAN.
Nonlinear wave interaction: Nonlinear wave interaction plays a very important role in wave deformation. In deep water, the fourth-order wave interaction determines the evolution of spectral frequency, which can transform wave energy from peak frequency to low-frequency and highfrequency. In shallow water, third-order wave interaction can transform wave energy from low-frequency to high-frequency to generate high-order harmonics. In the SWAN model, the discrete interaction approximation (DIA) and the centralized third-order approximation (LTA) are used to approximate. Figure 5 shows the disturbed area of numerical calculation results in the model calculation domain under wind and wave conditions. The rectangle represents the calculation water area, the thick solid line represents the wave front boundary with known boundary conditions, the thin solid line represents the lateral boundary of the unknown wave, and the shaded area represents the disturbed area of the model result due to the assumption. In the case of wind wave, the disturbed region propagates from the top of the known wave boundary conditions to the shore in the range of 30-45 • with the average wave direction. For the swell wave, the disturbed area will be smaller due to its smaller directional distribution range. A typical Typhoon Lily in 2007 was selected to calculate the typhoon wave. Its path is shown in Figure 6. The same water depth and shoreline data are used as that in calculation of storm surge, and the calculation method of wind field is the same. The verification of calculation results is a difficult problem in wave numerical simulation. Generally, there are three methods: the first is to verify with the measured data, the second is the data collected by the satellite altimeter, and the third is the data provided by other recognized and reliable numerical wave prediction models. Obviously, it is the most ideal if it can be verified by the measured data, but generally there is little chance to obtain the wave data of fixed area and fixed time period. For satellite observation, the reliability is greatly lower than the measured data, and it may not be able to capture the data we just A typical Typhoon Lily in 2007 was selected to calculate the typhoon wave. Its path is shown in Figure 6. The same water depth and shoreline data are used as that in calculation of storm surge, and the calculation method of wind field is the same. A typical Typhoon Lily in 2007 was selected to calculate the typhoon wave. Its path is shown in Figure 6. The same water depth and shoreline data are used as that in calculation of storm surge, and the calculation method of wind field is the same. The verification of calculation results is a difficult problem in wave numerical simulation. Generally, there are three methods: the first is to verify with the measured data, the second is the data collected by the satellite altimeter, and the third is the data provided by other recognized and reliable numerical wave prediction models. Obviously, it is the most ideal if it can be verified by the measured data, but generally there is little chance to obtain the wave data of fixed area and fixed time period. For satellite observation, the reliability is greatly lower than the measured data, and it may not be able to capture the data we just The verification of calculation results is a difficult problem in wave numerical simulation. Generally, there are three methods: the first is to verify with the measured data, the second is the data collected by the satellite altimeter, and the third is the data provided by other recognized and reliable numerical wave prediction models. Obviously, it is the most ideal if it can be verified by the measured data, but generally there is little chance to obtain the wave data of fixed area and fixed time period. For satellite observation, the reliability is greatly lower than the measured data, and it may not be able to capture the data we just need for the satellite may not be running in the wave calculation time. For the data Water 2021, 13, 1425 9 of 18 provided by other models, the biggest drawback is the reliability of the data, especially the nearshore terrain.
Due to the limitation of conditions, the data provided by the global wave prediction model WAVEWATCH III of NOAA is used for verification. The resolution of the model is 1 • × 1.25 • and there are special calculations for typhoons in some regions such as the Gulf of Mexico. The resolution of the model is relatively high because it is a targeted calculation, and its accuracy is relatively high. Unfortunately, in the East China Sea, this kind of study is not provided. Figure 7 shows the distribution of global ocean effective wave height predicted by WAVEWATCH III at 14:00 Beijing time on 15 July 2007. Near the East China Sea (the area in the red box), the annular wave field caused by typhoon Lily is clearly discernible. Effective wave height refers to the actual wave height calculated according to certain rules. Since sea waves are actually a random combination of waves with different wave heights, periods and directions, the wave height value of a wave is not representative. For this reason, in any wave group composed of n waves, the wave heights in the wave train are arranged in order from large to small, and the first 1/3 waves are determined as effective waves. The wave height and period of the effective wave are equal to the average wave height and average period of N over three waves [19]. need for the satellite may not be running in the wave calculation time. For the data provided by other models, the biggest drawback is the reliability of the data, especially the nearshore terrain. Due to the limitation of conditions, the data provided by the global wave prediction model WAVEWATCH III of NOAA is used for verification. The resolution of the model is 1° × 1.25° and there are special calculations for typhoons in some regions such as the Gulf of Mexico. The resolution of the model is relatively high because it is a targeted calculation, and its accuracy is relatively high. Unfortunately, in the East China Sea, this kind of study is not provided. Figure 7 shows the distribution of global ocean effective wave height predicted by WAVEWATCH III at 14:00 Beijing time on 15 July 2007. Near the East China Sea (the area in the red box), the annular wave field caused by typhoon Lily is clearly discernible. Effective wave height refers to the actual wave height calculated according to certain rules. Since sea waves are actually a random combination of waves with different wave heights, periods and directions, the wave height value of a wave is not representative. For this reason, in any wave group composed of n waves, the wave heights in the wave train are arranged in order from large to small, and the first 1/3 waves are determined as effective waves. The wave height and period of the effective wave are equal to the average wave height and average period of N over three waves [19]. In the calculation area of this paper, three effective points are selected for verification, and the verification results are shown in Figure 8. Although the two values are not in good agreement, the trend of significant wave height change is very consistent, and the peak values of three points appear at the same time. It shows that the model can reflect the response of waves in Zhejiang coastal area when typhoon passes by. The difference of the significant wave height between the two is caused by many reasons, some of which are even unavoidable. First of all, although they belong to the third generation of numerical wave models, there are differences in theory, especially in the nearshore, the processing of various physical processes is not the same, in addition, the input of wind energy is not consistent; secondly, the three boundaries of SWAN model in this paper cannot provide real boundary conditions, which will certainly affect the calculation accuracy; finally, the WAVEWATCH of NOAA in terms of the accuracy of typhoon analysis, is lower than that of the global wind field. In the calculation area of this paper, three effective points are selected for verification, and the verification results are shown in Figure 8. Although the two values are not in good agreement, the trend of significant wave height change is very consistent, and the peak values of three points appear at the same time. It shows that the model can reflect the response of waves in Zhejiang coastal area when typhoon passes by. The difference of the significant wave height between the two is caused by many reasons, some of which are even unavoidable. First of all, although they belong to the third generation of numerical wave models, there are differences in theory, especially in the nearshore, the processing of various physical processes is not the same, in addition, the input of wind energy is not consistent; secondly, the three boundaries of SWAN model in this paper cannot provide real boundary conditions, which will certainly affect the calculation accuracy; finally, the WAVEWATCH of NOAA in terms of the accuracy of typhoon analysis, is lower than that of the global wind field. Comparing the significant wave heights at different locations, it is obvious that the points in the open sea (with greater longitude) have larger peaks, and the points in the South have larger peaks than those in the north. It can be seen from Figure 8 that the two calculation results are consistent.
Although the two lines of effective height do not correspond well, they are very consistent in terms of the trend of effective wave heights, and the peaks of the three points all appear at the same time. For this research, the ultimate purpose is not to study the platform wind and wave itself, but to consider the influence of waves on the seabed shear stress and scour effect. Therefore, such a result is acceptable. This shows that the model can reflect the wave response in the coastal area of Zhejiang when a typhoon passes by. In addition, the difference of effective wave heights between the two is also caused by many reasons, some of which are even unavoidable. First of all, although both belong to the third-generation numerical wave model, there are theoretical differences, especially in the nearshore, the processing of various physical processes is not the same, and the input of wind energy is not consistent. Secondly, none of the three boundaries of the SWAN model in this paper can provide real boundary conditions, which will certainly affect the calculation accuracy. Finally, NOAA's WAVEWATCH Ⅲ forecast value is for the whole world, with a relatively low resolution. Due to the topography and other factors, the accuracy of the near shore is relatively poor. For the regions affected by the typhoon field, there is no detailed wind field, which is bound to affect the calculation accuracy. Figure 9 shows the distribution of significant wave height at four typical moments in the calculation area under Typhoon Lily. The four stages are: (a) when the center is far away from the calculation area, the wind field has little influence on the region, so the wave is mainly controlled by the open sea boundary; (b) when the typhoon is close to the calculation area, the regional wave is obviously affected by the wind field, thus the effective wave height of the region affected by strong wind began to increase significantly; (c) after the effect of the wind field on the region has lasted for appropriately 10 h, the wave is significantly increased, and the typhoon center is closest to the center of the region, which is also the time when the typhoon field has the strongest impact on the entire region. At this time, the wave of the whole region reaches its peak; (d) the typhoon leaves slowly to the north, the wind field begins to weaken, and the wave also weakens. Comparing the significant wave heights at different locations, it is obvious that the points in the open sea (with greater longitude) have larger peaks, and the points in the South have larger peaks than those in the north. It can be seen from Figure 8 that the two calculation results are consistent.
Although the two lines of effective height do not correspond well, they are very consistent in terms of the trend of effective wave heights, and the peaks of the three points all appear at the same time. For this research, the ultimate purpose is not to study the platform wind and wave itself, but to consider the influence of waves on the seabed shear stress and scour effect. Therefore, such a result is acceptable. This shows that the model can reflect the wave response in the coastal area of Zhejiang when a typhoon passes by. In addition, the difference of effective wave heights between the two is also caused by many reasons, some of which are even unavoidable. First of all, although both belong to the third-generation numerical wave model, there are theoretical differences, especially in the nearshore, the processing of various physical processes is not the same, and the input of wind energy is not consistent. Secondly, none of the three boundaries of the SWAN model in this paper can provide real boundary conditions, which will certainly affect the calculation accuracy. Finally, NOAA's WAVEWATCH III forecast value is for the whole world, with a relatively low resolution. Due to the topography and other factors, the accuracy of the near shore is relatively poor. For the regions affected by the typhoon field, there is no detailed wind field, which is bound to affect the calculation accuracy. Figure 9 shows the distribution of significant wave height at four typical moments in the calculation area under Typhoon Lily. The four stages are: (a) when the center is far away from the calculation area, the wind field has little influence on the region, so the wave is mainly controlled by the open sea boundary; (b) when the typhoon is close to the calculation area, the regional wave is obviously affected by the wind field, thus the effective wave height of the region affected by strong wind began to increase significantly; (c) after the effect of the wind field on the region has lasted for appropriately 10 h, the wave is significantly increased, and the typhoon center is closest to the center of the region, which is also the time when the typhoon field has the strongest impact on the entire region. At this time, the wave of the whole region reaches its peak; (d) the typhoon leaves slowly to the north, the wind field begins to weaken, and the wave also weakens. Water 2021, 13, x FOR PEER REVIEW 12 of 19 Figure 9. Distribution of effective wave height at four moments under Typhoon Lily.

Coupling Model
At present, the wave elements are considered in the most commonly used numerical model for calculating the interaction of two-dimensional waves and currents based on the N-S equation to deduce the mathematical model. This model is similar to the two-dimensional tidal current control equations, with radiation stress terms, where the wave current interaction is included.
However, due to the computational efficiency and complexity of the model, most of the models only consider the effect of unidirectional wave elements on the flow field and thus cannot consider the influence of the flow field on the wave field synchronously. Therefore, the wave current interaction is uncoupled or pseudo coupled, resulting in relatively lower simulation accuracy when the wave influence is large.
Driven by the boundary dynamic force, at first, the flow field calculation is carried out, outputting the elements affecting the wave to the wave calculation model, and then the wave field calculation of the first step is started. After that, the elements affecting the It can also be seen from Figure 9 that although the nearshore is not directly affected by a strong typhoon field, the effective wave height is about 3 m due to the wave conduction from the open sea. In Zhoushan Archipelago, due to the blocking effect of many islands on waves, the wave height in the Western sea area of the islands is obviously smaller.

Coupling Model
At present, the wave elements are considered in the most commonly used numerical model for calculating the interaction of two-dimensional waves and currents based on the N-S equation to deduce the mathematical model. This model is similar to the twodimensional tidal current control equations, with radiation stress terms, where the wave current interaction is included.
However, due to the computational efficiency and complexity of the model, most of the models only consider the effect of unidirectional wave elements on the flow field and thus cannot consider the influence of the flow field on the wave field synchronously. Therefore, the wave current interaction is uncoupled or pseudo coupled, resulting in relatively lower simulation accuracy when the wave influence is large.
Driven by the boundary dynamic force, at first, the flow field calculation is carried out, outputting the elements affecting the wave to the wave calculation model, and then the wave field calculation of the first step is started. After that, the elements affecting the flow field are delivered to the flow field model for the calculation of the next step. This cycle is repeated until the calculation reaches the specified time.

Scouring Model
A coupled large-area model is established with consideration of ocean dynamic elements under typhoon, such as the interaction of waves and currents. Meanwhile, the scouring and silting of the offshore areas with complicated topography during typhoon was investigated, which help understand the scouring and silting rules of the area where submarine pipelines are located.
The sediment transport module was used to simulate the transport of cohesive and non-cohesive sediment (such as the diffusion of dredged material) and to study the sediment deposition/erosion model.
The appropriate boundary conditions are used to solve the two-dimensional convectiondiffusion equation. In sediment calculation, the two-dimensional diffusion equation of suspended sediment is as follows: For the study area which is dominated by cohesive fine-grained sediment, full silting model was used to calculate the sediment erosion and deposition rate which were simulated by the Partheniades-Krone equation [18]: where E l represents the scour sediment flux; M l represents scour coefficient; D l represents deposition flux; w l s represents sediment settling velocity; c l b is average sediment concentration near bottom. S τ αv , τ l cr,e and S τ cw , τ l cr,d are respectively expressed in the following equation: where τ l cr,e represents the critical shear stress of sediment erosion, τ l cr,d denotes the critical shear stress of deposition, τ cv is the shear stress on the bed surface.
The bed sediment transportation will lead to the change of bed elevation which is directly determined by the sediment deposition flux and suspended flux, expressed as: where S x , S y are the sediment transport components, zb is the bed elevation, ε por is the bed porosity, which is generally 0.4. T d is the bed silting deposition or erosion rate, T d = D − E.

Results and Discussion
In local scour modeling, the calculation time was set from 8:00 September 10 to 23:00 16 September 2007, which fully included the whole time period of Typhoon Lily. Figure 10a shows the change of flow velocity over time around KP300. The velocity response is not very strong, which is related to a variety of factors, such as the intensity of the typhoon, the type of path, the relationship between its acting time and ebb and flow, and topographic factors [10]. However, the change of the wave is very obvious. During the period of typhoon transit, the effective wave height of the wave near KP300 increases obviously, and decreases with the distance of the typhoon. The extreme effective wave height reached near 3 m at around 18:00 on 15 September (Figure 10b).

Results and Discussion
In local scour modeling, the calculation time was set from 8:00 September 10 to 23:00 16 September 2007, which fully included the whole time period of Typhoon Lily. Figure  10a shows the change of flow velocity over time around KP300. The velocity response is not very strong, which is related to a variety of factors, such as the intensity of the typhoon, the type of path, the relationship between its acting time and ebb and flow, and topographic factors [10]. However, the change of the wave is very obvious. During the period of typhoon transit, the effective wave height of the wave near KP300 increases obviously, and decreases with the distance of the typhoon. The extreme effective wave height reached near 3 m at around 18:00 on 15 September (Figure 10b).
For the study area which is dominated by cohesive fine-grained sediment, the full silting model was used to calculate the sediment erosion and deposition rate which were simulated by the Partheniades-Krone equation [21]: The results of the local scour modeling are shown in Figures 11 and 12. As can be seen from Figure 11, four points (a, b, c, d) on the near-shore section of the pipeline route are selected to obtain their scour and deposition ranges during the typhoon. It can be intuitively seen from Figure 11 that an area centered on the location of the KP300 pipeline (shown by the red circle in Figure 11) experienced an obvious scouring process during the For the study area which is dominated by cohesive fine-grained sediment, the full silting model was used to calculate the sediment erosion and deposition rate which were simulated by the Partheniades-Krone equation [21]: The results of the local scour modeling are shown in Figures 11 and 12. As can be seen from Figure 11, four points (a, b, c, d) on the near-shore section of the pipeline route are selected to obtain their scour and deposition ranges during the typhoon. It can be intuitively seen from Figure 11 that an area centered on the location of the KP300 pipeline (shown by the red circle in Figure 11) experienced an obvious scouring process during the typhoon. However, in the area east of the pipeline route KP280, with the continuous increment of water depth, the impact of waves on the seabed bottom becomes smaller, and the topography becomes uncomplicated, so the flow structure becomes simple and the scour effect is not obvious.
Water 2021, 13, x FOR PEER REVIEW 15 typhoon. However, in the area east of the pipeline route KP280, with the continuou crement of water depth, the impact of waves on the seabed bottom becomes smaller, the topography becomes uncomplicated, so the flow structure becomes simple and scour effect is not obvious.   Figure 12 shows the scouring process at pipelines for four different points. The scour and deposition of points a, b, c and d can be summarized as silting, balance, scour and intense scour, showing a transitional trend. Explaining it in detail, point a is silted up during the typhoon, point b is basically in equilibrium, point c is scoured during the most violent period of the typhoon and then slowly silted up, and point d (near KP300) has a relatively violent scour process during the typhoon. After the typhoon, it begins to slowly accumulate so as to recover to the original state and reach equilibrium.
A strong current will be formed when the typhoon comes. As a result, the shear stress at the bed bottom will increase rapidly, resulting in a large number of sediment movement. Typhoons transfer energy to a water body and make it possible to scour violently in a short time by changing the shear stress of the bed and the sediment-carrying capacity of water. From point a to point d, these four observation areas are far away from the shore (Figure 11). The sediment cannot be supplied in time for areas such as point c and point d, which are relatively far away from the shore, resulting in severe regional scouring, exposure, and even suspension of pipelines. Restricted by the terrain, the flows and wind waves are roughly in the same direction. In areas (point c and point d) that are scoured during a typhoon, sediment from these areas is transported by the current and deposited in an adjacent area (point a) along the path of the typhoon. Meanwhile, in the middle area, a balance between scouring and deposition appears (point b). As a result, the original dynamic balance of scouring and silting was ruined which led to severe local scouring within a short period of time, resulting in critical suspension of pipelines and other dangerous conditions. Figure 11. Regional Erosion and Deposition after Typhoon Lily.   Figure 13 shows the shallow stratum profile of the suspended pipeline detected in 2007. The reflection wave of the suspended pipeline and the V-shaped trench scoured by the bottom current are clearly visible, and the V-shaped trench is asymmetric on both sides of the pipeline. The south slope is gentle, while the north slope is steep, which has to do with the direction of the bottom flow. From the scour pit shape, the scour is mainly caused by the south-to-north bottom flow, and a long gentle slope on the left side of the pipeline is formed under the influence of the wake flow. Seabed sediment near KP300 is unconsolidated clay soil, so the critical starting velocity is low, and the water depth is shallow. Therefore, the wind wave has a great influence on the bottom shear stress, providing conditions for the scouring of pipelines near KP300. Combined with mathematical model results and theoretical analysis, it is reasonable to believe that the main culprit behind submarine pipeline suspension danger in the Chunxiao gas field group was the typhoon weather with high density and intensity in 2007. pipeline is formed under the influence of the wake flow. Seabed sediment near KP300 is unconsolidated clay soil, so the critical starting velocity is low, and the water depth is shallow. Therefore, the wind wave has a great influence on the bottom shear stress, providing conditions for the scouring of pipelines near KP300. Combined with mathematical model results and theoretical analysis, it is reasonable to believe that the main culprit behind submarine pipeline suspension danger in the Chunxiao gas field group was the typhoon weather with high density and intensity in 2007. Under the influence of Typhoon Lily in 2007, the original dynamic equilibrium of scouring and deposition was ruined in the offshore section of submarine pipeline in the Chunxiao Gas Field Group, which led to severe local scouring in a short time, resulting in critical suspension and a suspension of the pipeline in this section. Energy was transferred Under the influence of Typhoon Lily in 2007, the original dynamic equilibrium of scouring and deposition was ruined in the offshore section of submarine pipeline in the Chunxiao Gas Field Group, which led to severe local scouring in a short time, resulting in critical suspension and a suspension of the pipeline in this section. Energy was transferred from the typhoons to the water body by changing the bed shear stress and sedimentcarrying capacity of water body, making it possible to scour intensely in a short time.
According to the monitoring data for many years, it is concluded that the main reason for the local suspension of the pipeline in 2007 is the high-intensity typhoon weather in 2007. Based on the wave current coupled numerical model, the erosion and deposition of the pipeline route area after Typhoon Lily in 2007 were calculated. The calculation results show that the area near KP300 where the suspended pipeline was found in 2007 is the most serious scour, which also verifies the above inference.

Conclusions
A numerical model of storm waves in the coastal area of Zhejiang is established. The calculated wind and wave fields caused by typhoon 0712 are in good agreement with the calculated results of NOAA WAVEWATCH III, which indicates that the model can reflect the distribution of wind and waves in this region under typhoon conditions. Analysis and calculation show that under extreme weather conditions, the contribution of strong waves to the bottom shear stress cannot be ignored. Only when the influence of waves is considered, the actual situation can be reflected in the calculation of local scour.
The contribution of waves (especially under extreme conditions) to the bottom shear stress is not to be ignored. In this paper, a numerical model of storm waves in the coastal area of Zhejiang is established. The calculated wind and wave fields caused by typhoon 0712 are in good agreement with the calculated results of NOAA WAVEWATCH III, which indicates that the model can reflect the distribution of wind and waves in this region under typhoon conditions. Because the area near the submarine pipeline KP300 of the Chunxiao gas field group is located at the entrance of tidal channel, when the typhoon storm comes, it will inevitably cause a strong current. Restricted by the terrain, the flow, wind and waves are roughly in the same direction. At this time, the shear stress at the bottom of the bed will increase sharply, and a large amount of sediment will be lifted. Because the area is relatively far from the shore, the sediment cannot be replenished in time, and the sediment deposition is far from offset by the amount of scouring under extreme conditions, which leads to regional scouring at a rapid rate, resulting in exposed or even suspended pipelines, endangering their safety.
Coupled with the hydrodynamic model and wave model under the action of the typhoon, and the calculations around typhoon 0712, the results obtained from the coupling model are provided to the sediment model to obtain the regional scour distribution under the typhoon. The results show that the area around KP300 shows strong scour under typhoon action, and the calculation results are consistent with the monitoring results of the pipeline. It shows that the typhoon has an impact on the safety of submarine pipelines in this area. It is the first time in China that a numerical model is used to study the safety condition of submarine pipelines in a specific area, specific environment and specific extreme natural conditions, considering the combined action of wave and current. This new method can provide a basis for the selection of submarine pipelines in complex terrain and hydrodynamic environment, and also provides a reference for the monitoring of submarine pipelines.
Based on the monitoring data of submarine pipelines in the Chunxiao gas field group after the typhoon in 2007, it was observed that the suspension section of the pipeline was near KP300. Based on the analysis of topography, marine hydrology, meteorological conditions and submarine sediment, it was inferred that the strong scour and pipeline suspension near KP300 were caused by the typhoon. Then, the pipeline of KP287.66~KP301.906 was selected as the study area, and a numerical model was established to calculate the scouring and deposition condition of the pipelines after Typhoon Lily, a typical typhoon in 2007. Simulation results show that the most serious scouring occurs around KP300, where suspended pipelines were found in 2007 monitoring. According to the monitoring data of the submarine pipelines of the Chunxiao Gas field group for many years combined with the comprehensive environmental factors such as hydrology, sediment, meteorology and topography, it is concluded that the main reason for the partial suspension of the pipeline is the typhoon weather with high intensity and density. On the one hand, this new method can provide a sound basis for line selection of submarine pipelines under complex conditions. On the other hand, it can also provide a reference for submarine pipeline monitoring.  Data Availability Statement: Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.