利用Argo数据分别计算温度、盐度和温盐所造成的比容海平面变化

本文所用到的温盐数据集：IPRC（美国夏威夷大学国际太平洋研究中心）

Argo data products | Argo (ucsd.edu)https://argo.ucsd.edu/data/argo-data-products/

理论知识（相关计算公式）：

代码和工具包准备：

参照论文 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和下面博文计算由于海洋温度和盐度变化产生的比容海平面变化-CSDN博客中的详细描述，需下载计算海底压力、海水密度以及其他有关海水信息的“seawater”代码包，还需下载Yan-Ning Kuo编写的用于计算比容海平面变化的代码包（“steric_height_calculation.m”）

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

② steric_height_calculation: ynkuo/TWS_emphysize_GMSL_difference_in_9798_1516_ElNino_paper_codes: v1.0.0 (zenodo.org)https://zenodo.org/records/5144491

下面重点以IPRC数据集为例，分别计算由温度、盐度和温盐度造成的比热容变化

1、首先展示steric_height_calculation.m计算比容海平面变化的核心代码：

按照上述盐容和热容海平面变化的计算公式，计算热容时，保持盐度为平均盐度不变；同理，计算盐容时，保持温度为平均温度不变；故得到如下热容和盐容的求解代码：

%% ii 取值为 1：热容  2：盐容  3：比容
rho = zeros(length(lon),length(lat),length(depth),time_step);
salinity_0(:,:,:,1)=mean(salinity,4,'omitnan');temperature_0(:,:,:,1)=mean(temperature,4,'omitnan');
for t = 1:time_step
    for k = 1:length(depth)
    for j=1:length(lat)
        if ii==1
               rho(:,j,k,t)=sw_dens(salinity_0(:,j,k,1),temperature(:,j,k,t)-273.15,pressure(j,k));
        elseif ii==2
               rho(:,j,k,t)=sw_dens(salinity(:,j,k,t),temperature_0(:,j,k,1)-273.15,pressure(j,k));
        else
               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
end

2、值得注意的是

① IPRC和SIO的温度单位均为摄氏度，EN4为开尔文；所以在计算前两者比容变化时，不需要再对温度进行单位转换。即下面三项不需要减去273.15。

② 虽然SIO 给出变量名为“压强/压力”的变量信息，但建议还是利用sw_pres函数将其转换为某深度处的压力大小。

3、读取IPRC温盐数据，并绘制比容海平面变化全球趋势图与盐容、热容以及比容在太平洋某处的时间序列图

% read data
clear;
addpath E:\Data\Temp_sal\Argo;
addpath E:\Data\Temp_sal\Argo\function;
full_path = 'E:\Data\Temp_sal\Argo\IPRC\argo_2005-2020_grd.nc\argo_2005-2020_grd.nc';
info_IPRC = ncinfo(full_path);
TEMP=ncread(full_path,'TEMP');
SLAT=ncread(full_path,'SALT');
dep=ncread(full_path,'LEVEL'); 
lat=ncread(full_path,'LATITUDE');
lon=ncread(full_path,'LONGITUDE');

 %  cal steric
 time_step = length(squeeze(TEMP(1,1,1,:)));
 type = {'_T';'_S';'_'};
 %%最后一位数字 1：热容  2：盐容  3：比容
 tic
 for i=1:3
 [steric_height] = steric_height_calculation(TEMP,SLAT,dep,lat,lon,time_step,i);
 steric_height_avg=mean(steric_height,3);
   for ii=1:time_step
      temp(:,:)=steric_height(:,:,ii);
      SH(:,:,ii)=flipud((temp-steric_height_avg)');
   end
   eval(['SH' type{i,1} '=SH;']);
 end
 toc

FileNameTime05={'2005-01','2005-02','2005-03','2005-04','2005-05','2005-06','2005-07','2005-08','2005-09','2005-10','2005-11','2005-12'};
FileNameTime06={'2006-01','2006-02','2006-03','2006-04','2006-05','2006-06','2006-07','2006-08','2006-09','2006-10','2006-11','2006-12'};
FileNameTime07={'2007-01','2007-02','2007-03','2007-04','2007-05','2007-06','2007-07','2007-08','2007-09','2007-10','2007-11','2007-12'};
FileNameTime08={'2008-01','2008-02','2008-03','2008-04','2008-05','2008-06','2008-07','2008-08','2008-09','2008-10','2008-11','2008-12'};
FileNameTime09={'2009-01','2009-02','2009-03','2009-04','2009-05','2009-06','2009-07','2009-08','2009-09','2009-10','2009-11','2009-12'};
FileNameTime10={'2010-01','2010-02','2010-03','2010-04','2010-05','2010-06','2010-07','2010-08','2010-09','2010-10','2010-11','2010-12'};
FileNameTime11={'2011-01','2011-02','2011-03','2011-04','2011-05','2011-06','2011-07','2011-08','2011-09','2011-10','2011-11','2011-12'};
FileNameTime12={'2012-01','2012-02','2012-03','2012-04','2012-05','2012-06','2012-07','2012-08','2012-09','2012-10','2012-11','2012-12'};
FileNameTime13={'2013-01','2013-02','2013-03','2013-04','2013-05','2013-06','2013-07','2013-08','2013-09','2013-10','2013-11','2013-12'};
FileNameTime14={'2014-01','2014-02','2014-03','2014-04','2014-05','2014-06','2014-07','2014-08','2014-09','2014-10','2014-11','2014-12'};
FileNameTime15={'2015-01','2015-02','2015-03','2015-04','2015-05','2015-06','2015-07','2015-08','2015-09','2015-10','2015-11','2015-12'};
FileNameTime16={'2016-01','2016-02','2016-03','2016-04','2016-05','2016-06','2016-07','2016-08','2016-09','2016-10','2016-11','2016-12'};
FileNameTime17={'2017-01','2017-02','2017-03','2017-04','2017-05','2017-06','2017-07','2017-08','2017-09','2017-10','2017-11','2017-12'};
FileNameTime18={'2018-01','2018-02','2018-03','2018-04','2018-05','2018-06','2018-07','2018-08','2018-09','2018-10','2018-11','2018-12'};
FileNameTime19={'2019-01','2019-02','2019-03','2019-04','2019-05','2019-06','2019-07','2019-08','2019-09','2019-10','2019-11','2019-12'};
FileNameTime20={'2020-01','2020-02','2020-03','2020-04'};

FileNameTime=[FileNameTime05,FileNameTime06,FileNameTime07,FileNameTime08,FileNameTime09,FileNameTime10,...
              FileNameTime11,FileNameTime12,FileNameTime13,FileNameTime14,FileNameTime15,FileNameTime16,FileNameTime17,...
              FileNameTime18,FileNameTime19,FileNameTime20];

num_file=size(FileNameTime,2);
FileNameTimeChar=char(FileNameTime);
int_year=zeros(num_file,1);int_month=zeros(num_file,1);
for i=1:num_file
    int_year(i)=str2double(FileNameTimeChar(i,1:4));
    int_month(i)=str2double(FileNameTimeChar(i,6:7));
end
time=int_year+(int_month-0.5)/12;

[ Amp_T, ~, ~, ~,  ~,  ~,  ~,  ~, Trend_T] = gmt_harmonic(time,[],SH_T);
[ Amp_S, ~, ~, ~,  ~,  ~,  ~,  ~, Trend_S] = gmt_harmonic(time,[],SH_S);
[ Amp_, ~, ~, ~,  ~,  ~,  ~,  ~, Trend_] = gmt_harmonic(time,[],SH_);

% time series
mask='E:\Data\Basin\Pacific_rectangle.vec';  rows=5;
Amp_T_serie=gmt_grid2serie(Amp_T,mask,'line',rows);
Amp_S_serie=gmt_grid2serie(Amp_S,mask,'line',rows);
Amp_serie=gmt_grid2serie(Amp_,mask,'line',rows);

Tre_T_serie=gmt_grid2serie(Trend_T,mask,'line',rows);
Tre_S_serie=gmt_grid2serie(Trend_S,mask,'line',rows);
Tre_serie=gmt_grid2serie(Trend_,mask,'line',rows);

SH_T_serie=gmt_grid2serie(SH_T,mask,'line',rows);
SH_S_serie=gmt_grid2serie(SH_S,mask,'line',rows);
SH_serie=gmt_grid2serie(SH_,mask,'line',rows);

close all;
grid_1=SH(:,:,10);
gmt_grid2map(Trend_.*1000,-10,10,1,0,'mm/year',['Trend IPRC'],20);

 close all;
 plot(time,SH_T_serie,'r','Linewidth',1);hold on; %% mm
 plot(time,SH_S_serie,'g','Linewidth',1);
 plot(time,SH_serie,'--','color',[0 0 1],'Linewidth',1);
l1=legend('热容','盐容','比热容');
set(l1,'box','off',"FontSize",10);

function [steric_height] = steric_height_calculation(temperature,salinity,depth,lat,lon,time_step,ii)
%--------------------------------------------------------------
%% ii 1:热容、2:盐容、3:比容

% 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
%--------------------------------------------------------------
addpath E:\Data\Temp_sal\matlab-master\matlab-master\external\seawater;
% 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);
salinity_0(:,:,:,1)=mean(salinity,4,'omitnan');temperature_0(:,:,:,1)=mean(temperature,4,'omitnan');
for t = 1:time_step
    for k = 1:length(depth)
    for j=1:length(lat)
        if ii==1
               rho(:,j,k,t)=sw_dens(salinity_0(:,j,k,1),temperature(:,j,k,t),pressure(j,k));
        elseif ii==2
               rho(:,j,k,t)=sw_dens(salinity(:,j,k,t),temperature_0(:,j,k,1),pressure(j,k));
        else
               rho(:,j,k,t)=sw_dens(salinity(:,j,k,t),temperature(:,j,k,t),pressure(j,k)); %[kg/m^3]
        end
    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

如有任何问题，尽情批评指正！

PS：感谢“我是水怪的哥”博文的启发，还有引领我进入GRACE海洋研究的贵人，再次一并致谢！

参考文献：

① 2021 杨元元（博）基于卫星重力、卫星测高和海洋温盐数据的海平面收支及深海冷暖研究

② 2021 Kuo Yan-Ning Terrestrial Water Storage Anomalies Emphasize Interannual Variations in Global Mean Sea Level During 1997–1998 and 2015–2016 El Niño Events

③ 2021 王奉伟（博）基于卫星重力的质量海平面变化及其闭合度研究

④ 2022 李杨（硕）联合时变重力数据和卫星测高数据反演地表质量迁移及全球海平面变化

⑤ 2013 江敏（博）重海平面变化及其成因的空间大地测量监测与分析

⑥ 2012 文江汉联合Argo浮标、卫星测高和GRACE数据研究海平面变化