%% Magma Ocean Solidification Example close all clear clc %% IMPORT PARTITIONING DATA, LMO BULK COMPOSITION AND CHONDRITIC NORMALIZATION DataImport=importdata('LMPC.txt'); LMPC=DataImport.data; %% Summary For Section #1 % In the section below we define the partitioning parameters for each % mineral, assumed starting composition of the magma ocean, and element % order for spider diagrams. %% Crystallization Variables % Trace Element Partitioning % Olivine Partitioning DOl = LMPC(1:16,1); %Olivine Coefficients Layer#1 DOl2 = LMPC(1:16,2); %Olivine Coefficients Layer#2 DOl3 = LMPC(1:16,3); %Olivine Coefficients Layer#3 DOl4 = LMPC(1:16,4); %Olivine Coefficients Layer#6 % Orthopyroxene Partitioning DOpx = LMPC(1:16,5); %Orthopyroxene Coefficients Layer#2 DOpx2 = LMPC(1:16,6); %Orthopyroxene Coefficients Layer#3 % Clinopyroxene Partitioning DCpx = LMPC(1:16,7); %Clinopyroxene Coefficients Layer#4 DCpx2 = LMPC(1:16,8); %Clinopyroxene Coefficients Layer#5 DCpx3 = LMPC(1:16,9); %Clinopyroxene Coefficients Layer#6 % Pigeonite Partitioning DPig = LMPC(1:16,10); % Plagioclase Partitioning DPlag = LMPC(1:16,11); %Plagioclase Coefficients Layer#3 DPlag2 = LMPC(1:16,12); %Plagioclase Coefficients Layer#4 DPlag3 = LMPC(1:16,13); %Plagioclase Coefficients Layer#5 DPlag4 = LMPC(1:16,14); %Plagioclase Coefficients Layer#6 % Ilmenite Partitioning DIlm = LMPC(1:16,15); % Quartz Partitioning DQtz = LMPC(1:16,16); % Interstitial Melt Partitioning DIntMlt = LMPC(1:16,17); % Starting Composition of The Magma Ocean LMO = LMPC(1:16,18); % Chondrite Composition CHON = LMPC(1:16,19); % Order of Elements for the Spider Diagrams ElementOrder = LMPC(1:16,20)'; AxisLabel = {'La' 'Ce' 'Pr' 'Nd' 'Sm' 'Eu' 'Gd' 'Tb' 'Dy' 'Y' 'Ho' ... 'Er' 'Tm' 'Yb' 'Lu' 'Hf'}; %% Summary of Section #2 % In the section below we build the melt fraction matrix for the magma % ocean solidification process. For this matrix each magma ocean cumulate % layer represents a 100 point portion of the total 600 point matrix. %% LMO Melt Fractions F1out = 0.49; % Melt fraction at which cumulate layer#1 ends F2out = 0.22; % Melt fraction at which cumulate layer#2 ends F3out = 0.19; % Melt fraction at which cumulate layer#3 ends F4out = 0.05; % Melt fraction at which cumulate layer#4 ends F5out = 0.02; % Melt fraction at which cumulate layer#5 ends % layer#6 ends at F=0 % Building a melt fraction Matrix Fs1 = linspace(1,F1out,100); % 100 point matix for cumulate layer#1 Fs2 = linspace(F1out,F2out,100); % 100 point matix for cumulate layer#2 Fs3 = linspace(F2out,F3out,100); % 100 point matix for cumulate layer#3 Fs4 = linspace(F3out,F4out,100); % 100 point matix for cumulate layer#4 Fs5 = linspace(F4out,F5out,100); % 100 point matix for cumulate layer#5 Fs6 = linspace(F5out,0,100); % 100 point matix for cumulate layer#6 F = [Fs1 Fs2 Fs3 Fs4 Fs5 Fs6]; % 600 point matix for LMO solidification %% Summary of Section #3 % In the section below we determine the modal mineralogy of the magma ocean % cumulates. For calculations of the residual liquid composition we % determine the modal mineralogy prior to floatation crust formation. For % calculations of the mantle cumulate composition we determine the modal % mineralogy of mantle cumulates after floatation. For calculations of % crustal cumulates compositions we determine the modal mineralogy of % crustal cumulates after floatation. For this example calculation we % assume 80% plagioclase floatation and 100% quartz flaotation. %% Modal Mineralogy of LMO Cumulates % The modal mineral proportions below will be used to calculate the % bulk partition coefficient of the magma ocean cumulates as a whole. % This bulk partition coefficient will be used to calculate the evolution % of the residual magma ocaean liquids. This means we leave the total % amount of crystallizing plagioclase and quartz in the total modal % mineralogy. The proportion of each mineral in each cumulate layer are % presented in the order below. % Mineralmode = [Layer#1 Layer#2 Layer#3 Layer#4 Layer#5 Layer#6] OlMode1 = [0.990 0.099 0.050 0.000 0.000 0.198]; %fraction of Olivine in LMO cumulates OpxMode1 = [0.000 0.891 0.198 0.000 0.000 0.000]; %fraction of Orthopyroxene in LMO cumulates CpxMode1 = [0.000 0.000 0.000 0.198 0.198 0.158]; %fraction of Clinopyroxene in LMO cumulates PigMode1 = [0.000 0.000 0.248 0.317 0.248 0.000]; %fraction of Pigeonite in LMO cumulates PlagMode1 = [0.000 0.000 0.495 0.475 0.446 0.416]; %fraction of Plagioclase in LMO cumulates IlmMode1 = [0.000 0.000 0.000 0.000 0.000 0.119]; %fraction of Ilmenite in LMO cumulates QtzMode1 = [0.000 0.000 0.000 0.000 0.099 0.099]; %fraction of Quartz in LMO cumulates IntMltMode1 = [0.010 0.010 0.010 0.010 0.010 0.010]; %fraction of Trapped Melt in LMO cumulates % The modal mineral proportions below will be used to calculate the % bulk partition coefficient of the magma ocean cumulates at each % calculation interval. This bulk partition coefficient will be used to % calculate the evolution of magma ocean mantle cumulates. This means we % remove a porportion of crystallizing plagioclase and quartz in the total % modal mineralogy since this plagioclase and quartz is removed to form % the crust. For this model we assume 80% plagioclase floatation and 100% % quartz floatation. The proportion of each mineral in each cumulate layer % are presented in the order below. % Mineralmode = [Layer#1 Layer#2 Layer#3 Layer#4 Layer#5 Layer#6] OlMode2 = [0.990 0.099 0.083 0.000 0.000 0.352]; %fraction of Olivine in mantle OpxMode2 = [0.000 0.891 0.331 0.000 0.000 0.000]; %fraction of Orthopyroxene in mantle CpxMode2 = [0.000 0.000 0.000 0.322 0.367 0.281]; %fraction of Clinopyroxene in mantle PigMode2 = [0.000 0.000 0.413 0.515 0.459 0.000]; %fraction of Pigeonite in mantle PlagMode2 = [0.000 0.000 0.164 0.153 0.164 0.146]; %fraction of Plagioclase in mantle IlmMode2 = [0.000 0.000 0.000 0.000 0.000 0.211]; %fraction of Ilmenite in mantle QtzMode2 = [0.000 0.000 0.000 0.000 0.000 0.000]; %fraction of Quartz in mantle IntMltMode2 = [0.010 0.010 0.010 0.010 0.010 0.010]; %fraction of Trapped Melt in mantle % The modal mineral proportions below will be used to calculate the % bulk partition coefficient of the crustal cumulates at calculation % intervals 201 to 600. This modal mineral proportion represents the total % amount of plagioclase and quartz that we previously removed from the % mantle cumulates plus 0.5% trapped melt. The proportion of each mineral % or phase in each cumulate layer are presented in the order below. % Phasemode = [Layer#3 Layer#4 Layer#5 Layer#6] PlagMode3 = [0.995 0.995 0.779 0.767]; %fraction of Plagioclase in crust QtzMode3 = [0.000 0.000 0.216 0.228]; %fraction of Quartz in crust IntMltMode3 = [0.005 0.005 0.005 0.005]; %fraction of trapped melt in crust %% Summary of Section #4 % The section below calculates the bulk partition coefficients fro the % total magma ocean cumulate pile, the mantle cumulates, and the crustal % cumulates. %% Bulk Partition Coefficient Calculations % Bulk Partition Coefficient for each step Ds11 = OlMode1(1,1).*DOl+OpxMode1(1,1).*DOpx+... CpxMode1(1,1).*DCpx+PigMode1(1,1).*DPig+PlagMode1(1,1).*DPlag+... IlmMode1(1,1).*DIlm+IntMltMode1(1,1).*DIntMlt+QtzMode1(1,1).*DQtz; %Layer #1 bulk partition coefficient without floatation Ds21 = OlMode1(1,2).*DOl2+OpxMode1(1,2).*DOpx+... CpxMode1(1,2).*DCpx+PigMode1(1,2).*DPig+PlagMode1(1,2).*DPlag+... IlmMode1(1,2).*DIlm+IntMltMode1(1,2).*DIntMlt+QtzMode1(1,2).*DQtz; %Layer #2 bulk partition coefficient without floatation Ds31 = OlMode1(1,3).*DOl3+OpxMode1(1,3).*DOpx2+... CpxMode1(1,3).*DCpx+PigMode1(1,3).*DPig+PlagMode1(1,3).*DPlag+... IlmMode1(1,3).*DIlm+IntMltMode1(1,3).*DIntMlt+QtzMode1(1,3).*DQtz; %Layer #3 bulk partition coefficient without floatation Ds41 = OlMode1(1,4).*DOl+OpxMode1(1,4).*DOpx+... CpxMode1(1,4).*DCpx+PigMode1(1,4).*DPig+PlagMode1(1,4).*DPlag2+... IlmMode1(1,4).*DIlm+IntMltMode1(1,4).*DIntMlt+QtzMode1(1,4).*DQtz; %Layer #4 bulk partition coefficient without floatation Ds51 = OlMode1(1,5).*DOl+OpxMode1(1,5).*DOpx+... CpxMode1(1,5).*DCpx2+PigMode1(1,5).*DPig+PlagMode1(1,5).*DPlag3+... IlmMode1(1,5).*DIlm+IntMltMode1(1,5).*DIntMlt+QtzMode1(1,5).*DQtz; %Layer #5 bulk partition coefficient without floatation Ds61 = OlMode1(1,6).*DOl4+OpxMode1(1,6).*DOpx+... CpxMode1(1,6).*DCpx3+PigMode1(1,6).*DPig+PlagMode1(1,6).*DPlag4+... IlmMode1(1,6).*DIlm+IntMltMode1(1,6).*DIntMlt+QtzMode1(1,6).*DQtz; %Layer #6 bulk partition coefficient without floatation % Total LMO Cumulate Bulk Partition Coefficient for i = 1:100 DBulkFractional1(1:16,i) = Ds11; end for i = 101:200 DBulkFractional1(1:16,i) = (Ds11.*(1-F1out)+Ds21.*(F1out-... F(i)))./(1-F(i)); end for i = 201:300 DBulkFractional1(1:16,i) = (Ds11.*(1-F1out)+Ds21.*(F1out-... F2out)+Ds31.*(F2out-F(i)))./(1-F(i)); end for i = 301:400 DBulkFractional1(1:16,i) = (Ds11.*(1-F1out)+Ds21.*(F1out-... F2out)+Ds31.*(F2out-F3out)... +Ds41.*(F3out-F(i)))./(1-F(i)); end for i = 401:500 DBulkFractional1(1:16,i) = (Ds11.*(1-F1out)+Ds21.*(F1out-... F2out)+Ds31.*(F2out-F3out)... +Ds41.*(F3out-F4out)+... Ds51.*(F4out-F(i)))./(1-F(i)); end for i = 501:600 DBulkFractional1(1:16,i) = (Ds11.*(1-F1out)+Ds21.*(F1out-... F2out)+Ds31.*(F2out-F3out)... +Ds41.*(F3out-F4out)+Ds51.*(F4out-F5out)+... Ds61.*(F5out-F(i)))./(1-F(i)); end % Bulk Partition Coefficient for each step, Mantle Cumulates Ds12 = OlMode2(1,1).*DOl+OpxMode2(1,1).*DOpx+... CpxMode2(1,1).*DCpx+PigMode2(1,1).*DPig+PlagMode2(1,1).*DPlag+... IlmMode2(1,1).*DIlm+IntMltMode2(1,1).*DIntMlt+QtzMode2(1,1).*DQtz; %Layer #1 bulk partition coefficient after floatation Ds22 = OlMode2(1,2).*DOl2+OpxMode2(1,2).*DOpx+... CpxMode2(1,2).*DCpx+PigMode2(1,2).*DPig+PlagMode2(1,2).*DPlag+... IlmMode2(1,2).*DIlm+IntMltMode2(1,2).*DIntMlt+QtzMode2(1,2).*DQtz; %Layer #2 bulk partition coefficient after floatation Ds32 = OlMode2(1,3).*DOl3+OpxMode2(1,3).*DOpx2+... CpxMode2(1,3).*DCpx+PigMode2(1,3).*DPig+PlagMode2(1,3).*DPlag+... IlmMode2(1,3).*DIlm+IntMltMode2(1,3).*DIntMlt+QtzMode2(1,3).*DQtz; %Layer #3 bulk partition coefficient after floatation Ds42 = OlMode2(1,4).*DOl+OpxMode2(1,4).*DOpx+... CpxMode2(1,4).*DCpx+PigMode2(1,4).*DPig+PlagMode2(1,4).*DPlag2+... IlmMode2(1,4).*DIlm+IntMltMode2(1,4).*DIntMlt+QtzMode2(1,4).*DQtz; %Layer #4 bulk partition coefficient after floatation Ds52 = OlMode2(1,5).*DOl+OpxMode2(1,5).*DOpx+... CpxMode2(1,5).*DCpx2+PigMode2(1,5).*DPig+PlagMode2(1,5).*DPlag3+... IlmMode2(1,5).*DIlm+IntMltMode2(1,5).*DIntMlt+QtzMode2(1,5).*DQtz; %Layer #5 bulk partition coefficient after floatation Ds62 = OlMode2(1,6).*DOl4+OpxMode2(1,6).*DOpx+... CpxMode2(1,6).*DCpx3+PigMode2(1,6).*DPig+PlagMode2(1,6).*DPlag4+... IlmMode2(1,6).*DIlm+IntMltMode2(1,6).*DIntMlt+QtzMode2(1,6).*DQtz; %Layer #6 bulk partition coefficient after floatation % Mantle Bulk Partition Coefficient for i = 1:100 DBulkFractional2(1:16,i) = Ds12; end for i = 101:200 DBulkFractional2(1:16,i) = Ds22; end for i = 201:300 DBulkFractional2(1:16,i) = Ds32; end for i = 301:400 DBulkFractional2(1:16,i) = Ds42; end for i = 401:500 DBulkFractional2(1:16,i) = Ds52; end for i = 501:600 DBulkFractional2(1:16,i) = Ds62; end % Crustal Bulk Partition Coefficient for each interval for i=1:100 DAnorth(1:16,i) = DPlag.*PlagMode3(1,1)+DIntMlt.*IntMltMode3(1,1)... +DQtz.*QtzMode3(1,1); end for i=101:200 DAnorth(1:16,i) = DPlag2.*PlagMode3(1,2)+DIntMlt.*IntMltMode3(1,2)... +DQtz.*QtzMode3(1,2); end for i=201:300 DAnorth(1:16,i) = DPlag3.*PlagMode3(1,3)+DIntMlt.*IntMltMode3(1,3)... +DQtz.*QtzMode3(1,3); end for i=301:400 DAnorth(1:16,i) = DPlag4.*PlagMode3(1,4)+DIntMlt.*IntMltMode3(1,4)... +DQtz.*QtzMode3(1,4); end %% Summary of Section #5 % In the following section we model the evolution of cumulate and melt % compositions as the magma ocean solidifies. %% Chemical Evolution Calculations ClFractional = LMO.*F.^(DBulkFractional1-1); % Residual Liquid CsFractional = (ClFractional.*DBulkFractional2); % Mantle Cumulates CsCrust = ClFractional(1:16,201:600).*DAnorth; % Crustal Cumulates %% Summary of Section #6 % The section below selects data from specific intervals of solidification % so the data can be plotted on spider diagrams. %% Matrix for the spider diagrams ClFractional0 = ClFractional(1:16,2); %Liquid @ 0.5% solidification ClFractional1 = ClFractional(1:16,98); %Liquid @ 50% solidification ClFractional2 = ClFractional(1:16,189); %Liquid @ 75% solidification ClFractional3 = ClFractional(1:16,364); %Liquid @ 90% solidification ClFractional4 = ClFractional(1:16,500); %Liquid @ 98% solidification ClFractional5 = ClFractional(1:16,595); %Liquid @ 99.9% solidification ClFractionalC = [ClFractional0 ClFractional1 ClFractional2 ClFractional3... ClFractional4 ClFractional5]; % Matrix for residual liquid compositions CsFractional0 = CsFractional(1:16,2); %Mantle @ 0.5% solidification CsFractional1 = CsFractional(1:16,98); %Mantle @ 50% solidification CsFractional2 = CsFractional(1:16,189); %Mantle @ 75% solidification CsFractional3 = CsFractional(1:16,364); %Mantle @ 90% solidification CsFractional4 = CsFractional(1:16,500); %Mantle @ 98% solidification CsFractional5 = CsFractional(1:16,595); %Mantle @ 99.9% solidification CsFractionalC = [CsFractional0 CsFractional1 CsFractional2 CsFractional3... CsFractional4 CsFractional5]; % Matrix for mantle cumulate compositions CsCrust0 = CsCrust(1:16,67); %Crust @ 80% solidification CsCrust1 = CsCrust(1:16,164); %Crust @ 90% solidification CsCrust2 = CsCrust(1:16,300); %Crust @ 98% solidification CsCrust3 = CsCrust(1:16,395); %Crust @ 99.9% solidification CsCrustC = [CsCrust0 CsCrust1 CsCrust2 CsCrust3]; % Matrix for Crustal cumulate compositions %% Plotted Data % Liquid Compositions figure plot(ElementOrder,(ClFractionalC(1:16,:)./CHON(1:16,:))... ,'LineWidth',2,'LineStyle','-') set(gca,'yscale','log') set(gca,'xlim',[0 17]) set(gca,'ylim',[1E-2 1E4]) set(gca,'xtick',ElementOrder) set(gca,'xticklabel',AxisLabel,'FontSize',12) xlabel('REEs') set(gca,'TickLength',[0.02 0.02]) ylabel('Log Concentration Cnorm') title('Residual Magma Ocean Liquids') legend('0.5 PCS','50 PCS','75 PCS','90 PCS','98 PCS','99.9 PCS',... 'Location','best') % Mantle Compositions figure plot(ElementOrder,(CsFractionalC(1:16,:)./CHON(1:16,:))... ,'LineWidth',2,'LineStyle','-') set(gca,'yscale','log') set(gca,'xlim',[0 17]) set(gca,'ylim',[1E-2 1E4]) set(gca,'xtick',ElementOrder) set(gca,'xticklabel',AxisLabel,'FontSize',12) xlabel('REEs') set(gca,'TickLength',[0.02 0.02]) ylabel('Log Concentration Cnorm') title('Mantle Cumulates') legend('0.5 PCS','50 PCS','75 PCS','90 PCS','98 PCS','99.9 PCS',... 'Location','best') % Crustal Compositions figure plot(ElementOrder,(CsCrustC(1:16,:)./CHON(1:16,:))... ,'LineWidth',2,'LineStyle','-') set(gca,'yscale','log') set(gca,'xlim',[0 17]) set(gca,'ylim',[1E-2 1E4]) set(gca,'xtick',ElementOrder) set(gca,'xticklabel',AxisLabel,'FontSize',12) xlabel('REEs') set(gca,'TickLength',[0.02 0.02]) ylabel('Log Concentration Cnorm') title('Crustal Cumulates') legend('80 PCS','90 PCS','98 PCS','99.9 PCS',... 'Location','best')