An electric potential gradient in a soil lets the cations in the pore water move to the cathode and the anions to the anode, counterbalancing the negative and positive charges on the electrodes. It can be used to obtain the diffusion rate of Cl^{} in concrete that is otherwise too slow to be measurable in laboratory experiments. And it can assist in the remediation of contaminated soils together with the oxidation/reduction reactions that accompany the application of an electric field. Details of the calculations can be found in Solute transport solved with the NernstPlanck equation for concrete pores with 'free' water and a double layer, Cem. Concr. Res. 101, 102113 (Ask Tony Appelo for a reprint)
The flux due to the gradients in concentration and electric potential is given by the NernstPlanck equation:
(1) 
Figure 1. The (negative) current calculated by PHREEQC over time, given by the special basic term 'current_A', from the H_{2} generated, and from the concentration changes in the cathode cell (PHREEQC input file current1.phr).
The current 'from H_{2} in cell 0' (green line, overlaps with current_A) is calculated from the moles of H_{2} that are captured in EQUILIBRIUM_PHASES 0; H2(g) 10 0 when the concentration is higher than 10^{13.2} M (SI = 10). The current 'from Na+ transport' is calculated from the concentration increase in the cathode cell (blue line). The orange line, marked with '+ Cl' gives the sum of Na^{+} and Cl^{} transport. After 10 hours, the sum is smaller than the total current (red and green lines), because OH^{}, that is generated by the reduction reaction contributes significantly to the current (see figure 3 below). The magenta crosses give the sum of Na^{+}, Cl^{} and OH^{} transport, which is the same as the actual current.  
At zero time, the current agrees with what was calculated before. But as time progresses, the ions move into the boundary reservoirs, and the smaller concentrations in the pore water of the sample can only conduct a smaller current with a constant potential difference over the column.
To find the details of the current decrease over time, we must look at the NernstPlanck equation. The electric flux of the ions is directly proportional to their concentration, and the concentration will be higher near the attracting electrode than at the other end where the ion is repelled. It means that the sample is depleted of the ion, or, in other words, that the concentration decreases in the column. Consequently, the current decreases, and it must do so equally, everywhere in a circuit that has no capacitor. But then, the potential gradient must be higher where the concentration is smaller. And this is the case indeed, as shown in Figure 2.
Figure 2. (Left) The progress of the Na^{+} concentration with time in the sample and the boundary cells where the potential is applied. The Cl^{} concentration is a mirror of the Na^{+} concentration and is shown when the PHREEQC file is run. (Right) The potential over the sample at the same timesteps. Initially, the gradient is linear when the concentrations are the same everywhere. But when the concentrations change, the gradient changes as well, and it becomes highest where the concentrations are smallest, maintaining so at all points in the column the same, but smaller than initial current.
The current due to Na^{+} and Cl^{}, divided by the total current, give the transport numbers of the ions (the fraction of the current transported by an ion). Clearly, it changes when the pH changes markedly in the cell. If one of the solute ions takes part in the electrode reactions, for example, Cl^{} in an Ag/AgCl electrode, the pH change is insignificant, since the reactions AgCl = Ag^{+} + Cl^{} and Ag^{+} + e^{} = Ag are not depending on H^{+} and OH^{}.
The charge balancing by OH^{} and H^{+} is evident in the pH change shown in Figure 3. The concentrations of the two ions are significant, but at the opposite ends of the column, and their transport must be included in the calculations, while the activities of OH^{} and H^{+} are connected by the law of mass action:
[OH^{}][H^{+}] = 10^{14} (at 25°C) (2).
Figure 3. The pH for a number of timesteps as a function of distance in the column. At high pH, OH^{} balances the surplus of positive charge at the cathode, at low pH, H^{+} balances the negative charge at the anode. The hydroxyls and protons are generated by reduction and oxidation of water:

We started with equal diffusion coefficient of 1e9 m^{2}/s for all solute species. But if we set D(Cl^{}) = 2e9 m^{2}/s, the initial current is 17.7 * (1 + 2) / 2 = 26.55 mA. You can check this, uncommenting the definition of the Clsolution species in the file. And, if phreeqc.dat is used, which contains definitions of individual diffusion coefficients of the species, the initial current will be 29.7 mA.
When phreeqc.dat is used, all the solutes have their own diffusion coefficient, and the protons, generated at the anode, diffuse 7 times faster than Na^{+}, and about twice faster than the hydroxyls. As a result, the pHjump, that is in the middle of the column in figure 3, will move towards the cathode, and the concentration patterns of the ions become more asymmetric, as shown in figure 4. The increasing gradient of the pH enhances diffusional transport of OH^{} from the cathode reservoir, and this flow becomes an even more important part of the current.
Thus, when the file is run with phreeqc.dat, the current initially increases with time, in contrast with the results shown in Figure 1. And the current, calculated from the transfer of Na^{+} and Cl^{} in the cathode reservoir, is now much smaller than the actual current. However, when the Na^{+} and Cl^{}concentrations become nearzero, the current drops like in Figure 1.
Figure 4. (Left) Timechanges of the Na^{+} and Cl^{}concentrations vs. distance when phreeqc.dat is used. Note that transport of Na^{+} is slower than of Cl^{} (the concentration changes of Na^{+} in the reservoirs are smaller), because the diffusion coefficient is smaller (compare with Figure 2). (Right) The current as a function of time. The model is run with a smaller grid for numerical stability (100 column cells of 0.3 mm length).
For checking the code, Figure 5 compares the steady state concentrations and potential with an analytical solution (Krabbenhøft and Krabbenhøft, 2008, Cement and Concrete Research 38, 7788), using input file current2.phr. The OH^{}concentration in the two boundary solutions is 108 mmol/kgw, and H^{+} is negligible according to eqn (2). In the model, there are 10 cells of 0.1 mm length at the ends, and 5 cells of 1.6 mm in the middle, thus, 25 cells for the 10 mm sample. It is of interest to note that in most of the column, the ratios of the concentrations are determined by the ratios in the source reservoir, thus, the anode reservoir for the cations (K^{+} / Na^{+} = 83 / 25 = 3.32), and the cathode reservoir for the anions (Cl^{} / OH^{} = 500 / 108). For example, at 5 mm, K^{+} / Na^{+} = 275 / 82.9 (= 3.32), and Cl^{} / OH^{} = 295 / 63.6.
The analytical solution is independent of the diffusion coefficients, and it does not matter whether wateq4f.dat is used as database, or phreeqc.dat (with individual diffusion coefficients for the solute species). However, with phreeqc.dat the runtime is longer because the diffusion coefficient of H^{+} is higher; the runtime can be reduced, redefining dw 1e9 for H^{+} since the proton concentration is everywhere negligible in this example.
Figure 5. (Left) Steady state concentrations of Na^{+}, K^{+}, Cl^{} and OH^{}concentrations with distance. (Right) The steady state potential over the sample. Lines from an analytical solution, + symbols from PHREEQC.
It is well known that the diffusion coefficients in an electric field vary with solution composition, but in the calculations for Figure 5 they are constant since the analytical solution cannot account for such variations. In PHREEQC, an option is available to let diffusion coefficients change similar to the calculation of the specific conductance of a solution (SC), see calculate the specific conductance of a solution with PHREEQC . The specific conductances calculated from the current with Ohm's law and from SC are then identical, as illustrated in Figure 6. The option is invoked when 'true' is added to the line with identifier multi_D:
see file SC_Ohm.phr.

with Ohm's law and from SC are equal when the diffusion coefficients are corrected for ionic strength. The red line with the XCross symbols shows the current`s specific conductance when the diffusion coefficients remain constant at their infinite dilution values. The magenta crosses show measured SC`s. 
Finally, can the code calculate uphill diffusion that results from the charge imbalance due to different diffusion coefficients, as shown in the upgradient example?
We can define in keyword TRANSPORT that the current is fixed, and PHREEQC calculates transport with the NernstPlanck equation, adapting the potentials in the cells to obtain the current. Fixing the current to a small value, say 1e10 Ampère, is identical with zerocharge transfer and can be used to compare the two calculation types. Alternatively, the implicit algorithm can be used, setting implicit true 5 in keyword TRANSPORT, with 'max_mixf' = D * Δt / Δx^{2} = 5. The time steps then are more than 10 times larger, and the calculations 10 times faster.
The results of running file uphill_NP.phr are shown in Figure 7.
Figure 7. (Left) Concentrations of H^{+}, NO_{3}^{}, Na^{+} and Cl^{} in a column showing uphill diffusion of Na^{+} and Cl^{}. (Right) The potential in the column becomes negative because H^{+} diffuses faster out than NO_{3}^{}.
Initially, the concentrations of Na^{+} and Cl^{} are the same in the column and the boundary solution, while the concentrations of H^{+} and NO_{3}^{} are higher in the column. H^{+} diffuses faster than NO_{3}^{} from the column, and a negative potential develops that pushes Cl^{} out, and attracts Na^{+} from the boundary solution. The calculated potential is as low as 23 mV after 20 minutes.
The concentration pattern shown in Figure 7 is identical with what is calculated by imposing zerocharge transfer (also shown when the file is run), although the formulas used in the calculations are different. When eqn (1) is solved with zerocharge transfer, the potential is omitted from the calculations. On the other hand, the complete calculation with the NernstPlanck equation shows that a potential of 23 mV develops in the column that may well be measurable experimentally.
The input files current1.phr, current2.phr, SC_Ohm.phr and uphill_NP.phr use features for calculating electromigration with PHREEQC: