Integrated Photo-Electrochemical Solar Fuel Generators under Concentrated Irradiation I. 2-D Non-Isothermal Multi-Physics Modeling

We investigate the direct conversion of solar energy and water into a storable fuel via integrated photo-electrochemical (IPEC) devices. Here we focus particularly on a device design which uses concentrated solar irradiation to reduce the use of rare and expensive components, such as light absorbers and catalysts. We present a 2-dimensional coupled multi-physics model using finite element and finite volume methods to predict the performance of the IPEC device. Our model accounts for charge generation and transport in the photoabsorber, charge transport in the membrane-separated catalysts, electrochemical reaction at the catalytic sites, fluid flow and species transport in the porous charge collectors and channels, and radiation absorption and heat transfer for all components. We then develop performance optimization strategies utilizing device design, component and material choice, and adaptation of operational conditions. Our model predicts that operation under high irradiation is possible and that dedicated thermal management can ensure high performant operation. The model shows to be a valuable tool for the design of IPEC devices under concentrated irradiation at elevated temperatures. To our knowledge, it is the most detailed yet computationally low-cost model of an IPEC device reported.

The solar energy received on earth's surface can meet mankind's current and future energy demand. 1,2 Direct conversion of solar energy and water into chemical energy via photo-electrochemical (PEC) processes is one viable route for renewable fuel processing and energy storage. Integrated photo-electrochemical (IPEC) devices, i.e. composed of an integrated buried photovoltaic (PV) component and an integrated electrochemical component (consisting of membraneseparated catalysts and porous charge collectors), allow circumvention of some of the challenges imposed by solid-liquid interfaces in traditional PEC devices, while also exhibiting potential to operate at higher efficiencies and lower cost than externally wired (non-integrated) PV plus electrolyzer (EC) devices. [3][4][5] We refer to the design as "integrated" to signify that the PV and electrolyzer components are area-matched and in direct contact, allowing heat transfer from one component to the other and for thermal management strategies to be applied. In order to increase the economic competitiveness of IPEC devices compared to conventional hydrogen generation pathways, we consider concentration of irradiation. [6][7][8] This leads to large driving current densities (approximately proportional to the concentration factor), and thus introduces larger overpotentials and potential mass transport limitations. 9 Concentration also decreases the performance of the photoactive components due to increased temperature. On the other hand, the kinetics are enhanced with increased temperature. Ionic transport in the solid electrolyte is also enhanced with increased temperature, but this increase stops abruptly and sharply drops at temperatures above 120 • C due to membrane dry out. 3,10,11 This competing and coupled behavior of the components requires a detailed understanding of the heat transfer, charge transport, fluid flow, and reaction kinetics in order to formulate performance optimization strategies for concentrated integrated photo-electrochemical (CIPEC) devices via device design and adaptation of the operational conditions.
Multi-physics computational models are a crucial support in device design and engineering. They allow in-depth analysis of conceptual designs, support feasibility investigations of devices and integrated systems, and permit the quantification of performance. Modeling efforts of PEC and specifically IPEC devices are limited. The earliest attempts used lumped-circuit models of a photocell in series with a current-dependent electrochemical load. 12 Berger et al. 13 presented a basic 1-dimensional model for light absorbers and electrolysis with applicability to both wired and wireless PEC systems. Gaudy et al. 14 extended the 1D model by adding detailed wave propagation modeling and an advanced semiconductor-electrolyte interface model accounting for pinning and unpinning of interface states. Haussener et al. 3,15 developed a 2-dimensional PEC model focusing on the charge transport in the electrolyte and reaction kinetics incorporating the idealized Shockley-Queisser limit for the photoabsorber approximation. In Haussener et al., 3 an isothermal model was used to provide predictions for the temperature-dependent performance behavior of PEC devices.
However, none of the modeling efforts accounted for the detailed solution of the energy conservation and the resulting corresponding spatial variations in temperature. Understanding these variations assists the development of thermal management strategies to benefit overall performance, i.e. guiding operation and design to maximize high-temperature advantages in kinetics and transport while minimizing high-temperature disadvantages in charge transport and recombination in the photoabsorber. Thermal management is therefore a rationale for CIPEC device designs. The formulation of detailed guidelines not only requires the solution of the energy conservation in a multi-dimensional design, but also the detailed modeling of charge transport in the semiconductor materials in order to reliably predict charge separation and recombination. Furthermore, the detailed prediction of the charge generation (the radiation absorption) in a nanostructured heterogeneous component is needed, requiring the solution of Maxwell equations instead of the usually applied Beer's law. None of the modelling efforts, up till now, have accounted for all of these phenomena, namely a coupled 2-dimensional multi-physics model for light propagation and absorption, semiconductor physics, fluid flow and reaction kinetics, and energy conservation. Such advanced multiphysics, multi-dimensional models require focusing on the accurate definition of the boundary conditions between the components and the consistency of the physical phenomena. The various conservation and transport equations must be solved with accurate interface conditions for component coupling. This coupling introduces additional complexity, as detailed component models accounting for a subset of physical phenomena rely on the solution of other subsets of equations. Consequently, such coupled modeling efforts require additional external iterative solution steps, which generally increases computational efforts and features challenges related to model robustness. We address these complexities in this study.
We have developed an advanced multi-physics and multi-scale numerical modeling tool to assist in designing and building an integrated (a) 3D schematic (not to scale) of the integrated PEC device depicting the incoming concentrated irradiation, cooling and preheating water channel, triple-junction solar cell (violet colors), and integrated electrolyzer consisting of anodic and cathodic channels, gas diffusion layers (GDL), catalyst layers, and polymeric electrolyte (Nafion). The 2D simulation domain in the xy-plane is indicated by the dashed orange rectangle. Electromagnetic wave propagation (EM), semiconductor charge transport (SC), heat transfer (HT), fluid flow and reacting fluid flow (FF, RFF), and electrochemical charge transport (EC) simulation domains are represented by their respective color labels in the schematic. (b) The energy re-partition diagram shows the distribution of energy at various stages of its utilization for C = 18. The width of the arrows shows approximate re-partitioning for the device utilizing a III-V based solar cell. The darkened areas of the arrows of the PV and EC losses are not considered in the total heat source of the energy equation, as discussed in section Heat transfer.
photo-electrochemical device using concentrated solar irradiation. The CIPEC device is shown in Fig. 1a. A variation of this design presented in Fig. S.2 integrates a novel self-tracking concentrator 16 and provides additional design considerations specific to that design. Concentrated solar irradiation irradiates the device from the x-direction where it is delivered and absorbed by a buried dual/triple-junction (e.g. Ga 0.51 In 0.49 P-GaAs or aSi-ucSi-ucSi) PV cell. The radiation arriving at the PV cell produces electron-hole-pairs (EHPs) if the radiation energy exceeds the bandgap energy of the absorber materials and is effectively absorbed. The charge carriers are then delivered to the electrochemical components, i.e. the solid electrolyte and the catalysts, building the integrated electrolyzer. The arrangement of the PV and EC can be done either in p-n-cathodic-anodic configuration, shown in Fig. 1a, or n-p-anodic-cathodic configuration. The holes are delivered to the anode causing oxidation of water and production of oxygen and protons at the catalytic sites. The protons travel to the cathode through the polymeric electrolyte where they are reduced by the electrons delivered from the PV's n-terminal to produce hydrogen at the catalytic sites. A water channel between the concentrator exit and the PV cell is introduced to cool the PV cell as well as to preheat the reactant (water) before it enters the anodic electrolyzer channels. The connection between the water channel and the anodic chamber is represented by black dots in the schematic of Fig. 1a. The reactant is therefore preheated by the energy which is rejected from the photoabsorbers.
The motivation for this CIPEC device design -in addition to the aforementioned economic advantage 7 -is an expected increase in efficiency as the rejected heat of the PV cell (energy above the band energy which is converted to heat), see Fig. 1b, is utilized for preheating of the reactant. In the case where a heat-driven, self-tracking concentrator is used, 16 the radiation with longer wavelength (above the smallest bandgap of the solar cell materials) can additionally be utilized for the tracking of the concentrator (as detailed in the ESI), further increasing the system efficiency.
We developed a 2D multi-physics model of the CIPEC device by customizing the combination of a commercial finite element/volume solver 17 and an open source finite difference Newton Raphson solver, 18 coupling local mass and heat transfers for the electrochemical component of the device to the detailed multi-physics model of the photoabsorber. The simulation domain consisted of the xy-plane as depicted in Fig. 1a. The model supports the development of design and operational guidelines to maximize hydrogen production, energetic efficiency, and device durability, and to minimize size and cost.

Governing Equations and Methodology
Electromagnetic wave propagation (EM).-The combined form, Eq. 3, of Maxwell curl Equations 1 and 2, 19 is solved using the MUltifrontal Massively Parallel Sparse direct Solver (MUMPS) 17 via the finite element method. The electrical field vector, E, and the magnetic field vector, H, are solved for a finite number of wavelengths spanning the entire spectral range of the incoming solar irradiation.
The real part of the calculated electric and magnetic vector fields, (E λ ) and (H λ ), are used to calculate the time-averaged Poynting vector, 19 S avλ , Eq. 4, and the corresponding optical generation rate, G optλ , in the semiconductor region, Eq. 5, However, instead of the divergence calculation we use the numerically to calculate the generation rate for each wavelength where {ε λ } represents the imaginary part of the complex permittivity of the material for the respective wavelength. The optical quantum yield, η opt , is assumed to be 1 for photons with energies larger than the bandgap of the material and zero for less energetic photons. The wavelength dependence of η opt is different for each of the three (two) semiconductors of the triple (dual) junctions of the solar cell and is determined by their respective bandgaps. The overall The EM simulation domain consists of the water channel and PV regions for the CIPEC device, as it is assumed that the light doesn't penetrate into the electrochemical cell. The irradiation flux at the top boundary of the device is provided by the concentrator's output, which is equal to C · I in , where C is the effective concentration ratio and I in is the solar intensity equal to the AM1.5G solar spectrum with a flux of 1000 W/m 2 . Floquet periodicity is used to account for realistic propagation of the plane wave.
Out of the total solar energy incident on the PV, part is absorbed according to the material's complex refractive index and the remainder, which is not absorbed, is lost. The absorbed energy comprises two parts: a part due to photons above the bandgap energy (E a > E g ), and a part due to photons below the bandgap energy (E b < E g ). It should be noted that photons with energies below the bandgap energy of the semiconductor can be absorbed according to the extinction coefficient of the material. However, these absorbed photons won't contribute to the EHP generation and instead result in lattice perturbations which contribute to the heating of the material as described by Q R and Q M , discussed below. Out of E a , one part leads to EHP generation (i.e. equal to E a (E g /hv) for each photon) and the other part (i.e. E a (1− E g /hv)) comprises the energy lost due to thermalization of photoexcited electrons and holes.
We define two fractions, f EHP = E g /hv and f TH = (1− f EHP ), such that f EHP (−∇ · S avλ ) gives the optical generation rate when divided by respective E g , consistent with the calculations in Eq. 5, and f TH (−∇ · S avλ ), the heat dissipation density due to the thermalization process (Q TH ).
Accordingly, the net heat source term coming from EM wave propagation is defined to be the sum of the individual source and sink terms: Q R are the electrical (resistive) losses, and Q M are the magnetic losses, given by: where J = σE and B = μH. For simulation domains not involving semiconductors, the net heat source is simply The overall Q EM is calculated by summing the individual wavelengths' net heat source terms. The data of the complex refractive index for Ga 0.51 In 0.49 P is from Schubert et al., 20 and data for water and GaAs from Palik. 21 The optical constants for various layers of the aSi-ucSi-ucSi triple junction solar cell were adopted from Collins et al. 22 The spectrally resolved data has been plotted in Fig

Semiconductor charge transport (SC).-The Poisson equation,
Eq. 9, and current conservation equations, Eqs. 10 and 11, are solved simultaneously, incorporating the definitions of both the electron current density vector, J n , and the hole current density vector, J p , from Eqs. 12 and 13, respectively. 23 ε s is the static electric permittivity, q is the electronic charge, and N + D , N − A are the concentrations of ionized donors and acceptors, respectively. R and G are the carrier recombination and optical generation rates (Eq. 5).
In Eqs. 12 and 13, n and p are the free electron and free hole densities, ∇E c and ∇E v are the gradients of the local conduction and valence band energies, T is the lattice temperature, μ n and μ p are the electron and hole mobilities, and D n.th and D p.th are the thermal diffusion coefficients. The function M is defined as M(α) = α/(F -1/2 (F 1/2 −1 (α))) with F -1/2 as the Fermi-Dirac integral. The total current, and thus its variation with voltage, is obtained from the sum of the hole and electron current densities. The recombination term in the current conservation equations considers Shockley-Read-Hall, Auger, and direct recombination phenomena. Temperature independent and temperature dependent baseline parameters are presented in Table I and section S.2.
The existence of continuous density of states (DOS) in the bandgap of a-Si:H/uc-Si:H, with no well-defined conduction band (CB) and valence band (VB) edges, is due to the spatial disorders in its atomic structure. The localized states in the mobility gap of a-Si:H are represented by the CB and VB tail states and the defect states. These states strongly influence the trapping and recombination phenomena. 24 The Urbach tails and defect states in amorphous/microcystalline semiconductors have been incorporated using wxAMPS, 18 and the underlying physics and methodology is discussed in section S.1 of the ESI.
The internal heat source term coming from ohmic losses and net recombination losses in the semiconductor is given by [14]  where R is the net recombination. We used Fermi-Dirac statistics and finite volume solvers 17 to solve Eqs. 9 to 11. The BCs for Eqs. 9 to 11 are depicted in Fig. 2b. They describe the different ohmic contacts used for each part of the PV. Each p-n junction is simulated separately and then the overall PV's current voltage characteristic is generated assuming series connection of these two or three p-n junctions. The tunnel diode connecting the two p-n junctions is not modeled, assuming that it has negligible optical and resistive losses. This assumption is justified by an experimental measurement (which includes all existing resistances) of the 1 cm 2 InGaP/InGaAs solar cells fabricated by CESI S.p.A. Italy, using a Pasan BV-81 AM1.5D concentration cell tester at 25 • C. The measured fill factor (FF) at C = 500 was 0.86 and the solar-to-electricity efficiency was 33.7%, showing that the optimized cells exhibited minimal losses even at high concentrations.

Electrolyzer charge transport (EC) and reacting fluid flow (RFF)
.-Charge transport in the electrode and electrolyte (subscript k = l for the ionic conductor and s for the electronic conductor) is given by: 25 [16] obeying electro-neutrality, i z i c i = 0. The electrochemical reaction at the electrode-electrolyte interface is accounted for via the reaction current, i loc , modeled via the Butler-Volmer expression, 25 RT , [17] where η act,m is the activation overpotential for reaction m, and the total overpotential at the electrode-electrolyte interface is given by for m = 1,2, accounting for the anodic one-step oxygen evolution reaction (OER), [19] and cathodic one-step hydrogen evolution reaction (HER), The charge conservation equations, Eqs. 15 and 16, are solved via finite element methods and a MUMPS solver. 17 The corresponding boundary conditions are: a positive electric potential applied to the anode side, a cathode maintained at zero potential, and insulated sidewalls of the simulation domain. The fluid flow and mass transport in the channels and the porous gas diffusion layers (GDLs) are modelled by the Navier-Stokes equation (Eq. 21) with Darcy extension, 26,27 [21] and species transport (Eq. 22) is modelled by the Maxwell-Stefan diffusion model 28 for a low density fluid mixture, with the diffusivities replaced by the binary diffusivities for the existing species pairs, where ω i is the mass fraction and j i is the mass flux vector relative to the mass average velocity vector u given by 29 R i is the rate expression describing production or consumption, D i T are the thermal diffusion coefficients, d k is the diffusional driving force acting on species k, and D ik are the multicomponent Fick diffusivities. Argon is incorporated as the sweeping gas for both anodic and cathodic chambers. The binary diffusivities and effective thermal conductivity, density, and specific heat capacities in each chamber vary with temperature as detailed in section S.2. The species transport equation, Eq. 22, is solved using a MUMPS solver 17 via finite element methods. BCs for Eq. 22 are depicted in Fig. 3a, describing the inlet and outlet conditions of the different channels. The water at the output of the water channel is fed to the anodic channel at its inlet with normal velocity, v water .
Heat transfer (HT).-We solve the steady-state energy conservation equation, 30 ρC p u · ∇T = ∇ · (k th ∇T ) + Q H , [24] in order to calculate the temperature field in the PEC device. Finite element methods and a PARDISO solver 17 are used for its solution. Q H = Q EM + Q PV + Q EC describes the heat source for the integrated PEC which includes the heat from electromagnetic heating (Eqs. [6][7][8], semiconductor transport (Eq. 14), and solid/liquid-phase charge transport and electrochemical reactions, 25 Eq. 25: The velocity vector, u, is zero for the solid components. The SC (section Semiconductor charge transport) and EC (section Electrolyzer charge trasnport and reacting fluid flow) simulations do not explicitly solve for energy conservation, but are based on charge and current conservation. This results in the appearance of two heat loss terms, Q PV_lost and Q EC_lost , which are defined as Q i_lost = P input,i -Q i -P output,i , for i = PV, EC. The total heat source in the integrated device domain is defined as Q EM + Q PV + Q EC . Q PV_lost and Q EC_lost are not considered in the heat source term. For porous media (e.g. GDL), the thermal conductivity, k th , is replaced by an effective conductivity, k eff , accounting for a volumeaveraged solid conductivity. The BCs for energy conservation are depicted in Fig. 3b, and consist of thermal insulation on the sidewalls and back side of the simulation domain, and heat flux (natural convection) on the top side. For the PEC device simulation, the concentrator is modelled by a thermal resistance approach. 30 Thus, we use an effective heat transfer coefficient of h eff = (1/h coeff + L conc /k conc ) −1 , with L conc representing the mean photon flux length in the concentrator. The temperature-dependent heat transfer coefficient for natural convection at the top of the water channel is calculated using h coeff = Nu · k air (T film )/δ where Nu is the Nusselt number, 31 k air is the temperature dependent thermal conductivity of the air, and δ is the characteristic length.
Meshing strategy and computational expense.-We adopted different meshing strategies for different physical simulation modules in order to minimize the overall solution time of the coupled multiphysics problem and assure mesh independence of the solution. The most restrictive requirement came from the EM simulations (Eq. 3), as minimum element size in the direction of incident radiation had to be at least 1/5 (or 1/10) of the incoming wavelength for second order elements (or first order elements). Additionally, this minimum element size was adjusted by a factor of 1/n ref accounting for the change in the refractive index as light travels from one material to another. A variable meshing approach was used for different wavelengths to minimize the computational time. The solution of the other transport and conservation equations, Eqs. 9-11, 15-18, and 24, were similarly optimized in order to allow for mesh convergence and to minimize the computational expense. A distributed mesh with arithmetic sequenced symmetric distribution and with an element ratio of 9 was used in the direction perpendicular to the flow in order to resolve the boundary layer for fluid flow, mass transport, and heat transfer calculations. A dense meshing approach in the catalyst layer of EC and around the junction interfaces in the PV was utilized. The mesh size was increased in the other domains for the solution of Eqs. 9-18 in order to reduce computational time. A workstation with 128 GB RAM and 12 cores was used to solve the coupled equations. Approximately 15 hours were required (with 150 wavelength bands of ∼7 nm and ∼41 nm in the above bandgap energy and below bandgap energy spectrum, respectively) to find a converged solution for reference case I (see following section). Typically, 3 global (blue loop in Fig. 4), 3 local (dark red loop in Fig. 4), and 20 internal (dotted box in Fig. 4) iterations were needed to achieve final convergence. The number of global and local iterations increased with increasing irradiation concentration and with decreasing water mass flow rate, as these parameters significantly influenced the device temperature profile. Fig. 4 is a simulation flow diagram of the coupled model containing six physical modules: EM, HT, FF, RFF, EC, and SC. The flow depicts the simulation flow for the PEC device without the concentrator, and takes as input the irradiation spectrum from any concentrator simulation.

Simulation flow.-
The EM simulation is performed for the given temperature field, starting with a constant initial temperature (T init = 293 K). The resulting heat dissipation density is input to the HT simulation module. The FF module provides the velocity of the water (flowing in the cooling water channel) to the RFF module and the HT module. The HT module provides the temperature input to all the other modules. The EC module -being fed by water mass fraction from the RFF module and temperature from the HT module -provides Q EC to HT and i loc to RFF. In parallel, the semiconductor device simulation is performed starting with T init (using the input of the defect/disorder recombination profile, Rec d, from AMPS). The resulting surface averaged current densities and electric potential from the semiconductor and electrolyzer charge transfer equations are input to the operating point calculation module, which calculates the operating current and voltage decided by the intersection of the PV and EC J-V curves. If the PV temperature distribution, T op (r), at the operating point is not equal to the initial temperature profile of the semiconductor simulation module (T semi (r)), a new T op (r) is provided to the corresponding modules (similar for T EM (r)) and the loop is repeated. This process is iteratively repeated until the temperature converges (T op,old (r) -T op,new (r) < ε rr ) where ε rr is the tolerable error in temperature, which equals 0.1 K in accordance with Eq. 26.
The AM1.5G solar irradiation spectrum was adapted to contain 150 wavelengths ( λ = 7-41 nm) between 280 nm-2500 nm in order to reduce the number of wavelengths simulated but to still capture the spectral variations in the visible and infrared parts of the spectrum. The change in absorbable spectrum with temperature requires the choice of an adaptive wavelength spectrum for each global iteration step. The minimum wavelength difference for a temperature change from T 1 to T 2 is given by [26]

Results and Discussion
Reference case and validation.-We defined two reference cases: reference case I with aSi-ucSi-ucSi and p-n-cathodic-anodic configuration, and reference case II with Ga 0.51 In 0.49 P-GaAs and n-p-anodiccathodic configuration. The dimensions and component characteristics used for the reference case, with aSi-ucSi-ucSi thin film triple junction (V oc ∼ 2.1 V, J sc ∼ 47 A/m 2 , FF ∼ 64%) or Ga 0.51 In 0.49 P-GaAs dual junction (V oc ∼ 2.35 V, J sc ∼ 100 A/m 2 , FF ∼ 93%) example solar cells, are given in Table I. The detailed modelling parameters for the aSi-ucSi-ucSi thin film cell, and temperature dependent baseline parameters used in the simulations, are given in Tables S.1 and section S.2 of the ESI, respectively. The operating normal mean flow velocity of water in the top cooling channel is 0.2 m/s for both reference cases. The choice of material for these two reference cases was motivated by the techno-economic analysis presented in 7 which highlights that both choices have the potential to achieve low cost production of hydrogen. Here, the F factor, defined as A EC /A PV , 7 is chosen to be 1. This choice is not optimal from a performance and cost point of view 7 but represents a closely integrated device. The solar to hydrogen efficiency (STH) is calculated using [27] assuming faradaic efficiencies of 1, negligible product crossover, and E 0 eq equals to 1.23 V (equilibrium potential for pH = 0, T = 298 K, and atmospheric pressure).
Non-integrated, externally wired concentrated PV and electrolysis devices have experimentally demonstrated efficiencies in the range of 17 to 24%. [32][33][34] Due to the lack of CIPEC device demonstration, we separately validate the PV and EC component. The experimental validation for PV and EC components is presented in Figs. 5a and 5b, respectively. The modelling parameters and operating conditions were chosen to be similar to those presented in other studies. [35][36][37][38] In order to validate the mass transport effects in the EC, a comparison was made with fuel cell (FC) data available from Kim et al. 38 assuming the electrolysis to be a perfectly reversible process of FC. The voltage scale was then adapted using V new = (E eq −V exp ) + E eq , where V exp is the experimentally measured voltage, V new is the adapted voltage, and E eq is the equilibrium voltage of 1.208 at 353 K, eq. S.8. This was due to limited availability of experimental data for current saturation and mass transport effects in PEM electrolysis. The materials and dimensional parameters used in the EC simulations are detailed in Table S.2. The simulated and experimental characteristic curves were in good agreement, i.e. maximal differences in the current of the ECs of 8.6% at a voltage of 1.78 V for the black curve and 12% at a voltage of 1.64 V for the green curve were observed. The FF of PV for the simulated case was 91% compared to the 75% for the experimental case. It should be noted that the PV model didn't take into account the contact resistance (instead we assumed an ohmic contact) and the anti-reflection coating used in Lueck et al., 35 explaining the observed fill factor difference between simulated and experimental curves. The calculated FF is within the physical feasibility limits, as the maximal FF for a V oc of 2.35 V is approximately 0.94 as reported by Green. 39 Irradiation concentration analysis.-The characteristic J-V curves (see Fig. 6) for both reference cases show that the STH efficiency, η STH , decreases with increasing irradiation concentration. J sc increases linearly with C. V oc is influenced by the combined effect of increasing temperature and increasing J sc , resulting in its logarithmic increase with increasing C. FF of the Si based PV (∼64%) is significantly lower than that of the III-V based PV (∼93%) at C = 1. The fill factor for both cases decreases (∼4% for reference case I, ∼1% for reference case II) when increasing C from 1 to 1000. J op, and hencė M H 2 (signifying the rate of hydrogen production in g/min/m 2 and calculated by Faraday's law of electrolysis), increase with increasing C, but at a decreasing rate for larger C. PV and EC average temperatures (i.e. T PV and T EC ) increase with C, and consequently the activation overpotentials for the EC decrease with C. These trends for η STH , T PV , and T EC are shown in Figs. 6c-6d for Si and III-V PV based devices.
The optimal irradiation concentration for the III-V material based case is around C opt = 180. Higher concentrations lead to current saturation in the EC and a further increase in C doesn't benefit J op . For the Si case, the optimal concentration is around C opt = 707. For both reference cases, the STH efficiency decreases when increasing C from 1 to C opt . For the III-V based PV cell this decrease is 0.3% while for the Si based PV cell this decrease is ∼38%. This behavior is due to the relatively low FF of the Si based case. Generally, for devices using low FF PVs, PV and EC are limiting for operation under high irradiation concentration. For devices using large FF PVs, the EC's saturation current is the limiting factor.
The heat source contribution (f Q i = Q i / Q i = Q i /P input ) for various sources including: absorption losses in the water channel (Q w ), EM resistive and magnetic losses (Q R+M ), thermalization losses (Q TH ), PV ohmic losses (Q ohm ), PV recombination losses (Q rec ), EC ohmic losses in electrolyte (Q ohm_l ) and electrodes (Q ohm_s ), and EC kinetic  losses (Q kinetic ), are plotted in Figs. 6c and 6d for Si and III-V based cases respectively. Q unabsorbed , Q PV_lost , and Q EC_lost are closing the energy balance. The total heat source in the integrated domain comes from f Q . For the Si based device, the f Q EM and f Q EC increases with increasing C while f Q PV decreases. For III-V material based device, the same trend is followed for f Q EM and f Q PV , but f Q EC initially increases up to C = 265 and then decreases. This f Q EC trend is the result of the C opt being 180 for the III-V based case, because after C = 180 the operating current is not increasing, resulting in a smaller relative contribution of Q EC . The f Q R+M is higher for reference case II because of the refractive index of the materials used which have significant absorption below the bandgap wavelength.
The sum of Q EM + Q PV + Q EC is higher (1.5-2 times) for the III-V based case than the Si one for the choice of design parameters presented in Table I. Despite this, the T PV and T EC were observed to be less than or comparable to reference case I. This implies that the n-p-anodic-cathodic configuration, used in reference case II, provides better cooling. This results from the direct feed of the water from the cooling channel above the n-side of the PV to the anodic channel lying between the p-side of the PV and the EC's anode, making the location central within the major heat sources in the device, resulting in better heat removal. Figs. 6c-6d also present f Q i for a lower mean flow velocity (0.03 m/s), implying that lower mass flow rates are less effective in heat removal, leading to higher device temperatures. The parasitic absorption in the top water channel is found to be minimal (see Figs. 6c-6d) and is in accordance with the findings reported by Döscher et al. 40 for a water channel of thickness of 0.2 mm with overpotentials in the range observed in our device.

Sensitivity analysis with material and dimensional parameters.-
The characteristic operating curves for EC and PV for reference case I, for the variation of parameters including exchange current density, membrane thickness, as well as catalyst and GDL thicknesses, are presented in Fig. 7. Only one parameter is varied at a time, holding all others at the reference case values. The intersection of the EC and PV curves represent the operating point of the integrated PEC device.
The operating current density and the STH efficiency, shown in Table II, increase with increasing exchange current density and catalyst thickness (leading to improved EC kinetics), Figs. 7a-7c. STH efficiency decreases with increasing membrane and GDL thickness as visible in Figs. 7b-7d. Increasing the membrane thickness results in increased mass transport limitations in the EC. The coupled model allows a temperature maintenance in the range of ∼300 K for all parameter variations, even at C = 450. A water inlet mean flow velocity of 0.2 m/s was found be sufficient in effectively maintaining the device temperature at ∼300 K, signifying the importance of thermal  Two-dimensional results.-The 2D modelling for the various physics involved allowed for the generation of an in-depth knowledge of the spatial variation of various physical parameters in the device. This 2D treatment is important because various physical phenomena have orthogonal main directions to one another, for example: the fluid flow's main direction is along the y-axis, perpendicular to the main direction of the charge transport, and absorption in the PV occurs along the x-direction.
A typical temperature distribution in the EC of the integrated PEC is shown in Figs. 8a-8b. The EC consists of the anodic and cathodic GDLs, catalyst layers, and Nafion membrane (in middle). The cathodic side operates at higher temperature than the anodic side. The temperature increases along the positive y-direction, the direction of the fluid flow. This results from the heating of the fluid once it enters the channel and moves toward the positive y-direction. The cathodic side shows higher temperatures because it hosts the exothermic reaction with a positive heat source, Q GDL+cat,c = 2.93 · 10 7 W/m 3 , while the anodic side is at lower temperatures as it hosts the endothermic reaction with a negative heat source, Q GDL+cat,a = −7.50 · 10 6 W/m 3 , both for C = 1000 and mean flow velocity of 0.2 m/s. The same trend at absolute increased temperatures is observed for lower water flow velocities, as depicted in Fig. 8b.
The H 2 mass fraction (wH 2 c) in the cathodic GDL and channel assembly increases relatively homogenously in the direction of sweep gas flow, because the generated H 2 gas has better diffusivity and is swept in flow direction. The input water diffuses through the anodic GDL to the catalytic site where it is oxidized. This results in a decrease of H 2 O mass fraction (wH 2 O) in the anodic GDL and also in the channel in the direction of water flow as shown in Fig. 8d. The ionic flux given by the Nernst-Plank equation decreases in the Nafion membrane from inlet to outlet, forcing the electrolyte current density to follow the same trend, shown in Fig. S.4(a). The local current density, given by the concentration dependent Butler-Volmer equation, is governed by the concentration of H 2 O at anodic side and H + ion on the cathodic side, giving rise to the current density profiles shown in Fig. S.4(c)-(d).
The 2D modelling enables knowledge of the full distribution profile of various physical parameters, which in turn helps in identifying the local maxima and minima. This is beneficial from the perspective of identifying and removing hot spots in the device which may lead to thermal stress and generally affect performance and operational uniformity. The 2D profile of wH 2 O shows that the channel length is optimized for the given conditions. Longer channels do not benefit the operation, as the H 2 O is fully consumed toward the end of the channel. Hence, only a simultaneous and corresponding increase of both the channel length and the mean flow velocity provide additional benefit to the performance.
As illustrated, the correct realization of complex geometries for electromagnetic propagation can only be done with 2D and 3D treatment. 2D modelling provides locally resolved information which is necessary for complete design guidance and optimization of PEC devices.

Conclusions
We developed a coupled 2D multi-physics model and solution methodology to simulate the performance of integrated photoelectrochemical devices using concentrated solar irradiation. The model couples electromagnetic wave propagation, semiconductor charge generation and transport, heat transfer, fluid flow, mass transport, electrolyte and electrode charge transport, and electrochemical reactions. Finite element and finite volume methods were used to solve the governing equations and the corresponding boundary conditions. Complex temperature dependencies were included in the model. The absorbable spectrum changes with temperature and therefore requires an adaptive spectrum, changing for each iteration step, giving rise to a trade-off between precision and computation time. The various heat source/sink term calculations were treated in detail to ensure accurate energy calculations and to allow for subsequent effective thermal management. The model and its simulation flow was fully automatized. Efficient computational power saving techniques (including variable meshing and adaptive radiation spectrum approaches) have been rigorously employed, making our model detailed while maintaining a low computational cost.
Two reference cases were defined, utilizing different photovoltaic (PV) components: i) a triple junction thin film aSi-ucSi-ucSi cell, and ii) a dual junction III-V based Ga 0.51 In 0.49 P-GaAs cell. These cells showed a decreasing trend in STH efficiency with increasing irradiation concentration. The low FF of the aSi-ucSi-ucSi cell resulted in a PEC device performance limited by both the electrocatalysts (i.e. the integrated electrolyzer) and the photoabsorbers (i.e. the PV). The high FF of the Ga 0.51 In 0.49 P-GaAs cell, on the other hand, was limited only by the electrolyzer's saturation current.
In order to maximize the produced amounts of H 2 , we recommend using high irradiation concentrations. Large concentrations additionally benefit the economic competitiveness of the device. 7 However, due to the limiting saturation current of the electrolyzer and the relatively small FF of PV, there exists an optimal irradiation concentration. Larger concentrations (for device configurations with equal EC and PV areas) do not increase performance. The water channel on top of the PV can effectively cool the device if a large enough mass flow rate is chosen. For example, the device could be maintained at around 300 K for C = 450. The optimal mean flow velocity of water was found to be 0.2 m/s (or 40 g/s/m) for all irradiation concentrations. Large water mass flow rates provide a greater benefit in terms of cooling capability compared to smaller rates, and additionally provide more reactant to the electrodes, alleviating mass transport limitations.
J op increases with increasing exchange current density, i.e. catalyst activity, at a particular C. A similar trend is observed for the STH efficiency. Changes in the exchange current density show minimal changes in device temperature due to optimized thermal management. The dimensional properties such as membrane thickness lead to significant changes in the operating points, highlighting the importance of the membrane for thermal management in the integrated device. With increasing membrane thickness, the mass transport limitations are instigated earlier, resulting in reduced saturation currents. Despite large operating point variations, temperature variations are small (a few Kelvin only), as the water mass flow rate ensures proper device cooling. In spite of significant increases of heat sources, these minimal temperature variations attest to functional and efficient device thermal management. STH efficiency and H 2 production increases with increasing catalyst thickness, in contrast to the decreasing efficiency with increasing GDL thickness resulting from increased electronic and diffusional resistance.
To our knowledge, the model we develop in this analysis is the most detailed yet computationally low-cost model reported. Our model allows for the investigation of any complex device design and geometry, and its simulation in fine physical detail. The model shows to be a valuable tool for the design of integrated PEC devices working with concentrated irradiation at elevated temperatures and illustrates that smart thermal management can assist in achieving efficient and low cost production of solar fuel at large volumes. Thermal hot spots in a device operating at high irradiation concentration can be reduced utilizing calculated, spatially resolved temperature profiles, reducing the thermal and operational stress on photoabsorbers or catalysts, and potentially slowing their degradation rate. A more detailed analysis and the quantitative/qualitative benefits of smart thermal management for the integrated design of IPECs will be detailed in follow-up studies. 41