近些年由于全球气候变暖，全球的海平面不断上升。目前的研究显示，造成海平面变化的原因主要有两个：一个是由于陆地质量的流入（比如两级冰川的融化，冰雪以径流的形式汇入海洋，总体上使得海洋的总质量产生变化），另一个是由于海洋面积巨大，会吸收全球范围内大量的热量，使得海水的温度和研究发生变化，进一步产生了比容海平面的变化。
计算比容海平面变化的基本原理见以下的专栏：
温盐海平面计算专题--01 - 哔哩哔哩 (bilibili.com)
温盐海平面计算专题--02 - 哔哩哔哩 (bilibili.com)
本专栏将进一步用matlab实现对比容海平面的计算，具体的程序代码如下：
% % % % % % % % % % read data
for i = 1:258SA(:,:,:,i) = salt{i,1};CT(:,:,:,i) = temp{i,1};disp(i)end% % % % % % % % %% cal sterictime_step = length(squeeze(CT(1,1,1,:)));[steric_height] = steric_height_calculation(CT,SA,dep,lat,lon,time_step);% % % % % % % % % % save dataload 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_stepfor 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]endend
endDEPTH = 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
endfor t = 1:time_stepsteric_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
欢迎交流学习！
