Evolution of causal relationships under climate change: controls of Net Primary Productivity in the North Altantic Subpolar Gyre

Bénard, Germain; Gehlen, Marion; Vrac, Mathieu

doi:10.5194/egusphere-2025-4680

Preprints

https://doi.org/10.5194/egusphere-2025-4680

Preprints

13 Nov 2025

| 13 Nov 2025

Status: this preprint is open for discussion and under review for Earth System Dynamics (ESD).

Evolution of causal relationships under climate change: controls of Net Primary Productivity in the North Altantic Subpolar Gyre

Germain Bénard, Marion Gehlen, and Mathieu Vrac

Abstract. Understanding how climate change affects marine primary productivity requires examining the evolving causal relationships between physical and biogeochemical processes. We applied the PCMCI+ causal discovery algorithm to investigate how the mechanisms controlling Net Primary Productivity (NPP) in the North Atlantic Subpolar Gyre evolve under different climate scenarios across five Earth System Models. Using 100-year sliding windows, we compare causal relationships in future scenarios against pre-industrial conditions, focusing on the roles of mixed layer nutrients, vertical mixing and horizontal transport. Our analysis reveals three main categories of relationship evolution: the emergence of links, the disappearance of links, and changes in link strengths. For example, while the link between stratification and NPP emerges under climate change in CanESM5-CanOE, it strengthens in CMCC-ESM2 and remains stable with moderate to high intensities in other models. At the end of the 21st century, the spread between models regarding the effect of stratification on NPP is reduced compared to pre- industrial conditions, suggesting a reduction in inter-model uncertainty. However, the transport and vertical mixing controls on the supply of nutrients to the mixed layer exhibit a more diverse evolution among the ESMs studied. The CMCC-ESM2 model has a strengthening of the relationships between winter vertical mixing and nutrients, while IPSL-CM6A-LR and CanESM5CanOE show weakening of these relationships. Furthermore, the evolution of the link between nutrient supply to the mixed layer for NPP exhibits a large variability between models. These divergent pathways reveal that the dynamics of nutrients has uncertain evolution between models. Lastly, model-specific dynamics are also observed, such as the strengthening of the link between horizontal transport and the nutrient content of the mixed layer in IPSL-CM6A-LR. Together with the decreasing strength of the vertical mixing/nutrients link, this suggests the presence of compensation mechanisms and a shift from vertical mixing dominance to enhanced horizontal transport control over the course of the scenario. These findings offer mechanistic insights into the dynamics of ESMs, specifically in the evolving relationships between physical and biogeochemical processes that shape the projections of NPP and nutrients. The causality-based approach identifies mechanisms that traditional analyses miss, offering a novel framework for model intercomparison and understanding ecosystem responses to climate change.

Received: 24 Sep 2025 – Discussion started: 13 Nov 2025

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.

Download & links

Preprint (PDF, 1518 KB)

Supplement (1146 KB)

Download & links

Germain Bénard, Marion Gehlen, and Mathieu Vrac

Status: open (until 05 Jan 2026)

Post a comment Subscribe to comment alert

RC1:
'Comment on egusphere-2025-4680', Anonymous Referee #1, 11 Dec 2025 reply
This paper aims to investigate the evolution of mechanisms that control NPP in the North Atlantic Subpolar Gyre under three emission scenarios for five ESMs using causal analysis in a sliding window approach. It builds on a previous study by the same authors in which model differences of said mechanisms were analysed with causal analysis on piControl runs, i.e. under absence of anthropogenic forcing.
The paper is well-structured and the question of causal stationarity of drivers of North Atlantic NPP under global warming in different ESMs is relevant and interesting. The sliding window approach is sensible to address this question, with possible shortcomings discussed, and the findings, although limited to a small number of ESMs, give insights into inter-model differences and non-stationarities in the relationships between physical and biogeochemical oceanic variables.
The Discussion section is well written, the Limitations as well as the advantages of the presented approach are contextualised appropriately, and the relationship to other methodologies such as emergent constraints are interesting and timely. The Conclusion offers a good summary and outlook.
However, I have a few general methodological concerns and questions that should be addressed:
According to your Methods description, your causal analysis consists in applying the causal discovery algorithm PCMCI+. To my knowledge, causal discovery is only really suited to obtain the qualitative causal model; the actual causal effect strength needs to be estimated in a subsequent causal effect estimation (see Runge et al. Nat Comms 2015). It is unclear to me, whether you actually perform such a causal effect estimation or whether you simply use the coefficients of the causal discovery / PCMCI+ algorithm as values to describe causal effect strength. The Tigramite package provides functions for the causal effect estimation as well, for example LinearMediation, in which actual path coefficients are derived through a linear regression of each link to all parent variables. If you did in fact use this in your analysis, please add a description of this to the Methods, so it is clear to the reader. If you did not, I strongly recommend doing so, as the coefficients taken from PCMCI+ directly may be misleading.

Your Methods are very brief on how you determine statistical significance of the Links and it is a bit unclear to me; significance of Links vs significance of changes in Links are also not clearly discernible to me. This makes it difficult for me to interpret Tables 1-3, and the lack of significance of the vast majority of points in Figure 3. Further, are the Intensity-values in the tables obtained as averages of all Link strengths or do you leave the non-significant links out of the calculation, treating them as having disappeared? In Figure 1, it is unclear to me what you mean with “the average minimum significant strength” of 0.31 - was this measured across the different links (some of which may be considerably stronger than others) and what is the reasoning behind using this as a threshold for all? And once more, those links that were not marked as significant, but still shown in the Figure, how would you interpret those? Please provide a clear description of your approach.

I would like to positively highlight the sensitivity tests you performed. The application of the sliding window analysis to the piControl data is important as a reference. Here I noticed that the spread of link strengths you find in the unforced data is considerable in Figures 2 and 3. In case of the link between stratification and NPP even the sign of the causal effect in the piControl runs is not consistent across windows. It would be good to address this in the Discussion and think through potential reasons and implications.

Your tables contain a lot of information and are generally difficult to read. Each Table would benefit from a title that makes it clear at first sight which link it describes, this should then easily be identifiable in the according Figures of the causal graphs. It is unclear to me why you provide nutrient-resolved information in Tables 1 and 3, but flatten the nutrient-level information in Table 2 and instead provide ranges. This inconsistency is a bit confusing. In general, I would encourage rethinking which information here is really the most important and how you can make the main takeaways in terms of changes in links under the different scenarios easily accessible to the reader; possibly you could move the full tables to the Supplementary and provide simplified versions for the main text.

Figures 2+3: Please add model and nutrient to the Figure, for example as a title, and please include the full link, not only one variable, so it is clear at first sight what the respective Figure is describing. There are no networks shown (neither in the main text nor in the Supplement) that contain the variable “no3”. This would be good to add. Further, it would be good to visually set aside also the points that are identical between all three rows, i.e. that span the historical runs only, like you did with the last 35 that only span scenarios. I also suggest adding a panel that shows GMT over time for the three scenarios as a reference.

Your causal networks consist of contemporaneous and, hence, undirected links; when discussing your results, however, you do describe causal influences in a directional manner. I assume that this interpretation of the links as oriented is based on expert knowledge; it would benefit the paper if you described the reasoning behind this. The same holds for the choice of timescales of the analysis; if all links you find are contemporaneous, this suggests that the time resolution chosen might be too coarse to resolve lagged links, which would more finely describe the actual physical processes. What was your reasoning behind this choice and might there be potential for more fine-grained analysis in the future? Analysing sub-annual data for this current paper certainly exceeds the scope, and you find novel and valuable information on the chosen timescale that is worth publishing. However, the paper would benefit from some reflections on the point of time resolution and future research potential in the Discussion/Perspective section.

Throughout, please unify nomenclature: You use “Intensity”, “Link intensity”, and “Link strength” interchangeably as far as I can tell. In Figure captions, please make sure to always add the abbreviations of variables that are used in the Figures to the variable names.

For the sake of reproducibility, please provide your code and the necessary aggregated time series to run it in a public repository.

Furthermore, please find below a few specific minor comments and suggestions:
L 7-9: I suggest to re-formulate this sentence from “Our analysis reveals three main categories of relationship evolution” to “Our analysis detects three main types of relationship evolution".
Section 2.3: Please add a short description of the piControl runs, specifically also their respective length.
L 187-188: “the proportion of points outside the pre-industrial distribution” would correctly be “the number of points…”; you could consider alternatively giving the proportion in %
L 345-347: Do I understand correctly that the external forcing narrows the uncertainty of this link? If this interpretation is accurate, it might make sense to specify this here, acknowledging that future ocean productivity might once more behave differently under emission scenarios with less forcing, such as Zero Emission Commitment scenarios, which could be interesting for future analysis.
L 409-411: As your final conclusion, you present the suggestion to conduct similar analysis on observational data, also stating that data availability may be limited. I agree that the study of observations is always desirable - however, often not possible for the reason you describe. If you could add a sentence or two on the actual situation, i.e. what variables would actually be available, how long and in what time resolution those records are available etc. and if you can derive any recommendations for measurement series from your study design and findings, that would be helpful.

Reply
Citation: https://doi.org/10.5194/egusphere-2025-4680-RC1
RC2:
'Comment on egusphere-2025-4680', Anonymous Referee #2, 12 Dec 2025 reply
The paper “Evolution of causal relationships under climate change: controls of Net Primary Productivity in the North Atlantic Subpolar Gyre” by Bénard et al. uses a causal discovery algorithm to investigate how causal drivers of the Net Primary Productivity in the North Atlantic Subpolar Gyre evolve with time and especially under climate change. It compares causal graphs and causal effects between pre industrial and future time periods, under three different emission scenarios and across five Earth System Models.
The paper is well-written and clearly structured. The findings (e.g., large variability between models for some links) support the relevance of the scientific question. The limitations and discussion section effectively highlights the significance and implications of the results.

I would first like to raise a few points, which are as follows:
I could not find a link to the code used for the analysis, only to the Tigramite repository. For reproducibility, it would be helpful to have access to this code.

It seems like the causal strength/intensity you are referring to is the value of the causal effect (l71: “PCMCI+ enables the quantification of causal relationships”). However, this is not what PCMCI(+) outputs with that color bar, but rather the cross-MCI value, which is a test statistic that sort of represents the algorithm's confidence that the link exists, but this is limited by the quality of the algorithm (see https://github.com/jakobrunge/tigramite/blob/master/tutorials/causal_discovery/tigramite_tutorial_pcmciplus.ipynb). It is not directly comparable between links because some links undergo more tests than others. To determine the value of a causal effect, one needs to perform causal effect estimation with a known graph (here it would be the one obtained using PCMCI+) (https://github.com/jakobrunge/tigramite/blob/master/tutorials/causal_effect_estimation/).

Only stationarity is discussed as an assumption required for the use of PCMCI+. A (short) description of other assumptions and of the conditional independence test used would be valuable (causal sufficiency, faithfulness, (non)linearity of relationships influencing the choice of test, etc.), as these can influence the results.

On Fig1 and the corresponding supplementary figures, the meaning of the double arrows ‘↔’ is not given, though it seems important. Does it mean that the algorithm found effects from X to Y and from Y to X? The lags then might differ? I strongly recommend not using double arrows in that case because in causal graphs they are usually used to represent hidden confounding. Moreover, the lags of causal effects are not discussed nor appear on the figures. This might be intentional, but should then be justified as it is also a parameter that could change with climate change, and because the defined maximum lag (3 years) sets a limit beyond which an effect can disappear from the analysis (but could for example exist at a longer lag).

A more extensive description of how you assess the significance of links would be very helpful. You say at l163 “The statistical significance of changes is assessed using both the Cramer-von Mises test (Anderson, 1962) for distributional changes and t-tests (Kim, 2015) for shifts in mean values.” Could you link this sentence to the figures? What are you measuring with distributional changes as opposed to shift in mean values? For example on Fig2 I am not sure to understand what the boxes with values of the statistic and pvalue refer to. You might be referring to the distributional change between the grey and red histograms. It could help the reader to have a description of what these values are and how they are computed in the caption.

I also listed some more specific comments:
l75 It is not clear to me what you refer to with “within a predefined conceptual scheme linking the physical and biochemicals variables”. You could clarify or reformulate.

l84 “base on Granger causality” is ambiguous, I would suggest to reformulate the sentence as: “PCMCI+ is a causal discovery algorithm for time series that relies on conditional independence tests grounded in Granger-causal principles. It identifies both contemporaneous and time-lagged relationships between multiple variables.”

l88 Suggestion for something a bit more accurate to describe PCMCI+: “The algorithm proceeds in two main steps: (1) the PC (Peter & Clark) step, which iteratively tests conditional independence between variables to eliminate spurious links and obtain a set of candidate parent variables; and (2) the MCI (Momentary Conditional Independence) step, which evaluates each candidate link by testing conditional independence while conditioning on the parents of the target variable and a relevant subset of the parents of the driver.”

l140 “The graphs of causal relationships will be explored separately for each nutrient (nitrate, dissolved iron, silicate and phosphate).” Why do you choose not to include all nutrients in the same causal graph? Depending on the relationships between these nutrients, one could identify effects that are actually mediated effects. The comparison between nutrients is then not straightforward. It may nonetheless be justifiable to separate nutrients depending on expert knowledge. It would be good to explain the rationale behind this choice in more detail.

l157 “We examined the performance of PCMCI+ according to the length of time-series with auto-regressive vectors.” What do you mean here? Did you simulate some data and applied PCMCI+ using different window lengths in your sensitivity analysis? If so, I would suggest to do a sensitivity analysis directly on the data used for the study, as it might have different characteristics than those of simulated auto-regressive vectors and the results might strongly depend on the window size.

l159 I don’t understand what you mean with 'estimation errors in link strength', especially what the ground truth is (i.e., what the error is measured against). You could add a short explanation.

l162 It would be very valuable to describe the pre-industrial control simulation, especially the years covered. For example, how do you obtain at l. 175 “one causal graph obtained under pre-industrial conditions” and “the mean causal graph representative for each SSP scenario”?

l176 “ …the mean causal graphs representative for each SSP scenario. The latter are derived from the 35 last sliding windows”. I would suggest to add a slightly more extensive explanation on how the “mean causal graph” is obtained. Are ”link strengths” simply averaged over each graph?

It would be clearer to use one terminology throughout the paper for the quantification of causal effects (and not strength, intensity, etc). Although at present it is not the value of the causal effect that is derived (see general comments), if you later perform causal effect estimation, I recommend defining the causal effect in your setting (see e.g. https://github.com/jakobrunge/tigramite/blob/master/tutorials/causal_effect_estimation/tigramite_tutorial_general_causal_effect_analysis.ipynb) and choosing only one word to refer to it.

And finally, some minor technical corrections:
typo in the title for Atlantic

l176 causal graphs —> causal graph

Fig1 The arrows are a little too small (especially arrow heads).

l252 5x10-3 —> 5x10^-3; same l256 to be consistent with the rest of the text

l355 point missing after ‘framework’

Reply
Citation: https://doi.org/10.5194/egusphere-2025-4680-RC2

Germain Bénard, Marion Gehlen, and Mathieu Vrac

Supplement

https://doi.org/10.5194/egusphere-2025-4680-supplement

Germain Bénard, Marion Gehlen, and Mathieu Vrac

Viewed

Total article views: 232 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	Supplement	BibTeX	EndNote
165	45	22	232	33	14	13

HTML: 165
PDF: 45
XML: 22
Total: 232
Supplement: 33
BibTeX: 14
EndNote: 13

Views and downloads (calculated since 13 Nov 2025)

Month	HTML	PDF	XML	Total
Nov 2025	123	18	10	151
Dec 2025	42	27	12	81

Cumulative views and downloads (calculated since 13 Nov 2025)

Month	HTML	PDF	XML	Total
Nov 2025	123	18	10	151
Dec 2025	42	27	12	81

Viewed (geographical distribution)

Total article views: 230 (including HTML, PDF, and XML) Thereof 230 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 22 Dec 2025

Short summary

Climate change is transforming ocean ecosystem dynamics. We used causality methods to study how the relationships between ocean physics and biogeochemistry will change in the North Atlantic over the next century. By analyzing five different climate models, we discovered that environmental drivers' influence on ocean productivity evolves in complex ways under global warming. One environmental driver can become of major importance while others can become irrelevant.


Total:	0
HTML:	0
PDF:	0
XML:	0