• Run the code in MATLAB

• Several errors come up which were expected, the first thing we've encountered is in this line of code

s_H(:,1) = (z_H^2)*v_H*G_H(:,1);
s_OH(:,1) = (z_OH^2)*v_OH*G_OH(:,1);
s_C(:,1) = (z_C^2)*v_C*G_C(:,1);
Sigma(:,1) = (F^2)*(s_H + s_OH + s_C);

• As we are only changing values for all x nodes in initial time node we need to add (:,1) to fix the error, ideally we will have the following code

s_H(:,1) = (z_H^2)*v_H*G_H(:,1);
s_OH(:,1) = (z_OH^2)*v_OH*G_OH(:,1);
s_C(:,1) = (z_C^2)*v_C*G_C(:,1);
Sigma(:,1) = (F^2)*(s_H(:,1) + s_OH(:,1) + s_C(:,1));

• In this line we are using alpha which is not defined in code and it should be 1alpha_Cinstead. not to mention all concentration should be for Carbon orG_Cin initial code I have usedG` which is not mentioned in this code.

G_C(i,m+1) = G_C(i,m) + sub(i,m)*(alpha*(G_C(i+1,m) -2*G_C(i,m) + G_C(i-1,m)) + beta_C*(G(i+1,m) - G(i-1,m)) + R_C(i,m)/R_D);

• So new code will look like this:

G_C(i,m+1) = G_C(i,m) + sub(i,m)*(alpha_C*(G_C(i+1,m) -2*G_C(i,m) + G_C(i-1,m)) + beta_C*(G_C(i+1,m) - G_C(i-1,m)) + R_C(i,m)/R_D);

• In line 250 we face this code:

G_HA(i,m+1) = G_HA(i,m) + sub(i,m)*(alpha_HA*(G_HA(i+1,m) -2*G_HA(i,m) + G_HA(i-1,m)) + beta_HA*(G_HA(i+1,m) - G_HA(i-1,m)) + R_HA(i,m)/R_D);

• This is mainly because we have not created any of active models arays, to do so these lines of codes were added

G_HA = zeros(nx,nt);
G_BOH = zeros(nx,nt);
G_A = zeros(nx,nt);
G_B = zeros(nx,nt);

• Following code returned an unknown parameter error

G_HA(i,m+1) = G_HA(i,m) + sub(i,m)*(alpha_HA*(G_HA(i+1,m) -2*G_HA(i,m) + G_HA(i-1,m)) + beta_HA*(G_HA(i+1,m) - G_HA(i-1,m)) + R_HA(i,m)/R_D);

• This is because reaction rates were not included for this part of model, after adding needed parameter problem fixed for both acid and base

• Since rate for both acid and base are only calculated separately at both x=0 and x=L there are two conditional if in code, though a mistake happened and made both R_prime for acid and base as array which is wrong. As result following line of code resulted in returning syntax error

R_OH(i,m) = -1*R_prime_OH(i,m);

and

R_H(i,m) = -1*R_prime_H(i,m);

• In order to fix the issue we need to remove (i,m) from prime values.

Pic1

How ever this is not the case in modified code and no plots appear, it is assumed to have an issue with constants but to make sure everything is fine mat data of original code and modified one were exported to make sure everything aligns within accepted margin.

So basically, there are two mat files, one is test.mat and the other one is final.mat which final is all about modified code as it is going to be our final code usable in this repo

• 1. Valency changes:

In this term all values related to this changed naming upon to what species we are dealing with, basically z1 z2 z3 changed to z_H1 z_OH z_C and so on. Does it seem to affect plots? No ❌

• 2. Physical Info

These values are mostly keeping their naming mechanism however some values were removed as instance there is no value after electric field in old code but there is one with phi except these changes epsilon for each specie has been removed, this needs some confirmation but epsilon seems to only applied to solution not each species in domain. Does it seem to affect plots? No ❌

• 3. Diffusion and primary advection

These values are essentially as same as old code so it isn't expected to consider these values essential for this error, although the most effective value in advection for plots are these values, unlike the old code in this one all parameters are in one block and are sorted based on what value and species we are dealing, that being said some names for these values might have change/removed or added. Even with these changes values within the parameters are same.

Does it seem to affect plots? No ❌ (Only for plot issue we have)

• 4. Velocity field

Same as diffusion, they might have several differences when it comes to value and names but are not considered too much of a change however a simple change might change the velocity compared to old code which is reflected in total velocity value, in this new code a new form of velocity has been added an it is convection velocity, since it is a big parameter, it leads into huge offsets which is not expected to see in this code

Does it seem to affect plots? Yes ✔

In order to check and validate differences in both working codes, mat files of both were created, with the help of compare tool within MATLAB we can locate what has been added or what has been removed. In this specific case most of variables changed their names but not values, so basically almost all the variables are same as before however some parameters have changed in a meaningful way so it is worth mentioning what are these variables.

Generally speaking, this code contains EKR and BKR however finished code I am sharing here does only have Bio-EKR, although it can be instituted with EKR only with changing diffusion and velocity advections with proper terms already calculated in the code