Luận án Nghiên cứu và ứng dụng chương trình Flexpart trong đánh giá phát tán phóng xạ tầm xa - Phạm Kim Long

à GFS và ECMWF (Hình 2.6), đây là hai mô hình khí tượng của Trung tâm Dự báo Môi trường Quốc gia Mỹ (NCEP) và Trung tâm Dự báo Thời tiết hạn vừa châu Âu (ECMWF). Trong nghiên cứu này, mô hình khí tượng CFSv2 của NCEP [56] sẽ được sử dụng chính làm dữ liệu khí tượng đầu vào cho chương trình FLEXPART trong bài toán mô phỏng phát tán phóng xạ từ Fukushima. Mô hình khí tượng GFS Analysis và Forecasts sẽ được sử dụng trong bài toán đánh 65 giá khả năng phát tán phóng xạ và dự báo phát tán phóng xạ trong trường hợp xảy ra sự cố từ các NMĐHN gần biên giới Việt Nam. Hình 2.6. Mô hình khí tượng đầu vào cho chương trình FLEXPART. Đối với bài toán phát tán phóng xạ tầm xa, hiện nay đang sử dụng dữ liệu khí tượng CFSv2 NCEP [56]. CFSv2 là phiên bản thứ hai của Hệ thống dự báo khí hậu của Trung tâm dự báo môi trường Hoa Kỳ NCEP. Hệ thống được vận hành từ tháng 3 năm 2011, cũng chính là khoảng thời gian bắt đầu xảy ra tai nạn hạt nhân Fukushima. Hệ thống cung cấp dữ liệu khí tượng bắt đầu từ 0h00 ngày 1/1/2011 cho đến thời điểm hiện tại với gần 50 thông số khí tượng. Dữ liệu được cập nhập liên tục theo ngày trên hệ thống, khởi tạo bốn lần trên ngày vào các khung giờ 0000, 0600, 1200 và 1800 UTC. Đây là mô hình tương tự như phiên bản đầu tiên CFSR (Climate Forecast System Reanalysis) [57]. Các dữ liệu phân tích và dự báo khí quyển, đại dương và đất liền trong 6 giờ, ở các độ phân giải ngang là 0.2, 0.5, 1.0, và 2.5 độ. Dữ liệu được cung cấp dưới định đạng GRIB2. Yêu cầu các thông số khí tượng đầu vào đối với chương trình FLEXPART như trong Bảng 2.6, bao gồm 12 thông số khí tượng, dữ liệu lấy theo độ phân giải 0.5x0.5 độ. Đối với bài toán phát tán phóng xạ tầm xa, độ phân giải 0.5x0.5 66 có thể đáp ứng được yêu cầu của bài toán. Các nghiên cứu của Stohl [5] và Drexler [58] cũng đánh giá cao các kết quả thu được từ hệ thống dữ liệu GFS so với dữ liệu khí tượng ECMWF của cơ quan khí tượng châu Âu. Bảng 2.6. Các thông số khí tượng đầu vào mô hình CFSv2 NCEP. Biến Thông số Kiểu Đơn vị DPT Dewpoint temperature Specified height level above ground (2m) K HGT Geopotential height Isobaric surface gpm HGT Geopotential height Ground or water surface gpm LAND Land cover (0 = sea, 1 = land) Ground or water surface - PRES Pressure Ground or water surface Pa PRES Pressure Mean sea level Pa PRMSL Pressure reduced to MSL Mean sea level Pa RH Relative humidity Isobaric surface % RH Relative humidity Specified height level above ground (2m) % RH Relative humidity Sigma level % TMP Temperature Specified height level above ground (2m) K TMP Temperature Sigma level K TMP Temperature Ground or water surface K TMP Temperature Isobaric surface K HPBL Planetary Boundary Layer Height Ground or water surface m VVEL Vertical velocity (pressure) Sigma level Pa/s VVEL Vertical velocity (pressure) Isobaric surface Pa/s WEASD Water equivalent of accumulated snow depth Ground or water surface kg.m-2 UGRD u-component of wind Sigma level m/s UGRD u-component of wind Specified height level above ground (10m) m/s UGRD u-component of wind Isobaric surface m/s VGRD v-component of wind Sigma level m/s VGRD v-component of wind Specified height level above ground (10m) m/s VGRD v-component of wind Isobaric surface m/s 67 2.5.2. Thu thập dữ liệu khí tượng thời điểm xảy ra tai nạn Trong nghiên cứu này, chúng tôi sử dụng dữ liệu khí tượng mô hình CFSv2 cung cấp bởi Trung tâm Dự báo Môi trường Quốc gia Mỹ NCEP [56], theo lưới chia 0.5ox0.5o từ 0oE đến 359.5oE và từ 90oN đến 90oS. Các bước thu thập dữ liệu như sau: Bước 1: Truy cập kho lưu trữ dữ liệu khí tượng "NCEP Climate Forecast System Version 2 (CFSv2) 6-hourly Products" (Hình 2.7) tại địa chỉ "https://rda.ucar.edu/datasets/ds094.0". Bước 2: Đăng ký tài khoản và đăng nhập để có quyền truy cập kho dữ liệu khí tượng. Bước 3: Truy cập vào mục "Data Access" để lấy dữ liệu. Bước 4: Trong mục "Customizable Data Requests" lựa chọn mục "Get a Subset". Hình 2.7. Truy cập kho dữ liệu khí tượng CFSv2 NCEP. Bước 5: Lựa chọn khoảng thời gian xảy ra tai nạn Fukushima (Hình 2.8) và các thông số khí tượng đầu vào phù hợp với chương trình FLEXPART (như trong Bảng 2.6). 68 Hình 2.8. Lựa chọn thời gian và thông số khí tượng đầu vào cho FLEXPART. Hình 2.9. Lựa chọn chi tiết các thông số khí tượng CFSv2 cho FLEXPART. Bước 6: Lựa chọn chi tiết các thông số khí tượng (Hình 2.9) với định dạng dữ liệu khí tượng GRIB2. 69 Bước 7: Tiến hành gửi yêu cầu cung cấp dữ liệu khí tượng và tải dữ liệu khí tượng đầu vào cho chương trình FLEXPART. Ngoài dữ liệu mô hình CFSv2 thu thập như trên, ta có thể sử dụng dữ liệu khí tượng mô hình GFS Reanalysis và Forecasts phục vụ cho bài toán dự báo phát tán phóng xạ. Truy cập vào kho dữ liệu khí tượng "Global Forecast System (GFS)" để thu thập tại địa chỉ "https://www.ncdc.noaa.gov/data-access/model- data/model-datasets/global-forcast-system-gfs". 2.5.3. Phân tích dữ liệu khí tượng Trong nghiên cứu mô hình đánh giá phát tán phóng xạ tầm xa, cần thiết phải phân tích dữ liệu khí tượng để có thể diễn giải được yếu tố khí tượng tác động vào quá trình lan truyền của phóng xạ trong bầu khí quyển (Mục 1.2). Kết hợp với kết quả thu được từ mô phỏng và quan trắc, qua đó giải thích chặt chẽ được các nguyên nhân ảnh hưởng trực tiếp đến qua trình lan truyền phóng xạ trong khí quyển. Các tệp dữ liệu khí tượng từ các mô hình CFSv2, GFS analysis và forecasts được cung cấp dưới định dạng GRIB2, có thể đọc trực tiếp bởi FLEXPART thông qua thư viện GRIB_API. Hình 2.10. Áp suất khí quyển (Pa) thời điểm xảy ra tai nạn NMĐHN Fukushima (11/3/2011). 70 Để phân tích dữ liệu khí tượng trong nghiên cứu này sử dụng phần mềm Panoply của NASA [59] để phân tích các dữ liệu khí tượng như áp suất, nhiệt độ, tốc độ gió, hướng gió (Hình 2.10). Qua đó phân tích diễn biến khí tượng để thấy được các ảnh hưởng chính của các vòng tuần hoàn khí quyển tác động lên quá trình lan truyền của phóng xạ trong bầu khí quyển. Qua đó kiểm chứng các kết quả thu được từ mô hình và thực tế. Ngoài ra, các phần mềm phân tích dữ liệu khí tượng khác cũng được sử dụng như NCL, WGRIB [60, 61]. 2.6. Phân tích và xử lý số liệu mô phỏng đầu ra 2.6.1. Giới thiệu Như nội dung nghiên cứu tổng quan về chương tình mô phỏng phát tán phóng xạ FLEXPART, định dạng đầu ra FLEXPART với tùy chọn nhị phân hoặc netCDF với các phiên bản 9.2 trở lên. Đối với định dạng nhị phân, các công cụ phân tích và xử lý số liệu viết dưới các ngôn ngữ như sau: - Ngôn ngữ Python: + Pflexible phát triển bởi J. Burkhart, là công cụ vẽ bản đồ phát tán phóng xạ kết hợp với quỹ đạo phát tán. Phiên bản tiếp theo hiện đang được phát triển là reflexible với nhiều tính năng mới hơn. + Quicklook phát triển bởi R. Hofman, là công cụ cũng dùng để vẽ bản đồ phát tán phóng xạ, rơi lắng khô và rơi lắng ướt. + Quickdose cũng được phát triển bởi R. Hofman, đây là chương trình cung cấp nhiều tùy chọn xử lý dữ liệu. Sử dụng trực tiếp FLEXPART làm nhân chạy mô phỏng, sau đó xử lý luôn kết quả đầu ra với các tùy chọn vẽ bản đồ như: nồng độ, rơi lắng khô, rơi lắng ướt, suất liều bức xạ qua chiếu xạ trong, chiều xạ ngoài (bề mặt và đám mây). - Ngôn ngữ Fortran: + Code chuyển dữ liệu đầu ra FLEXPART v6.2 sang định dạng CTBTO SRM phát triển bởi Petra Seibert. + Chương trình con để đọc FLEXPART 8.x/9.x phát triển bởi A. Stohl. 71 + Chương trình để đọc tệp nhị phân về hạt rơi lắng trong FLEXPART- WRF bởi Petra Seibert. - Các công cụ khác để đọc FLEXPART như IDL (Pasquale Sellitto), Matlab (S. Eckhardt, H. Sodemann, S. Raja). Trong đề tài nghiên cứu này, phần lớn các kết quả được phân tích và xử lý bằng chương trình Quicklook cho các định dạng dữ liệu đầu ra là nhị phân. 2.6.2. Chương trình Quicklook Chương trình Quicklook được xây dựng bởi Radek Holfman viết trên ngôn ngữ Python 2.7 [62]. Hiện tại Quicklook hỗ trợ FLEXPART các phiên bản từ 9.0.x trở về trước, các phiên bản từ 9.2 trở lại đây không được hỗ trợ do thay đổi về cấu trúc trong tệp đầu ra headers và grids. Phiên bản 9.3 của CTBTO vẫn hỗ trợ do phát triển lên từ phiên bản 9.02 vẫn giữ nguyên các cấu trúc tệp đầu ra trên. Quicklook là công cụ miễn phí sử dụng giao diện dòng lệnh cho việc vẽ nhanh bản đồ các dữ liệu đầu ra của FLEXPART (Hình 2.11). Hình 2.11. Bản đồ phát tán phóng xạ toàn cầu vẽ bằng chương trình Quicklook. 72 Bảng 2.7. Các tùy chọn xử lý số liệu đầu ra FLEXPART trong Quicklook. Tùy chọn Nội dung --version Xem phiên bản sử dụng -h, --help Xem trợ giúp tùy chọn -e EXP_PATH, --explore=EXP_PATH Thông tin về đầu vào tại EXP_PATH -i INPUT, --input=INPUT Đường dẫn INPUT đến thư mục đầu ra của FLEXPART -d DOMAIN, --domain=DOMAIN Tọa độ của lãnh thổ cần vẽ -m, --max_lat_lon Sử dụng tối đa tọa độ có sẵn (ghi đè lên tùy chọn -d) -l LEVELS, --levels=LEVELS Thiết lập khoảng độ cao, ví dụ -l 0 0 thiết lập cho lớp mặt đất. -t TYPE, --type = TYPE Thiết lập kiểu đầu ra Flexpart [mother | nested] -o OUTPUT, --output = OUTPUT Đường dẫn OUTPUT thư mục lưu kết quả xử lý -r RECEPTORS, --receptors=RECEPTORS Đánh dấu các điểm nhận trên bản đồ -f DATAFACTOR, --datafactor = DATAFACTOR Nhân dữ liệu với hệ số f -u UNITLABEL, --unitlabel = UNITLABEL Hiện đơn vị người dùng xác định -n FILENAME, --filename = FILENAME Đổi tên tệp hình động GIF -q, --reverse Vẽ bản đồ theo thời gian ngược lại -x TITLE, --title=TITLE Tiêu đề của hình ảnh -z PROJECTION, --projection=PROJECTION Tùy chọn định dạng bản đồ sử dụng -s SPECIES, --species=SPECIES Chỉ định nhân phóng xạ cần vẽ --conc (mặc định) Hiển thị kết quả nồng độ --drydepo Hiển thị kết quả rơi lắng khô --wetdepo Hiển thị kết quả rơi lắng ướt Quicklook được cung cấp miễn phí cho cộng đồng sử dụng FLEXPART dưới giấy phép công cộng GPL 2.0. Người dùng có thể dễ dàng sửa đổi lệnh trong Quicklook theo nhu cầu. Tải về Quicklook thông qua dòng lệnh: 73 git clone https://radekhofman@bitbucket.org/radekhofman/quicklook.git hoặc tải trực tiếp từ nguồn: https://bitbucket.org/radekhofman/quicklook/downloads/ Quicklook sử dụng mô đun flex_81.py phát triển bởi J. Brioude (NOAA) để đọc dữ liệu đầu ra FLEXPART. Các tùy chọn trong Quicklook thể hiện trong Bảng 2.7. Quicklook tương đối tiện lợi, không yêu cầu phải biên soạn và cài đặt, chỉ cần đáp ứng các thư viện cần thiết cho Python để chạy. Yêu cầu Python phiên bản 2.7 trở lên, bao gồm các thư viện cần thiết sau: Numpy, Matplotlib phiên bản 1.3 trở lên, Basemap, để tạo ảnh động (GIF) cần thêm ImageMagick. Đây là công cụ chính được sử dụng để vẽ bản đồ nồng độ, rơi lắng phóng xạ từ các kết quả thu được qua mô phỏng FLEXPART. Do cũng được cung cấp dưới dạng nguồn mở nên rất dễ dàng để chỉnh sửa các dòng lệnh trong Quicklook theo một số đòi hỏi đầu ra theo yêu cầu của bài toán. 2.6.3. Chương trình Quickdose Quickdose là chương trình tích hợp nhiều tính năng gồm cả vẽ bản đồ phát tán cũng như tính toán thêm các dữ liệu chiếu xạ trong và chiếu xạ ngoài thông qua dữ liệu rơi lắng và nồng độ phóng xạ thu được [63], tác giả Radek Hofman đã xây dựng chương trình trên ngôn ngữ Python, sử dụng FLEXPART viết trên Fortran làm nhân tính toán mô phỏng, các chương trình con khác làm tác vụ vẽ đồ thị và tính suất liều chiếu xạ. Tương tự như Quicklook, chương trình Quickdose tương đối tiện lợi, không yêu cầu phải biên soạn và cài đặt, chỉ cần đáp ứng các thư viện cần thiết cho Python để chạy. Bao gồm các thư viện: Numpy, Matplotlib, Basemap, Scipy, để tạo ảnh động (GIF) cần thêm ImageMagick. Quickdose được cung cấp miễn phí cho cộng đồng dưới giấy phép công cộng GPL 2.0. Người dùng có thể dễ dàng sửa đổi lệnh trong Quickdose theo nhu cầu. Tải về Quickdose thông qua dòng lệnh: git clone https://bitbucket.org/radekhofman/quickdose.git 74 hoặc tải trực tiếp từ nguồn: https://bitbucket.org/radekhofman/QuickDose/src Quickdose ngoài cung cấp bản đồ nồng độ phóng xạ phát tán trong khí quyển, còn cung cấp thêm bản đồ suất liều gamma do chiếu xạ ngoài từ đám mây hoặc từ phóng xạ rơi lắng trên mặt đất, hay chiếu xạ trong do hít phải các chất phóng xạ. Liều được tính toán như sau [64, 65]: Suất liều chiếu gamma do chiếu xạ ngoài từ đám mây hoặc bề mặt (EXT): EXT = C(t) * EC (2.19) Suất liều chiếu gamma do chiếu xạ trong do hít thở (INH): INH = C(t) * Br * DC (2.20) Trong đó: C(t) là nồng độ phóng xạ trong không khí (Bq/m3) hoặc rơi lắng dưới mặt đất (Bq/m2). EC là hệ số chuyển đổi cho đám mây hoặc rơi lắng. Br là tốc độ hít thở (m3/s). DC: là hệ số chuyển đổi liều của các nhân phóng xạ cho hít thở (Sv/Bq). Bảng 2.8. Tốc độ hít thở theo nhóm tuổi. Nhóm tuổi Tốc độ hít thở mặc định (m3/d) 0-1 tuổi 2.86 1-2 tuổi 5.16 2-7 tuổi 8.72 7-12 tuổi 15.3 12-17 tuổi 20.1 >17 tuổi 22.2 ICRP Publication 71 (1995b) Bảng 2.9. Hệ số chuyển đổi liều chiếu trong do hít thở cho 131I và 137Cs. Hệ số liều do hít thở (Sv/Bq) 0-1 tuổi 1-2 tuổi 2-7 tuổi 7-12 tuổi 12-17 tuổi >17 tuổi 131I 7.2E-8 7.2E-8 3.7E-8 1.9E-8 1.1E-8 7.4E-9 137Cs 8.8E-9 5.4E-9 3.6E-9 3.7E-9 4.4E-9 4.6E-9 ICRP (1996), Eckerman and Legett (1996) 75 Bảng 2.10. Hệ số chuyển đổi liều chiều ngoài cho 131I và 137Cs [63, 64]. Hệ số liều chiếu ngoài Đám mây (Sv s-1 Bq-1 m3) Mặt đất (Sv s-1 Bq-1 m2) 131I 1.69E-14 3.6E-16 137Cs 2.55E-14 5.51E-16 ICRP (1996), Eckerman and Legett (1996) 2.6.4. Chương trình Pflexible Pflexible phát triển bởi John F. Burkhart [66], đây là một công cụ Python cho phân tích dữ liệu từ mô hình vận chuyển FLEXPART. Pflexible ban đầu được phát triển để làm việc với FLEXPART v8.x, nó có một số tính năng mới để tạo dữ liệu đầu ra. Mô đun này được dựa trên các chương trình Matlab do Sabine Eckhardt và Harald Sodemann phát triển tại NILU. Chương trình con và các phát triển tiếp theo trên Python được thực hiện bởi John Burkhart. Pflexible được cung cấp miễn phí cho cộng đồng, tải về Pflexible thông qua dòng lệnh: git clone https://bitbucket.org/jfburkhart/pflexible pflexible hoặc tải trực tiếp từ nguồn: https://bitbucket.org/jfburkhart/pflexible/src Pflexible yêu cầu các thư viện sau: Numpy, Matplotlib, Basemap, f2py, netCDF, PIL. Tương tự như Quicklook và Quickdose, Pflexible cũng sử dụng các dữ liệu nồng độ phóng xạ thu được từ kết quả để vẽ bản đồ phát tán phóng xạ theo thời gian (Hình 2.12). Kết hợp một số tùy chọn như hiển thị quỹ đạo di chuyển của chùm phóng xạ. 76 Hình 2.12. Mô phỏng phát tán phóng xạ từ NMĐHN Fukushima sử dụng chương trình pflexible. 2.6.5.Chương trình Reflexible Reflexible được phát triển bởi John F. Burkhart [67], là bản thay thế cho Pflexible. Reflexible cũng được cung cấp miễn phí cho cộng đồng, tải về Reflexible thông qua dòng lệnh: git clone https://github.com/spectraphilic/reflexible.git hoặc tải trực tiếp từ nguồn: https://github.com/spectraphilic/reflexible Reflexible cũng tập trung xử lý dữ liệu đầu ra của FLEXPART để vẽ bản đồ phát tán phóng xạ (Hình 2.13). Hiện đang tiếp tục được phát triển và sử dụng chưa ổn định. Phiên đang vẫn đang dừng lại ở việc xử lý với định dạng nhị phân, tiến tới xử lý với các định dạng netCDF ở FLEXPART phiên bản mới. 77 Hình 2.13. Mô phỏng phát tán phóng xạ toàn cầu sử dụng chương trình reflexible. 2.7. Các phương pháp đánh giá thống kê mô hình phát tán Để đánh giá thống kê sự phù hợp giữa kết quả mô phỏng và quan trắc, sau khi thu thập các dữ liệu đo đạc tại các trạm quan trắc trong khu vực (Phụ lục 5). Nghiên cứu này sử dụng các phương pháp đánh giá các tính toán của một mô hình vận chuyển và khuếch tán trong khí quyển dựa trên các nghiên cứu trước đây [43, 68-72]. Tiếp nối các nghiên cứu trước, phương pháp xếp hạng dựa trên các thông số thống kê của Draxler [73] được sử dụng nhiều trong đánh giá chất lượng của một mô hình vận chuyển và phát tán trong khí quyển. Phương pháp xếp hạng này được định nghĩa bởi tổng của hệ số tương quan (R - correlation coefficient), phân số chênh lệnh (FB - fractional bias), hệ số phẩm chất trong không gian (FMS - figure-of-merit in space), và thông số Kolmogorov-Smirnov (KSP - Kolmogorov-Smirnov parameter). Cụ thể như sau: 78 𝑅𝑎𝑛𝑘⁡1 = 𝑅2 + (1 − | 𝐹𝐵 2 |) + ( 𝐹𝑀𝑆 100 ) + (1 − 𝐾𝑆𝑃 100 ) (2.21) Trong đó các thành phần chứa các thông số thông kê có giá trị từ 0 đến 1, nghĩa là tổng điểm xếp hạng một mô hình trong khoảng từ 0 đến 4 tương ứng với từ kém đến tốt. Hệ số tương quan (R), còn được gọi là hệ số tương quan Pearson (PCC), được sử dụng để biểu diễn sự phân tán giữa các cặp giá trị quan trắc (M) và dự đoán (P). R có giá trị trong khoảng -1 đến +1 và có giá trị lý tưởng là +1, R được tính như sau: 𝑅 = ∑(𝑀𝑖 − �̅�)(𝑃𝑖 − �̅�) √∑(𝑀𝑖 − �̅�)2(𝑃𝑖 − �̅�)2 (2.22) Phân số chênh lệnh (FB) giữa giá trị quan trắc và dự đoán nằm trong khoảng -2 đến +2, được xác định bởi biểu thức sau: 𝐹𝐵 = 2 (�̅� − �̅�) (�̅� + �̅�) (2.23) Kết quả tốt khi mà các cặp giá trị quan trắc và dự đoán nằm gần đường thẳng y = x. Nếu có N cặp (M,P) sử dụng thêm hệ số vượt (FOEX – factor of exceedance) để xác định phần trăm số cặp vượt trên hoặc dưới đường thẳng, cho bởi biểu thức sau: 𝐹𝑂𝐸𝑋 = [ 𝑁(𝑃𝑖>𝑀𝑖) 𝑁 − 0.5] × 100,−50 ≤ 𝐹𝑂𝐸𝑋 ≤ +50 (2.24) Giá trị FOEX bằng -50% nghĩa là tất cả các cặp giá trị nằm dưới đường thẳng y = x, và ngược lại giá trị bằng +50% nghĩa là tất cả các cặp giá trị nằm phía trên đường thẳng, còn với giá trị 0% có nghĩa tỉ lệ vượt trên và dưới bằng nhau không phải là N cặp nằm trên đường thẳng y = x. Sai số toàn phương trung bình (NMSD - Normalized Mean Square Deviation) được định nghĩa như sau: 𝑁𝑀𝑆𝐷 = 1 �̅� ∙ �̅� ∙ 1 𝑁 ∑(𝑀𝑖 − 𝑃𝑖) 2 (2.25) 79 NMSD cung cấp thông tin về độ lệnh không phải là độ vượt trên hay dưới. Đây là thông số rất nhạy với sự khác biệt giữa các giá trị quan trắc và dự đoán, kết quả mô hình hoàn hảo sẽ có giá trị NMSE bằng không. Tương tự với NMSE là thông số căn bậc hai của NMSE (NRMSD – normalized root mean square deviation, được tính bởi biểu thức sau: 𝑁𝑅𝑀𝑆𝐷 = 𝑅𝑀𝑆𝐷 𝑚𝑎𝑥(𝑂𝑖) − 𝑚𝑖𝑛⁡(𝑂𝑖) (2.26) 𝑅𝑀𝑆𝐷 = √ ∑(𝑀𝑖 − 𝑂𝑖)2 𝑛 (2.27) Hệ số phẩm chất trong không gian (FMS) được tính như sau: 𝐹𝑀𝑆 = 100 𝑁𝑃 ∩ 𝑁𝑀 𝑁𝑃 ∪ 𝑁𝑀 (2.28) Sự khác biệt giữa phân bố không ghép cặp của các giá trị quan trắc và dự đoán được thể hiện bằng thông số Kolmogorov-Smirnov (KSP), được định nghĩa là sự khác biệt lớn nhất giữa hai phân bố tích lũy khi Mk = Pk, trong đó: 𝐾𝑆𝑃 = 𝑀𝑎𝑥|𝐷(𝑀𝑘) − 𝐷(𝑃𝑘)| (2.29) Trong đó D là sự phân bố lũy tích của các nồng độ quan trắc và dự đoán trong khoảng k sao cho D là xác suất mà nồng độ sẽ không vượt quá Mk hoặc Pk. Đó là một thước đo cho mô hình tái tạo lại sự phân bố nồng độ đo bất kể khi nào và ở đâu. Sự khác biệt tối đa giữa hai phân bố không thể lớn hơn 100%. Hệ số FAx được định nghĩa là phần trăm của các kết quả mô phỏng nằm trong dải với hệ số 2, 5 và 10 các giá trị quan trắc, được cho bởi các biểu thức tương ứng FA2, FA5 và FA10 dưới đây: 0.5 ≤ 𝑃 𝑀 ≤ 2.0 (2.30) 0.2 ≤ 𝑃 𝑀 ≤ 5.0 (2.31) 0.1 ≤ 𝑃 𝑀 ≤ 10 (2.32) 80 Ngoài tính xếp hạng mô hình theo Rank 1 với các thông số thống kê R, FB, FMS và KSP, ta có thể đánh giá xếp hạng mô hình kết hợp các thông số thống kê như sau [33]: 𝑅𝑎𝑛𝑘⁡2 = 𝑅2 + (1 − | 𝐹𝐵 2 |) + ( 𝐹𝐴2 100 ) + (1 − 𝐾𝑆𝑃 100 ) (2.33) Trong đó Rank 2 trong khoảng từ 0 đến 4 tương ứng với từ kém đến tốt. 𝑅𝑎𝑛𝑘⁡3 = 𝑅2 + (1 − | 𝐹𝐵 2 |) + ( 𝐹𝑀𝑆 100 ) + (1 − | 𝐹𝑂𝐸𝑋 50 |) + (1 − 𝐾𝑆𝑃 100 ) (2.34) Trong đó Rank 3 trong khoảng từ 0 đến 5 tương ứng với từ kém đến tốt. 𝑅𝑎𝑛𝑘⁡4 = 𝑅2 + (1 − | 𝐹𝐵 2 |) + ( 𝐹𝑀𝑆 100 ) + ( 𝐹𝐴2 100 ) + (1 − | 𝐹𝑂𝐸𝑋 50 |) + (1 − 𝐾𝑆𝑃 100 ) (2.35) Trong đó Rank 4 trong khoảng từ 0 đến 6 tương ứng với từ kém đến tốt. 81 CHƯƠNG 3. KẾT QUẢ NGHIÊN CỨU VÀ BÀN LUẬN 3.1. Vận chuyển của các nhân phóng xạ trong khí quyển từ Fukushima đến Đông Nam Á 3.1.1. Thiết kế mô phỏng phát tán phóng xạ từ Fukushima đến Tây Thái Bình Dương và Đông Nam Á Tháng 3/2011, một loạt sự cố các lò phản ứng nghiêm trọng đã xảy ra tại NMĐHN Fukushima Daiichi (FNPP) phía Đông Bắc của Nhật Bản. Từ ngày 12/3/2011 đến hết tháng 4/2011, một lượng lớn các chất phóng xạ đã phát thải vào bầu khí quyển và phát tán trên bắc bán cầu do vận chuyển đối lưu nhanh bởi dòng xiết. Nhiều nguyên cứu đã thực hiện mô phỏng phát tán khí quyển và kiểm chứng các mô phỏng với dữ liệu quan trắc từ hàng trăm trạm quan trắc trên bắc bán cầu [5-9]. Trong khi đó, rất hiếm các mô phỏng vận chuyển nhân phóng xạ từ tai nạn Fukushima đến vùng nhiệt đới phía tây Thái Bình Dương (TWP) và Đông Nam Á (SEA) (khu vực rộng lớn trải dài từ phía Đông đến Tây Nam Nhật Bản). Huh và cộng sự (2013) đã thực hiện mô phỏng với mô hình đánh dấu WRF/Chem và cho thấy rằng các chùm phóng xạ được vận chuyển trong lớp biên hành tinh bởi gió mùa Đông Bắc hướng tới vùng TWP [39]. Tức là khu vực thấp của tầng đối lưu nơi bề mặt trái đất ảnh hưởng mạnh đến nhiệt độ, độ ẩm và gió thông qua sự di chuyển hỗn loạn của khối không khí. Trong nghiên cứu này, chúng tôi sử dụng mô hình phát tán hạt Lagrangian FLEXPART [33, 43] để mô phỏng vận chuyển trong khí quyển của các chất phóng xạ từ NMĐHN Fukushima đến TWP và SEA. Các kết quả được đánh giá với dữ liệu quan trắc trong khu vực. Chúng tôi đánh giá sự đóng góp của gió mùa Đông Bắc và các điều kiện khí tượng liên quan đến vận chuyển khu vực của chất phóng xạ vì điều này có thể cho biết khả năng tác động của phóng xạ phát thải từ các NMĐHN khác ở khu vực đông châu Á. Hiện có hơn một trăm lò phản ứng hạt nhân (NPR) đang hoạt động, đang xây dựng và trong kế hoạch xây 82 dựng ở phía Đông Bắc châu Á (PRIS, IAEA) [1]. Hàng chục NPR nằm cách biên giới phía bắc Việt Nam phạm vi 70 đến 500 km. Trong trường hợp xảy ra tai nạn hạt nhân, các nhân phóng xạ có thể được vận chuyển bởi gió Đông Bắc của gió mùa mùa đông Đông Á (EAWM) [74] đến các khu vực trong TWP và SEA. Do đó, một mô hình phát tán khí quyển được thử nghiệm thông qua tai nạn FNPP sẽ giúp thông báo mức độ sự cố của các tai nạn hạt nhân trong tương lai, cũng như cho biết cả việc theo dõi và biện pháp đối phó cho các tai nạn như vậy. Hình 3.1. Vị trí của các trạm quan trắc hoạt động trong thời gian tai nạn FNPP trong khu vực TWP và SEA (các điểm đỏ). Vị trí NMĐHN Fukushima điểm sao vàng. Mô hình phát tán hạt Lagrangian FLEXPART phiên bản 9.02 [33, 43] kết hớp với dữ liệu khí tượng CFSv2 từ Trung tâm Dự báo môi trường quốc gia Hoa Kỳ (NCEP) [56]. Dữ liệu khí tượng khả dụng ở 37 mức với độ phân giải không gian là 0.5°. Độ phân giải ngang của lưới đầu ra FLEXPART là 0,5×0,5° trên toàn cầu. 83 Các kết quả mô phỏng được so sánh với nồng độ phóng xạ của 131I và 137Cs quan trắc được tại 10 trạm quan trắc trong khu vực TWP và SEA (Hình 3.1), cụ thể là Fukuoka [16]; Okinawa; Guam; Tanay gần Manila (gọi là Manila); Tana Rata gần Kuala Lumpur (gọi là KL) [14]; Nankang ở mũi phía bắc của Đài Loan [39]; Hồng Kông [75]; Hà Nội; Đà Lạt; và Thành phố Hồ Chí Minh (TP.HCM) ở Việt Nam [15]. Các trạm quan trắc này trải dài hơn 5300 km về phía tây nam của FNPP, như trong Hình 3.1. Nồng độ phóng xạ của hai nhân phóng xạ này trải rộng trên ba bậc độ lớn, từ 103 µBq/m3 ở Guam và Fukuoka đến vài µBq/m3 tại Kuala Lumpur. Hình 3.2. Phát thải của 131I và 137Cs trong thời gian xảy ra tai nạn Fukushima theo Katata và cộng sự (2015). Đối với các kịch bản phát thải của 131I và 137Cs, số hạng nguồn từ nghiên cứu Katata và cộng sự (2015) [38] được sử dụng trong mô phỏng này (Hình 3.2). Các nhân phóng xạ được giả thiết phát tán đồng đều từ bề mặt lên độ cao 500m [5]. Tổng lượng phát thải của 131I và 137Cs ước tính bởi Katata và cộng sự (2015) lần lượt là 151 và 14.5 PBq. Tỷ lệ khí so với tổng của 131I là 0.5 theo Katata và cộng sự (2015). Các thông số đầu vào cho rơi lắng ướt của các hạt 131I và 137Cs như sau: hệ số scavenging A (1x10-4 s-1) tương ứng tốc độ mưa 1 mm/h và hệ số phụ thuộc mưa B (0.80) [5]. Đối với các thông số liên quan đến rơi lắng khô, có hai thiết lập đường kính trung bình hạt (dp) và mật độ hạt (ρp), cụ thể là dp = 0,6 µm, ρp = 2500 kg/m3 [76] trong mô phỏng I, và dp = 0,4 µm, ρp = 2500 kg/m3 trong mô phỏng II. 84 Các nhân phóng xạ phát tán trong khí quyển được mô phỏng từ 15h00 ngày 11/3/2011 đến 0h00 ngày 20/4/2011 UTC. Tổng gần 10 triệu hạt được phát thải trong mô phỏng. Mô phỏng được thiết kế chạy trên siêu máy tính PARAM- HUST của Trung tâm Tính toán hiệu năng cao, Trường Đại học Bách khoa Hà Nội. Siêu máy tính PARAM-HUST bao gồm 15 node với 2 chíp 8 core Intel Xenon-