% ============================================================== % Combined Figure: Volcanic Residual Pressure + Crustal Thickness Models % Left: Residual pressure contour + average crustal models % Right: thick_min / thick_ave / thick_max three-line models % Horizontal colorbar placed between two subplots at bottom clear; clc; close all; %% Part 1: Volcano Physics Model (Left Plot Data) % Fixed constants E_star = 23.7; % Deformation modulus E*, GPa g = 3.72; % Mars gravitational acceleration, m/s² gEarth = 9.8; % Earth gravitational acceleration, m/s² rhoEarth = 2830; % Average Earth crustal density, kg/m³ mu0 = 151.4; % Pure magma viscosity μ₀, Pa·s phi = 0.55; % Crystal volume fraction φ phi0 = 0.64; % Viscosity constant φ₀ C = 3; % Viscosity crystal constant C kappa = 3e-7; % Magma thermal diffusivity κ, m²/s c = 1600; % Magma heat capacity c, J/(kg·K) dT_dx = 13e-3; % Wall rock temperature gradient, K/m L = 414e3; % Magma latent heat L, J/kg H = 11.65 * 1000; % Vertical extent of deep magma chamber H, m rho_M = 3462; % Mantle density ρ_M, kg/m³ sigma_d = 0.11 * 1e6; % Differential stress σ_d, Pa q = 20e-3; % Reserved parameter % Adjustable parameter rho_m = 3000; % Magma-crystal density ρ_m, kg/m³ % Variable ranges D_km = linspace(20, 60, 100); % Crustal thickness D, 20~60 km rho = linspace(2500, 3100, 100);% Crustal density ρ, 2500~3100 kg/m³ %% Precompute constant terms exponent_mu = C * phi0 / (1 - phi0); mu_m = mu0 * ((1 - phi) / (1 - phi/phi0))^exponent_mu; P_e = (rho_M - rho_m) * g * H; %% 2D grid calculation [rho_grid, D_grid] = meshgrid(rho, D_km); % Step 1: Young's modulus E E_GPa = E_star + (g/gEarth) .* (rho_grid ./ rhoEarth) .* D_grid.^0.4; E_Pa = E_GPa * 1e9; % Constants for ΔP_crit formula term1 = sqrt(kappa * mu_m / pi); term2 = c * dT_dx / L; % Step 2: Magma chamber overpressure P_o D_m = D_grid * 1000; P_o = P_e + (rho_grid - rho_m) .* g .* D_m + sigma_d; % Step 3: Critical overpressure ΔP_crit inside = term1 * term2 * E_Pa.^2; deltaP_crit = 3.5 * inside.^(2/5); % Step 4: Residual pressure P_r, in MPa Pr_MPa = (P_o - deltaP_crit) / 1e6; %% Part 2: Common Model Definitions model_names = {'DWThotCrust1r','EH45ThotCrust2r','SANAK','DWAK','YOTHotRc1810',... 'EH45Tcold','LFAK','TAYAK','ZGDW','Khan2022'}; model_colors = { [ 0 114 189]/255, % 1 DWThotCrust1r [119 172 48]/255, % 2 EH45ThotCrust2r [126 47 142]/255, % 3 SANAK [216 86 155]/255, % 4 DWAK [237 177 32]/255, % 5 YOTHotRc1810 [217 83 25]/255, % 6 EH45Tcold [162 20 47]/255, % 7 LFAK [139 69 19]/255, % 8 TAYAK [112 112 112]/255, % 9 ZGDW [140 120 34]/255 % 10 Khan2022 }; fit_order = 2; n_fit_pts = 200; %% Part 3: Left Plot: Average Crustal Models % Model raw data rhoc_DWThotCrust1r = [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000,3050]; thick_DWThotCrust1r = [32.005,31.610,31.169,30.673,30.111,29.469,... 28.734,27.872,26.855,25.637,24.150]; rhoc_EH45ThotCrust2r = [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000,3050]; thick_EH45ThotCrust2r = [31.940,31.527,31.064,30.541,29.945,29.267,... 28.474,27.542,26.433,25.089,23.429]; rhoc_SANAK = [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000]; thick_SANAK = [31.849,31.425,30.949,30.410,29.795,29.093,... 28.271,27.302,26.144,24.734]; rhoc_DWAK = [2550,2600,2650,2700,2750,2800,2850,2900,2950]; thick_DWAK = [31.505,31.022,30.474,29.850,29.126,28.278,... 27.274,26.066,24.584]; rhoc_YOTHotRc1810 = [2550,2600,2650,2700,2750,2800,2850,2900]; thick_YOTHotRc1810 = [31.726,31.239,30.689,30.051,29.309,28.435,... 27.391,26.122]; rhoc_EH45Tcold = [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000,3050]; thick_EH45Tcold = [31.829,31.409,30.938,30.405,29.798,29.105,28.296,... 27.346,26.212,24.839,23.139]; rhoc_LFAK = [2550,2600,2650,2700,2750,2800,2850,2900,2950]; thick_LFAK = [31.467,30.980,30.426,29.796,29.063,28.205,... 27.187,25.960,24.454]; rhoc_TAYAK = [2550,2600,2650,2700,2750,2800,2850,2900,2950]; thick_TAYAK = [31.467,30.980,30.426,29.796,29.063,28.205,... 27.187,25.960,24.454]; rhoc_ZG_DW = [2550,2600,2650,2700,2750,2800,2850,2900,2950]; thick_ZG_DW = [31.622,31.152,30.618,30.009,29.311,28.492,... 27.523,26.361,24.941]; rhoc_Khan2022 = [2550,2600,2650,2700,2750,2800,2850,2900,2950]; thick_Khan2022 = [31.806,31.343,30.817,30.221,29.527,28.716,... 27.755,26.599,25.184]; models = { rhoc_DWThotCrust1r, thick_DWThotCrust1r; rhoc_EH45ThotCrust2r, thick_EH45ThotCrust2r; rhoc_SANAK, thick_SANAK; rhoc_DWAK, thick_DWAK; rhoc_YOTHotRc1810, thick_YOTHotRc1810; rhoc_EH45Tcold, thick_EH45Tcold; rhoc_LFAK, thick_LFAK; rhoc_TAYAK, thick_TAYAK; rhoc_ZG_DW, thick_ZG_DW; rhoc_Khan2022, thick_Khan2022; }; % Fit curves models_fit_left = cell(size(models,1),2); h_model = gobjects(size(models,1),1); for i = 1:size(models,1) rhoc_i = models{i,1}; thick_i = models{i,2}; p_coeff = polyfit(rhoc_i, thick_i, fit_order); rhoc_fit = linspace(min(rhoc_i), max(rhoc_i), n_fit_pts); thick_fit = polyval(p_coeff, rhoc_fit); models_fit_left{i,1} = rhoc_fit; models_fit_left{i,2} = thick_fit; end % Calculate intersection with Pr=0 line C0 = contourc(rho, D_km, Pr_MPa, [0 0]); npts0 = C0(2,1); rho0 = C0(1, 2:1+npts0); D0 = C0(2, 2:1+npts0); rho_intersect_all = []; D_intersect_all = []; for i = 1:size(models_fit_left,1) rhoc_i = models_fit_left{i,1}; thick_i = models_fit_left{i,2}; D_pr0_interp = interp1(rho0, D0, rhoc_i, 'linear', 'extrap'); diff_val = thick_i - D_pr0_interp; idx_cross = find(diff_val(1:end-1) .* diff_val(2:end) <= 0); if ~isempty(idx_cross) for k = 1:length(idx_cross) idx = idx_cross(k); rho1 = rhoc_i(idx); rho2 = rhoc_i(idx+1); d1 = diff_val(idx); d2 = diff_val(idx+1); rho_cross = rho1 - d1 * (rho2 - rho1) / (d2 - d1); D_cross = interp1(rho0, D0, rho_cross, 'linear','extrap'); rho_intersect_all = [rho_intersect_all, rho_cross]; D_intersect_all = [D_intersect_all, D_cross]; end end end [rho_min_perm,idx_min] = min(rho_intersect_all); [rho_max_perm,idx_max] = max(rho_intersect_all); D_at_rhoMin = D_intersect_all(idx_min); D_at_rhoMax = D_intersect_all(idx_max); %% Part 4: Right Plot: Three-Line Thickness Data line_style = {'-'}; % thick_min data models_min = { [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000,3050], ... [21.008334,20.05043,18.984803,17.792591,16.449564,14.925357,13.185777,11.170555,8.814933,6.024796,2.667763]; [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000,3050], ... [20.460974,19.439924,18.300242,17.019789,15.57074,13.922418,12.020333,9.806839,7.198766,4.080368,0.285617]; [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000], ... [20.244296,19.197812,18.028053,16.711622,15.219276,13.518152,11.550432,9.254264,6.540005,3.282119]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [18.840438,17.620985,16.244742,14.683762,12.889109,10.809213,8.370282,5.470673,1.966311]; [2550,2600,2650,2700,2750,2800,2850,2900], ... [18.213124,16.911532,15.440648,13.756153,11.812952,9.546652,6.869481,3.658689]; [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000,3050], ... [20.417242,19.391119,18.245439,16.957765,15.500139,13.841113,11.925911,9.695891,7.066594,3.920321,0.087953]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [19.262196,18.096211,16.784092,15.296508,13.600908,11.639588,9.350897,6.645536,3.398413]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [18.759974,17.530168,16.141489,14.565437,12.752184,10.649027,8.180486,5.242388,1.686729]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [19.031045,17.835827,16.488655,14.958775,13.211579,11.185772,8.81551,6.00491,2.61878]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [18.877662,17.662981,16.292484,14.738969,12.953081,10.884198,8.459372,5.578209,2.098615]; }; % thick_ave data models_ave = { [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000,3050], ... [47.51838,48.001594,48.543314,49.154836,49.850596,50.649281,51.575754,52.66297,53.95698,55.523034,57.457104]; [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000,3050], ... [47.909778,48.433721,49.023542,49.692505,50.457696,51.341648,52.373998,53.595762,55.064368,56.86311,59.117574]; [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000], ... [47.969612,48.504557,49.107769,49.793219,50.578951,51.488895,52.554675,53.820299,55.347841,57.22808]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [48.728368,49.361489,50.083247,50.913786,51.879435,53.016251,54.374199,56.02481,58.074242]; [2550,2600,2650,2700,2750,2800,2850,2900], ... [49.6115,50.327743,51.148377,52.097724,53.208826,54.526893,56.115738,58.068449]; [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000,3050], ... [47.875572,48.398855,48.988111,49.656669,50.421685,51.305836,52.338962,53.562395,55.03409,56.838227,59.101966]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [48.558504,49.164771,49.853626,50.643186,51.557446,52.628127,53.899361,55.433362,57.321081]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [48.747005,49.384444,50.11588,50.948939,51.923333,53.071582,54.444807,56.116366,58.195452]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [48.726877,49.352575,50.064796,50.882815,51.832301,52.947323,54.275563,55.88469,57.874449]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [49.128122,49.785374,50.534362,51.395879,52.39704,53.574981,54.98109,56.688838,58.807071]; }; % thick_max data models_max = { [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000,3050], ... [77.296046,79.613713,82.248263,85.239595,88.671193,92.686188,97.398102,103.018803,109.824868,118.211841,128.806991]; [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000,3050], ... [78.286886,80.777832,83.597747,86.816782,90.52679,94.878936,100.012588,106.165612,113.783284,123.32468,135.599298]; [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000], ... [78.693727,81.247017,84.14285,87.455842,91.282738,95.786004,101.114664,107.52596,115.445321,125.485678]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [81.346992,84.329068,87.753483,91.731922,96.402232,101.992019,108.766184,117.148193,127.837924]; [2550,2600,2650,2700,2750,2800,2850,2900], ... [82.515457,85.702599,89.389129,93.689226,98.780251,104.931419,112.470356,121.930023]; [2550,2600,2650,2700,2750,2800,2850,2900,2950,3000,3050], ... [78.364763,80.866979,83.700572,86.936721,90.667397,95.046755,100.214925,106.413624,114.08656,123.714836,136.117516]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [80.540446,83.388863,86.647333,90.411937,94.818735,100.057578,106.361638,114.093831,123.846807]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [81.502137,84.510226,87.966989,91.986342,96.709632,102.368785,109.236551,117.748943,128.628518]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [80.978576,83.899817,87.248748,91.126386,95.679067,101.107972,107.666434,115.749376,126.006426]; [2550,2600,2650,2700,2750,2800,2850,2900,2950], ... [81.242607,84.210939,87.618487,91.576592,96.217461,101.772581,108.500917,116.820961,127.41924]; }; % Fit right plot curves fit_min = cell(10, 2); fit_ave = cell(10, 2); fit_max = cell(10, 2); rho_intersect_right = zeros(10, 1); p_min = cell(10, 1); p_ave = cell(10, 1); p_max = cell(10, 1); for i = 1:10 rhoc = models_min{i,1}; thick = models_min{i,2}; p = polyfit(rhoc, thick, fit_order); p_min{i} = p; rhoc_fit = linspace(min(rhoc), max(rhoc)+150, n_fit_pts); thick_fit = polyval(p, rhoc_fit); fit_min{i,1} = rhoc_fit; fit_min{i,2} = thick_fit; roots_all = roots(p); roots_real = roots_all(imag(roots_all)==0); roots_valid = roots_real(roots_real > min(rhoc)); rho_intersect_right(i) = max(roots_valid); rhoc = models_ave{i,1}; thick = models_ave{i,2}; p = polyfit(rhoc, thick, fit_order); p_ave{i} = p; rhoc_fit = linspace(min(rhoc), max(rhoc)+150, n_fit_pts); thick_fit = polyval(p, rhoc_fit); fit_ave{i,1} = rhoc_fit; fit_ave{i,2} = thick_fit; rhoc = models_max{i,1}; thick = models_max{i,2}; p = polyfit(rhoc, thick, fit_order); p_max{i} = p; rhoc_fit = linspace(min(rhoc), max(rhoc)+150, n_fit_pts); thick_fit = polyval(p, rhoc_fit); fit_max{i,1} = rhoc_fit; fit_max{i,2} = thick_fit; end rho_band_left = min(rho_intersect_right); rho_band_right = max(rho_intersect_right); n_models = length(model_names); thick_min_rhoMin = zeros(n_models, 1); thick_ave_rhoMin = zeros(n_models, 1); thick_max_rhoMin = zeros(n_models, 1); thick_min_rhoMax = zeros(n_models, 1); thick_ave_rhoMax = zeros(n_models, 1); thick_max_rhoMax = zeros(n_models, 1); for i = 1:n_models rho_val = rho_intersect_all(i); thick_min_rhoMin(i) = polyval(p_min{i}, rho_val); thick_ave_rhoMin(i) = polyval(p_ave{i}, rho_val); thick_max_rhoMin(i) = polyval(p_max{i}, rho_val); rho_val = rho_intersect_right(i); thick_min_rhoMax(i) = polyval(p_min{i}, rho_val); thick_ave_rhoMax(i) = polyval(p_ave{i}, rho_val); thick_max_rhoMax(i) = polyval(p_max{i}, rho_val); end fprintf('\n'); fprintf('======================================================================================================================\n'); fprintf(' Crustal Thickness Summary Table (Two Characteristic Densities)\n'); fprintf('======================================================================================================================\n'); fprintf('%-22s | %12s | %10s %10s %10s | %12s | %10s %10s %10s\n', ... 'Model Name', 'Min Density', 'thick_min', 'thick_ave', 'thick_max', ... 'Max Density', 'thick_min', 'thick_ave', 'thick_max'); fprintf(' | (kg/m³) | (km) (km) (km) | (kg/m³) | (km) (km) (km) \n'); fprintf('----------------------|-------------|------------------------------|-------------|-----------------------------\n'); for i = 1:n_models fprintf('%-22s | %12.2f | %10.2f %10.2f %10.2f | %12.2f | %10.2f %10.2f %10.2f\n', ... model_names{i}, rho_intersect_all(i), ... thick_min_rhoMin(i), thick_ave_rhoMin(i), thick_max_rhoMin(i), ... rho_intersect_right(i), ... thick_min_rhoMax(i), thick_ave_rhoMax(i), thick_max_rhoMax(i)); end fprintf('======================================================================================================================\n\n'); T = table(); T.ModelName = model_names'; T.Min_Permissible_Crustal_Density = rho_intersect_all'; T.MinDensity_thick_min = thick_min_rhoMin; T.MinDensity_thick_ave = thick_ave_rhoMin; T.MinDensity_thick_max = thick_max_rhoMin; T.Max_Permissible_Crustal_Density = rho_intersect_right; T.MaxDensity_thick_min = thick_min_rhoMax; T.MaxDensity_thick_ave = thick_ave_rhoMax; T.MaxDensity_thick_max = thick_max_rhoMax; % Write the CSV file to the current working directory writetable(T, 'CrustalThickness_Summary.csv'); fprintf('>>> The table has been exported as a file：CrustalThickness_Summary.csv (Located in the current MATLAB working directory) \n\n'); %% Original console output (retained) fprintf('\nLeft plot shaded band bounds:\n'); fprintf(' Left boundary = %.2f kg/m³ (thickness = %.2f km)\n', rho_min_perm, D_at_rhoMin); fprintf(' Right boundary = %.2f kg/m³ (thickness = %.2f km)\n', rho_max_perm, D_at_rhoMax); fprintf('\n\n'); fprintf('\nRight plot shaded band bounds:\n'); fprintf(' Left boundary = %.2f kg/m³\n', rho_band_left); fprintf(' Right boundary = %.2f kg/m³\n', rho_band_right); %% Part 5: Combined Plotting % Create wide figure figure('Color','w','Position',[100 100 1300 800]); % --- Left Subplot: Residual Pressure Contour --- ax1 = subplot(1,2,1); hold on; grid on; box on; % Custom colormap n_col = 1024; c_blue = [0.55, 0.78, 1.00]; c_white = [1.00, 1.00, 1.00]; c_brown = [0.93, 0.82, 0.65]; cmap = [interp1([0 0.75 1],[c_blue(1) c_white(1) c_brown(1)],linspace(0,1,n_col))', ... interp1([0 0.75 1],[c_blue(2) c_white(2) c_brown(2)],linspace(0,1,n_col))', ... interp1([0 0.75 1],[c_blue(3) c_white(3) c_brown(3)],linspace(0,1,n_col))']; colormap(ax1, cmap); % Filled residual pressure contour contourf(rho_grid, D_grid, Pr_MPa, 200, 'LineStyle','none'); caxis(ax1, [-60 20]); % Plot Pr=0 equilibrium line h_pr0 = plot(rho0, D0, 'Color','k', 'LineWidth', 2); % Plot 10 average crustal model curves for i = 1:size(models_fit_left,1) h_model(i) = plot(models_fit_left{i,1}, models_fit_left{i,2}, ... 'Color', model_colors{i}, 'LineWidth', 1.5); end % Permissible density shaded band y_lim_left = [20 35]; patch([rho_min_perm, rho_max_perm, rho_max_perm, rho_min_perm], ... [y_lim_left(1), y_lim_left(1), y_lim_left(2), y_lim_left(2)], ... [0.3 0.7 1.0], 'FaceAlpha', 0.3, 'EdgeColor', 'none'); text(rho_min_perm -10, 29.3, 'Minimum permissible density','Color', [0 0 0], ... 'FontSize', 10, 'FontName', 'Arial','FontWeight', 'bold', 'Rotation', 90); % Axis settings xlim([2550 3100]); ylim(y_lim_left); xticks(2600:100:3000); yticks(20:5:35); xlabel('Crustal density (kg m^{-3})', 'FontSize', 16, 'FontName','Arial'); ylabel('Global crustal thickness at the nakhlite ejection site (km)', 'FontSize', 16, 'FontName','Arial'); pbaspect(ax1, [1 1 1]); % Legend at bottom-left corner legend([h_model; h_pr0], [model_names, 'P_r=0 equilibrium'], ... 'Location','southwest', 'FontSize', 11, 'FontName','Arial', 'Box','on'); % Axis styling set(ax1, 'FontSize', 12, 'FontName','Arial', 'LineWidth', 1.0, ... 'TickDir','out', 'TickLength',[0.01 0.01]); % Right Subplot: Three-Line Thickness Models ax2 = subplot(1,2,2); hold on; grid off; box on; % Shaded vertical band y_lim_right = [0 140]; patch([rho_band_left, rho_band_right, rho_band_right, rho_band_left], ... [y_lim_right(1), y_lim_right(1), y_lim_right(2), y_lim_right(2)], ... [0.93 0.82 0.65], 'FaceAlpha', 0.3, 'EdgeColor', 'none'); text(rho_band_right - 20, 70, 'Maximum permissible density','Color', [0 0 0], ... 'FontSize', 10, 'FontName', 'Arial','FontWeight', 'bold','Rotation', 90); % Plot three lines per mode for i = 1:10 color = model_colors{i}; plot(fit_min{i,1}, fit_min{i,2}, 'Color', color, 'LineStyle', line_style{1}, 'LineWidth', 1.5); plot(fit_ave{i,1}, fit_ave{i,2}, 'Color', color, 'LineStyle', line_style{1}, 'LineWidth', 1.5); plot(fit_max{i,1}, fit_max{i,2}, 'Color', color, 'LineStyle', line_style{1}, 'LineWidth', 1.5); end % Axis settings xlim([2550 3110]); ylim(y_lim_right); xticks(2600:100:3100); yticks(0:20:140); xlabel('Crustal density (kg m^{-3})', 'FontSize', 16, 'FontName','Arial'); ylabel('Global crustal thickness (km)', 'FontSize', 16, 'FontName','Arial'); pbaspect(ax2, [1 1 1]); text(2560, 93, 'Maximum', 'Color', 'k', ... 'FontSize', 13, 'FontName', 'Arial', 'FontWeight', 'bold'); text(2560, 55, 'Average', 'Color', 'k', ... 'FontSize', 13, 'FontName', 'Arial', 'FontWeight', 'bold'); text(2560, 25, 'Minimum', 'Color', 'k', ... 'FontSize', 13, 'FontName', 'Arial', 'FontWeight', 'bold'); % Axis styling set(ax2, 'FontSize', 12, 'FontName','Arial', 'LineWidth', 1.2, ... 'TickDir','out', 'TickLength',[0.01 0.01]); % Horizontal Colorbar: placed between two plots at bottom cb = colorbar(ax1, 'Position', [0.35 0.1 0.3 0.03], ... 'Orientation','horizontal', 'FontName','Arial'); cb.Label.String = 'Residual pressure (MPa)'; cb.Label.FontSize = 13; % Adjust subplot positions to fit colorbar set(ax1, 'Position', [0.08 0.15 0.38 0.75]); set(ax2, 'Position', [0.54 0.15 0.38 0.75]); drawnow; exportgraphics(gcf,'Volcphymodel.png','Resolution',600);