Phân tích xu thế quá trình vận chuyển trầm tích và biến đổi đường bờ, đáy khu vực cửa sông đáy bằng mô hình Mike

Giá trị tối đa cho phép đối với x phải được chọn trên cơ sở này. 39 - Điều kiện Courant dưới đây có thể dùng như một hướng dẫn để chọn bước thời gian sao cho đồng thời thoả mãn được các điều kiện trên. Điển hình, giá trị của Cr là 10 đến 15, nhưng các giá trị lớn hơn (lên đến 100) đã được sử dụng: x gyVt Cr    )( với V là vận tốc. Cr thể hiện tốc độ nhiễu động sóng tại nước nông (biên độ nhỏ). Số Courant biểu thị số các điểm lưới trong một bước sóng phát sinh từ một nhiễu động di chuyển trong một bước thời gian. Sơ đồ sai phân hữu hạn dùng trong MIKE 11 (sơ đồ 6 điểm Abbott), cho phép số Courant từ 10- 20 nếu dòng chảy dưới phân giới (số Froude nhỏ hơn 1). 2.5.1.3. Mô đun AD Bên cạnh mô đun thủy động lực HD nói trên, trong MIKE có thể lựa chọn thêm mô đun AD (mô đun khuyếch tán và lan truyền chất) để tính toán và dự báo xâm nhập mặn. Phương trình khuyếch tán Song song với việc sử dụng hệ phương trình thủy động lực nói trên, khi tính toán với mô đun khuyếch tán và lan truyền chất, trong mô hình MIKE 11 sử dụng thêm phương trình khuyếch tán có dạng như sau: qCAKC x C AD xx QC t AC 2)(            (2.5) Trong đó: C là nồng độ chất ô nhiễm (chất hòa tan); D là hệ số khuyếch tán; A là diện tích mặt cắt ngang; K là hệ số tự phân hủy tuyến tính; C2 là nồng độ của nguồn gia nhập/ra khỏi của hệ thống; q là lượng gia nhập khu giữa; x,t là tọa độ theo không gian và thời gian. Phương trình trên đây phản ánh hai quá trình diễn ra đồng thời: 1. Quá trình vận chuyển bình lưu theo dòng chảy 2. Quá trình khuyếch tán do sự chênh lệch về nồng độ chất hòa tan Phương trình (2.5) được xây dựng dựa trên cơ sở các giả thiết: chất hòa tan xáo trộn toàn bộ và xem là có mật độ đồng nhất trên tất cả mọi điểm thuộc mặt cắt 40 ngang, đồng nghĩa với việc xem rằng các nguồn vào hay ra đều ngay lập tức được xáo trộn đều trên toàn mặt cắt ngang, chất hòa tan đó được xem là bảo toàn hoặc nếu có tự phân hủy thì quá trình phân hủy đó là tuyến tính và định luật khuyếch tán của Fick được áp dụng ở đây, nghĩa là sự vận chuyển do khuyếch tán tỷ lệ với gradient về nồng độ. 2.5.2. Mô hình MIKE 21 2.5.2.1. Mô hình tính sóng MIKE 21 SW MIKE 21 SW là mô đun tính phổ sóng gió được tính toán dựa trên lưới phi cấu trúc. Mô đun này tính toán sự phát triển, suy giảm và truyền sóng tạo ra bởi gió và sóng lừng ở ngoài khơi và khu vực ven bờ. Động lực học của sóng trọng lực (the dynamics o the gravity wave) được mô phỏng dựa trên phương trình mật độ tác động sóng (wave action density). Khi áp dụng tính cho vùng nhỏ thì phương trình cơ bản được sử dụng trong hệ toạ độ Cartesian, còn khi áp dụng cho vùng lớn thì sử dụng hệ toạ độ cầu (spherical polar coordinates). Phổ mật độ tác động sóng thay đổi theo không gian và thời gian là một hàm của 2 tham số pha sóng. Hai tham số pha sóng là vevtor sóng k với độ lớn k và hướng θ. Ngoài ra, tham số pha sóng cũng có thể là hướng sóng θ và tấn suất góc trong tương đối σ hoặc tần suất góc tuyệt đối ω. Trong mô hình này thì hướng sóng θ và tấn suất góc tương đối σ được chọn để tính toán. Tần số phổ được giới hạn giữa giá trị tần số cực tiểu σ min và tấn số cực cực đại σ max Năng lượng mật độ khi đó là: mrEE  )(),(),( max2max   với m = 5. MIKE 21 SW bao gồm hai công thức khác nhau: - Công thức tham số tách hướng - Công thức phổ toàn phần Công thức tham số tách hướng được dựa trên việc tham số hoá phương trình bảo toàn hoạt động sóng. Việc tham số hoá được thực hiện theo miền tần số bằng cách đưa vào mô men bậc không và bậc một của phổ hoạt động sóng giống như các giá trị không phụ thuộc (theo Holtuijsen 1989). Xấp xỉ tương tự được sử dụng trong 41 mô đun phổ sóng gió ven bờ MIKE 21 NSW. Công thức phổ toàn phần được dựa trên phương trình bảo toàn hoạt động sóng, như được mô tả bởi Komen và cộng sự (1994) và Young (1999), tại đó phổ hướng sóng của sóng hoạt động là giá trị phụ thuộc. Các phương trình cơ bản được xây dựng trong cả hệ toạ độ Đề các với những áp dụng trong phạm vi nhỏ và hệ toạ độ cầu cho những áp dụng trong phạm vi lớn hơn. MIKE 21 SW bao gồm các hiện tượng vật lý sau: - Sóng phát triển bởi tác động của gió; - Tương tác sóng-sóng là phi tuyến; - Tiêu tán sóng do sự bạc đầu; - Tiêu tán sóng do ma sát đáy; - Tiêu tán sóng do sóng vỡ; - Khúc xạ và hiệu ứng nước nông do sự thay đổi độ sâu; - Tương tác sóng- dòng chảy; - Ảnh hưởng của thay đổi độ sâu theo thời gian. Việc rời rạc hoá phương trình trong không gian địa lý và không gian phổ được thực hiện bằng cách sử dụng phương pháp thể tích hữu hạn lưới trung tâm. Sử dụng kỹ thuật lưới phi cấu trúc trong miền tính địa lý. Việc tích phân theo thời gian được thực hiện bằng cách sử dụng xấp xỉ chia đoạn trong đó phương pháp hiện đa chuỗi được áp dụng để tính truyền sóng. Phương trình cơ bản chính là phương trình cân bằng tác động sóng được xây dựng cho cả hệ toạ độ Đề các và toạ độ cầu (xem Komen và cộng sự (1994) và Young (1999)). Phương trình cho tác động sóng được viết như sau: (2.6) trong đó N(σ,θ)là mật độ hoạt động; t là thời gian; x  =(x,y)là toạ độ Đề các đối với hệ toạ độ Đề các và x  =(  ,λ) là toạ độ cầu trong toạ độ cầu với  là vĩ độ và λ là kinh độ; ),,,(  ccccv yx  là vận tốc truyền nhóm sóng trong không gian bốn chiều 42 v  , σ và θ; và S là số hạng nguồn cho phương trình cân bằng năng lượng.  là toán tử sai phân bốn chiều trong không gian v  , σ và θ. Điều kiện biên Ở biên đất trong không gian địa lý, điều kiện biên trượt toàn phần được áp dụng. Các thành phần thông lượng đi vào được gán bằng không. Ở biên mở, thông lượng đi vào cần được biết. Do đó, phổ năng lượng phải được xác định ở biên mở. 2.5.2.2. Mô hình tính thủy lực Mike 21FM HD Mô đun dòng chảy được giải bằng phương pháp lưới phần tử hữu hạn. Mô đun này dựa trên nghiệm số của hệ các phương trình Navier-Stokes trung bình Reynolds cho chất lỏng không nén được 2 hoặc 3 chiều kết hợp với giả thiết Boussinesq và giả thiết áp suất thuỷ tĩnh. Do đó, mô đun bao gồm các phương trình: phương trình liên tục, động lượng, nhiệt độ, độ muối và mật độ và chúng được khép kín bởi sơ đồ khép kín rối. Với trường hợp ba chiều thì sử dụng xấp xỉ chuyển đổi hệ toạ độ sigma. Việc rời rạc hoá không gian của các phương trình cơ bản được thực hiện bằng việc sử dụng phương pháp thể tích hữu hạn trung tâm. Miền không gian được rời rạc hoá bằng việc chia nhỏ miền liên tục thành các ô lưới/phần tử không trùng nhau. Theo phương ngang thì lưới phi cấu trúc được sử dụng còn theo phương thẳng đứng trong trường hợp 3 chiều thì sử dụng lưới có cấu trúc. Trong trường hợp hai chiều các phần tử có thể là phần tử tam giác hoặc tứ giác. Trong trường hợp ba chiều các phần tử có thể là hình lăng trụ tam giác hoặc lăng trụ tứ giác với các phần tử trên mặt có dạng tam giác hoặc tứ giác. Phương trình cơ bản Phương trình liên tục (2.7) Phương trình động lượng theo phương x và y tương ứng 43 (2.8) (2.9) Trong đó, t là thời gian; x, y và z là toạ độ Đề các; η là dao động mực nước; d là độ sâu; h= η+d là độ sâu tổng cộng; u, v và w là thành phần vận tốc theo phương x, y và z; =2Ω sinΦ là tham số Coriolis; g là gia tốc trọng trường; ρ là mật độ nước; ν t là nhớt rối thẳng đứng; pa là áp suất khí quyển; ρ 0 là mật độ chuẩn;. S là độ lớn của lưu lượng do các điểm nguồn và (u S ,v S ) là vận tốc của dòng lưu lượng đi vào miền tính. F u , F v là các số hạng ứng suất theo phương ngang. Phương trình tải cho nhiệt và muối (2.10) (2.11) trong đó D v là hệ số khuếch tán rối thẳng đứng; Hˆ là số hạng nguồn do trao đổi nhiệt với khí quyển. T S và S S là nhiệt độ và độ muối của nguồn; F T và F S là các số hạng khuếch tán theo phương ngang. Phương trình tải cho đại lượng vô hướng (2.12) 44 trong đó C là nồng độ của đại lượng vô hướng; kp là tốc độ phân huỷ của đại lượng đó; C S là nồng độ của đại lượng vô hướng tại điểm nguồn; D v là hệ số khuếch tán thẳng đứng; và F C là số hạng khuếch tán ngang. Điều kiện biên Biên đất Dọc theo biên đất thông lượng được gán bằng không đối với tất cả các giá trị. Với phương trình động lượng điều này gây ra sự trượt toàn phần dọc theo biên đất. Biên mở Điều kiện biên mở có thể được xác định dưới cả dạng lưu lượng hoặc mực nước cho các phương trình thuỷ động lực. Với phương trình tải thì giá trị xác định hoặc chênh lệch xác định có thể được đưa vào. Phương pháp giải MIKE 21 áp dụng sơ đồ sai phân hữu hạn và phương pháp giải hiệu quả là kỹ thuật ADI (Alternating Direction Implicit) để giải các phương trình bảo toàn khối lượng và động lượng trong miền không gian và thời gian. Các ma trận phương trình kết quả đối với mỗi hướng và mỗi đường lưới tính toán được giải bằng thuật giải quét đúp (Double Sweep). Các phương trình trên được giải bằng phương pháp sai phân hữu hạn theo sơ đồ QUICKEST do Lars Ekebjerg và Peter Justesen đề xướng 1997. Để giải hệ phương trình trên, người ta đã sử dụng phương pháp ADI (Alternating Direction Implicit) để sai phân hoá theo lưới không gian - thời gian. Hệ phương trình theo từng phương và tại mỗi điểm trong lưới được giải theo phương pháp Double Sweep (DS). Biểu diễn các thành phần theo các phương được thể hiện trên hình 2.8. 45 Hình 2.8. Các thành phần theo phương x và y 2.5.2.3. Mô hình tính vận chuyển trầm tích MIKE 21 ST MIKE 21 ST là mô đun tính toán tốc độ vận chuyển trầm tích (cát) không kết dính dưới tác động của cả sóng và dòng chảy. Các thành phần vận chuyển trầm tích có thể gây ra biến đổi đáy. Việc tính toán được thực hiện dưới điều kiện thuỷ động lực cơ bản tương ứng với độ sâu đã cho. Không có sự tương tác trở lại của thay đổi độ sâu đến sóng và dòng chảy. Do đó, kết quả cung cấp bởi MIKE 21 ST có thể được sử dụng để xác định khu vực có khả năng xói hoặc bồi và để chỉ ra tốc độ biến đổi đáy nhưng không xác định được việc cập nhật độ sâu ở cuối mỗi chu kỳ tính toán. Đặc trưng chính của mô đun vận chuyển trầm tích không kết dính MIKE 21 ST được mô tả như sau: - Các đặc trưng của vật chất đáy có thể không đổi hoặc biến đổi theo không gian (ví dụ tỉ lệ và cỡ hạt trung bình) - Năm lý thuyết vận chuyển trầm tích khác nhau đều có giá trị cho việc tính toán tốc độ vận chuyển trầm tích trong điều kiện chỉ có dòng chảy: + Lý thuyết vận chuyển tổng tải Engelund và Hansen + Lý thuyết vận chuyển tổng tải (được xác định như tải đáy + tải lơ lửng) Engelund và Fredsoe + Công thức vận chuyển tổng tải (tải đáy + tải lơ lửng) Zyserman và Fredsoe + Lý thuyết vận chuyển tải đáy Meyer-Peter 46 + Công thức vận chuyển tổng tải Ackers và White - Hai phương pháp có giá trị để tính toán tốc độ vận chuyển trầm tích kết hợp giữa sóng và dòng chảy + Áp dụng mô đun vận chuyển trầm tích STP của DHI + Phương pháp vận chuyển tổng tải của Bijker - Phương pháp vận chuyển cát do người sử dụng xác định (2 chiều hoặc tựa 3 chiều) trong tính toán kết hợp sóng và dòng chảy khi mô đun STP được sử dụng. Tính toán tốc độ vân chuyển được đẩy mạnh thông qua việc sử dụng bảng vận chuyển trầm tích được tạo ra trước đó. - Sử dụng STP cho phép tính toán ảnh hưởng của hiện tượng sau đến tốc độ vận chuyển trầm tích: + Hướng truyền sóng bất kỳ tác động đến dòng chảy + Sóng vỡ hoặc sóng không vỡ + Đặc tính hình học của vật chất đáy được mô tả thông qua một loại cỡ hạt hoặc đường cong phân bố cỡ hạt + Đáy phẳng hoặc đáy gợn cát - Tính ổn định chuẩn Courant-Friedrich-Lewy. Phân bố thẳng đứng của trầm tích lơ lửng trong tính toán sóng kết hợp với dòng chảy dùng để đánh giá vận chuyển trầm tích trong biển. Cách thông thường để mô tả phân bố thẳng đứng của trầm tích lơ lửng đó là áp dụng phương trình khuếch tán: (2.13) trong đó c là nồng độ trầm tích; t là thời gian; w là tốc độ chìm lắng của trầm tích lơ lửng; y là toạ độ thẳng đứng; s là thừa sô trao đổi rối. 2.5.3. Mô hình LITPACK 2.5.3.1. Khái quát về mô hình Litpack Litpack (Littoral transport and coastline kinetics) nằm trong gói phần mềm MIKE của Viện Thủy lực và Môi trường Đan Mạch (DHI). Mô hình này có thể tính toán được 47 các quá trình vận chuyển trầm tích ven bờ và diễn biến đường bờ nhằm phục vụ các bài toán chỉnh trị cửa sông và kỹ thuật đường bờ. Trong mô hình này, có các mô đun mô phỏng các quá trình ven bờ riêng biệt và có liên kết động với nhau. Do đó, các quá trình biến đổi phức tạp của đường bờ có thể miêu tả một cách chi tiết thông qua các mô đun này. Mô hình Litpack bao gồm 5 mô đun. Trong đó có hai mô đun cơ sở: mô đun Litstp, mô đun Litdri t; 3 mô đun tính các đặc tính khác nhau của quá trình vận chuyển trầm tích: mô đun Litline, mô đun Litpro và mô đun Littren. Hình 2.9. Các mô đun trong mô hình Litpack 2.5.3.2. Các mô đun trong Litpack a) Mô đun tính dịch chuyển trầm tích không kết dính Litstp Mô đun này tính toán sự dịch chuyển của các hạt trầm tích dưới tác động của sóng và dòng chảy. Đây là môđun cơ sở cho các mô đun khác trong mô hình Litpack và mô hình MIKE 21 ST. Trong mô hình, quá trình vận chuyển trầm tích được tính toán theo lý thuyết của Engelund và Fredsoe (1976) và thông qua hàm số ứng suất đáy: Mô đun tính dịch chuyển trầm tích không kết dính Litstp Mô đun tính vận chuyển trầm tích dọc bờ Litdrift Mô đun tính quá trình diễn biến đường bờ Litline Mô đun tính quá trình bồi lấp luồng tàu Littren Mô đun tính quá trình biến đổi địa hình đáy Litpro 48 gds U f )1( 2   (2.14) Khối lượng trầm tích dịch chuyển q s , được tính toán thông qua nồng độ C. z C w zzt C s               (2.15) Tính toán vận chuyển trầm tích phụ thuộc vào vận tốc, nó tính bằng tổng các khối lượng vận chuyển trong một chu kỳ sóng T: dydtCU T q T D d meanS   0 2 1 (2.16) Trong đó: s là mật độ tương đối của trầm tích; g: gia tốc trọng trường; d: kích thước hạt; U f : vận tốc; t: thời gian; S : hệ số khuếch tán rối; w: vận tốc dòng theo phương thẳng đứng b) Mô đun tính toán vận chuyển trầm tích dọc bờ Litdrift Mô đun này kết hợp mô đun Litstp ở trên với các điều kiện thủy động lực để tính toán dòng chảy dọc bờ. Mô đun này có thể tính toán được lượng vận chuyển trầm tích, nó là cơ sở cho việc xem xét nghiên cứu hình thái đường bờ. Phương trình cơ bản của mô đun này xây dựng trên phương trình cân bằng lực theo phương song song với bờ. curw xy b dy ds dy du ED dy d         (2.17) Trong đó: b : ứng suất đáy do dòng dọc bờ;  : mật độ nước; E: hệ số trao đổi thông lượng; D: độ sâu khối nước; u: vận tốc dòng chảy dọc bờ; s xy : năng lượng bức xạ; w và cur : ứng suất gió và ứng suất do dòng chảy gây ra. c) Mô đun tính toán biến đổi địa hình đáy Litpro Mô đun này dựa trên cơ sở của Litstp mô tả các mặt cắt vuông góc với bờ, bằng phương trình liên tục đối với trầm tích đáy. x q nt h S       1 1 (2.18) 49 Các mô hình hình thái không thể tính vận chuyển trầm tích trực tiếp từ các thông số sóng, dòng chảy mà phải tính qua các hàm trung gian. L qq t q slssl    (2.19) Với h: độ sâu đáy; n: độ xốp của vật liệu đáy; L: chiều dài quy mô (tỉ lệ thuận với độ sâu tại nơi nghiên cứu). d) Mô đun tính diễn biến đường bờ Litline Đây là mô đun tính diễn biến đường bờ với các ứng dụng: nghiên cứu biến động đường bờ dưới ảnh hưởng của các yếu tố tự nhiên; nghiên cứu biến động đường bờ dưới ảnh hưởng của các công trình bảo vệ bờ; nghiên cứu biện pháp khôi phục đường bờ bằng phương pháp nuôi bãi nhân tạo. Mô đun tích hợp các kết quả của các mô đun Litstp và Litdri t để tính toán đưa ra kết quả. Litline tính toán diễn biến đường bờ căn cứ vào công thức liên tục phụ thuộc vào tổng lượng bùn cát vận chuyển dọc bờ Q(x): xxh xQ x xQ xht xy act sou act c        )( )()( )( 1)( (2.20) Trong đó: )(xyc : khoảng cách từ đường cơ sở đến đường bờ; T: thời gian tính; h act (x): chiều cao hữu hiệu của mặt cắt ngang tính toán; Q(x): lượng bùn cát vận chuyển dọc bờ; X: vị trí dọc bờ tính toán; x : độ lớn bước tính toán; Q sou (x): lượng bùn cát vận chuyển qua đơn vị tính một đoạn x e) Mô đun tính biến đổi địa hình lòng dẫn Littren Mô đun này kết hợp với các mô đun tính toán di chuyển hạt trầm tích Litstp để xác định mức độ thay đổi địa hình lòng dẫn vào luồng sông. Mô đun tính toán sự thay đổi địa hình dưới tác động đồng thời của cả sóng và dòng chảy. Mô đun được xây dựng trên hệ phương trình gồm: Phương trình nồng độ trầm tích: (2.21) Trong đó: c: nồng độ trầm tích; ω: tốc độ lắng đọng; t: thời gian; εs: hệ số nhớt rối 50 Độ cao địa hình đáy được tính theo công thức: (2.22) Trong đó: z: độ cao đáy; qs: lưu lượng trầm tích và n: độ xốp vật liệu đáy 51 CHƯƠNG 3. TÍNH TOÁN CÁC QUÁ TRÌNH VẬN CHUYỂN TRẦM TÍCH VÀ BỒI TỤ, XÓI LỞ 3.1. Đặt vấn đề Trong nghiên cứu này đã sử dụng bộ mô hình MIKE nghiên cứu tổng thể các quá trình vận chuyển bùn cát, xói lở và bồi tụ bờ biển khu vực cửa Đáy thông qua 3 nghiên cứu chính sau đây: - Nghiên cứu và phân tích xu thế vận chuyển trầm tích trong thời gian 1 năm bằng mô hình MIKE 21 ST. - Nghiên cứu tính toán biến đổi bùn cát dài hạn có xét đến dâng cao của mực nước biển và mô hình quá trình phát triển cửa Đáy trong thời gian 20 năm bằng mô hình MIKE 21 ST. - Nghiên cứu, tính toán mức độ thay đổi đường bờ trong thời gian 20 năm bằng mô hình LITPACK. Để giải quyết được các hướng nghiên cứu trên, báo cáo này đã sử dụng các một số mô hình cơ bản trong bộ mô hình MIKE theo sơ đồ nghiên cứu sau: Hình 3.1. Các bộ thông số và các mô hình toán cơ bản sử dụng trong nghiên cứu 52 3.2. Xây dựng bộ số liệu cơ sở cho mô hình 3.2.1. Địa hình, miền tính, lưới tính Địa hình miền tính cho mô hình MIKE 21 được lấy từ số liệu đo đạc của Bộ Tư lệnh Hải quân từ các bản đồ địa hình đáy biển với tỉ lệ khác nhau từ tỉ lệ 1:10,000 đến 1:1,000,000. Tọa độ miền tính từ 19046’N đến 20021’N và 105056’E đến 106045’E. Trong đó, các bản đồ tỉ lệ lớn được dùng cho khu vực ven bờ và các đảo; bản đồ tỉ lệ nhỏ dùng cho vùng ngoài khơi. Để thực hiện cho nhiệm vụ của nghiên cứu, trong nghiên cứu này đã lựa chọn lưới phần tử hữu hạn tăng dần từ ngoài biển vào trong sát bờ. Diện tích nhỏ nhất của 1 phần tử là 1250m2 ở khu vực các cửa sông như: cửa Đáy, cửa Ninh Cơ và cửa Ba Lạt. Diện tích lớn nhất là 25km 2 ở khu vực biên ngoài khơi. Miền tính có 2879 nút điểm, với độ phân giải thô nhất ở vùng ngoài khơi là 5000m, mịn nhất ở vùng bờ khu vực cửa sông là 50m. Hình 3.2. Địa hình khu vực nghiên cứu 53 Hình 3.3. Minh họa lưới tính sử dụng trong mô phỏng 3.2.2. Điều kiện biên Trong các mô hình thủy lực được sử dụng của nghiên cứu này, gồm 3 mô hình: MIKE 21 SW, MIKE 21 FM HD, MIKE 21 ST. Các mô hình có thể tính toán riêng biêt, nhưng trong nghiên cứu này sử dụng kết hợp liên hợp 3 mô hình trên. Trước tiên, tính toán mô hình MIKE 21 SW để tính toán phổ sóng, sau đó kết quả của mô hình SW làm số liệu đầu vào cho 2 mô hình còn lại. Tại biên ngoài khơi: Trong mô hình MIKE 21 SW các thông số đầu vào gồm: độ cao, chu kỳ và hướng sóng cho các biên ngoài khơi. Với mô hình MIKE 21 HD sử dụng mực nước từ phân tích điều hòa làm điều kiên biên ngoài khơi. Điều kiện biên cho mô hình MIKE 21 ST là cân bằng nồng độ. Tại biên cửa sông: Thông số đầu vào là lưu lượng dòng chảy từ sông ra, số liệu được tính toán từ mô hình MIKE 11. 3.2.3. Các thông số khác 54 Trường khí tượng được lấy theo số liệu thống kê trương khí tượng tại trạm đo Văn Lý trong “Sổ tay tra cứu các đặc trưng khí tượng thủy văn vùng thềm lục địa Việt Nam”. Các đặc trưng cấp hạt và nồng độ trầm tích ban đầu được lấy từ số liệu đo đặc cấp độ hạt tại khu vực Cửa Đáy trong 2 đọt khảo sát: tháng 7 và tháng 11 năm 2009. Số liệu được bổ sung thêm theo các nghiên cứu về địa chất, địa mạo khu vực này trước đây [3,7]. Thời gian mô phỏng: - Bài toán phân tích xu thế bồi tụ và xói lở: 1 năm bắt đầu từ 01/01/2010 tới ngày 31/12/2010. - Bài toán biến đổi bùn cát dài hạn cửa Đáy: 20 năm, bắt đầu từ ngày 01/01/2010 tới 31/12/2030. - Bài toán biến đổi đường bờ: 20 năm từ 1/1/1990 tới 1/1/2010. 3.3. Hiệu chỉnh và kiểm nghiệm mô hình 3.3.1. Mô hình MIKE 11 Mô hình MIKE 11 cho hệ thống sông Hồng – Thái Bình trong dự báo tác nghiệp được sử dụng để tính toán lưu lượng dòng chảy và nồng độ bùn cát tại các trạm thủy văn cửa sông Ba Lạt, Phú Lễ và Như Tân. Do không có số liệu bùn cát thực đo tại các sông Hồng và Ninh Cơ số liệu nồng độ bùn cát đo đạc tại mặt cắt trạm thủy văn Như Tân từ ngày 1/11/2009 đến ngày 30/11/2009 được sử dụng để hiệu chỉnh. Kết quả tính toán cho thấy, có sự tương đồng lớn về pha và độ lớn giữa nồng độ bùn cát tính toán và giá trị thực đo. 55 0 0.04 0.08 0.12 0.16 0.2 1/11/2009 0:00 7/11/2009 0:00 13-11-2009 00:00:00 19-11-2009 00:00:00 25-11-2009 00:00:00 Thời gian (giờ) Nồ ng độ (k g/m 3) Tính toán Thực đo Hình 3.4. So sánh nồng độ trầm tích tại mặt cắt trạm Như Tân trên sông Đáy 3.3.2. Mô hình tính sóng MIKE 21 SW Hiệu chỉnh mô hình bằng số liệu sóng trung bình đặc trưng mùa tại trạm hải văn Văn Lý (20007’N;106018’E) trong sổ tay tra cứa các đặc trưng khí tượng thủy văn thềm lục địa Việt Nam [9]. (a) Tháng 2 – 4 (b) Tháng 5 – 7 (c) Tháng 8 – 10 (d) Tháng 11 – 1 Hình 3.5. Độ cao và hướng sóng đặc trưng cho các tháng trong năm 56 Kết quả mô phỏng trường sóng đặc trưng cho từng giai đoạn bằng mô hình MIKE 21 SW tương đối phù hợp với số liệu thống kê. Kết quả này có thể sử dụng để tính toán các quá trình động lực cũng như vận chuyển bùn cát nhằm phục vụ các bài toán khác. 3.3.3. Mô hình thủy lực MIKE 21 FM Hiệu chỉnh mô hình dựa theo số liệu đo đạc mực nước tại cửa sông Ninh Cơ (vị trí HC, xem trong hình 3.6 trong đợt khảo sát thuộc đề tài: “Nghiên cứu phương pháp dự báo biên triều cho các cửa sông” của Trung tâm Nghiên cứu Biển và tương tác Biển – Khí quyển thực hiện tháng 8 năm 2010. Chuỗi số liệu 7 ngày, từ 12h ngày 19-08-2010 tới 12h ngày 26-08-2010. Hình 3.6. Vị trí các điểm hiệu chỉnh và kiểm nghiệm mô hình 57 Hình 3.7. So sánh mực nước thực đo và tính toán tại điểm HC (trạm đo Ninh Cơ;106012’7.14”E, 2001’26.49”N) Kiểm nghiện mô hình bằng cách so sánh số liệu mực nước tính toán chỉ có yếu tố triều bằng mô hình với số liệu mực nước tính toán từ bộ hằng số điều hòa tại 3 điểm gần bờ. Thời gian kiểm nghiệm 1 tháng, bắt đầu từ 7h ngày 01-09-2010 tới 7h ngày 01-10-2010. Vị trí 3 điểm kiểm nghiệm KN1, KN2, KN3 xem trong hình 3.6. Kết quả so sánh mực nước kiểm nghiệm được trình bày qua các hình ảnh dưới đây: m -1.5 -1 -0.5 0 0.5 1 1.5 9/1/2010 0:00 9/11/2010 0:00 9/21/2010 0:00 10/1/2010 0:00 Thời gian Mô hình HSĐH Hình 3.8. So sánh mực nước tính toán bằng mô hình và tính toán từ bộ hằng số điều hòa tại điểm KN1 (106035’E, 20013’N) 58 m -1.5 -1 -0.5 0 0.5 1 1.5 9/1/2010 0:00 9/11/2010 0:00 9/21/2010 0:00 10/1/2010 0:00 Thời gian Mô hình HSĐH Hình 3.9. So sánh mực nước tính toán bằng mô hình và tính toán từ bộ hằng số điều hòa tại điểm KN2 (106006’E, 19055’N) m -1.5 -1 -0.5 0 0.5 1 1.5 9/1/2010 0:00 9/11/2010 0:00 9/21/2010 0:00 10/1/2010 0:00 Thời gian Mô hình HSĐH Hình 3.10. So sánh mực nước tính toán bằng mô hình và tính toán từ bộ hằng số điều hòa tại điểm KN3 (105056’E, 19051’N) Nhìn vào các hình trên, cho thấy kết quả tính toán có sự tương đồng cao về pha và biên độ mực nước giữa kết quả mô phỏng của mô hình với số liệu thực đo/tính toán bằng hằng số điều hòa trong cả quá trình hiệu chỉnh mô hình (hệ số tương quan là 0,99548) và kiểm nghiệm mô hình (hệ số tương quan từ 0,88089 tới 0,9675). Do vậy, có thể sử dụng kết quả của mô hình thủy lực MIKE 21 FM để làm đầu cào phục vụ mô phỏng và tính toán vận chuyển trầm tích tại khu vực cửa sông Đáy bằng mô hình MIKE 21 ST. 3.4. Các kết quả trong nghiên cứu 3.4.1. Phân tích xu thế vận chuyển trầm tích 59 Quá trình vận chuyển trầm tích vùng cửa sông phụ thuộc vào các quá trình thủy động lực của khu vực, trong đó có 3 yếu tố quan trọng là: dòng chảy mang trầm tích từ sông đưa ra, chế độ sóng vùng ven bờ và dòng chảy khu vực ven bờ. Dòng chảy từ sông đưa ra có chu kỳ 1 năm với 2 mùa: mùa mưa và mùa khô. Chế độ sóng ven bờ chủ yếu phụ thuộc vào điều kiện gió, mà các đặc trưng gió có chu kỳ nhỏ là 1 năm với các mùa khác nhau. Khu vực này có chế độ thủy triều là nhật triều đều, với chu kỳ khoảng 30 ngày tương ứng hai con nước: con nước triều lên và con nước triều xuống với khoảng thời gian bằng nhau là 15 ngày. Dòng chảy khu vực ven bờ phụ thuộc chủ yếu vào chế độ thủy triều, do đó chu kỳ để tính toán các đăch trưng liên quan tới dòng chảy ở khu vực này khoảng 30 ngày. Do đó, trong quá trình nghiên cứu phân tích xu thế vận chuyển trầm tích cho khu vực đã tính toán với khoảng thời gian là 1 năm, bắt đầu từ ngày 01/01/2010 đến ngày 31/12/2010. 3.4.1.1. Mô phỏng thủy lực Lưu lượng -1000 -500 0 500 1000 1500 2000 1/1/ 10 2/1/ 10 3/1/ 10 4/1/ 10 5/1/ 10 6/1/ 10 7/1/ 10 8/1/ 10 9/1/ 10 10/1 /10 11/1 /10 12/1 /10 1/1/ 11 Thời gian m3/s Hình 3.11. Biến trình lưu lượng qua mặt cắt tại trạm Như Tân trên sông Đáy năm 2010 60 Lưu lượng -1000 -500 0 500 1000 1500 2000 1/1/ 10 2/1/ 10 3/1/ 10 4/1/ 10 5/1/ 10 6/1/ 10 7/1/ 10 8/1/ 10 9/1/ 10 10/1 /10 11/1 /10 12/1 /10 1/1/ 11 Thời gian m3/s Hình 3.12. Biến trình lưu lượng qua mặt cắt tại trạm Phú Lễ trên sông Ninh Cơ năm 2010 Lưu lượng -1000 -500 0 500 1000 1500 2000 1/1/ 10 2/1/ 10 3/1/ 10 4/1/ 10 5/1/ 10 6/1/ 10 7/1/ 10 8/1/ 10 9/1/ 10 10/1 /10 11/1 /10 12/1 /10 1/1/ 11 Thời gian m3/s Hình 3.13. Biến trình lưu lượng qua mặt cắt tại trạm Ba Lạt trên sông Hồng năm 2010 Kết quả tính toán lưu lượng dòng chảy từ 3 sông đổ ra khu vực nghiên cứu ta thấy: lưu lượng dòng chảy lớn nhất xảy ra vào tháng 8 và 9 tưng ứng với mùa mưa ở miền bắc. Lưu lượng dòng chảy lớn nhất tại mặt cắt trên sông Hồng là 2000 m3/s xảy ra vào ngày 1/9/2010. Tổng lượng nước chuyển qua 3 mặt cắt vào mùa mưa (từ tháng 6 tới tháng 10) chiếm khoảng ¾ tổng lượng nước năm. 61 Hình 3.14. Mực nước tính toán lúc 0h ngày 15 tháng 8 năm 2010 Hình 3.15. Mực nước tính toán lúc 0h ngày 15 tháng 1 năm 2010 3.4.1.2. Mô phỏng phân bố trầm tích 62 Nồng độ trầm tích 0 0.2 0.4 0.6 0.8 1/1/10 2/1/10 3/1/10 4/1/10 5/1/10 6/1/10 7/1/10 8/1/10 9/1/10 10/1/10 11/1/10 12/1/10 1/1/11 Thời gian mg/l Hình 3.16. Biến trình nồng độ trầm tích tại mặt cắt trạm Như Tân năm 2010 Nồng độ trầm tích 0 0.2 0.4 0.6 0.8 1/1/10 2/1/10 3/1/10 4/1/10 5/1/10 6/1/10 7/1/10 8/1/10 9/1/10 10/1/10 11/1/10 12/1/10 1/1/11 Thời gian mg/l Hình 3.17. Biến trình nồng độ trầm tích tại mặt cắt trạm Phú Lễ năm 2010 Nồng độ trầm tích 0 0.2 0.4 0.6 .8 1 1/1/10 2/1/10 3/1/10 4/1/10 5/1/10 6/1/10 7/1/10 8/1/10 9/1/10 10/1/10 11/1/10 12/1/10 1/1/11 Thời gian mg/l Hình 3.18. Biến trình nồng độ trầm tích tại mặt cắt trạm Ba Lạt năm 2010 63 Hình 3.19. Biến đổi đáy lúc 0h ngày 1-6-2010 (thời điểm đầu mùa mưa) Hình 3.20. Địa hình đáy biển tại thời điểm 0h ngày 1-6-2010 (thời điểm đầu mùa mưa) 64 Hình 3.21. Biến đổi đáy lúc 0h ngày 1-10-2011 (thời điểm cuối mùa mưa) Hình 3.22. Địa hình đáy biển tại thời điểm 0h ngày 1-10-2011(thời điểm cuối mùa mưa) 65 Hình 3.23. Biến đổi đáy lúc 0h ngày 31-12-2010 (sau 1 năm tính toán) Hình 3.24. Địa hình đáy biển tại thời điểm 0h ngày 31-12-2010 (sau 1 năm tính toán) 3.4.1.3. Nhận xét 66 Đây là khu vực có động lực sóng tương đối yếu ở gần bờ, dòng chảy dọc bờ và sông là yếu tố quyết định đến quá trình bồi lắng. Ở khu vực này, trầm tích chủ yếu là cát hạt nhỏ và mịn với đường kính hạt trung bình là 0,17mm. Theo tính toán, lưu lượng dòng chảy từ sông Đáy đưa ra khoảng 195.000 m3/ năm. Ngoài ra, lưu lượng dòng chảy tổng hợp cả năm qua mặt cắt tại bờ khu vực Nghĩa Hưng – Nam Định có hướng đi từ bắc xuống nam với lưu lượng khoảng 1.807.000 m3/năm. Lượng bùn cát của khu vực cửa Sông Đáy được tiếp nhận từ 2 nguồn lớn: nguồn thứ nhất là nguồn từ Sông Đáy đưa ra khoảng 34.000 tấn bùn cát và nguồn thứ hai là dòng dọc bờ mang bùn cát từ phía bắc trở xuống khoảng 220.000 tấn bùn cát. Xét về tính chất mùa ta thấy rõ: lượng bùn cát từ sông mang ra được vận chuyển chủ yếu vào mùa mưa, chiếm khoảng 80% tổng lượng bùn cát cả năm (xem các hình từ 3.16 tới 3.18). Ngoài ra, tốc độ dòng chảy sông lớn mang bùn cát ra xa bờ, nơi đó chịu động lực biển phân bố bồi tụ trầm tích theo diện rộng. Do đó, độ bồi tụ đáy vào mùa mưa nhanh hơn và làm cao thêm bề mặt đáy. Sau đó, giai đoạn từ cuối mùa mưa đến cuối năm 2010 ta thấy: lượng bùn cát vẫn được sông đưa ra nhưng với nồng độ và lưu lượng thấp, tốc độ dòng chảy yếu không thể đưa bùn cát ra xa bờ được gây tích tụ tại lòng dẫn sông, cho nên độ cao đáy ở đây được nâng lên. Còn khu vực ngoài, dưới tác động của thủy động lực biển đã làm thay đổi đáy. Quá trình thủy động lực mang dòng bùn cát đưa tới cửa sông, rồi lắng đọng tạo nên các khu vực bồi tụ trước cửa và hai phía của cửa Đáy. Ta thấy rõ khu vực cửa Sông Đáy có xu thế bồi ở cả 2 bên của cửa, độ cao bồi tụ khoảng 1m. Khu vực bồi mạnh nhất là khu vực trước cửa Đáy bồi cao đến gần 2m. Khu vực phía đông cửa Đáy thuộc huyện Nghĩa Hưng – tỉnh Nam Định có xu thế bồi và hình thành dọc bờ. Khu vực phía tây thuộc huyện Kim Sơn – Ninh Bình bồi rất mạnh nhưng có xu thế bồi lắng rộng ra ngoài xuống phía nam. Tóm lại: Khu vực cửa Sông Đáy được tiếp nhận 2 nguồn bổ sung trầm tích rất lớn từ sông Đáy đưa ra và từ phìa bắc đưa xuống. Các con sông đổ vào khu vực nghiên cứu đều có lưu lượng dòng chảy lớn mang theo nguồn trầm tích dồi dào. 67 Chế độ thủy động lực của khu vực làm lắng đọng trầm tích ở 2 phía của cửa Đáy. Dòng bùn cát này bị ngăn lại bởi dòng chảy của sông Đáy và lắng đọng lại ở ven bờ, lấp đầy khoảng trống giữa cửa Đáy và cồn ngầm ở phía ngoài. Vì vậy, bãi bồi có diện tích rất lớn và khoảng cách giữa cồn ngầm và bãi bồi phía trong hầu như không đáng kể. 3.4.2. Tính toán xu thế biến động bùn cát dài hạn có xét đến dâng cao mực nước biển và mô hình hóa quá trình phát triển cửa Đáy 3.4.2.1. Kịch bản nước biển dâng cho khu vực cửa Đáy Các kịch bản phát thải khí nhà kính được lựa chọn để tính toán, xây dựng kịch bản nước biển dâng cho Việt Nam là kịch bản phát thải thấp (kịch bản B1), kịch bản phát thải trung bình của nhóm các kịch bản phát thải trung bình (kịch bản B2) và kịch bản phát thải cao nhất của nhóm các kịch bản phát thải cao (kịch bản A1FI). Các kịch nước biển dâng được xây dựng cho bảy khu vực bờ biển của Việt Nam sử dụng trong nghiên cứu này được lấy từ Kịch bản biến đổi khí hậu, nước biển dâng do Bộ tài nguyên và Môi trường công bố năm 2012, theo đó, kịch bản nước biển dâng cho khu vực Cửa Đáy được đưa ra trong bảng 3.1 và hình 3.25. Bảng 3.1. Mực nước biển dâng theo các kịch bản phát thải cho khu vực Cửa Đáy (cm) Năm Kịch bản cao A1F1 Kịch bản trung bình B2 Kịch bản thấp B1 Cận trên Cận dưới Cận trên Cận dưới Cận trên Cận dưới 2020 9 8 8 7 9 8 2030 14 12 13 11 13 11 2040 19 16 18 15 17 15 2050 27 22 24 20 23 19 2060 36 30 32 25 30 24 2070 47 38 39 31 37 29 2080 59 47 48 37 44 34 2090 72 56 56 43 51 38 2100 86 66 65 49 58 42 68 0 20 40 60 80 100 120 2000 2010 2020 2030 2040 2050 2060 2070 2080 2090 2100 Năm M ức tă ng m ực n ướ c (c m ) Cận trên - Kịch bản A1FI Cận dưới - Kịch bản A1FI Cận trên - Kịch bản B2 Cận dưới - Kịch bản B2 Cận trên - Kịch bản B1 Cận dưới - Kịch bản B1 Hình 3.25. Kịch bản mực nước biển dâng tại khu vực cửa Đáy Theo kịch bản phát thải thấp (B1) thì vào giữa thế kỷ 21, tại khu vực Cửa Đáy, mực nước biển dâng trong khoảng từ 19-23cm, đến cuối thế kỷ 21, mực nước biển dâng trong khoảng từ 42-58cm; Theo kịch bản phát thải trung bình (B2) thì vào giữa thế kỷ 21, tại khu vực Cửa Đáy, mực nước biển dâng trong khoảng từ 20- 24cm, đến cuối thế kỷ 21, mực nước biển dâng trong khoảng từ 49-65cm; Theo kịch bản phát thải cao (A1FI) thì vào giữa thế kỷ 21, tại khu vực Cửa Đáy, mực nước biển dâng trong khoảng từ 22-27cm, đến cuối thế kỷ 21, mực nước biển dâng trong khoảng từ 66-86cm. 3.4.2.2. Cập nhật mực nước biển dâng trong mô hình Theo như đã trình bày ở phần 3.4.2.1, mực nước biển có xu hướng tăng lên trong thời gian từ năm 2010 tới 2030. Do đó, trong nghiên cứu này đã sử dụng số liệu dự báo mức độ tăng của mực nước biển cho khu vực cửa Đáy để thêm vào diễn biến mực nước (bằng tính toán từ hằng số điều hòa) tại các biên, làm đầu vào cho mô hình thủy lực MIKE 21 HD. 69 Theo như đã nêu ở trên, giá trị tăng mực nước biển năm 2020 khoảng từ 7 – 9 cm, năm 2030 tăng từ 11 – 14 cm tương ứng với từng kịch bản phát thải (so với mốc tính là mực nước trong thời kỳ 1980 – 1999). Theo đó, mực nước tăng trung bình trong giai đoạn từ năm 2010 tới 2030 khoảng 0,5cm/năm. Cho nên, nghiên cứu này đã sử dụng giái trị tăng 0,5cm/năm là giá trị thêm vào biến trình mực nước trong từng năm tại các biên để tính toán. 3.4.2.3. Lưu lượng dòng chảy sông Trong nghiên cứu này, sử dụng giả thiết diễn biến dòng chảy trên 3 sông Hồng, Ninh Cơ và Đáy trong giai đoạn từ năm 1990 tới năm 2010 tương đồng với giai đoạn từ năm 2010 tới 2030. Do đó, nghiên cứu đã sử dụng số liệu tính toán lưu lượng tại 3 cửa sông Hồng, Ninh Cơ và Đáy qua các mặt cắt tương ứng là Ba Lạt, Phú Lễ và Như Tân trong giai đoạn từ năm 1990 tới năm 2010 làm số liệu đầu vào. Số liệu được tính toán bằng mô hình MIKE 11 nhằm đơn giản hóa quá trình mô phỏng thủy động lực. -1000 0 1000 2000 3 0 4000 5000 60 7000 1/1/89 1/1/94 1/1/99 1/1/04 1/1/09 thời gian m3/s Lưu lượng (a) Trạm Ba Lạt trên sông Hồng 70 -1000 0 1000 2000 3000 4000 5000 1/1/89 1/1/94 1/1/99 1/1/04 1/1/09 thời gian m3/s Luu luong (b) Trạm Phú Lễ trên sông Ninh Cơ -1000 0 1000 2 300 4000 1/1/89 1/1/94 1/1/99 1/1/04 1/1/09 thời gian m3/s Lưu lượng (c) Trạm Như Tân trên sông Đáy Hình 3.26. Biến trình lưu lượng dòng chảy qua mặt cắt các trạm trong 20 năm 3.4.2.4. Kết quả a) Mô phỏng thủy lực 71 Mực nước -2 -1 0 1 2 1/1/10 1/1/15 1/1/20 1/1/25 1/1/30 thời gian (m) Hình 3.27. Biến trình mực nước tính toán tại điểm gần cửa Đáy (10605’E; 19 050’N) b) Kết quả phân bố bùn cát dài hạn Hình 3.28. Địa hình khu vực cửa Đáy ban đầu 72 (a) (b) Hình 3.29. Địa hình (a) và mức độ biến đổi địa hình (b) đáy sau 1 năm tính toán (a) (b) Hình 3.30. Địa hình (a) và mức độ biến đổi địa hình (b) đáy sau 2 năm tính toán (a) (b) Hình 3.31. Địa hình (a) và mức độ biến đổi địa hình (b) đáy sau 3 năm tính toán 73 (a) (b) Hình 3.32. Địa hình (a) và mức độ biến đổi địa hình (b) đáy sau 4 năm tính toán (a) (b) Hình 3.33. Địa hình (a) và mức độ biến đổi địa hình (b) đáy sau 5 năm tính toán (a) (b) Hình 3.34. Địa hình (a) và mức độ biến đổi địa hình (b) đáy sau 10 năm tính toán 74 (a) (b) Hình 3.35. Địa hình (a) và mức độ biến đổi địa hình (b) đáy sau 15 năm tính toán (a) (b) Hình 3.36. Địa hình (a) và mức độ biến đổi địa hình (b) đáy sau 20 năm tính toán c) Nhận xét Nhìn vào các kết quả tính toán ta có thể thấy rõ: - Hình 3.27: thể hiện kết quả tính toán dao động mực nước của 1 điểm trong khu vực nghiên cứu có xu hướng tăng, biểu hiện qua đường xu thế màu đỏ đi lên. Mức độ tăng của dao động mực nước tương ứng với mức độ tăng do kịch bản nước biển dâng đề xuất ở mục 3.4.2.1. - So sánh các hình từ 3.29 tới 3.36 với hình 3.28 cho ta thấy địa hình đáy biển của từng giai đoạn tính toán và mức độ thay đổi địa hình đáy biển của khu vực nghiên cứu. Có sự sai khác nhất định của địa hình đáy biển sau các khoảng thời gian tính toán so với địa hình lúc ban đầu. Khu vực có xu thế bồi mạnh. 75 - So sánh các hình từ 3.29 tới 3.36 với nhau cho ta thấy: địa hình khu vực nghiên cứu có xu thế bồi tụ mạnh trong 15 năm đầu tiên tính toán, quá trình bồi tụ mạnh làm bồi lấp luồng vào cửa Đáy. Hậu quả của nó làm thay đổi dòng chảy. Sau khoảng thời gian 15 năm đã bồi lấp mất luồng vào cửa Đáy, không còn dòng chảy từ cửa sông đưa ra. Trong khoảng thời gian sau 15 năm chỉ biến đổi địa hình tại các khu vực đáy ở vùng ngoài. - Nhìn vào các hình từ 3.29 tới 3.36 có thể nhận thấy các khu vực có biến động mạnh nhất là: khu vực luồng vào cửa Đáy và khu vực ngoài của cửa Ninh Cơ. Mức độ biến động lên tới trên 10m. - Ngoài ra, khu vực luôn luôn được tiếp nhận hàng năm lượng bùn cát rất lớn từ sông Đáy đưa ra và từ phía bắc đưa xuống. Lượng bùn cát từ phía bắc được cung cấp do sông Hồng, sông Ninh Cơ. Một phần lượng bùn cát này được đưa xuống khu vực cửa Đáy gây bồi lấp. Khu vực cửa Đáy luôn luôn được tiếp nhận một lượng rất lớn bùn cát từ sông Đáy đưa ra và từ phía bắc đưa xuống. Do đó, khu vực này có xu thế bồi lấp cửa Đáy, trong khoảng thời gian ngắn nữa cửa Đáy sẽ bị lấp. Nhằm phòng tránh các khó khăn do bồi lấp cửa Đáy gây ra cần phải có kế hoạch nạo vét cửa, mở rộng luồng sông. Tỉnh Ninh Bình đã và đang triển khai dự án “Nạo vét sông Đáy đoạn từ cầu Gián Khẩu đến cửa Đáy và các tuyến nhánh qua sông Đáy để thoát lũ cho sông Hoàng Long, tỉnh Ninh Bình” để tránh sự cố do bồi lấp luồng vào cửa Đáy. Quá trình tính toán đã cho thấy xu thế bồi tụ khu vực cửa Đáy, kết quả này tương đối chính xác với các nghiên cứu gần đây về khu vực này. Tuy nhiên, cần có những căn cứ khoa học hơn nữa ví dụ như kiểm nghiệm mô hình tính lan truyền bùn cát (MIKE 21 ST) để có thể đưa ra thời gian bồi lấp mất cửa Đáy. Đây cũng là lý do tại sao mà cửa Đáy thường xuyên phải nạo vét để có thể lưu thông hàng hải bình thường và hàng năm có nhiều diện tích đất được bồi tụ. Ngoài ra, qua nghiên cứu này cho thấy mức độ cấp thiết của các số liệu đo đạc về sóng, dòng chảy, lưu lượng để cho kết quả tốt hơn, phù hợp với thực tế hơn. 76 3.4.3. Tính toán biến đổi đường bờ có xét đến dâng cao mực nước biển do biến đổi khí hậu 3.4.3.1. Điều kiện tính toán Trong nghiên cứu biến động đường bờ có xét đến mực nước biển dâng do biến đổi khí hậu, để kiểm nghiệm lại độ đúng đắn của tính toán bằng cách so sánh đường bờ tính toán với đường bờ thu được từ ảnh vệ tinh. Do đã có bộ ảnh vệ tinh đường bờ khu vực ở các mốc thời gian là các năm 1990, 2001 và năm 2010 cho nên ở nghiên cứu này đã tính toán với quãng thời gian là 20 năm bắt đầu từ năm 1990. Tuy nhiên, khó khăn gặp phải là không có trường địa hình, không có bộ thông số lưu lượng cũng như dòng chảy từ các sông đưa ra từ năm 1990, cho nên để đơn giản tác giả đã sử dụng các điều kiện tính toán sau: - Trong mô hình Litpack các pro ile mặt cắt chỉ có sự thay đổi theo hướng tịnh tiến ra xa ngoài bờ được lấy tại mốc năm 2010 thay cho mốc năm 1990. - Quá trình biến đổi lưu lượng, tốc dộ dòng chảy từ năm 2010 tới năm 2030 chính là quá trình biến đổi từ năm 1990 tới năm 2010. - Tốc độ tăng của mực nước biển dâng do BĐKH trong giai đoạn từ năm 2010 tới năm 2030 bằng với giai đoạn từ năm 1990 tới năm 2010. 3.4.3.2. Bộ thông số đầu vào - Dữ liệu đầu vào + Dữ liệu về địa hình: được xác định thông qua số liệu địa hình trình bày ở phần trên. + Dữ liệu về trầm tích: Với kích thước hạt trung bình tại khu vực nghiên cứu được lấy là 0.17mm, độ chọn lọc 1.44 + Dữ liệu về mực nước: lấy mực nước trung bình là 0 + Dữ liệu về sóng: được lấy từ số liệu sóng tính toán bằng mô hình MIKE 21 SW cho khu vực nghiên cứu. + Dữ liệu về dòng chảy: được lấy từ bộ số liệu tính toán bằng mô hình MIKE 21 FM. 77 + Dữ liệu về gió: Bộ thông số về trường gió được lấy theo các đặc trưng gió thống kê tại trạm Văn Lý. * Xác định vị trí đường bờ Số liệu đường bờ được lấy từ ảnh vệ tinh năm 1989 làm số liệu phục vụ tính toán. Vị trí đường bờ được xác định là khoảng cách từ đường bờ tới đường cơ sở. Có thể định nghĩa đường cơ sở là một đường thẳng tương đối song song với đường bờ mà quá trình bồi xói không vượt quá ranh giới đó. Hình 3.37. Đường bờ và đường cơ sở Trong nghiên cứu này, khu vực cửa Đáy được mô phỏng dựa trên đường cơ sở có hệ số góc =1, tương đối song song với đoạn bờ cần tính toán, chi tiết biểu diễn qua hình ảnh dưới đây, đường bờ gồm 550 điểm, mỗi điểm cách nhau 30m. Ứng với mỗi điểm sẽ có một mặt cắt địa hình với 2 mặt cắt địa hình đặc trưng cho 2 phía của cửa sông. Mỗi mặt cắt địa hình chứa 600 nút điểm, mỗi nút điểm cách nhau 30m. 78 (a) (b) Hình 3.38. Đường cơ sở, khu vực nghiên cứu (a) và biểu diễn đường bờ năm 1989 (b) trên đường cơ sở 79 (a) Phân bố mặt cắt địa hình đặc trưng trong tính toán biến đổi đường bờ Mặt cắt địa hình MC1 -25 -20 -15 -10 -5 0 5 1 51 101 151 201 251 301 351 401 451 501 551 vị trí m (b) Mặt cắt địa hình MC1 Mặt cắt địa hình MC2 -25 -20 -15 -10 -5 0 5 1 51 101 151 201 251 301 351 401 451 501 551 vị trí m (c) Mặt cắt địa hình MC2 Hình 3.39. Phân bố mặt cắt địa hình và địa hình sử dụng trong nghiên cứu + Thời gian tính toán: 40 năm từ 01/01/1990 tới 01/01/2030. 80 3.4.3.3. Kết quả tính toán 0 5000 10000 15000 1 51 101 151 201 251 301 351 401 451 vị trí m ĐB tính toán 2001 ĐB tính từ ảnh vệ tinh 2001 Hình 3.40. So sánh đường bờ tính toán bằng mô hình và đường bờ trích ra từ số liệu vệ tinh năm 2001 0.00 5000.00 10000.00 15000.00 1 51 101 151 201 251 301 351 401 451 vị trí m ĐB tính toán 2010 ĐB tính từ ảnh vệ tinh 2010 Hình 3.41. So sánh đường bờ tính toán bằng mô hình và đường bờ trích ra từ số liệu vệ tinh năm 2010 81 Biến đổi ĐB giai đoạn 1990 - 2001 (ve tinh) 0 200 400 600 800 1000 1 51 101 151 201 251 301 351 401 451 vị trí m Hình 3.42. Mức độ biến đổi đường bờ khu vực Cửa Đáy giai đoạn từ năm 1990 tới năm 2001 (số liệu tính toán từ ảnh vệ tinh) Biến đổi ĐB giai đoạn 1990 - 2010 (ve tinh) 0 200 400 600 80 10 1200 1 51 101 151 201 251 301 351 401 451 vị trí m Hình 3.43. Mức độ biến đổi đường bờ khu vực Cửa Đáy giai đoạn từ năm 1990 tới năm 2010 (số liệu tính toán từ ảnh vệ tinh) Nhận xét: Nhìn vào số liệu tính toán đường bờ cho thấy: Kết quả tính toán tương đối phù hợp với số liệu được lấy theo đường bờ từ ảnh vệ tinh kể cả năm 2001 và năm 2010. Đường bờ đoạn phía huyện Kim Sơn tỉnh Ninh Bình có hiện tượng bồi tụ nhưng không nhiều, còn đoạn bờ phía bờ đông thuộc huyện Nghĩa Hưng tỉnh Nam Định thì đang bồi tụ rất mạnh. Khu vực giáp ranh với cửa Đáy lấn ra biển hàng năm khoảng gần 300m. Tuy nhiên, do khu vực phía đông của cửa Đáy có hiện tượng bồi 82 tụ rất mạnh, khi tham chiếu lên đường cơ sở đã không mô phỏng được rõ ràng nên kết quả tính toán còn nhỏ hơn so với số liệu đường bờ lấy từ ảnh khoảng 500m. Theo kết quả được trình bày trong các hình từ 3.44 tới 3.46 cho thấy: khu vực hai phía bờ của Cửa Đáy có xu thế bồi tụ lấn ra biển. Khu vực bờ phía đông cửa Đáy (huyện Nghĩa Hưng tỉnh Nam Định) bồi tụ mạnh hơn khu vực bờ phía tây cửa Đáy (huyện Kim Sơn tỉnh Ninh Bình). Khu vực có mức độ biến động mạnh nhất là hai bờ của cửa Đáy, tốc độ lên tới 80m/năm (trong mô hình là các vị trí từ 230 tới 400). Khu vực có tốc độ bồi tụ nhỏ nhất trong nghiên cứu này là khu vực giáp ranh với tỉnh Thanh Hóa (vị trí từ 0 cho tới 70), tốc độ chỉ khoảng 1m/năm. Nhìn chung kết quả trong nghiên cứu trong giai đoạn từ năm 1990 tới 2010 cho thấy hai phía bờ của cửa Đáy đang có xu thế bồi và đặc biệt là phía đông cửa Đáy đang bồi tụ rất mạnh. Kết quả này phù hợp với thực tế. 83 KẾT LUẬN Kết quả nghiên cứu cho thấy, khu vực cửa Đáy có xu thế bồi tụ là chủ yếu, tạo thành các cồn cát trước cửa Đáy và các doi cát dọc bờ thuộc địa phận huyện Nghĩa Hưng tỉnh Nam Định. Nguyên nhân do khu vực này được tiếp nhận 2 nguồn bổ sung trầm tích từ: i) sông Đáy đưa ra và lắng động dưới tác động của động lực biển; ii) lượng bùn cát được vận chuyển dọc theo bờ biển từ phía bắc (cửa Ba Lạt và Ninh Cơ). Ngoài ra, qua nghiên cứu này cho thấy, vào 4 tháng mùa mưa từ tháng 6 tới tháng 10 lượng trầm tích tích được các con sông đưa ra chiếm khoảng 80% lượng trầm tích cả năm. Khu vực cửa Đáy luôn luôn được tiếp nhận một lượng rất lớn bùn cát từ sông Đáy đưa ra và từ phía bắc đưa xuống. Do đó, khu vực này có xu thế bồi lấp cửa Đáy, trong khoảng thời gian ngắn nữa (theo tính toán có thể 15 năm) cửa Đáy sẽ bị lấp nếu không có kế hoạch nạo vét luồng vào cửa. Để tránh các khó khăn do bồi lấp cửa Đáy gây ra cần phải có kế hoạch nạo vét cửa, mở rộng luồng sông. Tỉnh Ninh Bình đã và đang triển khai dự án “Nạo vét sông Đáy đoạn từ cầu Gián Khẩu đến cửa Đáy và các tuyến nhánh qua sông Đáy để thoát lũ cho sông Hoàng Long, tỉnh Ninh Bình” để tránh sự cố do bồi lấp luồng vào cửa Đáy. Vấn đề vận chuyển trầm tích và biến động hình thái – địa mạo sông, cửa sông và ven biển vẫn là vấn đề khó khăn cả trong lý thuyết và thực tiễn. Nghiên cứu này đã ứng dụng bộ mô hình MIKE nhằm phân tích xu thế bồi tụ và xói lở cũng như các mối tương tác Biển – Lục địa khu vực cửa Đáy. Tuy nhiên, do điều kiện hạn chế về mặt thời gian và số liệu, nghiên cứu này chỉ xem xét và đánh giá được xu thế biến động địa hình cũng như đường bờ ở các điều kiện khí hậu đặc trưng. Do vậy, để đánh giá hiện tượng một cách đầy đủ và khách quan, trong các nghiên cứu tiếp theo, cần xem xét đánh giá bổ sung những hạn chế trên. Bộ mô hình MIKE đã được sử dụng trong nghiên cứu, đánh giá và phân tích xu thế bồi xói và vận chuyển trầm tích tại khu vực cửa sông Đáy và cho kết quả khả quan. Qua đó khẳng định, có thể sử dụng bộ mô hình MIKE trong nghiên cứu tương tác các quá trình thủy thạch động lực. 84 TÀI LIỆU THAM KHẢO Tiếng Việt 1. Bộ Tài nguyên và Môi trường (2011), Kịch bản biến đổi khí hậu và nước biển dâng cho Việt Nam. 2. Đặng Ngọc Thành, Phạm Văn Ninh, Mai Thanh Tân…(2005), Chuyên khảo Biển Đông, NXB Đại học Quốc Gia Hà Nội. 3. Đỗ Minh Đức (2004), Nghiên cứu sự hình thành và biến đổi quá trình bồi tụ- xói lở ở đới ven biển Thái Bình-Nam Định, Luận án Tiến sĩ chuyên ngành Địa chất công trình. 4. Lương Phương Hậu, Trịnh Việt An, Lương Phương Hợp (2002), Diễn biến cửa sông vùng đồng bằng Bắc Bộ, Nhà xuất bản Xây dựng. 5. Nguyễn Biểu, Vũ Trường Sơn, Dương Văn Hải và nnk (2001), Địa chất khoáng sản biển nông ven bờ (0-30 m nước) Việt Nam tỷ lệ 1/500.000. Lưu trữ ĐC, Hà Nội. 6. Nguyễn Thị Hải Lý, Lương Phương Hậu (2004), “Ứng dụng Litpack trong nghiên cứu diễn biến cửa sông và bờ biển”, Tạp chí GTVT số 10. 7. Nguyễn Thế Tưởng, Sổ tay tra cứu các đặc trưng khí tượng thủy văn vùng thềm lục địa Việt Nam, NXB Nông Nghiệp, 2000. 8. Nguyễn Thọ Sáo, Nguyễn Minh Huấn, Ngô Chí Tuấn, Đặng Đình Khá (2010), “Biến động trầm tích và diễn biến hình thái khu vực cửa sông ven bờ Cửa Tùng, Quảng Trị”, Tạp chí khoa học ĐHQGHN, khoa học tự nhiên và công nghệ, tập 26, số 3S, 427. 9. Nguyen Xuan Hien, Duong Ngoc Tien, Le Quoc Huy, Nguyen Tho Sao (2012), “Long-term sediment distribution calculation taking into account sea level rise and theb development o Day estuary”, Tạp chí khoa học ĐHQGHN, khoa học tự nhiên và công nghệ, tập 28, số 3S, tr.57 - 62. 10. Nguyễn Xuân Hiển, Dương Ngọc Tiến (2012) Nguyễn Thọ Sáo. “Tính toán và phân tích xu thế bồi tụ, xói lở khu vực Cửa Đáy”, Hội thảo khoa học Quốc Gia về KT, TV, MT & BĐKH, lần thứ XV,Viện KH KTTV&MT, Tập 2, tr.241-246. 11. Phạm Quang Sơn (2004), “Diễn biến lòng dẫn hạ lưu sông Hồng trong 15 năm vận hành khai thác nhà máy thủy điện Hoà Bình” TC Các Khoa học về Trái đất, 26/4: 520-531. Hà Nội. 12. Trần Việt Liễn, Hoàng Đức Cường, Trương Anh Sơn, Trần Trung Thành (2006), “Xây dựng các kịch bản (scenarios) về biến đổi khí hậu của thế kỷ 85 XXI cho các vùng thuộc lãnh thổ Việt Nam”, Tạp chí KTTV số 541 (tháng 1 năm 2006), Hà Nội. Tiếng Anh 13. Aslak Grinsted, J. C. Moore & S. Jevrejeva (2009), Reconstructing sea level from paleo and projected temperatures 2000 to 2100 ad. 14. A.M. prospathopoulos, A. Sotiropoulos, E. Chatziopoulos, C.H. Anagnostou (2004), “Cross-shore profile and coastline changes of a sandy beach in Pieria, Greece, based on measurements and numerical simulation”, Mediterranean Marine Science, vol 5/1. 15. Cazenave, A., and R. S. Nerem (2004), “Present-day sea level change: Observations and causes”, Rev. Geophys., 42, RG3001. 16. DHI (2003), An Intergrated Modelling System for Littoral Processes and Coastline Kinetics, Short Introduction and Tutorial, DHI Software. 17. Hidalgo, H. G., M. D. Dettinger, and D. R. Cayan (2008), “Downscaling with Constructed Analogues: Daily Precipitation and Temperature Fields over the United States”. California Energy Commission, PIER Project Report, 2008. 18. Houghton, J.T., Meira Filho, L.G., Callander, B.A., Harris, N., Kattenberg, A., Maskell, K. (1996), “Climate Change 1995: The Science o Climate Change”. Contribution of Working Group I to the Second Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge Univ. Press, Cambridge, UK, pp. 363-397. 19. Jianjun Yin (2009), Sea level rise due to global warming poses threat to New York City. 20. John A.Church, Neil J. White (2006), “A 20th century acceleration in global sea level rise”, Geophysical research letter. 21. MIKE 21 FLOW Model (2005), Hydrodynamic Module – Scientific Documentation, DHI Software. 22. MIKE 21 ST Non-Cohesive Sediment Transport Module (2005), User Guide, DHI Softwave.