A physics-guided Fourier neural operator to predict wind-driven water film/rivulet flows pertinent to aircraft icing phenomena

Jincheng Wang; Charlelie Laurent; Suhas S. Jain; Hui Hu

doi:10.1017/jfm.2026.12103

A physics-guided Fourier neural operator to predict wind-driven water film/rivulet flows pertinent to aircraft icing phenomena

Published online by Cambridge University Press: 12 October 2026

and

Jincheng Wang: Affiliation:
Department of Aerospace Engineering, Iowa State University , Ames, IA 50011-1096, USA
Charlelie Laurent: Affiliation:
Center for Turbulence Research, Stanford University, Stanford, CA 94305, USANVIDIA Corporation, Santa Clara, CA 95051, USA
Suhas S. Jain: Affiliation:
Center for Turbulence Research, Stanford University, Stanford, CA 94305, USAFlow Physics and Computational Science Lab, Georgia Institute of Technology, Atlanta, GA 30318, USA
Hui Hu*: Affiliation:
Department of Aerospace Engineering, Iowa State University , Ames, IA 50011-1096, USA
*: Corresponding author: Hui Hu, huhui@iastate.edu

Article contents

Abstract
Introduction
Experimental set-up and machine learning framework
Results and discussion
Conclusion
Funding
Declaration of interests.
Data availability statement
References

Rights & Permissions

Abstract

Content of image described in text.

Aircraft icing poses a critical threat to aviation safety in cold climates, with the dynamic glaze ice accretion process governed by the transient behaviour of wind-driven runback water film and rivulet (WDRWF/R) flows on wing surfaces. Despite its significance, accurately modelling this complex multiphase flow phenomenon remains a fundamental challenge for conventional numerical and theoretical approaches. This study introduces a physics-guided Fourier neural operator (PGFNO) that achieves unprecedented accuracy in predicting the spatiotemporal evolution of WDRWF/R flows relevant to aircraft icing. To establish a high-fidelity training dataset, we conducted an extensive experimental campaign utilising a novel digital image projection (DIP) technique, which enables temporally and spatially resolved measurements of water film thickness distributions over a flat plate under diverse operating conditions. The PGFNO architecture incorporates a composite loss function that enforces physical constraints, ensuring both accurate field reconstruction and mass conservation. The trained model demonstrates robust predictive capability across a wide range of flow conditions, including extrapolation to unseen wind speeds, water flow rates and initial conditions, while accurately capturing both the film thickness field and its spectral characteristics. A comprehensive ablation study systematically evaluates the contributions of key architectural components to model performance, including the physics-guided loss function, embedding-layer configuration and Fourier mode/block selection. Long-term predictions via recurrent forecasting achieve errors approaching the experimental uncertainty of DIP measurements, while robustness tests with spatially distributed random noise confirm exceptional model stability. These results establish PGFNO as a powerful predictive framework for complex multiphase flows, with direct implications for aircraft icing prediction and mitigation strategies.

JFM classification

Interfacial Flows (free surface): Interfacial Flows (free surface)Interfacial Flows (free surface): Thin films Multiphase and Particle-laden Flows: Multiphase flow

Information

Type: JFM Papers
Information: Journal of Fluid Mechanics , Volume 1044 , 10 October 2026, A39
DOI: https://doi.org/10.1017/jfm.2026.12103 [Opens in a new window]
Copyright: © The Author(s), 2026. Published by Cambridge University Press

1. Introduction

Aircraft icing represents one of the most critical meteorological hazards to aviation safety in cold climates, posing substantial risks to both commercial and military operations (Gent et al. Reference Gent, Dart and Cansdale2000; Cebeci & Kafyeke Reference Cebeci and Kafyeke2003; Cao et al. Reference Cao, Tan and Wu2018). When aircraft traverse icing clouds, supercooled water droplets impact and freeze upon airframe surfaces, with the freezing characteristics determined by the rate at which latent heat of fusion dissipates into the ambient flow. This process manifests in two distinct regimes: dry icing conditions produce rime ice through complete freezing of impinged water, whereas wet icing conditions generate glaze ice, i.e. a considerably more hazardous phenomenon wherein only partial freezing occurs at the impingement zone, with remaining liquid water flowing downstream to freeze beyond the initial impact area. The irregular morphology of glaze ice, characterised by pronounced ice horns and extensive feather formations, severely distorts aerofoil profiles and dramatically degrades aerodynamic performance through simultaneous drag increase and lift reduction (Bragg et al. Reference Bragg, Gregorek and Lee1986; Lynch & Khodadoust Reference Lynch and Khodadoust2001; Wang et al. Reference Wang, Dhulipalla and Hu2026a ). Under moderate to severe icing conditions, ice accumulation can cause premature stall at substantially elevated speeds and reduced angles of attack, potentially rendering continued flight impossible (Heinrich et al. Reference Heinrich1991; Cao & Yuan Reference Cao and Yuan2007; Ma et al. Reference Ma, Lin, Jin, Bu, Shen and Wen2025). The catastrophic consequences of inadequate ice protection are tragically illustrated by incidents such as Continental Flight 3407, which crashed in Buffalo, New York, claiming 50 lives due to wing ice accumulation (Mosher et al. Reference Mosher, Schaum, Herbster and Guinn2010), while military operations similarly suffer mission delays, cancellations and abortions attributable to icing hazards (Ryerson et al. Reference Ryerson, Peck and Martel2003).

Despite significant advances in computational ice accretion prediction through established engineering packages such as IGLOO, TURBICE, NASA’s LEWICE and GlennICE, and ANSYS FENSAP-ICE (Wright Reference Wright2008; Ozcer et al. Reference Ozcer, Baruzzi, Reid, Habashi, Fossati and Croce2011; Villedieu et al. Reference Villedieu, Trontin, Guffond and Bobo2012, Reference Villedieu, Trontin and Chauvin2014; Wright et al. Reference Wright, Porter, Galloway and Rigby2022), fundamental limitations persist in accurately forecasting glaze ice formation. These widely used frameworks combine droplet impingement models, thermodynamic freezing formulations and surface water transport models to simulate ice accretion on aerodynamic surfaces under atmospheric icing conditions. While these tools demonstrate satisfactory performance for rime ice scenarios and incorporate well-established classical models for multiphase transport phenomena, they systematically fail to capture the intricate microphysical processes governing glaze ice accretion. In particular, rime icing is characterised by immediate freezing of impinging droplets, whereas glaze icing involves significant liquid water transport along the surface prior to freezing, which introduces additional hydrodynamic complexity. Critical deficiencies include oversimplified treatment of transient surface water transport dynamics and surface roughness effects, which are often modelled heuristically and decoupled from the strongly interactive aerodynamic and thermodynamic environment (Hansman & Turnock Reference Hansman and Turnock1989; Anderson & Tsao Reference Anderson and Tsao2003). In most engineering ice accretion tools, the transport of liquid water along the surface is represented using simplified thin-film formulations combined with empirical correlations for freezing and roughness effects. This incomplete representation of surface physics constitutes a principal source of discrepancy between computational predictions and experimental observations in glaze icing conditions (Feo Reference Feo2000; Myers Reference Myers2001). Recent investigations have established that accurate glaze ice prediction fundamentally requires precise characterisation of wind-driven runback water film and rivulet (WDRWF/R) flow behaviour (Waldman & Hu Reference Waldman and Hu2015). Under glaze icing conditions, impinging supercooled droplets accumulate on the surface and merge into a thin shear-driven liquid film that is transported downstream by aerodynamic shear prior to freezing. The thickness of this film is typically of the order of tens to hundreds of micrometres, and the liquid layer often develops interfacial waves, rivulet structures and intermittent breakup, leading to highly heterogeneous liquid transport along the surface. The dynamics of such air-driven thin films interacting with aerodynamic boundary layers has been investigated through theoretical and numerical studies (Rothmayer et al. Reference Rothmayer, Matheis and Timoshin2002; Wang & Rothmayer Reference Wang and Rothmayer2009; Wang et al. Reference Wang, He and Hu2026b ), which demonstrate that the film transport behaviour is governed by a delicate balance between aerodynamic shear stress, viscous forces, surface tension and gravitational effects. A central difficulty in modelling WDRWF/R flows lies in accurately predicting the wetting-front or contact-line motion of the liquid film over the solid surface. In thin-film formulations, the motion of this contact line determines the spatial extent of the wetted region and therefore directly controls the downstream redistribution of liquid mass. However, the hydrodynamics near a moving contact line is inherently singular under the classical no-slip boundary condition, requiring additional modelling assumptions such as slip-length regularisation, precursor-film models or dynamic contact-angle relations (Huh & Scriven Reference Huh and Scriven1971; de Gennes Reference de Gennes1985; Cox Reference Cox1986). In practical icing environments, these modelling challenges are further amplified by microscale surface roughness and contact-angle hysteresis, which introduce pinning and depinning events that strongly influence the motion of the wetting front (Joanny & de Gennes Reference Joanny and de Gennes1984; Quéré Reference Quéré2008). Consequently, even when the bulk thin-film transport is captured reasonably well, numerical simulations often struggle to reproduce the experimentally observed location and evolution of the wetting front. The influence of roughness on air-driven water films has been examined in several studies investigating the interaction between thin liquid films and subgrid surface roughness. Wang & Rothmayer (Reference Wang and Rothmayer2009) demonstrated that unresolved roughness elements can significantly modify the transport behaviour of thin films driven by aerodynamic shear by altering the effective hydrodynamic resistance and water redistribution within the film. Subsequent work proposed low-order modelling approaches to describe the interaction between air-driven water films and unresolved roughness structures on aerodynamic surfaces (Rothmayer & Hu Reference Rothmayer and Hu2015). These studies suggest that even relatively small roughness elements can substantially influence the effective mobility of the liquid layer and introduce localised water trapping mechanisms that strongly affect the downstream evolution of runback flows. In addition, roughness elements may modify the external boundary-layer structure and the wall shear stress distribution acting on the liquid film, thereby introducing additional coupling between aerodynamic flow and surface water transport (Servini et al. Reference Servini, Smith and Rothmayer2018); it also enhances convective heat transfer at the interface, thereby modulating the local freezing rate and subsequent ice accretion. Although numerous numerical and experimental studies have attempted to simulate WDRWF/R flow dynamics (Myers & Thompson Reference Myers and Thompson1998; Myers et al. Reference Myers, Liang and Wetton2004; Zhang & Hu Reference Zhang and Hu2021; Chen et al. Reference Chen, Yi, Wang and Liu2023; Ma et al. Reference Ma, Lin, Jin, Bu, Shen and Wen2025; Wang et al. Reference Wang, Dhulipalla, Zhao, Nunes and Hu2025; Zhao et al. Reference Zhao, Wang and Hu2025), persistent challenges in capturing the complex multiphase physics (even excluding phase-change complications) continue to yield substantial disagreement between computational and experimental results. One important reason is the inherently three-dimensional and unsteady nature of runback transport, where thin surface films frequently transition into rivulet networks that transport water through localised streams separated by partially dry regions. These rivulet structures introduce complex three-dimensional contact-line dynamics as well as intermittent breakup and merging events that are difficult to capture using simplified thin-film models. Recent advances in experimental diagnostics have provided valuable insight into these processes. For example, Zhang et al. (Reference Zhang, Tian and Hu2015) employed a digital image projection (DIP) technique to obtain time-resolved measurements of surface water thickness on aerofoil surfaces, revealing the highly dynamic and spatially heterogeneous evolution of wind-driven runback flows.

Machine learning has emerged as a transformative paradigm for modelling complex physical systems through data-driven pattern recognition, offering compelling advantages over traditional physics-based approaches when confronted with phenomena that resist analytical or numerical treatment. The availability of high-fidelity experimental databases from recent icing tunnel campaigns provides an unprecedented opportunity to leverage machine learning for accurate, efficient prediction of aircraft icing phenomena. Several pioneering studies have explored neural network applications in icing analysis, predominantly focusing on ice shape prediction and severity assessment. Cao et al. (Reference Cao, Yuan and Li2011) developed neural network models to quantify ice horn geometric effects on aerofoil aerodynamic performance, while McCann (Reference McCann2005) predicted icing severity from cloud microphysical parameters. More recently, Yu et al. (Reference Yu, Han, Zhang, Zhang, Liu and Chen2022) introduced a multi-autoencoder fusion architecture for rapid ice profile prediction, and Wang et al. (Reference Wang, Hu, He and Hu2023a ) identified sparse governing equations for runback water dynamics using data-driven discovery techniques. However, these studies exhibit critical limitations: most employ machine learning to predict steady or periodic flows that inadequately represent actual icing transience, and most train models on synthetic data generated from simplified numerical simulations rather than physics-preserving real experimental measurements. Consequently, their predictive capability for real water runback film dynamics remains limited. Wang et al. (Reference Wang, Hu, He and Hu2023b ) made significant progress by proposing a ConvLSTM network for unsteady WDRWF/R prediction on flat plates, yet their framework’s capability of forecasting only 10 future time steps substantially constrains practical applicability.

Recent advances in scientific machine learning have introduced several neural-network frameworks for modelling systems governed by partial differential equations (PDEs). One widely used approach is the physics-informed neural network (PINN), which incorporates governing equations and boundary conditions directly into the training process by minimising the PDE residual together with available data (Raissi et al. Reference Raissi, Perdikaris and Karniadakis2019). Such PINNs have been applied to a range of forward and inverse PDE problems, particularly when labelled data are limited. However, they are typically trained for a specific configuration of the governing system and therefore generalise poorly across different operating conditions; in addition, residual-based optimisation can become challenging for strongly nonlinear or multiscale dynamics (Krishnapriyan et al. Reference Krishnapriyan, Gholami, Zhe, Kirby and Mahoney2021; Wang et al. Reference Wang, Teng and Perdikaris2021). Another related framework is the deep operator network (DeepONet), which learns nonlinear operators mapping between function spaces using a branch–trunk architecture (Lu et al. Reference Lu, Jin, Pang, Zhang and Karniadakis2021). Unlike PINNs, DeepONet learns a parametric solution operator and can therefore predict solutions for families of PDE problems under varying inputs. However, the approach can become computationally demanding for high-dimensional spatial fields and does not explicitly exploit spectral structures often present in fluid mechanics problems. More recently, neural operators have been proposed to directly learn mappings between function spaces using neural parametrisations of integral kernels. Among these approaches, the Fourier neural operator (FNO) represents the operator kernel in Fourier space and performs global spectral convolutions to learn the mapping between input and output fields (Li et al. Reference Li, Kovachki, Azizzadenesheli and Anandkumar2021; Kovachki et al. Reference Kovachki, Li, Liu and Anandkumar2023). Because the learned operator acts on entire fields rather than discrete grid values, neural operators exhibit generalisation across discretisations and parameter variations while efficiently capturing non-local spatial interactions and multiscale structures. One limitation of the original FNO formulation is that it is naturally defined on regular computational grids and is therefore less straightforward to apply to problems with complex or varying geometries. Recent developments, such as geometry-aware neural operators (Zhong & Meidani Reference Zhong and Meidani2025), extend the framework to irregular domains by incorporating geometric information into the operator representation. These properties make the FNO particularly suitable for modelling WDRWF/R flows, which exhibit strongly coupled spatial interactions and multiscale spatiotemporal dynamics associated with thin-film transport, rivulet formation and surface-wave evolution. In particular, the spectral convolution employed in the FNO captures non-local spatial interactions across the computational domain while efficiently representing multiscale flow structures spanning thin films, rivulets and droplet regimes. Moreover, the Fourier basis naturally decomposes the flow field into sinusoidal components, which provides an effective representation of the surface-wave dynamics commonly observed in experimental WDRWF/R flows.

In parallel with these developments, physics-guided neural networks have recently emerged as a powerful paradigm for improving both the physical realism and generalisation capability of data-driven flow prediction models. For example, Yousif et al. Reference Yousif, Yu and Lim2022) proposed a physics-guided deep learning framework for generating turbulent inflow conditions by combining a multiscale convolutional autoencoder with subpixel convolution and a long short-term memory model. By incorporating physical constraints and turbulence statistics into the training process, the model restricts the learned solution space to physically admissible states. This physics-based guidance not only enables accurate reproduction of instantaneous flow structures, turbulence statistics and spectral characteristics but also significantly improves model generalisation and transferability, allowing the trained network to adapt efficiently to new flow conditions or regimes with substantially reduced training data requirements.

Motivated by these recent advances, the present study develops a comprehensive physics-guided machine learning framework that synergistically integrates multiphase flow physics understanding, advanced neural operator architectures and high-fidelity experimental data to achieve accurate predictive modelling of WDRWF/R flow in aircraft icing contexts. We propose a physics-guided FNO (PGFNO), incorporating explicit physical constraints, to predict the spatiotemporal evolution of WDRWF/R flow across varying wind speeds and water flow rates on a flat plate. To generate a comprehensive training dataset spanning the relevant parameter space, systematic wind tunnel experiments were conducted using a novel DIP technique (Zhang et al. Reference Zhang, Tian and Hu2015) to measure three-dimensional spatiotemporal thickness fields with high precision. Such experimentally resolved datasets are particularly valuable because the coupled effects of aerodynamic forcing, thin-film transport and contact-line evolution are difficult to simulate simultaneously within conventional water runback modelling frameworks, which typically rely on simplified thin-film models (Rothmayer et al. Reference Rothmayer, Matheis and Timoshin2002; Wang & Rothmayer Reference Wang and Rothmayer2009). By learning from experimentally measured spatiotemporal WDRWF/R evolution, the proposed framework captures these strongly coupled dynamics that is otherwise difficult to resolve using existing modelling approaches. This physics-guided, data-driven approach addresses fundamental limitations of current methodologies by leveraging experimental observations while embedding known physical principles, thereby enabling accurate long-term prediction of the complex, transient behaviour of WDRWF/R flow pertinent to aircraft icing and providing insight into the underlying physical mechanisms.

The remainder of this paper proceeds as follows. Section 2 describes experimental methodology and PGFNO architectural design. Section 3 presents a comprehensive performance evaluation of the PGFNO model against experimental ground truth across multiple operating conditions. Section 4 concludes with a synthesis of key findings and implications for aircraft icing research and applications.

Figure 1.

(a) The principle of the DIP technique. (b) The experimental set-up for measuring WDRWF/R flow by the DIP technique.

A diagram illustrating the principle and experimental setup of the DIP technique for measuring water droplet flow.

Figure 1. Long description

Panel A: A schematic diagram showing the principle of the DIP technique. A high-speed camera and a pattern projector are positioned above a water rivulet. The high-speed camera captures images in the imaging direction, while the pattern projector projects patterns in the projecting direction. The water rivulet flows on a reference plane, forming a water head and a water film. Panel B: A diagram of the experimental setup for measuring wind-driven runback water film or rivulet flow using the DIP technique. The setup includes a high-speed camera, a pattern projector, a flat plate, a water tank, a digital water pump, and a water collector. Incoming airflow drives the water film or rivulet over the flat plate, which is connected to the water tank and pump system.

2. Experimental set-up and machine learning framework

2.1. Fundamentals of the DIP technique

The DIP technique is a non-intrusive, structured optical method for reconstructing three-dimensional features of objects, originally developed for measuring surface water transport over aerofoil surfaces (Zhang et al. Reference Zhang, Tian and Hu2015). It was subsequently optimised by Liu et al. (Reference Liu, Zhang, Tian and Hu2020) and Gao et al. (Reference Gao, Liu and Hu2020) to quantify the wind-driven runback of water–ice mixtures during the glaze icing process on aerofoil surfaces. As shown schematically in figure 1(a), the DIP technique projects patterns (e.g. maze patterns in figure 2 a) onto a reference plane using a digital projector while a high-speed camera records these patterns; when an object is placed on the reference plane, the projected patterns deform in a linear relationship with the object’s height (or thickness), expressed as

(2.1)

\begin{align} \frac{g}{q}=\frac{d-h}{h}, \end{align}

Figure 2.

(a) A schematic diagram of quantifying the thickness field of the WDRWF/R flow by the DIP technique. (b) The displacement-to-height conversion map of the DIP system calibrated for the present study. (c) Distribution of the experimental test conditions (i.e. the incoming wind speed, V _∞, and water flow rate, Q) used in this study.

Diagram and graphs depicting the study of water droplet flow and experimental conditions.

where h represents the object’s height (i.e. line $\overline{\textit{SO}}$ in figure 1 a), d is the vertical distance between the camera and the reference plane, q denotes the displacement of the image pattern (i.e. line $\overline{\textit{EF}}$ in figure 1 a) and g is the distance between the camera and the projector (i.e. line $\overline{\textit{AB}}$ in figure 1 a). This relationship can be rearranged and simplified to

(2.2)

\begin{align} h=\frac{d-h}{g}q\approx \frac{d}{g}q=k \cdot q, \end{align}

since the object’s height h is typically much smaller than the distance between the camera and the reference plane d. Here, $k=d/g$ denotes the displacement-to-height ratio, which remains constant for a fixed DIP system configuration. This parameter is determined through calibration by mounting a test plate on a vertical translation stage at predetermined heights. The recorded maze patterns are processed to acquire displacement vectors relative to a reference position (z = 0). The slope of the linear height–displacement relationship is then derived at each point on the plane, yielding the conversion factor for translating displacement to height measurements, which enables three-dimensional shape reconstruction. Figure 1(b) presents a schematic of the DIP system applied in the present WDRWF/R study. Figure 2(a) shows sample DIP maze patterns used in the experiments, where the wall and passage block sizes are 2 pixels and 3 pixels, respectively. Figure 2(b) presents the corresponding calibrated displacement-to-height conversion factor map (i.e. the k map). The maximum spatial variation of the factor k across the entire flat plate is about 0.12 %, indicating a quite uniform displacement-to-height distribution on the flat plate. This small variation in k is caused by the angle between the camera and projector, given the divergent nature of the light beam emitted from the projector. As an image-correlation-based technique, DIP provides spatially averaged measurements with spatial resolution primarily determined by the interrogation window size used in cross-correlation. The temporal resolution of the DIP technique is determined by the frame rate of the camera. The measurement uncertainty of the DIP system arises from several sources, including minor vibrations in the test section, the finite resolution of the imaging and projection hardware and the signal-to-noise ratio of the captured images (Zhang & Hu Reference Zhang and Hu2021). Following the analysis by Gendrich & Koochesfahani (Reference Gendrich and Koochesfahani1996), displacement measurements derived from cross-correlation algorithms exhibit an uncertainty of approximately ±0.2 pixels at the 95 % confidence level, assuming a Gaussian error distribution. When applied to the displacement-to-height conversion factor (i.e. k ≈ 88 µm px⁻¹ as shown in figure 2 b) employed in this study, this translates to a film-thickness measurement uncertainty of around ±20 μm. The measurement uncertainty of the DIP system was also assessed using 3D-printed hemispherical roughness elements mounted on a test plate, following the methodology described by Zhang et al. (Reference Zhang, Tian and Hu2015). Approximately 500 measurement points were obtained from hemispherical elements with a nominal height of 4.00 mm. Based on these measurements, the uncertainty of the DIP system was estimated to be about 20 μm, corresponding to roughly 0.5 % of the nominal height H = 4.00 mm. Further details on the theoretical foundation and implementation of the DIP technique can be found in Zhang et al. (Reference Zhang, Tian and Hu2015), Liu et al. (Reference Liu, Zhang, Tian and Hu2020) and Gao et al. (Reference Gao, Liu and Hu2020).

2.2. Experimental set-up for WDRWF/R flow

The experimental set-up for characterising WDRWF/R flows is illustrated in figure 1(b). A rectangular test plate, measuring 254 mm in length, 152 mm in width and 12 mm in thickness, was mounted horizontally at the base of the wind tunnel test section. The plate surface was first coated with a white hydrophilic protective enamel (Rustoleum™ Flat Protective Enamel) and then sequentially polished with fine-grit sandpapers up to 2000 grit to achieve a smooth and uniform finish. The WDRWF/R flows were generated by pumping water from a reservoir tank through a spanwise array of 13 precision-drilled holes (2 mm in diameter, spaced 5 mm apart, spanning a total width of D = 60 mm) located near the leading edge of the test plate. The water pump (Cole-Parmer 75211-30) employed in this study operates over a rotational speed range of 60–3600 rpm. To confine the water film and prevent its lateral spreading under conditions of low wind speed and high flow rate, two parallel side strips were installed on the test plate. These strips were aligned with the streamwise direction of the incoming airflow and positioned 60 mm apart. To facilitate DIP measurements, a digital light projector (DLP^® LightCrafter™) projected a predefined maze pattern onto the plate surface, while a high-speed CCD camera (Photron FASTCAM Mini WX™, 2048 × 2048 pixels) recorded the resulting deformations of the projected pattern as the WDRWF/R flow advanced. The camera, equipped with a 55 mm Micro-Nikkor f/5.6 lens (Nikon), captured an image of a 100 mm × 48 mm observation area (2048 × 960 pixels), where a 6 mm wide region adjacent to each side boundary was excluded to eliminate the side-strip boundary effect. Water film thickness was determined using a cross-correlation algorithm with a 64 × 64 pixel interrogation window, yielding an effective spatial resolution of approximately 1.6 mm in the x–y plane. Figure 2(a) presents the raw DIP images (i.e. the reference pattern and the pattern distorted by the water film) alongside the reconstructed water film thickness field obtained from the cross-correlation analysis. For each experimental condition, the camera recorded 1250 consecutive frames at a frame rate of 125 Hz, corresponding to 10 s of continuous acquisition.

To enhance the optical contrast and image quality for DIP measurements, the working fluid was prepared by adding a flat white latex dye to deionised water at a volume fraction of 0.5 %. This dye concentration was selected to optimise light reflection at the free surface of the water film/rivulet flow while maintaining the fluid’s physical properties. The viscosity and surface tension measurements indicate that the dye addition induces a negligible change ( $\leqslant$ 1 %) in viscosity relative to that of pure water (1.0 mPa s), while reducing the surface tension from 0.072 to 0.070 N m⁻¹ (Zhang Reference Zhang2021), ensuring that the hydrodynamic behaviour of the working fluid remained representative of water-based flows. To characterise the wettability of the dyed working fluid on the test surface, the contact angles were measured using a shadowgraph-based imaging technique on the flat test plate, yielding a static contact angle of $\theta _{s}$ = 65° and an advancing contact angle of $\theta _{\alpha }$ ≈ 80°. These fluid properties were consistently maintained throughout all experiments to ensure reproducibility and comparability of results.

The experimental matrix encompassed systematic variations of two governing parameters: the free-stream velocity V _∞ (bulk airflow velocity), which ranged from 5 to 26 m s⁻¹, and the water flow rate Q, spanning 50–800 ml min⁻¹. To achieve comprehensive parametric coverage while optimising experimental efficiency, 79 distinct combinations of (V _∞, Q) were strategically selected using Latin hypercube sampling, as illustrated in figure 2(c). This sampling strategy ensures uniform exploration of the parameter space and enables robust characterisation of the coupled wind–water interaction dynamics across the operational envelope.

2.3. Machine learning framework for WDRWF/R flow

A FNO architecture with modifications to incorporate physical constraints is employed, forming a PGFNO model for predicting WDRWF/R flow dynamics. The PGFNO model establishes a nonlinear mapping between an input functional space and an output functional space. The input comprises a user-defined number of consecutive temporal snapshots of thickness fields discretised on a regular spatial grid with a temporal resolution of $\Delta t =8$ ms, along with the free-stream air velocity V _∞ and volumetric water flow rate Q. The output consists of a corresponding sequence of thickness-field snapshots that temporally succeed the input horizon, maintaining identical spatial discretisation and temporal resolution. For brevity, these temporal samplings are referred to as ‘snapshots’ throughout this work.

The PGFNO model can be interpreted as an expected-value estimator $\mathrm{\mathbb{E}}$ that predicts the conditional mean of future thickness-field realisations given past observations and flow conditions. At a high level, the PGFNO model $F_{\theta }$ , parametrised by θ, aims to approximate the mapping

(2.3)

\begin{align} F_{\theta } \big(\boldsymbol{P}_{i\leq {N_{{in}}}},V_{\infty },Q,X,Y,t \big)\approx \mathrm{\mathbb{E}} \big(\boldsymbol{P}_{i\gt {N_{{in}}}}| \boldsymbol{P}_{i\leq {N_{{in}}}},V_{\infty },Q,X,Y,t \big) ,\end{align}

where $\boldsymbol{P}_{i\leq {N_{{in}}}}$ denotes the sequence of input snapshots and $\boldsymbol{P}_{i\gt {N_{{in}}}}$ denotes the subsequent output snapshots. Figure 3 illustrates the schematic architecture of the PGFNO model, depicting the input–output structure and internal components in detail. The cubic representations in figure 3 denote the last three dimensions of the corresponding tensors. All trainable parameters within the PGFNO model are collectively denoted by θ.

Figure 3.

The PGFNO architecture for predicting WDRWF/R flow. Here $\boldsymbol{P}_i$ represents a snapshot of the water film thickness field at time step i. The three-dimensional block representations illustrate the last three dimensions of each intermediate tensor. The dimensions are evaluated using the default values employed in this study.

A diagram illustrating the PGFNO architecture for predicting WDRWF/R flow.

Figure 3. Long description

Panel A: The input sequence includes various parameters such as Pi, Ni, Y, t, V, and Q. The diagram shows the initial stages of the architecture, starting with an input layer and progressing through a series of transformations. Panel B: The Fourier block is detailed, showing the process from FFT to complex linear transform, inverse FFT, and convolution layers. The diagram includes multiple layers such as linear layers, complex-valued layers, inverse FFT layers, convolutional layers, and summation layers. The output sequence is represented as a series of transformed data.

The detailed operations involved in the PGFNO model are formulated as follows. Define $\boldsymbol{P}_{\kern-1pt j}\in \mathrm{\mathbb{R}}^{W\times L}$ as the jth thickness snapshot of WDRWF/R flow with width W = 30 and length L = 45, where j is a positive integer indicating the time index of the thickness snapshot with $\Delta t =8$ ms. Noteworthy here is that 1/4 of the upstream thickness fields are cropped due to limited GPU memory, given that we aim to build a PGFNO model with high complexity, allowing for the learning of complex WDRWF/R flow features. Consider a three-dimensional basic feature tensor $\boldsymbol{B}\in \mathrm{\mathbb{R}}^{W\times L\times D}$ with a feature dimension $D=N_{{in}}+D_{1}$ . The feature tensor consists of $N_{{in}}$ consecutive snapshots of the thickness field and $D_{1}$ = 5 additional fields. These additional fields include the spatial coordinates (i.e. x and y), the time coordinate t associated with the last snapshot, V _∞ and Q. To ensure compatibility with the spatial dimensions of the thickness snapshot, the scalar quantities t, V _∞ and Q are broadcast into matrices of dimensions W × L:

(2.4)

\begin{align} \begin{array}{lll} \boldsymbol{t}=t\boldsymbol{\cdot }\boldsymbol{J}_{W\times L}, & \boldsymbol{V}_{\infty }=V_{\infty }\boldsymbol{\cdot }\boldsymbol{J}_{W\times L}, & \boldsymbol{Q}=Q\boldsymbol{\cdot }\boldsymbol{J}_{W\times L} \end{array}, \end{align}

where $\boldsymbol{J}_{W\times L}$ is a matrix of ones with dimensions W × L. Thus, the basic feature tensor $\boldsymbol{B}\in \mathrm{\mathbb{R}}^{W\times L\times D}$ can be expressed as

(2.5)

\begin{align} \boldsymbol{B}= \big[\boldsymbol{P}_{1},\boldsymbol{P}_{2},\ldots ,\boldsymbol{P}_{{N_{{in}}}},\boldsymbol{X},\boldsymbol{Y},\boldsymbol{t},\boldsymbol{V}_{\infty },\boldsymbol{Q} \big], \end{align}

where $\boldsymbol{X}\in \mathrm{\mathbb{R}}^{W\times L}$ and $\boldsymbol{Y}\in \mathrm{\mathbb{R}}^{W\times L}$ are the spatial coordinate matrices. To predict $N_{out}$ sequential thickness fields in the time horizon at one time, PGFNO requires $N_{out}$ initial conditions. Therefore, the basic feature tensor is replicated $N_{out}$ times to form the input tensor $\boldsymbol{I}\in \mathrm{\mathbb{R}}^{{N_{{out}}}\times W\times L\times D}$ :

(2.6)

StartLayout 1st Row bold italic upper I equals left bracket ModifyingBelow bold italic upper B comma bold italic upper B comma ellipsis comma bold italic upper B With bottom brace Underscript upper N Subscript o u t Baseline t i m e s Endscripts right bracket period EndLayout

\begin{align} \boldsymbol{I}=\left[\underset{N_{{out}}\;{times}}{\underbrace{\boldsymbol{B},\boldsymbol{B},\ldots ,\boldsymbol{B}}}\right] \!. \end{align}

The PGFNO model then transforms the input tensor $\boldsymbol{I}$ through a sequence of operations. The initial linear transformation $\Omega \colon \mathrm{\mathbb{R}}^{{N_{{out}}}\times W\times L\times D}\rightarrow \mathrm{\mathbb{R}}^{{N_{{out}}}\times W\times L\times D'}$ lifts the input to a higher-dimensional space:

(2.7)

\begin{align} \boldsymbol{D}=\Omega (\boldsymbol{I}), \end{align}

where $D'$ = 40 is called the number of embedding layers. The lifted tensor $\boldsymbol{D}\in \mathrm{\mathbb{R}}^{{N_{{out}}}\times W\times L\times D'}$ is then transposed to be $\boldsymbol{H}^{(0)}\in \mathrm{\mathbb{R}}^{D'\times W\times L\times k}$ , which is then processed through a sequence of N _b Fourier blocks, where each Fourier block $\Phi _{m}\colon \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}\rightarrow \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}$ implements a dual-path architecture, where $m\in \{1,\ldots ,N\}.$ For the mth Fourier block, the transformation can be expressed as

(2.8)

\begin{align} \boldsymbol{H}^{\left(m\right)}=\mathrm{GELU} \big(\Phi _{m} \big(\boldsymbol{H}^{(m-1)} \big) \big)\in \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}, \\[-28pt] \nonumber \end{align}

(2.9)

\begin{align} \Phi \big(\boldsymbol{H}^{(m-1)} \big) = \mathrm{R}_{1} \big(\boldsymbol{H}^{(m-1)} \big)+\mathrm{R}_{2} \big(\boldsymbol{H}^{(m-1)} \big)\in \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}, \\[0pt] \nonumber \end{align}

where $\mathrm{GELU}\colon \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}\rightarrow \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}$ is a nonlinear activation function and the spectral branch $R_{1}\colon \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}\rightarrow \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}$ performs frequency-domain processing:

(2.10)

\begin{align} R_{1} \big(\boldsymbol{H}^{(m-1)} \big) = \big(\mathrm{\Upsilon}_{1}\circ \mathrm{\Upsilon}_{1}\circ \Xi \circ \Gamma \circ \Psi \big) \big(\boldsymbol{H}^{(m-1)} \big), \end{align}

where $\Psi \colon \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}\rightarrow \mathrm{\mathbb{C}}^{D'\times W\times L\times ({N_{{out}}}/2+1)}$ denotes the fast Fourier transform (FFT) operation applied along the spatial and temporal dimensions, $\Gamma \colon \mathrm{\mathbb{C}}^{D'\times W\times L\times ({N_{{out}}}/2+1)}\rightarrow \mathrm{\mathbb{C}}^{D'\times W\times L\times ({N_{{out}}}/2+1)}$ represents a learnable complex-valued linear transformation in the spatiotemporal frequency domain acting on a learnable spatial–temporal Fourier mode $\xi \in \mathrm{\mathbb{C}}^{K \times K \times K } ,\ \Xi \colon \mathrm{\mathbb{C}}^{D'\times W\times L\times ({N_{{out}}}/2+1)}\rightarrow \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}$ denotes the inverse FFT operation and $\mathrm{\Upsilon}_{1}\colon \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}\rightarrow \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}$ represents a convolutional operation. The second path $R_{2}$ consists of a convolution operation $\mathrm{\Upsilon}_{2}\colon \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}\rightarrow \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}$ :

(2.11)

\begin{align} R_{2} \big(\boldsymbol{H}^{(m-1)} \big) = \mathrm{\Upsilon}\big(\boldsymbol{H}^{(m-1)}\big). \end{align}

The features $\boldsymbol{H}^{(N)}$ from the final Fourier block $\Phi _{N}$ are processed by two convolutional operations ( $\mathrm{\Upsilon}_{3}\colon \mathrm{\mathbb{R}}^{D'\times W\times L\times {N_{{out}}}}\rightarrow \mathrm{\mathbb{R}}^{4D'\times W\times L\times {N_{{out}}}}$ and $\mathrm{\Upsilon}_{4}\colon \mathrm{\mathbb{R}}^{4D'\times W\times L\times {N_{{out}}}}\rightarrow \mathrm{\mathbb{R}}^{W\times L\times {N_{{out}}}}$ ) to generate the output:

(2.12)

\begin{align} \boldsymbol{S}=(\mathrm{\Upsilon}_{4}\circ \mathrm{\Upsilon}_{3}) \big(\boldsymbol{H}^{\left(N\right)} \big), \end{align}

where $\boldsymbol{S}\in \mathrm{\mathbb{R}}^{W\times L\times {N_{{out}}}}$ represents the predicted sequence of thickness fields, i.e. the subsequent $N_{{out}}$ future time steps of $\boldsymbol{P}_{\kern-1pt j}\in \mathrm{\mathbb{R}}^{W\times L}$ .The loss function comprises two components: the data-driven loss (L _data) and the physics-guided loss (L _physics). The data-driven loss term is formulated as the pointwise discrepancy between the predicted thickness values and the ground-truth measurements:

(2.13)

StartLayout 1st Row upper L Subscript d a t a Baseline equals sigma summation Underscript i equals upper N Subscript i n Baseline plus 1 Overscript upper N Subscript i n Baseline plus upper N Subscript o u t Baseline Endscripts StartRoot sigma summation StartAbsoluteValue StartFraction ModifyingAbove bold italic upper P With caret Subscript i Baseline minus bold italic upper P Subscript i Baseline Over normal m normal e normal a normal n left parenthesis ModifyingAbove bold italic upper P With caret Subscript i Baseline right parenthesis EndFraction EndAbsoluteValue EndRoot comma EndLayout

\begin{align} L_{data}= \mathop{\sum }\limits_{i=N_{in}+1}^{N_{in}+N_{{out}}}\sqrt{\sum \left| \frac{\hat{\boldsymbol{P}}_{i}-\boldsymbol{P}_{i}}{\mathrm{mean}\left(\hat{\boldsymbol{P}}_{i}\right)}\right| }, \end{align}

where $\hat{\boldsymbol{P}}_{i}$ represents the ground truth of the thickness snapshot to be predicted. The physics loss L_physics is composed of mass conservation loss L _phy, _mass and the mass flow rate conservation loss L _phy, _{mass flow rate} :

(2.14)

StartLayout 1st Row 1st Column upper L Subscript italic physics 2nd Column equals alpha bold dot upper L Subscript p h y comma italic mass Baseline plus beta bold dot upper L Subscript p h y comma italic mass italic flow italic rate Baseline 2nd Row 1st Column Blank 2nd Column equals alpha bold dot ModifyingBelow sigma summation Underscript i equals upper N Subscript i n Baseline plus 1 Overscript upper N Subscript i n Baseline plus upper N Subscript o u t Baseline Endscripts StartAbsoluteValue StartFraction sigma summation ModifyingAbove bold italic upper P With caret Subscript i Baseline minus sigma summation bold italic upper P Subscript i Baseline Over sigma summation ModifyingAbove bold italic upper P With caret Subscript i Baseline EndFraction EndAbsoluteValue With bottom brace Underscript normal m normal a normal s normal s normal c normal o normal n normal s normal e normal r normal v normal a normal t normal i normal o normal n normal l normal o normal s normal s Endscripts 3rd Row 1st Column Blank 2nd Column plus beta bold dot ModifyingBelow sigma summation Underscript i equals upper N Subscript i n Baseline plus 1 Overscript upper N Subscript i n Baseline plus upper N Subscript o u t Baseline minus 1 Endscripts StartAbsoluteValue StartFraction left parenthesis sigma summation ModifyingAbove bold italic upper P With caret Subscript i plus 1 Baseline minus sigma summation ModifyingAbove bold italic upper P With caret Subscript i Baseline right parenthesis minus left parenthesis sigma summation bold italic upper P Subscript i plus 1 Baseline minus sigma summation bold italic upper P Subscript i Baseline right parenthesis Over sigma summation ModifyingAbove bold italic upper P With caret Subscript i plus 1 Baseline minus sigma summation ModifyingAbove bold italic upper P With caret Subscript i Baseline EndFraction EndAbsoluteValue With bottom brace Underscript normal m normal a normal s normal s flow normal r normal a normal t normal e normal c normal o normal n normal s normal e normal r normal v normal a normal t normal i normal o normal n normal l normal o normal s normal s Endscripts period EndLayout

\begin{align} L_{\textit{physics}} & =\alpha \boldsymbol{\cdot }L_{{phy},\textit{mass}}+\beta \boldsymbol{\cdot }L_{{phy},\textit{mass}\,\textit{flow}\,\textit{rate}} \nonumber\\ & =\alpha \boldsymbol{\cdot }\underset{\mathrm{mass\,conservation\,loss}}{\underbrace{ \mathop{\sum }\limits_{i=N_{in}+1}^{N_{in}+N_{{out}}}\left| \frac{\sum \hat{\boldsymbol{P}}_{i}-\sum \boldsymbol{P}_{i}}{\sum \hat{\boldsymbol{P}}_{i}}\right| }}\nonumber\\ & \quad +\beta \boldsymbol{\cdot }\underset{\mathrm{mass}\,\text{flow}\,\mathrm{rate\,conservation\,loss}}{\underbrace{\mathop{\sum }\limits_{i=N_{in}+1}^{N_{in}+N_{{out}}-1}\left| \frac{\left(\sum \hat{\boldsymbol{P}}_{i+1}-\sum \hat{\boldsymbol{P}}_{i}\right)-\left(\sum \boldsymbol{P}_{i+1}-\sum \boldsymbol{P}_{i}\right)}{\sum \hat{\boldsymbol{P}}_{i+1}-\sum \hat{\boldsymbol{P}}_{i}}\right| }}. \end{align}

The mass conservation loss quantifies the discrepancy in integrated water volume between the predicted and ground-truth fields. The mass flow rate conservation loss represents the temporal derivative of this integrated water volume discrepancy. A comprehensive grid search was employed to systematically evaluate the influence of the hyperparameters α and β on model performance. The optimal configuration was found at α = 10⁻² and β = 10⁻⁴, yielding the best performance for the PGFNO model.

The PGFNO model was implemented in PyTorch. The default prediction problem settings are $N_{{in}}$ = 10 and $N_{{out}}$ = 50. Default key architectural hyperparameters include K = 14 Fourier modes, $D'$ = 40 embedding layers, N _b = 8 Fourier blocks, which were optimised through ablation tests discussed in § 3 and by considering the 80 GB memory constraint of the NVIDIA A100 GPU used in the current study. The total number of parameters in the PGFNO model used in this study, with the above default settings, is 281 032 321. Training was conducted over 1000 epochs with an initial learning rate of 0.005 and a batch size of 64. The Adam optimiser was employed in conjunction with a cosine annealing learning rate scheduler. To mitigate overfitting, L ₂ regularisation with a weight decay coefficient of 10⁻⁴ was applied. Data augmentation was implemented through horizontal reflection of the thickness field about the midspan (i.e. the x-axis arrow at the bottom part of figure 1 b) to enhance the dataset size and improve model generalisation.

2.4. Data splitting method

As mentioned in § 2.2, the dataset consists of $N_{c}$ = 79 cases, each representing a unique combination of design variables (V _∞ and Q). For each case, there are $N_{\textit{total}}$ = 1250 snapshots of the thickness field from the WDRWF/R flow dataset. To focus on the transient development of runback water, a time range $\chi _{\textit{start}}=\{1,2,\ldots ,100\}\subset \{1,2,\ldots ,N_{\textit{total}}\}$ for selecting starting snapshots is defined. For each starting time index $t_{\textit{start}}\in \chi _{\textit{start}}$ , the corresponding end index of the input–output pair is determined by $t_{end}=t_{\textit{start}}+N_{in}+N_{out}-1$ . This ensures that each input–output pair spans a total of $N_{in}+N_{out}$ = 60 consecutive snapshots, with the input sequence covering time range $[t_{\textit{start}},t_{\textit{start}}+N_{in}-1]$ and the output sequence covering time range $[t_{\textit{start}}+N_{in},t_{end}]$ . The total number of unique experimental snapshots utilised in this study is 12 640.

To evaluate the generalisation capability of the PGFNO model, the WDRWF/R flow dataset was partitioned using two distinct strategies: pairwise splitting and casewise splitting. Pairwise splitting assesses the model’s ability to generalise across different initial time states, whereas casewise splitting evaluates generalisation performance across different test configurations (V _∞ and Q).

In the pairwise splitting approach, we randomly sample a set of $N_{t}$ = 20 starting time indices $\mathrm{T }_{i}=\{t_{i,k}{\}}_{k=0}^{N_{t}}\subset \chi _{\textit{start}}$ for each case $i\in \{1,\ldots ,N_{c}\}$ . For each starting time index $t_{i,k}\in \mathrm{T }_{i}$ , an input–output pair is constructed:

(2.15)

\begin{align} \boldsymbol{X}_{i,k}=\left\{{\left\{{\boldsymbol{P}}_{j+t_{i,k}}^{\left(i\right)}\right\}}_{j=0}^{N_{in}-1},\boldsymbol{X,Y,}\boldsymbol{t}_{i,k}\boldsymbol{,}{\boldsymbol{V}}_{\infty }^{\left(i\right)}\boldsymbol{,}\boldsymbol{Q}^{\left(i\right)}\right\} \!, \\[-28pt] \nonumber \end{align}

(2.16)

\begin{align} \boldsymbol{Y}_{i,k}= \left\{{P}_{N_{in}+j+t_{i,k}}^{\left(i\right)} \right\}_{j=0}^{N_{out}-1} \!, \\[0pt] \nonumber \end{align}

where ${\boldsymbol{P}_{m}}^{(i)}\in \mathrm{\mathbb{R}}^{W\times L}$ represents the thickness snapshot at time step m for case i. The complete set of input–output pairs is $C=\{(\boldsymbol{X}_{i,k},\boldsymbol{Y}_{i,k})| i\in \{1,\ldots ,N_{c}\},k\in \{1,\ldots ,N_{t}\}\}$ with the total number of pairs of $| C| =N_{c}\times N_{t}=79\times 20=1580.$ These pairs are randomly split into training (90 %) and test sets (10 %):

(2.17)

\begin{align} C_{{train}}\subset C,\left| C_{{train}}\right| =\left\lfloor 0.9\times \left| C\right| \right\rfloor =1422, \\[-28pt] \nonumber \end{align}

(2.18)

\begin{align} C_{{test}}=C\backslash C_{{train}},\left| C_{{test}}\right| =158, \\[0pt] \nonumber \end{align}

where $\lfloor \boldsymbol{\cdot }\rfloor$ represents the floor function, which returns the greatest integer less than or equal to a given real number. This splitting strategy ensures that both the training and test sets contain samples from all 79 cases, but with different temporal starting points, allowing for the evaluation of the model’s ability to generalise across different phases of the flow development process.

However, in the casewise splitting approach, we first split the parameter space by randomly dividing the cases into training (90 %) and test (10 %) subsets:

(2.19)

\begin{align} I_{{train}}\subset \{1,\ldots ,N_{c}\},\left| I_{{train}}\right| =\left\lfloor 0.9\times N_{c}\right\rfloor =71, \\[-28pt] \nonumber \end{align}

(2.20)

\begin{align} I_{{test}}=\{1,\ldots ,N_{c}\}\backslash I_{{train}},\left| I_{{test}}\right| =8, \\[0pt] \nonumber \end{align}

where I _train and I _test are index sets identifying which cases belong to the training and test subsets. For each case in the respective subsets, we construct input–output pairs following the same methodology as in pairwise splitting:

(2.21)

\begin{align} G_{{train}}= \big\{(\boldsymbol{X}_{i,k},\boldsymbol{Y}_{i,k})| i\in I_{{train}},k\in \{1,\ldots ,N_{t}\} \big\}, \\[-28pt] \nonumber \end{align}

(2.22)

\begin{align} G_{{test}}= \big\{(\boldsymbol{X}_{i,k},\boldsymbol{Y}_{i,k})| i\in I_{{test}},k\in \{1,\ldots ,N_{t}\} \big\}. \\[0pt] \nonumber \end{align}

This approach evaluates the ability of the model to extrapolate to previously unseen physical parameters (V _∞ and Q), as the test set contains entirely different flow conditions from those present in the training set.

3. Results and discussion

3.1. Evaluation of PGFNO generalisability across different initial states

Figure 4.

Selected thickness-field snapshots predicted by the PGFNO model for the test case (V _∞ = 13 m s⁻¹, Q = 327 ml min⁻¹) under the pairwise splitting scheme at three time instances: (a) time = 0.032 s; (b) time = 0.24 s; (c) time = 0.40 s. Each row presents (left to right): ground-truth thickness field, PGFNO prediction and absolute error between prediction and ground truth.

Three sets of 3D surface plots compare ground truth, predictions, and errors over time.

Figure 4 presents the performance of the PGFNO model trained using the pairwise splitting methodology. Representative results are shown for a test case with a free-stream velocity of V _∞ = 13 m s⁻¹ and a flow rate of Q = 327 ml min⁻¹. Each row in the figure corresponds to a distinct temporal state, thereby visualising the temporal evolution of the WDRWF/R flow field. Absolute error maps are presented to provide a quantitative assessment of model accuracy.

The results demonstrate that the PGFNO model achieves high predictive accuracy across varying initial states, confirming its generalisability with respect to temporal initialisation. Across all 158 input–output pairs in the test set, the average L ₁ norm error is 4 μm, which falls well below the measurement uncertainty inherent to the experimental DIP system (ε ≈ 20 μm). This performance metric indicates that the model’s predictions are within the bounds of experimental precision, demonstrating not only temporal generalisability but also practical reliability for reconstructing the spatiotemporal dynamics of the WDRWF/R flow field.

3.2. Evaluation of PGFNO generalisability across different test conditions

Figures 5 and 6 present the performance of the PGFNO model trained using the casewise splitting methodology. The representative test cases shown correspond to two different flow regimes: a film flow case (V _∞ = 5 m s⁻¹ and Q = 196 ml min⁻¹; figure 5) and a rivulet flow case (V _∞ = 25 m s⁻¹ and Q = 527 ml min⁻¹; figure 6). Notably, the thickness field in figure 6 is thinner than that in figure 5, owing to the increased wind shear imposed by the elevated free-stream velocity. The figure layout follows the same convention as figure 4 to facilitate direct comparison.

Figure 5.

Selected thickness-field snapshots predicted by the PGFNO model for the test case (V _∞ = 5 m s⁻¹, Q = 196 ml min⁻¹) under the casewise splitting scheme at three time instances: (a) t = 0.064 s; (b) t = 0.20 s; (c) t = 0.40 s. Each row presents (left to right): ground-truth thickness field, PGFNO prediction and absolute error between prediction and ground truth.

Three sets of 3D surface plots compare ground truth, predictions, and absolute errors of ice accretion over time.

Figure 5. Long description

Panel A: At t = 0.064 s, three 3D surface plots are shown. The first plot, labeled Ground Truth, displays a color gradient representing ice thickness in millimeters, with a scale ranging from 0 to 2.5 mm. The second plot, labeled Prediction, shows a similar color gradient. The third plot, labeled Absolute Error, uses a grayscale to indicate the error magnitude, ranging from 0 to 0.3 mm. Panel B: At t = 0.20 s, the same three types of plots are presented. The Ground Truth and Prediction plots show increased ice thickness compared to Panel A, while the Absolute Error plot indicates the error distribution. Panel C: At t = 0.40 s, the plots again show the Ground Truth, Prediction, and Absolute Error. The ice thickness in the Ground Truth and Prediction plots is further increased, and the Absolute Error plot shows the corresponding error distribution.

Figure 6.

Selected thickness-field snapshots predicted by the PGFNO model for the test case (V _∞ = 25 m s⁻¹, Q = 527 ml min⁻¹) under the casewise splitting scheme at three time instances: (a) t = 0.032 s; (b) t = 0.096 s; (c) t = 0.32 s. Each row presents (left to right): ground-truth thickness field, PGFNO prediction and absolute error between prediction and ground truth.

Three sets of 3D surface plots compare ground truth, predictions, and errors for thickness fields at different times.

The results demonstrate that the PGFNO model successfully generalises across distinct flow regimes/conditions (i.e. film and rivulet flow regimes) defined by varying free-stream velocities (V _∞) and flow rates (Q). Across all eight test cases, the average L ₁ norm error is 15 μm. While this represents a modest increase relative to the pairwise split model (4 μm), it remains below the experimental DIP measurement uncertainty (ε ≈ 20 μm), confirming that predictions maintain physical fidelity within the bounds of measurement precision. This performance degradation is expected given the increased complexity of generalising across entirely unseen flow regimes, yet the model’s ability to maintain sub-uncertainty-level errors underscores its robustness for interpolation across the operational parameter space of V _∞ and Q.

3.3. The robustness of the PGFNO against noise

For the trained PGFNO model to be viable in experimental applications, robustness to input noise is essential. To assess this robustness, Gaussian white noise with varying amplitudes was added to the input thickness snapshots, and the resulting output error (i.e. average L ₁ error) was quantified. Given its broader applicability, the model trained on the case-split dataset (i.e. the model used in § 3.2) was selected for this evaluation, using the same eight test cases. As illustrated in figure 7, noise amplitudes ranging from 0 % to 200 % of ε were examined, where ε = 20 μm represents the DIP measurement uncertainty. The results demonstrate that a 200 % of ε increase in input noise leads to only an approximately 5 % of ε absolute increase in the output error (i.e. from about 70 % of ε to 75 % of ε), indicating strong noise tolerance. This robust performance is particularly significant when noise magnitudes are comparable to typical measurement uncertainties, suggesting that the model maintains high fidelity under realistic experimental conditions. These findings support the model’s potential for reliable deployment in practical measurement scenarios, where measurement noise embedded in input snapshots is inevitable.

Figure 7.

Evaluation of the trained PGFNO model’s robustness to experimental white noise.

A line graph showing the relationship between noise added to input and error in output.

3.4. Film-scale distribution spectrum

The WDRWF/R flow exhibits multiscale characteristics, encompassing features of varying spatial scales (e.g. films, rivulets and droplets) under different combinations of V _∞ and Q. Consequently, evaluating the trained PGFNO model’s capability to predict these diverse spatial scales is essential. To facilitate this assessment, we introduce a film-scale distribution spectrum, analogous to the energy spectrum used for characterising velocity fields in turbulent flows.

Let $P(\boldsymbol{x})$ denote the film thickness at the position $\boldsymbol{x}_{\kern-1pt j}=(x_{\kern-1pt j},y_{\kern-1pt j})$ at a given time moment, where the subscript $j\in \{1,2,\ldots ,W\times L\}$ is an index used to denote individual grid points in the discretised domain. First, we define the fluctuating component of the film thickness:

(3.1)

\begin{align} p'(\boldsymbol{x}_{\kern-1pt j})=P(\boldsymbol{x}_{\kern-1pt j})-\left\langle P(\boldsymbol{x}_{\kern-1pt j})\right\rangle , \end{align}

where $\langle P(\boldsymbol{x}_{\kern-1pt j})\rangle$ represents the spatial average of the film thickness over the entire spatial domain:

(3.2)

\begin{align} \left\langle P\!\left(\boldsymbol{x}_{\kern-1pt j}\right)\right\rangle =\frac{1}{W\times L} \mathop{\sum }\limits_{j=1}^{W\times L}P\left(\boldsymbol{x}_{\kern-1pt j}\right). \end{align}

The two-point thickness autocorrelation function $R_{P}$ is then defined as

(3.3)

\begin{align} R_{P}(\boldsymbol{r}_{m})=\left\langle p'(\boldsymbol{x}_{\kern-1pt j})\boldsymbol{\cdot }p'(\boldsymbol{x}_{\kern-1pt j}+\boldsymbol{r}_{m})\right\rangle , \end{align}

where $\boldsymbol{r}_{m}$ is the displacement vector between two points and $m\in \{1,2,\ldots ,W\times L\}$ . The film-scale distribution spectrum $S_{P}$ is then computed by applying the two-dimensional FFT to the autocorrelation function:

(3.4)

\begin{align} S_{P}\!\left(\boldsymbol{k}\right)= \mathop{\sum }\limits_{m=1}^{W\times L}R_{P}\left(\boldsymbol{r}_{m}\right){\rm e}^{-\textrm{i}\boldsymbol{k}\boldsymbol{\cdot }{\boldsymbol{r}_{m}}}\Delta x\Delta y, \end{align}

where $\boldsymbol{k}=(k_{x},k_{y})$ represents discrete wave vectors in Fourier space and

(3.5)

\begin{align} k_{x}=\frac{2\pi n_{x}}{L},\quad k_{y}=\frac{2\pi n_{y}}{W}, \end{align}

for $n_{x}=\{0,1,\ldots ,\lfloor L/2\rfloor -1\}$ and $n_{y}=\{0,1,\ldots ,\lfloor W/2\rfloor -1\}$ , representing the wavenumber indices. Although the Fourier kernel in (3.4) is complex, S _p( k ) is real because the autocorrelation R _P ( r _m) is real and even; the imaginary sine terms cancel over the overall displacement domain, leaving no phase information.

Figure 8 presents the two-dimensional spatial wavenumber spectra of the film thickness distribution for the WDRWF/R flow field. In each panel, the x and y axes represent the wavenumbers in the streamwise (k _x) and spanwise (k _y) directions, respectively. The contour magnitude quantifies the spectral energy at each wavenumber pair, providing a measure of the prevalence of geometric features at corresponding spatial scales. The spectra shown in figure 8 correspond to the thickness fields presented in figure 6. From figures 8(a) and 8(b), the spectral contour is stretched along the y direction, i.e. the spectral energy is concentrated more prominently along the spanwise wavenumber axis (k _y) than along the streamwise axis (k _x). This anisotropy reflects the dominance of small-scale spanwise structures, i.e. the varying widths of individual rivulets visible in figure 6(b), whereas the streamwise direction is characterised by larger-scale features corresponding to rivulet lengths. From figures 8(b) and 8(c), the contour shrinks in the y direction due to the merging of the adjacent rivulets under the perturbation of aerodynamic forces on the surface waves on the water film.

Figure 8.

Film-scale distribution spectra of selected thickness-field snapshots predicted by the PGFNO model for the test case (V _∞ = 25 m s⁻¹, Q = 527 ml min⁻¹) under the casewise splitting scheme. Each row corresponds to a different time instance: (a) time = 0.032 s; (b) time = 0.096 s; (c) time = 0.32 s. Each row presents (left to right): ground-truth spectrum, predicted spectrum and absolute error spectrum.

Three sets of contour plots compare experimental, predicted, and error data for different time instances.

Regarding the prediction accuracy in film scale, qualitative comparison between predicted and ground-truth spectra in figure 8 reveals strong morphological agreement across all temporal states. To provide quantitative validation, the two-dimensional correlation coefficient between predicted and ground-truth spectra was computed at each time instant and subsequently averaged over the entire temporal domain. For the case illustrated in figure 8, the time-averaged correlation coefficient is 0.98, demonstrating that the PGFNO model accurately captures the multiscale spatial characteristics of the WDRWF/R flow field across the full range of relevant length scales.

3.5. Effect of the number of embedding layers

As described in § 2.3, the number of embedding layers in the PGFNO architecture directly governs the model’s capacity to resolve intricate features of the WDRWF/R flow field. This subsection systematically investigates the influence of embedding-layer depth on the predictive performance of PGFNO for WDRWF/R flow dynamics. Figure 9(d) presents the prediction error as a function of the number of embedding layers ( $D'$ ), evaluated across the test dataset. The results reveal a rapid decrease in prediction error with increasing $D'$ , confirming that a deeper embedding latent space enhances feature extraction and representational capacity. Due to GPU memory limitations, the maximum number of embedding layers examined in this study was constrained to 50. To elucidate the qualitative impact of embedding-layer depth, instantaneous predictions from two representative configurations ( $D'$ = 5 and $D'$ = 40) are compared against ground-truth data in figures 9(a–c) and 9(e–f). Note that $D'$ = 40 is also the default setting for the proposed PGFNO model in this study. The instantaneous thickness-field comparison is illustrated using the case at V _∞ = 15 m s⁻¹ and Q = 357 ml min⁻¹at t = 0.288 s. The absolute error maps reveal that the shallow architecture ( $D'$ = 5) exhibits discernible errors in capturing the topological complexity of the WDRWF/R flow, particularly in regions with fine-scale structures. In contrast, the deeper architecture ( $D'$ = 40) demonstrates negligible deviations, indicating substantially improved spatial fidelity. The computational cost associated with increasing the number of embedding layers is presented in figure 9(d), where the training time, obtained using a single NVIDIA A100 GPU, increases with the number of embedding layers.

Figure 9.

The effect of the number of embedding layers on the prediction accuracy of the PGFNO model. (a) The ground truth of the representative instantaneous thickness field of case V _∞ = 15 m s⁻¹, Q = 357 ml min⁻¹ at t = 0.288 s. (b) The predicted instantaneous thickness field of the model with D^′ = 5. (c) The predicted instantaneous thickness field of the model with D^′ = 40. (d) PGFNO model’s L ₁ error on the test datasets versus the number of embedding layers and the corresponding training time of the model. (e) The absolute error of the prediction of the model with D^′ = 5. (f) The absolute error of the prediction of the model with D^′ = 40. (g) Evolution of the autocorrelation coefficient. (h) The PDF of the ground truth and predicted thickness fields for the model with D^′ = 5. (i) The PDF of the ground truth and predicted thickness fields for the model with D^′ = 40.

Multiple graphs depict the effect of embedding layers on prediction accuracy in a model.

Figure 9. Long description

Panel A: A 3D surface plot shows the ground truth of the instantaneous thickness field for a specific case at a given time. The axes are labeled with X/D, Y/D, and H (mm). Panel B: A 3D surface plot shows the predicted instantaneous thickness field of the model with D = 5. The axes are labeled with X/D, Y/D, and H (mm). Panel C: A 3D surface plot shows the predicted instantaneous thickness field of the model with D = 40. The axes are labeled with X/D, Y/D, and H (mm). Panel D: A line graph shows the L1 error on the test datasets versus the number of embedding layers and the corresponding training time of the model. The x-axis is labeled Number of embedding layers, the left y-axis is labeled L1 Error, and the right y-axis is labeled Training time (Hour). Panel E: A 3D surface plot shows the absolute error of the prediction of the model with D = 5. The axes are labeled with X/D, Y/D, and Abs. Err. (mm). Panel F: A 3D surface plot shows the absolute error of the prediction of the model with D = 40. The axes are labeled with X/D, Y/D, and Abs. Err. (mm). Panel G: A line graph shows the evolution of the autocorrelation coefficient over time lag (s). The x-axis is labeled Time lag (s) and the y-axis is labeled Correlation coefficient. Panel H: A histogram shows the probability density function (PDF) of the ground truth and predicted thickness fields for the model with D = 5. The x-axis is labeled Thickness (mm) and the y-axis is labeled P.D.F. Panel I: A histogram shows the PDF of the ground truth and predicted thickness fields for the model with D = 40. The x-axis is labeled Thickness (mm) and the y-axis is labeled P.D.F.

To assess temporal predictive accuracy, the spatial autocorrelation function of the thickness field is presented in figure 9(g). While the $D'$ = 5 model shows a modest deviation from the ground truth, the $D'$ = 40 model exhibits excellent agreement throughout the correlation length scales, demonstrating superior capability in preserving spatial coherence during temporal evolution. To identify the thickness regimes most susceptible to prediction error, the probability density function (PDF) of the thickness field was computed over the entire temporal domain (figures 9 h and 9 i), comparing ground-truth distributions with predictions from both architectures. This analysis reveals that errors in the shallow model ( $D'$ = 5) are concentrated predominantly in two critical regimes: the thin-film region (i.e. H ≈ 0.1 mm) and the thick-film region (i.e. H > 1.5 mm), corresponding physically to the contact-line vicinity and surface-wave-crest regions, respectively. These regions are characterised by steep gradients and complex interfacial dynamics, requiring higher model capacity for accurate resolution. The deeper architecture ( $D'$ = 40) successfully mitigates these deficiencies, achieving close agreement with the ground-truth PDF across the full thickness range.

3.6. Effect of physical loss on model performance

As described in § 2.3, the PGFNO incorporates the principles of mass conservation and mass flow rate conservation for WDRWF/R flows directly into the loss function of the standard FNO architecture. The total loss comprises three components corresponding to the relative errors in pointwise thickness prediction (L_data ), total mass (L_phy,mass ) and mass flow rate (L_{phy,mass flow rate} ). As shown in figure 10, all unweighted loss components for test datasets decrease consistently during training and remain within approximately one order of magnitude at convergence, indicating a well-balanced formulation in which no term dominates or becomes negligible. Notably, the mass flow rate loss remains higher than the other components because it is computed from the difference in total mass between consecutive time steps and normalised by the corresponding true mass change, which is typically much smaller than the total mass. As a result, the mass flow rate loss depends on the accuracy of two neighbouring predictions and is more sensitive to small inconsistencies between them. To assess the impact of these physical constraints on predictive performance, an ablation study was conducted using test data at V _∞ = 15 m s⁻¹ and Q = 357 ml min⁻¹, comparing models trained with and without the physics-guided loss terms.

Figure 10.

Evolution of the unweighted loss components of the PGFNO model during training.

A line graph showing the evolution of loss components during training.

Figure 11(a–d) presents four diagnostic metrics that characterise different aspects of the flow field: temporal evolution of the total runback water mass, temporal evolution of the total runback water mass flow rate, streamwise spectral energy distribution of the thickness field and spanwise spectral energy distribution of the thickness field, respectively. These metrics collectively evaluate the model’s ability to capture both global integral quantities and multiscale spatial structures.

Figure 11.

The effect of the physical loss term on the PGFNO model’s performance. Taking the case of V _∞ = 15 m s⁻¹, Q = 357 ml min⁻¹ as an example: (a) the mass evolution, (b) the mass flow rate evolution, (c) the spectrum of the thickness field along the streamwise direction and (d) the spectrum of the thickness field along the spanwise direction.

Four line graphs compare the performance of PGFNO models with and without a physical loss term.

The results demonstrate that the PGFNO model with embedded physical constraints achieves markedly superior accuracy across all four metrics when benchmarked against the ground-truth thickness field. In contrast, the baseline FNO model without physical constraints exhibits substantially larger deviations. This performance disparity highlights the crucial role of physics-guided regularisation in ensuring physical consistency and enhancing model fidelity, particularly for quantities directly governed by the incorporated conservation laws. These findings validate the efficacy of the physics-guided approach in enhancing both the physical plausibility and quantitative accuracy of data-driven predictions for complex interfacial flows.

3.7. Effect of the number of Fourier blocks

As noted in § 2.3, the Fourier block represents a modular component whose depth can be adjusted to enhance PGFNO’s representational capability. Here, we examine how varying the number of stacked Fourier blocks (N _b) influences predictive fidelity for WDRWF/R flows. The relationship between Fourier block depth and prediction accuracy is shown in figure 12(a). A monotonic reduction in error is observed as N _b increases from 1 to 8, with the most pronounced improvements occurring at lower values. Hardware memory constraints limited our investigation to a maximum of eight blocks. This trend indicates that iterative spectral processing through multiple Fourier blocks progressively refines the model’s ability to capture multifrequency flow features. The computational cost associated with increasing the number of Fourier blocks is illustrated in figure 12(a), which shows that the model training time (on a single NVIDIA A100 GPU) scales approximately linearly with the number of Fourier blocks.

Figure 12.

Effect of the number of Fourier blocks on the prediction accuracy of the PGFNO model. (a) The L ₁ error of the PGFNO model on test datasets as a function of the number of Fourier blocks, and the corresponding training time of the model. (b) The PDF of the ground truth and predicted thickness fields for the model with N _b = 1, shown for a representative case with V _∞ = 15 m s⁻¹ and Q = 357 ml min⁻¹ at t = 0.288 s. (c) The PDF of the ground truth and predicted thickness fields for the model with N _b = 8. (d) Spatial correlation map of the ground-truth thickness field. (e) Spatial correlation map of the predicted thickness field for the model with N _b = 1. (f) Absolute spatial correlation error map for the predicted thickness field of the model with N _b = 1. (g) Evolution of the autocorrelation coefficient. (h) Spatial correlation map of the predicted thickness field for the model with N _b = 8. (i) Absolute spatial correlation error map for the predicted thickness field of the model with N _b = 8.

Multiple graphs depict the effect of Fourier blocks on prediction accuracy and training time, including error analysis and spatial correlation maps.

Figure 12. Long description

Panel A: A line graph shows the L1 error on the left y-axis and training time on the right y-axis as a function of the number of Fourier blocks on the x-axis. Panel B: A histogram compares the probability density function (PDF) of ground truth and predicted thickness fields for a model with 1 Fourier block. Panel C: Another histogram compares the PDF of ground truth and predicted thickness fields for a model with 8 Fourier blocks. Panel D: A contour plot displays the spatial correlation map of the ground-truth thickness field. Panel E: A contour plot shows the spatial correlation map of the predicted thickness field for the model with 1 Fourier block. Panel F: A contour plot illustrates the absolute spatial correlation error map for the predicted thickness field of the model with 1 Fourier block. Panel G: A line graph depicts the evolution of the autocorrelation coefficient over time. Panel H: A contour plot shows the spatial correlation map of the predicted thickness field for the model with 8 Fourier blocks. Panel I: A contour plot illustrates the absolute spatial correlation error map for the predicted thickness field of the model with 8 Fourier blocks.

Taking the same test case utilised in § 3.6 as an example, figures 12(b) and 12(c) reveal where residual errors manifest within the thickness distribution. A PDF analysis demonstrates that the single-block architecture struggles disproportionately at distribution extremes: the near-contact-line regime (i.e. H = 0.1–0.2 mm) and the wave-crest regime (i.e. H > 1.3 mm). Both regions involve sharp interfacial transitions and localised flow accelerations that require multiscale spectral resolution. The eight-block model substantially reduces these biases, producing PDF in close statistical agreement with experimental measurements across the full thickness spectrum.

Although the film-scale spectrum $S_{P}$ and film spatial correlation $R_{P}$ introduced in § 3.4 are mathematically related through the Fourier transformation, they provide complementary perspectives for evaluating model performance in multiphase WDRWF/R flows. The energy spectrum quantifies variance distribution across spatial scales, revealing whether the predicted film dynamics accurately captures scale-dependent energy content. However, spectral analysis is inherently phase-insensitive and thus cannot detect spatial misalignment of coherent liquid structures or distortions in rivulet patterning. In contrast, spatial correlation analysis directly quantifies spatial coherence and directional organisation, maintaining sensitivity to structural alignment, anisotropic features and characteristic correlation lengths, which are attributes that are critical for characterising streamwise rivulet development and spanwise film spreading under aerodynamic forcing. Consequently, spatial correlation evaluation complements spectral analysis by revealing deviations in multiphase structural organisation that remain obscured in the frequency domain, thereby providing essential validation of the model’s capability to reproduce WDRWF/R flow topology.

To illustrate this complementarity, figures 12(d), 12(e) and 12(h) present the film spatial correlation of the instantaneous thickness field for the test case at V _∞ = 15 m s⁻¹ and Q = 357 ml min⁻¹ at t = 0.288 s. The ground-truth spatial correlation is compared against predictions from models with N _b = 1 and N _b = 8. Notably, the spatial correlation maps directly reveal the anisotropic characteristics of the flow, with pronounced elongation in the streamwise direction reflecting the rivulets’ large-scale structure in this direction, while the correlation decay in the spanwise direction quantifies the extent of rivulet spreading laterally. To facilitate quantitative assessment of structural alignment accuracy, the absolute error maps of spatial correlation are shown in figures 12(f) and 12(i). The model with N _b = 1 exhibits a maximum error of approximately 10 %. In contrast, the model with N _b = 8 achieves substantially improved accuracy, with maximum errors below 1 %, demonstrating enhanced fidelity not only in capturing the thickness distribution of the multiphase flow field but also in preserving the directional coherence and characteristic correlation lengths that govern rivulet topology. This superior performance confirms the model’s ability to reproduce the precise spatial arrangement of liquid structures under aerodynamic forcing, validating its predictive capability for WDRWF/R flow.

Temporal consistency is evaluated through spatial autocorrelation analysis in figure 12(g). The shallow model exhibits detectable decorrelation relative to the ground truth, suggesting limitations in preserving spatial structure during time marching. The deep model maintains high correlation fidelity across all lag distances, confirming robust temporal prediction capability.

3.8. Effect of the number of Fourier modes

As outlined in § 2.3, the Fourier mode truncation in PGFNO spans three dimensions: streamwise, spanwise and temporal. While our baseline implementation employs equal mode counts across all dimensions (parametrised by K), preliminary investigations have revealed that streamwise, spanwise and temporal modes exhibit comparable influences on model performance. Consequently, this subsection focuses specifically on the sensitivity to streamwise Fourier modes (N _x), with results generalised to the spanwise and temporal dimensions.

The evaluation protocol involves training independent PGFNO models with varying N _x while holding spanwise and temporal mode counts fixed at their baseline value of 14. This isolates the contribution of streamwise spectral resolution to overall predictive accuracy. Figure 13(a) demonstrates that prediction error exhibits a sharp initial decline with increasing N _x, asymptotically approaching a performance plateau at N _x = 14. Beyond this threshold, additional streamwise modes yield negligible accuracy gains, suggesting that 14 modes are sufficient to capture the dominant streamwise wavelengths present in the WDRWF/R flow within the experimental parameter range investigated. The computational cost associated with increasing the number of Fourier modes is illustrated in figure 13(a), where the model training time, obtained using a single NVIDIA A100 GPU, exhibits an approximately linear increase with the number of Fourier modes.

Figure 13.

Impact of streamwise Fourier modes on PGFNO model accuracy. (a) The L ₁ test error versus number of streamwise Fourier modes, and the corresponding model training time. The RMS of (b) streamwise and (c) spanwise gradient thickness fields for predicted instantaneous thickness (case: V _∞ = 15 m s⁻¹, Q = 357 ml min⁻¹, t = 0.288 s) comparing N _x = 1 and N _x = 14 models.

Three line graphs depict the impact of streamwise Fourier modes on model accuracy and training time.

To rigorously assess the spatial resolution requirements for capturing interfacial complexity, we conducted a detailed gradient analysis using the representative test case from § 3.6. The interfacial thickness gradients in both streamwise (∂H/∂x) and spanwise (∂H/∂y) directions were computed throughout the entire spatiotemporal domain, providing a comprehensive measure of the model’s capability to resolve steep interfacial transitions and small-scale topological features inherent to the WDRWF/R flow. The gradient fluctuation statistics were processed as follows: for streamwise gradients ∂H/∂x, root-mean-square (RMS) values were computed, spatially averaged in the streamwise direction and temporally averaged to yield spanwise RMS profiles (figure 13 b); conversely, for spanwise gradients ∂H/∂y, RMS values were spatially averaged in the spanwise direction and temporally averaged to produce streamwise RMS profiles (figure 13 c). These gradient-based metrics exhibit high sensitivity to localised interfacial features, including rivulet edges, surface waves and irregular contact lines.

Comparative analysis reveals distinct resolution-dependent behaviour across directional gradients. In the streamwise direction (figure 13 b), the coarse-resolution configuration (N _x = 1) underestimates gradient magnitudes by approximately 20 % at peak locations, with the most pronounced discrepancies occurring in regions of high interfacial curvature where capillary and inertial forces interplay. This substantial underprediction directly reflects the model’s inability to resolve highly fluctuating interfacial wave structures. In contrast, the fine-resolution model (N _x = 14) accurately reproduces the reference RMS profile, demonstrating adequate capability in capturing both the steep interfacial fluctuations on the surface of rivulets and the fine-scale geometric complexity associated with surface wave trains propagating along the WDRWF/R flow surface.

The spanwise gradient analysis (figure 13 c) presents a fundamentally different resolution sensitivity pattern. While the N _x = 1 configuration again exhibits reduced accuracy relative to the N _x = 14 model, the magnitude of prediction error is notably smaller, i.e. approximately 2 % versus the 20 % observed for streamwise gradients. This directional asymmetry in resolution requirements can be attributed to the underlying physics governing WDRWF/R flows. The primary driving mechanism (i.e. aerodynamic shear stress exerted by the gas phase) acts predominantly in the streamwise direction, generating instabilities that manifest as streamwise-propagating surface waves. Conversely, spanwise rivulet widths are shorter than streamwise lengths and remain relatively invariant once the flow is fully developed, resulting in reduced fluctuation and geometric complexity compared with streamwise surface waves.

This anisotropic resolution behaviour has important implications for optimising the PGFNO model architecture. The results suggest that adaptive Fourier mode approaches, which prioritise streamwise refinement in regions of high gradient activity while maintaining moderate spanwise resolution, may offer an efficient pathway for accurate WDRWF/R prediction without the computational burden of uniformly high Fourier mode numbers in all three dimensions.

3.9. Long-term prediction performance

While the PGFNO architecture is trained on fixed-length temporal sequences, it can achieve extended-horizon forecasting through autoregressive deployment, where predicted states are iteratively fed as initial conditions for subsequent prediction windows. This section evaluates the model’s performance in this recurrent prediction mode using a PGFNO trained with N _out = 50, corresponding to a training time horizon from t = 0 to t = 0.4 s. As illustrated in figure 14(a), the forecasting protocol employs a sliding window approach in which the final 10 snapshots from each prediction window provide initial conditions for the next, enabling sequential forecasting across time steps 51–100 (i.e. t = 0.4–0.8 s) and 101–150 (i.e. t = 0.8–1.2 s) through two successive autoregressive iterations.

Figure 14.

Long-term prediction performance of PGFNO for the test case with V _∞ = 15 m s⁻¹ and Q = 357 ml min⁻¹. (a) Schematic of autoregressive deployment for extended-horizon forecasting, where the time step is 0.008 s given experimental data acquisition frequency of 125 Hz. The temporal evolution of (b) total water mass, (c) water mass flow rate and (d) L ₁ error.

Multiple graphs depict the performance of PGFNO for predicting water mass, flow rate, and error over time.

Figure 14. Long description

Panel A: A schematic diagram showing the autoregressive deployment for extended-horizon forecasting. The diagram includes multiple 3D surface plots representing input and output data over time. Panel B: A line graph depicting the total runback water mass in grams over time in seconds. The graph shows a rapid increase in water mass followed by a steady state. Panel C: A line graph showing the mass flow rate in grams per second over time in seconds. The graph displays fluctuating values with peaks and troughs. Panel D: A line graph illustrating the L1 error in millimeters over time in seconds. The graph indicates a low error initially, followed by a significant increase.

To assess autoregressive performance, we examine the representative test case of V _∞ = 15 m s⁻¹ and Q = 357 ml min⁻¹, consistent with previous analyses. Figure 14(b) presents the temporal evolution of total water mass under autoregressive operation. The predicted trajectory maintains close fidelity to ground truth throughout the extended forecast horizon, demonstrating robust conservation of this global integral quantity despite cumulative error propagation. Similarly, the mass flow rate evolution in figure 14(c) exhibits excellent agreement with reference data for extended forecast time steps 51–100 (i.e. t = 0.4–0.8 s). While predictions for steps 101–150 (i.e. t = 0.8–1.2 s) show modest quantitative deviations, the model critically preserves the characteristic oscillation frequency, indicating that dominant temporal dynamics remains accurately captured even as local errors accumulate.

The corresponding L ₁ error evolution quantified in figure 14(d) reveals the expected compounding behaviour inherent to autoregressive forecasting. Notably, however, the maximum accumulated error remains approximately 0.02 mm, which is comparable to the experimental measurement uncertainty (ε ≈ 0.02 mm). This level of accuracy is particularly significant given that it represents prediction fidelity after multiple autoregressive cycles, where small inaccuracies can potentially amplify through iterative propagation.

To contextualise these results, we compare the PGFNO’s predictive horizon against a previously reported state-of-the-art ConvLSTM model (Wang et al. Reference Wang, Hu, He and Hu2023b ) for WDRWF/R forecasting. The ConvLSTM architecture maintained sub-uncertainty prediction accuracy for only 10 time steps at a 30 Hz sampling, corresponding to approximately 0.33 s of physical time. In contrast, the present PGFNO model achieves comparable accuracy over 150 time steps at a sampling rate of 125 Hz, representing 1.2 s of physical time, which is nearly a fourfold improvement in temporal extrapolation capability. This substantial advancement demonstrates the superior capacity of PGFNO to capture and propagate the complex spatiotemporal dynamics governing WDRWF/R flows.

Collectively, these results establish PGFNO’s capability for extended temporal extrapolation of WDRWF/R dynamics, maintaining predictive capability across time scales substantially exceeding the training sequence length. The simultaneous preservation of integral conservation properties (i.e. mass and mass flow rate) and characteristic temporal frequencies under autoregressive operation provides compelling evidence that the model has internalised the fundamental governing physics of multiphase flow, rather than merely interpolating temporal patterns. This physics-guided learning is essential for reliable long-horizon forecasting in complex fluid dynamics applications. In the context of aircraft icing, these capabilities are particularly significant, as accurate prediction of WDRWF/R dynamics is critical to ice accretion modelling. In particular, the ability of PGFNO to capture and propagate the coupled spatiotemporal evolution of surface water enables improved characterisation of water redistribution, which governs ice growth, shape formation and accretion rate under glaze icing conditions. The demonstrated long-horizon predictive capability further provides a foundation for integrating data-driven flow evolution models into existing icing simulation frameworks, thereby enhancing their ability to resolve transient surface water transport processes that are not adequately captured by conventional modelling approaches.

4. Conclusion

This study presents a PGFNO model for predicting the transient evolution of WDRWF/R flows. By integrating pointwise data loss with integral physical constraints, i.e. mass conservation and mass flow rate conservation, the PGFNO framework achieves enhanced accuracy and generalisability. The model was trained on a high-fidelity experimental dataset acquired through advanced DIP techniques, which captures the intricate physics of WDRWF/R flows without simplifying assumptions. Notably, the trained model demonstrates robust generalisation across varying initial conditions, including different free-stream velocities (V _∞) and flow rates (Q). To rigorously assess the model’s capability in resolving multiscale flow features, a novel film-scale distribution spectrum metric was introduced. The PGFNO model exhibited exceptional resilience to experimental noise and maintained reliable long-term predictive performance. Comprehensive ablation studies elucidated the influence of physical loss constraints and architectural parameters (including the number of embedding layers, Fourier blocks and Fourier modes) on model performance, thereby identifying critical components essential for capturing complex interfacial flow physics. These findings underscore the significant promise of the PGFNO framework for advancing the understanding and prediction of challenging multiphase WDRWF/R flows.

Future work will extend the model to incorporate additional interfacial phenomena, including contact-line dynamics, surface roughness effects and surface energy variations, thereby enabling more refined predictions of coating-induced modifications to film thickness and runback characteristics. The model will be validated under more demanding operational conditions, including high-Weber-number regimes, high liquid water content and geometrically complex curved surfaces, through integration with high-fidelity physics-based simulations. Furthermore, this computational framework will be generalised to address other critical runback multiphase flow challenges, such as evaporative runback on heated substrates and ice accretion processes on supercooled surfaces, thereby broadening its applicability to diverse engineering scenarios of practical importance.

Acknowledgements

The authors would like to thank the 2024 Center for Turbulence Research (CTR) Summer Program at Stanford University for their support.

Funding

This research was partially supported by the Iowa Energy Center under IEC Competitive Grant no. 312350 and the National Science Foundation (NSF) under Award nos. CBET-1935363, TIP-2140489, and CBET-2313310.

Declaration of interests.

The authors have no conflicts of interest to disclose.

Data availability statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

Anderson, D. & Tsao, J. 2003 Overview of icing physics relevant to scaling. In FAA In-flight Icing / Ground De-icing International Conference & Exhibition, Chicago, Illinois, United States, June 16, 2003. SAE International.10.4271/2003-01-2130CrossRef Google Scholar

Bragg, M.B., Gregorek, G.M. & Lee, J.D. 1986 Airfoil aerodynamics in icing conditions. J. Aircraft 23 (1), 76–81.10.2514/3.45269CrossRef Google Scholar

Cao, Y., Tan, W. & Wu, Z. 2018 Aircraft icing: an ongoing threat to aviation safety. Aerosp. Sci. Technol. 75, 353–385.10.1016/j.ast.2017.12.028CrossRef Google Scholar

Cao, Y. & Yuan, K. 2007 Aircraft flight characteristics in conditions of windshear and icing. Aeronaut. J. 111 (1115), 41–49.10.1017/S0001924000001743CrossRef Google Scholar

Cao, Y., Yuan, K. & Li, G. 2011 Effects of ice geometry on airfoil performance using neural networks prediction. Aircr. Engng Aerosp. Tech. 83 (5), 266–274.10.1108/00022661111159870CrossRef Google Scholar

Cebeci, T. & Kafyeke, F. 2003 AIRCRAFT ICING. Annu. Rev. Fluid Mech. 35, 11–21.10.1146/annurev.fluid.35.101101.161217CrossRef Google Scholar

Chen, N., Yi, X., Wang, Q. & Liu, Y. 2023 Numerical study on wind-driven thin water film runback on an airfoil. AIAA J. 61 (6), 2517–2525.10.2514/1.J062618CrossRef Google Scholar

Cox, R.G. 1986 The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow. J. Fluid Mech. 168, 169–194.10.1017/S0022112086000332CrossRef Google Scholar

Feo, A. 2000 Icing Scaling with Surface Film Thickness Similarity for High LWC Conditions. Instituto Nac. Tec. Aeroesp. Rep., Madrid INT, p. AE/PRO/4420/184/INTA/00. Available at: https://Scholar.Google.Com/Scholar.Google Scholar

Gao, L., Liu, Y. & Hu, H. 2020 An experimental investigation on the dynamic glaze ice accretion process over a wind turbine airfoil surface. Intl J. Heat Mass Transfer 149, 119120.10.1016/j.ijheatmasstransfer.2019.119120CrossRef Google Scholar

Gendrich, C.P. & Koochesfahani, M.M. 1996 A spatial correlation technique for estimating velocity fields using molecular tagging velocimetry (MTV). Exp. Fluids 22 (1), 67–77.10.1007/BF01893307CrossRef Google Scholar

de Gennes, P.G. 1985 Wetting: statics and dynamics. Rev. Mod. Phys. 57 (3), 827–863.10.1103/RevModPhys.57.827CrossRef Google Scholar

Gent, R.W., Dart, N.P. & Cansdale, J.T. 2000 Aircraft icing. Phil. Trans. R. Soc. A: Math. Phys. Engng Sci. 358 (1776), 2873–2911.10.1098/rsta.2000.0689CrossRef Google Scholar

Hansman, R.J. & Turnock, S.R. 1989 Investigation of surface water behavior during glaze ice accretion. J. Aircraft 26 (2), 140–147.10.2514/3.45735CrossRef Google Scholar

Heinrich, A., et al. 1991 Aircraft Icing Handbook, Volume 1 of 3, vol. 1, p. ADA238040. Gates Learjet Corporation.Google Scholar

Huh, C. & Scriven, L.E. 1971 Hydrodynamic model of steady movement of a solid/liquid/fluid contact line. J. Colloid Interface Sci. 35 (1), 85–101.10.1016/0021-9797(71)90188-3CrossRef Google Scholar

Joanny, J.F. & de Gennes, P.G. 1984 A model for contact angle hysteresis. J. Chem. Phys. 81 (1), 552–562.10.1063/1.447337CrossRef Google Scholar

Kovachki, N., Li, Z., Liu, B., Anandkumar, A., et al. 2023 Neural operator: learning maps between function spaces with applications to PDEs. J. Mach. Learn. Res. 24 (89), 1–97.Google Scholar

Krishnapriyan, A., Gholami, A., Zhe, S., Kirby, R. & Mahoney, M.W. 2021 Characterizing possible failure modes in physics-informed neural networks. In Advances in Neural Information Processing Systems (ed. M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang & J.W. Vaughan), vol. 34, pp. 26548–26560. Curran Associates, Inc.Google Scholar

Li, Z., Kovachki, N., Azizzadenesheli, K., Anandkumar, A., et al. 2021 Fourier neural operator for parametric partial differential equations. Available at: https://arxiv.org/abs/2010.08895.Google Scholar

Liu, Y., Zhang, K., Tian, W. & Hu, H. 2020 An experimental study to characterize the effects of initial ice roughness on the wind-driven water runback over an airfoil surface. Intl J. Multiphase Flow 126, 103254.10.1016/j.ijmultiphaseflow.2020.103254CrossRef Google Scholar

Lu, L., Jin, P., Pang, G., Zhang, Z. & Karniadakis, G.E. 2021 Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nat. Mach. Intell. 3 (3), 218–229.10.1038/s42256-021-00302-5CrossRef Google Scholar

Lynch, F.T. & Khodadoust, A. 2001 Effects of ice accretions on aircraft aerodynamics. Prog. Aerosp. Sci. 37 (8), 669–767.10.1016/S0376-0421(01)00018-5CrossRef Google Scholar

Ma, K., Lin, G., Jin, H., Bu, X., Shen, X. & Wen, D. 2025 3D numerical investigation of surface wettability induced runback water flow behavior. Renewab. Energy 246, 122937.10.1016/j.renene.2025.122937CrossRef Google Scholar

McCann, D.W. 2005 NNICE – a neural network aircraft icing algorithm. Environ. Modell. Softw. 20 (10), 1335–1342.10.1016/j.envsoft.2004.09.027CrossRef Google Scholar

Mosher, F., Schaum, D., Herbster, C. & Guinn, T. 2010 Analysis of causes of icing conditions which contributed to the crash of Continental flight 3407. In 14th Conference on Aviation, Range, and Aerospace Meteorology. American Meteorological Society (AMS).Google Scholar

Myers, T.G. 2001 Extension to the Messinger model for aircraft icing. AIAA J. 39 (2), 211–218.10.2514/2.1312CrossRef Google Scholar

Myers, T.G., Liang, H.X. & Wetton, B. 2004 The stability and flow of a rivulet driven by interfacial shear and gravity. Intl J. Nonlinear Mech. 39 (8), 1239–1249.10.1016/j.ijnonlinmec.2003.08.001CrossRef Google Scholar

Myers, T.G. & Thompson, C.P. 1998 Modeling the flow of water on aircraft in icing conditions. AIAA J. 36 (6), 1010–1013.10.2514/2.472CrossRef Google Scholar

Ozcer, I.A., Baruzzi, G.S., Reid, T., Habashi, W.G., Fossati, M. & Croce, G. 2011 FENSAP-ICE: numerical prediction of ice roughness evolution, and its effects on ice shapes. In SAE Technical Paper. SAE International.10.4271/2011-38-0024CrossRef Google Scholar

Quéré, D. 2008 Wetting and roughness. Ann. Rev. Mater. Res. 38, 71–99.10.1146/annurev.matsci.38.060407.132434CrossRef Google Scholar

Raissi, M., Perdikaris, P. & Karniadakis, G.E. 2019 Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J. Comput. Phys. 378, 686–707.10.1016/j.jcp.2018.10.045CrossRef Google Scholar

Rothmayer, A. & Hu, H. 2015 A low-order model of the interaction between air-driven water films and sub-grid surface roughness. In SAE 2015 International Conference on Icing of Aircraft, Engines, and Structures, Prague, Czech Republic, June 22, 2015. SAE International.10.4271/2015-01-2078CrossRef Google Scholar

Rothmayer, A.P., Matheis, B.D. & Timoshin, S.N. 2002 Thin liquid films flowing over external aerodynamic surfaces. J. Engng Math. 42 (3), 341–357.10.1023/A:1016152918706CrossRef Google Scholar

Ryerson, C., Peck, L. & Martel, C. 2003 Army aviation operations in icing conditions. In FAA In-flight Icing / Ground De-icing International Conference & Exhibition, Chicago, Illinois, United States, June 16, 2003. SAE International.10.4271/2003-01-2094CrossRef Google Scholar

Servini, P., Smith, F.T. & Rothmayer, A.P. 2018 The impact of dynamic roughness elements on marginally separated boundary layers. J. Fluid Mech. 855, 351–370.10.1017/jfm.2018.589CrossRef Google Scholar

Villedieu, P., Trontin, P. & Chauvin, R. 2014 Glaciated and mixed phase ice accretion modeling using ONERA 2D icing suite. In Proceedings of the 6th AIAA Atmospheric and Space Environments Conference, Atlanta, GA, USA, June 2014, pp. 16–20. American Institute of Aeronautics and Astronautics (AIAA).Google Scholar

Villedieu, P., Trontin, P., Guffond, D. & Bobo, D. 2012 SLD Lagrangian modeling and capability assessment in the frame of ONERA 3D icing suite. In Proceedings of the 4th AIAA Atmospheric and Space Environments Conference, New Orleans, LA, USA, 25–28 January 2012, p. 3132. American Institute of Aeronautics and Astronautics (AIAA).Google Scholar

Waldman, R.M. & Hu, H. 2015 High-speed imaging to quantify transient ice accretion process over an airfoil. J. Aircraft 53 (2), 369–377.10.2514/1.C033367CrossRef Google Scholar

Wang, J., Dhulipalla, A. & Hu, H. 2026 a An experimental study on the aerodynamic performance degradation of a high-lift UAV airfoil induced by ice accretion. In AIAA SCITECH 2026 Forum. American Institute of Aeronautics and Astronautics.Google Scholar

Wang, J., Dhulipalla, A., Zhao, Y., Nunes, A. & Hu, H. 2025 An experimental investigation on dynamic water runback process over the blade surface of a rotating drone propeller. In AIAA AVIATION FORUM and ASCEND 2025. American Institute of Aeronautics and Astronautics.Google Scholar

Wang, J., He, P. & Hu, H. 2026 b A numerical study on wind-driven runback characteristics of a thin water film flow over a solid surface. Acta Mechanica Sin. 42 (1), 725073.10.1007/s10409-025-25073-xCrossRef Google Scholar

Wang, J., Hu, H., He, P. & Hu, H. 2023 a A data-driven approach to study nonlinear dynamics of wind-driven water runback flows pertinent to aircraft icing phenomena. In AIAA AVIATION 2023 Forum. American Institute of Aeronautics and Astronautics.Google Scholar

Wang, J., Hu, H., He, P. & Hu, H. 2023 b A machine learning study to predict wind-driven water runback characteristics. Phys. Fluids 35 (10), 102104.10.1063/5.0167545CrossRef Google Scholar

Wang, G. & Rothmayer, A.P. 2009 Thin water films driven by air shear stress through roughness. Comput. Fluids 38 (2), 235–246.10.1016/j.compfluid.2008.02.009CrossRef Google Scholar

Wang, S., Teng, Y. & Perdikaris, P. 2021 Understanding and mitigating gradient flow pathologies in physics-informed neural networks. SIAM J. Sci. Comput. 43 (5), A3055–A3081.10.1137/20M1318043CrossRef Google Scholar

Wright, W.B. 2008 User’s manual for LEWICE version 3.2. CR-2008-214255. National Aeronautics and Space Administration (NASA).Google Scholar

Wright, W.B., Porter, C.E., Galloway, E.T. & Rigby, D.L. 2022 GlennICE 2.1 capabilities and results. In AIAA AVIATION 2022 Forum, AIAA Paper 2022-3309. American Institute of Aeronautics and Astronautics.Google Scholar

Yousif, M.Z., Yu, L. & Lim, H. 2022 Physics-guided deep learning for generating turbulent inflow conditions. J. Fluid Mech. 936, A21.10.1017/jfm.2022.61CrossRef Google Scholar

Yu, D., Han, Z., Zhang, B., Zhang, M., Liu, H. & Chen, Y. 2022 A multi-autoencoder fusion network for fast image prediction of aircraft ice accretion. Phys. Fluids 34 (7), 76107.10.1063/5.0091068CrossRef Google Scholar

Zhang, Z. 2021 Experimental Investigations on the Interactions Between Water/Deicing Fluids and Icephobic Coatings Pertinent to Aircraft Anti-/de-icing. Iowa State University.10.31274/td-20240329-319CrossRef Google Scholar

Zhang, K. & Hu, H. 2021 An experimental study on the transient runback characteristics of wind-driven film/rivulet flows. Phys. Fluids 33 (11), 112104.10.1063/5.0067672CrossRef Google Scholar

Zhang, K., Tian, W. & Hu, H. 2015 An experimental investigation on the surface water transport process over an airfoil by using a digital image projection technique. Exp. Fluids 56 (9), 173.10.1007/s00348-015-2046-zCrossRef Google Scholar

Zhao, Y., Wang, J. & Hu, H. 2025 An experimental study on transient wind-driven water droplet runback and dynamic ice accretion process over a frozen cold surface. In AIAA AVIATION FORUM and ASCEND 2025, AIAA paper 2025-3438. American Institute of Aeronautics and Astronautics.Google Scholar

Zhong, W. & Meidani, H. 2025 Physics-informed geometry-aware neural operator. Comput. Meth. Appl. Mech. Engng 434, 117540.10.1016/j.cma.2024.117540CrossRef Google Scholar

Figure 1.Figure 1 long description.(a) The principle of the DIP technique. (b) The experimental set-up for measuring WDRWF/R flow by the DIP technique.

Figure 2.(a) A schematic diagram of quantifying the thickness field of the WDRWF/R flow by the DIP technique. (b) The displacement-to-height conversion map of the DIP system calibrated for the present study. (c) Distribution of the experimental test conditions (i.e. the incoming wind speed, V∞, and water flow rate, Q) used in this study.

Figure 3.Figure 3 long description.The PGFNO architecture for predicting WDRWF/R flow. Here Pi$\boldsymbol{P}_i$ represents a snapshot of the water film thickness field at time step i. The three-dimensional block representations illustrate the last three dimensions of each intermediate tensor. The dimensions are evaluated using the default values employed in this study.

Figure 4.Selected thickness-field snapshots predicted by the PGFNO model for the test case (V∞ = 13 m s−1, Q = 327 ml min−1) under the pairwise splitting scheme at three time instances: (a) time = 0.032 s; (b) time = 0.24 s; (c) time = 0.40 s. Each row presents (left to right): ground-truth thickness field, PGFNO prediction and absolute error between prediction and ground truth.

Figure 5.Figure 5 long description.Selected thickness-field snapshots predicted by the PGFNO model for the test case (V∞ = 5 m s−1, Q = 196 ml min−1) under the casewise splitting scheme at three time instances: (a) t = 0.064 s; (b) t = 0.20 s; (c) t = 0.40 s. Each row presents (left to right): ground-truth thickness field, PGFNO prediction and absolute error between prediction and ground truth.

Figure 6.Selected thickness-field snapshots predicted by the PGFNO model for the test case (V∞ = 25 m s−1, Q = 527 ml min−1) under the casewise splitting scheme at three time instances: (a) t = 0.032 s; (b) t = 0.096 s; (c) t = 0.32 s. Each row presents (left to right): ground-truth thickness field, PGFNO prediction and absolute error between prediction and ground truth.

Figure 7.Evaluation of the trained PGFNO model’s robustness to experimental white noise.

Figure 8.Film-scale distribution spectra of selected thickness-field snapshots predicted by the PGFNO model for the test case (V∞ = 25 m s−1, Q = 527 ml min−1) under the casewise splitting scheme. Each row corresponds to a different time instance: (a) time = 0.032 s; (b) time = 0.096 s; (c) time = 0.32 s. Each row presents (left to right): ground-truth spectrum, predicted spectrum and absolute error spectrum.

Figure 9.Figure 9 long description.The effect of the number of embedding layers on the prediction accuracy of the PGFNO model. (a) The ground truth of the representative instantaneous thickness field of case V∞ = 15 m s−1, Q = 357 ml min−1 at t = 0.288 s. (b) The predicted instantaneous thickness field of the model with D′ = 5. (c) The predicted instantaneous thickness field of the model with D′ = 40. (d) PGFNO model’s L1 error on the test datasets versus the number of embedding layers and the corresponding training time of the model. (e) The absolute error of the prediction of the model with D′ = 5. (f) The absolute error of the prediction of the model with D′ = 40. (g) Evolution of the autocorrelation coefficient. (h) The PDF of the ground truth and predicted thickness fields for the model with D′ = 5. (i) The PDF of the ground truth and predicted thickness fields for the model with D′ = 40.

Figure 10.Evolution of the unweighted loss components of the PGFNO model during training.

Figure 11.The effect of the physical loss term on the PGFNO model’s performance. Taking the case of V∞ = 15 m s−1, Q = 357 ml min−1 as an example: (a) the mass evolution, (b) the mass flow rate evolution, (c) the spectrum of the thickness field along the streamwise direction and (d) the spectrum of the thickness field along the spanwise direction.

Figure 12.Figure 12 long description.Effect of the number of Fourier blocks on the prediction accuracy of the PGFNO model. (a) The L1 error of the PGFNO model on test datasets as a function of the number of Fourier blocks, and the corresponding training time of the model. (b) The PDF of the ground truth and predicted thickness fields for the model with Nb = 1, shown for a representative case with V∞ = 15 m s−1 and Q = 357 ml min−1 at t = 0.288 s. (c) The PDF of the ground truth and predicted thickness fields for the model with Nb = 8. (d) Spatial correlation map of the ground-truth thickness field. (e) Spatial correlation map of the predicted thickness field for the model with Nb = 1. (f) Absolute spatial correlation error map for the predicted thickness field of the model with Nb = 1. (g) Evolution of the autocorrelation coefficient. (h) Spatial correlation map of the predicted thickness field for the model with Nb = 8. (i) Absolute spatial correlation error map for the predicted thickness field of the model with Nb = 8.

Figure 13.Impact of streamwise Fourier modes on PGFNO model accuracy. (a) The L1 test error versus number of streamwise Fourier modes, and the corresponding model training time. The RMS of (b) streamwise and (c) spanwise gradient thickness fields for predicted instantaneous thickness (case: V∞ = 15 m s−1, Q = 357 ml min−1, t = 0.288 s) comparing Nx = 1 and Nx = 14 models.

Figure 14.Figure 14 long description.Long-term prediction performance of PGFNO for the test case with V∞ = 15 m s−1 and Q = 357 ml min−1. (a) Schematic of autoregressive deployment for extended-horizon forecasting, where the time step is 0.008 s given experimental data acquisition frequency of 125 Hz. The temporal evolution of (b) total water mass, (c) water mass flow rate and (d) L1 error.

Article contents

A physics-guided Fourier neural operator to predict wind-driven water film/rivulet flows pertinent to aircraft icing phenomena

Abstract

JFM classification

Information

1. Introduction

2. Experimental set-up and machine learning framework

2.1. Fundamentals of the DIP technique

2.2. Experimental set-up for WDRWF/R flow

2.3. Machine learning framework for WDRWF/R flow

2.4. Data splitting method

3. Results and discussion

3.1. Evaluation of PGFNO generalisability across different initial states

3.2. Evaluation of PGFNO generalisability across different test conditions

3.3. The robustness of the PGFNO against noise

3.4. Film-scale distribution spectrum

3.5. Effect of the number of embedding layers

3.6. Effect of physical loss on model performance

3.7. Effect of the number of Fourier blocks

3.8. Effect of the number of Fourier modes

3.9. Long-term prediction performance

4. Conclusion

Acknowledgements

Funding

Declaration of interests.

Data availability statement

References

Save article to Kindle

Save article to Dropbox

Save article to Google Drive

Reply to: Submit a response

Your details

You have entered the maximum number of contributors

Conflicting interests