Report of Reviewer #1

[MarieEtienne]

Associate Editor: MP Etienne
Reviewer : DqP1 (chose to lift his/her anonymity / remain anonymous)

Reviewer: Reviewing history

• Paper submitted Dec, 19th 2023
• Reviewer invited Feb, 14th 2024
• Review 1 received Apr, 10 2024,
• Paper revised July, 8th 2024
• Paper conditionally accepted July, 10th

First Round

The comment from reviewer 1 are publicly available on Openreview

Summary Of Contributions:

The article is very well-written and well-organized. I thoroughly enjoyed reading it. The author provides a comprehensive description of spatiotemporal modeling of high-dimensional fire occurrences and sizes within the Bayesian inferential paradigm, utilizing INLA for inference. Fire occurrences are modeled using the Poisson response with a random intensity measure of log-Gaussian Cox types, while mark sizes are modeled using the gamma distribution with spatially varying means incorporating spatiotemporal structures. The author also suggests a posterior simulation-based approach for future projections under two commonly known climate scenarios. Overall, the flow of the manuscript is handled well, and the content is well-presented.

Strengths And Weaknesses:

I had some difficulty understanding the novelty of the paper compared to Koh et al. (2023) and Pimout et al. (2021). I believe that the model proposed herein could be considered a sub-model of Koh et al. (2023). The one novelty that I could identify was in simulating future wildfires while combining the posteriors of model parameters with climate model output. This is a great idea, though I feel that this can be further improved, especially in the way the predictor variables FWI and FR are created for future projections. I'm wondering whether FWI and FR cannot be predicted or calculated based on reanalysis data or using available historical data, and then using these predicted FWI and FR when simulating future wildfires? Please disregard if this involves too many technicalities. However, assuming FWI and FR surfaces remain constant above historical ranges while projecting into the future might not be the most accurate approach, especially when using scenarios like RCP8.5. These clarifications should be clearly stated in the manuscript and/or in the discussion section if your goal is to do it as future research.

Thank you for these comments. We have added some text in the revised version to better explain the differences with the paper of Pimont et al (2021) and the paper of Koh et al (2023).

You are totally correct that our model can be viewed as a submodel of Koh et al. (2023). However, we here perform modeling for a study region that has never been modeled before and for which climate and also land-use conditions are quite different from the classical study area in the historical wildfire-prone area in Southeastern France used for the original Firelihood model and also the model of Koh et al (2012). A particularity is that there are generally fewer wildfire occurrences in the region studied here.

Moreover, and as well identified in your remark, we explore the combination of posterior simulation and climate model output to provide projections of wildfire risk under climate change for this region.
Also, a methodological novelty is that we show how these projections can be smoothed using R-INLA.

Regarding the extrapolation of the predictor variables, we fully agree with you that this might not be the most accurate approach, but for lack of better solutions, we consider it the best possible approach.

However, regarding the FWI, various results in already published work have shown that the FWI effect becomes saturated at very high levels, i.e., it does not further increase even if the FWI values increase, as for example discussed in the Firelihood papers Pimont et al (2021) or Koh et al (2023). The reasons for this can be found in the construction of the FWI as an index for Canadian forests, where the influence of variables such as relatively high wind speeds is stronger than for France. Therefore, a constant extrapolation does not seem problematic in this case. There certainly is higher variance in estimated spline functions for the high FWI values, and constructing novel models that appropriately address this issues is still an open research question and is too complex to be addressed here; see also the written comment on this problem published in the context the JRSS ``The First Discussion Meeting on Statistical aspects of climate change" by Legrand and Opitz (2023), that is written by some of the authors of the current manuscript. To better highlight this issue and give some directions for future work as suggested, we have added the sentence in the discussion: One limitation of our approach for the simulation of future wildfire activities is the constant extrapolation of the function $f_{FWI}$ .* Various studies have demonstrated that the FWI effect becomes saturated at very high levels, suggesting that a constant extrapolation may not be problematic. However, there is a higher variance in the estimated spline functions for the high values of FWI, and developing novel models that appropriately address this issue is the subject of future research Legrand and Opitz (2023).*

For the extrapolation of the FA, as of today, there do not yet exist any reliable projections of future FA in the study area, and therefore our solution is the best one we could think of in these circumstances. In fact, projections of future FA (or, at least, the creation of some plausible "storylines" how future FA could look like) is one of the goals of the FIRE-RES project, under which this study was conducted.}

Another concern of mine is related to jointly modeling counts and size processes by the use of sharing random effects. You did mention the possibility of sharing some of the spatial random effects between the occurrence model and the size model, similar to Koh et al. (2023). While computational cost is one consideration, I believe you might see improved performance due to strong dependencies between the two processes and enhanced uncertainty estimations. This is something I will recommend to worth give a try to a datasets that you can handle computationally. Again, this is the trade-off between model flexibility and computational costs.

Indeed, linking the two components, as done in \cite{Koh2021}, would result in a more sophisticated model, which would be more complex to construct and more difficult to estimate. This would be the price to pay in order to obtain possibly more accurate (and narrower) uncertainty estimations. However, in line with one of the scope of the journal, the aim of this study was rather to present a simple flexible model that can be easily applied with other datasets. Moreover, even if the estimation is conducted independently, the simulation of sizes is conditioned on the simulation of the occurrences, thereby introducing a certain degree of dependence between the two.

Furthermore, the model proposed here includes spatial and temporal random effects in an additive fashion, hence lacking space-time interaction. I think this effect is evident from Figure 8, where the spatial pattern of estimated counts and burnt areas for future projections appears very similar to those in historical periods. I believe this could be an artifact of using a model that lacks space-time interactions. For instance, the models proposed herein account for temporal dependence, which is independent a priori from spatial effects. A proper spatiotemporal model would incorporate non-separable space-time structures to effectively capture space-time interactions. Therefore, I will be careful while drawing any such conclusions in future projections with models without space-time interaction. It's also possible that the observed similarity in spatial patterns is influenced by the local climate and specific to the study region.

The purpose of the spatial and temporal effects in the model is there to capture certain spatial and temporal patterns not well captured by the information provided by FWI and FA covariates. Note that the values of these covariates are not "separable" in space and time. Climatologists and fire experts are usually quite skeptical about including relatively complex random effects into models, since the physical covariates (FA, FWI...) should already be able to explain spatiotemporal patterns as far as possible. The current model formulation is a good compromise since clear explanations can be given to the spatial and temporal effect, for example in terms of different land-use, land-cover and wildfire management for the spatial effect, and in terms of biases of the FWI for appropriately representing fuel moisture properties across the season for the temporal effect. Such explanations were given in various other papers, such as \cite{Castel2023,Koh2021,Pimont2021}. In the paper, we have added some text to explain this in the Discussion section.

Requested Changes:

Major comments

Section 3.1: Do you have any insight into how these subsampling schemes affect spatial dependencies in your model? I agree that the use of the weights you introduced doesn't introduce bias, given how you've defined them, which is justified by the conditional independence assumption at the data level, making convolution of Poisson acceptable. It's a clever approach. However, I'm concerned that at the latent process level, subsampling might alter the original spatial dependence pattern in the zeros as well as the transition from zeros to non-zero occurrences. I understand that it's costly to account for all observed locations that indeed have large proportions of zeros. But it would be beneficial to provide a brief overview of how this subsampling of zeros might impacts spatial dependencies and final conclusions. Have you attempted to analyze the sensitivity of subsampling different proportions of zeros in your final assessments? For example, changing ( p ) and ( q ) in Section 3.1 that you can easily handle computationally?

Subsampling introduces a greater variability in the model a posteriori as we use fewer data but still everything works well. Such subsampling schemes were first proposed by \cite{Koh2021}, and they also showed through a simulation study that estimation results are reliable for choices of $p$ and $q$ as in our case. We did not check in depth for the potential impacts on spatial dependencies, but if the reviewer wishes, we could do so. Since another referee also had a comment but more about our choices of $p$ and $q$, we added the following sentence to the manuscript at the end of Section 3.1: In our study, using $p = 0.5$ allows us to put the same weight on the two subsets in a partition of the zero counts of the dataset, and then setting $q = 0.9$ allows us to relatively strongly overweight
relatively large values of FWI (note that the proportion of fires occurring for FWI values greater than its $0.9$-quantile is more than $35%$ of the observed fire occurrences).

When defining the Bayesian hierarchical models in Equations (1) and (4), I found that using the notation and at the data level is confusing and not in-line with traditional Bayesian hierarchical models. According to your notation, the top layer (data level) also depends on hyperparameters that are not related to the data level. For example, conditional on the log-intensity in (1), the distribution of does not depend on the hyperparameters . Similarly, in equation (4), hyperparameters only appear at the process level, and is the only hyperparmeter that is related to the data level. I have similar comments regarding equation (5) as well where is the hyperparameter related to data level and not . I will try to remove these thetas from the data level to avoid any confusion for readers.

Thank you for pointing out that these notations were confusing. Following your recommendations, we have changed the notations in equations (1), (4) and (5).

I'm curious about the estimated spatial random effects in Figure 4; they appear very small, almost close to zero. You should provide some uncertainty estimates to determine if these effects are truly significant, and whether spatial effects are needed in the size model. Furthermore, there seems to be a dip in the estimated FA and FWI for larger FA and FWIs, which appears counterintuitive. Do you have any insights into why this might be the case?

All credibility intervals (i.e. 0.025-0.975 quantile intervals) contain zero for the estimated spatial effect, and therefore not significantly different from 0. For information, we report below these estimates. This table is not reported in the main paper but if the reviewer wishes to add it, this could be done. However, even if the spatial effect is almost zero, we still want to keep this effect in the model because it reflects a spatial variability, which we think is an important aspect to take into account in our modelling framework.

Dep	mean	0.025quant	0.975quant
1	-0.0136	-0.0909	0.0131
2	-0.0003	-0.0359	0.0345
3	0.0044	-0.0248	0.0511
4	0.0090	-0.0178	0.0732

The dip in the estimated FWI for high FWIs reflects the uncertainties associated with the estimation procedure and the choice of the spline knots. ~\
On the other hand, the dip for the estimated FA is a fairly classic feature \citep[e.g.][]{serra2014spatio,Pimont2021,Koh2021}: for very dense or large forests, fewer wildfires are observed, in particular because there is less human activity in such forests.

Minor comments

In Eqn (1), "size" is used as the superscript of $\theta^\text{year}$ which might be a typo?

Indeed, thank you for the typo.
In Figure 1, you might want to complete the sentence after "over"?
This problem is due to the conversion from the html version to the pdf version, and also certainly to the fact that a punctuation mark was missing at the end of the second subtitle. In the html version, the two graphic columns are better separated, so the first sub-heading does not appear incomplete. We will check with the editors how to solve this issue.
In line 204, $B^\text{eff}$ should be in bold.
Thank you, we also made the change in lines 173 and 230.
In line 514 references, capitalize "Safran" to "SAFRAN."
Done, thank you.

Broader Impact Concerns:

N/A
Claims And Evidence: Yes
Audience: Yes