En route to a unified model for photo-electrochemical reactor optimisation. I - Photocurrent and H 2 yield predictions †

A semi-empirical model was developed for prediction of photocurrent densities and implemented to predict the performance of a photo-electrochemical reactor for water splitting in alkaline solutions, using Sn IV -doped α -Fe 2 O 3 photo-anodes produced by spray pyrolysis. Photo-anodes annealed at different temperatures were characterized using photoelectrochemical impedance spectroscopy, cyclic voltammetry in the presence and absence of a hole scavenger and also the open circuit potential under high intensity illumination. Mott-Schottky analysis was used cautiously to estimate charge carrier concentration and the flat band potential. In addition to overpotential/current distribution and ohmic potential losses, the model also accounts for absorbed photon flux, surface and bulk electron-hole recombination rates, gas desorption, bubble formation and (H 2 -O 2 ) cross-over losses. This allows the model to estimate the total yield of hydrogen, charge and gas collection efficiencies. A methodology is presented here in order to evaluate parameters required to assess the performance of a photo-electrochemical reactor in 1D and 2D geometries. The importance of taking into account bubble generation and gas desorption is discussed, together with the difficulties of measuring charge carrier concentration and electron-hole recombination in the bulk of the semiconductor, which are of major importance in the prediction of photocurrent densities.


Introduction
The well-known impact of CO2 emissions on the climate change, human health and local environments are driving the search for more benign sources of energy. The sun is the earth's ultimate energy source, though subject to diurnal and seasonal intermittency. This calls for a technology to store the harvested solar energy, and make it readily available when needed. Batteries have been proposed as a practical solution. Their cost, potential toxicity and limited specific energies, contrast with their portability and versatility. As an alternative, solar fuels, such as hydrogen, have been proposed as a suitable energy carrier. Recent advances in technologies such as water splitting via PV+electrolyser, 1 integrated photo-electrochemical, 2 photoelectrochemical 3 and particulate photocatalysts 4 are approaching a level of readiness for up-scaling. However, commercial deployment is still constrained by the degradation of the systems with time and the cost of the produced hydrogen, which must be competitive with current technologies, e.g. steam reforming. Recently, most of the efforts have been directed towards the development and prediction of performance of photo-anodes 5 and in a lesser extent towards photo-cathode development. 6,7 Hence, technoeconomic modelling and engineering optimisation need to be tackled in parallel with addressing scientific and engineering challenges. The seminal review of Gerisher 8 mapped the path towards a fundamental understanding and modelling of semiconductors coupled with electrochemical devices. Carver et al. 9 and Haussener 10 started proper electrochemical modelling of such devices, often operating as an electrolyser at relatively low current densities of 2 -200 A m -2 , c.f. 1000 -3000 A m -2 for conventional water electrolysers. 11 Berger et al. 12,13 expanded this work for membrane based systems. However, most of the parameters that defined material performances were from the literature. The effects of bubbles and desorption of gas were not quantified and there was a limited comparison of predictions with experimental data. A more complete integrated solar driven water-splitting device operating at neutral pH was modelled by Jin et al. 14 This model considered pH gradients and cross-over losses in presence and absence of a membrane. Little attention was paid to the effects of bubble generation, causing light reflection/scattering, decreasing conductivities of electrolyte solutions and spatial distribution of current densities of larger electrodes. However, it showed clearly the issues related to current density distribution even at small size electrodes (16 mm). More recently, this model was expanded to study the effect of electrolyte flow on membrane based systems. 15 Recent works by Haussener et al. [16][17][18] extended their model for temperature dependence under concentrated irradiation and statistical analysis validated against experimental data. The effect of perforations on the performance of integrated solar devices 18 and photoelectrochemical reactors has also been reported recently. 19 However, hitherto a model had yet to be developed for photoelectrochemical reactors, in which the photo-absorber and electrode are the same material, coupled with bubble generation at the electrode and bulk of electrolyte. Incorporating bubble formation at the photo-electrode | electrolyte interface is important as: (i) it enables quantification of the amount of gas produced near the surface and the actual dissolved gas flux diffused into the bulk of the electrolyte that could take part in cross-over losses; 20 (ii) rates of generation of bubbles in photo-electrochemical reactors should be minimised by judicious design, to minimise resistive losses 21,22 and the reflection and scattering of incoming light 23 that otherwise would decrease photo-electro-active areas. The latter effect will be addressed in more detail in a subsequent paper.
Validation against experimental results obtained with different photo-electrodes materials is often bypassed. In this paper, we extend previous work, 19 aiming to predict the performance of up-scaled photo-electrochemical reactors based on easily measurable characteristics of components of the system at small scale, specifically photo-anodes. Sn IV -doped α-Fe2O3 photo-anodes treated at different annealing temperatures were used to obtain parameter values as inputs to a 1D model and validate it against experimental data measured in small scale cells. The novelty of the work presented here rests on the use of a single equation to predict photo-anode photocurrent densities, considering electron-hole recombination rates in the bulk semiconductor, the effects of bubble generation and diffusion of dissolved gases in the electrolyte, as outlined in (i) in the previous paragraph. Next, a hypothetical and enhanced photo-anode was studied in an 2D model. This will be followed in part II of this series of papers by 3D geometrical optimisation of perforated photo-anodes.

Definition of bulk and surface electron-hole recombination
Definitions of bulk and surface electron-hole recombination rates in semiconductors are widely spread and vary depending on the framework in which they are discussed. For clarity, we define bulk recombination as the annihilation of generated holes and electrons in the bulk of the semiconductor, outside the depletion region or space charge layer. This recombination accounts for low conductivity or sluggish charge transfer into the substrate, as well as for absorbed photons in the bulk that ultimately do not generate charge. Surface recombination is closely related to the rate of charge carrier extraction at the surface of the semiconductor and the rate at which the charge is migrating and diffusing from the depletion region to and across the semiconductor | electrolyte interface. Generally, surface electron-hole recombination depends on the chemistry of the surface and is a function of the applied potential. A schematic view of these definitions is shown in Fig. 1.

Photocurrent density prediction
The photocurrent density produced in the anode was predicted using a modified version of the Gärtner-Butler equation, 19,24,25 full derivation in supporting information, which does not account for charge migration and diffusion rates into the semiconductor bulk, where electron-hole recombination may occur.
where ( 0 − ) is the actual absorbed light after attenuation e.g. by electrolyte solutions, membranes, bubbles, etc. The photocurrent density should also be limited by the photon flux ( 0 − ) as: This formulation was used previously and explained in more detail by Hankin et al. 19 According to Dotan et al., 26 photocurrent densities should also be corrected for surface and bulk recombination rates as: absorbed is a value that does not depend on the potential and is the maximum achievable current density that corresponds to × ∑ ( 0 − ) . Hence,  bulk is a function of potential. We propose here a modified version of this equation to decouple the dependence of  bulk from the potential and, instead, account for this dependence in the photocurrent as predicted by the Gärtner-Butler equation: In which only ph(GB3) and surface are a function of the band bending Δ SC ., whereas bulk is independent of potential, as will be described below.
Surface recombination. The interfacial charge transfer efficiency,  surface , can be estimated from photoelectrochemical impedance spectroscopy (PEIS) as a function of potential 27-29 as described previously for Sn IV -doped α-Fe2O3. 19 Then,  surface can be estimated as a function of the potential after fitting the impedance spectra to a two-time constant equivalent electrical circuit, Rs(QSC[Rt(QrRr)]). 27,30 This simplified circuit removes the charge transfer resistance in the bulk, which cannot be discriminated easily from the spectra, and assumes the charge transfer is dominated by surface states. 30 This circuit explains well the presence of two processes: interfacial electron-hole recombination and charge transfer at the photoanode | electrolyte interface under illuminated conditions. The fitting also accounts for surface distribution after replacing capacitances by constant phase elements. 31 The interfacial charge transfer efficiency then can be written as: where Rt is the charge transfer resistance and Rr the recombination resistance. In order to express  surface as a continuous function of potential, the values were fitted to an exponential equation in terms of constants A and B, and subsequently combined into the model.
where, Δ SC is the driving force for photocurrent defined here as the difference between the flat band potential and the applied electrode potential, Δ SC = applied − fb .
Bulk recombination. As defined here, the recombination in the bulk should not be a function of applied potential. If surface recombination is avoided, that means  surface = 1 when using a fast hole scavenger such as H2O2, see reaction (7),  bulk can be determined as the ratio between the experimental photocurrent ph,H 2 O 2 and the expected photo-current from the modified Gärtner-Butler equation, ph(GB3) .
Then, the rate of bulk recombination should be determined by the Shockley-Read-Hall effect, which does not depend on the applied potential, but on the charge carrier concentration in the bulk of the semiconductor. As shown in Fig. 1, photons absorbed deeper into the semiconductor than the depletion layer tend to recombine, due to the absence of an electric field, so produce no net charge. The Gärtner-Butler equation already considers the depletion layer thickness as a function of the band bending, so the ratio of photocurrents in absence of electronhole recombination and expected photocurrents from the modified Gärtner-Butler equation, should be constant. This can be tested by comparing the absorbed photons in the depletion layer and total number of photons absorbed by the whole depth of the semiconductor, ( 0 − ) . The former photon flux can be estimated using the semiconductor Debye length as an approximation for the depletion layer width in a nanostructured electrode, , eqn (9). 32 Then, the absorbed photons inside the region of the Debye length can be estimated by a modified Gartner relationship, eqn (10). 25 Using previously reported parameters for charge carrier concentration and absorptivity of annealed Sn IV -doped α-Fe2O3, 33 the number of absorbed photons inside the Debye length is only 7.3% of the total absorbed photons. This value is comparable to the values measured here for  bulk of annealed samples, ca. 7%.

Dark (anodic and cathodic) current
Dark anodic current densities were approximated by the Tafel equation, as dark overpotentials are usually high for oxygen evolution on hematite electrodes (>100 mV), so the rate of the reverse reaction could be neglected: ,dark = 0, 10 / The overpotential for oxygen evolution reaction (OER) was computed against its equilibrium potential: Assuming small changes of local pH in the electrolyte and insignificant change in the partial pressure of oxygen, it can be assumed that O 2 |H 2 O (RHE) = 1.229 V. Likewise, cathodic currents were computed with the Butler-Volmer equation for hydrogen evolution on platinized titanium.
The hydrogen overpotential was calculated from the equilibrium potential for the hydrogen evolution reaction (HER) at 298 K: Again, assuming small changes of local pH and negligible change in partial pressure of hydrogen, it was assumed that

Ionic current
As reported previously, 19 the ionic current in the electrolyte for (infinitely dilute) charged electro-active species was modelled after substituting Faraday's law into the Nernst-Planck equation: Assuming electroneutrality conditions: and in the absence of concentration gradients, the distribution of potential satisfies Laplace's equation: The mass transport limited current density for a quiescent 1 M NaOH electrolyte is L,OH − 10 4 A m −2 . 19 All experiments and predicted current densities in this study were < 50 A m -2 . Hence, the OHgradients and ion depletion at the photo-anode | electrolyte interface were such that ∇ = 0 was a suitable assumption.
Assuming negligible resistance for the electric contacts and wiring in the cell, a potential balance across the cell gives: 9 Where O 2 |H 2 O − H 2 O|H 2 is the required thermodynamic potential difference; + | c | the overpotentials at the anode and cathode, respectively, and the last term is the sum of the ohmic potential losses in the phases with conductivity , path length and cross section area .

Bubble formation and gas desorption
The evolution of gas at the surface of anode and cathode was modelled in terms of efficiency of gas evolution, . 20 where is the flux density of hydrogen or oxygen produced by the faradic current at the surface of the electrode.
is the fraction of converted into bubbles at the electrode surface, desorption is the gas evolving from natural desorption under supersaturated conditions and bulk is the diffusive flux density. A schematic representation of these processes is given in Fig. 2, following definitions proposed by Vogt. 20 The effects of bubble generation on light scattering and reflection has yet to be implemented. Also, the present model does not consider the effect of flow, which can decrease gas evolution efficiencies substantially 20 and modify the micro-and macro-convection effects due to the removal of bubbles at the electrode surface.
Efficiency of gas evolution: several models have been proposed to estimate in terms of measurable quantities such as current density and fractional surface coverage. 20,[34][35][36] However, these models are difficult to implement when coupled to other physics, such as electrochemical reaction and transport of dissolved species. Root squared and exponential expressions often result in complex values when evaluated under certain conditions. Vogt was wary to use simplified and ad hoc models, but also was aware of difficulties found when using fundamental models. 35 A phenomenological model for gas evolution proposed by Janssen et al. 37 relates closely to the system studied here. Hydrogen and oxygen bubble efficiencies at nickel and platinum-wire electrodes in alkaline solutions over current densities of 10 to 10 4 A m -2 were predicted by a power function: where , and , are fitted parameters from experimental data for a given system. The measurements by Janssen et al. differ mainly in two aspects: the electrode material was a nickel wire and the fluid was not quiescent, but flowing at 0.12 m s -1 (Re ≈ 67). Also, eqn (20) does not account directly for bubble coverage which could have a significant effect, as demonstrated previously for electrolysers operating at current densities > 100 A m -2 . 38 As explained in the supporting information, relatively low current densities (0 -100 A m -2 ), as reported here, produced fractional bubble coverages < 0.1 that do not affect gas evolution efficiencies significantly. 20 For increased current densities, i.e. higher photon flux intensities or improved photoelectrodes, the effect of bubble coverage should be taken into account. However, eqn (20) is still a reasonable approximation that enables estimation of gas evolution efficiencies as a function of local current density. A more detailed comparison between Vogt and Janssen models can be found in the ESI, Fig.  S2. The semi-empirical eqn (20), 37 for the estimation of , was used in the model described in here.
Volumetric mass transfer coefficient: gas desorption was modelled by a simple mechanism of volumetric desorption when the specie is above the saturated concentration.
where, is the mass transport coefficient and the area of transfer per unit volume. Usually, these values are reported as , the volumetric mass transfer coefficient. The seminal work by Frossling 39 has been used widely 40 to estimate as a function of Re and Sc numbers for bubbles with diameters 0.1 to 2 mm.
The minimum average bubble diameter for mass transport in quasi-static conditions in a variety of electrolytes, including 1 M KOH, was found to be = 4.5 × 10 −4 m. 41 This value can be used to estimate the bubble slip velocity (relative to the liquid) using the Wüest relationship. 42,43 = 4474 ( 2  ) 1.357 Then, Re and Sc are calculated as usual for a spherical geometry: Where is the kinematic viscosity of the fluid and the diffusion coefficient of the oxygen or hydrogen in the liquid, see a detailed list of used parameters in the ESI,  (26) where is the gas holdup, which can be estimated graphically as a function of 45 or using the correlation by Joshi and Sharma: 40,46 = 0.3+2 (27) Table S2 in the ESI provides results of a detailed numerical calculation of , which has values of 1.94×10 -4 s -1 for oxygen and 3.02×10 -4 s -1 for hydrogen at 298 K and in quiescent solutions. These values are similar to the predictions by Wüest 43 (2.76×10 -4 s -1 ), but are significantly lower than those predicted by Deckwer 45 for oxygen absorption and desorption in tap water for a well stirred reactor (3.58×10 -2 and 9.96×10 -2 s -1 , respectively). The differences with the latter values were due to the absence of electrolyte flow and convection effects in the reactors.

Diffusion of dissolved species
Dissolved gases, oxygen and hydrogen, were modelled by the convection-diffusion equation: where is the concentration of the specie , is the diffusion coefficient, is the velocity of the fluid, and is the reaction Please do not adjust margins Please do not adjust margins rate or in this case the rate of desorption. Assuming absence of homogeneous reactions and convection effects, the equation simplifies to: The flux density of species at the surface of the electrode was modelled according to Faraday's law of electrolysis and reaction stoichiometries for water splitting in alkaline conditions, hydrogen (HER) by reaction (33) and oxygen evolution (OER) by reaction (34), which can also be expressed in terms of photogenerated oxidising species, ℎ + , by re-allocating the charge on the right-hand side of the reaction.
In the case of electrolyte | membrane interfaces, the concentrations of oxygen and hydrogen at the interface were modelled using a partition coefficient.
The partition coefficient was estimated from the saturated concentrations in the membrane and electrolyte, as follows: In some cases, OHions were also included with an initial concentration of 1 M. However, the depletion and increase in concentration at the anode and cathode, respectively, prevented the model converging to a steady state solution during long term simulations, specifically when a cationpermeable membrane was used to separate the catholyte and anolyte. A quasi-steady state was considered, i.e. depletion rates of OHin the bulk was very slow compared to rates of other processes in the model.
Lastly, the flux densities of hydrogen and oxygen, , is defined as the total amount of the gas coming from bubble formation at the surface of the electrode and desorption in the electrolyte, or catholyte and anolyte when a membrane is used, according to eqn (37). Fluxes were normalized by the area of the electrodes, defined as hydrogen flux density and reported as mol m 2 s -1 . The faradaic efficiency, Φ , for hydrogen and oxygen evolution is defined as usual: the ratio between the charge transferred for hydrogen evolution (HER) or oxygen evolution reaction (OER), and the total charge transferred at the electrode, according to eqn (38). The collection efficiency,  collection, , is defined here as the ratio between the flux density and the maximum achievable gas evolution rate given by the total current at the electrode and Faraday's law, according to eqn (39). The collection efficiency is intrinsically related to the cross-over of products between anolyte and catholyte, via the membrane if present.
An overall efficiency for a photo-electrochemical device was estimated following the definition proposed previously 19 for photo-electrodes operating under an electrical bias: assisted solar to hydrogen efficiency (ASTH). This efficiency is defined as the power associated to the rate of hydrogen evolution (hydrogen flux density) divided by the power input to the system (combination of photons and electrical bias in this case):

Boundary conditions
Boundary conditions for the 1D model are described in Fig. 3.
Photo-and dark current densities were calculated at the boundary electrode | electrolyte. The potential applied through the system was set by the potential at the anode vs. the reference potential (RHE at the cathode). The spatial distribution of electric potential, eqn (17), is then solved between the electrodes in such a way that the cell potential difference is resolved for Δ bias = − . Similarly, the flux density of species, associated with currents by the Faraday's law, was coupled to transport of species at the electrode | electrolyte interface. Bubble evolution and gas desorption rates were estimated after volumetric integration along the electrolyte to compute the rate of oxygen and hydrogen collection in the anolyte and catholyte, respectively. The oxidation of hydrogen and reduction of oxygen at the anode and cathode, respectively, were assumed to be controlled by mass transport, i.e. assuming concentrations of these species at the electrode is zero. The concentration in the membrane was limited by the solubility of gas in the membrane (saturation concentration).
The same boundary conditions can be applied for 2D and 3D geometries, with different distances between electrodes, membrane and perforations and different electrode arrangements. For a simplified 1D model, the electrode surfaces were facing each other and 0.02 m apart, the membrane was 330 µm thick and positioned equidistant from the 740 µm thick electrodes.

Model summary
The prediction of photocurrent densities, accounting for interfacial and bulk electron-hole recombination rates, and the efficiency of gas evolution were modelled by semi-empirical correlations. The prediction of dark current densities, diffusion and volumetric desorption of dissolved gases were modelled using well-established functions, based on fundamental properties of the system.

Photo-anode fabrication
Sn IV -doped α-Fe2O3 photo-anode films deposited on titanium and FTO were fabricated by spray pyrolysis, following a similar procedure reported previously 19,53 . In brief, the precursor solution comprised 0.1 M FeCl3.6H2O (99.99%, Sigma Aldrich, UK) and 6×10 -4 M SnCl4 (99.995 Sigma Aldrich, UK) dissolved in ethanol absolute (AnalaR Normapur, VWR), with Sn IV concentrations corresponding to 1.3% mass doping relative to Fe III . The precursor was nebulized 40 times with a quartz spray nozzle (Meinhard, USA) controlled by WinPC-NC CNC Software (BobCad-CAM, USA) on the substrates pre-heated at 480 °C. Post-deposition annealing was done in a pre-heated oven (Elite Thermal Systems, Ltd.) in air at atmospheric pressure at 400 and 500 °C over an hour.

Photo-anode characterization
A photo-electrochemical cell of 0.06 dm 3 (PVC body with a quartz window) was used for electrochemical characterization of photo-anodes. A potentiostat/galvanostat (Autolab PGSTAT 30) was used to control the three-electrode cell. Hematite samples acted as working electrodes, a platinised titanium mesh (Expanded Metal Company, UK) as counter electrode and HgO|Hg as reference electrode (0.913 V vs. RHE and 0.109 V vs. SHE). 1 M NaOH was used as electrolyte (pH 13.6) at room temperature, ca. 25 °C. Hydrogen peroxide, H2O2 (30%, AnalaR Normapur, VWR) was added in some cases as a sacrificial electron donor (hole scavenger) at a concentration of 0.5 M. A Xe arc lamp (LOT-Oriel) of 300 W with Fresnel lenses was used to obtain a beam of 3 mm × 6 mm and a maximum power intensity of 3646 W m -2 . The intensity of light was corrected by electrolyte and quartz attenuation to give an effective power intensity of 3000 W m -2 on the photo-anode. The light source was characterised and calibrated with a UV-Vis spectrophotometer coupled to a CR2 cosine receptor (Black-Comet CXR-25, StellarNet, USA). Absorption spectra of FTO samples were obtained using a Shimadzu UV-2600 UV-Vis spectrophotometer with an integrating sphere and corrected by baseline subtraction assuming reflection and scattering are wavelength independent. 54 Thickness of the films deposited on FTO was measured using a stylus profilometer (TencorAlphastep200 Automatic Step Profiler). For the photoelectrochemical cell described here, i.e. small beam size to volume ratio, there was no significant increase in temperature after long term exposure to illumination.

Simulation parameters
Comsol Multiphysics 5.2a with a Batteries and Fuel Cells module was used to solve the system using the finite element method. Secondary current distribution (siec) and transport of diluted species (tds) physics were coupled; gas desorption was treated as a homogeneous reaction in the electrolyte. A stationary nonlinear solver with a relative tolerance of 0.001 and a maximum of 100 iterations was typically used, as with previous works. 19 The mesh was geometry dependent with a minimum element size of 4 × 10 -7 m and maximum size of 2 × 10 -2 m; the mesh was optimized for higher resolutions in narrow regions. The model converged typically in less than 20 seconds (1D) or one minute (2D) for a given set of conditions. A parametric sweep was used to evaluate a wide range of settings and geometries. The suitability of the mesh size and tolerance is demonstrated in the ESI.

Results and discussion
Photocurrent response of photo-anodes: Fig. 4 shows linear sweep voltammograms of Ti|Sn IV -doped α-Fe2O3 samples. Samples were characterised in alkaline conditions in presence and absence of a hole scavenger, 0.5 M H2O2. In the absence of H2O2, the effect of annealing on the surface electron-hole recombination was evident in the shift of photocurrents to lower potentials. This effect disappeared once photo-anodes were immersed in electrolyte containing hydrogen peroxide. Maximum photocurrents were also affected by the annealing process and were comparable in the presence and absence of the hole scavenger. The decrease in the maximum photocurrent with annealing could be associated with an increase of electronhole recombination rates in the bulk of the semiconductor. The migration of the dopant to the surface after heat treatment seemed to have been responsible for this behaviour, as discussed in more detail previously. 33 Mott-Schottky analysis was used to estimate the charge carrier concentration of the semiconductor, as shown in Fig. S3  The flat band potential was estimated using several methods. As shown in Fig. S4 and Table S5 in the ESI, good estimations can be obtained from open circuit potential (OCP) measurements at high intensity or from the linear portion of the potential dependence of the squared photocurrent in presence of a hole scavenger. As spray-pyrolysed hematite samples were nanostructured, 33 they violated one of the assumptions of Mott-Schottky analysis: the electrode | electrolyte interface should be perfectly planar. 56,57 Therefore, Mott-Schottky analysis should not be relied upon implicitly for determination of flat band potentials of such materials. The absorbance and absorptivity were calculated using a film thickness of 47 ± 6 nm and were very similar for all samples. The photon flux from the light source was integrated between wavelengths of 280 and 800 nm to obtain 0 . Next, the absorbed photon flux 0 − and the product between this and the absorptivity ( 0 − ) were also integrated to be used in the Gärtner-Butler equation. A typical band gap of 2.11 eV was found for all hematite samples and determined from absorbance measurements using Tauc plot analysis. 54 All the spectra are reported in Fig. S5 to S7 in the ESI.
Electron-hole recombination rates at Ti|Sn IV -doped α-Fe2O3 | 1 M NaOH interfaces were extracted and analysed from photoelectrochemical impedance spectroscopy (PEIS) measurements, following the methodology proposed by Peter et al. [27][28][29] Fig. 5 shows typical impedance spectra, in which the charge transfer and recombination processes are evident at ca. 1 kHz and 1 Hz, respectively. After extracting the resistances corresponding to charge transfer, , and recombination processes, ,  surface was fitted to eqn (6) as a function of the potential and fitting parameters A and B, see Fig. 6.

Electron-hole recombination rates in bulk Sn IV -doped α-Fe2O3
were derived using eqn (8) and the charge transfer efficiency ( bulk ) plotted in Fig. 7 as a function of the applied potential.  bulk was not wholly independent of applied potential, though the dependence was considerably less pronounced compared with  surface . At lower potentials,  bulk could have been affected by residual surface recombination rates; hence, it is proposed that the values of  bulk asymptote to their actual values as  surface asymptote to their maximum values. In order to estimate a constant and potential independent  bulk for each sample, an average value was calculated in the range between 1.4 and 2 V vs. RHE.  bulk can also be estimated from the photo-response in the absence of a hole scavenger (1 M NaOH), but after correction for surface recombination, i.e. dividing by  surface . Similar values were found with this approach; data from the two methods are compared in the ESI ( Fig. S8 and Table S7). The model uses the former values, in presence of H2O2, as it is a more direct method, and does not require mathematical corrections due to interfacial electronhole recombination rates in absence of H2O2.

1D Model
A summary of all measured parameters input into the model and the techniques used for their determination is presented in the ESI, Table S3 and S4.
Photocurrent prediction. Fig. 8 shows the effect of electrode potential on predicted photocurrent densities for different versions of the Gärtner-Butler equation: eqn (1), (2) and (4). The maximum achievable photocurrent density was limited significantly by absorbed photon fluxes, as expected. After correction for surface and bulk electron-hole recombination rates, j-E data exhibited behaviour typical of hematite photoanodes. Predicted dark current densities were very low for potentials < 1.6 V vs. RHE, in agreement with experimental data shown in Fig. 4. Fig. 9 shows good agreement between predicted and experimental current densities (photo + dark), with regression analysis giving R 2 values of 0.987, 0.994 and 0.935 for Ti|Sn IVdoped α-Fe2O3 samples unannealed and annealed at 400 and 500 °C, respectively. All values input to the model were obtained from independent experimental measurements.
Surface recombination, 27 absorbed photon flux 54 and flat band potential 58 were obtained using well documented methodologies. However, charge carrier concentration values estimated from Mott-Schottky analysis can vary by several orders of magnitude from sample to sample, especially when the semiconductor surface is contaminated and nanostructured. The estimation of bulk recombination rates from measurements in the presence of a hole scavenger, supressing surface recombination, has a known origin as discussed previously in the theoretical treatment section, but still lacks a complete fundamental explanation. For modelling purposes and assuming these properties do not change drastically with the size of the electrode, predicted photocurrent densities can be used to predicted the behaviour of up-scaled reactors using photo-electrochemical results obtained at small scale. In practice, larger solar simulators may produce photon flux inhomogeneities.
The difference between predicted photocurrent densities and experimental results can be associated with several causes: (i) uncertainties in the measurement of photo-anode properties. Difficulties in the fitting of the impedance spectra might alter the estimated interfacial electron-hole recombination rates; this seems to be the main reason of disagreement between predicted and experimental values for the sample annealed at 500 °C. The higher the recombination rate, the easier to discriminate it from the charge transfer rate. This would also explain why the predictions for the other samples were more accurate. (ii) Irradiated photon flux can be affected by small changes in the position of the cell and decrease of intensity with time due to degradation of the lamp. (iii) Dark current densities, predominant only at higher current densities, were dependent on the electro-active area of the electrode, which could also vary from sample to sample.

Hydrogen yield predictions in membrane-less configurations.
Hydrogen yield, faradaic and collection efficiencies were predicted for all samples in a 1D membrane-less model. Experimental data for comparison with these predictions are not available, since anode and cathode were immersed in the same electrolyte with no gas separation and the small-scale reactor for photo-electrochemical characterisation was not gas tight. Fig. 10 shows that there was a potential for each sample at which efficiencies started to increase towards unity. These potentials correspond to the current density ( total ≈ 0.04 A m −2 and hydrogen flux densities of 1 × 10 −9 mol m −2 s −1 for conditions used) at which bubble generation rates became significant at the surface ( > 0.01) and desorption of gas due to the saturation near the electrode started to contribute to the hydrogen flux density, as shown in Fig. S9 in the ESI. Hence, H2 gas is not effectively desorbed from the electrolyte and collection efficiencies (continuous lines in Fig. 10) were near zero. Faradaic efficiencies (dashed lines) were never below 0.4, which is the predicted value for the configuration at which reduction of oxygen occurred at the same rate as reduction of water on the cathode; this value can be increased almost to unity if a membrane is used. The small lumps before the efficiencies start increasing are artefacts associated with a discontinuity (additional hydrogen flux density due to the sudden formation of bubbles) in the model.
With increasing photo-anode potential, all efficiency values approached unity asymptotically. Hence, the contribution of bubble generation and desorption rates should be considered when operating at low current densities. However, once a certain photocurrent density threshold is reached, depending on the geometry and photo-electrode properties, the effects of potential on the faradaic and collection efficiency are minimal. Also, it is noteworthy that the collection efficiency estimated for a membrane-less reactor was estimated assuming that oxygen and hydrogen are separated effectively in a subsequent process and that the explosion limit of > 4 vol% O2 in H2 is not exceeded. 14  Hydrogen flux densities followed a similar pattern, as shown in Fig. 11, which plots fluxes against both current densities (bottom x-axis) and electrode potential (upper x-axis). As expected, the dependence of hydrogen flux densities on current density were independent of the photo-anode, the properties of which determined its current density-potential relationship, so hydrogen fluxes were ultimately highly dependent on the potential. At current densities less than 0.04 A m -2 , bubble generation and desorption rates were predicted not to contribute to hydrogen fluxes, as most of the hydrogen was notionally oxidised at the anode, resulting in hydrogen fluxes close to zero.
Hydrogen yield predictions using a membrane. Faradaic and collection efficiencies for hydrogen evolution were estimated for the case of a cation-permeable membrane (330 µm) positioned between the electrodes. Nafion was chosen due to its chemical stability, durability and superior mechanical properties over an anion-permeable counterpart, which otherwise would have been more suitable for water splitting under alkaline conditions. However, OHconcentrations were high enough to ensure slow depletion and negligible pH changes with time. As shown in Fig. 12, faradaic efficiencies were higher Please do not adjust margins Please do not adjust margins and almost unity at any given potential for all samples. This demonstrates the utility of using a membrane or micro-porous separator to separate gases when operating at low current densities. On the other hand, collection efficiencies were predicted to be lower when compared to the membrane-less model, because predicted hydrogen fluxes came only from the catholyte, which was then only half of the electrolyte volume. However, the use of a membrane or micro-porous separator is a more practical, as it obviates the need of further gas separation and circumvents safety issues related to hydrogen and oxygen gas mixtures. 14 The discontinuities in collection efficiencies at ca. 1.3 and 1.4 V was due to gas saturation of the membrane, in which transport of hydrogen and oxygen occurs by diffusion only. Hydrogen fluxes were predicted to be slightly greater for a membrane reactor, i.e. shifted to lower current densities, but followed a similar behaviour to that shown in Fig.  11 for a membrane-less model. In this case, the tipping point of current density for gas evolution at the cathode and desorption rate was 0.025 A m -2 . Hydrogen fluxes at greater current densities were the same as predicted in a membrane-less model. This means that once the current setup is operating at current densities at which gas is evolving via bubbles or desorption for saturated electrolytes, then a membrane or microporous separator is needed only to separate the products and prevent H2-O2 explosions, but otherwise is predicted not to affect the reactor performance significantly. As proposed by others, 10, 13 instead of using membranes, a material with pores smaller than the Sauter diameter of bubbles, usually between 50 59, 60 and 450 µm 41 , may be sufficient to separate H2 and O2 gases. However, this will depend heavily on the (i) thickness of the microporous separator; (ii) the superficial and morphological properties of the electrodes and its effect on the size and life time of the adhered bubbles; (iii) the flow rate of the electrolyte, which affects the flow regime and physical separation of the products, which in turn is also affected by the distance and characteristic length of the electrodes.

2D Model
A simple monolithic photo-electrode was modelled in order to confirm the importance of geometrical optimisation, as reported previously. 19 A membrane was absent for this case. The current density distribution and hydrogen fluxes of a hypothetical enhanced photo-anode material and a platinized cathode were predicted. The properties of the photo-anode were notionally enhanced in the model, to demonstrate the importance of geometrical optimisation, as current density distributions become less homogeneous and product crossover rates are exacerbated on scale-up. The properties of this hypothetical photo-anode are described in the ESI, Table S8. Surface and bulk recombination were minimized; its band gap was deceased to 1.7 eV, so that the wavelength corresponding to the photon absorption edge was increased to 729 nm. Photon fluxes were taken from the solar AM1.5D spectrum. A membrane was not included in the model, but hydrogen fluxes were integrated only over the 1 M NaOH catholyte. The electrolyte volume (projected as 2D) was proportional to the electrode length on both dimensions. Fig. 12. Predicted effects of electrode potential faradaic (dashed) and collection efficiencies (continuous) as a function of potential for hydrogen evolution on Ti|Pt when using a photo-anode Ti|Sn IV -doped α-Fe2O3 | 1 M NaOH and separated by Nafion membrane. Fig. 13 shows the resulting current density distribution for 1 V vs. RHE electrical bias as a function of the relative distance to the edge of the electrode. A significant decrease of local current densities (mostly photocurrent densities at this potential) was predicted along the photo-electrode, especially for electrodes with characteristic lengths > 0.02 m. These results are not as severe as the reported current distributions for neutral pH. 14 After integration of current densities on the anode and hydrogen fluxes in the catholyte, collection and faradaic efficiencies were also predicted, as shown in Fig. 14. Total current densities (normalized by the length of the electrode) were greater for shorter electrodes, as implied by the current density distributions shown in Fig. 13. However, hydrogen fluxes were smaller for shorter electrodes and reached a maximum at ca. 0.01 m, with a subsequent decrease for longer electrodes. A similar behaviour was predicted in ASTH efficiencies as a function of length. This means that due to crossover losses and current density distributions, there is an optimized length at which photo-electrodes should be fabricated. Faradaic and collection efficiencies are also reported in Fig. 14.
As expected for electrodes shorter than 0.03 m, the cross-over losses and reduction of oxygen at the cathode, causing loss of faradaic efficiencies, were significant. For longer electrodes, both efficiency values approached unity asymptotically. Electrochemical reduction of oxygen at the cathode occurred close to the cathode edge, at which overpotentials were highest, decreasing towards the centre. Hence, longer cathodes will have less area at which oxygen reduction occurs, so increasing average faradaic efficiencies and gas collection efficiencies. The consequences of current density distributions and cross-over losses for ASTH efficiencies and hydrogen flux data in Fig.14 emphasise the importance of geometrical optimisation.
possible to use this model to define optimized conditions to achieve the highest ASTH and hydrogen fluxes. More detailed results and comparisons with other efficiency definitions (ABPE and AB-STH) can be found in the ESI, Fig. S11.  For the case presented above, product-gas cross-over, i.e. H2 flux in the anolyte and O2 in the catholyte, was not as severe as reported by Jin. 14 For small electrodes (< 0.01 m), product-gas cross-overs were ca. 2% to 10% molar for oxygen and hydrogen in the catholyte and anolyte, respectively. This still represents a safety issue and increased risk associated with ignition of stored products. 10,14 For electrodes longer than 0.1 m, the predicted product-gas cross-over was < 2% in both cases. If a membrane is used, the amount of hydrogen collected in the anolyte reached zero values, as predicted by the 1D model, due to the low fluxes of hydrogen and oxygen crossing between anolyte and catholyte are effectively anodically oxidised and cathodically reduced, respectively. A more detailed discussion of the effects of membranes on the geometrical optimisation of photo-electrodes will be reported in part II of this series of papers.

Conclusions
A model for photo-electrochemical reactors was developed to describe the behaviour of Sn IV -doped α-Fe2O3 photo-anodes deposited on titanium and annealed in air at different temperatures. Photo-anodes were characterised in terms of photon absorption, electron-hole recombination rates in the bulk semiconductor and at the semiconductor | electrolyte interface, and other electrochemical parameters obtained experimentally. Charge carrier concentration and bulk recombination rates were found to require special care and theoretical treatment for their estimations. Predicted values showed good agreement with experimental data. For the first time, bubble formation rates at the electrodes and in the bulk of electrolyte solution for supersaturated conditions, were predicted as a function of current density; hence, reactor performance did not increase drastically by including a membrane or microporous separator when operating at current densities > 0.04 A m −2 , as shown in Fig. 11. This critical current density depended mainly on the geometry and properties of the photo-anode. Faradaic efficiencies for hydrogen formation reached unity when the photo-anode was operating at a sufficiently high potential or at all potentials when the reactor incorporated a membrane separator. In addition, product collection efficiencies and current density distributions were predicted successfully for a hypothetical photo-anode with 'enhanced' behaviour and are decisive parameters for further geometrical optimisation of the electrodes. The electrolyte flow, light scattering and reflection due to bubbles and other optical phenomena, pH and temperature changes were not included in the present model, which will be extended to include them, based on the works of Singh, 15 Haussener 16 and Holmes-Gentle. 61