计算由于海洋温度和盐度变化产生的比容海平面变化

近些年由于全球气候变暖，全球的海平面不断上升。目前的研究显示，造成海平面变化的原因主要有两个：一个是由于陆地质量的流入（比如两级冰川的融化，冰雪以径流的形式汇入海洋，总体上使得海洋的总质量产生变化），另一个是由于海洋面积巨大，会吸收全球范围内大量的热量，使得海水的温度和盐度发生变化，进一步产生了比容海平面的变化。

计算比容海平面变化的基本原理见以下的专栏：

温盐海平面计算专题--01 - 哔哩哔哩 (bilibili.com)

温盐海平面计算专题--02 - 哔哩哔哩 (bilibili.com)

本专栏将进一步用matlab实现对比容海平面的计算，具体的程序代码如下：

% % % % % % % % % % read data
for i = 1:258
     SA(:,:,:,i) = salt{i,1};
     CT(:,:,:,i) = temp{i,1};
     disp(i)
 end
 % % % % % % % % %% cal steric
 time_step = length(squeeze(CT(1,1,1,:)));
 [steric_height] = steric_height_calculation(CT,SA,dep,lat,lon,time_step);
 % % % % % % % % % % save data
 load lon_lat_dep.mat
 [lon,lat] = meshgrid(lon,lat);
 ocean_height.lon = lon;
 ocean_height.lat = lat;
 ocean_height.rg  = steric_height;
 ocean_height.tt  = tt;

% % % % % % % % % % % time series
tp_gind=load('BS.txt');  
Tibet_ind=inpolygon(interpn(ocean_height.lon,1),interpn(ocean_height.lat,1),tp_gind(:,1),tp_gind(:,2));
area_scale=cal_grid_region(ocean_height);
for ii=1:size(ocean_height.rg,3)
    A = ocean_height.rg(:,:,ii);A(isnan(A))=0;
    CSR_rg(ii)=sum(sum(interpn(A',1).*Tibet_ind.*interpn(area_scale,1)))/sum(sum(Tibet_ind.*interpn(area_scale,1)));
    jgc_tt(ii)=ocean_height.tt(ii);
end
plot(jgc_tt,CSR_rg*1000) %% mm

其中使用到的函数 steric_height_calculation.m，需要注意的是，一般数据读取的温度是以开尔文为单位，而程序计算的是以摄氏度为单位，因此需要进行转换。本专栏使用的温度和盐度数据是MetOffice提供的再分析数据，需要自己提前读取。

function [steric_height] = steric_height_calculation(temperature,salinity,depth,lat,lon,time_step)
%--------------------------------------------------------------
% [steric_height] = steric_height_calculation(temperature,salinity,depth,lat,lon,time_step)
% This function is used to calculate the steric height.
% Note that the SEAWATER linrary version 3.2 by Lindsay Pender is 
% used in the code. 
%--------------------------------------------------------------
% input:
%  temperature(lon,lat,depth,time): temperature, unit: degree C
%  salinity(lon,lat,depth,time): salinity, unit: psu (PSS-78)
%  depth: depth of the ocean layer, unit: m
%  lat: latitude
%  lon: longitude
%  time_step: the number of time step of temperature/salinity 
%             ***  time_step = length(squeeze(temperature(1,1,1,:)))
%--------------------------------------------------------------
% output:
%  steric_height(lon,lat,time): steric height
%--------------------------------------------------------------
% calculate pressure from depth
pressure = zeros(length(depth),length(lat));
for k=1:length(depth)
    for j=1:length(lat)
        pressure(k,j)=sw_pres(depth(k),lat(j));%[db]
    end
end
pressure=pressure';
clear k j
rho = zeros(length(lon),length(lat),length(depth),time_step);
for t = 1:time_step
    for k = 1:length(depth)
    for j=1:length(lat)
    rho(:,j,k,t)=sw_dens(salinity(:,j,k,t),temperature(:,j,k,t)-273.15,pressure(j,k)); %[kg/m^3]
    end
    end
end

DEPTH = repmat(depth',length(lat),1);
steric_height = NaN(length(lon),length(lat),time_step);

rhobar = mean(rho,4,'omitnan'); % time-meaned rho
rho0_dep = squeeze(mean(mean(rhobar,1,'omitnan'),2,'omitnan')); % rho0 of each depth
dz =NaN(length(depth),1); 
dz(1) = abs(DEPTH(1,1)-0); 
dz(2:length(depth)) = abs(DEPTH(1,2:length(depth))-DEPTH(1,1:length(depth)-1));
DZ = NaN(length(lon),length(lat),length(depth)); rho0 = DZ;
for i = 1:length(lon)
    for j = 1:length(lat)
        DZ(i,j,:) = dz; %create DZ(lon,lat,depth) from dz(depth)
        rho0(i,j,:) = rho0_dep; % create rho0(lon,lat,depth) from rho0_dep(depth) 
    end
end

for t = 1:time_step
    steric_height(:,:,t) = -sum(DZ.*((squeeze(rho(:,:,1:length(depth),t))-rhobar)./rho0),3,'omitnan');
    disp(t)
end

需要注意的是，需要下载一个计算seawater各种物理量的库，下载链接如下：

https://github.com/ashao/matlab/tree/master/external/seawater

程序运行的结果：credit:Yuan et al., 2017.

全球某个月的比容海平面的空间分布，单位：m

渤海湾12个月的比容海平面空间分布：[-0.1,0.1]，单位：m

渤海比容海平面的时间变化：单位：mm

参考资料：

https://github.com/ashao/matlab/tree/master/external/seawater

Kuo, Y.-N., Lo, M.-H., Liang, Y.-C.,Tseng, Y.-H., & Hsu, C.-W. (2021). Terrestrial water storage anomalies emphasize interannual variations in global mean sea level during 1997–1998 and 2015–2016 El Niño events. Geophysical Research Letters, 48, e2021GL094104. https://doi.org/10.1029/2021GL094104

欢迎交流学习！