Advanced X-ray radiographic imaging analysis for identifying the water dynamics regime of polymer electrolyte fuel cells

Abstract

Operando X-ray radiography is widely used to evaluate liquid water distributions in polymer electrolyte fuel cells, yet conventional analyses typically condense pixel-resolved time series into spatiotemporal averages, potentially obscuring dynamically fluctuating water states. In this study, we present a computational spatiotemporal statistical framework to extract fluctuating information from sequential radiographic datasets of a cathode gas diffusion layer operating under steady conditions. Beyond mean saturation index, pixel-wise time-series descriptors, including skewness, excess kurtosis, Shannon entropy, frequency, autocorrelation and mutual information, were systematically evaluated to characterise distributional variability, periodicity, temporal persistence and spatial coordination of water fluctuations. The analysis revealed that regions with similar mean saturation can exhibit distinct fluctuation structures. By combining Shannon entropy and temporal autocorrelation, a two-dimensional statistical regime map was established to classify fluctuation patterns into four states. The proposed framework does not directly establish universal water-transport mechanisms but instead introduces a fluctuation-oriented statistical perspective complementary to traditional mean-value analyses. This approach provides a reproducible and extensible methodology for comparative interpretation of operando radiographic datasets and may contribute to future data-driven diagnostic analyses.

---

Introduction

Water management in polymer electrolyte fuel cells (PEFCs) remains challenging because water simultaneously enables membrane hydration and obstructs gas transport when it accumulates as liquid water.¹ Adequate water content in a hydrated membrane is desired to maintain high proton conductivity. Conversely, excessive liquid water accumulation, often within the cathode gas diffusion layer (GDL), restricts oxygen transport to the catalyst layer and results in increased mass-transport overpotential. Balancing ohmic losses against mass-transport losses is therefore important for the development of next-generation high-performance PEFCs.^2,3 To directly observe water distributions and their evolution during fuel cell operation, researchers increasingly use operando X-ray imaging.^4–7 Imaging techniques, whether tomography or radiography, are typically chosen based on the required spatial and temporal resolution.

Tomography has been widely applied to study three-dimensional liquid water behaviours in membrane electrode assemblies (MEAs), particularly within the cathode GDL. By reconstructing three-dimensional fields (and, when repeated, time-resolved four-dimensional datasets), tomography enables detailed analysis of the morphology of water droplets and their connectivity in a cell. For example, Liu et al.⁸ developed a computational algorithm to measure internal water contact angles at fibre–water–air boundaries, demonstrating that GDLs exhibit a state of mixed wettability with a broad distribution of local contact angles. Kato et al.⁹ performed network analysis to classify liquid water distributions into three distinct cluster types, revealing that most internal water paths result in dead ends due to local pore narrowing within the complicated GDL structure. Nagai et al.¹⁰ classified water droplets into four distinct categories, full-connected, top-connected, bottom-connected and non-connected pathways from the catalyst layer to gas channels, demonstrating that GDLs attached to the microporous layers (MPLs) with large pores stabilise the formation of fully connected clusters to minimise overall water saturation and enhance cell performance.

The limitations of tomography are its temporal resolution and data volume. A typical minute-level temporal resolution can obscure fast liquid water transients and make it difficult to interpret dynamic water redistribution, such as flushing, evaporation and condensation, on (sub-)second timescales. Recent advances have improved scan rates (e.g., 4.2 s per scan for ex situ water injection experiments¹¹ and 3 s per scan for operando experiments¹²). In such cases, terabyte-scale data handling and analysis become major bottlenecks. AI-based segmentation and labelling are being developed to reduce the analysis burden.¹³ Recent deep-learning approaches have improved multi-class segmentation of MPL/GDL structures by training with measured porosity data.¹⁴ Furthermore, Shum et al.¹⁵ compared machine learning algorithms, including decision tree learning and convolutional neural networks, for the automated three-dimensional segmentation of liquid water in GDLs to improve accuracy over traditional thresholding methods.

Radiography offers substantial advantages in acquisition speed and data handling compared with tomography. Whereas operando X-ray tomography requires extensive reconstruction and post-processing, operando X-ray radiography can often be processed and interpreted on much shorter timescales. The trade-off is that radiography provides a line-of-sight projection, integrating attenuation along the X-ray irradiation direction. As a result, details of water droplet morphology and their connectivity are missing. Despite this dimensional constraint, radiography effectively provides time-averaged water saturation within user-defined regions of interest (ROI). For example, Banerjee et al.¹⁶ demonstrated the utility of X-ray radiography for tracking real-time changes in membrane hydration. Yoshimune et al. established a diagnostic method using a specially designed cell that visualised water accumulation within the catalyst layer¹⁷ and simultaneously evaluated the GDL degradation rate using a water saturation indicator across all MEA components during accelerated stress testing.¹⁸ Furthermore, Kato et al.¹⁹ developed a fuel cell design to simultaneously quantify cross-flow rates and liquid water distribution in GDLs, demonstrating that thicker GDL substrates can effectively reduce liquid water accumulation when combined with a cross-flow-inducing parallel flow field. Recently, Chadwick et al.²⁰ revealed that the distribution of water droplets is influenced by gravity when using a cell with a rotating mechanism, and that this effect is critical in anode GDL water management.

A critical limitation remains in conventional radiographic studies. Frame-by-frame pixel-level fluctuations are typically condensed into spatiotemporal averages to obtain a single representative water saturation value. This averaging process is robust for materials screening, but it inherently masks fluctuation patterns, such as variability, persistence, intermittency and coordination, that reflect a dynamic water state. Despite advances in the radiographic imaging hardware, less attention has been paid to how radiographic time-series data are analysed beyond spatiotemporal averaging. Regions exhibiting similar mean saturation may nonetheless display different temporal fluctuation characteristics. During PEFC operation, water configurations likely rearrange at microscales even when the cell voltage appears stable, and these hidden dynamics may modulate local oxygen transport resistance and trigger flooding onset. As these variations are not directly reflected in a single scalar metric, a multi-dimensional statistical description of water dynamics is potentially beneficial. Thus, we here propose a method of extracting dynamic information from sequential operando X-ray radiographic images beyond conventional spatiotemporal-averaged saturation metrics. Using an available radiography dataset, we first reviewed standard mean-value analyses and then computed pixel-resolved time-series indexes, including skewness, excess kurtosis, Shannon entropy, Fourier spectral features, temporal/spatial autocorrelation and mutual information, to elucidate intermittency, persistence, and spatial coordination of water saturation in the cathode GDL. Finally, we integrated these metrics to classify local water behaviour into dynamic regimes associated with water drainage and stagnation.

Methodology

Dataset

The operando X-ray radiography dataset analysed in the study was previously reported in ref. 21. The radiographic measurements were conducted at an X-ray energy of 11.4 keV with a pixel resolution of 1.3 µm. Data acquisition was conducted after stabilising the fuel cell for 3 min under operating conditions of 40 °C and 100% relative humidity at a constant cell voltage of 0.2 V. Gas flow rates were maintained at 100 cm³ min⁻¹ for synthetic air at the cathode and 100 cm³ min⁻¹ for hydrogen at the anode, without applying additional back pressure at the outlets. The dataset consisted of N = 128 radiographic frames acquired with a frame period Δt = 1.5 s, including a 0.5 s inter-frame interval. A pixel-resolved water thickness map at a specific frame number t(x,y,n) was obtained by normalising wet radiographs with a time-averaged dry reference radiograph acquired under identical operating parameters except under the open-circuit condition. The equivalent water thickness was calculated using the Beer–Lambert law:
where I_wet(x,y,n) and I_dry(x,y) denote the pixel-wise X-ray intensities in the wet and dry states, and μ is the attenuation coefficient of liquid water (0.344 mm⁻¹ at 11.4 keV taken from ref. 22). The resulting dataset was stored as pixel-resolved water thickness fields for each frame.

Conventional data analysis

The local water thickness was converted into water saturation s(x,y,n) using the effective X-ray path length l (6 mm) and the GDL porosity ε:

Fig. 1 presents representative saturation maps illustrating the apparent spatiotemporal evolution within the cathode GDL (see Video S1 in the SI). Within the observed domain, the upper boundary interfaces with the cathode flow field, while the lower boundary is in contact with the cathode catalyst layer. Specifically, the upper-central region (x = 400–800 px) is adjacent to the gas flow land (rib), whereas the upper-left and upper-right regions are connected to the gas flow channels. Higher saturation levels are observed in the central region beneath the gas flow land. The liquid water distribution appears visually quasi-stable across sequential snapshots (Fig. 1a–h).

Fig. 1 Spatiotemporal evolution of water saturation in a cathode GDL under a steady-state fuel cell operation. (a)–(h) Sequential maps captured from t = 0 s to 190.5 s illustrate the temporally dependent liquid water distribution. The upper boundary of the domain (y = 0 px) is in contact with the cathode flow field, while the lower boundary (y = 110 px) interfaces with the cathode catalyst layer. The top central region (x = 400–800 px) is adjacent to the gas flow land, whereas the upper-left and upper-right regions are in contact with the gas flow channels (x = 0–400, 800–1100 px). These saturation fields serve as the dataset for subsequent analyses.

The average water saturation at a given frame s̄(n) within an arbitrary ROI with a width W and height H is calculated as follows:

Fig. 2 shows the temporal evolution of the averaged water saturation. The corresponding current density profile is also presented to indicate the electrochemical stability during operando radiographic acquisition. Although small temporal fluctuations in current density were observed around approximately 1 A cm⁻², the overall electrochemical response remained relatively stable throughout the measurement period. As indicated by both the spatial maps (Fig. 1) and the corresponding time-series profile (Fig. 2), the spatially averaged saturation remains nearly constant under steady-state operating conditions. Temporal fluctuations are typically reduced to a single representative mean saturation value s̄ under steady-state conditions:

Fig. 2 Temporal evolution of the rib-region averaged water saturation within the cathode GDL domain presented in Fig. 1, together with the corresponding current density profile.

In this case, the resulting spatiotemporally averaged water saturation over the analysis period is 0.412 ± 0.007 (mean ± standard deviation). Such averaged saturation values are routinely reported in operando X-ray radiography studies to characterise cell performance under varying key operating parameters, such as temperature, relative humidity and gas flow rate.

Proposed data processing

Analytical framework overview. Conventional radiographic analyses typically reduce pixel-resolved saturation time series to a single spatiotemporal mean value under steady-state conditions. While such averaging provides a robust scalar descriptor, it suppresses the temporal variability retained in the underlying image stack.

To characterise these suppressed dynamics, we developed a multi-metric analytical framework operating directly on the pixel-resolved saturation time series s(x,y,n). The framework groups metrics into three categories: distribution metrics, temporal metrics, and spatial metrics. All analyses were implemented in Python and operate on water saturation image stacks with an arbitrary spatial resolution and size.

Distribution metrics. These metrics quantify the statistical properties of saturation fluctuations without considering temporal ordering at each pixel.

Skewness S(x,y) is calculated as follows:²³

where s̄(x,y) and σ(x,y) are the temporal mean and standard deviation. Skewness quantifies distribution asymmetry. As a standardised statistic, it is scale invariant and therefore not further normalised. A positive or negative S value indicates a distribution biased toward lower or higher saturation regions, respectively, with a long tail in the opposite direction.

Excess kurtosis K(x,y) is described as follows:²³

Excess kurtosis measures tail heaviness and intermittency. As with skewness, it is inherently normalised and therefore not rescaled. In this definition, K = 0 corresponds to a standard normal distribution. A positive K value indicates a leptokurtic distribution with a sharp peak and heavy tails. Conversely, a negative K value represents a flatter distribution with thinner tails.

Shannon entropy²⁴ is computed based on the saturation distribution. Importantly, this metric evaluates the statistical distribution of saturation values and does not account for temporal ordering, such as persistence or periodicity. For each pixel, the saturation time series is discretised into M equally spaced bins spanning the global range of saturation values within the dataset. Here, c_i(x,y) denotes the number of frames in which the saturation value falls within bin i, and

denotes the corresponding empirical probability. The Shannon entropy H(x,y) is then defined as follows:
where bins with p_i = 0 are excluded from the summation. As the Shannon entropy is estimated from finite samples, its magnitude depends on both the time-series length N and the bin number M. To mitigate estimation bias, the bin number was chosen such that M ≪ N (here M = 32, N = 128). All pixels were evaluated using identical bin edges to ensure spatial comparability. To allow comparison across datasets and parameter choices, the entropy was normalised using the theoretical maximum entropy H_max = log₂ M, corresponding to a uniform distribution:²⁵

The resulting normalised entropy E(x,y) provides a dimensionless measure of distributional diversity. A high E indicates that the saturation values span a broad range of states over the analysed period, whereas a low E indicates that the saturation remains confined to a narrow range of values.

Temporal metrics. These metrics evaluate the temporal fluctuations of the saturation time series at each pixel.

Temporal saturation signals were transformed into the frequency domain using the fast Fourier transform (FFT) algorithm²⁶ to identify dominant periodic components X(v) and the power spectral density P(v). As the dataset consists of uniformly sampled discrete observations, the discrete Fourier transform (DFT)²⁷ was applied:

The power spectrum was calculated as follows:

P(v) = |X(v)|²

and
where v(= 0, 1, K, N − 1) is the frequency index and f_v represents the frequency.

In this study, the local water saturation was first averaged along the horizontal direction (x-axis):

The DFT was applied to s̄(y,n) to obtain frequency components for each y-coordinate. In addition, pixel-resolved spectral analysis was performed on s(x,y,n) to construct spatial maps of the spectral power at selected frequencies. No detrending was applied prior to the FFT to retain and analyse the direct current (DC) component at 0 Hz.

Temporal autocorrelation TA(x,y,τ) was employed to evaluate the self-similarity over a specific frame lag τ.²⁸ It is written as follows:

Higher TA(x,y,τ) values indicate stronger temporal persistence of the saturation fluctuations.

Spatial metrics. These metrics calculate the spatial connectivity of saturation fluctuations between pixel pairs.

To quantify the spatial distribution and the characteristic length scale of saturation patterns, spatial autocorrelation²⁹ was employed. The two-dimensional spatial correlation C(Δx,Δy,n) is defined as follows:

where Δx and Δy represent spatial lags in the horizontal and vertical directions. Then, these two-dimensional maps are averaged to yield the time-averaged spatial correlation map C̄(Δx,Δy):

To obtain an orientation-independent representation, radial averaging was performed to derive the spatial autocorrelation function SA(r) as follows, using a radial distance :

SA(r) = 〈C̄(Δx,Δy)〉_r

The analysis is performed up to a maximum spatial lag distance r_max = 50 pixels. The characteristic correlation length is defined as the distance at which SA(r) decays below a selected threshold (typically 0.5).

To quantify the shared information between adjacent pixels, mutual information analysis was employed. For two saturation time series observed at pixel locations i and j,

s_i = {s_i(n)}^N_n=1, s_i = {s_i(n)}^N_n=1

The mutual information³⁰ is defined as

where p(x,y) is the joint probability distribution of the paired samples (s_i(n), s_j(n)), and p(x) and p(y) are the marginal distributions. Although this expression assumes discrete variables, mutual information is estimated from continuous-valued time series using the Kraskov k-nearest-neighbour estimator³¹ with k = 5. To characterise spatial coupling, the mutual information connectivity at each pixel MC(x,y) was defined as the average mutual information shared with neighbouring pixels located within a specific search radius r as shown below:
where N_r denotes the number of neighbouring pixels within the search radius and d_i,j is the Euclidean distance between pixels i and j.

High MC(x,y) values indicate stronger non-linear temporal coupling between neighbouring pixels.

To evaluate the decay of dynamic coordination over distance, a mutual information decay profile was constructed by averaging mutual information values for pixel pairs separated by a radial distance. The mutual information decay function MD(r) is given by

where
denotes the set of pixel pairs whose separation lies within the radial bin centred at distance r.

Regime classification

The computed statistical metrics was integrated to categorise the spatiotemporal states of liquid water into operational regimes within a two-dimensional statistical space. The regime map was constructed combining values of Shannon entropy, representing distributional variability, and the normalised temporal autocorrelation, representing temporal persistence. For each metric, a threshold value th is defined as the median of its spatial distribution over the analysed GDL area. The median was selected as a robust, data-driven threshold that does not assume a specific parametric distribution and ensures a balanced partitioning of the spatial domain. The resulting classification was therefore relative within the present dataset and did not imply universal physical thresholds. Each pixel was assigned to one of four regimes according to its position in the variability–persistence plane (Table 1).

Table 1 Regime classification criteria

Classification	Definition
Low variability/low persistence	E(x,y) ≤ th
	TA(x,y,τ) ≤ th
Low variability/high persistence	E(x,y) ≤ th
	TA(x,y,τ) > th
High variability/low persistence	E(x,y) > th
	TA(x,y,τ) ≤ th
High variability/high persistence	E(x,y) > th
	TA(x,y,τ) > th

---

Results

This section presents the application of the proposed spatiotemporal diagnostic framework to the operando X-ray radiography dataset introduced in Fig. 1. Conventional analysis based on the time-averaged saturation (Fig. 2) indicates quasi-steady behaviour. In contrast, the following results demonstrate spatial and temporal heterogeneity revealed by the distributional, temporal, and spatial metrics. In this section, we focus on presenting the analysis results. Detailed physical interpretation and mechanistic discussion are provided in a separate Discussion section.

Fig. 3 presents the spatial maps of (a) skewness, (b) excess kurtosis and (c) Shannon entropy obtained from the pixel-wise time-series analysis. Each statistical metric was calculated independently at every pixel based on the temporal distribution of saturation values. The geometric structure is symmetric with respect to x ≈ 600 px. The regions x = 0–400 px and x = 800–1100 px correspond to channel-connected areas of identical width, while x = 400–800 px represents the land region. The interfaces between these regions are located at approximately x ≈ 400 px and x ≈ 800 px. No pronounced systematic variation was observed along the vertical (through-plane) direction; therefore, the following discussion focuses primarily on the horizontal (in-plane) distributions.

Fig. 3 Spatial mapping of temporal dynamics: (a) skewness, (b) excess kurtosis, and (c) Shannon entropy of water saturation fluctuations.

The skewness map (Fig. 3a) exhibits variations, but the spatial pattern is not perfectly symmetric about the geometric center. In the channel-connected region at the right side, positive values are locally observed, although negative values are also present. In contrast, the land region shows values predominantly distributed around −0.2–0.2, without extensive areas of strongly positive or negative skewness. Near the interfaces, sign changes in skewness are detected: positive values are observed at the left side, while negative values are also found at the opposite side.

The excess kurtosis distribution (Fig. 3b) does not display a clear separation compared to the skewness map. Across the entire domain, values are generally centred around zero, with localized high-value spots reaching up to approximately 1.6 near the right channel-connected and channel–land interface regions. These appear as scattered, high-intensity regions rather than continuous bands.

Among the three metrics, the Shannon entropy exhibits the most distinct in-plane contrast. The land region (400–800 px) shows consistently higher values (approximately 0.5–0.8), with peaks near the interfaces at x ≈ 400 px and x ≈ 800 px. In contrast, the channel-connected regions are characterized by lower entropy values (approximately 0.2–0.4), resulting in a clear lateral differentiation.

The three basic statistical descriptors do not simply reproduce the same spatial pattern. The Shannon entropy most clearly delineates the land region, skewness highlights sign asymmetry and local reversals near the interfaces, and excess kurtosis emphasizes spatially localized high-intensity events. These differences indicate that each metric captures a distinct statistical aspect of the saturation time-series distribution, as described in detail in the Discussion section.

The application of FFT analysis characterizes the temporal dynamics of water saturation averaged across the in-plane width of the cathode GDL. Fig. 4a presents the spatial temporal map along the through-plane, from the flow field interface (y = 0 px) to the catalyst layer side (y = 110 px). An overall decreasing trend in saturation is observed from the flow field toward the catalyst layer. The highest values are located near the flow field interface (approximately y ≈ 0–30 px, s ≈ 0.28–0.30), while the lowest values appear near the catalyst layer side (y ≈ 80–110 px, s ≈ 0.05–0.10). This lower region corresponds to the MPL. Although the profile exhibits a downward trend toward the catalyst side, the decrease is not monotonic. Instead, the saturation profile shows an undulating structure. Along the temporal axis, no pronounced long-term drift is observed, and the system remains quasi-steady over the measurement period. Temporal fluctuations are present throughout the thickness, but their amplitudes are small relative to the mean saturation gradient in the through-plane axis.

Fig. 4 FFT analysis of water saturation in the cathode GDL. (a) Spatial temporal map showing the saturation distribution from interface between the flow field toward the catalyst layer. (b) Spatial frequency map representing the power spectral density.

The corresponding frequency-domain representation is shown in Fig. 4b as a spatial frequency map of the power spectral density. Across the entire GDL thickness, spectral power is predominantly concentrated at 0 Hz, indicating that slow temporal variations dominate the dynamics. However, the spectral distribution becomes broader in the intermediate and lower regions of the GDL substrate (approximately y ≈ 40–90 px). In these areas, weak spectral components extend toward higher frequencies compared to the flow field side. Near the flow field interface, the spectrum remains more strongly confined to the lowest frequencies.

To examine frequency-specific behavior, Fig. 5 presents spatial distributions of the spectral power density at 0 Hz, 0.026 Hz, and 0.104 Hz. The 0 Hz component is proportional to the time-averaged saturation at each pixel, and thus this map reflects the spatial distribution of saturation averaged over the entire measurement period, equivalent to the field obtained by time-averaging the full image stack (see Video S1), of which Fig. 1a–h provide representative snapshots. At 0.026 Hz (period ≈ 38.4 s), localized high-intensity regions emerge near the channel–land interfaces (x ≈ 400 px and x ≈ 800 px). These features are most pronounced within the GDL substrate region. At 0.104 Hz (period ≈ 9.6 s), the intensity is nearly zero across the entire range, indicating that the intensity observed at this corresponding frequency in Fig. 4b is attributed to noise. This spatial pattern corresponds to the enhanced spectral intensity observed at the same frequency range in the spatial frequency map (Fig. 4b). While Fig. 4b presents the frequency characteristics averaged across the in-plane width, Fig. 5 spatially resolves those frequency-specific features, revealing their localization along the channel–land boundaries.

Fig. 5 Frequency-specific power spectral density maps at 0, 0.026 and 0.104 Hz.

Fig. 6 shows the spatial distributions of the temporal correlation coefficient of saturation at different time lags. At Δt = 1.5 s (which correspond to 1 frame lag), relatively high correlation values (0.6–1.0) are observed across most of the domain. Elevated correlation bands appear within the land region (x ≈ 400–800 px) and near the channel–land interfaces. In contrast, the channel-connected regions exhibit locally reduced correlation values. At Δt = 7.5 s, the overall correlation decreases. The reduction is particularly pronounced in the central part of the land region (x ≈ 500–700 px), where the previously continuous high-correlation areas become fragmented and weaker. In contrast, relatively higher correlation values persist near the interfaces (x ≈ 400 px and x ≈ 800 px). The spatial pattern becomes increasingly patchy, and the overall contrast is moderated. At Δt = 15.0 s, the correlation further decreases, with values predominantly distributed in the range of 0.2–0.5. In the central land region, high-correlation structures disappear. However, localized regions of elevated correlation remain near the channel–land boundaries, especially on the right side. Overall, the spatial decay of temporal correlation with increasing time lag is non-uniform. While correlation diminishes rapidly in the central land region, it persists for longer durations near the interfaces.

Fig. 6 Spatial mapping of temporal autocorrelation at multiple time lags at 1.5, 7.5 and 15.0 s.

Fig. 7 shows the spatial autocorrelation coefficient of the saturation as a function of radial lag distance. The correlation is unity at r = 0 by definition and decreases monotonically with increasing r. The decay is gradual rather than abrupt, and the correlation remains above 0.5 even at r = 50 px. The curve does not exhibit a pronounced inflection point within the measured range, indicating a smooth reduction in spatial similarity with distance. These results demonstrate that the saturation field retains substantial in-plane spatial correlation over length scales of several tens of pixels.

Fig. 7 Spatial autocorrelation profile as a function of lag distance.

Fig. 8a presents the spatial distribution of mutual information connectivity. High mutual information values (>0.5 bits) are primarily localized near the channel–land interfaces (x ≈ 400 px and x ≈ 800 px). The maximum values reach approximately 0.7 bits, demonstrating strong spatially heterogeneous coupling across the domain. Fig. 8b shows the decay of mutual information as a function of radial separation distance. At shorter distances, mutual information is approximately 0.32 bits and decreases over the first few pixels. Beyond this initial drop, the decay becomes more gradual. Even at r = 30 px, this metric remains around 0.12 bits, indicating that spatial dependence persists over length scales of several tens of pixels rather than vanishing within the measured range.

Fig. 8 Mutual information analysis of spatial coordination. (a) Spatial map of mutual information connectivity. (b) Mutual information decay as a function of spatial distance.

Fig. 9 presents the statistical regime map constructed using Shannon entropy (variability) and temporal autocorrelation (persistence) values, with thresholds defined by the median values of each metric. Each pixel was successfully classified into one of four operational regimes in the entropy–persistence plane. Areas previously appearing homogeneous in the time-averaged saturation distribution are separated into distinct statistical states based on variability and persistence.

Fig. 9 Statistical regime map constructed using Shannon entropy (variability) and normalised temporal autocorrelation (persistence). Thresholds for each metric were determined using median values over the analysed dataset. Each pixel was classified into one of four regimes in the entropy–persistence plane, defined in Table 1.

Discussion

Projection constraints and statistical interpretation

In this study, a time-series water saturation dataset obtained from operando X-ray radiography was processed computationally using statistical descriptors capturing distributional, temporal, and spatial characteristics. As transmission images provide a depth-integrated projection of three-dimensional structures, multiple three-dimensional liquid configurations are folded into the recorded two-dimensional images. This intrinsic limitation may affect the interpretation of spatial dependence metrics.

Table 2 summarises the sensitivity of each statistical indicator to projection effects. The mean saturation metric represents an integrated physical measure and its interpretation is straightforward. The distribution metrics (here, skewness, excess kurtosis and Shannon entropy), defined from the marginal distribution of each pixel time series, do not preserve temporal ordering and primarily reflect the diversity of visited saturation states. While projection may smooth amplitude variations, these metrics are comparatively robust for assessing distributional complexity. Temporal metrics (here, FFT and autocorrelation) include the superposition of independent structures along the depth direction; however, the temporal axis itself remains preserved. Consequently, relative comparisons of persistence across spatial locations remain meaningful within the projection framework. In contrast, spatial metrics (here, spatial autocorrelation and mutual information) are more strongly influenced by projection. Depth integration may artificially enhance apparent spatial coupling. Such effects can lead to overestimation of three-dimensional connectivity, as demonstrated in Fig. 7 and 8b. In particular, the mutual connectivity map (Fig. 8a) displays a plausible distribution, despite overestimation. Spatial maps can offer valuable insights, but they also carry the potential to lead to misleading interpretations. This study employs Shannon entropy (variability) and temporal autocorrelation (persistence) for regime classification (Fig. 9) to ensure robustness against artifacts arising from prediction. While various indicators can be integrated computationally, these specific metrics were prioritized to minimize potential risks. The resulting statistical regime map reveals spatial segmentation that is not discernible from a single statistical metric alone. The complexed indexed metric by combining variability with persistence demonstrates that dynamic behavior cannot be inferred solely from average values. Complexed indicators combining volatility and persistence (or other metrics) demonstrate new possibilities for describing the water dynamic equilibrium.

Table 2 Influence of radiographic projection on statistical metrics

Metrics	Sensitivity to depth integration	Potential bias introduced by projection
Mean value	Low	Depth integration introduces minimal interpretational ambiguity
Distribution	Moderate	Depth averaging may flatten the distribution
Temporal	Moderate	Depth superposition may modify apparent persistence
Spatial	High	Correlation length may be overestimated by projection

---

Implications for water dynamics in PEFCs

The statistical indicators introduced in this study reveal spatial heterogeneity in water distribution that is not captured by conventional spatiotemporal averaging analysis. The discussion presented in this subsection provides physically possible interpretations of water dynamics in an operating PEFC, based on the statistical maps shown in this article. For this reason, it is intentionally presented separately from other sections.

The skewness (Fig. 3a) map reflects the asymmetry of saturation distribution, which suggests the presence of water intrusion and drainage events. Positive skewness observed near the left gas channel suggests that dry GDL pores sometimes experience water intrusion events for continuing water drainage toward the gas channel. In contrast, negative skewness near the channel–land interface on the right side indicates that regions maintaining relatively high saturation still undergo drainage events.

The excess kurtosis distribution (Fig. 3b) serves as an indicator of rare events. Regions with large positive kurtosis exhibit strongly intermittent behavior rather than smooth, continuous variation. The elevated kurtosis observed near the right gas channel suggests that water intrusion may occur in short, localized bursts driven by transient pressure conditions or shifts in capillary balance, rather than through steady, continuous inflow.

The Shannon entropy (Fig. 3c) reflects the diversity of saturation states. High-entropy regions indicate that water is not pinned to a single saturation level. Conversely, low-entropy regions remain within a restricted range of perfectly dry or wet states.

FFT analysis (Fig. 5) of the fluctuation component, except for the DC component, provides complementary information regarding the distribution of dominant periodic time-scales in water fluctuations. Regions where low-frequency components dominate are associated with slower accumulation–release processes, whereas regions with stronger high-frequency components exhibit more rapid fluctuation behavior. In the vicinity of the channel–rib transition, relatively stronger low-frequency components were observed, suggesting that slower accumulation–drainage dynamics contribute in these areas.

Temporal autocorrelation characterizes the persistence of water states (Fig. 6). High autocorrelation indicates that saturation conditions are maintained over extended time intervals, suggesting slow evolution and possible stabilization by capillary retention or structural confinement. Low autocorrelation, in contrast, implies frequent state changes and more dynamically evolving water transport.

Hereafter, the discussion combines these indicators. The dominance of low-frequency components (Fig. 5) may appear inconsistent with the event-dominated behavior implied by elevated excess kurtosis (Fig. 3b). This apparent discrepancy could explain that water dynamics involve a multi-scale structure in which processes with different intermittency and temporal characteristics coexist. Intermittent capillary-driven water liquid drainage and fluctuation-induced evaporation/condensation processes associated with external variations in temperature, pressure, or gas flow rate may simultaneously contribute to the observed behavior. In such cases, composite indicators may contribute to the in-depth understanding of coexisting water dynamics regimes. As an example, the combination of variability and persistence is presented in Fig. 9. As discussed in the previous subsection, this pairing was selected to demonstrate a statistically robust combination for comparative classification, not representing a universal water dynamics mechanism. Regions characterized by high variability with low persistence, observed particularly beneath the land, correspond to areas, where water transitions among multiple states with rapid replacement, suggesting a relatively mobile liquid water regime. In contrast, regions exhibiting both high variability and high persistence represent distinctive water dynamics behavior in which water fluctuates within a limited range without being completely expelled, consistent with a dynamically maintained balance that may involve repeated condensation and evaporation. Low variability combined with low persistence corresponds to persistently dry conditions, whereas low variability with high persistence indicates that a relatively constant amount of liquid water remains pinned over time. This two-dimensional classification framework enables water distributions to be positioned not by a single scalar quantity but within a space defined by water state diversity.

Conclusions

In this study, we established a spatiotemporal analytical framework based on a dataset of operando X-ray radiography to characterise the water dynamics within the gas diffusion layer in polymer electrolyte fuel cells. Moving beyond conventional evaluation methods based on spatiotemporal averaging, we applied time series-based indicators, including skewness, excess kurtosis, Shannon entropy, frequency and temporal autocorrelation, to demonstrate that quasi-steady-state water distributions exhibit hidden features with dynamically maintained equilibrium states.

The results demonstrate that even when the average saturation appears stationary, the underlying behaviour of water is governed by spatially heterogeneous and multi-scale temporal dynamics. These findings indicate that water in the gas diffusion layer should not be regarded as a static distribution, but rather as a spatially structured dynamic system characterised by intermittency, persistence, and distributed fluctuation time scales.

Although the current framework does not directly establish universal water transport mechanisms under practical operating conditions, it provides a new diagnostic perspective that cannot be captured by traditional mean-value analyses. We expect that this framework may support future comparative analyses of GDL structures and contribute to the further development of data-driven diagnostic methodologies for PEFC water management.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6an00365f.

The analysis scripts and processed dataset used in this study are available from the corresponding author upon reasonable request for academic research purposes.

References

1. W. Yoshimune, Cell Rep. Phys. Sci., 2026, 7, 103146.
2. K. Kodama, T. Nagai, A. Kuwaki, R. Jinnouchi and Y. Morimoto, Nat. Nanotechnol., 2021, 16, 140–147.
3. D. A. Cullen, K. C. Neyerlin, R. K. Ahluwalia, R. Mukundan, K. L. More, R. L. Borup, A. Z. Weber, D. J. Myers and A. Kusoglu, Nat. Energy, 2021, 6, 462–474.
4. Y. Chen, G. Stelmacovich, A. Mularczyk, D. Parkinson, S. K. Babu, A. Forner-Cuenca, S. Pylypenko and I. V. Zenyuk, ACS Catal., 2023, 13, 10010–10025.
5. J. T. Lang, D. Kulkarni, C. W. Foster, Y. Huang, M. A. Sepe, S. Shimpalee, D. Y. Parkinson and I. V. Zenyuk, Chem. Rev., 2023, 123, 9880–9914.
6. W. Yoshimune, Bull. Chem. Soc. Jpn., 2024, 97, uoae046.
7. Y. P. Wijaya, F. P. Orfino and E. Kjeang, ACS Appl. Energy Mater., 2025, 8, 12997–13019.
8. C. P. Liu, P. Saha, Y. Huang, S. Shimpalee, P. Satjaritanun and I. V. Zenyuk, ACS Appl. Mater. Interfaces, 2021, 13, 20002–20013.
9. S. Kato, S. Yamaguchi, A. Kato, Y. Matsuoka, Y. Nagai and T. Suzuki, J. Chem. Eng. Jpn., 2022, 55, 71–76.
10. Y. Nagai, J. Eller, T. Hatanaka, S. Yamaguchi, S. Kato, A. Kato, F. Marone, H. Xu and F. N. Büchi, J. Power Sources, 2019, 435, 226809.
11. S. Yamaguchi, S. Kato, A. Kato, Y. Matsuoka, Y. Nagai and T. Suzuki, Electrochem. Commun., 2021, 128, 107059.
12. J. Eller, F. Marone and F. N. Büchi, ECS Trans., 2015, 69, 523–531.
13. A. Colliard-Granero, K. Duan, R. Zeis, M. H. Eikerling, K. Malek and M. J. Eslamibidgoli, Digital Discovery, 2025, 4, 2724–2736.
14. A. Colliard-Granero, S. Ranieri, A. Bazylak, T. Morawietz, K. A. Friedrich, J. Jankovic, M. H. Eikerling, K. Malek and M. J. Eslamibidgoli, J. Electrochem. Soc., 2025, 172, 074515.
15. A. D. Shum, C. P. Liu, W. H. Lim, D. Y. Parkinson and I. V. Zenyuk, Transp. Porous Media, 2022, 144, 715–737.
16. R. Banerjee, N. Ge, C. Han, J. Lee, M. G. George, H. Liu, D. Muirhead, P. Shrestha and A. Bazylak, Int. J. Hydrogen Energy, 2018, 43, 9757–9769.
17. W. Yoshimune, A. Kato, T. Hayakawa, S. Yamaguchi and S. Kato, Adv. Energy Sustainability Res., 2024, 5, 2400126.
18. W. Yoshimune, A. Kato, T. Hayakawa, S. Yamaguchi and S. Kato, npj Mater. Degrad., 2024, 8, 106.
19. A. Kato, S. Kato, S. Yamaguchi, T. Suzuki and Y. Nagai, Int. J. Hydrogen Energy, 2024, 50, 1218–1227.
20. E. A. Chadwick, B. Derebaşi, V. P. Schulz and A. Bazylak, Sci. Rep., 2025, 15, 39380.
21. W. Yoshimune, A. Kato, T. Hayakawa, S. Yamaguchi and S. Kato, J. Power Sources, 2025, 655, 237896.
22. A. Kato, S. Kato, S. Yamaguchi, T. Suzuki and Y. Nagai, J. Power Sources, 2022, 521, 230951.
23. D. N. Joanes and C. A. Gill, The Statistician, 1998, 47, 183.
24. C. E. Shannon and W. Weaver, in A mathematical theory of communication, University Illinois Press, Urbana, 1949.
25. E. C. Pielou, J. Theor. Biol., 1966, 13, 131–144.
26. J. W. Cooley and J. W. Tukey, Math. Comput., 1965, 19, 297–301.
27. W. Press, S. Teukolsky, W. Vetterling and B. Flannery, Second-order conservative equations, in Numerical Recipes: The Art of Scientific Computing, 2007.
28. G. E. P. Box and G. M. Jenkins, Time Series Analysis: Forecasting and Control, Holden-Day, San Francisco, 1976.
29. S. Torquato, in Random Heterogeneous Materials: Microstructure and Macroscopic Properties, Springer, 2002.
30. M. Thomas and J. A. T. Cover, Elements of Information Theory, Wiley, 2005.
31. A. Kraskov, H. Stögbauer and P. Grassberger, Phys. Rev. E: Stat. Nonlin. Soft Matter. Phys., 2004, 69, 066138.

---