En Route to a Unified Model for Photoelectrochemical Reactor Optimization. II–Geometric Optimization of Perforated Photoelectrodes

Results have been reported previously of a model describing the performance of photoelectrochemical reactors, which utilize semiconductor | liquid junctions. This model was developed and verified using SnIV-doped α-Fe2O3 as photoanodes. Hematite films were fully characterized to obtain parameter inputs to a model predicting photocurrent densities. Thus, measured photocurrents were described and validated by the model in terms of measurable quantities. The complete reactor model, developed in COMSOL Multiphysics, accounted for gas evolution and desorption in the system. Hydrogen fluxes, charge yields and gas collection efficiencies in a photoelectrochemical reactor were estimated, revealing a critical need for geometric optimization to minimize H2-O2 product recombination as well as undesirable spatial distributions of current densities and “overpotentials” across the electrodes. Herein, the model was implemented in a 3D geometry and validated using solid and perforated 0.1 × 0.1 m2 planar photoanodes in an up-scaled photoelectrochemical reactor of 2 dm3. The same model was then applied to a set of simulated electrode geometries and electrode configurations to identify the electrode design that would maximize current densities and H2 fluxes. The electrode geometry was modified by introducing circular perforations of different sizes, relative separations and arrangements into an otherwise solid planar sheet for the purpose of providing ionic shortcuts. We report the simulated effects of electrode thickness and the presence or absence of a membrane to separate oxygen and hydrogen gases. In a reactor incorporating a membrane and a photoanode at 1.51 V vs RHE and pH 13.6, an optimized hydrogen flux was predicted for a perforation geometry with a separation-to-diameter ratio of 4.5 ± 0.5; the optimal perforation diameter was 50 µm. For reactors without a membrane, this ratio was 6.5 and 8.5 for a photoanode in a “wired” (monopolar) and “wireless” (photo-bipolar) design, respectively. The results and methodologies presented here will serve as a framework to optimize composite photoelectrodes (semiconductor | membrane | electrolyte), and photoelectrochemical reactors in general, for the production of hydrogen (and oxygen) from water using solar energy.


INTRODUCTION
Thus far, the bulk of the research on photoelectrochemical systems for splitting water with solar energy has focused on material developments (Chen et al., 2010;Tachibana et al., 2012;Moss et al., 2021), aiming to conceive materials that would be: 1) efficient in converting photons to chemical product(s), 2) economical to fabricate and 3) chemically and mechanically durable. These materials are typically synthesized and tested at small scale, with electro-active areas ≤1 cm 2 (Khaselev and Turner, 1998;Rocheleau et al., 1998;Kelly and Gibson, 2006;Jia et al., 2016;Bedoya-Lora et al., 2021). Interest in identifying and resolving the many engineering challenges associated with photoelectrochemical device scale-up has been developing through modelling and experiments at a comparatively slower pace (Carver et al., 2012;Haussener et al., 2012;Lopes et al., 2014;Turan et al., 2016;Hankin et al., 2017;Vilanova et al., 2018;Ahmet et al., 2019;Tembhurne et al., 2019;Vilanova et al., 2020).
Engineering challenges apply to four broad types of photoelectrolysis systems, which have emerged: 1. Photovoltaics (PVs) + electrolysers (Kelly et al., 2011;Jia et al., 2016), which are not thermally or chemically integrated with each other and are connected only electronically. 2. Integrated photoelectrochemical devices (IPECs), in which the electrolysis is powered by a solar cell/PV embedded in the device, but in which the PV is either protected from the electrolyte by an interposed layer to prevent photoelectrode corrosion (Reece et al., 2011;Turan et al., 2016) or is integrated with the electrolyser by some other means, such as thermally (Tembhurne et al., 2019). 3. Photoelectrochemical devices (PECs), with semiconductor | liquid junction(s) in which the semiconductor(s) often also function(s) as catalyst (Brillet et al., 2012;Liu et al., 2013;Shaner et al., 2013;Vilanova et al., 2018). 4. Particle-based photocatalytic water splitting devices (PCWS), which are based on a dispersion of one type of particle, or two particles operating in tandem, connected by a redox mediator (Kuang et al., 2016;Wang et al., 2016Wang et al., , 2020Hisatomi and Domen, 2019).
These systems have sub-classifications based on the finer details of their conceptual design (Nielander et al., 2014). Device designs for systems 1-3 have also been explored and classified schematically to illustrate the wide variety of geometric, electronic and optical configurations that have been conceived to date (Jacobsson et al., 2014;Holmes-Gentle et al., 2018;Moss et al., 2021).
Looking deeper into the designs of systems 2 and 3, which are the focus of the present study, the electrode arrangements may be "wireless" or "wired" bipolar monoliths or "wired" monopolar photo electrodes in terms electric configurations (Newman, 2013;Holmes-Gentle et al., 2018). However, while the focus has been on improving electron/hole transport, the ion transport profiles between electroactive surfaces are also of paramount importance and, if non-uniform, can profoundly affect the reactor performance. Additionally, the usual requirement for "front illumination" (especially for cases when the photoelectrode or PV substrate is not transparent to light or when the charge mobility in the constituent semiconductors are too poor for back-illumination to be appropriate (Eichhorn et al., 2018)) of an embedded PV or semiconductor liquid | junction to maximize the charge separation efficiency can force the hydrogen-and oxygenevolving electro-active surfaces to face away from each other, also leading to non-uniform electric potential distributions and large ionic transport distances (Hankin et al., 2017;Moss et al., 2021). Two-or three-dimensional (c.f. 1-D) electrical potential distributions in the electrolyte may not affect device performance significantly at the millimeter scale that is relevant to most present spontaneous water splitting tests. However, they become critically important on up-scaling (Newman, 2013;Moss et al., 2021).
Regarding perforated electrodes, minimization of inhomogeneities in spatial distributions of current densities between opposing electroactive surfaces in systems 2 and 3 by means of macro-or microscopic perforations has been proposed conceptually (Orazem and Newman, 1984a;1984b) experimentally, and modelled for wired monopolar PECs (Hankin et al., 2017) and wireless bipolar IPECs (Bosserez et al., 2016;Trompoukis et al., 2018;Vijselaar et al., 2019). Example structures of the perforations, as well as their purpose, are illustrated in Figure 1. It has been confirmed through experimental data that electrode perforations offer ionic shortcuts through which current density distributions are minimized, but at the same time they also result in H 2 -O 2 cross-over, thereby decreasing hydrogen collection efficiencies (Bosserez et al., 2016). The task at hand now is to demonstrate using a robust photo-electro-kinetic model exactly what the ideal perforation geometry should be that simultaneously 1) maximizes the homogeneity in the current density distribution across the electrode surface and 2) minimizes losses due to H 2 -O 2 crossover. Such a model should necessarily account for efficiencies of gas evolution from the liquid phase (Vogt, 1984) and effects of diffusion of dissolved gas species within a photoelectrochemical reactor. Moreover, the model should enable the prediction of photo-electro-kinetic current densities as a function of operating conditions, such that the current density distributions over electrode surfaces can be modelled; when photocurrents are not predicted in the reactor models, but for simplicity are assigned fixed values, such distributions cannot be modelled accurately. Some of these aspects were addressed in our 2D model of electrode geometry effects in photoelectrochemical reactors (Bedoya-Lora et al., 2017b), as were the effects of photoanode properties on photocurrent densities (Bedoya-Lora et al., 2017a) and theoretical predictions of photocurrent densities using Sn IVdoped α-Fe 2 O 3 as an example of photoanodes. In this study we focus specifically on the influence of photoelectrode geometry on photoelectrochemical reactor performance.
Modelling of photoelectrochemical reactors involves a wide range of physics: charge, photon, mass, heat and momentum transfer. Electronic charge and photon transfer have been the focus of most of the models, as these processes are inherently related with the performances and efficiencies of photoelectrode materials (Bedoya-Lora et al., 2021). Electron-hole recombination, absorptivity and band gap are material properties related to both processes. The remaining physics, especially mass and momentum transfer in the context of PEC device design, have been studied to a lesser degree, but are most important for up-scaling of these devices.
The geometrical optimization of perforations using appropriate models for PEC and IPEC systems will be of major importance for the design and up-scaling of photoelectrochemical reactors. Differences in ionic pathlengths across photoelectrode surfaces tend to increase with electrode size, depending on geometry, so giving rise to increasingly inhomogeneous electric potential and current density distributions, that also increase with the magnitudes of (mean) current densities and with ionic resistivities. Hence, it has been suggested that large-scale monolithic photoabsorbers are perhaps not practical and that electrodes should be fabricated from Below, we report the performance of a planar and perforated 0.1 × 0.1 m 2 Ti | Sn IV -doped α-Fe 2 O 3 photoanodes (shown in Figure 1B), accompanied by a simulation-based geometrical optimization of different perforation configurations in the absence and presence of a membrane. In addition to coupling gas bubble (c.f. dissolved gas) formation efficiency and product gas (H 2 -O 2 ) cross-over, the objective of this paper was to use our model to optimize several PEC geometries and arrangements, including the effects of photon flux densities and electrolyte conductivities. The model can be used to design composite photoelectrodes, fabricated with µm resolution and, for given conditions and photoelectrode properties, optimize geometries of perforations, which cause loss of photoabsorber area, but homogenize spatial distributions of reaction rates. In principle, perforations can be achieved by using an expanded mesh substrate for photoelectrode deposition (Hankin et al., 2017) or, alternatively, laser ablation (Bosserez et al., 2016) or deep reactive etching (Vijselaar et al., 2019) may be used to introduce perforations into pre-formed (photo-)electrodes.
As discussed previously (Bedoya-Lora et al., 2017b), the effects of bubble generation on modelled photocurrent densities have yet to be implemented, incorporating the consequences of light reflection and scattering (Stevens, 2012), coverage of electro-active area (Vogt, 2011;Hernández et al., 2015), ohmic potential losses in the electrolyte (Vogt and Thonstad, 2017), hydrodynamic (Taqieddin et al., 2017) and local convection (Boissonneau and Byrne, 2000). The complex interactions of these factors on the performance of photoelectrochemical reactors is not yet understood adequately (McKone and Lewis, 2013), so their engineering challenges will be discussed in a future publication.

Geometry Definitions
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 2 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 cationpermeable 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 2 and H 2 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 ) 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.

FIGURE 3 | (A)
Square and hexagonal distributions of perforations in the photo-anode, with lengths and areas defined. Electrode designs for geometrical optimisation: (B) "wired" (photo-anode facing away from cathode) and (C) "wireless" (monolithic photo-bipole) with a single photo-anode.

Target Definitions
As defined in part I, system performances were evaluated in terms of four parameters: current density (j i ), hydrogen flux (N i ), charge yield (Φ e i ) and gas collection efficiency (Φ collection,i ) of species i.
EXPERIMENTAL Photoanode Fabrication 0.1 × 0.1 m 2 Sn IV -doped α-Fe 2 O 3 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 Reactor
Photoelectrochemical measurements were carried out in a photoelectrochemical reactor designed for operation with 0.1 × 0.1 m 2 electrodes. The reactor comprised two chambers, each with a volume ≈1 dm 3 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 −1 , and a twoelectrode 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 −2 was achieved over an area of 0.12 × 0.12 m 2 and an estimated 497 ± 12 W m −2 over 0.1 × 0.1 m 2 , 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 −2 over the photoanode surface.
Photoelectrochemical properties of un-annealed Ti | Sn IV -doped α-Fe 2 O 3 have been measured previously using a small photoelectrochemical cell (0.06 dm 3 ) and an electroactive area of 3 × 6 mm 2 . 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.

Software
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 −3 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 −7 m and maximum size of 8 × 10 −4 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.

PRELIMINARY CONSIDERATIONS
Effects of the Temperature on the Performance of the Cell Efficiencies of hydrogen bubble formation, f g,H2 , 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 2 than for H 2 (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,O2 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 1 and T 2 : 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,O2 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 2 O 3 ) 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.

Photon Flux Distribution of Solar Simulator
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 2 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 , 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.

3D Up-Scaled Reactor Model
The predicted photocurrent densities as a function of potential in the up-scaled reactor using 0.1 × 0.1 m 2 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 2 and O 2 fluxes were not stable and produced nonstoichiometric 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 2 :O 2 , oscillated between 2.3 and 2.0, in agreement with electron stoichiometries of the cathodic and anodic reactions (Supplementary Figure S9); values of H 2 :O 2 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 2 and H 2 collection efficiencies

3D (Symmetrical) Model Optimization
An un-annealed Ti | Sn IV -doped α-Fe 2 O 3 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 FIGURE 5 | Effect of inter-perforation distance (l d ) on current densities for "wired" reactor: Ti | Pt | 1 M NaOH | Sn IV -doped α-Fe 2 O 3 | Ti, (A) with and (B) without a membrane; photoanode at 1.51 V vs RHE with a perforation size of l p 50 μm, l d /l p 2 to 10; and concentration profiles of oxygen and hydrogen for the same 'wired' reactor (C) with and (D) without a membrane in the z-direction in the center of the perforation (x 0, y 0) for l d /l p 2 and 10.
Frontiers in Chemical Engineering | www.frontiersin.org September 2021 | Volume 3 | Article 749058 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 −5 mol m −2 s −1 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 2 -O 2 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 2 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 2 -O 2 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 + 2H 2 O + 4e − → 4OH − . 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 photoabsorber and Pt and IrO 2 as catalyst films for H 2 and O 2 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 FIGURE 7 | Effect of inter-perforation distance (l d ) on current densities for "wireless" (photo-bipolar) reactor (A) with and (B) without a membrane, with 1 M NaOH | Pt | Ti | Sn IV -doped α-Fe 2 O 3 | 1 M NaOH, photoanode at 1.51 V vs RHE with a perforation size of l p 50 μm, l d /l p 2 to 10; and concentration profiles of oxygen and hydrogen for the same "wireless" reactor (C) with and (D) without a membrane in the z-direction (i.e. through catholyte | membrane (z 0.02 m)| anolyte) in the center of the perforation (x 0, y 0) for l d /l p 2 and 10.
Frontiers in Chemical Engineering | www.frontiersin.org September 2021 | Volume 3 | Article 749058 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 2 -O 2 crossover 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 2 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.

Sensitivity Analysis
The effects on hydrogen fluxes were predicted by varying the following parameters: semiconductor donor density (10 24 -10 27 m −3 ), irradiance (10 −1 -10 2 suns), electrode thickness (0.1-3 mm) and electrolyte conductivity (10 −1 -10 2 S m −1 ). 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 −2 ) 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.

CONCLUSION
A model was developed in terms of measurable properties of the photoelectrode (donor density, photon absorption, bulk and interfacial electron-hole efficiencies, flat band, thickness, exchange current density and Tafel slope), spectrally resolved photon flux of the light source, volumetric mass transfer for gas desorption, bubble formation of hydrogen and oxygen at the electrode surface, and the presence or absence of a membrane between electrodes. The model was validated successfully against experimental data obtained from an upscaled reactor using a 0.1 × 0.1 m 2 photoanode, Sn IV -doped α-Fe 2 O 3 spray pyrolysed on titanium. Different reactor designs, parameter values and geometries were studied and optimized as a function of current densities, hydrogen fluxes and collection efficiencies, accounting for H 2 -O 2 cross-over, ohmic potential losses and current density distributions. Optimal ratios of separation to size of perforations were found for a wide range of geometrical configurations and photoanode properties. The effects of temperature were considered to a limited extent, but a deeper understanding of its effects on the photoanode performance should be addressed, as well as the inclusion of convection effects (or flow) on bubble efficiency and H 2 -O 2 crossover rates.
The model predictions emphasized the importance of using a membrane for quiescent electrolytes, especially when a "wireless" photo-bipolar reactor design is used. Hydrogen fluxes for "wired" and "wireless" reactor configurations were comparable when a membrane was incorporated. By contrast, in the absence of a membrane, the model predicted significant H 2 and O 2 cross-over for a 'wired' system, while "wireless" systems exhibited low charge yields for the range of parameters evaluated in the photoelectrochemical model.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
FB made substantial contributions to the conception of the work, to the acquisition, analysis and interpretation of data for the work, as well as writing the manuscript. AH contributed significantly to the model formulation and to the production of the manuscript. GK was the project leader.