The dimensions and geometry of an existing up-scaled reactor (Hankin et al., 2017) was replicated in a 3D model to predict photocurrent densities and hydrogen fluxes and to enable their comparison with values determined experimentally. All flux densities and current densities were normalized by the geometrical area of the photoanode. A planar photoanode with dimensions 0.1 × 0.1 m² was exposed to illumination with (photo-electro-) active surface facing away from the cathode, as shown in Figure 2. Hydrogen and oxygen gases were collected in two chambers separated by a cation-permeable Nafion® membrane. This reactor geometry, with illumination of the front of the photoelectrode through the electrolyte, was selected following previous results with better performance than its counterparts with back-illumination and front-illumination through the mesh counter electrode. All the nomenclature referred to in this manuscript is listed in Table 1.

A symmetrical 3D model was developed to optimize the photoelectrode geometry with the aim of reducing losses due to ohmic potential drops and cross-over of O₂ and H₂ products, via perforations and/or membrane. As shown in Figure 3, square and hexagonal distributions of circular perforations were introduced in the planar photoanode structure. The ratio between the separation distance (l _d) and diameter (l _p) of the perforations, l _d/l _p, was optimized for several perforation sizes and electrode thicknesses, l _e. As mentioned previously, to enable rational comparison between the performances with different configurations, all the fluxes and currents densities were calculated with respect to the geometrical area (including perforated area) and not only the electro-active area, as shown in Figure 3A. Furthermore, two electrode arrangements were studied as shown in Figures 3B,C: spatially separated monopolar electrodes (‘wired’) with the photoanode facing away from the cathode and a monolithic photo-bipole (“wireless”) (Newman, 2013). The model with “wired,” front-illuminated photoanodes, Figure 3B, was selected due to its similarities with existing photoelectrochemical cells for material characterisation and scaled-up reactors (Figure 2).

The distance between electrodes for the wired configuration, l _a-c, was kept constant at 0.01 m. Other possible configurations, e.g. photoelectrodes separated by perpendicular separators (Haussener et al., 2012) and integrated light-absorber between electrodes (Berger and Newman, 2014) have been investigated. The “wireless” configuration, Figures 3C is possibly the most promising and studied arrangement to date. Several materials and models have been reported using this configuration (Reece et al., 2011; Jin et al., 2014; Bosserez et al., 2016; Vijselaar et al., 2019), although up-scaled reactors have yet to be developed.

As defined in part I, system performances were evaluated in terms of four parameters: current density ( j i ), hydrogen flux ( N i ), charge yield ( Φ i e ) and gas collection efficiency ( Φ collection , i ) of species i. N i = N g , i + N desorption , i (1) Φ i e = j i j total (2) Φ collection , i = N i j total / ν e . i F (3)

0.1 × 0.1 m² Sn^IV-doped α-Fe₂O₃ photoanodes were deposited on titanium by spray pyrolysis, following the procedure reported previously (Bedoya-Lora et al., 2016; Hankin et al., 2017). The electrodes were produced on planar Ti sheets and also on perforated Ti sheets, with rectangular perforations of ca. 3 × 1.5 mm, spaced 5.8 mm horizontally and 2 mm vertically apart. In the present case, photoanodes were not annealed after the deposition of hematite films.

Photoelectrochemical measurements were carried out in a photoelectrochemical reactor designed for operation with 0.1 × 0.1 m² electrodes. The reactor comprised two chambers, each with a volume ≈1 dm³ and separated by a cation-permeable membrane (PTFE-reinforced Nafion® 424, DuPont Inc.), selected for its chemical stability, durability and superior mechanical properties over an anion-permeable counterpart, which otherwise would have been more appropriate for alkaline water splitting. Both compartments were filled with 1 M NaOH (pH 13.6). A potentiostat/galvanostat (Autolab PGSTAT 30) was used to control the reactor with a three-electrode configuration for voltammetry on the photoanode at 50 mV s⁻¹, and a two-electrode configuration for chronoamperometry at an applied cell potential difference of 1.6 V (resulting in a photoanode potential of ≈1.51 V vs RHE) for 24 h. The reason for the two-electrode configuration in the latter case is that the use of a reference electrode did not permit adequate sealing of the reactor for gas evolution rate measurements to be made. Hematite photoanodes acted as working electrodes, platinized titanium mesh (Expanded Metal Company, United Kingdom) as counter electrode and saturated AgCl|Ag as reference electrode (1.001 V vs RHE and 0.197 V vs SHE). The reactor was operated in batch mode without recirculation of electrolyte solution, the hydrodynamic effects of which have yet to be implemented in the model.

A solar simulator (Sun 2000; Abet technologies, United States) with a 550 W Xe arc lamp was used to irradiate the photoanode. The light source was calibrated and mapped at the reactor working distance using a UV-Vis spectrophotometer coupled to a CR2 cosine receptor (Black-Comet CXR-25, StellarNet, United States). An average power density of 433 ± 82 W m⁻² was achieved over an area of 0.12 × 0.12 m² and an estimated 497 ± 12 W m⁻² over 0.1 × 0.1 m², which corresponds to the photoanode area, as shown in Supplementary Figure S1 and S2 in the supporting documentation. An unavoidable degree of non-uniformity in the spatial distribution of the irradiance was evident. The intensity of incoming light was corrected for attenuation by the electrolyte and quartz to give an effective power density of 461 W m⁻² over the photoanode surface.

Photoelectrochemical properties of un-annealed Ti | Sn^IV-doped α-Fe₂O₃ have been measured previously using a small photoelectrochemical cell (0.06 dm³) and an electroactive area of 3 × 6 mm². As discussed in the next section, parameters from experimental results reported previously (Bedoya-Lora et al., 2017b) were also employed in the model described below after correction for any changes in the power density delivered by the light source. The average absorbed photon flux for the photoanode is presented in Supplementary Figure S3 in the supporting information.

Voltammograms were obtained over a sufficiently short time period for electrolyte temperatures to remain constant at 25 C. As shown in Supplementary Figure S4 in the supporting documentation, prolonged illumination caused electrolyte temperatures to rise to ca. 40 C within the first ca. 6 h, after which thermal equilibrium with surroundings was achieved. Hence, chronoamperometry results are not shown for this initial 6-h period. Gas flow rates for cathode (hydrogen) and anode (oxygen) compartments were recorded by gas flow meters (MilliGascounter MGC-1, Ritter, Germany) during longer term (24 h) chronoamperometry; the steady state temperature was used to correct molar flux densities, assuming ideal gas behavior.

Comsol Multiphysics 5.2a with a Batteries and Fuel Cells module was used to solve the system of equations 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 bulk of the electrolyte. A stationary non-linear solver was used, typically with a relative tolerance of 0.001 and a maximum of 50 iterations; the suitability of these parameters has been reported previously (Hankin et al., 2017).

The mesh of the up-scaled reactor was geometry-dependent with a minimum element size of 3.33 × 10⁻³ m and maximum size of 0.02 m. The mesh was optimized for higher resolutions in narrow regions. The model converged typically in under 3 min for a given set of conditions. In the case of models for geometric optimization, a minimum element size of 8 × 10⁻⁷ m and maximum size of 8 × 10⁻⁴ m was used. The model converged usually in under 10 s. A parametric sweep was used to evaluate a wide range of settings and geometries.

Efficiencies of hydrogen bubble formation, f g , H 2 , do not depend significantly on the temperature, whereas oxygen bubble formation efficiencies increase rapidly with temperature. This could be a result of the effect of temperature on saturated concentration of dissolved gas, which is stronger for O₂ than for H₂ (Shoor, 1968). Conforming with the extensive data reported by Chin Kwie et al. (2020) and Bedoya-Lora et al. (2017b) and given the lack of empirical correlations to calculate efficiencies of oxygen bubble formation at higher temperatures (>25 C), it was assumed that during chronoamperometric operation at thermal equilibrium (≈40 C under our experimental conditions), f g , O 2 was approximately double that calculated for 25 C. Bubble formation coefficients for hydrogen gas evolution remained unchanged (Joe et al., 1988; Bedoya-Lora et al., 2017b). Diffusion coefficients were corrected for higher temperatures using the (ratio of the) Stokes-Einstein equation for temperatures T ₁ and T ₂: D 1 D 2 = T 1 T 2 μ 2 μ 1 (4)

Concentrations of dissolved oxygen and hydrogen at saturation were interpolated from published experimental data (Shoor, 1968) for KOH solutions at different temperatures and concentrations. Values estimated at 25 C using these experimental data were remarkably similar to more recent values determined for this temperature (Davis et al., 1967; Joe et al., 1988; Tromans, 1998), so it was assumed that interpolations at 40 C and 1 M KOH were sufficiently reliable to be used in the model. In contrast to oxygen, it was found that hydrogen solubility decreased only slightly with temperature, offering an explanation of why f g , H 2 do not increase significantly with temperature, while f g , O 2 is double at 40 C (Joe et al., 1988), Furthermore, this is in agreement with Vogt’s model, in which the saturation concentration plays an important role in the efficiency of gas evolution (Vogt, 1984). Due to changes in diffusion coefficients and kinematic viscosity of water with temperature, volumetric mass transfer coefficients, k L a , were recalculated as a function of these parameters. A numerical description of these corrections can be found in Supplementary Tables S1–S3,S6. (Frossling, 1938; Marrucci and Nicodemo, 1967; Joshi and Sharma, 1979; Deckwer, 1992; Wüest et al., 1992; McGinnis and Little, 2002; Painmanakul et al., 2009).

The effects of the change in temperature, from 25 to 40 C, on the properties of the semiconductor (Sn^IV-doped α-Fe₂O₃) were considered too small to merit consideration in this particular study. Also, the effect of temperature on membrane properties were disregarded due to the lack of data available for the conditions and materials assessed in this work.

As discussed and demonstrated earlier, the photon flux delivered by the solar simulator was not radially uniform across the photoanode surface. This non-uniformity was a potential concern, as it can cause a significant distribution in “overpotentials” and photocurrent densities across a 0.1 × 0.1 m² electrode. The photon flux density was evaluated as a function of position on the photoanode surface, I 0 = f ( x , y ) (Supplementary Figure S5). Measurements were corrected for attenuation by electrolyte and the quartz window, and were fitted by regression analysis into a spatially rotated exponential function that was subsequently fed into the Comsol Multiphysics model: I 0 =   D 1 +   D 2 ⁡ exp [ − ( x ⁡ cos π 4 − y ⁡ sin π 4 ) 2 ]   +   D 3 ⁡ exp [ − ( x ⁡ sin π 4 − y ⁡ cos π 4 ) 2 ] (5)

I 0 , I 0 − I x and ( I 0 − I x ) α functions were fed into the model and used to estimate the spatial distribution of photons absorbed locally by the photoanode. A more detailed numerical description of these functions can be found in Supplementary Table S4 and Supplementary Figure S5.

The predicted photocurrent densities as a function of potential in the up-scaled reactor using 0.1 × 0.1 m² planar and perforated electrodes were within the error of experimental measurement (Supplementary Figure S6), as confirmed in an earlier version of the model (Hankin et al., 2017).

Hydrogen and oxygen fluxes were measured and estimated after operating the reactor at 1.6 V cell potential difference, which corresponded to ca. +1.51 V vs RHE at the anode. During the first ca. 6 h, H₂ and O₂ fluxes were not stable and produced non-stoichiometric yields due to a combination of heating of the system and gradual saturation of catholyte and anolyte with gases. The generated volumes of both gases were recorded over time during these experiments and are reported in Supplementary Figure S7. Reactor performance was relatively stable after 10 h of operation. As shown in Supplementary Figure S8, predicted current densities and fluxes of hydrogen and oxygen agreed with experimental data after 10 h of photoelectrolysis. The molar ratio between evolving gases, H₂:O₂, oscillated between 2.3 and 2.0, in agreement with electron stoichiometries of the cathodic and anodic reactions (Supplementary Figure S9); values of H₂:O₂ above and below 2.0 may have been caused by oxygen leaking from the anolyte, due to imperfect sealing around the membrane and quartz window, and lower precision of the gas flow meters in measuring slower rates of oxygen evolution compared to hydrogen evolution. Charge yields for the planar electrode were ca. 0.9 and 0.8 for hydrogen and oxygen evolution, respectively. In the case of the perforated photo anode, charge yields oscillated between 0.9 and 1.0 for both oxygen and hydrogen. However, the predicted value for this parameter according to the model should have been close to unity (0.99). Again, the deviation of experimental data from model predictions was probably associated with non-steady state operation and leaking of gases during prolonged experiments. However, the agreement of the model with experimental results was still very good.

Current densities at the center of a planar photo anode (Supplementary Figure S10) were higher than at the edges due to the non-uniform illumination associated with the light source (Supplementary Figure S5). However, the predicted “overpotential” ( Δ ϕ SC = U applied − U fb ) followed the opposite trend and was higher at the edges of the electrode. The mismatch of these spatial distributions was attributed to the dependence of photocurrent densities on two distinct variables: the incident photon flux, I 0 , a function of the light source and the “overpotential,” Δ ϕ SC , affected by the potential distribution between anode and cathode. In this case, the distribution of the photon flux overrode the effect of the “overpotential” distribution over the photoanode. Hence, distributions of current densities were not as pronounced due to the relatively poor performance of the hematite photoanode used to exemplify the effects. In contrast, the predicted “overpotential” for a perforated photoanode was more homogeneous and slightly higher current densities are achieved in the center of the electrode surface compared to a planar photoanode (Supplementary Figure S10). These findings highlight the importance of accounting for the various unavoidable imperfections in the synthesis of materials and conditions in up-scaled reactors in order to draw correct conclusions from the experimental data and the model. More heterogeneous distributions were predicted for photoelectrodes with notionally enhanced catalytic properties for oxygen evolution (Bedoya-Lora et al., 2017b).

According to model predictions, steady state oxygen concentrations reached supersaturation in the anolyte chamber adjacent to the photoanode surface, with analogous behaviour for hydrogen in the catholyte. Consequently, the flux of dissolved gases through the compartments was small and was minimized effectively when a membrane was present (Supplementary Figure S11). This also implies that significant gas desorption rates were operative only near electrode surfaces (<10 mm), as shown in Supplementary Figure S12. Although small fluxes of dissolved gas were predicted to occur through the perforations, even in the absence of a membrane, their magnitude was not sufficient to affect model predictions, which corresponded to O₂ and H₂ collection efficiencies of ca. 0.99. In contrast, experimentally determined values of 0.8 suggested potential leakage of gases through the membrane seals (nitrile rubber) of the reactor.

An un-annealed Ti | Sn^IV-doped α-Fe₂O₃ photo-anode in 1 M NaOH at 25 C and at a potential of 1.51 V vs RHE was modelled under illumination by 1.5AM light; properties are presented in Supplementary Tables S5, S7 and S8 (Frossling, 1938; Powell and Tye, 1961; Marrucci and Nicodemo, 1967; Ogumi et al., 1985; Joe et al., 1988; Sone et al., 1996; Haug and White, 2000; Painmanakul et al., 2009; Sheng et al., 2010; ASTM G173-03, 2012; Lee, 2013). A Ti | Pt cathode was used in 1 M NaOH electrolyte solution. Nafion® 117 was used in some cases as membrane and for filling the perforations in the photoanode throughout its thickness, which was assumed initially as 500 µm. As shown schematically in Figure 3, this model was used to predict values for the ratio between perforation separation/size, l d / l p , that correspond to current density and hydrogen flux maxima. As reported previously (Bedoya-Lora et al., 2017b), the optimum l d / l p values for hydrogen fluxes were affected by the collection efficiency, which was lower for shorter photoanodes and reached almost unity when the distance (l _d) between perforations was increased. On the other hand, inhomogeneous current density distributions were exacerbated for higher values of l d / l p . Furthermore, the total current was also affected by the inactive area due to perforations, so for very low values of l d / l p , low current densities and hence low rates of hydrogen evolution were achieved.

Hexagonal vs square perforation distributions: the effects of spatial distributions of perforations on current densities and hydrogen fluxes for a “wired” reactor configuration without the perforations being filled with a membrane, are presented in Supplementary Figures S13 and S14. There was no difference between squared and hexagonal distributions, when current densities and hydrogen fluxes were plotted as functions of the ratio between the active and geometrical area. When plotted as a function of l d / l p , there was a noticeable shift towards higher values of current densities and hydrogen fluxes for the hexagonal geometry. However, the optimum value was very similar for both geometries, with small changes in a range of l d / l p = 5.0 ± 1 for maximum current density and l d / l p = 6.5 ± 2 for maximum hydrogen fluxes. Evidently, current densities may not be the key performance indictor to target for optimization, since higher current densities do not translate into higher collection efficiencies, particularly in the absence of a membrane.

For l p < 50 μm, there was no significant increase in the resulting current densities and hydrogen fluxes, as shown in Supplementary Figure S14B. Hence, 50 µm was selected as an appropriate perforation size in a hexagonal geometry for the simulations reported below.

Optimized l d / l p   for a ‘wired’ reactor configuration: Figure 4A compares the model reactor’s performance with and without a membrane. As shown in Figure 4B, when a membrane was used, collection efficiencies were almost unity for all geometries, so that the optimum point for current densities and hydrogen fluxes occurred at the same l d / l p = 4.5 ± 0.5 . Figure 4C shows the maximum current densities and hydrogen fluxes for different perforation sizes with electrodes separated by a membrane; both outputs were lower for larger perforations. This figure also confirms that perforations with diameters in the range of 200 to 50 µm did not affect the photoanode’s performance significantly at optimised conditions. When a membrane was used to separate the hydrogen and oxygen, the optimum geometries for hydrogen evolution rates and current densities were essentially the same, due to the collection efficiencies having been maximised. However, current densities were higher in the absence of a membrane, due to decreased ohmic potential losses across the perforation, dependent on its thickness.

In contrast, Figure 4A also shows that optimum hydrogen fluxes reached similar values of ca. 2.45 × 10⁻⁵ mol m⁻² s⁻¹ with or without a membrane at different l d / l p ratios. Hence, the decision of whether to incorporate a membrane rests on the extent of H₂-O₂ cross-over between anolyte and catholyte. As shown in Figure 4B, in the absence of a membrane and for l d / l p = 6.5 , 3 mol% of oxygen was predicted in the H₂ gas collected from the catholyte, increasing to 14% for l d / l p = 2 . As the explosive limit of oxygen in hydrogen is ca. 4%, membranes are recommended for safety reasons.

Figure 5 shows that higher current densities can be achieved without using a membrane, but a higher concentration of oxygen (above saturation concentration) was found in the catholyte. The membrane effectively separated hydrogen and oxygen gases, and concentrations were almost identical for low and higher values of l d / l p . However, current densities for l d / l p = 10 decreased by almost 10% compared to the optimum.

Optimized l d / l p for a “wireless” configuration: Predicted current densities and hydrogen fluxes differed greatly in the presence and absence of a membrane filling the perforations in a “wireless” system. In absence of a membrane, current densities and hydrogen fluxes were at their maximum when l d / l p = 2.5 and 8.5 , respectively, as shown in Figure 6A. Again, this was explained by the increase of H₂-O₂ cross-over rates in the absence of a membrane. However, the amount of oxygen collected in the catholyte was 2 orders of magnitude lower (<0.18% molar) when compared to a “wired” configuration (<14% molar). In a wireless design, oxygen reached the cathode surface faster and most of it was reduced near the perforation, while in a “wired” reactor, oxygen had to diffuse through the membrane, into the catholyte and then to the cathode surface, at which it was reduced under transport control according to the reaction: O 2 + 2 H 2 O + 4 e − → 4 OH − .

This was confirmed by predictions shown in Figure 7, in which concentrations of oxygen and hydrogen in the catholyte and anolyte, respectively, were near zero. Despite the lack of oxygen found in the gas collected from the cathode compartment, collection efficiencies of a membrane-less and “wireless” system decreased greatly for low values of l d / l p , as shown in Figure 6B. Consequently, hydrogen fluxes were 7% lower when compared to the system with perforations filled with ion-permeable membrane material. In the latter case, current densities and hydrogen fluxes were again found to be optimal at l d / l p = 4.5 . Hydrogen fluxes were greater than for “wired” reactors, but comparable current densities were predicted for “wired” and “wireless” reactors. It has been predicted that up-scaled “wired” reactors will exhibit superior performance, but the effects of perforations were not considered and hydrogen fluxes were not reported (Newman, 2013).

The optimum geometry did not change significantly on increasing perforation sizes (l _p), but the total current density and net hydrogen fluxes decreased, due to loss of photoactive area, as implied by data in Figure 6C; hence, perforations <200 µm are preferred. Bosserez et al. performed a similar experimental study using perforated monoliths (“wireless” configuration without a membrane, using n-Si as photo-absorber and Pt and IrO₂ as catalyst films for H₂ and O₂ evolution, respectively (Bosserez et al., 2016)) aiming to minimize ohmic potential losses using different sizes and distances between perforations, for 3.2 < l d / l p < 4.3 . It was concluded that for 1 M KOH and ohmic potential losses <100 mV, l d should be <1,000 μm, which corresponds to l d / l p < 4.3 . These results are consistent with our predicted current densities presented in Figure 6A for the same system (“wireless” photo-bipolar reactor without a membrane). However, higher current densities were achieved at the cost of decreasing hydrogen collection efficiencies, as suggested by Figure 6B, in which Φ collection < 0.63 for l d / l p < 4.3 . Experimental collection efficiencies were also estimated by Bosserez as Φ collection ≈ 0.78 ± 0.1 (Bosserez et al., 2016), which agrees with our predictions ( Φ collection = 0.76 ) using a modified model geometry that resembles Bosserez’s system ( l e = 550   μm , l p = 232   μm , l d / l p = 4.3 and j = 79   A m − 2 ). Our results are also in broad agreement with those of (Vijselaar et al., 2019).

As shown in Figure 7, in the presence of a membrane, there was no difference between the hydrogen and oxygen concentration profiles at low and high l d / l p values. In the absence of a membrane, concentration profiles were higher for l d / l p = 10 due to an increased diffusion resistance through the perforation. Nevertheless, concentrations of oxygen and hydrogen in the catholyte and anolyte, respectively, were significantly lower than for a membrane-less “wired” system, compare Figure 5D and Figure 7D, implying that H₂-O₂ cross-over rates were minimized, but charge yields decreased. However, Bosserez et al. (Bosserez et al., 2016) reported experimental charge yields for hydrogen of >0.97 in all cases, because on their IrO₂ anode, hydrogen oxidation would have been kinetically controlled, while in our model, it was assumed to be controlled by mass transport (as very reasonably was oxygen reduction at the Ti/Pt cathode).

Figures 7A,B also show that current densities for different l d and a perforation size of l p = 50   μm were only slightly higher than those for a “wired” system (Figure 5A,B, with a small difference at the optimum l d / l p = 4.5 . The reasons for differences between the predictions for a “wireless” system with and without a membrane are the same as discussed in the section above for a “wired” system.

The effects on hydrogen fluxes were predicted by varying the following parameters: semiconductor donor density (10²⁴–10²⁷ m⁻³), irradiance (10⁻¹–10² suns), electrode thickness (0.1–3 mm) and electrolyte conductivity (10⁻¹–10² S m⁻¹). For each range of variables, the relative optimum hydrogen fluxes were defined as a function of l d / l p . Figure 8 and Figure 9 illustrate these predictions for “wired” and “wireless” reactors, respectively, with perforations filled with ion-permeable membrane material. Dashed lines represent the typical values used in the model (also listed in Supplementary Tables S7 and S8); under these conditions an optimum hydrogen flux was reached at ca. l d / l p =   4.5 . These figures also show the significant effects of electrode thickness and irradiance on the optimum configuration.

The effect of electrode thickness, increasing ohmic potential losses, was expected and comparable to results reported previously (Bosserez et al., 2016). When the perforation diameter/area was smaller while maintaining the same electrode thickness, l e / A perforation , increased as did the ohmic resistances. Hence, for thicker electrodes, higher perforation areas should be used (lower l d / l p ). When concentration radiation (>1000 W m⁻²) is to be used, it is especially important to define an optimum l d / l p , as hydrogen fluxes were more highly dependent on photon fluxes reaching the photoanode, compared to the effects of the other three parameters.

On the other hand, optima for donor density and electrolyte conductivity were near the same l d / l p ≈   5 ± 1 for a wide range of values for those properties. However, Figure 8 and Figure 9 report hydrogen fluxes relative to that at the local maximum; current densities and hydrogen fluxes were affected greatly by changes in donor densities and, to a limited extent, by conductivities of electrolyte solutions, as reported in Supplementary Figures S15 and S16.