Fusion of multi-source precipitation records via coordinate-based generative models

Domain: earth_science
Published: 2025-12-29
Editorial outcome: major_revision

Abstract
Precipitation remains one of the most challenging climate variables to observe and predict. Existing datasets face intricate trade-offs: gauges are relatively trustworthy but sparse, satellites provide near-global coverage with retrieval uncertainties, and numerical models offer physical consistency but are biased. Here we introduce PRIMER (Precipitation Records Infinite MERging), a framework that fuses these complementary sources. PRIMER employs a coordinate-based diffusion model that learns from arbitrary spatial locations and associated intensity values, enabling seamless integration of gridded data and irregular gauge observations. Through two-stage training—first learning large-scale patterns, then refining with gauge measurements—PRIMER captures both large-scale structure and local precision. Once trained, it can correct biases in existing datasets—yielding significant error reductions at most gauge sites—and downscale reanalysis. In addition, by combining background estimates with extra gauges, it produces analysis fields that further reduce errors. All tasks are achieved through posterior sampling utilizing the prior obtained by fusing multi-source records. Crucially, it generalizes without retraining, correcting biases in operational forecasts and downscaling future scenario precipitation fields. This demonstrates how PRIMER can transform imperfect data into a source of strength.

Introduction
Precipitation–when, where, and how much water falls from the sky to the earth's surface–governs freshwater availability, agricultural productivity, flood hazards, and ecosystem health across the globe1. Despite its significance, precipitation remains one of the most challenging climate variables to observe and predict accurately. This challenge stems from precipitation’s fundamental nature: unlike most climate variables that vary smoothly across space and time, precipitation manifests as discrete, intermittent pulses with striking discontinuities2,3. These complex spatiotemporal organization depends crucially on small-scale cloud microphysics processes4 that remain poorly understood or simulated. Besides, these processes are highly sensitive to environmental conditions: small perturbations in temperature, humidity, or aerosol concentrations can determine whether clouds produce no rain, light drizzle, or torrential downpours5,6. Furthermore, the triggering and organization of convection–the primary mechanism for intense precipitation–depends on complex interactions between boundary layer turbulence7, atmospheric stability8, and mesoscale circulations9,10 that remain computationally prohibitive to simulate explicitly. These complexities create fundamental observational and predictive challenges.

Currently, we rely upon three sources to derive precipitation information, namely in situ gauge observations, remote sensing, and numerical simulation that often assimilate in situ and remote-sensed data11. Each of these three sources comes with inherent limitations regarding their accuracy, coverage, and resolution. Ground-based observations from rain gauges provide the most direct and accurate measurements at point locations. However, gauge networks exhibit severe spatial limitations: even 2.5∘×2.5∘ grid cells contain less than two gauges on average12, let alone the oceanic and remote regions. Satellite remote sensing offers near-global coverage, but measures precipitation indirectly. Passive microwave sensors on polar-orbiting satellites detect emission and scattering signatures from hydrometeors, providing relatively direct estimates but with limited temporal sampling13. Infrared sensors on geostationary satellites offer frequent observations (every 10-30 minutes) but only measure cloud-top temperatures, requiring empirical relationships to infer surface precipitation–a particularly poor assumption for shallow, warm clouds that produce significant precipitation in tropical and maritime regions14. Numerical weather prediction and reanalysis products provide physically consistent, complete spatiotemporal coverage by assimilating available observations into dynamical models15. However, precipitation in these systems emerges as the end result of a complex chain of parameterized processes—radiation, convection, cloud microphysics, and boundary layer turbulence—each contributing its own errors16, with their errors compounding multiplicatively. The consequence of these observational and simulational limitations is profound: current precipitation datasets often disagree by as much as the signal itself11,16. In tropical regions, the spread among different products can exceed 300 mm/hr of the mean precipitation11, fundamentally limiting our ability to close the global water budget, validate climate models, or provide reliable information for water resources management.

A promising solution to these challenges lies in data fusion–leveraging the complementary strengths of multiple data sources to produce precipitation estimates that surpass any individual source in accuracy, resolution, and coverage17–27. Among data-fusion approaches, Bayesian methods offer a coherent and probabilistically grounded solution. The key insight is elegant: by deriving an informative prior from all available sources, we can encode existing knowledge in a statistically coherent form. Once established, this prior can be updated via Bayes’ theorem with any new observation–accounting for each source’s unique error characteristics and observational modalities through tailored likelihood functions.28–31. The framework naturally weights observations by their reliability and propagates uncertainties to yield full posterior distributions32, essential for risk assessment.

Recent advances in deep generative models, particularly probabilistic diffusion models33,34, offer a transformative opportunity for implementing the above Bayesian framework. Diffusion models approximate target distributions by learning to reverse a gradual noising process. The forward process progressively perturbs data with Gaussian noise, while the reverse process, parameterized by a neural network, learns to invert this corruption to recover samples. This iterative denoising procedure enables diffusion models to generate high-quality and diverse samples, spanning domains from natural images35 to protein structures36, while also serving as priors for Bayesian inference, making them particularly well-suited for capturing the intricate patterns of precipitation. Once trained, they function as “plug-and-play” priors37–41: the same learned distribution can be applied to diverse inference tasks–bias correction, downscaling, or gap-filling–by simply changing the likelihood function without retraining. Despite the promises, implementing this framework for precipitation faces three fundamental challenges. First, precipitation’s extreme spatiotemporal variability—from localized convective cells to continental-scale fronts—makes it extraordinarily difficult to be captured in a single prior distribution. Second, constructing an informative prior becomes paradoxical when no individual data source is trustworthy or comprehensive. Each source captures different aspects of precipitation across mismatched scales, creating a circular dependency where we need accurate data to build a prior, yet need a prior to evaluate data accuracy. Third, even with a reasonable prior, posterior sampling remains challenging because, from a machine learning perspective, a precipitation field is high-dimensional and the associated observation likelihood is complex. These barriers define the frontier for deploying generative AI in Earth-system science, demanding innovations that transcend conventional approaches.

To address these challenges, we introduce PRIMER (Precipitation Record Infinite MERging), a general framework that reconceptualizes how diffusion models can learn from imperfect, heterogeneous precipitation records. Our key insight is that probabilistic diffusion models need not be trained on perfect samples—instead, they can be viewed as spectral regression models that progressively learn from low-frequency structures to high-frequency details as we gradually corrupt the target distribution using Gaussian noise42. This property enables us to construct an informative prior by learning conditional distributions of precipitation patterns for each data source, where the conditioning explicitly captures each dataset’s characteristic biases. As emphasized by ref. 43, integrating data with varying degrees of sparsity—from sparse grids to dense fields—poses a major machine learning challenge. We acknowledge this issue and propose an approach to better handle such heterogeneity (see SI Section 2.8 for a comparison with ref. 43). Conventionally, diffusion models work on samples residing on fixed-resolution grids44, forcing us to interpolate heterogeneous observations to common resolutions. This interpolation is particularly destructive for precipitation: it smooths sharp gradients at convective boundaries, introduces artificial correlations between sparse gauge points, and—most critically—destroys the very precision that makes gauges valuable. For sparse gauge networks covering less than 1% of the domain, interpolation essentially fabricates information that doesn’t exist. We therefore require an architecture that can learn priors directly from each source’s native sampling structure. This necessity drives our adoption of coordinate-based diffusion models, which represent precipitation as spatial fields x: R2 -> R rather than tensors. In this formulation, both dense grids and sparse gauge observations are simply different sampling patterns of the same underlying field. PRIMER directly learns from arbitrarily and sparsely distributed points—each defined by its location and precipitation intensity—without relying on spatial interpolation (see Fig. 1b)—gauge observations influence the function locally while gridded data constrain large-scale structure. Our two-stage training strategy is thus a natural choice: we first learn the baseline priors P_ERA5(x) and P_IMERG(x), which represent the climatological distributions of precipitation fields x derived from climate reanalysis, i.e., fifth-generation ECMWF atmospheric reanalysis (ERA5), and satellite-based retrieval dataset, i.e., Integrated Multi-satellitE Retrievals for GPM (IMERG). We then fine-tune the model using gauge observational information at sparse grid locations (we refer hereafter to these densely observed grid cells as “gauge observations”; see Method 4.6 for data sources and detailed descriptions), so as to incorporate local accuracy, yielding the updated prior P⋆(x) (Fig. 1b; star indicates that it is supposed to be a better prior). The coordinate-based representation ensures that gauge information enhances rather than corrupts the prior, as each source contributes at its natural scale. Once trained, PRIMER supports diverse applications through principled posterior sampling: given observations O—whether from biased satellites, sparse gauges, or coarse forecasts—we can sample from posterior P(x|O) to produce improved ensemble estimates (Fig. 1a). Empirical evaluations demonstrate the effectiveness of our approach: it achieves statistically significant error reductions for grids which are densely observed with gauges, supplements high-frequency details through downscaling, and further reduces errors by merging gauge observations with the background, in a way similar to optimal interpolation, underscoring its potential for operational use. It also generalizes to unseen operational forecasts without retraining and extends to downscaling future scenario precipitation fields in CMIP6. By transforming the challenge of heterogeneous, imperfect data from a limitation into a strength, PRIMER establishes a paradigm for precipitation data fusion that extends naturally to other Earth-system variables plagued by observational trade-offs.

Results
Reproducing climatological distributions
The gist of the PRIMER is to learn a trustworthy prior of precipitation fields, thereafter applying it for a broad range of relevant probabilistic inference tasks. Before verifying the probabilistic inference results, we should ensure the accuracy of the learned prior distribution. As directly evaluating such high-dimensional priors is intractable, we instead assess their statistical properties as proxies45–47. We compare unconditionally generated samples from P_IMERG(x), P_ERA5(x), and the updated prior P⋆(x) against their respective referenced datasets. In particular, we focus on the climatological mean and standard deviation of precipitation (Fig. 2). At the grid-point level, the agreement is clear. For mean precipitation (Fig. 2a–f), both P_IMERG(x) and P_ERA5(x) exhibit strong spatial correspondence with IMERG and ERA5, achieving Pearson correlation coefficients (PCCs) of 0.85 and 0.97, respectively. The standard deviation fields (Fig. 2g–l) are likewise well reproduced (PCC=0.75 and 0.86), highlighting PRIMER’s capacity to represent not just the average precipitation spatial structure but also its variance. Notably, we also introduce P⋆(x), constructed by fine-tuning PRIMER using sparse but reliable gauge observational information. Despite the limited spatial coverage of gauge observations, this calibration yields a climatologically consistent prior that preserves spatial structures learned from the gridded products while injecting localized realism. This “climatological jailbreak” illustrates how PRIMER can adapt to sparse gauge records without compromising coherence across scales. To further evaluate spatial structure, we perform a radially averaged power spectral density (RAPSD) analysis (Fig. 2m), which confirms that the learned priors accurately recover the multiscale spectral characteristics of the referenced datasets, especially across mesoscale wavelengths, which are crucial for convective processes (see also Supplementary Information (SI) Fig. 9). Additional statistical evaluations—including precipitation frequency, extremes, skewness, and empirical orthogonal function (EOF) modes—are provided in the SI Fig. 10.

Case study on high-impact events
The previous section evaluated PRIMER’s ability to match climatology. After Stage 2 fine-tuning, the updated prior P⋆(x) is expected to align more closely with gauge observations; however, its actual skill remains to be validated through posterior sampling experiments. To this end, we perform posterior sampling using different priors while conditioning on the same observations O. By comparing the posterior samples against the held-out gauge observations, we directly assess the impact of the prior on posterior accuracy, thereby quantifying how much fine-tuning improves alignment with real-world observations. We examine three representative high-impact events. These events were selected to span a wide range of precipitation regimes, including prolonged precipitation associated with the Meiyu front, heavy precipitation driven by landfalling typhoons, and localized convective extremes. The primary case, which occurred over Hubei Province, China, during the East Asian summer monsoon on 2 July 2016, is shown in Fig. 3; additional examples are provided in SI Figs. 12, 13.

To evaluate the effectiveness of PRIMER, we employ two standard performance metrics: the mean absolute error (MAE) and the continuous ranked probability score (CRPS), with the latter providing a probabilistic measure of an ensemble system’s accuracy (see Method “Evaluation metrics”). For each metric, we define a relative skill score, ΔM, as the non-negative error measure of the original precipitation dataset (ERA5 or IMERG) minus that of the posterior sample, so that positive values indicate reduced error and thus enhanced skill. All evaluations are conducted at a spatial resolution of 0.1∘, where ERA5, IMERG, and posterior samples are compared against gauge observations treated as ground truth.

As shown in Fig. 3c, f, the updated prior P⋆(x) substantially outperforms baseline priors derived from ERA5 and IMERG. The ensemble-mean ΔMAE decreases from 0.46 mm/hr for P(x|O_ERA5) to 0.14 mm/hr for P(x|O_ERA5); a similar improvement is observed in the IMERG case, where the ΔMAE decreases from 0.29 mm/hr to 0.14 mm/hr. These gains extend beyond ensemble means: across individual samples, ΔMAE values for P(x|O_ERA5) are consistently lower than those for P(x|O_ERA5). Figures 4a, b reveal widespread reductions in MAE, highlighting PRIMER’s ability to systematically correct biases inherent in the original datasets. Figure 4c, d shows pronounced reductions in CRPS, with deeper blue tones indicating substantial gains in probabilistic estimates. These results demonstrate that PRIMER captures the posterior distribution accurately, with the improvements confirmed as statistically significant by t-tests. In addition to PRIMER, we also evaluated baseline priors (P_ERA5(x) and P_IMERG(x)) as well as baseline methods BCSD-EQM (Bias correction and spatial disaggregation–equitable quantile mapping)48 and RM (random mixing)49 (for notes of two methods, refer to SI Section 4). These baselines were subjected to the same application. As detailed in the SI Figs. 6–8, PRIMER generally outperforms these baselines. What’s more, the largest improvements are observed in the Sichuan Basin and the Pearl River Delta—regions with dense populations and strong economic activity. We further analyzed the correlation between gauge density and performance improvement (SI Section 5). Although a positive trend is apparent, the correlation is not statistically significant, indicating that PRIMER delivers relatively spatially consistent improvements irrespective of local gauge density.

Beyond reducing pointwise error, PRIMER also enhances the physical realism of existing precipitation datasets. To comprehensively evaluate the performance of PRIMER, we adopt two complementary perspectives: the member view and the envelope view. The member view analyzes statistics from a single sample, representing one physically plausible realization. In contrast, the envelope is constructed by selecting, at each gauge location, for a given event, the maximum precipitation value across 50 posterior samples. As illustrated in Fig. 5a, both P(x|O_ERA5) and P(x|O_IMERG) more accurately reproduce the frequency distribution of precipitation, particularly at higher intensities. Both perspectives reveal improvements in the representation of heavy precipitation tails compared to the existing datasets, underscoring PRIMER’s capacity to detect high-impact precipitation events that are often underrepresented in original products. Improvements in spatial structure are further quantified using PCCs with respect to gauge observations (Fig. 5b). P(x|O_ERA5) and P(x|O_IMERG) show markedly enhanced structural agreement relative to existing datasets, suggesting that PRIMER not only reduces local biases but also restores spatial coherence. While various methods have been proposed to assess spatial organization and feature propagation50,51, we employ a simplified yet informative diagnostic based on two-dimensional spatial lagged correlation coefficient (Method “Evaluation tool”, Fig. 5c). Physically, this correlation characterizes how anomalies at a reference point are spatially linked to those at surrounding locations, thereby revealing key features of precipitation system organization. We approximate the 0.6 correlation contour with an ellipse and extract two geometric descriptors: the focal length (F), indicative of spatial extent, and the orientation (O), which capture the dominant directional alignment. Results show that both P(x|O_ERA5) and P(x|O_IMERG) produce orientations that are more consistent with precipitation reference orientations derived from gauge observations, indicating improved spatial alignment. In terms of focal length, P(x|O_ERA5) exhibits a clear reduction, while P(x|O_IMERG) shows no substantial improvement. These results demonstrate PRIMER’s effectiveness in correcting the spatial anisotropy of precipitation systems.

Generalization test
PRIMER is not only effective for existing precipitation datasets, but also exhibits a certain degree of generalization. Figure 6 illustrates PRIMER’s ability to correct biases in previously unseen operational precipitation forecasts, using the ECMWF High-Resolution Forecast (HRES) as a representative example52. Despite never being trained on HRES, PRIMER successfully corrects systematic biases in a typical precipitation event caused by typhoon landing (Fig. 6a, e). The ensemble mean of P(x|O_HRES) (Fig. 6b, f) aligns with HRES, while each member (Fig. 6c, g) captures a diverse range of physically plausible scenarios, reflecting the model’s ability to encode meaningful uncertainty. Maps of ΔCRPS (Fig. 6d, h) with widespread positive values (blue dots) indicate that PRIMER produces a reliable ensemble system for HRES. These improvements arise from the Bayesian posterior sampling mechanism. By drawing samples from P(x|O_HRES), we effectively use the learned prior distribution P⋆(x)—which has been calibrated to match gauge statistics—to adjust the original HRES forecasts. To illustrate these benefits more intuitively, we present time series at two representative gauge locations (Fig. 6i, j). The ensemble envelopes generated by PRIMER closely track observed precipitation peaks, confirming that the HRES guidance is effectively incorporated. Occasional deviations arise when the HRES forecasts and the local gauge observations exhibit divergent trends—for instance, a slower decrease in HRES versus a sharper observed decline (panel i at +27h), or mismatched peaks (panel j at +27h and +42h). In such cases, the posterior may not fully capture the observed variability.

To assess PRIMER under a future scenario, we selected a model simulation output from CMIP6 for 205053, when elevated CO2 forcing is expected to alter precipitation regimes. Hourly precipitation was downscaled to 0.1∘ with PRIMER and compared against the raw model output. As shown in SI Fig. 17, the domain-mean precipitation curves from the raw output and the downscaled fields remain closely aligned, indicating that PRIMER preserves large-scale variability while adding fine-scale structure under a shifted climate state. Taken together, these results empirically suggest the broader utility of PRIMER as a foundation model for downstream applications without additional retraining (zero-shot adaptation).

Methods
Problem formulation
A general formulation of the precipitation data fusion task involves two key components: (1) constructing an informative prior distribution, and (2) performing posterior inference given new observations. Let x denote the target precipitation field. Different data sources—including gridded products such as satellite retrievals and reanalysis, as well as sparse gauge observations—provide multiple versions of x, each with varying spatial coverage and accuracy. Our goal is to effectively leverage these heterogeneous sources to construct a unified prior P(x). This prior plays a central role, as it is expected to integrate statistical characteristics of each source through a balanced fusion. A key innovation of this work lies in the design of a principled framework for modeling such a prior.

Once an informative prior is established, posterior inference is conducted as new observational evidence O becomes available. Posterior distribution P(x|O) can be factored into two components: the prior distribution P(x), and the likelihood P(O|x). Another innovation of our work is the effective implementation of posterior inference that balances the prior and the observations, ensuring the inferred precipitation field reflects both the climatological variability and the specific constraints provided by O. This Bayesian framework naturally enables various downstream applications, such as super-resolution by conditioning on coarse data, bias correction by conditioning on biased estimates, and optimal interpolation by jointly conditioning on observations and background (Fig. 1a).

Preliminary on diffusion models
To construct a prior, we employ score-based diffusion models. To enable the model to distinguish between sources during training, we associate each sample with a corresponding entity embedding e_i (e_1=(1,0,0), e_2=(0,1,0), e_3=(0,0,1))54, which is injected into the model. This embedding functions as a source identifier, enabling the model to learn distinct priors for different data sources. Specifically, e_1 corresponds to ERA5, e_2 to IMERG, and e_3 to gauge observations. Here, we first outline the foundations of the traditional diffusion framework before extending its conceptual scope. The forward diffusion process evolves the data distribution into a tractable Gaussian through a stochastic differential equation (SDE)33,34,55,56:
dx_t = f(x_t, t)dt + g(t)dW_t, (1)
where x_t is the state at time t, f(x_t, t) is the drift function, and W_t is a standard Wiener process. To generate samples, we solve the reverse-time SDE55,57:
dx_t = [f(x_t, t) - g^2(t) ∇_x_t log P_θ(x_t|e_i)] dt + g(t)dW_t, (2)
where the score function ∇_x_t log P_θ(x_t|e_i) denotes the gradient of the log-density. Since this score is intractable, we approximate it using a neural network f_θ.

PRIMER
Traditional diffusion models typically rely heavily on U-Net architectures44, which require inputs and outputs to be uniformly gridded data with fixed resolution. This architectural constraint limits their flexibility, particularly when processing discrete, sparse gauge observations. PRIMER utilizes a framework inspired by recent theoretical advances58–61, which generalizes diffusion models from finite-dimensional Euclidean space to an infinite-dimensional Hilbert space H, as illustrated in SI Fig. 1. In this setting, each element x is a function x: R^n -> R^d, where R^n denotes coordinates and R^d represents physical quantities. Both dense gridded data and sparse gauge observations are treated as partial realizations of an underlying function, allowing PRIMER to natively integrate heterogeneous records. Following ref. 58, we define H as L^2([0,1]^n -> R^d), where L^2 denotes the space of functions f such that integral |f(x)|^2 dx < infinity. The rationale for the name PRIMER is discussed in SI Section 1.

Mollification. While tempting, using white noise in the forward diffusion process poses a fundamental issue. Let ϵ(c) be a white noise where each c is sampled independently from N(0,1). For ϵ to lie in the Hilbert space H, it must be square-integrable. However, ϵ(c) violates this, as its norm diverges. To address this, PRIMER applies a Gaussian kernel k to mollify the noise: ξ(c) = (k*ϵ)(c) = integral k(c - c')ϵ(c')dc'. The resulting smoothed noise is square-integrable and thus belongs to H, as rigorously proven in SI Section 2.2. Similarly, PRIMER also mollifies x_0, which ensures that L_x_0 inherits the same smoothness properties. In practice, this operation is implemented efficiently using discrete Fourier transformations (DFT). In Fourier space, mollification corresponds to: ϵ(ω) = e^kωk^2t ξ(ω), where ω denotes the frequency, and t = σ^2/2, with σ being the standard deviation of kernel k (a detailed derivation is provided in SI Section 2.3). Directly applying the inverse transformation is often numerically unstable; thus, we employ the Wiener filter, defined as58,62: ϵ(ω) = e^-kωk^2t ξ(ω) / (e^-2kωk^2t + δ^2), where δ is a small positive regularization parameter.

Network architecture. Neural Operators are capable of learning a map between two functional spaces61,63–65. Neural operators achieve discretization invariance by learning integral kernels parameterized via neural networks. Specifically, for an input function x: R^n -> R^d, with observations at m distinct spatial locations, the operator K(x; θ) is defined as:
(K(x; θ)x)(c) = integral κ_θ(c, b, x(c), x(b))x(b)db, (3)
where κ_θ is a neural network that captures complex non-local dependencies. PRIMER implements a hybrid multiscale architecture that synthesizes the strengths of Neural Operators and convolutional networks. PRIMER first processes the input x, together with their corresponding locations c using a series of SparseConvResBlocks, which primarily employ sparse depthwise convolutions66, producing updated features, where D >> d. This embedding step projects low-dimensional input features into a higher-dimensional space, a crucial operation that enables the model to capture richer representations. For the motivation behind SparseConvResBlock, see SI Section 2.6. Since the features lie on an irregular set of discrete locations, we project them onto a coarse regular grid based on their spatial coordinates (see SI code 1). This transformation aligns the features to a regular grid, capturing the multiscale structure of the gridded data. After the SparseConvResBlocks, the processed grid features are reprojected to these coordinates via bilinear interpolation, yielding a feature tensor. Finally, a subsequent series of SparseConvResBlocks are applied to produce the final output tensor. Details of the network are provided in SI Section 2.5.

Model training. The model is optimized by minimizing a simplified denoising objective33,55,58 (derivation provided in SI Section 2.4):
L = E_t [|f_θ(x_t, t, e_i) - ξ|^2_H], (4)
where x_t denotes the noisy input at time step t, e_i represents the entity embedding, ξ is the ground-truth noise, and |.|_H denotes the loss norm defined in Hilbert space H. We adopt a two-stage training procedure. In Stage 1, the model is jointly trained on ERA5 (e_1) and IMERG (e_2). In Stage 2, we specialize the pretrained model to sparse gauge observations (e_3), following a strategy akin to DreamBooth67. Specifically, we fine-tune the model using a shared-weight strategy, where training samples are proportionally drawn from multiple data sources. The total loss is computed as:
L_fine-tuning = α_1 L_ERA5 + α_2 L_IMERG + α_3 L_gauge, (5)
with weights α_1 = 0.1, α_2 = 0.4, and α_3 = 0.5. Assigning α_3 = 0.5 prevents catastrophic forgetting of ERA5 and IMERG knowledge while ensuring strong gauge influence68. Among the remaining weights, IMERG (α_2 = 0.4) is favored over ERA5 (α_1 = 0.1) given its finer resolution. Although not optimized through exhaustive search, this empirical configuration preserves climatological priors while adapting to high-fidelity signals, thereby grounding the generative manifold in real-world observations.

The full training and inference pipelines are summarized in SI Algorithm 1 and SI Algorithm 2, with an overview schematic shown in SI Fig. 2. For the configuration of the hyperparameters, see SI Section 2.7.

Posterior sampling
In tasks such as bias correction, downscaling, and optimal interpolation, the objective is to infer an unknown target state x given observations O. PRIMER enables the incorporation of prior knowledge through a prior P(x), facilitating posterior inference via Bayes’ theorem: P(x|O) ∝ P(O|x)P(x). The standard reverse-time SDE can be modified to sample from the posterior distribution, yielding the following reverse diffusion process:
dx_t = [f(x_t, t) - g^2(t) (∇_x_t log P_θ(x_t|e_i) + ∇_x_t log P_θ(O|e_i, x_t))] dt + g(t)dW_t, (6)
This formulation requires two key components: the time-dependent score of the prior, ∇_x_t log P_θ(x_t|e_i), and the gradient of the log-likelihood, ∇_x_t log P_θ(O|e_i, x_t), which captures the observational constraints. Several studies have proposed various strategies to address posterior sampling within the diffusion framework37,38,69. In light of the characteristics of our problem setting, we adopt two representative approaches: Inpainting70–72 and SDEdit73.

Inpainting reconstructs unobserved regions by conditioning on partial observations O. A binary mask m indicates observed entries (m=1 if observed). At each reverse-time step t, a denoised estimate x_t is first computed. To enforce consistency with known observations, we blend the latent state using
x_t = m ⊙ q(x_t|O) + (1-m) ⊙ x_t,
where ⊙ denotes element-wise multiplication. The term q(x_t|O) is constructed by applying the same forward noise process to O; that is, for each observed entry, we simulate its noisy counterpart at step t under the forward SDE. This blending operation preserves observed values while allowing the model to impute missing regions, approximating the posterior distribution p(x|O). SDEdit can be viewed as a special case of inpainting where the entire input field is treated as observed, i.e., m=1. However, a key distinction lies in its use of a noise level parameter τ, which determines the strength of forward noise applied to the input before denoising. This parameter controls the extent to which the model is allowed to deviate from the original input, balancing fidelity and diversity. To select an appropriate τ, we conduct a sensitivity analysis on IMERG for 13 June 2016 at 23:00 UTC. For each noise level from 0.1 to 0.9 in steps of 0.1, we generate an ensemble of 50 samples from posterior P(x|O_IMERG) and compute both the RMSE and CRPS over 50 repeated subsampling trials, each selecting 10 members randomly. As shown in SI Fig. 4, performance improves with increasing τ up to around 0.6, beyond which both RMSE and CRPS begin to deteriorate. This suggests an optimal trade-off at 0.6 noise levels, where PRIMER maintains consistency with observations while capturing high-frequency variability.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Dear Sun et al.,\nCongratulations on writing this high-quality manuscript. The results are certainly noteworthy, interesting, and potentially\nuseful to the progression of more reliable weather forecasts. Please consider the comments below to further improve the\nwork.\nLine 77: “potentially” seems like an odd word choice here. Consider replacing with “often”\nLine 106: “Bayesian methods provide the most principled solution” Unless you can empirically show this or provide a source\nthat does so, I suggest softening this claim.\nLine 115: The text could benefit from a brief introduction of what diffusion models are and how they work before describing\nhow they have been applied to different tasks\nFigure 3: It is unclear in panel b and e what the line contours represent and what the shaded contours represent (The same\nissue applies to Figure 6)\nFigure 4: It would be useful to run significance tests to see if the changes in CRPS and MAE are significant\nFigure 6 and D15: The rain gauge measurements appear to fall outside of the PRIMER envelope several times. This should\nbe noted in the text.\nThe text could benefit from describing the noise scheduler used for the diffusion model and a brief discussion that justifies\nthe choice in noise scheduler\nLine 819: It would be useful to know the reasoning behind the selected alpha values (.1, .4, .5)\nLine 969: Should be “… we use a dataset from Shen et al. …”\nLine 1294: There is a typo here. “denosing” should be “denoising”\nIt could also be worth discussing why the experiments only examine lead times out to two days. Would this work at longer\nlead times and at what lead time would this no longer be effective?\nIt could be useful to provide the number of training epochs that were used and the possible variation of performance across\nepochs\nFigure 4: There seems to be a correlation between the change in MAE and the change in CRPS. This indicates that\nPRIMER is more effective in some stations compared to others. I suggest including some discussion that could hypothesize\npossible reasons for that happening.\nBest regards,\nTim Higgins\n(Remarks on code availability)\nThe code is clear and well-organized. The experiments are reproducible and I was able to run the code.", "reviewer_2": "Review of manuscript NCOMMS-25-51999-T entitled \"Fusion of multi-source precipitation records via coordinate-based\ngenerative models\".\nThe authors present a novel generative AI model that is capable of fusing multi-source precipitation data, including gridded\nsatellite and reanalysis data as well as irregular point observations from automatic weather stations. The model is designed\nto produce continuous precipitation fields that can be used for downscaling, interpolation, and bias correction tasks. The\ngeneral principle is to leverage score-based diffusion to learn a prior distribution of precipitation fields and then condition this\nprior on the available data to produce potentially higher-resolution precipitation fields. The model is evaluated on several\ndatasets, including ERA5 and IMERG, and shows promising results in terms of accuracy and flexibility.\nThe manuscript is well-written and structured, with a clear explanation of the model architecture in the method section. The\nmain document shows detailed results, but is hard to follow in a first read due to the complex methodology that is not\nsufficiently broken down to a digestible level while keeping the most important details intact.\nThe display material is generally of high quality, although some improvements could be made to the figures for better clarity\nand intuitiveness, especially for the general reader interested in applying improved precipitation observations in a practical\nsetting.\nThe methodology is novel and the results are promising. The methodical framework provides a high level of flexibility and\nadaptability, but the paper lacks a thorough comparison to existing methods (even though these may lack some key features\nlike producing continuous fields), which makes it difficult to assess the real performance gain of the proposed approach. The\nauthors should also clarify some methodological details, such as the independence of the training and validation data, and\nprovide more context on how their model compares to established precipitation fusion, downscaling and bias correction\nmethods.\nMy overall recommendation is to accept the paper after major revisions. The results are very promising and the method has\nthe potential to be extremely impactful, but the authors should address the comments below to improve the clarity and\nrobustness of their claims.\nMajor comments:\n- The results look promising, but it is hard to assess the real performance gain with a baseline other than the raw input\ndatasets. The authors should compare and or discuss their results with respect to baseline models, e.g. statistical\ndownscaling and bias correction or other statistical reconstruction and fusion methods [Lorenz et al., 2021; Yan et al., 2021].\nThere are deterministic AI models like MetNet3 that also make use of sparse observations, both as input and target\n[Andrychowicz et al., 2023]. MetNet3 successfully integrates sparse observations using coordinate-based encoding and\ncross-attention mechanisms, demonstrating how to maintain both local accuracy from point observations and large-scale\nspatial coherence. At least a discussion of such a similar approach, would be particularly valuable because it processes the\nsame types of multi-source data (gridded data + irregular point observations) and has established performance benchmarks,\nalthough for forecasting. It will allow readers to assess whether PRIMER's generative approach producing continuous fields\noffers clear advantages over MetNet3's deterministic neural network architecture for precipitation data fusion tasks on fine\ngrids.\n- The authors claim that PRIMER can be used for downscaling, interpolation, and bias correction, yet the results mainly\nfocus on bias correction. While it is clear that the model can be used for producing higher resolution data, downscaling\nresults are shown only for ERA5 enhancing the spatial resolution by a factor of 2.5. It would be good to see how the model\nperforms for different spatial resolutions, e.g. 2.5km, 5km, 10km, 20km, etc. and how the performance changes with respect\nto the resolution and training data resolution to support the claim that PRIMER performs well for downscaling. Similarly, the\nauthors should show how the model performs for interpolation tasks, e.g. by conditioning only on a hold-out set of AWS data.\n- The generalization claims are based on testing only one forecast model (HRES) from the same organization that produces\nERA5, and only one typhoon case study. This is evidence for claiming \"zero-shot generalization\", but it is not sufficient to\ndemonstrate robustness across different models and scenarios.\n- Since the model is conditioned on instantaneous snapshots of precipitation and produces random realizations of the\nposterior distribution, it is not clear how the model can ensure temporal consistency in the output. This could be improved by\nenforcing it in the training process and it should at least be analyzed in more detail. This is a crucial aspect for applications\nin hydrology and meteorology, where temporal coherence of single ensemble members is often required.\n- The feature that PRIMER is producing continuous fields should be stressed and analyzed more, e.g. by producing different\nresolution outputs and comparing them in terms of the correct space-time integration at this resolution. Here, I only see point\nmeasurements as a reference, but I also have to acknowledge that it is hard to understand the precise utilization of the\ngauge data. Especially since AWS data is referenced as both a gridded dataset and point observations. The authors should\nclarify in a more transparent manner how the AWS data is used in the model and how it is utilized for training and validation.\n- The authors should clarify the independence of the IMERG data used for training and validation. Since IMERG\nincorporates rain gauge observations in its calibration process, it is crucial to ensure that the AWS stations used for fine-\ntuning were not part of IMERG's production. This may lead to data leakage and potentially inflated performance metrics. I\nsee that IMERG is compared to AWS with imperfect scores, but this does not completely rule out the possibility of data\nleakage.\n- Figure C5 shows highly uneven gauge distribution, but there's no analysis of how this spatial bias affects validation results.\nPerformance improvements may simply reflect dense gauge coverage in certain regions rather than model skill. It seems like\nthe model performs better in regions with more gauges.\n- The validation is mainly based on point measurements, which may not fully capture the model's performance in producing\ncontinuous fields. It would be beneficial to generate outputs at multiple resolutions and show consistent performance across\nthese resolutions. This would validate the model's ability to produce continuous fields that maintain realistic precipitation\nscaling laws across all spatial scales [Fabry, 1996].\n- In terms of the general applicability of the framework I also have a circular validation concern. The model is fine-tuned on\ngauge observations and then validated against gauge observations. This circular validation approach likely inflates\nperformance metrics and limits demonstration of true generalization capability. Independent validation datasets not used in\nany training stage, such as weather radar data, would be essential for proper validation.\n- The authors could show how the produced ensembles are calibrated apart from a mean and standard deviation\ncomparison. A good comparison where the authors show how the ensemble spread is calibrated with respect to the\nobserved precipitation distribution is in Harris et al. (2023) or Glawion et al. (2023).\nMinor comments:\n- The authors should mention how PRIMER compares to an explicit traditional Bayesian inference approach, especially in\nterms of computational efficiency and scalability. This would provide a clearer context for the advantages of the proposed\ngenerative model.\n- Line 133: It might be unclear to readers what \"high dimensionality of precipitation fields\" means. While fields are 2-3\ndimensional for most applications, I think you want to refer to the number of dimensions needed to represent the data\naccurately, or you want to refer to the number of grid points in the fields each generating a dimension when looking at the\nproblem from a machine learning perspective. Please clarify.\n- Fig 1: Be aware that the first figure readers see should be extremely intuitive and easy to understand. Many readers might\neven first look at Fig. 1 and then read the model framework description. I think there is a large potential to improve this figure.\nSome suggestions:\n- a):\n- I think a top view of the discrete stations would be more intuitive and would emphasize the irregular distribution.\n- Visualize how PRIMER produces continuous fields and how observations can be queried from the continuous field.\n- It is still unclear to me where AWS data is used as sparse gridded data and where is is used with precise coordinates.\nMaybe you can clarify this in the figure using appropriate visualization.\n- b):\n- It takes some time to decode the color scheme. Maybe use a more direct connection between stages and gridded/irregular\ndata.\n- The training stages are important, but it would be easier to understand the model if inference was visualized first. Maybe\nyou can show the data characteristics and inference combined in a)?\n- Section 2.2: The description of the relative error improvement is hard to grasp. It requires re-reading sections several times\nto get the performances in order. I also think that absolute errors should be reported in addition to relative errors, as they are\nmore intuitive to understand. Please clarify the description and consider adding absolute errors.\n- Fig. 4:\n- Use histograms to visualize the distribution of the relative errors.\n- unclear what baseline means in this figure.\n- Line 648 f.: The outlook mainly talks about future, unexplored model development. It would be good to stress the\napplication potential of the current model more, e.g. by discussing how the model can be used in practice.\n- Animations are a powerful tool to demonstrate spatial and temporal coherence of rainfall fields\nTechnical corrections:\nLine 82: \"avearge\" should be \"average\"\nLine 140: \"here after\" should be \"hereafter\"\nLine 144: \"Gaussain noise\" should be \"Gaussian noise\"\nFig D7: Some colorbar labels are missing\nReferences:\nYan, J., Li, F., Bárdossy, A., and Tao, T.: Conditional simulation of spatial rainfall fields using random mixing: a study that\nimplements full control over the stochastic process, Hydrol. Earth Syst. Sci., 25, 3819–3835"}, "author_response": ""}]

Decision letter