4 Discussion 4.1 Interpretation In summary, conditioning a flow-matching model on PM 2.5 and meteorology yields spectra whose distribution, not only whose mean, tracks the observations, and a reduced GDE residual improves them further. Since this study is a benchmark on a two-station testbed, the value of these results lies less in the absolute scores than in what they reveal about model families, constraints and failure modes. Three results call for interpretation: the different verdicts of RMSE and CRPS, the contrast between constraint types, and the dependence of transfer on quantity and direction. Why ensembles beat regression on CRPS but not always on RMSE. Under the mass–number ambiguity, several distinct spectra are consistent with the same input. A squared-error regression returns their average, which RMSE rewards but which lacks the structure and spread of any single realization (Fig. 2 ). A generative ensemble spreads probability across the plausible spectra, which CRPS rewards. The two metrics therefore answer different questions. Gradient boosting’s lower RMSE at Hyytiälä and the generative models’ lower CRPS at both sites are consistent: the former is the better point estimate at that site, the latter the better description of what remains uncertain. Why the dynamical constraint helps and the integral constraints do not. An integral constraint fixes one summary number and leaves many spectra that satisfy it. In the nowcast, the mass reference also duplicates information already given as input. When the resolved spectrum and the instrument do not sample the same size range, forcing agreement can distort the spectrum. The GDE residual instead constrains how neighbouring bins and hours relate to each other, and it also widens the too-narrow FM ensemble. These observations are consistent with a regularizing effect of the dynamics, but the evidence is indirect: the constrained runs also have more parameters and a longer training budget. Why transferability depends on quantity and direction. Number and condensation sink weight the same spectrum differently. The sink is weighted towards larger particles, which tend to be regionally mixed and therefore shared between the two sites, whereas ultrafine particle number is dominated by the smallest particles, which reflect more local sources such as traffic at Kumpula and particle formation at Hyytiälä. We hypothesize that this explains why the sink transfers better from Kumpula to Hyytiälä. The boreal-to-urban failure is consistent with a missing source: the Hyytiälä record contains no analogue of the urban daytime small-particle enhancement, and a soft dynamical constraint cannot supply one. The present two-site comparison cannot separate the roles of emissions, formation, growth, transport and mixing, so these explanations remain hypotheses. A practical implication follows regardless: a PSD model should be judged by the quantity and uncertainty that its application needs, because full-spectrum correlation can remain high while the relevant integral degrades. 4.2 Limitations Several limitations follow from the deliberate choice of a small, controlled testbed, and should be read in that light. The evidence rests on two stations, one late-year holdout of about 60 days and one training run per configuration, so the bootstrap intervals do not include training-run variability; repeated training and holdouts from other seasons would address this. Two stations are enough to compare model families under identical conditions and to expose failure under transfer, but not to show how far neural networks can generalize when trained on many aerosol regimes. That question requires large, multi-station training data and is the natural next step. Every learned ensemble is too narrow, with 80 % intervals containing about half of the observations locally and fewer after transfer, and eight members limit how finely quantiles can be resolved. Post-hoc calibration, larger ensembles and guidance selected on validation data rather than the exploratory test-set sweep (Appendix D.4 ) are natural remedies. The fitting period of the shared PSD normalization has not been verified, so a scaling leak from test data is not yet ruled out, and quality flags were not applied (Appendix A ). The high-number episode analysis contains only 0, 3, 29 and 0 qualifying days in the four settings, too few for reliable probability verification (Appendix F.1 ). The physical interpretation is also incomplete. The GDE gain cannot yet be attributed to the residual itself, because the constrained runs have an auxiliary head and a longer training budget, and the sulfuric-acid proxy in the growth term depends on SO 2 , which is not among the stated inputs (Appendix C ); capacity- and budget-matched controls are needed. The integral constraints use references with different measurement support, and in the nowcast the pseudo-mass reference is the input itself, so a size-resolved conversion between DMPS, PM 2.5 and CPC support would be required to make them physically meaningful. Finally, the derived quantities are computed from the same DMPS spectra used for evaluation and so do not validate process physics independently, and the SOSAA-FP comparison uses windows that overlap training data and a single model configuration, so it cannot rank data-driven and process-based modelling in general. \conclusions We asked how much of the aerosol particle size distribution can be recovered from the concurrent PM 2.5 and meteorological record. Treating this nowcast as conditional generation, we benchmarked six generative models and three physical constraints on a controlled two-station testbed, with the aim of identifying which modelling choices matter rather than proposing a single best model. Four findings stand out. (1) Generative ensembles give better distributional skill than deterministic and climatological estimates: all six generative models have lower CRPS than gradient boosting at both sites, although gradient boosting remains competitive on point error. (2) The form of a physical constraint matters: a reduced GDE lowers CRPS by 11.5 % and 5.5 % and widens the ensemble, whereas integral constraints on mass or number alone degrade the spectra. (3) The routine record adds information beyond historical spectra: FM+GDE has the lowest CRPS of the learned models for ultrafine particle number and the condensation sink at both sites, with skill scores of 0.27–0.83 against the historical ensemble, but its intervals contain only about half of the observations. (4) Transferability depends on the quantity and the direction: the condensation sink transfers better than ultrafine particle number from the urban to the boreal site, whereas boreal-to-urban transfer fails for ultrafine particle number. Local PSD nowcasting from routine observations is therefore feasible today, and the testbed has shown which ingredients are worth carrying forward. These results raise the question of what would be required to nowcast PSDs at a station that measures only PM 2.5 . First, and most importantly, the method should be scaled up: trained on large data sets from many stations that cover contrasting aerosol regimes, neural networks can be tested where their main strength lies, in learning from diverse data and generalizing to new domains and out-of-distribution conditions. The boreal-to-urban failure observed here, after training on a single station, is exactly the situation that multi-station training is designed to overcome. Second, calibrated uncertainty is a prerequisite for any use of the intervals. Third, training-derived climatology could be represented explicitly, for example by starting generation from a historical PSD template instead of pure noise. Fourth, source-sensitive predictors such as traffic or emission indicators could reduce the mass–number ambiguity directly. Finally, capacity-matched ablations would clarify how much of the constraint’s benefit is physical. These improvements are complementary to process models: data-driven nowcasts can extend the observational record, while process models remain necessary to explain it. The evaluation framework presented here can be used to measure progress along each of these axes.
Appendix A Data preparation and evaluation splits A.1 Size grid and data assembly The DMPS archives at both stations report 61 diameter bins between 1.00 and 1000 nm. We use a fixed 50-bin subset between 2.82 and 794 nm, excluding the nine smallest archived bins (1.00–2.51 nm) and the two largest (891 and 1000 nm), and we keep the archived bin centres without interpolation. Hourly values are means of the native 10 min measurements, each weighted equally. ERA5 variables are taken from the grid cell nearest to each station. After assembly, each station has 61 368 hourly records with complete PM 2.5 , PSD and ERA5 values on a common time axis. A.2 Transformations Particle concentrations span several orders of magnitude across bins and over time, so the model operates on log-transformed values. For each bin, the concentration n is transformed to ℓ=log(1+n) and then scaled linearly to the interval [−1,1] , ℓ~=2ℓmax−ℓminℓ−ℓmin−1, (13) with fixed bin-wise constants ℓmin and ℓmax shared by both stations. The same scaling therefore applies at the training site and at the transfer target. Values outside the scaled range are kept rather than clipped. To return to physical units, the scaling is inverted and n=exp(ℓ)−1 is applied after negative values of ℓ are set to zero, which guarantees non-negative concentrations. PM 2.5 is transformed with log(1+x) , and positive, strongly skewed meteorological variables are also log-transformed. A.3 Examples and splits An example is a 24 h window of inputs together with the PSD of the same 24 hours (Sect. 1 ). Windows start at every hour, so consecutive windows overlap by 23 hours, and a window is used only if all 24 hours are valid. Missing hours are kept as gaps on the time axis rather than removed, so that every window spans exactly 24 consecutive clock hours. The data are split in time. All training windows end before 2 November 2024, 00:00 UTC, and the temporal holdout consists of the 1417 windows that start on or after this cutoff (60 days, 2 November–31 December 2024). Before the cutoff, 10 % of the eligible windows (6048 per station), stratified by month and hour of day, form an interpolation split, and the remaining 53 857 windows per station are used for training. Because the interpolation windows overlap with neighbouring training windows by up to 23 hours, the interpolation split measures how well a model reproduces the training period, not its skill on unseen periods. We use it only for the diagnostics in Appendix D.2 and for the matched comparison with SOSAA-FP (Sect. 3.5 ); all skill statements in the main text refer to the temporal holdout. No separate validation set or early stopping is used. Appendix C Physical constraints C.1 Resolved mass and number The resolved mass and number (Eqs. 3 and 5 ) are sums over the 50 bins of the predicted PSD, using logarithmic bin widths Δi=Δlog10Dp,i and a constant particle density ρ=1.5 g cm -3 . PM 2.5 is defined by the aerodynamic diameter, whereas the DMPS sizes particles by their electrical mobility. For spherical particles, the two diameters are approximately related by Dm≈Da(ρ/ρ0)−1/2 with ρ0=1 g cm -3 , neglecting shape and slip corrections. The PM 2.5 cut therefore corresponds to a mobility diameter of about 2.5μm×1.5−1/2≈2.0μ m, well above the 794 nm upper limit of the grid. All predicted bins thus lie inside the PM 2.5 size range, but the mass between 794 nm and about 2 μ m is missing from the resolved mass, which is why the mass term is used only as a soft reference. Both integral constraints compare logarithms of the resolved and the measured quantity through the Huber function, Huberδ(u)={u2/2,δ(∣u∣−δ/2),∣u∣≤δ,∣u∣>δ,δ=0.5, (17) where u is the difference of the two logarithms. In logarithmic space, a constant multiplicative mismatch, such as the missing mass above 794 nm, becomes a constant offset, and the linear branch of the Huber function prevents such systematic offsets from dominating the gradient. The reference scales M0 and N0 in Eqs. ( 4 ) and ( 5 ) make the logarithms well defined at zero concentration, and negative reference values, which can arise from instrument noise, are set to zero. The number constraint is evaluated only in hours with a valid CPC measurement. Each physical term is weighted so that its initial magnitude is comparable to that of the FM loss; for the mass term, this gives λm=0.05 , and the number and GDE weights were set by the same rule. A mass weight of 1.0 led to a degenerate solution that satisfied the mass constraint at the expense of a realistic spectrum, which illustrates why the mass term must remain soft. The CPC and the DMPS measure closely related but not identical quantities. Over the full record, their log-transformed total number concentrations correlate with Pearson r=0.988 at Hyytiälä (59 689 hours with both measurements) and r=0.963 at Kumpula (57 424 hours), with median DMPS/CPC ratios of 0.87 and 1.03, respectively. Possible contributors to the differences include particles below the DMPS lower limit, which the CPC counts, and differences in detection efficiency and sampling losses ( Wiedensohler et al., 2012 ) , which the number constraint does not model. C.2 Condensation and coagulation sinks The condensation sink (Eq. 12 ) is the rate at which a condensable vapor is lost to the particle population ( Kulmala et al., 2001 ) . The Fuchs–Sutugin correction accounts for the transition between the continuum regime, where vapor diffuses to large particles, and the free-molecular regime of small particles, βm(Dp)=1+(3α4+0.377)Kn+3α4Kn21+Kn,Kn=Dp2λv, (18) where λv is the mean free path of the vapor molecules and α the accommodation coefficient, set to one. For small particles, βmDp grows approximately with Dp2 , and for large particles approximately with Dp , so the sink is weighted towards larger particles. The coagulation sink of a particle of diameter d , CoagS(d)=∑iK(d,Dp,i)Ni, (19) is its loss rate by coagulation with the existing particles, where K is the Brownian coagulation kernel ( Kulmala et al., 2001 ; Seinfeld and Pandis, 2016 ) and Ni=niΔlog10Dp,i is the number concentration in bin i . Both sinks are evaluated at fixed reference conditions (293.15 K, 101325 Pa) for predicted and observed spectra alike, so that differences between them reflect the size distribution only. C.3 Reduced general dynamic equation For a size distribution n=dN/dlog10Dp at a fixed station, the continuity equation underlying the GDE constraint is ∂τ∂n+∂lnDp∂[DpGRn]=S−CoagSn−λdiln, (20) where τ is physical (clock) time. The second term on the left describes condensational growth as a flux through size space with velocity GR/Dp in lnDp . On the right, S collects sources and transport, the coagulation term removes particles that collide with the existing population, and λdil=(dH/dτ)+/H dilutes the population when the boundary layer of height H grows. The equation is reduced in three respects: coagulation enters only as a loss, not as a gain at larger sizes; deposition is not represented explicitly; and new particle formation appears only as the growth flux into the smallest resolved bin. The source term S is unconstrained in sign and can absorb these omitted processes, as well as the advection of different air masses past the station. The residual of Eq. ( 20 ) is evaluated on the predicted spectrum (Eq. 6 ). Time derivatives are computed by finite differences over the 24 hourly values, central in the interior and one-sided at the first and last hour, and size derivatives by upwind differences in lnDp , which follow the direction of growth. At the lower edge of the grid, the incoming flux is set by the growth velocity of the first bin; it represents growth into the observed range rather than a nucleation rate. All rates are expressed per hour. The growth rate and the source term are predicted by an auxiliary network head from the conditioning only, so that they do not inherit noise from the intermediate generative state. The growth rate is constrained to be non-negative and is informed by a sulfuric-acid proxy proportional to SO 2 concentration times radiation divided by the condensation sink ( Petäjä et al., 2009 ) ; the dependence of this proxy on SO 2 , which is not among the model inputs, is discussed in Sect. 4 . The head also predicts the log residual variance sh,i of Eq. ( 7 ), which is bounded to [−5,5] , i.e. to variances between 6.7×10−3 and 1.5×102 , so that the model can neither ignore a residual entirely nor overweight it without limit. Because S and sh,i are flexible, a low residual shows consistency with Eq. ( 20 ) but does not uniquely identify the growth or source rates. C.4 Coupling the constraints to flow matching During training, the network sees partly noised states xt , whereas the physical relationships apply to clean spectra. Inverting the straight path of Eq. ( 14 ) gives the clean-sample estimate implied by the current velocity prediction, x^1=xt+(1−t)vθ, (21) which equals x1 exactly when the predicted velocity is correct. The constraints are evaluated on x^1 after conversion to physical units, and their gradients flow back into vθ . Near the noise end ( t→0 ), x^1 is an extrapolation from an almost random state, so the weights t2 and t3 of Eq. ( 2 ) suppress the constraints there; the GDE term, which contains derivatives that amplify noise, receives the steeper weight. The constraints act only on the training objective; generation at inference proceeds exactly as for unconstrained FM. All physical ablations use FM as the base model, and every constrained variant uses this time weighting.