本项目利用 PLASIM-GENIE 气候系统模型输出的三组不同 CO₂ 浓度(控制实验 + 两倍/四倍 CO₂)全球表面温度数据,开展以下工作:
大气二氧化碳(CO₂)浓度的持续升高是全球变暖的核心驱动力之一(IPCC, 2021)。深入理解地表气温在不同 CO₂ 浓度条件下的变化规律与响应机制,对预测未来气候变化具有重要意义。本案例基于 Planet Simulator(PLASIM)模式模拟了多种 CO₂ 浓度情境,并通过合成分析、t 检验和蒙特卡洛检验评估地表气温的变化特征与统计显著性,并进一步探究温度变化背后的物理机制。
import xarray as xrimport pandas as pdimport numpy as npimport globimport osimport numpy as npimport matplotlib.pyplot as plt本研究基于 PlaSim–GENIE 耦合气候模型(中分辨率大气环流模型 PlaSim 与地球系统模型 GENIE 耦合)设计了以下 CO₂ 浓度情景实验:
说明:上述低浓度偏差不影响本研究主要结论(以 560 ppm 与 1120 ppm 相对于 280 ppm 的差值分析为主),但已在结果讨论中予以说明并进行稳健性检验。
def get_dataset(path,file_count,attribute_name): file = glob.glob(os.path.join(path, "*.nc")) nc=xr.open_dataset(file[0]) attribute_dataset= np.zeros((file_count,nc.latitude.shape[0], nc.longitude.shape[0])) for i in range(0,file_count): nc_name=file[i] nc=xr.open_dataset(nc_name) attribute_dataset[i]=nc[attribute_name].values mean_attribute=attribute_dataset.mean(0) return attribute_dataset,mean_attributepath='/home/mw/input/PLASIM_tem7051/PLASIM_surface_tem (2)/PLASIM_surface_tem/280ppm'file = glob.glob(os.path.join(path, "*.nc"))CO2_280ppm,mean_tem_280ppm=get_dataset(r'/home/mw/input/PLASIMsimple8604/PLASIM_surface_tem/PLASIM_surface_tem/280ppm',10,'puma_temperature_surface')CO2_560ppm,mean_tem_560ppm=get_dataset(r'/home/mw/input/PLASIMsimple8604/PLASIM_surface_tem/PLASIM_surface_tem/560ppm',10,'puma_temperature_surface')CO2_1120ppm,mean_tem_1120ppm=get_dataset(r'/home/mw/input/PLASIMsimple8604/PLASIM_surface_tem/PLASIM_surface_tem/1120ppm',10,'puma_temperature_surface')CO2_2805ppm,mean_tem_2805ppm=get_dataset(r'/home/mw/input/PLASIMsimple8604/PLASIM_surface_tem/PLASIM_surface_tem/280.5ppm',10,'puma_temperature_surface')CO2_5605ppm,mean_tem_5605ppm=get_dataset(r'/home/mw/input/PLASIMsimple8604/PLASIM_surface_tem/PLASIM_surface_tem/560.5ppm',10,'puma_temperature_surface')CO2_11205ppm,mean_tem_11205ppm=get_dataset(r'/home/mw/input/PLASIMsimple8604/PLASIM_surface_tem/PLASIM_surface_tem/1120.5ppm',10,'puma_temperature_surface')#查看返回数据集的维度np.mean(CO2_280ppm, axis=(1,2)).shapeimport matplotlib.pyplot as pltimport cartopy.crs as ccrsimport cartopy.feature as cfeatureimport cartopy.mpl.ticker as ctickerimport cmaps#计算温差tem_ano_560ppm=mean_tem_560ppm-mean_tem_280ppmtem_ano_1120ppm=mean_tem_1120ppm-mean_tem_280ppmprint('{}ppmCO2浓度下地表温差范围:min:{} max:{}'.format(560,tem_ano_560ppm.min(),tem_ano_560ppm.max()))print('{}ppmCO2浓度下地表温差范围:min:{} max:{}'.format(1120,tem_ano_1120ppm.min(),tem_ano_1120ppm.max()))print('{}ppmCO2浓度下全球增温:{:.2f}°'.format(560,tem_ano_560ppm.mean()))print('{}ppmCO2浓度下全球增温:{:.2f}°'.format(1120,tem_ano_1120ppm.mean()))#读取数据经纬度nc=xr.open_dataset(r'/home/mw/input/PLASIMsimple8604/PLASIM_surface_tem/PLASIM_surface_tem/280ppm/0051-12-30_plasim.nc')lat = nc.latitudelon = nc.longitude#绘图fig = plt.figure(figsize=(12,8))ax1 = fig.add_axes([0.1, 0.1, 0.8, 0.4],projection = ccrs.PlateCarree(central_longitude=180))ax1.set_extent([-180,180,-90,90], crs=ccrs.PlateCarree())ax1.add_feature(cfeature.LAKES, alpha=0.5)ax1.set_xticks(np.arange(-180,210,30), crs=ccrs.PlateCarree())ax1.set_yticks(np.arange(-90,120,30), crs=ccrs.PlateCarree())lon_formatter = cticker.LongitudeFormatter()lat_formatter = cticker.LatitudeFormatter()ax1.xaxis.set_major_formatter(lon_formatter)ax1.yaxis.set_major_formatter(lat_formatter)ax1.set_title('(a) Surface Temperature. anomaly in 580ppm ',loc='left',fontsize=18)c1 = ax1.contourf(lon,lat, tem_ano_560ppm,levels =np.arange(0,12,0.5) , extend = 'both', transform=ccrs.PlateCarree(), cmap=cmaps.NCV_jaisnd)ax1.add_feature(cfeature.COASTLINE.with_scale('50m')) ax2 = fig.add_axes([0.7, 0.1, 0.8, 0.4],projection = ccrs.PlateCarree(central_longitude=180))ax2.set_extent([-180,180,-90,90], crs=ccrs.PlateCarree())ax2.add_feature(cfeature.COASTLINE.with_scale('50m')) ax2.add_feature(cfeature.LAKES, alpha=0.5)ax2.set_xticks(np.arange(-180,210,30), crs=ccrs.PlateCarree())ax2.set_yticks(np.arange(-90,120,30), crs=ccrs.PlateCarree())ax2.xaxis.set_major_formatter(lon_formatter)ax2.yaxis.set_major_formatter(lat_formatter)ax2.set_title('(b) Surface Temperature. anomaly in 1120ppm ',loc='left',fontsize=18)c2 = ax2.contourf(lon,lat, tem_ano_1120ppm, zorder=0,levels =np.arange(0,12,0.5) , extend = 'both', transform=ccrs.PlateCarree(), cmap=cmaps.NCV_jaisnd)position=fig.add_axes([0.3, 0.02, 0.35, 0.025])fig.colorbar(c1,cax=position,orientation='horizontal',format='%d',) position=fig.add_axes([0.9, 0.02, 0.35, 0.025])fig.colorbar(c2,cax=position,orientation='horizontal',format='%d',) from scipy.stats.mstats import ttest_ind_,p_560ppm = ttest_ind(CO2_560ppm,CO2_280ppm)_,p_1120ppm = ttest_ind(CO2_1120ppm,CO2_280ppm)#查看不显著点个数print((p_560ppm>0.01).flatten().sum(0))print((p_1120ppm>0.01).flatten().sum(0))在合成分析结果上叠加显著性检验结果
fig = plt.figure(figsize=(12,8))ax1 = fig.add_axes([0.1, 0.1, 0.8, 0.4],projection = ccrs.PlateCarree(central_longitude=180))ax1.set_extent([-180,180,-90,90], crs=ccrs.PlateCarree())ax1.add_feature(cfeature.LAKES, alpha=0.5)ax1.set_xticks(np.arange(-180,210,30), crs=ccrs.PlateCarree())ax1.set_yticks(np.arange(-90,120,30), crs=ccrs.PlateCarree())lon_formatter = cticker.LongitudeFormatter()lat_formatter = cticker.LatitudeFormatter()ax1.xaxis.set_major_formatter(lon_formatter)ax1.yaxis.set_major_formatter(lat_formatter)ax1.set_title('(a) Surface Temperature. anomaly in 580ppm ',loc='left',fontsize=18)c1 = ax1.contourf(lon,lat, tem_ano_560ppm,levels =np.arange(0,12,0.5) , extend = 'both', transform=ccrs.PlateCarree(), cmap=cmaps.NCV_jaisnd)c1p = ax1.contourf(lon,lat,p_560ppm, levels =[0,0.05,1],hatches=['...', None],zorder=1,colors="none", transform=ccrs.PlateCarree())ax1.add_feature(cfeature.COASTLINE.with_scale('50m')) ax2 = fig.add_axes([0.7, 0.1, 0.8, 0.4],projection = ccrs.PlateCarree(central_longitude=180))ax2.set_extent([-180,180,-90,90], crs=ccrs.PlateCarree())ax2.add_feature(cfeature.COASTLINE.with_scale('50m')) ax2.add_feature(cfeature.LAKES, alpha=0.5)ax2.set_xticks(np.arange(-180,210,30), crs=ccrs.PlateCarree())ax2.set_yticks(np.arange(-90,120,30), crs=ccrs.PlateCarree())ax2.xaxis.set_major_formatter(lon_formatter)ax2.yaxis.set_major_formatter(lat_formatter)ax2.set_title('(b) Surface Temperature. anomaly in 1120ppm ',loc='left',fontsize=18)c2 = ax2.contourf(lon,lat, tem_ano_1120ppm, zorder=0,levels =np.arange(0,12,0.5) , extend = 'both', transform=ccrs.PlateCarree(), cmap=cmaps.NCV_jaisnd)c2p = ax2.contourf(lon,lat,p_1120ppm, levels =[0,0.05,1],hatches=['...', None],zorder=1,colors="none", transform=ccrs.PlateCarree())position=fig.add_axes([0.3, 0.02, 0.35, 0.025])fig.colorbar(c1,cax=position,orientation='horizontal',format='%d',) position=fig.add_axes([0.9, 0.02, 0.35, 0.025])fig.colorbar(c2,cax=position,orientation='horizontal',format='%d',) import numpy as npfrom scipy.stats import percentileofscore蒙特卡洛检验,基于矩阵打乱
def matrix_permutation_test(A, B, n_permutations=1000): """ 蒙特卡罗方法对两个矩阵进行检验,返回每个元素对应的p值 :param A: 第一个矩阵 :param B: 第二个矩阵 :param n_permutations: 随机次数 :return: 每个元素的p值 """ n_rows, n_cols = A.shape D = A - B # 计算差异矩阵 A=A.flatten() B=B.flatten() permuted_ranks = np.zeros((n_permutations, n_rows, n_cols)) for i in range(n_permutations): #生成随机序列,对A和B进行重排列 A_index=np.random.permutation(range(len(A.flatten()))) B_index=np.random.permutation(range(len(A.flatten()))) D_P =A[A_index]-B[B_index] # 计算每个置换下的差异向量 permuted_ranks[i] = D_P.reshape(n_rows, n_cols) # p_values = np.zeros((n_rows, n_cols)) for i in range(n_rows): for j in range(n_cols): p_values[i, j] =1-percentileofscore(permuted_ranks[:, i, j].flatten(),D[i, j]) / 100 # 计算p值 return p_values蒙特卡洛检验,基于序列打乱
def MK_test(A, B, n_permutations=1000): """ 蒙特卡罗方法对两个矩阵进行检验,返回每个元素对应的p值 :param A: 第一个时间序列 :param B: 第二个时间序列 :param n_permutations: 随机次数 :return: 每个元素的p值 """ n_length,n_rows, n_cols = A.shape # D = A - B # 计算差异矩阵 # A=A.flatten() # B=B.flatten() D=A.mean(0)-B.mean(0) permuted_ranks = np.zeros((n_permutations*n_permutations, n_rows, n_cols)) for i in range(n_rows): for j in range(n_cols): # A_P=np.tile(A[:,i,j],(1,n_permutations)) # B_P=np.tile(B[:,i,j],(1,n_permutations)) # D_P =np.random.permutation( A_P)-np.random.permutation( B_P) # permuted_ranks[:, i, j]= D_P for n in range(n_permutations): #生成随机序列,对A和B进行重排列 D_P =np.random.permutation(A[:,i,j])-np.random.permutation(B[:,i,j]) # 计算每个置换下的差异向量 permuted_ranks[n*n_length:(n+1)*n_length, i, j]= D_P p_values = np.zeros((n_rows, n_cols)) for i in range(n_rows): for j in range(n_cols): p_values[i, j] =1-percentileofscore(permuted_ranks[:, i, j],D[i, j]) / 100 # 计算p值 return p_values蒙特卡洛检验,基于序列打乱
def MK_test(A, B, n_permutations=1000): """ 蒙特卡罗方法对两个矩阵进行检验,返回每个元素对应的p值 :param A: 第一个时间序列 :param B: 第二个时间序列 :param n_permutations: 随机次数 :return: 每个元素的p值 """ n_length,n_rows, n_cols = A.shape # D = A - B # 计算差异矩阵 # A=A.flatten() # B=B.flatten() D=A.mean(0)-B.mean(0) permuted_ranks = np.zeros((n_permutations*n_permutations, n_rows, n_cols)) for i in range(n_rows): for j in range(n_cols): for n in range(n_permutations): #构造样本值数组 values=np.concatenate((A[:,i,j],B[:,i,j]),axis=0) #生成随机序列,对A和B进行重排列 D_P =np.random.permutation(values)[:10]-np.random.permutation(values)[:10] # 计算每个置换下的差异向量 permuted_ranks[n*n_length:(n+1)*n_length, i, j]= D_P p_values = np.zeros((n_rows, n_cols)) for i in range(n_rows): for j in range(n_cols): p_values[i, j] =1-percentileofscore(permuted_ranks[:, i, j],D[i, j]) / 100 # 计算p值 return p_valuesA=np.concatenate((CO2_560ppm[:,1,1],CO2_560ppm[:,1,1]),axis=0)np.random.permutation(A)[:10]p_560ppm=MK_test(CO2_560ppm,CO2_280ppm,n_permutations=500)p_1120ppm=MK_test(CO2_1120ppm,CO2_280ppm,n_permutations=500)sum(p_1120ppm<0.05).sum()sum(p_560ppm<0.05).sum()fig = plt.figure(figsize=(12,8))ax1 = fig.add_axes([0.1, 0.1, 0.8, 0.4],projection = ccrs.PlateCarree(central_longitude=180))ax1.set_extent([-180,180,-90,90], crs=ccrs.PlateCarree())ax1.add_feature(cfeature.LAKES, alpha=0.5)ax1.set_xticks(np.arange(-180,210,30), crs=ccrs.PlateCarree())ax1.set_yticks(np.arange(-90,120,30), crs=ccrs.PlateCarree())lon_formatter = cticker.LongitudeFormatter()lat_formatter = cticker.LatitudeFormatter()ax1.xaxis.set_major_formatter(lon_formatter)ax1.yaxis.set_major_formatter(lat_formatter)ax1.set_title('(a) Surface Temperature. anomaly in 580ppm ',loc='left',fontsize=18)c1 = ax1.contourf(lon,lat, tem_ano_560ppm,levels =np.arange(0,12,0.5) , extend = 'both', transform=ccrs.PlateCarree(), cmap=cmaps.NCV_jaisnd)c1p = ax1.contourf(lon,lat,p_560ppm, levels =[0,0.05,1],hatches=['...', None],zorder=1,colors="none", transform=ccrs.PlateCarree())ax1.add_feature(cfeature.COASTLINE.with_scale('50m')) ax2 = fig.add_axes([0.7, 0.1, 0.8, 0.4],projection = ccrs.PlateCarree(central_longitude=180))ax2.set_extent([-180,180,-90,90], crs=ccrs.PlateCarree())ax2.add_feature(cfeature.COASTLINE.with_scale('50m')) ax2.add_feature(cfeature.LAKES, alpha=0.5)ax2.set_xticks(np.arange(-180,210,30), crs=ccrs.PlateCarree())ax2.set_yticks(np.arange(-90,120,30), crs=ccrs.PlateCarree())ax2.xaxis.set_major_formatter(lon_formatter)ax2.yaxis.set_major_formatter(lat_formatter)ax2.set_title('(b) Surface Temperature. anomaly in 1120ppm ',loc='left',fontsize=18)c2 = ax2.contourf(lon,lat, tem_ano_1120ppm, zorder=0,levels =np.arange(0,12,0.5) , extend = 'both', transform=ccrs.PlateCarree(), cmap=cmaps.NCV_jaisnd)c2p = ax2.contourf(lon,lat,p_1120ppm, levels =[0,0.05,1],hatches=['...', None],zorder=1,colors="none", transform=ccrs.PlateCarree())position=fig.add_axes([0.3, 0.02, 0.35, 0.025])fig.colorbar(c1,cax=position,orientation='horizontal',format='%d',) position=fig.add_axes([0.9, 0.02, 0.35, 0.025])fig.colorbar(c2,cax=position,orientation='horizontal',format='%d',) 显著性检验结果表明,不同 CO₂ 浓度情景下的地表温度合成结果均具有统计显著性。t 检验显示,在 560 ppm 情景下共有 7 个格点未通过显著性检验,而在 1120 ppm 情景下,仅有 2 个格点不显著。进一步采用蒙特卡洛检验进行验证后发现,1120 ppm 情景下所有格点均通过显著性检验,而 560 ppm 情景下仍存在 2 个格点未达显著水平。
合成分析结果显示,相较于 280 ppm(工业化前基准情景),在 560 ppm 情景下全球地表年均温升为 2.98°C,1120 ppm 情景下则升高至 5.74°C。空间分布上,极地和亚洲地区是主要的增温热点;非洲、北美、南美和澳大利亚亦表现出显著的气温上升。整体来看,陆地和高纬地区的增温幅度明显高于海洋区域,体现出“极地放大”与“陆地-海洋温差”的典型特征。

#读取数据经纬度nc=xr.open_dataset(r'/home/mw/input/PLASIM_tem7051/PLASIM_surface_tem/PLASIM_surface_tem/1120.5ppm/0051-12-30_plasim.nc')lat = nc.latitudelon = nc.longitudeimport matplotlib.pyplot as pltfig=plt.figure(figsize=(4,4),dpi=150)ax=fig.add_axes([0,0,1,1])x=lon.valuesy1=tem_ano_560ppm.mean(0)y2=tem_ano_1120ppm.mean(0)y3=y2-y1ax.plot(x,y1,ls='-.',c='blue',fillstyle='none',label='560ppm')ax.plot(x,y2,ls=':',c='red',label='1120ppm')ax.set_ylabel('warming amplitude ($^\circ$C)')ax.set_xticks(np.arange(0, 420, 60))lon_formatter = cticker.LongitudeFormatter()ax.xaxis.set_major_formatter(lon_formatter)ax.set_xlabel('Longitude')ax.fill_between(x,y1,y2,where=(y1>y2),interpolate=True, facecolor='tab:blue',alpha=0.8)ax.fill_between(x,y1,y2,where=(y2>y1),interpolate=True, facecolor='tab:orange',alpha=0.6)plt.legend()fig=plt.figure(figsize=(4,4),dpi=150)ax=fig.add_axes([0,0,1,1])x=lat.valuesy1=tem_ano_560ppm.mean(1)y2=tem_ano_1120ppm.mean(1)ax.plot(x,y1,ls='-.',c='blue',fillstyle='none',label='560ppm')ax.plot(x,y2,ls=':',c='red',label='1120ppm')ax.set_ylabel('warming amplitude ($^\circ$C)')ax.set_xlabel('Latitude')ax.set_xticks(np.arange(-90, 120, 30))lat_formatter = cticker.LatitudeFormatter()ax.xaxis.set_major_formatter(lat_formatter)ax.fill_between(x,y1,y2,where=(y2>y1),interpolate=True, facecolor='tab:orange',alpha=0.6)plt.legend()


from scipy.stats.mstats import ttest_ind#计算不同浓度下的温差tem_ano_560ppm=mean_tem_560ppm-mean_tem_280ppmtem_ano_1120ppm=mean_tem_1120ppm-mean_tem_280ppmtem_ano_5605ppm=mean_tem_5605ppm-mean_tem_280ppmtem_ano_11205ppm=mean_tem_11205ppm-mean_tem_280ppm#显著性检验_,p_560ppm = ttest_ind(CO2_560ppm,CO2_280ppm)_,p_1120ppm = ttest_ind(CO2_1120ppm,CO2_280ppm)_,p_5605ppm = ttest_ind(CO2_5605ppm,CO2_280ppm)_,p_11205ppm = ttest_ind(CO2_11205ppm,CO2_280ppm)fig = plt.figure(figsize=(12,8))ax1 = fig.add_axes([0.1, 0.1, 0.8, 0.4],projection = ccrs.PlateCarree(central_longitude=180))ax1.set_extent([-180,180,-90,90], crs=ccrs.PlateCarree())ax1.add_feature(cfeature.LAKES, alpha=0.5)ax1.set_xticks(np.arange(-180,210,30), crs=ccrs.PlateCarree())ax1.set_yticks(np.arange(-90,120,30), crs=ccrs.PlateCarree())lon_formatter = cticker.LongitudeFormatter()lat_formatter = cticker.LatitudeFormatter()ax1.xaxis.set_major_formatter(lon_formatter)ax1.yaxis.set_major_formatter(lat_formatter)ax1.set_title('(c) Surface Temperature. anomaly in 580.5ppm ',loc='left',fontsize=18)c1 = ax1.contourf(lon,lat, tem_ano_5605ppm,levels =np.arange(0,12,0.5) , extend = 'both', transform=ccrs.PlateCarree(), cmap=cmaps.NCV_jaisnd)c1p = ax1.contourf(lon,lat,p_5605ppm, levels =[0,0.05,1],hatches=['...', None],zorder=1,colors="none", transform=ccrs.PlateCarree())ax1.add_feature(cfeature.COASTLINE.with_scale('50m')) ax2 = fig.add_axes([0.7, 0.1, 0.8, 0.4],projection = ccrs.PlateCarree(central_longitude=180))ax2.set_extent([-180,180,-90,90], crs=ccrs.PlateCarree())ax2.add_feature(cfeature.COASTLINE.with_scale('50m')) ax2.add_feature(cfeature.LAKES, alpha=0.5)ax2.set_xticks(np.arange(-180,210,30), crs=ccrs.PlateCarree())ax2.set_yticks(np.arange(-90,120,30), crs=ccrs.PlateCarree())ax2.xaxis.set_major_formatter(lon_formatter)ax2.yaxis.set_major_formatter(lat_formatter)ax2.set_title('(d) Surface Temperature. anomaly in 1120.5ppm ',loc='left',fontsize=18)c2 = ax2.contourf(lon,lat, tem_ano_11205ppm, zorder=0,levels =np.arange(0,12,0.5) , extend = 'both', transform=ccrs.PlateCarree(), cmap=cmaps.NCV_jaisnd)c2p = ax2.contourf(lon,lat,p_11205ppm, levels =[0,0.05,1],hatches=['...', None],zorder=1,colors="none", transform=ccrs.PlateCarree())position=fig.add_axes([0.3, 0.02, 0.35, 0.025])fig.colorbar(c1,cax=position,orientation='horizontal',format='%d',) position=fig.add_axes([0.9, 0.02, 0.35, 0.025])fig.colorbar(c2,cax=position,orientation='horizontal',format='%d',) print('{}ppmCO2浓度下全球增温:{:.2f}°'.format(560,tem_ano_5605ppm.mean()))print('{}ppmCO2浓度下全球增温:{:.2f}°'.format(1120,tem_ano_11205ppm.mean()))