the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 31 Aug 2026
MAROR: a multipoint autoregressive network for forcing-driven, uncertainty-calibrated reconstruction of hourly coastal sea level from 1940 to the present
Abstract. Tide gauge records are the primary source for coastal sea level research, yet they are short, gappy, and affected by instrumental errors, which limits the study of extreme events and their long-term changes. We present MAROR (Multipoint Autoregressive Reconstruction Of sea-level Records), a neural network trained on tide gauge records that reconstructs gap-free, hourly coastal sea level from 1940 to the present, driven only by atmospheric forcing and, when available, satellite altimetry. Building on the HIDRA3, a convolutional neural network designed for short-term predictions, MAROR turns a short-range forecast model into a long-range autoregressive reconstruction engine, with an exposure-bias correction that removes drift and a generative ensemble that yields a calibrated, per-tide-gauge uncertainty band. We evaluate it on 373 tide gauges across eleven regions spanning a wide range of surge regimes, from extratropical shelves to tropical cyclones and western boundary currents. MAROR reconstructs the hourly storm surge with a median correlation of 0.88, reproduces the most extreme events, and outperforms both a hydrodynamic hindcast and a state-of-the-art deep-learning model.
- Preprint
(5119 KB) - Metadata XML
-
Supplement
(376 KB) - BibTeX
- EndNote
Status: final response (author comments only)
- RC1: 'Comment on egusphere-2026-5073', Anonymous Referee #1, 09 Sep 2026
-
RC2: 'Comment on egusphere-2026-5073', Anonymous Referee #2, 22 Sep 2026
The paper proposes a neural network MAROR for reconstructing gap-free hourly coastal sea surface heights over very long periods, possibly spanning decades. The method is a very interesting reformulation of the HIDRA3 model, which has been initially designed for short-horizon SSH forecasting. While HIDRA3 relies on past SSH and predicts future 72h, the novel auto-regressive formulation is trained to recursively predict SSH values one-hour-at-a-time, solely by astronomic tide, atmospheric forcing and satellite altimerty where available. Bias correction is proposed to mitigate inevitable drifts in long-range forecasting setups. The experimental analysis is extensive and the demonstrated example of tide gauge anomalous performance detection is very interesting. I liked the experimental evaluation very much. If authors are planning to make the data formatted for the specific setup published in this paper freely available to the community, it would be great to include a link to the dataset bulk. I think the paper very good and brings a strong contribution. Nevertheless, I have several comments for the authors below ordered as they appear in the paper.
Page 3: There is a bit of confusion in the second paragraph. It stars by saying "...given the recent sea level at each tide gauge..." and later end "...so no observed sea level is ever ingested.." -- these are contradictory statements.
Page 5, Section 2.3.1.: The text implies general WxH values may be used for atmospheric encoding, but in fact they have to be quite specific, since the computational blocks reduce it to a 1x1 vector, right?
Page 5: "The number of output kernels..." --> I suggest to use a more descriptive term, "output feature channels". It tells the reader the output dimension and number of kernles is implied for the readers used to NN terminology.
Page 5: "...axis gives 2 · 512+2 · 64 = 1152 channels over 36 time steps". I suggest explictly stating that the operations produce a 36x1152 tensor.
Page 5: "feature vector a ∈ R^8192" --> To be consistent with the prior notation of strides and tensors, it might be worth considering the forms of type "a ∈ R^1x8192". This might be a bit unusual for readers coming from non-machine learning background, but on the other hand the machine-learning part of the community might be more used to this kind of notation, as it maps to actual datastructure sizes used in implementations. I let it to the authors to decide.
Page 5: "...SLA is enabled, the altimetry window ℓi ∈ R30 are concatenated and encoded by..." --> I suggest to write the dimensions of the tensor after concatenation? Afterall if the SLA is unavailable, a zero vector of equal size is concatenated? SLAs colocated with the tide gauge must be quite gappy, right? Could you give some statistics about the gappiness or percentage of missing data, i.e., the usual portion of hourly data available in he 30h window?
Page 7: "it produces many alternative predictions (ensemble members), all equally plausible, ...": Are these really equally plausible? I would think those closer to the distribution mode would be more plausible than those generated at the tails.
Line 185: The ensemble-generation head with random noise sampling is interesting and reminiscent of the GAN-based generative models. I suggest giving some additional insights to this in the experiments section. Does the r_i(,) change the spread depending on the atmosphere, or has it learned a more or less independent constant variance growth? An empirical insight would be valuable. Also, the paper mentions that non Gaussian pdfs can be learned by r_i. Could you please show whether these are in practice indeed substantially non-Gaussian? I would also suggest testing against simply predicting the variance as in HIDRA3 to validate the need for the sampling-based estimate.Section 2.4: The section is missing explicit definition of the loss functions used. I suggest writing them out at each phase. I also urge the authors to denote losses in Figure 1 if possible.
Section 2.4.2: The argument for stop gradient after each inferred hour goes against the established practice (and theory) of the recurrent networks. In fact, the theory shows that it would be beneficial to backprop over the entire chain, as it puts the network closer to test-time performance. However, the memory might become the issue given the large architecture (though, again, since the atmospheric encoder is likely well trained at Phase2, that part of the network might be kept frozen, which then reduces the memory consumption). If the primary reason for stop gradient is memory consumption, the paper should state that clearly and not place it as a convenient property of stop gradient decision. I also suggest avoiding the "stochastic spread", since it enters only at Phase 3. Again, a clear definition of the loss by an equation will clear things up.
Section 2.4.3: Please write out the CRPS loss and expose the parameters. Also some argument for using \alpha=0.95 should be provided (why not 0.5 or 0.9 or any other?).
Section 2.4.4: Are separate factors used on each station? It would be interesting to see their values and whether they vary substantially between the stations. Since a post-hoc inflation is trained on the held-out set, would it not be simpler to have the network estimate simply the variance, which should be easy to train on the held-out set? Also, how exactly is it trained? On just 72h reconstruction setup as the previous stages? I suggest writing out the cost/loss function for this training phase as well.
Section 2.4.5: "The columns of the first per-gauge dense layer ...". The term "columns" is very implementation-oriented, i..e, how data is structured for implementation and does not mean much to the reader. The "column" is where the SLA's are kept. So I'd suggest saying something like "the SLA" measurements input tensor...
Section 2.5.1, Cold start: The cold start is very interesting. Could you quantify how well it usually fits the tide-gauge measurements and how the reconstruction quality might change if tide-gauge data was used instead (assuming it's available)? This would give some practical insight into network's robustness to imperfect initialization.
Section 2.5.2: At inference time, 50 members are run to evolve in parallel. I understand that's useful for estimating the spread and also the mode of the distribution, i.e., most likely outcome. But I also think it might be very interesting to run the "mean member" -- right after the cold start, the network will predict a common mean SSH for the next hour, along with 50 members for that hour. If I understand correctly, only the 50 members continue to be propagated by their own predictions. But why not do the same for the prediction of the "mean"? And I would be curious to know how this single "mean member" compares in accuracy after long unrolling to the average computed over the 50 members.
Section 2.5.3: The calibration factor \beta is estimated on the held out set. It makes sense. Nevertheless, what is the reason for underestimation? Since the training data is rather huge, it seems overfitting should not be causing it. In that sense, it might be something to do with the CRPS loss. What happens when the \beta is estimated on training data instead of the held out set? Do values differ much?
Page 12: I suggest detailing exactly from which data the standard deviation is computed for outlier rejection at tide gauge.
Section 3.3: Given the satellite data might be quite sparse, could you plot a histogram of the measurements over the region of interest? What is the typical number of SLA satellite observations over the cell occupied by a tide gauge?
Page 13: Tide gauges are typically referenced to a local zero and not globally calibrated, e.g. among themselves nor with SLA. I wonder if the authors calibrated the tide gauges with the SLAs, or is the network expected to do that internally? Also, since a max radius of 300km was used for pairing SLA with a tide-gauge, what is the expected SSH difference typically observed for two points at 300km distance? Is the difference negligible?
Results: I generally like the experimental evaluation. Nevertheless, it is stated at several points that the tide gauge measurements do not enter the reconstruction. But I do not see why it would be beneficial to "always ignore measurements". Surely the best setup is that the reconstruction network ingests the measurements when available and uses its own predictions when they are unavailable. So I suggest to take the best MAROR version, then just ingest the SSH measurement in the past buffer when available. Does reconstruction performance further improve?
Citation: https://doi.org/10.5194/egusphere-2026-5073-RC2 - RC3: 'Comment on egusphere-2026-5073', Anonymous Referee #3, 27 Sep 2026
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 418 | 241 | 64 | 723 | 75 | 54 | 59 |
- HTML: 418
- PDF: 241
- XML: 64
- Total: 723
- Supplement: 75
- BibTeX: 54
- EndNote: 59
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Dear Authors,
My detailed comments and revision requests are provided in the attached PDF file.
Regards