diff --git a/grassland/elem_ratios.csv b/grassland/elem_ratios.csv new file mode 100644 index 0000000..deb1e69 --- /dev/null +++ b/grassland/elem_ratios.csv @@ -0,0 +1,13 @@ +,C,N,P +DeadMic,0.0,0.0,0.0 +DeadEnz,0.0,0.0,0.0 +Cellulose,1.00,0.0,0.0 +Hemicellulose,1.00,0.0,0.0 +Starch,1.00,0.0,0.0 +Chitin,6.01,1.00,0.0 +Lignin,118.58,1.00,0.0 +Protein1,5.05,1.00,0.0 +Protein2,5.05,1.00,0.0 +Protein3,5.05,1.00,0.0 +OrgP1,26.08,0.0,1.00 +OrgP2,3.8,1.67,1.00 diff --git a/src/grid.py b/src/grid.py index a1d1367..2511ca4 100644 --- a/src/grid.py +++ b/src/grid.py @@ -67,13 +67,18 @@ def __init__(self, runtime, data_init): ] # Rate for change of moisture reaction in decay module # Uptake - # self.Microbes_init = data_init['Microbes_pp'] # microbial community before placement + self.Microbes_init = data_init[ + "Microbes_pp" + ] # microbial community before placement self.Microbes = data_init["Microbes"].copy( deep=True ) # microbial community after placement - # self.Monomers_init = data_init['Monomers'] # Monomers initialized + self.Monomers_init = data_init["Monomers"] # Monomers initialized self.Monomers = data_init["Monomers"].copy(deep=True) # Monomers self.MonInput = data_init["MonInput"] # Inputs of monomers + self.C_frac_org = np.tile(data_init["elem_fracs"]["C"].values, self.gridsize) + self.N_frac_org = np.tile(data_init["elem_fracs"]["N"].values, self.gridsize) + self.P_frac_org = np.tile(data_init["elem_fracs"]["P"].values, self.gridsize) self.Uptake_Ea = data_init["Uptake_Ea"] # transporter enzyme Ea self.Uptake_Vmax0 = data_init["Uptake_Vmax0"] # transporter Vmax self.Uptake_Km0 = data_init["Uptake_Km0"] # transporter Km @@ -235,7 +240,11 @@ def degradation(self, day): ) # Update Substrates Pool by removing decayed C, N, & P. Depending on specific needs, adding inputs of substrates can be done here - self.Substrates -= SubstrateRatios.mul(DecayRates, axis=0) # + self.SubInput + self.Substrates -= SubstrateRatios.mul( + DecayRates, axis=0 + ) # First remove decayed matter + self.Substrates += self.SubInput # then add daily substrates + self.Substrates[self.Substrates < 0] = np.float32(0) # Pass these two back to the global variables to be used in the next method self.SubstrateRatios = SubstrateRatios @@ -264,6 +273,8 @@ def uptake(self, day): self.Monomers.index != "PO4" ) # organic monomers # is_mineral = (Monomers.index == "NH4") | (Monomers.index == "PO4") + is_NH4 = self.Monomers.index == "NH4" + is_PO4 = self.Monomers.index == "PO4" # Update monomer ratios in each time step with organic monomers following the substrates self.Monomer_ratios[is_org] = self.SubstrateRatios.values @@ -271,11 +282,32 @@ def uptake(self, day): # Determine monomer pool from decay and input # Organic monomers derived from substrate-decomposition Decay_Org = self.Monomer_ratios[is_org].mul(self.DecayRates.values, axis=0) + # Monomer pool determined + self.Monomers.loc[is_org] += Decay_Org # + Input_Org + + # --- Route monomer inputs by pool --- + input_vals = self.MonInput.values # positionally aligned with self.Monomers + + # Organic monomers: input is carbon, goes to the C column + org_input = input_vals[is_org] + self.Monomers.loc[is_org, "C"] += org_input * self.C_frac_org + self.Monomers.loc[is_org, "N"] += org_input * self.N_frac_org + self.Monomers.loc[is_org, "P"] += org_input * self.P_frac_org + + # Mineral N: NH4 input goes to the N column + self.Monomers.loc[is_NH4, "N"] += input_vals[is_NH4] + + # Mineral P: PO4 input goes to the P column + self.Monomers.loc[is_PO4, "P"] += input_vals[is_PO4] + + self.Monomers = self.Monomers.fillna(0.0) + self.Monomers[self.Monomers < 0] = np.float32(0) + # inputs of organic and mineral monomers # Input_Org = MR_transition[is_org].mul(self.MonInput[is_org].tolist(),axis=0) # Input_Mineral = MR_transition[is_mineral].mul((self.MonInput[is_mineral]).tolist(),axis=0) # Monomer pool determined - self.Monomers.loc[is_org] += Decay_Org # + Input_Org + # self.Monomers.loc[is_org] += Decay_Org # + Input_Org # self.Monomers.loc[is_mineral] += Input_Mineral # Get the total mass of each monomer: C+N+P @@ -330,6 +362,7 @@ def uptake(self, day): # Update Monomers # By monomer: total uptake (monomer*gridsize) * 3(C-N-P) self.Monomers -= self.Monomer_ratios.mul(Uptake.sum(axis=1), axis=0) + self.Monomers[self.Monomers < 0] = np.float32(0) # Derive Taxon-specific total uptake of C, N, & P # By taxon: total uptake; (monomer*gridsize) * taxon @@ -744,6 +777,7 @@ def mortality(self, day): # Update monomer pools self.Monomers.loc[is_NH4, "N"] += sum(MicLoss["N"]) / self.gridsize self.Monomers.loc[is_PO4, "P"] += sum(MicLoss["P"]) / self.gridsize + self.Monomers[self.Monomers < 0] = np.float32(0) # Update Substrates pool by adding dead microbial biomass self.Substrates.loc[is_DeadMic] += Death_gridcell.values @@ -911,6 +945,15 @@ def reinitialization(self, initialization, microbes_pp, output, mode, pulse, *ar self.Substrates = initialization["Substrates"].copy(deep=True) self.Monomers = initialization["Monomers"].copy(deep=True) self.Enzymes = initialization["Enzymes"].copy(deep=True) + self.SubstrateRatios = pd.DataFrame( + np.zeros_like(initialization["Substrates"]), + index=initialization["Substrates"].index, + columns=initialization["Substrates"].columns, + ) + self.DecayRates = pd.Series( + np.zeros(len(initialization["Substrates"])), + index=initialization["Substrates"].index, + ) # reinitialize microbial community in a new pulse as per the mode in three steps # first: retrieve the microbial pool; NOTE: copy() @@ -940,23 +983,12 @@ def reinitialization(self, initialization, microbes_pp, output, mode, pulse, *ar # calculate frequency of every taxon frequencies = cum_abundance / cum_abundance.sum() frequencies = frequencies.fillna(0) - frequencies = frequencies.clip( - lower=0, upper=1 - ) # guard against any stray negative/over-1 values (sum of 1) - frequencies = ( - frequencies / frequencies.sum() if frequencies.sum() > 0 else frequencies - ) # renormalize to exactly 1 # last: assign microbes to each grid box randomly based on prior densities choose_taxa = np.zeros((self.n_taxa, self.gridsize), dtype="int8") for i in range(self.n_taxa): - freq = np.float64( - frequencies.iloc[1] - ) # working independent of indexing (based on location not index labels) - p_vec = np.array([freq, 1.0 - freq], dtype=np.float64) # solve numoy issues - p_vec = p_vec / p_vec.sum() - choose_taxa[i, :] = np.random.choice( - [1, 0], self.gridsize, replace=True, p=p_vec + choose_taxa[i, :] = np.random.binomial( + 1, frequencies.iloc[i], self.gridsize ) self.Microbes.loc[np.ravel(choose_taxa, order="F") == 0] = np.float32( 0 diff --git a/src/initialization.py b/src/initialization.py index 6f11437..cfb4847 100644 --- a/src/initialization.py +++ b/src/initialization.py @@ -45,12 +45,19 @@ def initialize_data(runtime_parameters, site): Ea_input = pd.read_csv(site + "/" + "enzyme_ea.csv", header=0, index_col=0).astype( "float32" ) # enzyme activation energy + elem_ratios = pd.read_csv( + site + "/" + "elem_ratios.csv", header=0, index_col=0 + ).astype("float32") # elemental (C:N:P) ratios of organic monomer inputs # climate forcings climate = pd.read_csv(site + "/" + "climate.csv", header=0, index_col=0) # daily temperature and water potential daily_temp = climate["Temp"].to_numpy(dtype="float32") # temperaure series daily_psi = climate["Psi"].to_numpy(dtype="float32") # water potential series + # Convert C:N:P ratios into fractions of total (each row sums to 1.0) + row_totals = elem_ratios.sum(axis=1) + elem_fracs = elem_ratios.div(row_totals, axis=0).fillna(0.0) + # ...an instance of Substrate class Substrates = Substrate(runtime_parameters, parameters, substrates_init) # ...substrate initial pool size @@ -138,6 +145,7 @@ def initialize_data(runtime_parameters, site): "Monomers": expand(monomers_initial_pool, gridsize), "Monomer_ratio": expand(monomer_ratio_inital, gridsize), "MonInput": expand(monomers_input_rate, gridsize), + "elem_fracs": elem_fracs, "Uptake_ReqEnz": expand(monomers_uptake_reqenzyme, gridsize), "Enzymes": expand(enzymes_initial_pool, gridsize), "Km0": expand(enzymes_Km, gridsize), # enzyme half-saturation constant diff --git a/src/substrate.py b/src/substrate.py index 30b663e..7ac4dc5 100644 --- a/src/substrate.py +++ b/src/substrate.py @@ -56,7 +56,7 @@ def substrate_input(self, sub_mon_input): SubInputN.name = "N" SubInputP.name = "P" SubInput_df = pd.concat([SubInputC, SubInputN, SubInputP], axis=1, sort=False) - SubInput_df["DeadMic"] = SubInput_df["DeadEnz"] = 0 # Change NAs to 0 + SubInput_df = SubInput_df.fillna(0) SubInput_df = SubInput_df.astype("float32") return SubInput_df