Attributing daily ocean temperatures to anthropogenic climate change

Ocean temperatures are rising and hit record levels around the world in 2023. While trends are clear and likely strongly connected to human-caused climate change, the oceans also exhibit variability on the daily level, leading to local extremes such as marine heatwaves. We present an operational system to estimate the impact of human-caused climate change on daily sea surface temperatures anywhere in the ocean. This system uses a multi-method approach combining observed trends and paired control/forced climate model runs from CMIP6. Our approach is novel in its flexibility and ease of application for global, daily use for any day since the beginning of the satellite era (1982–2023). The system allows for rapid evaluation for further study of attributable ocean temperatures and real-time communications of attributable ongoing events. We apply the system to well-documented heatwaves in the Tasman Sea, Gulf of Maine, and Mediterranean Sea over the past decade, as well as global conditions in July 2023, to confirm that the system produces estimates consistent with other attribution methods, and to simulate how our system handles interesting events as they are occurring. Each of these events strongly reflected impacts of climate change: their temperatures were consistently made at least four times as likely to occur in our human-influenced climate than in a world without climate change. Meanwhile, in July 2023, almost all ( >70 %) of the ocean’s temperatures were made at least twice as likely to occur on any given day. Rapid attribution of daily ocean temperatures provides a pathway for quantifying the influence of climate change on ecological impacts like coral bleaching and on ocean-generated/influenced storms like tropical cyclones.


Introduction
Ocean temperatures have been rising steadily, both at the surface and at depth.In 2022, the 0-2000 m Ocean Heat Content was at a record high, having warmed a rate of 5.5 zetta-Joules per year since 1985 (Cheng et al 2023).This observed increase is one of the strongest and most robust signals of human-caused climate change (IPCC 2021).While trends in mean ocean temperatures are clearly linked to climate change, attribution studies can provide a stronger, quantitative assessment of the climate change contributions to particular events.The science of attributing changes in weather events to climate change has made considerable progress.It is now common for individual extreme heat events in air temperature to be followed by rapid analyses on the degree to which climate change changed their likelihoods (Otto 2017).Gilford et al (2022) published an approach for automating multi-method calculations underlying standard rapid attribution studies (e.g.Philip et al 2020).In this study we apply their daily attribution system to observed present-day sea surface temperatures (SSTs) globally.
Consistent with other attribution approaches (Attribution of Extreme Weather Events in the Context of Climate Change 2016, Otto 2017), our system considers temperature events as a statistical process.Statistical attribution does not evaluate the specific dynamical processes such as a shift in a front or changes in surface heat flux that give rise to an event.Rather, it calculates the likelihood of the event, regardless of the process.
Thus, it provides a consistent approach that can be used anywhere, one that can highlight conditions that have a strong climate influence and may warrant further study.The ability to make these calculations every day on an operational basis creates important opportunities for climate communication, allowing scientists or the media to draw connections between ocean temperature events and climate change.
Marine heatwaves are one of the most visible categories of extreme heat events in the ocean (Hobday et al 2016, Oliver et al 2021).These events can lead to impacts on marine ecosystems (Harvell et al 2019, Strydom et al 2020) and on human activities like fisheries (Mills et al 2013, Pershing et al 2018a, 2018b, Chen et al 2023).When these events occur, a key question is how climate change increased the likelihood of the event.Xu et al (2022) used an ensemble of climate models to show that the increase in the frequency, duration, and intensity of marine heatwaves (MHWs) around the world over the past three decades are a direct consequence of climate change.The likelihoods and intensities of individual MHWs have also been attributed to climate change.For example, Laufkötter et al (2020) calculated that climate change made the anomalously warm temperatures of the 2012 Northwest Atlantic heatwave (Mills et al 2013) at least 33 times more likely.They also calculated that the 2015-2016 Tasman Sea event (e.g.Oliver et al 2017), e.g., had temperatures that were made 50 times more likely.
For this study, we use a single methodology to analyze well-documented MHWs over the past decade in order to demonstrate how our system handles complex and extreme temperature anomalies.Furthermore, a global analysis of July 2023 (the hottest month on record to date for air temperature) demonstrates the system's ability to rapidly quantify the influence of climate change on the ocean every day, everywhere at high resolution.We use a flexible metric called the change in information due to perspective (ChIP-defined in more detail in section 2) to quantify the attributable influence of climate change on present-day MHW events.ChIP enables aggregation and averaging of attributable ocean temperatures over time and space (Gilford et al 2024).

Data and methodology
We attribute human-driven changes in daily SST likelihoods using an automated framework developed in Gilford et al (2022) [hereafter G22] based on peer-reviewed attribution protocols (Philip et al 2020).We quantify attributable SST changes using a metric called the 'Change in Information due to Perspective' (ChIP) defined as, where f (SST) are the probability distribution functions evaluating the likelihood of a SST occurring in our present day climate (mod; which is forced by approximately +1.27We draw daily historical SST observations from the National Oceanic and Atmospheric Administration's Optimum Interpolation SST version 2 (OISST) over 1981-2023 at 0.25×0.25 degree resolution (Reynolds et al 2002).OISST data provides a balance between spatial resolution of data and length of period of record, allowing us to characterize the daily variability of each grid cell, as well as its response to climate change over the past ∼40 years.We restrict our analyses to between 60 • S and 75 • N to avoid the potentially confounding effects of historical sea ice changes (Meredith et al 2019).For a given observed daily SST, ChIP is determined using SST distributions found from two distinct methods: (1) an observation-based approach that uses the linear relationships between attributable GMT changes and local observed f mod (SST) to estimate f cf (SST) without climate change, and (2) a model-based approach that directly simulates the present-day climate system with (f mod (SST)) and without (f cf (SST)) human greenhouse gas emissions.
The observation-based approach uses the empirical relationships between SSTs and global mean temperature (GMT) evolution to project a counterfactual of what today's SST distribution would be in the absence of human-caused climate change.At each ocean grid cell for the 1st and 15th day of each month we compile annually-resolved SST distributions composed of the surrounding 31 days (for a total of 24 reference distributions per year to account for seasonal cycles) and take the median and 21 evenly-spaced quantiles (between 0.01 and 0.99) of these distributions.The resulting local annual time series of median and quantile Correlation coefficients from this analysis are referred to as median and quantile 'scale factors' β q50 and β qi , respectively, which describe the linear changes in the mean and shape of the local SST distribution as a function of greenhouse gas driven global warming (figure 1, 2).
We next fit a skew-normal distribution, f reference (SST), to each of the 24 underlying reference distributions over 1991-2020 (Azzalini and Capitanio 1999) (figure 1, 1).Following the observation-based 'median method' in G22, for any given period, the attributable shift in these distributions due to global warming is found by multiplying the change in GMT between that year and the reference period by the scale factors, i.e. ∆GMT period→reference * β.This value is then added to the mean of f reference to form f mod .In this study we estimate f cf and f mod by shifting the reference period distribution (with GMT +0.9 • C relative to 1850-1899) by ∆GMT period→reference * β q50 using attributable GMT differences between the world without global warming (defined as the period from 1850-1900, with a GMT +0.0 • C relative to 1850-1899) and the climate of the analyzed period with human greenhouse gas emissions (with a GMT for a given year predicted from a regression based on the previous 30 years of GMT data for years), respectively (Masson-Delmotte et al 2021) (figure 1, 3).
We also employ the G22 'quantile method' , extending the above approach by individually shifting each of the 21 evenly-spaced quantile reference SST distributions by ∆GMT period→reference * β qi .Skew-normal distributions are then directly fit on the shifted quantiles to obtain an additional set of f cf and f mod estimates that can be used to evaluate ChIP (equation (1); figure 1, 4).ChIP estimates from the median and quantile  Eyring et al (2016), CMIP6.For each of the relevant CMIP6 models with available daily SSTs and a pre-industrial ('piControl') control run, we concatenate data from 'historical' and greenhouse gas emission forced ('ssp370' , or 'ssp585' if 'ssp370' is unavailable) runs to form continuous time series of simulated historical+forced SSTs andGMT over 1900-2050 (table 1).Model data are regridded to a common 1.5×1.5 degree grid (Zhuang et al 2023) and each simulation is individually debiased against OISST observations over a 1985-2015 calibration period (G22, see their appendix) using the method of Lange (2019).The last 30 years of the pre-industrial control experiments are taken to represent a climate without attributable human influences (GMT +0 • C), forming each model's f cf (figure 1, 5).We next compute 'crossing times' for each model by determining when the ∆GMT period→control of centered 30 year periods first reach a certain threshold of GMT change compared with its control GMT.We then computed 30 year climatologies for each 1st and 15th day of each month to form 24 distributions representing SSTs across the annual cycle (including the counterfactual), and fit skew-normal distributions to each of these.To enable flexible attribution analysis, we compute model crossing times for each tenth of a degree of ∆GMT from +0.0 to +1.5 • C.These crossing times varied significantly across models.For example, the crossing time for +1.3 • C ranged from 2002 to 2042 across the 13 models.However, because we are attempting to capture possible climates at a given ∆GMT, the difference between the rate of warming of the models matters less than the spread of temperatures observed within a 30 year period of a model.By using models with different warming rates and warming behaviors, we are able to capture a wide variety of possible 30 year climates centered on these ∆GMT.Matching the observation-based method, f mod is drawn and fit with a skew-normal distribution over each model's 30 year distribution centered around each tenth of a degree of ∆GMT from +0.0 to +1.5 • C. Then we evaluate equation (1) for each model individually and average them together to compute a final model-based ChIP estimate for each observed SST (figure 1, 6).
In each grid cell and for each daily observed SST we synthesize observation-based and model-based ChIP estimates for each observed daily SST over 2000-2023.For each target day analyzed, we determine the closest two calendar days from among the 24 analyzed annual periods (i.e. the closest two 1st and 15th days of each month) and calculate the average of their ChIP scores weighted by their relative distance from the target day.The model-based and observation-based empirical ChIP estimates are averaged together to form a single synthesis ChIP estimate (figure 1, 7).Finally, values where there has been observed cooling on the median of the climatologies (i.e.where the scale factors are negative) are filtered out to account for the presence of 'cooling holes' .This is done for ease of communication-it would be confusing in an operational system for negative temperature anomalies to be associated with positive ChIP values, and vice versa, without knowing the underlying reason.We have included a sample error calculation for July 15, 2023 in appendix.
Computing daily ChIP values from temperature anomalies and skew-normal distributions is relatively quick when computed in parallel over the world's oceans.However, calculating skew-normal distributions and scale factors is a computationally laborious process.
To demonstrate the attribution system's performance, we applied these methods to four well-documented MHW events: the Tasman Sea 2015-2016 (Oliver et al 2017), the Gulf of Maine 2012 (Mills et al 2013) and 2016   event averaged over bounding boxes used in the original studies (table 2), as well as the spatial structures of each event on the highest attributable days.The purpose of analyzing these MHWs is twofold.First, it enables us to evaluate the ability of our statistically-driven system to handle potentially complex, well-established and previously attributed events.Second, it demonstrates how daily ChIP data could be used to evaluate how individual days, regions, or seasons in a complex event are more or less attributable.Rather than focus on the specific processes associated with MHWs, we use the documented MHWs to identify persistent, strongly anomalous warm SSTs, and understand how our system characterizes these events on a day-by-day basis.We compare our results to the attribution study by Laufkötter et al (2020), which used the FAR.FAR is a function of the probability that an event occurs in the counterfactual climate, divided by the probability that it occurs in the modern climate: As such, it is easy to convert FAR values to ChIP: ) . (3) For direct comparison and where necessary we convert FAR values from other studies to the equivalent ChIP values.
We also consider the influence of climate change on the global SST patterns during July 2023, the record hottest month on Earth to date (Copernicus 2023).We conclude with an illustration of the flexibility and sensitivity of our automated approach by exploring how ChIP estimates would change if the counterfactual is instead defined by the climate (i.e. the GMT) of 2000.

Tasman 2015-2016
The 2015/2016 Tasman Sea Marine Heatwave (Oliver et al 2017) lasted 251 days from 9 September 2015 to 16 May 2016 (Oliver et al 2017).Daily SST anomalies were above 0 • C for most of 2015-2016 (figure 2, red).During the heatwave event, SST anomalies increased, reaching a maximum of 2.60 • C on December 19, 2016.All of the top 100 anomalously high days in 2016 occurred within the defined MHW period.
ChIP estimates show the influence of climate change on the daily observed SSTs (figure 2, blue), and the ChIP time series is similar to that of the SST anomaly.Anomalies in early 2015 are low and sometimes negative.However, because the 1991-2020 reference period is already warm with regards to the counterfactual period, negative temperature anomalies can correspond to positive ChIP values.Prior to the event, the ChIP was consistently >2 indicating that baseline conditions in this region were already four times more likely because of climate change, even though temperatures anomalies were negative with respect to the 1991-2020 period.During the event, ChIP increased and reached a maximum of 7.48 on 12 February 2016.This is higher than the value of 5.64 calculated by Laufkötter et al (2020), and slightly higher than the calculated upper bound of 6.64.While ChIP increased during the MHW, 15 of the 100 days with the highest ChIP over the two-year period occurred outside of the MHW.Notably, in September 2016, ChIP reached levels comparable to those during the MHW, yet the anomaly was only 1.5 • C.
The difference in behavior between anomaly and ChIP time series can be explained by looking at the probability distribution functions underpinning the attribution calculations.The difference between the mean of f mod (averaged across all methods, latitudes, and longitudes over the Tasman Sea) and the mean of f cf shows a strong seasonal pattern (figure 3, left), with a peak in the summer/early fall months (December-March), and a minimum in the Winter (May-September).A large difference between these two mean likelihoods drives higher ChIP values, as it indicates that the difference between the region before and after attributable climate change is strong.The variance meanwhile (figure 3, right), shows a relatively  modest seasonal pattern, with a maximum in the summer and a minimum in the winter.Higher variances dampen attributability by directly accounting for extreme temperature anomalies.
Taken together, the high difference in the mean likelihoods early in the year (January-April/May) boosts the ChIP values during this time of the year.Conversely, the lower difference in the mean likelihoods later in the year (September-November), paired with the higher variance, dampens the ChIP, even though temperature anomalies are similar between the two periods.
On February 12, 2016, the date with the highest average ChIP value, the SST anomaly was slightly stronger inside the main event box (figure 4   The lack of consistent correspondence between these anomaly and ChIP time series can once again be explained by looking at the probability distribution functions associated with the Gulf of Maine throughout the year.The difference between the mean of f mod and the mean of f cf shows a strong seasonal pattern (figure 6, left), with a peak from July to September, and a minimum in the spring.The variance (figure 6, right) meanwhile further drives this pattern, with both the variance of f mod and f cf peaking from May to mid July.This suggests that the lower ChIP values in early spring were driven by relatively low mean warming during this season, while the lower ChIP values in July were mostly driven by high variance in both climatologies.Meanwhile, between late May and June, both mechanisms helped contribute to the lower ChIP values, even in the presence of high anomalies.
Nonetheless, the entire event, lasting all of 2012, showed strong attributability: the spatially-averaged daily ChIP is >1 for 359 days of the year, and >2 for 223 days.
The SST anomaly on August 23 was above 2 • C across most of the Gulf of Maine region (figure 7).A maximum was apparent on the eastern edge of the bounding box, while relatively low values were found on the southern edge of the bounding box.The highest ChIP values in the Gulf of Maine bounding box were in

2016
Compared with 2012, the 30 day-average SST anomalies in the Gulf of Maine in 2016 (figure 8, red) were relatively stable and never fell below 0.5 • C. 30 day-average anomalies increased in autumn, reaching a maximum of 2.07 • C on 11 September.Daily temperatures anomalies have higher variance during the spring and summer.The daily maximum temperature anomaly of 3.11 • C occurred on 22 September The 30 day-average ChIP values closely tracks the 30 day-average temperature anomalies.Single-day ChIP values began at just above 2, steadily decreased until reaching a minimum of 0.38 on 17 May, before they increased to a maximum of 6.54 on 22 September, finally followed by a gradual decrease to the end of the year.Over the year, 321 days had a ChIP >1, and 191 days had a ChIP >2.
The spatial pattern of SST anomalies from 22 September showed a strong warming signal in the southeast of the bounding box, with some pixels having temperature anomalies of >5 • C (figure 9(a)).The rest of the bounding box contained largely positive anomalies, with a minimum off the coast of Nova Scotia.The spatial distribution of ChIP resembled that of 2012, with elevated values in the northwest portion of the domain (figure 9(b)).Compared with 2012, the pattern more closely resembled observed regional trends in temperature; the median scale factor values were >4 • CSST • CGMT in this region (figure 9(c)).

Mediterranean 2022
The 2022 Mediterranean MHW began in May and persisted through at least the end of the year (Martínez et al 2023).SST anomalies (figure 10, red) increased markedly in late spring and remained elevated through the remainder of 2022.The average SST anomaly prior to 15th May was 0.06 • C, while the average after 15 May was 1.13 • C. The minimum SST anomaly of −0.50 • C occurred on 5 April, shortly before the onset of  warming.The 30 day maximum of 1.45 • C occurred on 31 July, while the daily maximum of 1.77 • C occurred on 26 July Prior to the start of the event, ChIP values (figure 10, blue) closely followed the SST anomalies.After this, the 30 day-average ChIP increased steadily throughout the summer to a maximum of 4.5 on 3 September.The daily maximum of 5.01 occurred on 12 September.While SST anomalies were steady and elevated, the ChIP time series decreased to around 2 throughout October and November.ChIP rarely fell below 2 during the event, indicating a high baseline of attributability throughout the documented MHW.
The difference between the mean of f mod and the mean of f cf shows a strong seasonal pattern (figure 11, left), with a peak in the summer months, and a minimum in the spring.The variance of both the probability distribution functions follows a similar pattern, with a high level of variance across both climates in the summer, and a minimum in March.However, the variance values drop off steeply after July, while the difference in mean values drops off much more slowly.Higher difference in mean in the summer drives higher ChIP values, while higher variance drives lower ChIP values.This explains the increase in ChIP values between June and September paired with relatively constant SST anomalies: as the variance decreases-while the difference between the likelihood means stays relatively elevated-similar SST anomalies become more attributable to climate change.
On the 12 September, the SST anomaly was positive throughout most of the Mediterranean (figure 12).A maximum of 3 • C occurred in the western Mediterranean, near the Balearic islands.The only negative SST anomalies occurred in the Ionian Basin.The eastern half of the Mediterranean generally experienced lower anomalies than the western half.
The map of ChIP strongly resembles the anomaly map (figure 12).The same maximum occurred near the Balearic islands, with values of ChIP reaching 15, while a minimum occurred in the Ionian basin.Although negative temperature anomalies occurred in this portion of the domain, ChIP values were above   zero indicating a certain but limited level of attributability.This suggests that the Ionian basin was not in a MHW at this point in time, but the averaged ChIP value over the entire region (figure 10) were still elevated, indicating that even outside the MHW period, regional conditions remained meaningfully impacted by climate change.

July 2023
To date, July 2023 was the hottest recorded month on Earth (for surface air temperature) (Copernicus 2023), and elevated ocean temperatures were a major driver of record global temperatures.Positive anomalies occurred over most of the world's oceans (figure 13).Strong positive values associated with a growing El Niño were apparent in the upwelling region west of Peru and along the equator.Prominent positive temperature anomalies were also evident in the Northeast Pacific, off of Japan, in and to the south of the Labrador Sea.Areas with negative anomalies were apparent to the west of the of the US and Mexico, in portions of the Arctic, and in several places in the southern hemisphere around Australia and southern Africa.
Several of the areas with prominent SST anomalies showed a strong signal in ChIP.These include the positive anomalies in the Labrador Sea and the northeast Pacific, both with ChIP >5, and the negative anomalies in the eastern North Pacific and around Australia.In contrast, several areas with strong positive anomalies had only modest ChIP values.These include the eastern equatorial pacific.This corresponds to the El Niño upwelling region.The strong temperature anomalies associated with the onset of 2023's El Niño did not correspond to proportionally elevated ChIP values (consistent with a low signal to variance in this region).Finally, some of the highest ChIP values occurred in areas with only modest SST anomalies.These include the Caribbean, the western Indian Ocean, and the western equatorial Pacific.Additionally, there were a handful of regional swaths of uncomputed ChIP, most strikingly off the west coast of South America, which had some of the strongest SST anomalies.This is mostly caused by the prevalence of negative scale factors in those areas, which are filtered out.The average spatially-weighted daily ChIP in July was 2.42-this indicates that the average daily SST for a given location was made 5.35 times as likely in a world with human-caused climate change.Furthermore, only 32% of the ocean had a ChIP <1.
Additionally, the anomalous heat of July 2023 shows attributability even when computing ChIP substituting a counterfactual based on the GMT from the year 2000 (+ • 0.79C, figure 14).There was still widespread signal apparent, as the high ChIP values northeast Pacific, the areas to the East of Papua New Guinea, the Caribbean, and the Labrador Sea stand out.Patches of the southern Pacific and Indian Oceans also emerge as notable regions.The average area-weighted ChIP for the month (with a counterfactual baseline centered on 2000) was 0.25, with a much higher portion of the oceans showing less or little attributability.

Discussion
Rapid attribution of climate-influenced events over land have become more common.We apply an analogous approach to attribute SSTs at a daily, global scale to human-caused climate change.We illustrate this approach with well-studied marine heatwaves to validate the effectiveness of our approach and to complement studies such as (Frölicher et al 2018) that consider the changing frequency of these events at a global scale.That study found that anthropogenic climate change has doubled the frequency of marine heatwaves.Consistent with Laufkötter et al (2020), our study indicates that the influence of climate change on particular events was even stronger, with ChIP values exceeding 6 during each of the MHWs we considered (implying they were made 64 times more likely because of climate change).This was even the case when ChIP values were averaged over regions not necessarily in a MHW (such as the Ionian basin in the Mediterranean MHW).Our results emphasize the robustness of the ChIP metric for quantifying attributability over time and space on daily timescales.
Daily attribution of both extreme events (figures 2-10) and global patterns in July 2023 (figure 13) underscores the high level of attributable warming observed across the world's oceans.Daily ocean ChIP values outside the defined marine heatwaves were rarely below 1 (with an exception during a brief period of negative ChIP and anomalies in the Mediterranean Sea in 2022).In July 2023, only 32% of the ocean had an average monthly ChIP of less than one.
The strong levels of attribution in the ocean reflects two of its properties.First, ocean warming is one of the strongest and most reliable signals of climate change, with ocean heat content rising steadily (Cheng et al 2023).Second, the high heat capacity of the ocean smooths out high frequency signals coming from the atmosphere.The consistent trends and reduced variability increase the detectability of the ocean's climate change signal.
Our analysis also highlights another property of short-term attribution studies: attributability does not strictly scale with temperature anomalies.This is most striking in the Gulf of Maine 2012 example (figure 5): whereas the most intense anomaly (2.41 • C on 21 May) was attributable (ChIP >2), less intense anomalies later in the summer had much higher ChIP values.The stronger attributability in the late summer and fall is driven by the stronger historical rates of warming during this season compared with other times of the year (Thomas et al 2017), paired with a higher variance both before and after global warming.The impact of variance on the difference between anomaly and attributability is further seen in results from July 2023 (figure 13).Due to the emerging El Niño, the eastern equatorial Pacific had one of the most prominent temperature anomalies in July 2023 (figure 13).However, due to the El Niño Southern Oscillation, this region has very large observed interannual variability (see Witman et al 2023).This means that the ChIP values in this region are modest relative to other places like the Caribbean, where temperature variability is lower.
This effect can also be seen when comparing ChIP and temperature anomaly values across MHWs.For example, when comparing the two Gulf of Maine MHWs, one might look at the temperature anomalies across both events, and assume that the elevated SST anomalies from April to July 2012 might be more attributable to climate change than the cooler SSTs in 2016 (figures 5 and 8).However daily ChIP values across both events generally stay between 1 and 2 between April and July due to elevated variance associated with this time period (figure 6).
ChIP values can also be compared across locations.Comparing the two Gulf of Maine events to the Tasman Sea MHW, the Gulf of Maine events tended to have higher temperature anomalies, but the Tasman Sea MHW had higher ChIP values (figures 2, 5 and 8).This illustrates how the ChIP (in concurrence with SST anomaly) can be quickly used to understand how different regions have reacted differently to climate change.
The Gilford et al (2024) daily attribution approach applied in this study builds on the extreme event method developed by Philip et al (2020).These methods depend in part on the relationship between local temperature and warming.Over the 42 year OISST record, GMTs have risen +0.9 • C (70% of the total warming since 1900) suggesting that OISST should be able to resolve much of the local relationships with global warming.However, the relatively short time series creates challenges for the empirical part of the attribution estimates.The relevance of interdecadal variability creates a challenge for ocean temperature attributability, a problem which is less prevalent over land.For example, the region off of the coast of Chile that led us to mask out these pixels in the July 2023 map (figure 13) are likely due to the relatively cool conditions in recent years related to the Pacific Decadal Oscillation.It is simultaneously possible that natural climate variability has amplified the warming trends in some locations.Incorporating longer time series like NOAA's ERSST (Huang et al 2015) would reduce the impact of decadal variability on our attribution quantification.However, empirical approaches will always have some sensitivity to natural variability; one strategy to reduce the influence of this sensitivity would be to use large ensembles of climate models (e.g.Alexander et al 2018, Xu et al 2022) to estimate distributions of local warming rates and then blend these estimates with those from OISST using Bayesian approaches (e.g.Robin and Ribes 2020).
A future challenge for our system, and for attribution studies generally, is the concept of emergence (Henson et al 2017).When a climate warms to such an extent that it no longer overlaps with the counterfactual, the resulting ChIP values become infinite, making them functionally meaningless.Two solutions present themselves using our system: first, the counterfactual climate can be readjusted to a different time period (as in figure 14).Second, the methods outlined in this paper could be adapted to compute counterfactual temperatures, representing what a location's temperatures would have been in the counterfactual climate.This can be done by recomputing equivalent quantiles from current climatologies on counterfactual climatologies.
Computing counterfactual temperatures creates new opportunities to attribute a variety of more complex but equally human-relevant processes.Each of the marine heatwaves considered above had notable impacts on marine ecosystems, and coral bleaching in the Caribbean has been documented during the exceptionally warm conditions in 2023 (NOAA 2023).Coral bleaching risk is calculated using a degree day model based on a threshold defined by the maximum climatological monthly temperature (Liu et al 2014).Our approach could be used to estimate the ChIP for each day above the threshold and then comparing the number of degree days in the current and counterfactual climates, providing an index of coral bleaching attributability.Our approach could also be extended to estimate the influence of human-caused climate change on tropical cyclones-several hurricanes during 2023 underwent rapid intensification off American coastlines while passing over waters that were strongly influenced by human-caused climate change (i.e.ChIP >3), consistent with known links between anomalous ocean temperatures and intensification (e.g.Kaplan andDeMaria 2003, Bhatia et al 2019).While tropical cyclone attribution is complex and will require additional analyses, our work highlights potential paths forward for considering how attributably warmer ocean temperatures could affect these storms.

Appendix. Error calculations
Here, we calculate the standard error for a sample day (15 July 2023), to help convey the order of magnitude associated with the errors for ChIP calculations.
First, we select the median scale factors associated with 15th July, and retrieve the error around the linear fit performed to compute the median scale factors.Using this, we are able to find the two-sigma confidence interval around the scale factors.From this, we are able to find the confidence interval for the ChIP values by computing ChIP values using the maximum and minimum values in the scale factor confidence interval, before converting to a Standard Error.
Next, for the model-based method, we find the standard error by dividing the mean ChIP by the square root of 13 (the number of models used to find the mean model ChIP).
We combined the two results, and found that for 15th July, the area-weighted average error was 0.66, while the area-weighted ChIP was 2.42.Below are the global maps associated with each value for 15th July (figure A1)

Figure 1 .
Figure 1.Implementation of the Gilford et al (2022) multi-method approach.(1) The reference climatology of daily temperatures at a location and time of year is represented by a skew norm distribution fit to the daily data from 1991-2020.(2) The scale factor (β) for the location is calculated by regressing the local SSTs against the global mean temperature for 1982-2020.(3) The scale factor is used to shift the reference distribution to the modern (red) and counterfactual (blue) temperatures.(4) ChIP(SST) is calculated from the resulting empirically-derived likelihoods of the observed SST.(5) An alternative modern climate is constructed from distribution of temperatures surrounding the year when a climate model's global mean temperature crosses the desired value (+1.3 • C in this example).A corresponding counterfactual climate is built from the model's control run.(6) ChIP(SST) is calculated from the model-based likelihoods.(7) The final ChIP(SST) is computed as the mean of two observation-based empirical estimates and estimates from the 13 paired-models.
(Pershing et al 2018b), and the Mediterranean Sea 2022(Martínez et al 2023).For each event, we consider the time series of SST anomalies (relative to 1991-2020) and ChIP during the documented

Figure 2 .
Figure 2. Time series of spatially-averaged temperature anomalies (relative to 1991-2020, red) and ChIP (blue) during the 2015-2016 Tasman Sea MHW.The thick lines are the daily averages, the thin lines are a 30 day running mean.The dashed vertical lines indicate the beginning and end of the event, as described in Oliver et al (2017).

Figure 3 .
Figure 3.The difference between the mean values of the modern probability distribution functions f mod and the counterfactual probability distribution functions f cf , averaged over the Tasman Sea (Left), and the variance of f mod and f cf (Right).
(a)).It featured a prominent maximum in the southern part of the domain and lower values in the north.ChIP was more clearly elevated in the event box, exceeding 15 in some locations in the north (figure 4(b)).In general, the spatial structure within the domain differed between ChIP and SST anomalies.Note that the spatial structure in ChIP is a combination of both the signal (i.e. the anomaly) and the strength of the associations with climate change (i.e.changes in likelihood).The median scale factor (a measure of how reactive a location has been to climate change for that time of year) is one component of this sensitivity (figure 4(c)); it had higher values in the eastern edge of the domain associated with a local ChIP peak.3.2.Gulf of Maine 2012, 2016The Gulf of Maine experienced two well-documented and similar MHWs(Pershing et al 2018b) in 2012(Mills et  al 2013) and 2016.

Figure 4 .
Figure 4. Temperature anomaly, ChIP, and scale factor for the day with the peak average ChIP (2 February 2012).Bounding box described in table 2.

Figure 5 .
Figure 5.Time series of spatially-averaged temperature anomalies (relative to 1991-2020, red) and ChIP (blue) during the 2012 Gulf of Maine MHW.The thick lines are the daily averages, the thin lines are a 30 day running mean.

Figure 6 .
Figure 6.The difference between the mean values of the modern probability distribution functions f mod and the counterfactual probability distribution functions f cf , averaged over the gulf of Maine (Left), and the difference between the variance values of f mod and f cf (Right).

Figure 7 .
Figure 7. Temperature anomaly (relative to 1991-2020), ChIP, and scale factor for the day with the peak spatially-averaged ChIP during 2012 in the Gulf of Maine (23 August 2012).Bounding box described in table 2.

Figure 8 .
Figure 8.Time series of spatially-averaged temperature anomalies (relative to 1991-2020, red) and ChIP (blue) during the 2016 Gulf of Maine MHW.The thick lines are the daily averages, the thin lines are a 30 day running mean.Bounding box described in table 2.

Figure 9 .
Figure 9. ChIP, temperature anomaly (relative to 1991-2020), and scale factor for the day with the peak spatially-averaged ChIP during 2016 in the Gulf of Maine (22 September 2016).Bounding box described in table 2.

Figure 10 .
Figure 10.Time series of spatially-averaged temperature anomalies (relative to 1991-2020, red) and ChIP (blue) during the 2022 Mediterranean MHW.The thick lines are the daily averages, the thin lines are a 30 day running mean.The dashed vertical line indicates 1 May.

Figure 11 .
Figure 11.The difference between the mean values of the modern probability distribution functions f mod and the counterfactual probability distribution functions f cf , averaged over the Mediterranean (Left), and the difference between the variance values of f mod and f cf (Right).

Figure 12 .
Figure 12.Temperature anomaly (relative to 1991-2020), ChIP, and scale factor for the day with the peak spatially-averaged ChIP during 2022 in the Mediterranean (12 September 2022).

Figure 13 .
Figure 13.SST anomaly (top) and Average ChIP (bottom) for July 2023.Dark Blue indicates regions where the median scale factor is negative (cooling has occurred while GMT has increased).

Figure 14 .
Figure 14.SST anomaly (top) and mean ChIP (bottom) for July 2023, with a counterfactual baseline of 2000.Dark blue indicates regions where the median scale factor is negative (cooling has occurred while GMT has increased).

Table 1 .
List of model names and forced Shared Socioeconomic Pathway (SSP) used.

Table 2 .
Details of the four marine heatwaves considered in this study.